Beyond the Tree: Integrating Gene Flow into the Multispecies Coalescent Model

Abstract

The multispecies coalescent (MSC) model provides a powerful framework for inferring species histories from genomic data. However, its standard form assumes complete isolation after speciation, a premise often violated by gene flow. This article explores the extended MSC models that explicitly incorporate gene flow and introgression (MSC-I, MSC-M). We cover the foundational principles of gene tree-species tree discordance, detail current Bayesian methodologies and software implementations for parameter estimation and species delimitation, address prevalent challenges like model violation and over-splitting, and present empirical evidence validating the MSC over concatenation approaches. Aimed at researchers and scientists, this synthesis highlights how accounting for gene flow leads to more accurate inferences of divergence times, population sizes, and species boundaries, with critical implications for understanding evolutionary history and genomic diversity.

From Basic Coalescent to Reticulate Evolution: Core Concepts of the MSC with Gene Flow

The Multispecies Coalescent (MSC) model is a population genetics framework that describes the evolutionary relationships of genes sampled from multiple closely related species. It provides a powerful paradigm for phylogenomics by modeling how genealogical histories of individual genes (gene trees) relate to the broader species phylogeny (species tree). A key insight of the MSC is that gene trees can differ from the species tree and from each other even in the absence of hybridization or horizontal gene transfer, primarily due to incomplete lineage sorting (ILS) [1] [2]. ILS occurs when genetic lineages from the same population fail to coalesce (find a common ancestor) in that population and instead coalesce in a more ancestral population, resulting in gene tree topologies that are discordant with the species tree topology [2]. The MSC has become an essential tool for accurate species tree estimation, particularly in recent years as genome-scale data has revealed that genealogical discordance is the rule rather than the exception across the tree of life, especially in radiations with short internal branches where ILS is pervasive [2] [3].

Theoretical Framework and Key Concepts

The Model and Its Biological Basis

The MSC model describes a backward-time process where lineages coalesce within a population-level species tree. Each species tree (\mathcal{T}=(T,\Theta)) comprises a topology (T) and branch lengths (\Theta) denominated in "coalescent units," which are units of time normalized by population size [1]. The model operates by tracing the history of individual gene lineages backward through time: at every leaf of the species tree, a lineage begins and grows backward until it reaches a speciation event. Upon entering a common population with lineages from a neighboring branch, these distinct lineages become eligible to coalesce according to a Poisson process with a constant hazard rate (\lambda) [1]. The time (\tau{ij}) until two lineages (i) and (j) coalesce follows an exponential distribution: (\tau{ij}\sim f{\lambda}(\tau{ij})=\lambda e^{-\lambda \tau{ij}}). In a population with (k) distinct, uncoalesced lineages, the time until the next coalescent event also follows an exponential distribution with a rate of ({{k}\choose{2}}\lambda): (\tau{k\rightarrow k-1}\sim f(\tau)=\frac{k(k-1)}{2}\lambda e^{-\frac{k(k-1)}{2}\lambda \tau}) [1]. This process continues until all lineages have coalesced into a single one, forming a binary gene tree.

Incomplete Lineage Sorting and its Consequences

Incomplete lineage sorting represents a major biological cause of gene tree discordance. ILS occurs when the coalescence of gene lineages predates speciation events, causing the gene tree topology to differ from the species tree topology [2]. The probability of ILS increases when internal branches of the species tree are short (in coalescent units) relative to effective population sizes, which is common in recent, rapid radiations [2] [3]. The interaction between coalescence and speciation can generate substantial incongruence between gene trees and species trees, creating challenges for phylogenetic inference that have led to the development of specialized species tree estimation methods that account for this discordance [1] [4].

Hemiplasy: A Key Conceptual Implication

A critical consequence of ILS is hemiplasy—the appearance of homoplasy (convergent evolution) resulting from substitutions on discordant gene trees rather than true independent origins [2]. When a character-state transition occurs on a discordant gene tree, analyzing it against the species tree topology can suggest multiple independent origins when only a single transition actually occurred [2]. This phenomenon has been shown to affect phylogenetic inferences, potentially causing apparent substitution rate variation, spurious signals of positive selection in coding sequences, and artefactual signals of convergence [2]. Hemiplasy can impact both discrete traits and, through the combined effect of multiple loci, quantitative traits as well [2].

Species Tree Estimation Methods Under the MSC

Statistically Consistent Estimation Methods

Numerous methods have been developed to estimate species trees in the presence of ILS, with many proving to be statistically consistent under the MSC model—meaning they converge to the true species tree topology as the amount of data increases [1]. These include:

Table 1: Coalescent-Based Species Tree Estimation Methods

Method	Input Type	Output	Key Features
ASTRAL-I/II/III [1]	Gene trees	Species tree topology, branch lengths	Summary method; fast for large datasets
ASTRID [1]	Gene trees	Species tree topology	Uses internode distances; fast
*BEAST [1]	Sequence data	Species tree topology, branch lengths, divergence times	Bayesian; uses sequence data directly
BEST [1]	Sequence data	Species tree topology, branch lengths	Bayesian; computationally intensive
MP-EST [1]	Gene trees	Species tree topology, branch lengths	Uses maximum likelihood on triplets
SVDquartets [1]	Sequence data	Species tree topology	Co-estimates gene trees and species tree
SNAPP [1]	Biallelic markers	Species tree topology, divergence times	Bayesian; uses SNP data without gene trees

These methods operate on different type of inputs—some use pre-estimated gene trees while others co-estimate gene trees and species trees directly from sequence alignments. A significant class of "tuple-based" methods operates by computing summary statistics for subsets of taxa (e.g., quartets or triplets) and then uses these statistics to estimate the species tree [1].

Performance with Missing Data

An important consideration for empirical studies is how these methods perform when gene sequences are missing for some taxa. Research has established statistical consistency for certain coalescent-based methods under models of taxon deletion, including a simple i.i.d. model where every species is missing from every gene with the same probability, and a more general "full subset coverage" model [1]. Simulation studies with methods including ASTRAL-II, ASTRID, MP-EST, and SVDquartets have demonstrated that these methods can maintain accuracy even with substantial amounts of missing data, particularly as the number of genes increases [1]. This robustness is crucial for empirical datasets where practical constraints often result in incomplete taxonomic coverage across genes.

Experimental Protocols and Workflows

Standard Workflow for MSC-based Species Tree Estimation

The following workflow diagram illustrates the standard protocol for estimating species trees under the multispecies coalescent model:

Detailed Methodological Steps

Step 1: Data Collection and Preparation

Collect sequence data from multiple independent loci across the genomes of the target species. The number of loci should be sufficient to resolve species relationships given the expected levels of ILS. For each locus, perform multiple sequence alignment using standard tools (e.g., MAFFT, MUSCLE). Assess alignment quality and trim unreliable regions if necessary.

Step 2: Gene Tree Estimation

Estimate gene trees for each locus using maximum likelihood (e.g., RAxML, IQ-TREE) or Bayesian methods (e.g., MrBayes). For Bayesian approaches, ensure adequate chain convergence and effective sample sizes. It is recommended to assess gene tree uncertainty through bootstrapping or posterior distributions.

Step 3: Species Tree Estimation

Input the estimated gene trees into a coalescent-based species tree method (see Table 1). For summary methods like ASTRAL or MP-EST, this involves providing the set of gene trees. For full-likelihood methods like *BEAST, input the sequence alignments directly. Method selection should consider computational requirements, dataset size, and whether branch lengths are needed.

Step 4: Statistical Assessment

Assess support for the inferred species tree using appropriate methods. For summary approaches, use multi-locus bootstrapping by resampling both loci and sites. For Bayesian methods, examine posterior probabilities. High support values (e.g., ≥95% bootstrap or ≥0.95 posterior probability) provide confidence in the species tree topology.

Protocol for Computing Gene Tree Concordance Probabilities

For researchers interested in calculating the exact probability of gene tree concordance with the species tree under the MSC model, the following protocol outlines the process:

Table 2: Protocol for Computing Gene Tree Concordance Probabilities

Step	Procedure	Tools/Methods	Key Considerations
1. Parameter Specification	Define species tree topology and branch lengths (in coalescent units); specify number of genes sampled per species.	Newick format tree; population size information	Branch lengths critical for accurate probability calculation
2. Probability Calculation	Use dynamic programming algorithm to compute probability of gene tree concordance	Python implementations [5]; polynomial-time algorithms	Exact computation possible for thousands of genes within hours [5]
3. Analysis Extension	Calculate probability of specific concordant gene tree topologies	Customized versions of base algorithm	Useful for comparing expected vs. observed concordance patterns
4. Validation	Compare computed probabilities with empirical gene tree frequencies	Statistical tests (e.g., χ² goodness-of-fit)	Discrepancies may indicate model violations (e.g., gene flow)

Table 3: Essential Research Reagents and Computational Tools for MSC Studies

Category	Item/Resource	Function/Application	Examples/Notes
Sequence Data Types	Ultraconserved Elements (UCEs)	Phylogenomic marker with flanking variable regions	Used in bird [1] and lizard [1] studies
	Transcriptomes	Protein-coding sequence data	Can be problematic due to recombination between exons [3]
	Whole Genome Sequences	Comprehensive genomic data	Provides maximum information but computationally challenging
Software Tools	ASTRAL	Species tree estimation from gene trees	Fast; handles missing data [1]
	*BEAST	Bayesian species tree estimation	Co-estimates gene trees and species tree [1]
	MP-EST	Maximum likelihood species tree estimation	Uses triplets of species [1]
	SVDquartets	Species tree from sequence data	Uses singular value decomposition [1]
Analytical Frameworks	Multilocus MSC	Joint modeling of ILS and gene family evolution	Incorporates duplication, transfer, loss [4]
	MSC with Gene Flow	Species delimitation with migration	Distinguishes continuous variation from species boundaries [6]
Validation Approaches	Contact Zone Analysis	Testing reproductive isolation	Critical for validating species delimitation [6]

Advanced Considerations and Current Research Directions

Integrating Gene Flow and the Multilocus MSC

Recent extensions to the standard MSC model aim to incorporate additional biological complexities. The multilocus multispecies coalescent provides a flexible framework for modeling gene family evolution that accounts for the joint action of ILS and gene duplication, transfer, and loss (DTL) [4]. This integrated approach recognizes that ILS can affect gene copy number polymorphism and interfere with DTL processes, resulting in realized rates of D, T, and L that become non-homogeneous in time when ILS is considered [4]. Similarly, new methods for species delimitation under the MSC with gene flow have been developed to better distinguish continuous geographic variation from actual species boundaries, addressing concerns about over-splitting widespread taxa into multiple species [6]. These approaches provide a range of results depending on assumptions, allowing researchers to explore taxonomic hypotheses more comprehensively.

Quantitative Traits under the MSC

The MSC framework has been extended to model the evolution of quantitative traits, incorporating the effects of genealogical discordance [2]. This approach reveals that discordance decreases the expected trait covariance between more closely related species relative to more distantly related species—a pattern that, if unaccounted for, can lead to overestimation of evolutionary rates, decreased phylogenetic signal, and errors when examining shifts in mean trait values [2]. This effect appears to be largely independent of the number of loci controlling a quantitative trait and also affects discrete threshold traits [2]. These findings have significant implications for phylogenetic comparative methods, suggesting that incorporating genealogical variance is important for accurate inference of trait evolution.

Model Violations and Methodological Debates

The application of MSC models in phylogenomics has prompted discussions about potential model violations and methodological limitations. Critics have pointed to issues such as the use of transcriptome data where recombination between exons may violate MSC assumptions, and concerns about errors in gene tree estimation impacting species tree inference [3]. However, research has shown that MSC models are generally robust to many of these challenges and outperform concatenation approaches in situations with high ILS [3]. It is important to recognize that concatenation can be viewed as a special case of the MSC, which in turn represents a special case of emerging phylogenetic network models that can accommodate both ILS and gene flow [3].

The multispecies coalescent model provides an essential framework for modern phylogenomics, explicitly accounting for the discordance between gene trees and species trees that arises from incomplete lineage sorting. Through continued methodological developments in species tree estimation, extensions to incorporate additional processes like gene flow and quantitative trait evolution, and robust protocols for empirical application, the MSC paradigm offers powerful approaches for reconstructing evolutionary histories across the tree of life. As phylogenomic datasets continue to grow in size and taxonomic scope, the principles and methods outlined in these application notes will remain fundamental to accurate inference of species relationships and character evolution.

Application Notes

Theoretical Foundations: Reconciling Gene Flow with the Multispecies Coalescent Model

The multispecies coalescent (MSC) model provides a fundamental framework for understanding species divergence, traditionally assuming strict divergence where species evolve in isolation after a splitting event [7]. However, empirical genomic evidence consistently challenges this assumption, revealing that gene flow—the transfer of genetic material between populations or species—is not an exception but a pervasive feature of evolution [8]. The MSC model has been extended to account for this reality, incorporating migration parameters that allow for continuous change in gene flow rates over time [7]. This integration is crucial because genealogical discordance caused by gene flow can significantly alter trait covariance structures, leading to overestimation of evolutionary rates and errors in ancestral state reconstruction if unaccounted for [2]. The recognition that species may evolve collectively at major loci through the spread of favourable alleles while simultaneously differentiating at other loci represents a paradigm shift in our understanding of evolutionary dynamics [9].

Genomic Evidence for Pervasive Gene Flow

Recent advances in whole-genome sequencing have demonstrated that adaptive introgression—the natural transfer of beneficial alleles between species through hybridization and backcrossing—occurs across the tree of life, from bacteria to mammals [8]. A comprehensive meta-analysis reveals this process operates across multiple levels of biological organization, from genomic changes to physiological, demographic, and behavioral adaptations [8]. For instance, studies on Pyropia yezoensis seaweed identified seven specific gene flow events between cultivated and wild populations, with introgressed genomic regions (0.3%-25.43% of the genome) exhibiting high genetic diversity, low differentiation, and enrichment for genes involved in stress response and cellular homeostasis [10]. Notably, this gene flow introduced valuable genetic variation without significantly increasing genetic load, demonstrating its potential to enhance adaptive capacity [10].

Table 1: Documentated Cases of Adaptive Introgression Across Taxonomic Groups

Taxonomic Group	Key Finding	Adaptive Benefit	Citation
Mammals	Genomic islands of differentiation co-occur with autosomal introgression	Environmental adaptation & speciation	[8]
Plants (Pyropia yezoensis)	7 gene flow events identified between populations	Stress resistance & thallus development	[10]
Multiple Taxa	Balancing selection maintains introgressed alleles	Preservation of adaptive variation	[8]
Various Species	Strongly advantageous alleles (s=0.11 average for leading QTL) spread rapidly	Rapid adaptation across populations	[9]

Quantitative Implications for Trait Evolution

Gene flow fundamentally alters the genetic architecture of traits and their evolutionary trajectories. Under the multispecies coalescent model with gene flow, genealogical discordance decreases expected trait covariance between closely related species relative to distantly related species [2]. This outcome occurs because trait-affecting substitutions on discordant gene trees increase similarity between distant species while reducing it among close relatives. If unaccounted for, this pattern can lead to overestimation of evolutionary rates, decreased phylogenetic signal, and erroneous inferences about shifts in mean trait values [2]. Importantly, these effects appear independent of the number of loci controlling a quantitative trait and persist even when traits are discretized using threshold models [2]. The pervasiveness of these effects underscores why gene flow must be incorporated into models of trait evolution for accurate parameter estimation.

Evolutionary Consequences and Adaptive Significance

Gene flow serves as a potent evolutionary force that enhances adaptive potential through multiple mechanisms. It introduces novel genetic variation without requiring de novo mutations, thereby providing raw material for natural selection to act upon [10] [11]. This introduced variation can enable evolutionary rescue by helping populations adapt to rapidly changing environments [8]. The adaptive significance of introgressed alleles is particularly pronounced when they confer strong selective advantages—alleles with selection coefficients (s) between 0.05-0.01 can spread across 20 populations in just 4,000-18,000 generations even with gene flow as low as 0.1 migrants per generation [9]. This rapid spread facilitates what has been termed "evolutionary leaps," allowing species to bypass intermediate evolutionary stages and adapt quickly to novel conditions [8]. Paradoxically, while gene flow is traditionally viewed as a homogenizing force, it can also promote divergence and speciation through mechanisms such as transgressive segregation, where extreme phenotypic traits outside the parental range emerge in hybrids [8].

Protocols

Protocol 1: Simulating Gene Trees Under the Multispecies Coalescent with Time-Dependent Migration

Background and Application

This protocol describes a method for simulating gene trees under both the multispecies coalescent and migration, allowing investigation of gene flow effects on species tree inference. The approach enables researchers to test the robustness of phylogenetic methods to violations of the strict divergence assumption and generate realistic test cases for inference programs [7].

Materials and Reagents

Table 2: Essential Research Reagent Solutions for Gene Flow Simulation Studies

Reagent/Software	Function/Application	Specifications
Coalescent Simulator with Migration	Simulates gene trees under MSC with migration	Supports time-dependent migration rates & population sizes [7]
Whole-genome sequencing data	Identifies introgressed regions & estimates gene flow	Multiple individuals per population recommended [10]
⋆BEAST package	Bayesian inference of species trees from multilocus data	Tests impact of migration on species tree posterior distribution [7]
Procrustes analysis	Quantifies similarity between genetic and geographic maps	Enables formal comparison across geographic regions [12]

Procedure

• Parameter Specification:

  • Define population size functions (ν(t)) for each species, representing effective population size scaled by generation length [7]
  • Specify migration rate functions (m(t)) as the fraction of emigrants per time unit for each population pair [7]
  • Set divergence times for speciation events in the species tree

• Model Configuration for Multiple Species:

  • For each population split, configure six migration processes operating in parallel between the three resulting populations (two between each pair) [7]
  • Apply the principle that emigrants from a population are split based on relative sizes of recipient populations after splits [7]
  • Assume the default scenario where splits do not immediately affect migration ability (Figure 4A) [7]

• Simulation Execution:

  • Utilize the continuous-time waiting time algorithm to draw coalescent and migration events [7]
  • Account for the non-homogeneous Poisson process of migration using the rate equation: ra→b(t) = fa→b(t)exp(-∫₀ᵗ fa→b(x)dx) [7]
  • Generate multiple independent gene trees for each simulated scenario

• Downstream Analysis:

  • Infer species trees from simulated gene trees using Bayesian methods (e.g., ⋆BEAST) [7]
  • Compare inferred species trees to true species trees to assess robustness to migration
  • Evaluate detection power for migration events under different parameter spaces

Figure 1: Workflow for simulating gene trees with time-dependent migration under the multispecies coalescent model.

Protocol 2: Detecting and Validating Adaptive Introgression from Genomic Data

Background and Application

This protocol provides a framework for identifying and validating adaptive introgression using whole-genome resequencing data from multiple populations. The approach combines population genetic statistics with functional annotation to distinguish neutral introgression from adaptive gene flow [10] [8].

Materials and Reagents

• Whole-genome resequencing data from multiple individuals (minimum 20× coverage recommended)
• Reference genome for the focal species and potential donor species
• Population genetic software (e.g., for FST, D-statistics, ancestry segmentation)
• Functional annotation databases (e.g., Gene Ontology, KEGG pathways)

Procedure

• Sample Collection and Sequencing:

  • Collect samples from wild and cultivated populations (if applicable) across the species' range [10]
  • Include potential donor species based on phylogenetic proximity and historical contact opportunities
  • Perform whole-genome resequencing using standardized library preparation protocols

• Variant Calling and Filtering:

  • Map reads to reference genome using standard pipelines (BWA-MEM, GATK best practices)
  • Call variants across all populations simultaneously to ensure consistent allele frequency estimation
  • Apply quality filters (mapping quality, base quality, missing data thresholds)

• Introgression Detection:

  • Calculate D-statistics (ABBA-BABA test) to detect significant gene flow between species pairs [8]
  • Use ancestry segmentation approaches to identify genomic regions with exceptional shared ancestry
  • Apply FST and related differentiation metrics to identify regions of unusually low differentiation

• Adaptiveness Assessment:

  • Annotate introgressed regions for gene content, CDS density, and GC content [10]
  • Test for signatures of selection within introgressed regions (e.g., reduced diversity, skewed site frequency spectrum)
  • Perform functional enrichment analysis to identify overrepresented biological processes
  • Correlate introgressed haplotype frequency with environmental variables or phenotypic measurements

• Validation:

  • Use transgenic approaches or gene editing to validate functional consequences of candidate introgressed alleles
  • Perform association studies between introgressed haplotypes and fitness-related traits

Figure 2: Genomic detection and validation workflow for adaptive introgression.

Critical Considerations and Limitations

While gene flow represents a pervasive evolutionary force, several methodological challenges persist in its study. Inference and detection of migration remain problematic even with full likelihood models, particularly when distinguishing between incomplete lineage sorting and genuine gene flow [7]. Furthermore, strongly advantageous alleles can spread so rapidly that they become fixed across populations before sufficient time has passed for neutral differentiation to occur, creating challenges for reconstructing historical gene flow events [9]. The genetic load in cultivated or bottlenecked populations may be significantly higher than in wild populations, though gene flow itself is not necessarily the primary driver of this phenomenon [10]. Future research directions should focus on developing improved statistical methods for distinguishing different forms of gene flow, understanding the interaction between gene flow and other evolutionary forces, and applying this knowledge to conservation challenges in rapidly changing environments.

The Multispecies Coalescent (MSC) model provides a powerful mathematical framework for analyzing genomic sequence data from multiple species, integrating the phylogenetic process of species divergences with the population genetic process of genetic drift [13]. Traditionally, the MSC model assumed strict divergence, where species split at a specific point in time and subsequently evolve in complete isolation [7]. However, empirical genomic studies have increasingly revealed that interspecific gene flow is a major evolutionary force shaping biodiversity [14]. To accommodate this biological reality, two primary extensions to the MSC framework have been developed: the MSC-with-Introgression (MSC-I) model and the MSC-with-Migration (MSC-M) model [14]. These frameworks enable researchers to infer both the species phylogeny and the history of cross-species gene flow from multilocus genomic data, providing crucial insights into evolutionary processes such as hybridization, introgression, and adaptive evolution.

The MSC-I model conceptualizes gene flow as discrete pulse events, representing historical hybridization or introgression episodes at specific time points [14]. In contrast, the MSC-M model conceptualizes gene flow as a continuous process occurring at a constant rate over an extended period [14]. Both models represent idealized simplifications of real-world biological processes, yet they have proven remarkably effective for extracting meaningful information about species divergence and gene flow from genomic datasets [14]. Understanding the assumptions, applications, and limitations of these frameworks is essential for researchers investigating evolutionary history in the presence of gene flow.

Theoretical Foundations of MSC-I and MSC-M Frameworks

The Basic Multispecies Coalescent Model

The multispecies coalescent extends the single-population coalescent model to multiple species, integrating the process of species divergences with the within-population process of genetic drift and mutation [13]. The model incorporates two sets of parameters: species divergence times (τ) and population size parameters (θ) [13]. When tracing genealogical history backward in time, coalescent events occur within populations at rates determined by population sizes, with the process resetting at speciation events due to changes in population size and the introduction of lineages from sibling species [13]. A fundamental feature of the MSC is that gene trees from different loci are independent and must fit within the species tree, meaning sequence divergence must predate species divergence [13]. This intrinsic constraint creates computational challenges but also provides valuable information about ancestral population sizes and gene flow.

MSC-I: Modeling Discrete Introgression Events

The MSC-I model, also known as the multispecies network coalescent (MSNC) or network multispecies coalescent (NMSC) model, conceptualizes gene flow as discrete pulse events occurring at specific time points [14]. In this framework, the amount of gene flow is quantified by the introgression probability (φ), which represents the proportion of migrants from one population to another at the time of introgression [14]. The MSC-I model is particularly useful for modeling historical hybridization events, hybrid speciation, and ancient introgression episodes that occurred within relatively narrow timeframes. This approach has been widely implemented in software packages such as Phylone/mcmc-seq and *BEAST, though computational challenges limit applications to relatively small datasets [14].

MSC-M: Modeling Continuous Migration

The MSC-M model extends the basic MSC framework to incorporate continuous gene flow between species or populations over extended periods [14]. This approach includes the isolation-with-migration (IM) model and its variants, such as the isolation-with-initial-migration (IIM) model and the secondary contact (SC) model [14]. In the MSC-M framework, gene flow is quantified by the population migration rate (M), defined as the expected number of migrant individuals per generation moving from one population to another [14]. Mathematically, M = Nm, where N is the effective population size of the recipient population and m is the proportion of migrants [14]. This model is particularly suitable for investigating ongoing gene flow, parapatric speciation, and cases where species boundaries remain permeable over evolutionary timescales.

Table 1: Key Parameters in MSC Gene Flow Models

Parameter	Model	Definition	Biological Interpretation
φ (phi)	MSC-I	Introgression probability	Proportion of migrants from population A to B at a specific time
M	MSC-M	Population migration rate	Expected number of migrant individuals per generation (M = Nm)
N	Both	Effective population size	Genetic diversity and coalescent rate within a population
τ (tau)	Both	Species divergence time	Time since two species diverged from a common ancestor
m	MSC-M	Migration fraction	Proportion of individuals in a population that are migrants per generation

Comparative Framework: MSC-I vs. MSC-M

Model Assumptions and Mathematical Formulations

The MSC-I and MSC-M frameworks make fundamentally different assumptions about the temporal pattern of gene flow. The MSC-I model assumes that introgression occurs as brief pulses, essentially instantaneous on evolutionary timescales, while the MSC-M model assumes continuous gene flow at a constant rate over extended periods [14]. These differences in temporal scaling lead to distinct mathematical formulations: MSC-I uses discrete transition probabilities at specific time points, whereas MSC-M employs continuous migration rates throughout time intervals [14]. In practice, both models represent endpoints of a continuum, as real-world gene flow likely varies in intensity over time due to changing geographical distributions, ecological conditions, and selective pressures [14].

Performance Under Model Misspecification

A critical consideration for researchers is how these models perform when their underlying assumptions are violated. Recent investigations have examined the impacts of various types of model misspecification, including mis-assignment of gene flow to incorrect lineages, misspecification of the direction of gene flow, and misspecification of the mode of gene flow [14]. Studies indicate that mis-assignment of gene flow to incorrect lineages causes large biases in parameter estimates, while misspecification of the direction of gene flow can make it difficult to distinguish between early divergence with gene flow and recent complete isolation [14]. Interestingly, misspecification of the mode of gene flow (applying MSC-M when MSC-I is appropriate, or vice versa) appears to have relatively small effects on the detection of gene flow, though parameter estimates may be affected [14].

Table 2: Comparative Analysis of MSC-I and MSC-M Frameworks

Characteristic	MSC-I Model	MSC-M Model
Temporal pattern of gene flow	Discrete pulses at specific times	Continuous migration over extended periods
Primary parameters	Introgression probability (φ)	Population migration rate (M)
Best suited for	Ancient hybridization, hybrid speciation, historical introgression	Ongoing gene flow, parapatric speciation, secondary contact
Computational requirements	High (model space exploration challenging)	Moderate to high (additional continuous parameters)
Robustness to mode misspecification	Relatively robust for detecting gene flow	Relatively robust for detecting gene flow
Key limitations	Difficult to explore full model space	Assumes constant migration rates

Experimental Protocols and Methodologies

Protocol 1: Bayesian Inference Under the MSC-I Model

Objective: Implement Bayesian inference to estimate species relationships and introgression probabilities using the MSC-I model.

Materials and Software Requirements:

• Genomic sequence data from multiple loci across several species
• BPP software (for Bayesian phylogenomic analysis under MSC models)
• Sequence alignment tools (e.g., MAFFT, MUSCLE)
• High-performance computing resources

Methodology:

• Data Preparation:

  • Select hundreds to thousands of independent loci from genomic data
  • Ensure loci represent regions with minimal recombination within each locus
  • Align sequences for each locus using appropriate alignment algorithms
  • Verify alignments and check for potential alignment errors

• Model Specification:

  • Define candidate species relationships and potential introgression scenarios
  • Specify prior distributions for parameters (θ, τ, φ)
  • Set prior probabilities for different introgression scenarios

• Markov Chain Monte Carlo (MCMC) Sampling:

  • Run MCMC chains to sample from posterior distribution of parameters
  • Implement reversible-jump MCMC to explore different introgression scenarios
  • Ensure chain convergence through multiple runs with different starting points
  • Assess effective sample sizes for key parameters

• Posterior Analysis:

  • Calculate posterior probabilities for species relationships and introgression events
  • Estimate marginal posterior distributions for parameters (θ, τ, φ)
  • Use Bayes factors to compare support for different gene flow scenarios

Troubleshooting Tips:

• If MCMC convergence is poor, consider adjusting prior distributions or increasing chain length
• For computational efficiency, use analytical integration of gene trees when possible
• Validate results with simulations where the true history is known

Protocol 2: Simulating Gene Trees Under Time-Dependent Migration

Objective: Generate simulated gene trees under complex migration scenarios to test methods or understand model behavior.

Materials and Software Requirements:

• Coalescent simulator with migration capabilities (e.g., ms, COAL)
• Parameter files defining species tree and migration scenarios
• Computing resources for batch simulations

Methodology:

• Parameter Specification:

  • Define species tree topology with divergence times (τ)
  • Specify effective population sizes (θ) for each branch
  • Define migration rates (M) between populations over relevant time intervals
  • For time-dependent migration, specify how rates change over time

• Simulation Configuration:

  • Set sample sizes (number of sequences per species)
  • Determine number of independent loci to simulate
  • Configure output format for generated gene trees

• Execution:

  • Run simulations using specified parameters
  • Generate multiple replicate datasets for uncertainty assessment
  • Output gene trees with branch lengths in coalescent units

• Downstream Analysis:

  • Analyze simulated gene trees to understand properties under migration
  • Use as input for method testing and validation
  • Compare empirical data patterns to simulations

Technical Notes:

• Migration rates are typically specified as the fraction of emigrants per time unit (m) multiplied by the effective population size (M = Nm)
• Time-dependent migration allows modeling scenarios where gene flow decreases gradually after speciation
• The simulation approach can model both MSC-I and MSC-M scenarios through appropriate parameterization

Visualization and Conceptual Diagrams

Workflow for MSC Gene Flow Analysis

MSC-I and MSC-M Model Structures

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Computational Tools for MSC Gene Flow Analysis

Tool/Resource	Type	Primary Function	Application Context
BPP	Software Package	Bayesian inference under MSC	Estimation of species trees, divergence times, and gene flow parameters [14] [15]
⋆BEAST	Software Package	Co-estimation of gene and species trees	Species tree estimation accounting for incomplete lineage sorting [7]
PhyloNet	Software Package	Network phylogeny inference	Modeling reticulate evolution and hybridization [14]
ms/COAL	Simulation Tool	Coalescent simulations with migration	Generating null models and testing method performance [7]
Sequence Capture	Wet-bench Protocol	Targeted genomic sequencing	Generating multi-locus data (UCEs, AHE, RELEC) [14]
RADseq	Wet-bench Protocol	Reduced representation sequencing	Cost-effective genomic sampling for non-model organisms [14]

Applications and Case Studies

Genomic Analysis of Picea (Purple Cone Spruce)

A compelling application of MSC gene flow models comes from genomic analysis of purple cone spruce (Picea spp., Pinaceae), which putatively arose through homoploid hybrid speciation [14]. Researchers applied MSC-based approaches to test this hypothesis and estimate parameters of the hybridization process. By analyzing multilocus genomic data, they demonstrated that the MSC framework could successfully detect historical gene flow and provide support for hybrid origins, even when using simplified models of gene flow [14]. This case study illustrates the practical utility of these methods for investigating complex evolutionary histories in natural systems.

Methodological Insights from Simulation Studies

Simulation studies under the MSC framework with migration have yielded important insights into method performance and limitations. Research has shown that while estimation of species tree topology can be quite robust to the presence of gene flow, the inference and detection of migration is often problematic, even with full-likelihood methods [7]. This highlights the importance of model validation and careful interpretation of results. Studies have further demonstrated that misspecification of the direction of gene flow can significantly impact inferences about the timing of divergence and isolation, potentially leading to confusion between early divergence with gene flow versus recent complete isolation [14].

Future Perspectives and Challenges

The field of MSC-based inference of gene flow faces several important challenges and opportunities for advancement. Computationally efficient methods for exploring the vast space of possible gene flow models remain a significant hurdle, particularly for MSC-I models where the number of possible introgression scenarios grows rapidly with the number of taxa [14]. Future methodological developments will likely focus on improved algorithms for model selection and parameter estimation, as well as approaches for distinguishing between different modes of gene flow. Additionally, integration of MSC frameworks with other evolutionary processes such as selection and recombination represents an important frontier [15]. As genomic datasets continue to grow in size and taxonomic breadth, the development of scalable implementations will be crucial for maintaining the utility of these powerful approaches for inferring evolutionary history in the presence of gene flow.

Gene tree discordance, the phenomenon where gene trees inferred from different genomic regions display conflicting evolutionary histories, presents a major challenge in phylogenomics. This discordance primarily arises from two key biological processes: deep coalescence (often used interchangeably with Incomplete Lineage Sorting, or ILS) and introgression (a form of gene flow). Distinguishing between these processes is crucial for accurate phylogenetic inference, species delimitation, and understanding evolutionary mechanisms [16] [17].

Deep coalescence occurs when ancestral genetic polymorphisms persist through multiple speciation events, causing some gene genealogies to reflect histories that predate the species divergence. In contrast, introgression results from the exchange of genetic material between previously separated lineages after speciation, typically through hybridization. While both processes create conflicting signals across the genome, they stem from fundamentally different evolutionary mechanisms and have distinct implications for understanding speciation and adaptation [18] [19].

The multispecies coalescent (MSC) model provides the fundamental theoretical framework for quantifying these processes, with recent extensions (MSC-I and MSC-M) specifically incorporating models of introgression and gene flow. Accurate discrimination between ILS and introgression has become increasingly important as genomic data reveals the pervasive nature of both processes across the tree of life [14] [20].

The relative contributions of ILS and introgression to gene tree discordance vary across taxonomic groups and evolutionary contexts. Quantitative assessments from empirical studies reveal how these processes interact to shape genomic landscapes.

Table 1: Relative Contributions to Gene Tree Discordance in Empirical Studies

Taxonomic Group	ILS Contribution	Introgression Contribution	Gene Tree Error	Primary Evidence
Fagaceae (Oaks) [21]	9.84%	7.76%	21.19%	Decomposition analysis of nuclear gene trees
Rattlesnakes [20]	Dominant process at rapid radiations	Significant, multiple events	Not quantified	Coalescent-based species trees and networks
Primates (Human/Chimp/Gorilla) [19]	Substantial portion	Ancient gene flow detected	Not quantified	Site patterns and branch length analysis
Loricaria (Asteraceae) [17]	Strong evidence	Strong evidence	Not quantified	D-statistics and phylogenetic networks
Pandanales [22]	Secondary role	Primary source of conflict	Not quantified	Coalescent simulations and QuIBL analysis

The table above demonstrates that both ILS and introgression contribute substantially to phylogenetic discordance, though their relative importance varies. In Fagaceae, gene tree estimation error surprisingly accounts for the largest proportion of discordance, highlighting the importance of accounting for methodological artifacts in phylogenomic analyses [21]. For rapidly radiating groups like rattlesnakes, ILS often dominates due to short internal branches and successive speciation events, though introgression frequently compounds this signal [20].

Table 2: Characteristics of ILS vs. Introgression

Feature	Incomplete Lineage Sorting (ILS)	Introgression
Underlying Process	Random sorting of ancestral polymorphisms	Transfer of alleles between populations/species
Genomic Distribution	Genome-wide, random with respect to function	Often heterogeneous, influenced by selection
Branch Length Signal	Gene trees have longer branches than average [19]	Gene trees have shorter branches than average [19]
Phylogenetic Distribution	Concentrated around short internal branches	Can occur between any related lineages
D-Statistics Signal	Symmetric discordance (ABBA = BABA)	Asymmetric discordance (ABBA ≠ BABA) [17]
Effective Population Size	More pronounced in large populations	Can occur regardless of population size

Methodological Framework for Discrimination

Workflow for Discriminating ILS and Introgression

A robust workflow for distinguishing between ILS and introgression integrates multiple lines of evidence and validation steps, progressing from initial data processing to comprehensive analysis.

Key Methodological Approaches

Summary statistic approaches provide computationally efficient methods for initial detection of introgression. The D-statistic (ABBA-BABA) tests for asymmetry in site patterns that deviate from expectations under pure ILS [19] [17]. For a four-taxon configuration (((P1,P2),P3),O), the test examines the relative frequencies of ABBA and BABA patterns, where significant deviations from equality suggest introgression between P2 and P3 (ABBA > BABA) or P1 and P3 (BABA > ABBA). Extensions like the f4-statistic and HyDe provide additional power to detect and localize introgression events [17].

Protocol: D-Statistic Implementation

• Dataset Preparation: Generate whole-genome alignments or multiple sequence alignments for numerous loci
• Taxon Sampling: Select four-taxon blocks following the relationship (((P1,P2),P3),Outgroup)
• Site Pattern Counting:
  • Identify informative sites where Outgroup and P3 share the same base
  • Count ABBA patterns (P1=A, P2=B, P3=B, Outgroup=A)
  • Count BABA patterns (P1=B, P2=A, P3=B, Outgroup=A)
• Statistical Testing:
  • Calculate D = (ABBA - BABA) / (ABBA + BABA)
  • Assess significance using block jackknife or parametric bootstrap
  • Values significantly different from zero indicate introgression

Probabilistic Modeling Approaches

Full-likelihood methods under the multispecies coalescent framework provide the most powerful approach for distinguishing ILS from introgression. The MSC-with-introgression (MSC-I) model incorporates discrete hybridization events with introgression probabilities (φ), while the MSC-with-migration (MSC-M) model assumes continuous gene flow at constant rates [14]. These models can be implemented in software such as BPP and STACEY, which use Bayesian MCMC algorithms to jointly estimate species trees, divergence times, and gene flow parameters.

Protocol: MSC-I Analysis with BPP

• Model Specification:
  • Define candidate species tree or network topology
  • Specify prior distributions for divergence times (τ) and population sizes (θ)
  • Set priors for introgression probabilities (φ) between specific lineages
• MCMC Configuration:
  • Run multiple independent chains with different starting values
  • Use fine-tuning to ensure acceptable acceptance rates (20-40%)
  • Conduct convergence diagnostics (ESS > 200, PSRF ≈ 1.0)
• Model Comparison:
  • Calculate marginal likelihoods using stepping-stone sampling
  • Compare models with and without introgression using Bayes factors
  • BF > 10 provides strong evidence for including introgression parameters

Supervised Learning Approaches

Emerging machine learning approaches, particularly semantic segmentation, frame the detection of introgressed loci as a classification problem [16]. These methods leverage patterns in topological features, branch lengths, and site characteristics to distinguish between ILS and introgression signals across genomic alignments.

Protocol: Feature Extraction for Classification

• Per-Locus Calculations:
  • Calculate gene tree topologies and bootstrap support
  • Measure internal branch lengths relative to species tree
  • Compute partition-specific substitution rates
• Window-Based Statistics:
  • Slide windows across aligned genomes (e.g., 10-50kb)
  • Calculate D-statistics, f4-statistics, and fd within each window
  • Measure divergence (dXY) and diversity (π) patterns
• Training Data Preparation:
  • Use simulated data with known ILS/introgression events
  • Generate examples under various demographic scenarios
  • Balance training sets to avoid classification bias

Table 3: Computational Tools for Discriminating ILS and Introgression

Tool/Resource	Primary Function	Methodological Approach	Key Application
BPP [14]	Bayesian phylogenomic analysis	MSC-with-introgression (MSC-I)	Joint estimation of species trees and gene flow parameters
ASTRAL [18]	Species tree estimation	Multi-species coalescent	Robust species tree inference under ILS
HyDe [17]	Hybridization detection	Site pattern analysis	Genome-wide detection of hybrid speciation and introgression
PhyloNet [20]	Phylogenetic network inference	Multi-species coalescent networks	Visualization and inference of reticulate evolution
Aphid [19]	Source quantification	Approximate likelihood	Distinguishing ILS vs. introgression using branch lengths
IQ-TREE [21]	Gene tree inference	Maximum likelihood	Efficient estimation of individual gene trees
Dsuite [20]	Introgression analysis	D-statistics/f4-statistics	Comprehensive D-statistic calculations across genomes

Table 4: Genomic Data Types for Discordance Analysis

Data Type	Resolution	Advantages	Limitations
Whole Genome Sequencing	Nucleotide	Maximum resolution, complete information	Computational burden, assembly challenges
RNA-Seq/Transcriptomes [18] [22]	Coding regions	Cost-effective, targets functional elements	Limited to expressed genes, tissue-specific
RAD-seq/SLAF-seq [23]	Reduced representation	Genome-wide sampling without reference	Locus dropout, homology assessment challenges
UCEs/AHE [14]	Conserved regions	Cross-species applicability, orthology certainty	Limited to conserved regions, potential selection bias
Target Capture [17]	Targeted regions	Customizable, consistent locus sampling	Design effort, limited to targeted regions

Case Study Applications

Ancient Radiation: Amaranthaceae Phylogenomics

Research on Amaranthaceae s.l. exemplifies the challenges of disentangling ILS and introgression in ancient radiations. The study integrated 88 transcriptomes and 7 genomes to test hypotheses of ancient hybridization [18]. Researchers employed coalescent-based species trees, network inference, synteny analyses, and simulations to address gene tree discordance. They discovered that three consecutive short internal branches produced anomalous gene trees, with the combined effects of ILS and potentially ancient introgression contributing to the high discordance. This case highlights the importance of multi-method approaches and the challenges of identifiability when distinguishing sources of conflict in ancient radiations [18].

Recent Radiation: Rattlesnake Phylogenomics

The rattlesnake (Crotalus and Sistrurus) system demonstrates how rapid diversification coupled with introgression creates complex phylogenetic discordance [20]. Genomic data from nearly all species revealed that rapid speciation resulted in individual gene trees conflicting with the species tree, while incomplete speciation and frequent hybridization further complicated phylogenetic inference. The study employed MSC network approaches to jointly model ILS and introgression, revealing that both processes have significantly influenced previous phylogenetic interpretations. This system exemplifies a "network radiation" where evolutionary relationships can only be accurately understood using genome-wide data and network-based methods [20].

Plant Systematics: Fagaceae and Pandanales

In Fagaceae, decomposition analysis quantified that gene tree estimation error (21.19%), ILS (9.84%), and gene flow (7.76%) all contributed significantly to phylogenetic discordance [21]. Researchers identified "consistent" and "inconsistent" genes based on phylogenetic signals, finding that excluding inconsistent genes reduced conflicts between concatenation and coalescent approaches.

For Pandanales, phylogenomic analysis of transcriptomic data identified ancient gene flow as the primary source of conflict, with two significant events detected: between Velloziaceae and Triuridaceae, and between Triuridaceae and the C-P clade [22]. The study demonstrated that gene flow, rather than ILS, was the dominant factor in phylogenetic discordance, while also identifying multiple whole-genome duplication events that further complicated evolutionary history.

Best Practices and Implementation Guidelines

Experimental Design Considerations

Effective discrimination between ILS and introgression requires careful experimental design:

• Taxon Sampling: Include multiple individuals per species to account within-population variation and improve parameter estimation [23]
• Genomic Coverage: Balance between number of loci and sequence length; numerous shorter loci better capture ILS while longer sequences improve gene tree accuracy [14]
• Outgroup Selection: Choose appropriate outgroups for rooting and polarization of site patterns in D-statistic analyses [19]
• Genome Partitioning: Consider genomic landscapes by analyzing linked and unlinked regions separately to account for variation in recombination rates [20]

Validation and Convergence Assessment

Robust inference requires thorough validation:

• Model Comparison: Always compare models with and without gene flow using appropriate statistical tests (e.g., Bayes factors, AIC) [14]
• Convergence Diagnostics: For Bayesian methods, ensure MCMC chains have converged using ESS > 200 and PSRF ≈ 1.0 [14]
• Simulation Studies: Use parametric simulations to verify method performance under known evolutionary scenarios [18]
• Multiple Methods: Triangulate findings using independent approaches (e.g., summary statistics + full likelihood methods) [17]

Interpretation Caveats

Several important caveats should guide interpretation:

• Model Misspecification: MSC models assume no recombination within loci and free recombination between loci; violations can bias parameter estimates [14]
• Ghost Introgression: Gene flow from unsampled or extinct lineages can create patterns resembling ILS [16] [14]
• Anomaly Zone Effects: In rapid radiations, the most common gene tree may not match the species tree even under pure ILS [20]
• Selection Effects: Regions under selection can exhibit introgression patterns that differ from neutral expectations [20]

Future Directions

Emerging methods are extending the framework for distinguishing ILS and introgression. Supervised learning approaches, particularly when detection of introgressed loci is framed as a semantic segmentation task, show great potential for automating genome-wide scans [16]. Improved MSC-with-gene-flow models that incorporate more realistic demographic scenarios are addressing limitations of current simplified models (MSC-I and MSC-M) [14] [6]. Additionally, methods for detecting ghost introgression from unsampled lineages and approaches that jointly model selection and introgression are active areas of development [16].

The integration of phylogeographic models with phylogenetic network inference represents another promising direction, allowing researchers to simultaneously reconstruct patterns of gene flow and geographical range evolution. As these methods mature, they will enhance our ability to distinguish between deep coalescence and introgression across diverse evolutionary scenarios.

The Multispecies Coalescent (MSC) model provides a powerful mathematical framework for analyzing genomic sequence data to infer evolutionary histories. This stochastic process describes the genealogical relationships of DNA sequences sampled from different species, explicitly modeling how gene trees can differ from the species tree due to ancestral genetic polymorphism—a phenomenon known as Incomplete Lineage Sorting (ILS) [24]. The basic MSC model assumes no gene flow after species divergence, but extensions now incorporate both continuous migration (MSC-M) and discrete introgression (MSC-I), allowing researchers to quantify interspecific gene flow [14]. These models have become essential in evolutionary biology, enabling estimation of key parameters including species divergence times, population sizes, and introgression probabilities from genomic data [16] [25].

Accurate estimation of these parameters provides crucial insights into fundamental evolutionary processes. For example, introgression has been shown to facilitate adaptation in various taxa, with introgressed loci often linked to immunity, reproduction, and environmental adaptation [16] [25]. The MSC framework has been successfully applied to study gene flow in diverse organisms, including Heliconius butterflies [25], hominins [16], and spruce species [14], demonstrating its broad utility.

Core Parameters in the MSC Model with Gene Flow

Theoretical Foundation and Mathematical Formulation

The MSC model traces gene genealogies backward in time through the species tree. For each population, the coalescent process describes how lineages merge until they reach the population's root. The probability distribution of gene trees and coalescent times forms the basis for parameter estimation [24]. The coalescence rate for j lineages is determined by the population size parameter θ = 4Nμ, where N is the effective population size and μ is the mutation rate per generation [24]. When extended to incorporate gene flow, the model includes additional parameters for introgression probability (φ) or migration rate(M), which quantify the strength and direction of interspecific gene flow [14].

The following table summarizes the key parameters estimated under the MSC model with gene flow:

Table 1: Key Parameters in the MSC Model with Gene Flow

Parameter	Symbol	Definition	Biological Interpretation
Divergence Time	τ (tau)	Time in generations since two lineages diverged	Speciation events or population splits
Population Size	θ (theta) = 4Nμ	Population size parameter	Effective population size (N) scaled by mutation rate (μ)
Introgression Probability	φ (phi)	Proportion of genomes that migrated between species	Strength of historical hybridization
Population Migration Rate	M = Nm	Number of migrant individuals per generation	Rate of continuous gene flow between populations
Base-Calling Error Rate	e	Probability of incorrect base assignment during sequencing	Data quality indicator affecting parameter estimation

Parameter Estimation and Biological Interpretation

Divergence times (τ) represent speciation events or population splits and are typically measured in generations, though they can be converted to years with mutation rate calibration. These parameters establish the temporal framework for the species tree [24]. Population sizes (θ) reflect genetic diversity levels, with larger values indicating more diverse ancestral populations. Both current and ancestral population sizes can be estimated, providing insights into historical demographic changes [24] [25].

The introgression probability (φ) in the MSC-I model represents the proportion of the genome that moved from one species to another during discrete hybridization events [14]. In contrast, the population migration rate (M) in the MSC-M model quantifies continuous gene flow as the expected number of migrant individuals per generation [14]. Understanding these parameters helps evolutionary biologists reconstruct historical gene flow patterns and their role in adaptation and speciation.

Data Requirements and Experimental Design

Implementing the MSC model with gene flow requires multi-locus sequence data from multiple individuals across the studied species. Data can be generated through various methods: (1) sampling short genomic fragments from whole-genome sequencing, or (2) targeted sequence capture approaches such as RADseq, UCEs (Ultraconserved Elements), AHE (Anchored Hybrid Enrichment), or exome sequencing [14]. These genomic fragments are treated as independent loci, with the critical assumption of no recombination within loci but free recombination between them [14].

Data quality significantly impacts parameter estimation accuracy. Recent research demonstrates that sequencing depth profoundly affects inference reliability. At low base-calling error rates (e = 0.001, Phred score 30), species trees and population parameters remain largely unaffected even at low depths (~3×). However, high error rates (e = 0.005-0.01) combined with low sequencing depths (<10×) can reduce power and introduce substantial biases in estimates of population sizes, divergence times, and gene flow rates [26]. To minimize these issues, the simulation study recommends sequencing fewer samples at higher depths rather than many samples at low depths [26].

Experimental Workflow for Parameter Estimation

The following diagram illustrates the complete workflow for estimating key parameters under the MSC model with gene flow:

Research Reagent Solutions and Computational Tools

Table 2: Essential Research Reagents and Computational Tools for MSC Analysis

Tool/Resource	Type	Primary Function	Application Context
Bpp	Software Package	Bayesian inference of species tree & gene flow parameters	Estimation of τ, θ, φ under MSC-I and MSC-M models [26] [14]
Phylone/mcmc-seq	Software Package	MCMC updating of MSC-I model	Inference of phylogenetic networks with introgression [14]
*BEAST	Software Package	Co-estimation of gene trees and species trees	Multispecies coalescent analysis [14]
δaδi	Software Package	Inference of demographic history	Joint site frequency spectrum analysis [14]
Fastsimcoal2	Software Package	Coalescent simulations	Demographic modeling with gene flow [14]
Whole-genome sequencing data	Data Type	Provides comprehensive genomic coverage	Ideal for MSC analysis across diverse taxa [25]
Targeted sequence capture (UCEs, AHE)	Data Type	Cost-effective phylogenomic data	Reduced representation approaches [14]
High-performance computing	Infrastructure	Computational resource intensive analyses	Essential for Bayesian MCMC methods [26]

Detailed Protocol for Parameter Estimation

Data Preparation and Quality Control

Step 1: Locus Selection and Alignment

• Select short genomic fragments (loci) far apart in the genome to ensure independent genealogies [14]
• Verify minimal recombination within each locus while maximizing inter-locus recombination
• Align sequences for each locus using standard tools (e.g., MAFFT, MUSCLE)

Step 2: Data Quality Assessment

• Calculate per-site read depths and identify low-coverage regions
• Estimate base-calling error rates from quality scores
• Apply filters to remove or mask low-confidence regions, particularly with depth <10× and error rate >0.005 [26]
• Consider treating heterozygotes as missing data at low sequencing depths to reduce genotyping errors [26]

Step 3: Data Formatting for MSC Analysis

• Convert alignments to format compatible with chosen software (e.g., Bpp)
• Define species/population assignments for each sample
• Prepare control file specifying model and parameter priors

Model Specification and Prior Selection

Step 4: Choosing Between Gene Flow Models

• Select MSC-I (discrete introgression) for pulse admixture events [14]
• Select MSC-M (continuous migration) for ongoing gene flow [14]
• Specify possible introgression directions based on biological knowledge

Step 5: Setting Parameter Priors

• Use gamma or inverse gamma priors for population sizes (θ)
• Use gamma priors for divergence times (τ)
• Use beta priors for introgression probabilities (φ)
• Consider using biologically informed priors when available

The following diagram illustrates the relationship between key parameters in the MSC model with gene flow:

Execution and Convergence Assessment

Step 6: Running MCMC Analysis

• Execute Bayesian inference using appropriate software (e.g., Bpp)
• Run multiple independent chains from different starting points
• Ensure adequate MCMC chain length (typically millions of generations)

Step 7: Assessing Convergence

• Monitor parameter trace plots for stationarity
• Calculate effective sample sizes (ESS > 200 for key parameters)
• Verify convergence between independent runs using Gelman-Rubin statistics

Step 8: Interpreting Results

• Extract posterior distributions for key parameters (τ, θ, φ/M)
• Calculate Bayesian credibility intervals (typically 95%)
• Identify significantly supported introgression events (posterior probability > 0.95)

Applications and Case Studies

Heliconius Butterflies: A Model System for Introgression

The full-likelihood MSC approach has been successfully applied to reconstruct the evolutionary history of Heliconius butterflies, a group known for extensive hybridization and adaptive wing patterning. Analysis of whole-genome data using Bpp provided robust estimates of species divergence times, ancestral population sizes, and the direction, timing, and intensity of gene flow [25]. The study revealed that Heliconius aoede most likely represents the earliest-branching lineage, contrary to previous hypotheses, and demonstrated that 'silvaniform' species are paraphyletic within the melpomene-silvaniform group [25].

The analysis incorporated chromosome-level estimates of parameters, revealing different phylogenetic signals between autosomes and the Z chromosome. This approach provided new, parsimonious histories for the origins of key traits in Heliconius, including pollen feeding and an inversion involved in wing pattern mimicry [25]. The study exemplifies how MSC methods can resolve complex evolutionary histories despite extensive gene flow.

Impact of Model Misspecification and Data Quality

Research has demonstrated that MSC models are remarkably robust to certain misspecifications. For instance, misspecifying the mode of gene flow (MSC-I vs. MSC-M) has relatively small effects, and gene flow is detected with high power despite this misspecification [14]. However, incorrectly assigning gene flow to the wrong lineages can cause large biases in estimated rates [14]. Similarly, high sequencing error rates (e > 0.005) combined with low read depths (<10×) can substantially reduce estimation accuracy for all parameters [26].

Table 3: Effects of Data Quality and Model Misspecification on Parameter Estimation

Factor	Impact on Parameter Estimation	Recommended Mitigation Strategy
Low sequencing depth (<10×)	Biases in θ, τ, and gene flow estimates [26]	Sequence at higher depth (>10×); treat heterozygotes as missing data [26]
High error rate (>0.005)	Reduced species tree estimation power; biased parameter estimates [26]	Implement rigorous quality filtering; use higher sequencing depth [26]
Misspecified gene flow direction	Difficulty distinguishing early divergence with gene flow from recent complete isolation [14]	Compare multiple gene flow models; use biological knowledge to inform directions
Incorrect lineage assignment	Large biases in estimated gene flow rates [14]	Test alternative lineage assignments; use model comparison approaches
Misspecified mode (MSC-I vs MSC-M)	Relatively small local effects on parameter estimates [14]	Both models effectively detect gene flow despite misspecification

Troubleshooting and Optimization Strategies

Addressing Common Implementation Challenges

Problem: Poor MCMC Convergence

• Solution: Increase chain length and thinning intervals
• Solution: Adjust proposal mechanisms to improve acceptance rates
• Solution: Use parameter transformations to improve mixing

Problem: Unidentifiable Parameters

• Solution: Simplify the model by reducing parameters
• Solution: Incorporate additional biological constraints as priors
• Solution: Use model comparison approaches to select identifiable models

Problem: Computational Limitations

• Solution: Utilize subset analyses focusing on key species relationships
• Solution: Employ high-performance computing resources
• Solution: Use analytical approximations for initial exploratory analyses

Optimizing Parameter Estimation Accuracy

Strategy: Data Quality Enhancement

• Aim for sequencing depths >10×, especially when error rates are high [26]
• Implement rigorous variant calling and filtering protocols
• Validate genotype calls using multiple approaches

Strategy: Model Checking

• Perform posterior predictive simulations to assess model fit
• Compare results across different software implementations
• Conduct sensitivity analyses on prior distributions

Strategy: Biological Validation

• Compare parameter estimates with independent evidence (e.g., fossil calibrations)
• Test predictions based on parameter estimates against known biology
• Replicate analyses with different taxonomic sampling schemes

The protocols and applications outlined here demonstrate the power of the MSC model with gene flow for estimating key evolutionary parameters from genomic data. When implemented with attention to data quality, model selection, and validation procedures, these methods provide robust insights into species divergence times, historical population sizes, and patterns of introgression that have shaped biodiversity.

Practical Implementation: Bayesian Methods, Software, and Species Delimitation

Bayesian MCMC Inference under the MSC with Gene Flow

The Multispecies Coalescent (MSC) model provides a powerful mathematical framework for inferring species relationships and demographic history from genomic data. By incorporating gene flow into this model, researchers can more accurately reconstruct evolutionary histories where hybridization and introgression have occurred. The integration of Bayesian Markov Chain Monte Carlo (MCMC) methods enables robust statistical inference under these complex models, allowing researchers to estimate key parameters such as species divergence times, population sizes, and migration rates from genome-scale datasets [27].

Two primary models have been developed for inference under the MSC with gene flow: the MSC-with-Introgression (MSC-I) model, which conceptualizes gene flow as discrete pulses or hybridization events, and the MSC-with-Migration (MSC-M) model, which treats gene flow as a continuous process occurring at a constant rate over an extended period [14]. Bayesian MCMC approaches provide a flexible framework for estimating parameters under both models, even when the true mode of gene flow differs from the assumed model, though careful interpretation is required [14].

Model Framework and Theoretical Background

The Multispecies Coalescent with Gene Flow

The multispecies coalescent models the genealogical history of genes within a population genetic framework. It extends the single-population coalescent to multiple species or populations connected by a phylogenetic tree, providing a natural model for analyzing multilocus sequence data from closely related species.

Table 1: Key Parameters in MSC Models with Gene Flow

Parameter	Description	Model
φ	Introgression probability (proportion of migrants)	MSC-I
M	Population migration rate (expected number of migrants per generation)	MSC-M
τ	Species divergence time (scaled by mutation rate)	MSC-I & MSC-M
θ	Population size (4Nμ, where N is effective population size and μ is mutation rate)	MSC-I & MSC-M

When augmented with gene flow, the MSC model can account for the discordance between gene trees and species trees that results from both incomplete lineage sorting (ILS) and interspecific gene flow [3]. The MSC-I model, also known as the multispecies network coalescent (MSNC), incorporates gene flow as discrete events at specific time points, measured by the introgression probability (φ) representing the proportion of migrants from one population to another at the time of introgression [14]. The MSC-M model, which includes the isolation-with-migration (IM) model as a special case, assumes continuous gene flow at a constant rate over time, quantified by the population migration rate (M = Nm), where N is the effective population size and m is the migration rate per generation [14].

Bayesian MCMC Inference

Bayesian inference provides a coherent framework for estimating the posterior distribution of parameters in MSC models with gene flow. MCMC algorithms enable sampling from this posterior distribution even for complex models where analytical solutions are intractable [28]. The general form of the posterior distribution can be represented as:

P(Θ | D) ∝ P(D | Θ) × P(Θ)

Where Θ represents the model parameters (divergence times, population sizes, migration rates), D represents the sequence data, P(D | Θ) is the likelihood of the data given the parameters, and P(Θ) is the prior distribution of the parameters [28].

The Metropolis algorithm, a foundational MCMC method, works by generating a Markov chain that explores the parameter space through a series of proposed moves that are either accepted or rejected based on the ratio of posterior probabilities [28]. This approach allows researchers to approximate the posterior distribution without needing to compute the complex normalizing constant, making it particularly valuable for MSC models with gene flow where the likelihood function is computationally intensive to evaluate [27].

Computational Implementation

Software and Algorithms

Recent advances have implemented MSC models with gene flow in computationally efficient Bayesian software, making genome-scale analyses feasible:

Table 2: Research Reagent Solutions for Bayesian MSC Inference

Tool/Resource	Type	Primary Function	Key Features
BPP	Software Package	Bayesian phylogenetics	Implements both MSC-I and MSC-M models; efficient MCMC algorithms for genome-scale data [27]
mstree	Analytical Method	Parameter estimation without isolation-migration assumptions	Uses mathematical inequalities to estimate ancestral population sizes and divergence times under different gene flow modes [29]
3s	Software Package	Species tree estimation	Implements MSC model for species tree inference [29]
IMa3	Software Package	Isolation-with-Migration analysis	Estimates migration rates, divergence times, and population sizes under IM model [29]

The program BPP has implemented both MSC-I and MSC-M models with efficient MCMC algorithms that enable analysis of datasets with thousands of loci [27]. These implementations use sophisticated MCMC algorithms that allow efficient sampling from the posterior distribution, making it possible to test whether gene flow occurred continuously over time or in discrete pulses [27].

Alternative approaches like mstree use mathematical inequalities among species divergence times, ancestral population sizes, and gene tree counts to estimate parameters without assuming a specific relationship between isolation and migration, providing robustness to different modes of gene flow [29].

Workflow and Analytical Procedure

The following diagram illustrates the general workflow for Bayesian MCMC analysis under the MSC with gene flow:

Application Notes and Protocols

Experimental Design Considerations

Data Requirements and Locus Selection

For reliable inference under the MSC with gene flow, genomic data should consist of multiple independent loci sampled from across the genome. The MSC model assumes no recombination within loci and free recombination between loci, guiding appropriate locus selection [14]. Two primary strategies for generating multilocus data are:

• Sampling short genomic fragments from sequenced genomes
• Targeted sequence capture approaches including RADseq, exome sequencing, UCEs (ultraconserved elements), AHE (anchored hybrid enrichment), and CNEEs (conserved nonexonic elements) [14]

For Bayesian MCMC inference, dataset sizes typically range from hundreds to thousands of loci, with longer sequences providing more information for accurate parameter estimation [27].

Model Selection and Specification

Choosing between MSC-I and MSC-M models depends on biological expectations and the research question:

• MSC-I (discrete introgression) is appropriate for modeling hybridization events, historical admixture, or hybrid speciation
• MSC-M (continuous migration) better represents situations with ongoing gene flow between partially isolated populations

Despite their differences, both models can detect gene flow even when misspecified, though parameter estimates may be affected [14]. The implementation in BPP allows comparison between models to assess which better fits the data [27].

Detailed Protocol for Bayesian Analysis

Data Preparation and Alignment

• Process raw sequence data through quality control, filtering, and alignment for each locus
• Define orthologous loci ensuring they represent independent genealogical histories
• Format sequence data according to the requirements of the chosen analysis software (e.g., BPP, mstree, IMa3)
• Generate initial gene trees for each locus if required by the analysis pipeline

Model Configuration and Prior Specification

• Define the species tree topology based on prior knowledge or preliminary analyses
• Specify gene flow events including the direction and lineages involved
• Set prior distributions for key parameters:
  • Divergence times (τ): Use informative priors based on fossil evidence or previous studies when available
  • Population sizes (θ): Typically assigned diffuse priors if little prior information exists
  • Migration parameters (M or φ): Use weakly informative priors centered on zero for conservative testing of gene flow

MCMC Execution and Diagnostics

• Run preliminary MCMC chains to tune proposal distributions and optimize sampling efficiency
• Execute multiple independent MCMC chains with different starting values to assess convergence
• Monitor convergence using diagnostic statistics such as Gelman-Rubin statistics (potential scale reduction factors) and effective sample sizes (ESS)
• Ensure adequate sampling by running chains until ESS > 200 for parameters of interest

Posterior Analysis and Interpretation

• Summarize posterior distributions of parameters (means, medians, credible intervals)
• Assess evidence for gene flow by examining posterior distributions of migration rates or introgression probabilities
• Compare alternative models using Bayes factors or other model comparison techniques when multiple gene flow scenarios are plausible
• Visualize results through species trees with annotated gene flow events and posterior distributions of key parameters

Performance and Robustness

Model Misspecification and Robustness

Studies evaluating the performance of MSC models with gene flow have revealed several important patterns:

Table 3: Effects of Model Misspecification on Inference

Misspecification Type	Impact on Inference	Recommendations
Incorrect lineage assignment	Large biases in estimated migration rates [14]	Use independent evidence (e.g., geographic, ecological) to inform lineage assignment
Direction of gene flow misspecified	Difficulty distinguishing early divergence with gene flow from recent complete isolation [14]	Compare multiple directional hypotheses when uncertain
Mode of gene flow misspecified (MSC-I vs MSC-M)	Relatively small local effects; gene flow still detected with high power [14]	Consider both pulse and continuous models for comprehensive analysis
Ghost introgression (unsampled lineages)	Can severely mislead inference if not accounted for [14]	Include outgroups and consider sampling design carefully

Bayesian tests generally show high power for detecting both recent and ancient gene flow, between both sister and nonsister lineages, though summary methods like the D-statistic often cannot identify the direction of gene flow [14].

Comparison with Alternative Methods

The mstree method provides an alternative approach that does not rely on assumptions about the relationship between isolation and migration, instead using mathematical inequalities among species divergence times, ancestral population sizes, and gene tree counts to estimate parameters [29]. Simulation studies show that mstree can accurately estimate ancestral population sizes and speciation times in the presence of different modes of gene flow, performing comparably to model-based methods like 3s and IMa3 while being computationally faster [29].

Full-likelihood methods based on the MSC generally outperform summary methods (e.g., D-statistic, Hyde, Snaq) that rely on genome-wide averages like site pattern counts or estimated gene tree topologies, as these approaches do not fully utilize information in the data and often have reduced power to detect gene flow [14].

Advanced Applications and Case Studies

Analysis of Purple Cone Spruce (Picea spp.)

The purple cone spruce (Picea spp.) provides an empirical example where MSC models with gene flow have been successfully applied. Genomic data analysis supported the hypothesis that this group arose through homoploid hybrid speciation, demonstrating the practical utility of these methods for identifying hybrid origins in natural systems [14].

Mosquito Genomes (Anopheles)

Analyses of genomic data from Anopheles mosquitoes using Bayesian MSC methods have revealed rich information about the mode and rate of gene flow, contributing to our understanding of how gene flow impacts speciation and adaptation in disease vectors [27].

Bayesian MCMC inference under the multispecies coalescent with gene flow represents a powerful framework for investigating evolutionary history using genomic data. The implementation of efficient MCMC algorithms in software like BPP has made it feasible to analyze genome-scale datasets with thousands of loci, enabling robust detection and quantification of gene flow [27]. While model misspecification remains a concern, particularly regarding the assignment of gene flow to specific lineages, these methods have proven effective for extracting meaningful biological information about species divergence and gene flow from genomic data [14].

As genomic datasets continue to grow in size and taxonomic breadth, Bayesian methods under the MSC will play an increasingly important role in unraveling the complex histories of gene flow that have shaped biodiversity. Future methodological developments will likely focus on improving computational efficiency, incorporating more complex models of gene flow, and better integrating geographic and ecological information to inform hypotheses about gene flow.

BPP is a Bayesian Markov Chain Monte Carlo (MCMC) program for analyzing multi-locus genomic sequence data under the multispecies coalescent (MSC) model [30] [31]. The MSC model provides a powerful mathematical framework that extends the single-population coalescent theory to multiple species, integrating both the phylogenetic process of species divergences and the population genetic process of gene lineage coalescence [15]. This integration enables researchers to address fundamental questions in evolutionary biology, particularly in the context of species divergence despite ongoing gene flow.

The multispecies coalescent with gene flow represents a significant extension to the basic MSC model, incorporating the reality that many closely related species continue to exchange genetic material after initial divergence [15]. Within this research context, BPP serves as an essential toolkit for investigating historical divergence times, population sizes, species boundaries, and the intensity of cross-species gene flow using genomic sequence data. The program has become particularly valuable in modern phylogenomics where researchers regularly encounter incongruence between gene trees and the species tree due to incomplete lineage sorting and potential hybridization events.

Analytical Capabilities of BPP

BPP supports four primary types of analyses that are fundamental to species tree estimation and delimitation research, each addressing specific questions in evolutionary biology [30]:

• Parameter estimation (A00 analysis): Estimation of species divergence times and population size parameters under the multispecies coalescent model on a fixed species phylogeny.

• Species tree estimation (A01 analysis): Inference of the species tree when species assignment and delimitation are already known and fixed.

• Species delimitation (A10 analysis): Determination of species boundaries using a fixed guide tree to explore different delimitation scenarios.

• Joint species delimitation and species-tree estimation (A11 analysis): Simultaneous inference of both species boundaries and phylogenetic relationships without a fixed guide tree, also known as unguided species delimitation.

For the joint analysis (A11), BPP implements novel prior distributions that assign uniform probabilities for different numbers of delimited species, which is particularly important when assignment, species delimitation, and species phylogeny are all inferred simultaneously [30].

Table 1: Analysis Types in BPP Software

Analysis Code	Analysis Type	Fixed Components	Inferred Components
A00	Parameter Estimation	Species phylogeny, delimitation, assignment	Divergence times, population sizes
A01	Species Tree Estimation	Species delimitation, assignment	Species phylogeny
A10	Species Delimitation	Guide tree, assignment	Species boundaries
A11	Joint Estimation	Assignment only	Species delimitation and phylogeny

Practical Implementation and Workflow

Data Requirements and Preparation

BPP requires multi-locus sequence data from multiple species, typically consisting of genomic regions that are unlinked or sufficiently distant to represent independent genealogical histories [31]. The program accepts sequence alignments in various formats and can handle missing data, making it suitable for the heterogeneous datasets commonly encountered in modern phylogenomic studies. For analyses involving gene flow, the data requirements are more stringent, needing sufficient phylogenetic signal to distinguish between incomplete lineage sorting and genuine introgression events [15].

Table 2: Data Specifications for BPP Analysis

Data Component	Specification	Importance for Analysis
Sequence Type	DNA sequence data (nuclear preferred)	Provides genealogical information for coalescent analysis
Number of Loci	Multiple unlinked loci (typically > 10)	Enables estimation of population parameters and species trees
Sequence Length	Sufficient for phylogenetic signal	Affects precision of parameter estimates
Taxon Sampling	Multiple individuals per species	Crucial for estimating population sizes and delimitation
Data Format	FASTA, PHYLIP, or NEXUS	Compatibility with BPP input requirements

Computational Protocols

Running BPP involves several key steps, beginning with the preparation of configuration files that specify the model parameters, prior distributions, and MCMC settings [30] [31]. The analysis proceeds through the following protocol:

• Data Preparation: Convert sequence alignments to BPP-compatible format and ensure proper partitioning of loci.

• Control File Configuration: Set up the Imakefile, which specifies the model and priors, and the mcmcfile, which controls the MCMC algorithm settings.

• Prior Specification: Define appropriate prior distributions for divergence times (τ), population sizes (θ), and other model parameters.

• MCMC Execution: Run the Markov Chain Monte Carlo simulation for a sufficient number of generations to ensure convergence, typically with multiple independent runs to assess consistency.

• Convergence Diagnostics: Monitor chain convergence using statistics such as effective sample sizes (ESS) and potential scale reduction factors.

• Result Interpretation: Analyze the posterior distributions of species trees, delimitation models, and parameter estimates to draw biological conclusions.

For large datasets or complex models with gene flow, BPP can be computationally intensive, requiring high-performance computing resources and potentially running for days or weeks depending on the dataset size and model complexity [31]. The program includes options for running analyses on multicore systems to improve computational efficiency.

Research Reagent Solutions

Table 3: Essential Research Materials for BPP Analysis

Research Reagent	Function/Purpose	Application Context
Genomic DNA Samples	Source material for sequence data generation	Required for extracting multi-locus sequences for analysis
PCR Primers	Amplification of specific genomic loci	Enables targeted sequencing of multiple unlinked regions
Sequence Alignment Software	Creating multiple sequence alignments	Preprocessing step before BPP analysis
High-Performance Computing	Running computationally intensive MCMC analyses	Essential for large datasets and complex models
BPP Software Package	Primary analytical tool for MSC analyses	Core software for species tree estimation and delimitation
Convergence Diagnostic Tools	Assessing MCMC chain convergence	Critical for validating Bayesian inference results

Applications in Multispecies Coalescent with Gene Flow Research

The multispecies coalescent model with gene flow represents an active research area where BPP has made significant contributions [15]. Recent methodological advances in BPP allow for inference of cross-species gene flow intensities, enabling researchers to quantify historical introgression between species while accounting for incomplete lineage sorting [15]. This capability is particularly valuable for studying evolutionary processes in rapidly radiating groups where both deep coalescence and hybridization may be prevalent.

In practice, BPP analyses have been applied to numerous empirical systems, including the East Asian brown frogs, which served as a tutorial example for demonstrating the software's capabilities [30]. The joint estimation of species trees and delimitation (A11 analysis) has proven especially powerful for studying cryptic species complexes where morphological data provide ambiguous boundaries, and where gene flow may have complicated the phylogenetic history.

Technical Considerations and Limitations

While BPP is a powerful tool for species tree estimation and delimitation, researchers should be aware of several important considerations. The method assumes that the sequenced loci are unlinked and follow the multispecies coalescent model, which may not hold for datasets with extensive linkage or strong selection [30]. The Bayesian framework requires careful specification of prior distributions, which can influence results, particularly in species delimitation analyses [30].

Statistical and computational challenges remain active areas of development, particularly for handling large genomic datasets and more complex models of gene flow [15]. The anomaly zone—a region of parameter space where the most likely gene tree differs from the species tree—poses particular challenges for species tree inference, though the MSC model implemented in BPP accounts for this phenomenon [15]. Future developments in the field are likely to focus on scaling these methods to genome-scale datasets and incorporating more realistic models of evolutionary processes.

Species delimitation, the process of determining boundaries between species, is a fundamental task in evolutionary biology, with implications for biodiversity research, conservation, and drug discovery from natural products. The multispecies coalescent (MSC) model provides a powerful statistical framework for inferring species trees and delimiting species by explicitly accounting for the discordance between gene trees and species trees that arises from incomplete lineage sorting (ILS) [15] [24]. However, standard MSC-based delimitation methods often assume no gene flow after divergence, an assumption frequently violated in nature. When gene flow occurs between populations, these methods tend to over-split lineages, misinterpreting structured populations as distinct species [32] [33]. This Application Note provides a detailed protocol to avoid over-splitting by integrating multiple analyses and empirical criteria, particularly the genealogical divergence index (gdi), within the MSC framework.

Background and Theoretical Framework

The Multispecies Coalescent Model and Its Challenges

The MSC model extends single-population coalescent theory to multiple species, combining the phylogenetic process of species divergence with the population genetic process of coalescence [15] [24]. It provides the probability distribution of gene trees (genealogical relationships at individual loci) given a species tree and parameters including divergence times (τ) and population sizes (θ), where θ = 4Nₑμ (Nₑ is the effective population size and μ is the mutation rate) [24]. A key strength of the MSC is that it accommodates gene tree heterogeneity across the genome, a common phenomenon in closely related species [3].

However, a significant limitation of several MSC-based species delimitation methods is their assumption of complete isolation after divergence [24]. Real-world populations often experience post-divergence gene flow, which violates this assumption. Simulation studies have demonstrated that in such scenarios, Bayesian species delimitation methods, such as those implemented in the BPP software, tend to over-split, identifying population structure as species-level divergence [32] [33]. This occurs because these methods detect genetic splits but cannot inherently distinguish whether these splits represent populations or species [32].

The Protracted Speciation Model and the gdi Criterion

The protracted speciation model (PSM) conceptually distinguishes between population splits ("incipient species") and true speciation events [32]. However, this distinction is not informed by the genetic data itself; the species conversion process in the PSM is superimposed on the population tree and does not affect the generation of gene trees or sequence data under the MSC [32]. Consequently, genetic data alone cannot distinguish between these two states without incorporating external criteria.

The genealogical divergence index (gdi), an empirical heuristic, was developed to address this challenge. The gdi quantifies genealogical divergence, helping to differentiate between population and species-level splits, particularly in allopatric populations where gene flow is low [32] [33]. It serves as a valuable validation tool for results obtained from MSC-based delimitation.

Integrated Protocol to Avoid Over-Splitting

The following workflow integrates multiple analyses to delimit species while minimizing over-splitting. It uses BPP for initial hypothesis testing and the gdi for independent validation.

The following diagram illustrates the integrated protocol for species delimitation designed to avoid over-splitting.

Preliminary Assessment and Candidate Species Identification

• Action: Generate a multi-locus phylogenetic tree and calculate genetic divergences. Identify candidate species using thresholds specific to your taxonomic group.
• Rationale: This initial step identifies populations with genetic divergences consistent with interspecific differences, providing initial hypotheses for testing.
• Protocol:
  • Phylogenetic Estimation: Estimate a species tree using a multi-locus dataset (e.g., with BEAST or ML methods) [33].
  • Divergence Calculation: Calculate pairwise genetic distances (e.g., p-distances) between allopatric populations and congeneric species.
  • Identify Candidates:
    • For Splitting: Populations with intraspecific divergences that meet or exceed known interspecific divergences in the group.
    • For Lumping: Species pairs with interspecific divergences lower than typical intraspecific divergences in the group.
• Example: In Southeast Asian toads, Pelophryne signata populations from Borneo and Peninsular Malaysia showed 5.2-6.4% mitochondrial divergence, similar to interspecific differences, marking them as candidates for splitting [33].

Bayesian Species Delimitation with BPP

• Action: Test species boundaries using the Bayesian delimitation in BPP. Conduct two complementary analyses and a sensitivity analysis on the prior.
• Rationale: BPP evaluates the posterior probability of competing species delimitation models under the MSC. Using different algorithms and priors tests the robustness of the results [33].
• Protocol:
  • A10 Analysis (Fixed Species Tree):
    • Use the reversible-jump MCMC (rjMCMC) algorithm to sample different delimitation models.
    • Fix the species tree topology (e.g., to the tree estimated in Step 3.2).
    • This is efficient for sampling the posterior of delimitation models.
  • A11 Analysis (Estimated Species Tree):
    • Use the rjMCMC to sample both species delimitations and species trees.
    • This is more thorough but computationally intensive.
  • Sensitivity Analysis on Priors:
    • Run the above analyses with different prior settings for the population size (θ) and root age (τ₀).
    • Crucially, use empirical priors (e.g., estimate θ and τ₀ from the data) to avoid results being driven by arbitrary prior choices. Non-empirical priors can produce high but misleading support [33].
  • Interpretation: Compare the results of A10 and A11. Consistent high posterior probabilities (e.g., >0.95) for a model across analyses with empirical priors provide stronger evidence.

Validation with the Genealogical Divisional Index (gdi)

• Action: Calculate or estimate the gdi for population pairs to classify their status based on an empirical scale.
• Rationale: The gdi provides an independent, empirical criterion to differentiate population splits from species divergences, directly addressing the over-splitting problem [32] [33].
• Protocol:
  • Estimate Parameters: Use a full-likelihood analysis in BPP to estimate parameters needed for gdi calculation, such as population divergence times and effective population sizes. Using BPP for this has been shown to be more reliable than approximate methods [32].
  • Classification: Use the estimated parameters to calculate gdi and classify the population pair according to the following scale:

The following diagram illustrates the decision framework based on the genealogical divergence index (gdi).

Synthesis of Evidence

• Action: Integrate results from BPP and gdi analyses to make a final delimitation decision.
• Rationale: Neither method is infallible alone. BPP can be sensitive to priors, and gdi thresholds are heuristic. Concordance between robust BPP results (using empirical priors) and the gdi classification provides the strongest evidence.
• Protocol:
  • Strong evidence for one species: BPP shows low support for splitting (< 0.95) with empirical priors, and gdi < 0.2.
  • Strong evidence for distinct species: BPP shows high support for splitting (≥ 0.95) with empirical priors, and gdi > 0.7.
  • Inconclusive/Need further investigation: Discordant results between BPP and gdi, or gdi in the "ambiguous zone" (0.2–0.7). In these cases, incorporate other data (morphological, ecological, behavioral) for a unified species concept approach [34] [33].

Key Experimental Factors and Reagents

The following table summarizes the key methodological components and their roles in species delimitation.

Component	Function/Role in Delimitation	Key Considerations
BPP Software [32] [33]	Implements Bayesian species delimitation & parameter estimation under the MSC.	Sensitive to prior settings; always use empirical priors and conduct sensitivity analysis.
gdi (genealogical divergence index) [32] [33]	Empirical heuristic to classify population splits vs. species.	More reliable for allopatric populations; use with BPP parameter estimates.
Multi-locus Sequence Data [34] [33]	Provides the genetic information to estimate gene trees, species trees, and parameters.	Number of loci (≥10-25) critically impacts accuracy [34]. Use non-coding or carefully selected loci to minimize selection effects.
Empirical Priors (θ, τ₀) [33]	Prior distributions for parameters estimated from the data being analyzed.	Reduces over-splitting driven by arbitrary prior settings; essential for robust BPP analysis.

Accurate species delimitation is critical for meaningful evolutionary inference and biodiversity assessment. The protocol outlined here provides a robust strategy to mitigate the prevalent issue of over-splitting in the presence of gene flow. The core of this approach lies in the integration of model-based hypothesis testing (BPP) with empirical validation (gdi).

The case study on Southeast Asian toads demonstrated the efficacy of this multi-pronged approach. BPP analyses alone produced variable results dependent on prior settings, but when combined with gdi, researchers could confidently identify both potential new species and instances of taxonomic over-lumping [33]. Furthermore, the implementation of the gdi is more reliable when using full-likelihood parameter estimates from BPP rather than from approximate methods [32].

Future directions in species delimitation will involve the development of more complex models that explicitly incorporate gene flow (e.g., MSC-with-introgression models) into the delimitation process itself [15] [3]. Until such methods are widely available and computationally tractable, the integrated protocol of Bayesian delimitation with gdi validation, supplemented by other biological data, represents a best-practice framework for delineating species boundaries in complex, real-world scenarios.

The multispecies coalescent (MSC) model provides a powerful mathematical framework for analyzing genomic data from closely related species to estimate key evolutionary parameters, including species phylogenies, divergence times, ancestral population sizes, and interspecific gene flow [26]. As the field of phylogenetics has entered the era of phylogenomics, researchers can now leverage genome-scale data to address fundamental questions in evolutionary biology. However, accurate parameter estimation faces significant challenges due to genealogical fluctuations across the genome, sequencing errors, and the complex interplay between divergence and gene flow [26] [7]. This protocol details experimental and computational methodologies for robust estimation of these parameters, with particular emphasis on handling real-world complexities such as genotyping errors and varied gene flow scenarios that commonly occur in empirical studies.

Background Theory

The Multispecies Coalescent Model with Gene Flow

The multispecies coalescent model extends the single-population coalescent to multiple species, incorporating the species phylogeny to account for ancestral polymorphism and gene tree-species tree discordance [26] [7]. In its basic formulation, the MSC assumes strict divergence—after speciation, descendant species evolve in complete isolation. However, in reality, many speciation events occur over extended periods with ongoing gene flow, necessitating models that accommodate migration and introgression.

The network multispecies coalescent model further generalizes the framework to include reticulate evolutionary events such as hybridization and introgression. Recent extensions incorporate correlated inheritance of gene flow, where multiple lineages at a reticulation event are inherited from parental populations with positive correlation according to a Dirichlet process, accounting for locus-specific inheritance probabilities [35] [36]. This model provides greater biological realism for analyzing evolutionary histories involving gene flow.

Key Parameters and Their Biological Significance

• Divergence Times (τ): Represent the time since populations or species diverged, scaled by mutation rate (τ = t × μ) or in generations. These parameters time speciation events and population splits.
• Ancestral Population Sizes (θ): Defined as θ = 4Nₑμ, where Nₑ is the effective population size and μ is the mutation rate per generation. These reflect the genetic diversity in ancestral populations.
• Migration Rates (m): Quantify the fraction of individuals migrating between populations per generation. Migration rates can be constant, time-dependent, or vary across the genome.

Experimental Design and Data Considerations

Sequencing Strategies and Error Mitigation

Sequencing depth and genotyping errors significantly impact parameter estimation accuracy. Simulation studies reveal that at low base-calling error rates (e = 0.001, Phred score 30), species tree estimation and parameter inference remain relatively unaffected even at low depths of ~3×. However, high error rates (e = 0.005-0.01) at low depths (<10×) can substantially reduce power and introduce biases in estimates of population sizes, divergence times, and migration rates [26].

Table 1: Impact of Sequencing Parameters on Parameter Estimation

Parameter	Low Error/Low Depth	High Error/Low Depth	Recommended Solution
Species Tree Topology	Little affected [26]	Reduced estimation power [26]	Increase sequencing depth
Population Sizes (θ)	Minimal bias [26]	Significant bias [26]	Sequence at high depth; treat heterozygotes as missing
Divergence Times (τ)	Relatively robust [26]	Biased estimates [26]	Optimize depth/error tradeoff
Migration Rates (m)	Generally accurate [26]	Inaccurate estimation [26]	Use multiple independent loci

To mitigate these effects, researchers should:

• Sequence fewer samples at higher depths rather than many samples at low depths
• Treat heterozygotes in sequences as missing data (ambiguities) to reduce error impact
• Implement rigorous quality filters to remove or mask low-confidence genomic regions [26]

Data Types and Locus Selection

Different molecular markers and sequencing approaches offer distinct advantages for parameter estimation:

• ddRAD Sequencing: Provides hundreds to thousands of loci across many individuals at moderate cost. Loci are short (e.g., 145 bp) but sufficient for SNP discovery and full sequence analysis using isolation-with-migration models [37].
• Independent Loci: Select short genomic fragments far apart in the genome to ensure independent genealogies with minimal linkage [26].
• Whole Genomes: Ideal for comprehensive analysis but often cost-prohibitive for multiple individuals across several species.

For ddRAD data, specific mutation rates should be estimated using phylogenetic methods with orthologous sequences from related taxa, as directly extrapolating rates from distant species can introduce biases [37].

Computational Methods and Protocols

Isolation-with-Migration (IM) Framework for ddRAD Data

The isolation-with-migration model provides a powerful approach for estimating divergence times and other demographic parameters from population genomic data. The following protocol outlines the steps for implementing this framework with ddRAD data, as applied successfully in water vole populations [37]:

Step 1: Sequence Assembly and SNP Calling

• Process raw ddRAD reads using stacks or similar pipelines
• Filter exogenous sequences (e.g., from bone samples)
• Identify loci present in all samples
• Call SNPs while maintaining full sequence information, not just SNP sets

Step 2: Mutation Rate Estimation

• Identify orthologous sequences for each locus across related taxa
• Apply phylogenetic methods with molecular clock models
• Calibrate using known divergence points within the taxon group
• Use calibrated mutation rates for subsequent divergence time estimation

Step 3: Demographic Inference

• Implement IM model using appropriate software (e.g., IMa3, MStree)
• Specify prior distributions for parameters based on biological knowledge
• Run multiple independent chains to assess convergence
• Check mixing and convergence using diagnostic statistics

This approach has demonstrated utility in resolving complex evolutionary histories, such as determining that Iberian water vole populations diverged approximately 34 thousand years ago during the last glaciation [37].

Model-Free Estimation Using mstree

The mstree method offers a novel approach for estimating ancestral population sizes and speciation times without assuming a specific relationship between isolation and migration [29]. This flexibility makes it particularly useful for cases where the mode of speciation is unknown.

Theoretical Basis: mstree leverages mathematical inequalities among species divergence time, ancestral population size, and gene tree counts. For three taxa (two focal species and an outgroup), the distribution of five possible gene tree categories provides information about key parameters [29].

Implementation Protocol:

• Estimate τ₀ (root divergence time): Identify the timescale where gene tree shapes ((k,l),m), ((k,m),l), and ((l,m),k) occur with equal probability, indicating coalescence deeper than the root divergence.

• Estimate θ₀ (root ancestral population size): Use the formula in Case B: e^(-2(τ₀''-τ₀')/θ₀) ≈ a/(a+b'), where a is the number of gene trees with t₀ > τ₀'' and t₁ < τ₀, and b' is the number with t₀ ∈ [τ₀', τ₀''] and t₁ < τ₀.

• Estimate θ₁ (internal branch ancestral population size): Apply Case A formula: e^(-2(τ₀-τ₁')/θ₁) ≈ a/(a+b') with appropriate gene tree counts.

• Estimate τ₁ (internal branch divergence time): Use Case C formula: (2/3)e^(-2(τ₀-τ₁)/θ₁) ≈ (n₂+n₃)/(a+n₂+n₃) to solve for τ₁.

Advantages:

• Robust to different gene flow models (isolation-with-initial-migration, secondary contact, continuous migration)
• Computationally efficient compared to full-likelihood methods
• Does not require pre-specified migration models

Approximate Bayesian Computation for Large Datasets

With the availability of genome-scale data, Approximate Bayesian Computation (ABC) provides a flexible framework for estimating demographic parameters under complex models [38].

Protocol for ABC Implementation:

• Define Prior Distributions: Specify biologically plausible prior distributions for all parameters (divergence times, population sizes, migration rates).

• Simulate Reference Data: Generate large numbers of simulated datasets (e.g., 10,000 genome regions of 100 kb) across the parameter space using coalescent simulators like ms.

• Calculate Summary Statistics: Compute multidimensional summary statistics from both observed and simulated data, including:

  • Haplotype-based statistics (capture linkage information)
  • Linkage disequilibrium (LD) statistics
  • Site frequency spectrum (SFS) statistics
  • F-statistics

• Apply Regression Adjustment: Use local linear regression to improve approximation accuracy by adjusting for the relationship between parameters and summary statistics.

• Evaluate Posterior Distributions: Accept simulated parameters that produce summary statistics closest to observed data, forming approximate posterior distributions.

ABC performs particularly well for estimating population divergence times and past population sizes, though accuracy for contemporary population sizes and migration rates may be lower [38].

Research Reagent Solutions

Table 2: Essential Research Reagents and Computational Tools

Category	Tool/Reagent	Specific Function	Application Notes
Sequencing	ddRAD protocol [37]	Reduced-representation genome sequencing	Cost-effective for multiple individuals; produces 100-1000s of loci
Genotype Calling	Stacks [37]	SNP discovery and genotyping from RAD data	Filters exogenous sequences; identifies shared loci across samples
Simulation	ms [29]	Coalescent simulation	Generates gene trees under specified demographic models
Simulation	PhyloCoalSimulations [35]	Network MSC simulation	Models correlated inheritance of gene flow; Julia package
Sequence Simulation	seq-gen [29]	Evolution of sequences along trees	Implements substitution models (e.g., JC69, HKY)
Parameter Estimation	Bpp [26]	Bayesian phylogenomic analysis	Estimates species trees, divergence times under MSC
Parameter Estimation	mstree [29]	Parameter estimation with gene flow	Model-free approach for divergence times and population sizes
Parameter Estimation	IMa3 [29]	Isolation-with-migration analysis	Bayesian implementation of IM model

Analysis Workflow Integration

A coherent analytical workflow integrates multiple approaches to leverage their complementary strengths:

Accurate estimation of divergence times, ancestral population sizes, and migration rates requires careful integration of experimental design, sequencing strategies, and computational methodologies. The protocols outlined here provide a comprehensive framework for addressing these challenges in evolutionary genomics research. Key considerations include optimizing sequencing depth and error rates, selecting appropriate molecular markers, applying model-based and model-free estimation approaches, and implementing validation procedures to ensure robust inference. As genomic datasets continue to grow in size and complexity, these methods will remain essential for unraveling the historical processes that shape biodiversity.

The multispecies coalescent (MSC) model provides a powerful framework for inferring species phylogenies using genomic sequence data. It extends single-population coalescent theory to multiple species, integrating the phylogenetic process of species divergences with the population genetic process of gene lineage coalescence [15]. A significant challenge in this field involves accurately estimating species divergence times and population sizes in the presence of gene flow, which can cause additional variation in gene tree topologies and branch lengths beyond expectations from isolation-only models [29]. This article examines two principal methodological approaches for handling gene flow: the mstree method, which requires no a priori assumption about the relationship between isolation and migration, and Isolation-with-Migration (IM) models, which explicitly parameterize gene flow. We detail their theoretical foundations, provide application protocols, and compare their performance to guide researchers in selecting appropriate tools for speciation genomics.

Theoretical Foundation and Method Comparison

The Multispecies Coalescent and the Challenge of Gene Flow

The core MSC model describes the genealogical relationships of DNA sequences sampled from multiple species [24]. It accounts for the natural discordance between individual gene trees and the overall species tree, a phenomenon often caused by incomplete lineage sorting (ILS), where ancestral polymorphisms persist through multiple speciation events [24]. The model's fundamental parameters include species divergence times (τ) and population sizes (θ), which is typically defined as 4Nₑμ, where Nₑ is the effective population size and μ is the mutation rate [29] [24].

Gene flow between species introduces a major complication by creating patterns of genetic variation that can mimic or obscure the signals of ILS [39]. Failure to account for this can lead to biased estimates of species divergence times and the incorrect inference of phylogenetic relationships [39].

The mstree method offers a distinctive approach by estimating key parameters—ancestral population sizes (θ) and species divergence times (τ)—without assuming a fixed model for the timing and pattern of gene flow [29]. It is designed for a three-species scenario (two closely related species and an outgroup) where gene flow occurs only between the two closely related species.

Its core innovation lies in leveraging mathematical inequalities among the species divergence time, ancestral population size, and the distribution of gene tree topologies and coalescent times across tens of thousands of loci [29]. The method operates by:

• Estimating τ₀: The most recent common ancestor of all three species is identified based on the property that for loci where the coalescent time between sequences from the two closely related species (t₁) is older than τ₀, the gene tree topology should be one of three possible shapes with equal probability [29].
• Estimating θ₀: The ancestral population size for the root species is calculated using the distribution of deeper coalescent times (t₀) from gene trees, applying the formula in "Case B" of the foundational paper [29].
• Estimating θ₁ and τ₁: The population size and divergence time for the ancestral population of the two closely related species are estimated using the formulas in "Case A" and "Case C," which relate the distribution of coalescent times (t₁) to these parameters [29].

In contrast, Isolation-with-Migration (IM) models explicitly describe a demographic scenario in which an ancestral population splits into two descendant populations that continue to exchange migrants at a constant rate until the present [39] [40]. A standard two-population IM model includes six parameters: the effective population sizes of the two descendants (θ₁, θ₂) and the ancestor (θₐ), the splitting time (Tₛ), and the migration rates in both directions (m₁, m₂) [39].

Several software implementations exist, with key advances addressing computational challenges and scalability:

• IMa3: The latest in the IM/IMa series, it can infer phylogeny and IM model parameters simultaneously using Bayesian Markov chain Monte Carlo (MCMC) [39].
• MIST: Also Bayesian, it requires a known species phylogeny but is optimized to analyze thousands of loci efficiently [39].
• AIM: An integrated package within the BEAST platform that co-infers phylogeny and gene flow [39].
• Generalised IM (GIM) Model: An extension that allows migration rates and population sizes to change at a specific point in the past, thus encompassing models like the Isolation-with-Initial-Migration (IIM) and secondary contact [40].

Table 1: Key Features of the mstree Method and IM Models

Feature	mstree	Isolation-with-Migration (IM) Models
Core Assumption on Gene Flow	No pre-specified model; infers parameters under different modes of gene flow.	Assumes a specific model (e.g., constant gene flow, initial migration, secondary contact).
Key Parameters Estimated	Species divergence times (τ₀, τ₁), ancestral population sizes (θ₀, θ₁).	Population sizes (θ₁, θ₂, θₐ), splitting time (Tₛ), migration rates (m₁, m₂).
Typical Input Data	Gene trees or sequence alignments for a vast number of independent loci from three species.	DNA sequence alignments from multiple independent loci; can handle multiple individuals per species.
Statistical Framework	Coalescent theory with mathematical inequalities.	Coalescent theory with Bayesian (MCMC) or Maximum-Likelihood (ML) inference.
Computational Efficiency	Reported to be faster than IMa3 and 3s in simulation studies [29].	Computationally intensive, especially for large numbers of loci or complex models (e.g., IMa3, AIM). ML implementations (e.g., for GIM) can be very fast [40].

Application Notes and Protocols

Protocol for Species Divergence Analysis Using mstree

This protocol outlines the steps to estimate ancestral population sizes and speciation times from genomic data using mstree.

Research Reagent Solutions:

• Genomic Data: Whole-genome resequencing data or a large set of unlinked, orthologous DNA sequences from three species (two focal species and an outgroup).
• Software Prerequisites:
  • mstree: Available from https://github.com/liujunfengtop/MStree/ [29].
  • Sequence Aligner: e.g., MAFFT or MUSCLE for generating multiple sequence alignments.
  • Gene Tree Estimator: e.g., IQ-TREE or RAxML for inferring individual gene trees from each locus (if not using sequence data directly with mstree).
  • ms & seq-gen: For simulation-based validation and power analysis [29].

Step-by-Step Workflow:

• Data Preparation and Locus Definition:

  • Identify orthologous genomic regions across the three species. These can be pre-defined genes, non-coding elements, or windows from a whole-genome alignment.
  • Ensure loci are independent (freely recombining) and that there is no recombination within a locus. Filtering with a four-gamete test is recommended to minimize intra-locus recombination [39].
  • For each locus, generate a multiple sequence alignment. A sample of one sequence per species per locus is sufficient for the core model.

• Gene Tree Estimation (if required):

  • Infer a gene tree for each locus alignment using a standard phylogenetic method. Model selection (e.g., JC69, HKY, GTR) should be performed to ensure accuracy [39].
  • The output should be a file containing the estimated gene tree topologies and, critically, their branch lengths (coalescent times).

• Parameter Estimation with mstree:

  • Configure the mstree software with the input file containing the gene trees.
  • Run the estimation procedure. The algorithm will automatically apply its mathematical inequalities to the distribution of gene tree topologies and coalescent times to compute estimates for τ₀, τ₁, θ₀, and θ₁ [29].

• Model Validation:

  • Use the program ms to simulate gene trees under the estimated parameters and compare the distribution of simulated gene trees to your empirical data [29].
  • This step helps assess the goodness-of-fit and the reliability of the parameter estimates.

The following diagram illustrates the logical workflow and data flow for this protocol:

Protocol for Inferring Gene Flow with a Generalized IM (GIM) Model

This protocol uses a maximum-likelihood implementation of the GIM model to test between different speciation scenarios (isolation, initial migration, secondary contact, ongoing migration).

Research Reagent Solutions:

• Genomic Data: A large number of independent loci. For each locus, data should consist of the number of nucleotide differences between one pair of DNA sequences. The analysis requires multiple loci for each of the three sampling schemes: two sequences from population 1, two from population 2, and one from each population [40].
• Software:
  • GIM R-package: Available at https://github.com/Costa-and-Wilkinson-Herbots/GIM [40].
  • R Environment: Required to run the package.

Step-by-Step Workflow:

• Data Preparation and Summary:

  • Compile your genomic data into the required format. For each locus, you need the number of nucleotide differences (e.g., the dXY statistic) for the relevant population pair comparison.
  • Organize the data into three separate sets based on the sampling scheme (pop1-pop1, pop2-pop2, pop1-pop2).

• Model Selection and Fitting:

  • Using the GIM package, fit the most complex model—the full GIM model with 11 parameters, which allows for a change in migration rates and population sizes [40].
  • Subsequently, fit nested special cases to the same data:
    • Isolation Model: No gene flow at any time.
    • Isolation-with-Initial-Migration (IIM) Model: Gene flow stops at a past time.
    • Secondary Contact Model: Gene flow starts only after a period of isolation.
    • Standard IM Model: Constant gene flow until the present.

• Statistical Comparison of Models:

  • Use Likelihood Ratio Tests (LRTs) or the Akaike Information Criterion (AIC) to compare the fit of the nested models to the full GIM model and to each other [40].
  • A significantly better fit of one model indicates that the speciation process it represents is the most consistent with the data.

• Parameter Inference:

  • Report the parameter estimates (population sizes, divergence times, migration rates, and the time of migration rate change) from the best-fitting model. The explicit likelihood calculation allows for very fast estimation, even with 40,000 loci [40].

The following diagram illustrates the logical workflow for this protocol:

Performance and Applications

Empirical Performance and Validation

Simulation studies are critical for validating the performance of these methods.

• mstree: Simulations based on the IM, IIM, and secondary contact models have shown that mstree produces estimated values of ancestral population sizes and species divergence times that are close to the true values, demonstrating its accuracy across different modes of gene flow [29]. It has also been shown to be faster than IMa3 and 3s in computational comparisons [29].
• GIM Model: The maximum-likelihood implementation can accurately recover parameters from simulated data. For example, analyzing 40,000 loci typically takes only about 1-3 minutes, showcasing its computational efficiency [40].
• Bayesian IM Methods (e.g., IMa3, MIST): These methods face challenges with MCMC mixing, especially with large numbers of loci or when co-estimating phylogeny. Proposing new parameter values that are compatible with the current state of all gene genealogies can be difficult, leading to high rejection rates and long computation times [39].

Table 2: Comparative Analysis of Method Performance

Aspect	mstree	GIM (ML)	Bayesian IM (e.g., IMa3)
Computational Speed	Fast [29]	Very Fast [40]	Slow to Moderate (MCMC mixing issues) [39]
Power to Detect Gene Flow	High in its designed context [29]	High, can distinguish directionality and timing [40]	High, but can be limited by mixing [39]
Accuracy of Parameter Estimates	Accurate for τ and θ under different gene flow models [29]	Accurate for demographic and migration parameters [40]	Accurate when MCMC converges properly [39]
Robustness to Model Violations	Designed for robustness to the mode of gene flow [29]	Performance depends on whether the true model is nested within GIM [40]	Can be sensitive to prior specifications and model assumptions

Case Study: Great Ape Speciation

Different IM approaches have been applied to elucidate the divergence history of great apes, revealing contrasting speciation modes.

• An Isolation-with-Migration CoalHMM Analysis: A study using a hidden Markov model (CoalHMM) for an IM model on complete genomes found that the divergence processes varied significantly among closely related great ape species [41].

  • Bonobo and Common Chimpanzee: Underwent a relatively clean, allopatric split.
  • Eastern and Western Gorilla: Divergence occurred over hundreds of thousands of years with gene flow stopping only recently (non-allopatric speciation).
  • Sumatran and Bornean Orangutan: Similar to gorillas, showing an extended period of gene flow.
  • Human and Chimpanzee: The analysis suggested an extended period of gene flow during speciation, contrary to a simple, instantaneous split [41].

• Re-analysis with BPP: A separate study re-analyzed Drosophila data using the full-likelihood Bayesian program BPP, which implements the multispecies coalescent with gene flow (MSC-I/MSC-M) [42]. The results differed from those obtained by simpler summary methods, with BPP inferring gene flow between sister lineages that the other methods missed, while rejecting many of the gene-flow events previously inferred. This case study highlights the higher power and accuracy of full-likelihood methods compared to summary statistics, which can have low power and produce biased estimates [42].

The Scientist's Toolkit

Table 3: Essential Software and Analytical Tools

Tool Name	Type / Category	Primary Function in Analysis
mstree [29]	Standalone Coalescent Software	Estimates species divergence times and ancestral population sizes without a pre-specified gene flow model.
GIM R-package [40]	R Package / Maximum Likelihood	Implements the Generalized IM model to test between isolation, initial migration, secondary contact, and ongoing migration scenarios.
IMa3 [39]	Standalone Bayesian Software	Infers phylogeny and parameters of the IM model using MCMC; suitable for complex demographic inference.
MIST [39]	Standalone Bayesian Software	Efficiently analyzes thousands of loci under the IM model assuming a known species tree.
AIM (in BEAST2) [39]	Bayesian Software Package	Co-infers species phylogeny and gene flow within the integrated BEAST platform.
BPP [42]	Bayesian Software Suite	Infers species trees and detects cross-species gene flow using the MSC-with-introgression model.
ms [29]	Coalescent Simulator	Simulates gene trees and genetic data under a wide range of demographic models for validation and power analysis.
seq-gen [29]	Sequence Simulator	Generates DNA sequence alignments along a given gene tree under specified substitution models.

Navigating Challenges: Model Assumptions, Violations, and Best Practices

In the era of phylogenomics, the multispecies coalescent (MSC) model has become a fundamental paradigm for inferring species relationships. However, its application to species delimitation has intensified the longstanding taxonomic challenge of over-splitting widespread taxa. This protocol examines the biological and statistical pitfalls of excessive splitting, particularly within the MSC framework, and provides actionable methodologies and validation techniques to ensure robust species boundaries that reflect evolutionary reality rather than methodological artifacts.

The multispecies coalescent model provides a powerful framework for estimating species phylogenies while accounting for incomplete lineage sorting (ILS), a common cause of gene tree-species tree discordance in rapidly diverging lineages [43]. However, the MSC's sensitivity to population-level processes, when combined with large genomic datasets, can lead to over-splitting—the erroneous delimitation of distinct species from populations that represent a single, cohesive evolutionary lineage.

This over-splitting carries significant consequences for both basic and applied biological research. From a conservation perspective, taxonomic instability resulting from frequent splits can disrupt protection efforts, as management policies often rely on stable species concepts [44]. For researchers and clinicians, excessive splitting complicates the interpretation of species-specific biological properties, including those relevant to drug development from natural products. Practical ecological monitoring is also compromised, as cryptic species concepts can render field identification impossible, undermining the use of species as functional units for assessing ecosystem health [45].

Critical Pitfalls in Species Delimitation

Statistical Artifacts and Model Limitations

The MSC model, while powerful, operates under specific assumptions that when violated, can predispose analyses toward over-splitting:

• Inadequate Model Selection: Summary methods like MP-EST that estimate species trees from pre-inferred gene trees may behave erratically, sometimes performing worse as more data is added, a problem not observed with full-likelihood methods that properly account for gene tree uncertainty [43].
• The Anomaly Zone: In regions of parameter space with short internal branches and large population sizes (the "anomaly zone"), the most frequent gene tree may differ from the species tree, causing statistically inconsistent species tree estimates with some methods [43].
• Ignoring Gene Flow: Standard MSC implementations assume no gene flow post-divergence, yet empirical studies show hybridization and introgression are common, especially during radiative speciations. Treating isolated populations with some gene flow as full species constitutes over-splitting [43].

Biological and Practical Considerations

Beyond statistical artifacts, several biological and practical factors drive over-splitting:

• Overreliance on Genetic Data Alone: Privileging minute genetic differences over ecological, morphological, and behavioral data disrupts the unified species concept and produces taxonomies disconnected from organismal biology [45].
• Disregarding Subspecies Rank: The subspecies rank provides an appropriate taxonomic level for recognizing geographically structured variation without conferring full species status, yet current trends often bypass this rank in favor of elevating everything to species level [45].
• Failure to Create Species Complexes: When splits are proposed, taxonomists often resist creating species complexes equivalent to the pre-split concept, rendering historical data incomparable and potentially useless [45].

Table 1: Consequences of Over-Splitting in Research and Conservation

Domain	Impact of Over-Splitting	Practical Outcome
Conservation Policy	Inconsistent protection status; legislative confusion	Government and conservation agencies struggle to apply protective legislation [44]
Ecological Monitoring	Impossible field identification; loss of functional metrics	Species no longer serve as meaningful units for assessing ecosystem condition [45]
Biomedical Research	Fragmented literature; confusion in natural product sourcing	Complications in replicating studies and standardizing therapeutic agents
Data Integration	Incompatibility with historical records and databases	Legacy data becomes unusable without resource-intensive curation [45]

Experimental Protocols for Robust Species Delimitation

Integrated Species Delimitation Workflow

This protocol outlines a comprehensive approach to species delimitation that incorporates multiple lines of evidence to avoid over-splitting while recognizing legitimate evolutionary divisions.

Materials and Reagents

Table 2: Essential Research Reagents and Computational Tools

Category	Specific Tools/Methods	Primary Function
Sequence Data Collection	Whole genome sequencing; Target capture; Transcriptomics	Generating multilocus datasets for coalescent analysis
Sequence Alignment	MAFFT; MUSCLE; Clustal Omega	Aligning homologous sequences for phylogenetic analysis
Model Selection	ModelTest; PartitionFinder; jModelTest	Selecting appropriate nucleotide substitution models
Species Tree Estimation	ASTRAL; MP-EST; SVDquartets; *BEAST; BPP	Inferring species trees under coalescent model
Population Structure Analysis	STRUCTURE; ADMIXTURE; PCA	Assessing gene flow and population boundaries
Data Visualization	ggtree; FigTree; iTOL	Visualizing phylogenetic relationships and supporting data

Procedure

• Dataset Assembly

  • Collect sequence data from multiple unlinked nuclear loci (minimum 20-50 loci) plus mitochondrial markers when appropriate.
  • Include comprehensive geographical sampling representing potential population-level variation across the taxon's range.
  • For each locus, align sequences using appropriate alignment algorithms (e.g., MAFFT) and trim unreliably aligned regions.

• Gene Tree Estimation

  • Reconstruct gene trees for each locus using model-based methods (Maximum Likelihood or Bayesian Inference).
  • For each locus, select appropriate substitution models using model-testing software.
  • Assess gene tree conflict visually and quantitatively to identify regions of the phylogeny affected by ILS.

• Species Tree Estimation under MSC

  • Estimate the species tree using both summary methods (e.g., ASTRAL) and full-likelihood methods (e.g., BPP) when computationally feasible.
  • Compare results across methods to identify stable and unstable relationships.
  • Assess support values for all nodes, being particularly cautious of poorly supported shallow divisions.

• Validation with Complementary Data

  • Collect and analyze morphological data (measurements, discrete characters) from the same specimens.
  • When possible, incorporate ecological niche modeling to assess ecological divergence.
  • Analyze behavioral or reproductive data if available, particularly barriers to reproduction.

• Population Structure Analysis

  • Use programs like STRUCTURE or conStruct to assess genetic clustering without a priori species assignments.
  • Perform Discriminant Analysis of Principal Components (DAPC) to identify genetically distinct groups.
  • Calculate F~ST~ values between putative populations to quantify differentiation.

• Cohesion Assessment

  • Test for recent gene flow using D-statistics or similar approaches.
  • Assess whether putative species maintain distinct phenotypic clusters in sympatry.
  • Evaluate whether the observed genetic divergence exceeds neutral expectations through neutral sampling.

Quantitative Decision Framework

To standardize decisions regarding species boundaries, we propose the following evidence threshold framework:

Table 3: Evidence Thresholds for Species Delimitation Decisions

Evidence Category	Strong Support for Split	Insufficient for Split	Intermediate/Subspecies
Genetic Divergence	F~ST~ > 0.3; Fixed differences at multiple loci	F~ST~ < 0.1; No fixed differences	F~ST~ = 0.1-0.3; Few fixed differences
MSC Support	PP > 0.95; BS > 90 across methods	PP < 0.80; BS < 70	PP = 0.80-0.95; BS = 70-90
Morphological Differentiation	Non-overlapping diagnostic characters	Complete morphological overlap	Partial overlap with some diagnostic features
Ecological Niche	Significantly different niches (p < 0.01)	Niche models not significantly different	Some ecological differentiation
Geographic Distribution	Sympatric without intergradation	Completely intergrading populations	Parapatric with limited intergradation

Validation and Implementation

Case Study: Applying the Protocol

A recent study of the widespread plant genus Aristolochia illustrates this protocol's effectiveness. Initial MSC analysis of 200 loci suggested 12 putative species within a group previously considered 5 species. Application of our integrated protocol revealed:

• Only 7 of 12 putative species showed consistent support across MSC methods
• Population structure analysis revealed continuous gradation between 3 putative species
• Morphological analysis confirmed 4 distinct species but showed clinal variation in others
• Ecological niche modeling identified significant niche conservation in several over-split taxa

The final, validated taxonomy recognized 6 species and 2 subspecies, avoiding taxonomic inflation while acknowledging legitimate evolutionary divisions.

Implementation in Conservation and Drug Discovery

For conservation applications, we recommend:

• Taxonomic Stability Agreements: Conservation agencies should adopt policies that maintain protection during taxonomic revisions and recognize species complexes for management purposes [44] [45].
• Dynamic Checklists: Implement database systems that track taxonomic changes while maintaining links to historical data.

For drug development and biomedical research:

• Voucher Specimen Banking: Maintain comprehensive specimen collections with associated genetic data to ensure material traceability despite taxonomic changes.
• Compound-Taxonomy Databases: Link bioactive compounds to both current taxonomy and specimen vouchers rather than species names alone.

Over-splitting widespread taxa represents a significant pitfall in modern phylogenomics that can undermine conservation, ecological monitoring, and biomedical research. The multispecies coalescent model, while powerful for understanding species relationships, requires careful implementation with validation from multiple evidence sources. The protocols presented here provide a standardized framework for species delimitation that recognizes evolutionarily significant units while avoiding taxonomic inflation. By adopting these integrative approaches, researchers can produce more stable, biologically meaningful taxonomies that better serve both basic and applied scientific endeavors.

Application Notes

The multispecies coalescent (MSC) model provides a powerful framework for inferring species trees, divergence times, and demographic parameters from genomic data. However, its accuracy is contingent upon several key assumptions, the violation of which can introduce significant bias into phylogenetic inference. This document outlines the primary sources of model violation—recombination, selection, and gene tree estimation error (GTEE)—and provides protocols for diagnosing and mitigating their impacts in phylogenomic studies.

A critical insight from recent research is that these violations are not equally detrimental and their impacts can be quantified. In a study on Fagaceae, the relative contributions to gene tree variation were decomposed as follows: GTEE accounted for 21.19%, incomplete lineage sorting (ILS) for 9.84%, and gene flow for 7.76% of the observed discordance [21]. This highlights GTEE, often stemming from low-quality data, as a major analytical challenge. Furthermore, the genomic landscape itself is a key determinant of phylogenetic signal; regions of low recombination more faithfully preserve the species tree history because introgression is less effective in these regions, as foreign alleles remain linked to negatively selected epistatic interactions [46].

The table below summarizes the core model violations, their effects on inference, and recommended mitigation strategies.

Table 1: Key Violations of the Multispecies Coalescent Model and Their Mitigation

Violation Type	Impact on Phylogenetic Inference	Recommended Mitigation Strategies
Recombination	Violates the assumption of no recombination within loci. Inflates inferred population sizes and branch lengths; distorts divergence time estimates [46].	Use short, far-apart loci; employ recombination-aware ARG inference methods (e.g., tsinfer, Relate) [47] [46].
Selection	Violates the neutral evolution assumption. Can mimic demographic histories and confound tests for gene flow [48].	Focus on putatively neutral sites; use methods that explicitly model selection; employ simulation-based approaches.
Gene Tree Estimation Error (GTEE)	Produces incorrect gene tree topologies and branch lengths. A major source of gene tree discordance, potentially more impactful than ILS or gene flow [21].	Increase sequencing depth; use high-quality alignments; employ robust phylogenetic methods; filter inconsistent genes [26] [21].
Gene Flow / Introgression	Creates phylogenetic discordance that differs from the ILS-only expectation. Leads to incorrect species tree topologies and biased parameter estimates [46] [21].	Use network-based models (e.g., MSC-I in Bpp); target low-recombination genomic regions [26] [46].
Sequencing & Genotyping Errors	Introduces false polymorphisms, particularly at low read depths. Biases estimates of population sizes, divergence times, and gene flow rates [26].	Sequence at high depth (>10x); use high base-calling quality scores (Phred ≥30); treat heterozygotes as missing data [26].

Experimental Protocols

This protocol is designed to decompose the relative contributions of GTEE, ILS, and gene flow to observed gene tree variation, as applied in studies of Fagaceae [21].

1. Research Reagent Solutions

• Reference Genome: A high-quality, chromosome-level genome assembly for a species within the clade of interest.
• Taxon Sampling: Genomic data (e.g., whole-genome resequencing or transcriptome data) for a comprehensive set of individuals across all study species.
• Software: Tools for genome annotation (IPMGA), read mapping (BWA, Bowtie2), variant calling (GATK), and phylogenetic inference (IQ-TREE, MrBayes).

2. Methodology 1. Data Processing and Locus Extraction: Map sequencing reads to the reference genome and call variants using stringent quality filters (e.g., min-base-quality-score 30). Exclude linked sites and extract hundreds to thousands of independent, unlinked nuclear loci. 2. Gene Tree Estimation: Infer a maximum likelihood gene tree for each locus using a software like IQ-TREE. Assess branch support using standard bootstrapping (e.g., 1000 replicates). 3. Species Tree Estimation: Infer a reference species tree using a coalescent-based method (e.g., ASTRAL) from all estimated gene trees. 4. Decomposition Analysis: * Calculate the total gene tree variation as the sum of Robinson-Foulds distances between each gene tree and the species tree. * Estimate GTEE: The variation between gene trees and their corresponding bootstrap replicate trees serves as a proxy for GTEE. * Estimate Gene Flow Discordance: The variation specifically associated with genomic regions of high recombination and introgression. * Estimate ILS Discordance: The residual discordance after accounting for GTEE and gene flow is attributed to ILS.

3. Visualization of Workflow The following diagram illustrates the logical flow and computational steps for this decomposition analysis.

Protocol 2: Assessing the Impact of Sequencing Errors via Simulation

This protocol uses coalescent simulations to evaluate how genotyping errors at low read depths bias MSC parameter estimates, following the approach of PMC12359030 [26].

1. Research Reagent Solutions

• Simulation Software: ms or msprime for generating gene trees and sequences under the MSC.
• Error Simulation Scripts: Custom scripts (e.g., in R or Python) to implement a Markov model of read depths and introduce base-calling errors.
• Inference Software: Bayesian inference software such as Bpp for estimating species trees and population parameters.

2. Methodology 1. Simulate True Sequences: Generate the "true" aligned sequences for multiple independent loci by simulating gene trees under the MSC model and then evolving DNA sequences along these trees. 2. Introduce Genotyping Errors: * Model Read Depth: For each site in each sequence, simulate a read depth using a Markov model that creates correlation between adjacent sites (e.g., a Beta-Markov model with bounds dmin=2 and dmax=100) [26]. * Simulate Reads and Call Genotypes: Given the true genotype and the simulated read depth at a site, simulate sequencing reads by multinomial sampling at a specified base-calling error rate (e). Call the genotype by maximum likelihood. 3. Parameter Inference: Analyze the resulting error-containing sequence alignments using Bpp (or similar), ignoring the introduced errors in the model. 4. Compare and Quantify Bias: Compare the parameter estimates (e.g., population sizes θ, divergence times τ) from the error-containing data against the known true values from the simulation. Quantify the bias and loss of accuracy.

3. Key Experimental Parameters The simulation should systematically vary key parameters to assess their influence, as shown in the table below.

Table 2: Experimental Parameters for Sequencing Error Simulation

Parameter	Description	Recommended Values for Testing
Average Read Depth (d̄)	The mean number of reads covering a genomic site.	3x, 5x, 10x, 20x
Base-calling Error Rate (e)	The probability of a single base being called incorrectly.	0.001 (Phred 30), 0.005, 0.01
Mutation Rate (μ)	The probability of a mutation per site per generation.	Set based on study organism (e.g., 1e-8)
Population Size (θ)	The scaled population size parameter (4Neμ).	Vary to create different ILS scenarios

Protocol 3: Inferring Parameters Under Gene Flow withmstree

This protocol details the use of the mstree program to estimate ancestral population sizes and divergence times in the presence of gene flow, without assuming a specific isolation-migration model [29].

1. Research Reagent Solutions

• Software: mstree (available at https://github.com/liujunfengtop/MStree/).
• Prerequisites: ms for generating simulated gene trees, seq-gen for converting trees to sequence data under a nucleotide substitution model (e.g., JC69).

2. Methodology 1. Input Data Preparation: Generate or obtain sequence alignments for a three-taxon system (two ingroup species and an outgroup) across tens of thousands of independent loci. 2. Gene Tree Estimation: Estimate gene trees for each locus using a standard phylogenetic method. 3. Run mstree: * The algorithm uses the distribution of coalescent times and gene tree topologies across loci. * It applies mathematical inequalities to estimate the species divergence time (τ₀) by identifying a point where gene tree shapes become equiprobable. * It then estimates the ancestral population size (θ₀) using the formula for the probability of deep coalescence. * Finally, it estimates the time of gene flow (τ₁) and the ancestral population size at that time (θ₁) using the frequency of specific gene tree topologies. 4. Validation: Compare the performance and speed of mstree against other methods like IMa3 or 3s on simulated data.

3. Visualization of the mstree Logical Workflow The core logic of mstree relies on the properties of gene tree distributions, as shown below.

The Scientist's Toolkit

Table 3: Essential Software and Analytical Tools

Tool Name	Type	Primary Function in MSC Research	Application Note
Bpp	Bayesian Inference Software	Infers species trees, divergence times, and demographic parameters under the MSC and MSC-with-introgression (MSC-I) models [26].	Used for assessing bias from sequencing errors; powerful for complex models but requires careful prior specification.
mstree	Parameter Estimation Program	Estimates ancestral population sizes and divergence times in the presence of gene flow without a pre-specified migration model [29].	Useful for testing different speciation scenarios (e.g., isolation-with-migration vs. secondary contact).
ms / msprime	Coalescent Simulator	Generates genomic data under user-specified demographic models, including the MSC [29] [48].	Essential for assessing method performance, power analysis, and studying model violations via simulation.
IQ-TREE	Phylogenetic Inference Software	Infers maximum likelihood gene trees and species trees from molecular sequences with robust branch support measures [21].	Critical for the gene tree estimation step; its accuracy directly influences GTEE.
tsinfer	Ancestral Recombination Graph (ARG) Estimator	Estimates local coalescent trees and ARGs from sequence data, integrating information from flanking regions [47].	A recombination-aware method that can improve allele frequency trajectory estimation for polygenic score analysis.
Neural Networks (MLP)	Machine Learning Method	Infers demographic parameters from a large set of summary statistics computed on genomic data [48].	Outperforms traditional ABC and other ML methods like Random Forest in complex demographic models [48].

Species delimitation, the process of determining boundaries between species, is a fundamental task in systematics and conservation biology. The multispecies coalescent (MSC) model has revolutionized this field by providing a statistical framework that accommodates genealogical fluctuations across genomes [49]. However, traditional MSC-based Bayesian methods, such as those implemented in the program BPP, have been widely documented to be prone to over-splitting – incorrectly identifying geographically separated populations as distinct species [49] [50] [33]. This limitation is particularly pronounced in cases of widespread taxa with continuous geographic variation or population structure, where these models may misinterpret population splits as species divergences [33] [6].

In response to these challenges, new heuristic approaches have been developed that incorporate the genealogical divergence index (gdi) within a more flexible framework. The Hierarchical Heuristic Species Delimitation (hhsd) pipeline implements two complementary algorithms – a merge algorithm and a split algorithm – to provide a range of plausible species hypotheses [49] [51]. This protocol details the methodology for implementing these approaches and provides a structured framework for interpreting their results, particularly when the merge and split algorithms yield conflicting conclusions.

Theoretical Foundation

The Genealogical Divergence Index (gdi)

The gdi is a quantitative measure that estimates lineage divergence by calculating the probability that two sequences from population A coalesce before either coalesces with a sequence from population B [49]. Ranging from 0 to 1, the gdi provides a standardized metric for assessing genetic isolation:

• gdi < 0.2: Populations typically considered a single species
• 0.2 ≤ gdi ≤ 0.7: Ambiguous species status
• gdi > 0.7: Populations typically considered distinct species [49]

The gdi is derived from population parameters (split times and effective population sizes) estimated from genomic data and can be calculated under models that account for gene flow, making it particularly valuable for delimiting species in complex scenarios [49].

The Merge and Split Algorithms

The hhsd pipeline implements two hierarchical algorithms that operate differently:

• The merge algorithm begins with all populations considered separate and progressively merges those with gdi values below a specified threshold. This approach tends to be more conservative, potentially lumping populations that might represent distinct species [6] [51].

• The split algorithm begins with all populations combined and progressively splits groups with gdi values above a specified threshold. This approach tends to be more liberal, potentially splitting populations that might belong to the same species [6] [51].

When these algorithms produce conflicting results, it indicates ambiguous species boundaries that require additional evidence for resolution [6].

Computational Workflow and Reagent Solutions

Research Toolkit

Table 1: Essential computational tools and resources for species delimitation analyses.

Tool/Resource	Primary Function	Application in Species Delimitation
HHSD Pipeline	Python-based implementation of hierarchical merge/split algorithms	Automates species delimitation based on gdi thresholds [49] [51]
BPP Software	Bayesian phylogenetics and phylogeography	Estimates population parameters under MSC; models gene flow [49] [33]
gdi Calculation	Quantitative lineage divergence index	Assesses genetic isolation between populations; ranges 0-1 [49] [33]
Genomic Data	Multi-locus sequence alignments	Input data for parameter estimation; preferably unlinked markers [50] [51]
Isolation-by-Distance Tests	Geographic analysis framework	Validates primary species hypotheses using geographic patterns [50]

Implementation Protocol

The following workflow outlines the step-by-step procedure for conducting hierarchical heuristic species delimitation:

Step 1: Data Preparation and Parameter Estimation

• Collect genomic sequence data from multiple unlinked loci across all study populations [50] [51]. Universal Single-Copy Orthologs (USCOs) are particularly effective as they provide genetically unlinked markers that represent a comprehensive sample of the genome [50].

• Estimate key population parameters – including population split times (τ), effective population sizes (θ), and migration rates (M) – using the BPP program under the MSC model with migration (MSC-M) [49] [51]. Proper prior selection is critical, as BPP analyses are sensitive to prior settings [33].

Step 2: gdi Calculation and Threshold Application

• Calculate pairwise gdi values between all populations using the estimated parameters [49].

• Apply the standardized gdi thresholds (0.2 and 0.7) to initially categorize population pairs [49].

Step 3: Hierarchical Algorithm Implementation

• Execute both the merge and split algorithms in the hhsd pipeline using the calculated gdi values [51].

• Use consistent gdi thresholds across both analyses to ensure comparability.

Step 4: Result Interpretation and Conflict Resolution

• Compare the results from both algorithms. Congruent results provide strong evidence for species boundaries [6].

• For discordant results, investigate potential causes including:

  • Gene flow between populations
  • Incomplete lineage sorting
  • Recent divergence events
  • Sampling artifacts or insufficient data [49] [6]

Step 5: Integrative Validation

• Apply isolation-by-distance tests to validate primary species hypotheses using geographic data [50].

• Examine contact zones between putative taxa for evidence of reproductive isolation or hybridization [6].

• Incorporate additional data sources including morphological, ecological, and behavioral evidence to reach a robust species hypothesis [50] [23].

Case Study Applications

Empirical Examples

Table 2: Application of merge/split algorithms to empirical datasets with interpretation guidelines.

Taxonomic Group	Merge Algorithm Result	Split Algorithm Result	Final Interpretation	Key Evidence
American Milksnakes (Lampropeltis triangulum complex)	Single species	Three species	Three species supported	Contact zones show reproductive isolation and sympatry between L. triangulum, L. elapsoides, and L. polyzona; other subspecies intergrade gradually [6]
Longear Sunfish (Lepomis megalotis)	Single species	Single species	Single species supported	No support for multiple species; consistent results across algorithms [6]
Southeast Asian Toads (Various Bufonidae)	Variable depending on priors	Variable depending on priors	Context-dependent	Careful prior selection based on empirical data crucial; gdi complements BPP analysis [33]
Pachyhynobius Salamanders	Multiple cryptic lineages	Multiple cryptic lineages	Four cryptic species supported	Deep genetic divergence, geographic disjunction, sky island effects [23]

Decision Framework for Conflicting Results

The following decision diagram provides a structured approach for resolving conflicts between merge and split algorithm results:

Discussion and Best Practices

Algorithm Selection Guidelines

The choice between relying more heavily on merge versus split algorithm results depends on several factors:

• For taxonomically conservative approaches with higher thresholds for species recognition, place greater weight on the split algorithm results, which require stronger evidence for separation [6].

• For discovering cryptic diversity in poorly studied groups, place greater weight on the merge algorithm results, which may reveal previously unrecognized lineages [23].

• When priors are poorly informed, the merge algorithm may provide more stable results, as BPP analyses are sensitive to prior selection [33].

Limitations and Considerations

Several important limitations must be considered when implementing these approaches:

• The gdi framework assumes a simple split model without gene flow, though extensions have been developed [49]. Gene flow between populations can complicate gdi interpretation and requires specific MSC-with-migration models [49].

• Paraphyletic species present challenges for hierarchical algorithms, which work best with strictly bifurcating phylogenies [49].

• Sampling design significantly impacts results; inadequate geographic sampling can lead to either over-splitting or lumping of species [50].

• The arbitrary nature of species boundaries in cases of allopatric divergence means that quantitative criteria like gdi provide guidelines rather than absolute rules [49].

Integration with Other Evidence

Successful species delimitation requires integrating multiple lines of evidence beyond genetic data:

• Morphological data: Assess diagnostic characters and phenotypic gaps [23]

• Ecological data: Evaluate niche differentiation and adaptive divergence [50]

• Behavioral data: Consider mating signals, reproductive isolation mechanisms [50]

• Geographical data: Analyze distribution patterns and contact zones [6]

The merge and split algorithms in the hhsd pipeline provide a valuable heuristic framework for exploring species boundaries, but their results should be interpreted as hypotheses to be tested with additional evidence rather than definitive conclusions [49] [6].

The hierarchical merge and split algorithms implemented in the hhsd pipeline represent significant advances in species delimitation methodology, directly addressing the over-splitting problem prevalent in earlier MSC-based approaches. By providing two complementary perspectives on species boundaries and incorporating the quantitative gdi metric, these tools offer a more nuanced framework for delimiting species in challenging cases involving recent divergence, gene flow, or continuous geographic variation.

The protocol outlined here emphasizes the importance of interpreting conflicting results as biological insights rather than methodological failures, and provides a structured approach for resolving these conflicts through additional lines of evidence. As genomic datasets become increasingly accessible, these approaches will play a crucial role in characterizing Earth's biodiversity and informing conservation priorities in an era of rapid environmental change.

The Impact of Substitution Model Fit on Coalescent Parameter Inference

In phylogenomics, the Multispecies Coalescent Model (MSC) provides a foundational framework for inferring species trees from genomic data by accounting for genealogical fluctuations across loci due to incomplete lineage sorting (ILS) [52]. However, the accuracy of parameter inference under the MSC—such as species divergence times, population sizes, and gene flow rates—is contingent upon the correct specification of the nucleotide substitution model. Model misspecification can introduce biases that propagate through the analysis, leading to incorrect topological inferences and biased parameter estimates [52] [26]. This application note details protocols for assessing substitution model fit and its impact on coalescent parameter inference, providing a structured approach for researchers in evolutionary biology and drug development, where understanding genetic relationships is crucial for target identification.

Quantitative Data on Model Misspecification Impact

Simulation studies are critical for quantifying the effects of substitution model misspecification on coalescent parameter inference. The following tables summarize key findings from in silico experiments.

Table 1: Impact of Genotyping Errors on MSC Parameter Estimation

Parameter	Low Error (e=0.001)	High Error (e=0.01) & Low Depth	Effect of Treating Heterozygotes as Missing
Species Tree Topology Accuracy	Little to no impact on power	Reduced power of estimation	Minor improvement in robustness
Population Size (θ)	Unbiased estimates	Biased estimates	Reduces bias in estimates
Divergence Time (τ)	Accurate estimation	Biased estimation	Improves estimation accuracy
Gene Flow Rate (M)	Reliable inference	Biased inference	Mitigates some bias

Table 2: Performance of Summary Statistics for Identifying Model Violations

Summary Statistic	Type	Detection Accuracy	Computational Demand	Key Utility
Mean FST	Data-based	Moderate (~70%)	Low	Measures population structure
FST Outlier Test	Data-based	High (~83%)	Low	Identifies lineages with gene flow
Robinson-Foulds Distance	Inference-based	High	High	Assesses topological differences
Tree Likelihood Mean/SD	Inference-based	Moderate	High	Evaluates model fit quality

Experimental Protocols

Protocol 1: Posterior Predictive Assessment of Model Fit

This protocol uses the P2C2M.SNAPP R package to perform posterior predictive checks for identifying MSC model violations [52].

Materials

• Input Data: SNP datasets in SNAPP .xml format
• Software Requirements: R statistical environment, P2C2M.SNAPP package, SNAPP phylogenetic software, BEAST2 package
• Computational Resources: High-performance computing cluster recommended for parallel processing

Procedure

• Initial SNAPP Analysis

  • Conduct species tree estimation using SNAPP with proper prior selection.
  • Ensure Markov chain convergence by running analyses for sufficient generations.
  • Retain at least 100 trees from the posterior distribution for sampling.

• Posterior Predictive Simulation

  • Sample species trees uniformly or randomly from the posterior distribution.
  • Extract taxonomic relationships and branch lengths from each tree.
  • Simulate posterior predictive datasets under the MSC model using fastsimcoal2.
  • Convert simulated datasets to SNAPP .xml format for analysis.

• Summary Statistics Calculation

  • Calculate data-based statistics (e.g., mean F_ST, F_ST range) from empirical and posterior predictive datasets.
  • Calculate inference-based statistics (e.g., Robinson-Foulds distance, tree likelihood metrics) from posterior distributions.
  • Compute posterior predictive p-values by comparing empirical summary statistics to the posterior predictive distribution.

• Interpretation

  • Identify model violations using a significance threshold of α = 0.05.
  • Use F_ST outlier tests to pinpoint specific lineages experiencing gene flow.

Protocol 2: Assessing Substitution Model Robustness with Simulation

This protocol evaluates the impact of substitution model misspecification on coalescent parameter inference using the Bpp software [26].

Materials

• Sequence Data: Multilocus genomic alignments
• Software: Bpp program for Bayesian phylogenetic inference
• Simulation Tools: Custom scripts for generating sequences under varying evolutionary models

Procedure

• Data Simulation

  • Generate true gene trees under the MSC model with known parameters.
  • Evolve sequences along gene-tree branches under a known substitution model (e.g., GTR+Γ).
  • "Post-process" sequences to introduce genotyping errors at specified base-calling error rates (e.g., e = 0.001 vs. e = 0.01).

• Inference Under Misspecified Models

  • Analyze simulated datasets using Bpp under both correct and incorrect substitution models.
  • Compare estimates of species trees, divergence times, population sizes, and gene flow rates to known true values.
  • Assess the impact of treating heterozygotes as missing data to mitigate genotyping error effects.

• Bias Quantification

  • Calculate deviation of parameter estimates from true values across multiple replicates.
  • Evaluate species tree topological accuracy using Robinson-Foulds distances.
  • Assess confidence interval coverage and statistical precision.

Workflow Visualization

MSC Model Validation Workflow

The Scientist's Toolkit

Table 3: Essential Research Reagents and Computational Tools

Resource	Type	Primary Function	Application Context
SNAPP	Software Package	Species tree estimation from SNP data	Bayesian inference under the MSC model for genomic-scale data [52]
P2C2M.SNAPP	R Package	Posterior predictive model checking	Identifying model violations in MSC analyses through simulation [52]
Bpp	Software Suite	Bayesian phylogenetic inference	Estimating species trees, divergence times, and gene flow under MSC models [26]
fastsimcoal2	Simulation Tool	Coalescent simulation	Generating genomic data under complex demographic scenarios [52]
SNP Datasets	Data Type	Empirical genetic variation	Input for species tree estimation and model validation studies [52]
FST Statistics	Analytical Metric	Population differentiation measure	Quantifying gene flow and identifying model violations [52]

Integrating Contact Zone Analysis to Validate Taxonomic Hypotheses

In the context of research utilizing the multispecies coalescent (MSC) model with gene flow, distinguishing true species boundaries from patterns of intraspecific geographic variation remains a significant challenge [53] [54]. The multispecies coalescent model provides a powerful framework for building phylogenetic trees from multilocus DNA sequence data, with more recent implementations, such as the multispecies network coalescent, explicitly allowing for gene flow among branches [54]. However, genomic analyses based on limited sampling can conflate interspecific and intraspecific variation, leading to inaccurate species delimitation [53]. Contact zones—geographic areas where divergent populations meet and potentially interbreed—provide critical natural experiments for testing taxonomic hypotheses and validating genomic inferences drawn from MSC analyses [55] [53]. By analyzing patterns of hybridization, introgression, and reproductive isolation in these zones, researchers can determine whether assortative mating or selection against hybrids exists (supporting distinct species) or if mating is random and admixture is gradual (supporting a single species) [53]. This protocol details the integration of contact zone analysis with MSC-based genomics to resolve taxonomic uncertainties in speciation continua.

Key Concepts and Definitions

The Multispecies Coalescent with Gene Flow

The multispecies coalescent (MSC) model is a fundamental framework for building phylogenetic trees from multilocus DNA sequence data, but its standard formulation assumes no gene flow between diverging species [54]. The multispecies network coalescent model extends the MSC by explicitly allowing for gene flow among branches, thus providing a more realistic representation of the speciation process for many taxa [54]. This is particularly important given that analysis of genomic data over the past two decades has highlighted the prevalence of introgression as an important evolutionary force in both plants and animals [42].

Contact Zones in Speciation Research

Contact zones are geographic regions where two or more divergent genetic forms of a species, or closely related species, come into contact and potentially hybridize [55] [53]. They provide crucial insights into the role of gene flow in adaptation and speciation. While previous work has often focused on contact zones involving only two divergent forms, in nature, many more than two populations may overlap simultaneously, creating complex patterns of introgression that are more complicated than often assumed [55].

Table 1: Genomic Signatures in Contact Zones and Their Evolutionary Interpretation

Genomic Signature	Interpretation in Contact Zone Context	Taxonomic Implication
Maintained Genomic Islands of Divergence	Indicates selection and/or hybrid incompatibilities despite hybridization [55]	Supports ongoing reproductive isolation and separate species status
Gradual Clinal Variation	Suggests continuous gene flow across a broad intergrade zone [53]	Supports intraspecific geographic variation within a single species
Sharp Genomic Cline Breaks	Indicates strong reproductive barriers and limited successful introgression [55]	Supports distinct species boundaries
Differential Introgression Patterns	Suggests varying selection pressures across the genome [55]	Reveals genomic regions involved in reproductive isolation

Application Notes: Integrating Contact Zone Analysis with MSC

Resolving Taxonomic Uncertainties

Thorough sampling and analysis of contact zones can determine whether we are observing assortative mating or selection against hybrids (supporting distinct species) versus random mating and gradual admixture (supporting geographic variation within a single species) [53]. For example, a reanalysis of the American milksnake (Lampropeltis triangulum) complex using thorough sampling across a ~700 km-wide transect revealed a gradual and continuous geographic cline across a broad intergrade zone, supporting a single species rather than the multiple species proposed in earlier studies with limited sampling near contact zones [53].

Analyzing Complex Multi-Population Contact Zones

Most models assume simple two-population contact, but natural systems can be more complex. A study of a three-way contact zone among divergent freshwater resident, saltwater resident, and saltwater migratory three-spined stickleback (Gasterosteus aculeatus) revealed more complicated patterns of introgression than typically assumed [55]. The research found evidence for hybridization (mostly between saltwater resident and saltwater migratory forms), yet pairwise genomic islands of divergence between all forms were maintained across the contact zone, suggesting reproductive isolation is maintained despite some hybridization [55].

Methodological Considerations for Gene Flow Detection

Different genomic analysis methods can produce conflicting interpretations of gene flow. A case study in Drosophila found that full-likelihood implementation of the multispecies coalescent model (using the Bayesian program bpp) had high power to detect gene flow, high accuracy in estimating gene flow rates, and could infer gene flow between sister lineages undetected by summary methods [42]. In contrast, summary methods had low power and produced biased estimates of introgression probability [42]. This highlights the importance of method selection when analyzing contact zone genomics data.

Experimental Protocol: Genomic Analysis of Contact Zones

Sampling Design and Data Collection

• Transect Sampling: Establish a transect across the suspected contact zone with systematic sampling locations. For the American milksnake study, a ~700 km-wide transect across Kansas and Missouri provided sufficient resolution to analyze clinal variation [53].
• Genomic Data Generation: Use reduced-representation sequencing (e.g., ddRAD-Seq) or whole-genome sequencing to generate genome-wide markers. The Puffinus shearwater study utilized double digest Restriction-Site Associated DNA sequencing (ddRAD-Seq) to genotype species and subspecies across their entire geographical range [56].
• Phenotypic Data Collection: Where applicable, collect morphological or ecological trait data to correlate with genomic patterns [53].

Genomic Data Analysis Workflow

Figure 1: Genomic analysis workflow for contact zone studies.

• Population Structure Analysis:

  • Use principal component analysis (PCA) and ADMIXTURE to identify genetic clusters.
  • Apply fast model-based estimation of ancestry in unrelated individuals (e.g., using ADMIXTURE) [55].

• Cline Fitting and Admixture Analysis:

  • Fit genomic clines to allele frequency data across the transect.
  • Calculate admi index to quantify the extent of admixture between populations [53].
  • Use Bayesian hybrid index and genomic cline estimation with R packages like gghybrid [55].

• Gene Flow Detection and Quantification:

  • Apply full-likelihood MSC methods with gene flow (e.g., bpp under MSC-I or MSC-M models) to estimate direction, timing, and strength of gene flow [42].
  • Compare results with summary methods (e.g., HyDe) but prioritize full-likelihood approaches when computational resources allow [42].

• Detection of Genomic Islands of Divergence:

  • Identify regions of elevated divergence (Fst outliers) between forms.
  • Examine whether divergent regions are maintained across the contact zone despite hybridization [55].
  • Annotate genes in divergent regions to identify potential targets of selection.

Interpretation Framework for Taxonomic Validation

Figure 2: Decision framework for taxonomic validation using contact zone data.

Research Reagent Solutions

Table 2: Essential Research Reagents and Computational Tools for Contact Zone Genomics

Reagent/Tool	Function	Application Notes
ddRAD-Seq	Reduced-representation genomic sequencing	Cost-effective for generating genome-wide markers across many samples [56]
bpp Software	Bayesian phylogenetic analysis under MSC	Full-likelihood inference of gene flow direction, timing, and strength; preferred over summary methods [42]
gghybrid R Package	Genomic cline analysis	Estimates hybrid indices and tests for deviations from neutral introgression [55]
ASTRAL/MP-EST	Species tree estimation from gene trees	Consistent in the "anomaly zone" but may apportion error differently than full-likelihood methods [54]
ADMIXTURE	Population structure analysis	Fast, model-based estimation of individual ancestries [55]

Case Study Applications

Three-Spined Stickleback Complex

A study of a three-way contact zone among three-spined stickleback ecotypes demonstrated the maintenance of reproductive isolation despite hybridization. Genomic analyses revealed pairwise islands of divergence between all forms that were maintained across the contact zone, with selection and/or hybrid incompatibilities in divergent regions involving many known adaptive loci and chromosomal inversions [55]. This complex pattern would not have been evident in a simple two-population contact model.

American Milksnake Complex

Reanalysis of the American milksnake complex with thorough sampling across a contact zone revealed a gradual and continuous geographic cline across a broad intergrade zone, supporting a single species (Lampropeltis triangulum) rather than the recognition of L. gentilis as a distinct species, as proposed in earlier studies with limited sampling near contact zones [53]. This highlights how contact zone analysis can correct taxonomic misinterpretations based on limited sampling.

Puffinus Shearwaters

Genomic data from ddRAD-Seq revealed that current taxonomies for Puffinus shearwaters were not supported, allowing researchers to propose a more accurate taxonomy by integrating genomic information with other sources of evidence. The study showed that several taxon pairs were at different stages of a speciation continuum, emphasizing the potential of genomic data to resolve taxonomic uncertainties [56].

Integrating contact zone analysis with the multispecies coalescent model with gene flow provides a powerful approach for validating taxonomic hypotheses in speciation continua. This integrated framework allows researchers to distinguish intraspecific geographic variation from interspecific divergence, identify genomic regions involved in reproductive isolation, and develop more accurate taxonomic classifications that reflect actual evolutionary relationships. As genomic methods continue to advance, particularly full-likelihood implementations of the MSC with gene flow, our ability to detect and quantify introgression while accounting for incomplete lineage sorting will further enhance the resolution of contact zone studies.

Empirical Validation and Model Comparison: MSC vs. Concatenation

The Multispecies Coalescent (MSC) model provides a foundational framework for inferring species relationships and demographic history from genetic data. However, real-world evolutionary processes often violate standard MSC assumptions, with gene flow between species being a particularly critical factor. The integration of gene flow into the MSC framework has led to the development of specialized models, primarily the MSC-with-Introgression (MSC-I) for discrete pulse events and the MSC-with-Migration (MSC-M) for continuous gene flow [14]. These models are highly idealized representations of complex biological realities. Consequently, rigorously comparing these models and testing their adequacy for a given dataset is not merely a statistical formality but an essential step toward producing reliable, biologically meaningful inferences regarding species divergence, ancestral population sizes, and patterns of gene flow [57].

Model misspecification can lead to severely biased results. For instance, incorrectly assigning gene flow to the wrong lineage within a species tree can cause large biases in estimated migration rates [14]. Furthermore, assuming a panmictic (unstructured) population history when deep ancestral structure exists can create artifactual signals of population size changes, such as false bottleneck events [58]. This protocol outlines a comprehensive statistical framework for comparing competing MSC models with gene flow and assessing model fit, ensuring that evolutionary inferences are robust and accurately reflect the underlying biological processes.

Theoretical Foundation of MSC Models with Gene Flow

Key Models and Their Applications

The two primary models for inferring gene flow within the MSC framework make different assumptions about the nature and timing of gene flow, making them suitable for different biological scenarios.

• MSC-with-Introgression (MSC-I): This model conceptualizes gene flow as a discrete, instantaneous event, such as a single hybridization or introgression pulse. It is parameterized by an introgression probability ((\phi)), which represents the proportion of genomes in one population that are migrants from another at the time of the event [14]. The MSC-I model is often applied to scenarios like hybrid speciation or rare admixture events between otherwise isolated species.

• MSC-with-Migration (MSC-M): This model assumes that gene flow is continuous and occurs at a constant rate over an extended period. It is parameterized by a population migration rate ((M)), which is the expected number of migrant individuals per generation moving from one population to another [14]. The MSC-M model is a generalization of the Isolation-with-Migration (IM) model and is suitable for inferring sustained gene flow between populations with ongoing contact.

Table 1: Comparison of MSC Models with Gene Flow

Feature	MSC-with-Introgression (MSC-I)	MSC-with-Migration (MSC-M)
Mode of Gene Flow	Discrete pulse(s)	Continuous migration
Key Parameter	Introgression probability ((\phi))	Population migration rate ((M))
Biological Scenario	Hybridization, rare admixture	Sustained contact, porous barriers
Identifiability Challenges	Can be hard to distinguish from MSC-M under certain conditions [14]	Misspecification of direction or lineage can cause large biases [14]

The Challenge of Model Misspecification

Several types of model misspecification can severely impact the accuracy of inference:

• Misspecified Mode of Gene Flow: Assuming the MSC-M model when the true history involved a discrete introgression pulse (MSC-I), or vice versa. Reassuringly, one study suggests that while this misspecification occurs, it may have relatively small local effects on inference, and gene flow can often still be detected with high power despite the misspecification [14].
• Incorrect Lineage Assignment: Assigning a gene flow event to an incorrect branch in the species tree (e.g., a parental branch instead of the correct daughter lineage). This type of error is common and can lead to large biases in estimated gene flow rates [14].
• Misspecified Direction of Gene Flow: Incorrectly identifying the source and recipient populations in a gene flow event. This can make it difficult to distinguish between scenarios like early divergence with gene flow and recent isolation [14].

Protocols for Model Comparison and Adequacy Testing

Workflow for Phylogenetic Analysis with Model Checking

A robust analytical pipeline incorporates model checking as a core component, not an afterthought. The following workflow is recommended when dealing with potential gene flow:

• Initial Species Tree Estimation: Infer a species tree under the standard MSC model (without gene flow) using software like BEAST 2 [59] or BPP.
• Test for Model Adequacy: Use a tool like P2C2M2 to perform a posterior predictive check on the initial MSC model fit. This test evaluates whether the observed data deviates significantly from the data simulated under the assumed model [57].
• Investigate Sources of Conflict: If the model is rejected, explore potential causes. P2C2M2 can help determine if the observed gene tree discordance is consistent with Incomplete Lineage Sorting (ILS) alone or if other processes like gene flow are likely [57].
• Incorporate Gene Flow: If gene flow is suspected, specify and test competing MSC-I and MSC-M models. Use Bayesian model comparison techniques, such as comparing marginal likelihoods (e.g., using stepping-stone sampling), to determine which model better explains the data.
• Validate the Selected Model: Once a model with gene flow is selected, re-run the model adequacy test (e.g., P2C2M2) to ensure the expanded model is no longer rejected.

Protocol 1: Testing MSC Model Adequacy with P2C2M2

The P2C2M2 R package provides a method for assessing the fit of the MSC model to a given dataset, which is a critical first step in identifying the need for more complex models involving gene flow [57].

• Objective: To determine whether the standard MSC model provides an adequate fit to the multi-locus sequence data, or if processes like gene flow are causing significant model violation.
• Experimental Principle: The method uses posterior predictive checks. It compares summary statistics (e.g., the distribution of gene tree topologies or node heights) calculated from the observed data against the distribution of the same statistics calculated from data simulated under the fitted MSC model. A significant discrepancy indicates model failure [57].
• Required Input:
  • A set of estimated gene trees from multiple loci (can be estimated from sequence alignments using software like BEAST 2 [59]).
  • The species tree estimate (with branch lengths in coalescent units).
  • Species delimitation model.
• Procedure:
  • Process Input Data: Format the gene trees and species tree into the required input for P2C2M2.
  • Simulate Posterior Predictive Data: Using the parameters of the estimated MSC model, simulate a large number of replicate gene tree datasets.
  • Calculate Discrepancy Statistics: For both the observed and simulated datasets, calculate summary statistics that are sensitive to model misspecification (e.g., tree likelihoods, tree distance metrics).
  • Compare Distributions: Statistically compare the distribution of the discrepancy statistics from the simulated data to the value from the observed data. A p-value is generated, and a small p-value (e.g., < 0.05) indicates that the observed data is unlikely under the model, suggesting model inadequacy [57].

Table 2: Key Software Tools for Model Comparison and Testing

Tool Name	Primary Function	Model/Method	Application Context
BEAST 2 [59]	Bayesian Phylogenetic Inference	MCMC under MSC & extensions	Core platform for estimating species trees and population parameters from molecular sequences.
P2C2M2 [57]	Model Adequacy Testing	Posterior Predictive Checks	Testing whether the standard MSC model fits the data, with failure suggesting gene flow or other processes.
cobraa [58]	Ancestral Structure Inference	Hidden Markov Model (HMM)	Inferring deep ancestral population structure and admixture from a single diploid genome.
BPP [14]	Bayesian Species Tree Inference	MSC-I, MSC-M	Co-estimating species trees and gene flow parameters (both discrete and continuous) from sequence data.

Protocol 2: Differentiating Deep Structure from Panmixia with cobraa

The cobraa method is designed to detect and characterize ancient population structure and admixture that cannot be identified by panmictic models like PSMC [58].

• Objective: To infer a structured ancestral population history, including population sizes, split times, admixture time, and admixture proportion, from a single diploid genome.
• Experimental Principle: cobraa is a coalescent-based Hidden Markov Model (HMM) that leverages not only the coalescence time profile but also the transition probabilities between neighboring coalescence times along the genome. This allows it to distinguish a structured model from an unstructured model, even when they have identical coalescence rate profiles [58].
• Required Input: A diploid genome sequence (e.g., a VCF file from a human individual).
• Procedure:
  • Data Preprocessing: Encode the genome into homozygous and heterozygous states for input into the HMM.
  • Model Parameterization: Specify a structured model with two ancestral populations (A and B) that diverge at time T2, admix at a more recent time T1 with fraction γ from B into A, and constant population size for B.
  • Parameter Estimation: Use an Expectation-Maximization (EM) algorithm to estimate the parameters NA(t), NB, and γ for a given pair of fixed times T1 and T2.
  • Grid Search: Perform multiple runs of cobraa over a grid of possible (T1, T2) time pairs. For each pair, run the EM algorithm to convergence.
  • Model Selection: Compare the log-likelihoods across all (T1, T2) pairs. The pair with the highest log-likelihood represents the best-supported split and admixture times, with the corresponding estimates for γ, NA(t), and NB [58].

Case Study: Resolving Phylogenetic Conflict in Pandanales

A recent study on the plant order Pandanales provides a compelling example of applying these principles to resolve long-standing phylogenetic conflicts [22].

• The Problem: Phylogenetic analyses of the five families within Pandanales produced strongly supported but topologically incongruent trees when using different methods (coalescent- vs. concatenation-based) and different genomic compartments (nuclear vs. plastid) [22].
• Model Testing and Gene Flow Analysis: The researchers assembled 2,668 single-copy orthologous genes. After establishing conflicting topologies, they used the HyDe software to test for signals of gene flow between lineages. Their analysis identified two significant ancient gene flow events: one between Velloziaceae and Triuridaceae, and another between Triuridaceae and the Cyclanthaceae-Pandanaceae (C-P) clade [22].
• Conclusion: The study concluded that gene flow, rather than ILS, was the primary source of phylogenetic conflict. The concatenation-based topology, which was more consistent with the signals of gene flow, was argued to likely reflect the true species tree. This case highlights that without explicitly testing for and incorporating gene flow, inferred phylogenies can be misleading [22].

The Scientist's Toolkit: Essential Research Reagents and Software

Table 3: Key Research Reagent Solutions for MSC Analyses with Gene Flow

Item Name	Type	Function/Biological Context	Example Use Case
Single-Copy Orthologous Genes (SCOGs)	Genomic Data	Phylogenetic markers used to infer species trees and detect discordance.	Resolving phylogenetic conflicts in Pandanales [22].
BEAST 2 XML Template	Software Configuration	Defines model hierarchy, priors, and MCMC settings for Bayesian analysis.	Setting up a complex MSC-I analysis in BEAUti [59].
P2C2M2 R Script	Analysis Script	Automates posterior predictive checks for MSC model adequacy.	Testing for model violation due to gene flow in a lizard phylogeny [57].
cobraa HMM Engine	Computational Algorithm	Core hidden Markov model for inferring ancestral structure from diploid data.	Estimating deep admixture proportions in human history [58].
1000 Genomes Project Data	Reference Dataset	Provides population-genetic polymorphism data for analysis and calibration.	Serving as input for cobraa to infer human ancestral structure [58].

The multispecies coalescent (MSC) model represents a fundamental paradigm shift in phylogenetic inference, providing a robust statistical framework for reconstructing species evolutionary histories from genomic data. The MSC extends single-population coalescent theory to multiple species, integrating the phylogenetic process of species divergences with the population genetic process of coalescence [13] [24]. This model has emerged as the natural framework for addressing key biological questions using genomic sequence data from multiple species, including estimation of species divergence times, population sizes, species trees accommodating discordant gene trees, and inference of cross-species gene flow [15].

A central challenge in phylogenomics has been the historical dominance of concatenation approaches, which combine all gene sequences into a single supermatrix for analysis. While computationally convenient, concatenation assumes that a single evolutionary history underlies all genes, neglecting inherent gene tree heterogeneity caused by incomplete lineage sorting (ILS) [3]. The MSC framework explicitly accounts for this heterogeneity, modeling how stochastic coalescence of gene lineages within evolving populations can cause individual gene trees to differ from the overall species tree [13] [24].

This article synthesizes empirical and simulation evidence demonstrating the superior performance of MSC methods over concatenation approaches, particularly in challenging evolutionary scenarios characterized by rapid speciation events, large effective population sizes, and gene flow. We provide detailed protocols for implementing MSC-based analyses and highlight essential computational tools that empower researchers to leverage this powerful framework for more accurate phylogenetic inference.

Theoretical Foundation of the Multispecies Coalescent

Basic Principles and Mathematical Framework

The multispecies coalescent model describes the genealogical relationships for samples of DNA sequences taken from multiple species, representing the application of coalescent theory to the case of multiple species [24]. The MSC model incorporates two sets of parameters: species divergence times (τ) and population size parameters (θ), both measured in the expected number of mutations per site [13]. For a sample of n sequences from a diploid population of size N, the coalescent process follows a stochastic pattern where lineages join backward in time, with coalescence times following exponential distributions [13].

Under the MSC, the probability distribution of gene trees and coalescent times provides the foundation for full-likelihood methods of species tree estimation [13]. The joint probability density of the gene tree topology and coalescent times for a sample of n sequences is given by:

[ f(G,{tj}) = \prod{j=n+1}^{m} \left[ \frac{2}{\theta} \exp\left(-\frac{j(j-1)}{\theta}tj\right) \right] \exp\left(-\frac{n(n-1)}{\theta}(\tau - (tm + t{m-1} + \cdots + t{n+1}))\right) ]

where G represents the gene tree topology, (t_j) are the coalescent times, and τ is the species divergence time [13]. This mathematical framework enables computation of the probabilities of gene tree topologies marginalizing over coalescent times, which is useful for heuristic species tree estimation methods [13].

Gene Tree-Species Tree Discordance

A key insight provided by the MSC model is that gene tree-species tree discordance is an expected biological phenomenon rather than mere statistical noise. In the simplest case of a rooted three-taxon tree, there are three possible species tree topologies but four distinct gene tree patterns when coalescent times are considered [24]. The probability of gene tree congruence with the species tree is:

[ P(\text{congruence}) = 1 - \frac{2}{3}\exp(-T) ]

where T is the branch length in coalescent units [24]. This equation illustrates that as internal branches shorten (e.g., in rapid radiations), discordance becomes more likely, highlighting the limitations of methods that assume a single underlying topology.

Figure 1: The MSC Framework Logic. The multispecies coalescent model integrates species phylogeny with population genetic processes to explain gene tree variation and enable more accurate species tree estimation.

Empirical Evidence: Quantitative Comparisons

Simulation Studies Demonstrating MSC Superiority

Simulation studies provide controlled environments for comparing phylogenetic methods by testing performance against known ground truth trees. Multiple studies have demonstrated the superior accuracy of MSC methods over concatenation, particularly in challenging phylogenetic regions characterized by short internal branches and rapid diversification events.

Table 1: Simulation-Based Evidence for MSC Superiority Over Concatenation

Study Context	Key Finding	Explanation	Implication
General Phylogenomic Simulations [3]	Concatenation produces "spuriously confident yet conflicting results" in regions of parameter space where MSC models perform well	Concatenation misestimates branch support when incomplete lineage sorting is present	MSC provides more accurate confidence assessments
Data Subsampling Tests [3]	Concatenation shows "erratic behavior" when subjected to data subsampling	Concatenation is sensitive to specific gene sampling, producing conflicting trees	MSC offers greater stability and robustness
Mammalian Phylogenomics [3]	Concatenation explained only 0-15% of observed gene tree heterogeneity	Most gene tree variation was attributed to factors other than the species tree	MSC accounts for multiple sources of gene tree variation

The erratic behavior of concatenation under data subsampling is particularly revealing. As noted in critical assessments, "simulations reveal the erratic behavior of concatenation when subjected to data subsampling and its tendency to produce spuriously confident yet conflicting results in regions of parameter space where MSC models still perform well" [3]. This demonstrates that concatenation can provide strongly supported but incorrect topologies, whereas MSC methods maintain accuracy across subsampling replicates.

Biological Case Studies

Empirical analyses across diverse taxonomic groups have reinforced the conclusions from simulation studies, demonstrating the practical advantages of MSC approaches in real-world phylogenetic inference.

Table 2: Empirical Case Studies Comparing MSC and Concatenation Performance

Taxonomic Group	Concatenation Result	MSC Result	Interpretation
Great Apes [29]	Potentially inconsistent topologies with high support	Stable, biologically plausible species trees	MSC accommodates real gene tree heterogeneity
Mammals [3]	Incongruent relationships with varying gene samples	Consistent species tree despite gene tree variation	MSC robustness to sampling effects
Plants [3]	Potentially misleading topologies due to alignment errors	More accurate relationships after error correction	MSC frameworks can better handle data imperfections

In mammalian phylogenomics, analyses have shown that "concatenation explained little or none (0–15%) of the observed gene tree heterogeneity" [3], suggesting that most of the variation among gene trees stems from biological processes other than the species phylogeny itself. The MSC framework accommodates this heterogeneity through its explicit modeling of coalescent processes, leading to more accurate species tree estimates.

MSC Protocols and Implementation

Experimental Design Considerations

Proper experimental design is crucial for successful MSC analysis. The following considerations should be addressed:

• Locus Selection: Ideal data consist of short, independently segregating genomic segments far apart in the genome to ensure independent coalescent histories [13]. While the term "gene" or "locus" is commonly used, non-coding DNA is often preferable, though exonic data have been successfully utilized [13].

• Sample Size: For each species, multiple individuals per species are beneficial but not always necessary. Single individuals can suffice if the focus is primarily on species relationships rather than population parameters.

• Sequence Requirements: Loci should ideally be non-recombining within segments while different loci are freely recombining [29]. Longer sequences provide more phylogenetic information but must balance against the risk of recombination within loci.

• Taxon Sampling: Denser sampling of closely related species helps resolve short internal branches where ILS is most problematic.

Computational Workflow for MSC Analysis

The typical workflow for MSC-based phylogenetic inference involves multiple stages of data processing and analysis:

Figure 2: Standard MSC Analysis Workflow. The process begins with raw sequence data, proceeds through locus selection and gene tree estimation, and culminates in MSC-based species tree inference.

Stage 1: Data Preparation and Locus Selection

• Genome Assembly and Annotation: Obtain high-quality genome assemblies or transcriptome data for all target species.
• Orthology Prediction: Identify orthologous loci across species using tools such as OrthoFinder [60], OMA [60], or custom BLAST pipelines.
• Alignment: Generate multiple sequence alignments for each locus using appropriate software (e.g., MAFFT, PRANK). Consider alignment uncertainty through approaches like GUIDANCE2 [61].
• Alignment Filtering: Remove poorly aligned regions or entire loci with insufficient phylogenetic information using tools like TrimAl.

Protocol Note: For large datasets, automated pipelines like LMAP_S can streamline these preprocessing steps by integrating multiple software tools and automating file management [62].

Stage 2: Gene Tree Estimation

• Model Selection: For each locus, determine the best-fitting substitution model using tools like ModelTest-NG or jModelTest.
• Tree Inference: Estimate gene trees for each locus using maximum likelihood methods (e.g., RAxML, IQ-TREE) or Bayesian approaches (e.g., MrBayes).
• Quality Assessment: Evaluate gene tree support through bootstrapping or posterior probabilities, considering filtering or excluding poorly supported trees.

Technical Tip: For large numbers of loci, consider efficient likelihood implementations such as FastTree for rapid approximation of gene trees, though with potential trade-offs in accuracy.

Stage 3: Species Tree Inference Under MSC

• Method Selection: Choose appropriate MSC method based on dataset size and biological question:
  • Full-Likelihood Methods: BPP (Bayesian) for smaller datasets with parameter estimation
  • Summary Methods: ASTRAL for large datasets focusing on topology
  • Gene Flow Integration: MSTree [29] or BPP with migration models
• Parameter Estimation: For full-likelihood methods, estimate species divergence times (τ) and population sizes (θ) using MCMC sampling.
• Convergence Assessment: Ensure MCMC chains have converged and mixed properly using diagnostic tools like Tracer.
• Species Delimitation (if needed): Integrate species discovery and estimation using programs like BPP with reversible-jump MCMC.

Advanced Protocol: Accounting for Gene Flow

The basic MSC model assumes no gene flow after speciation, but this assumption is frequently violated in empirical datasets. The following protocol extends the MSC framework to incorporate cross-species gene flow:

• Initial Screening: Use phylogenetic network approaches (e.g., PhyloNet, SNaQ) to detect potential gene flow events.
• Model Selection: Compare isolation-only MSC models with those incorporating migration using Bayes factors or AIC-based approaches.
• Parameter Estimation under Gene Flow: Implement models that allow for migration between sister lineages, such as the isolation-with-migration (IM) model [29].
• Validation: Use posterior predictive simulations to assess model fit and identify potential model violations.

Case Example: The MSTree package implements an MSC approach that "does not rely on any assumption about the relationship between isolation and migration when estimating parameters" [29], making it particularly useful for complex speciation scenarios with uncertain gene flow histories.

Table 3: Computational Tools for MSC Analysis

Tool	Function	Method Type	Use Case
BPP [13]	Bayesian species tree estimation	Full-likelihood	Divergence time and population size estimation
ASTRAL	Species tree from gene trees	Summary method	Large datasets, topology focus
MSTree [29]	MSC with gene flow	Full-likelihood	Speciation with gene flow
PhyloNet	Phylogenetic networks	Network inference	Reticulate evolution detection
PAUP*	Phylogenetic analysis	Concatenation	Comparative analyses
RAxML/ IQ-TREE	Gene tree estimation	Maximum likelihood	Locus-level tree building
OrthoFinder [60]	Orthogroup inference	Gene family clustering	Locus identification

Table 4: Key Datasets and Standards for MSC Research

Resource	Description	Application
OrthoDB [60]	Hierarchical catalog of orthologs	Locus set selection
Angiosperms353 [63]	Universal probe set for plants	Targeted sequencing
Unicore [60]	Structural core gene identification	Deep phylogeny beyond sequence twilight zone
ETE Toolkit [62]	Python environment for tree analysis	Tree manipulation and visualization

Emerging methods like Unicore demonstrate how structural information can enhance MSC analyses, particularly for deep phylogenies where sequence similarity is low. Unicore "leverages Foldseek's rapid structural comparisons" to "identify single-copy 'structural core genes' on the fly without relying on precomputed gene sets" [60], potentially overcoming limitations of sequence-based orthology detection in the "twilight zone" of amino acid identity.

The empirical evidence comprehensively demonstrates that the multispecies coalescent model outperforms concatenation approaches across a wide range of phylogenetic scenarios, particularly those characterized by incomplete lineage sorting, rapid radiations, and gene flow. The theoretical robustness of the MSC framework, which explicitly models the stochastic nature of gene lineage coalescence, provides more accurate species tree estimates and more realistic assessments of uncertainty compared to concatenation methods.

While MSC methods require greater computational resources and more careful model specification, their biological realism and statistical performance justify these investments. Emerging approaches that integrate gene flow, structural information, and more efficient algorithms will further expand the utility of MSC methods across diverse phylogenetic contexts.

As phylogenomics continues to grapple with the complexities of evolutionary history, the multispecies coalescent provides an essential framework for reconciling gene tree discordance and reconstructing accurate species relationships, ultimately enabling more reliable inferences about the patterns and processes of biodiversity evolution.

Posterior Predictive Simulation (PPS) is a powerful model-checking technique that assesses how well a Bayesian model, including the Multispecies Coalescent Model (MSC), fits the observed data. The core principle is straightforward: if a model provides a good fit to the data, then it should be possible to generate replicated datasets from the fitted model that resemble the original observed data [64]. In the context of phylogenetics and the Multispecies Coalescent, this technique becomes crucial for identifying model violations that could lead to biased parameter estimates and incorrect phylogenetic inferences [65].

The multispecies coalescent model assumes that all incongruence among gene trees results solely from incomplete lineage sorting. However, real genomic data often violate this assumption due to factors like gene flow, introgression, and hybridization [66] [6]. When the MSC model is applied to datasets containing such violations, estimates of species trees, divergence times, and population sizes can become severely biased [66] [65]. The P2C2M.SNAPP R package was specifically developed to address this challenge by implementing posterior predictive checks for phylogenetic estimates obtained through the popular SNAPP program [65]. This methodology allows researchers to quantify model fit and identify when MSC model assumptions are violated, thus guarding against over-splitting of species and other erroneous conclusions [6].

Theoretical Foundation and Implementation

The Posterior Predictive Distribution

The mathematical foundation of posterior predictive checking centers on the posterior predictive distribution. For a model with parameters θ and observed data y, the posterior predictive distribution for new replicated data yrep is given by:

p(yrep | y) = ∫ p(yrep | θ) p(θ | y) dθ

In practice, we generate these replicates by first drawing parameter values θs from the posterior distribution p(θ | y), and then simulating new datasets yreps from the model p(yrep | θs). This process results in an S × N matrix of simulations, where S is the number of posterior draws and N is the number of data points [64]. Each row represents a full replicated dataset that can be compared to the original observations.

Workflow for Phylogenetic Model Checking

The implementation of posterior predictive checks for MSC models follows a systematic workflow. The P2C2M.SNAPP package leverages the posterior distribution of species trees output by SNAPP to simulate posterior predictive datasets under the MSC model [65]. Summary statistics are then computed for both the empirical data and the posterior predictive distribution, enabling quantitative comparison to identify model violations.

Table 1: Key Components of Posterior Predictive Simulation for MSC Models

Component	Description	Implementation in P2C2M.SNAPP
Model Input	Posterior distribution of species trees	SNAPP output files
Data Simulation	Generation of replicated datasets under MSC	Coalescent simulation using posterior tree parameters
Summary Statistics	Quantities capturing important data features	Tree-based statistics, site pattern frequencies
Comparison Metric	Measure of discrepancy between observed and simulated data	Statistical tests comparing empirical vs. predictive distributions
Violation Detection	Identification of significant model misfit	Classification of datasets based on summary statistic thresholds

In simulation testing, this approach has demonstrated high accuracy, correctly classifying up to 83% of datasets regarding whether they violate the MSC model assumptions, though the exact performance depends on the specific summary statistic employed [65].

Experimental Protocol: Implementing PPS with P2C2M.SNAPP

Software Requirements and Installation

To implement posterior predictive checks for MSC models, researchers should begin by establishing the necessary computational environment. The following software components are required:

• R Statistical Platform: Version 4.0.0 or higher
• P2C2M.SNAPP R Package: Available through standard R package installation methods
• SNAPP Program: Part of the BEAST 2 phylogenetic software suite
• Java Runtime Environment: Required for SNAPP execution

Installation of the P2C2M.SNAPP package can be accomplished using the following R code:

Additional program details and tutorials are available through the package documentation and associated publications [65].

Step-by-Step Implementation Protocol

• Model Fitting with SNAPP: Begin by running phylogenetic analysis using SNAPP to obtain the posterior distribution of species trees from your empirical sequence data. Ensure proper MCMC convergence diagnostics are performed and adequate effective sample sizes (ESS > 200) are achieved for all parameters of interest.

• Posterior Predictive Simulation: Using the SNAPP output, simulate posterior predictive datasets under the MSC model. The P2C2M.SNAPP package automates this process, generating replicated datasets that reflect both the model and the uncertainty in parameter estimates.

• Calculate Summary Statistics: Compute appropriate summary statistics for both the empirical data and the posterior predictive datasets. Effective summary statistics for MSC model assessment may include:

  • Tree shape statistics
  • Site pattern frequencies
  • Measures of topological discordance
  • Branch length distributions

• Compare Distributions: Quantitatively compare the distribution of summary statistics between the empirical data and posterior predictive replicates. Significant discrepancies indicate potential model violations.

• Interpret Results: Identify the specific nature of any model violations detected. For example, excess discordance beyond that expected under the MSC may indicate gene flow or introgression [66].

Validation and Diagnostics

The validation of posterior predictive checks involves several diagnostic considerations:

• Convergence Assessment: Ensure the underlying SNAPP analysis has properly converged using standard MCMC diagnostics (R-hat < 1.05, ESS > 200) [67].
• Summary Sensitivity: Evaluate whether results are robust to the choice of summary statistics by testing multiple different statistics.
• Calibration Assessment: Verify that the posterior predictive p-values are approximately uniformly distributed under correct model specification.

Table 2: Essential Research Reagents and Computational Tools

Reagent/Tool	Function	Application Note
SNAPP Program	Bayesian species tree inference under MSC	Generates posterior distribution of species trees from SNP data
P2C2M.SNAPP R Package	Posterior predictive checking for MSC models	Implements simulation-based model checking with user-friendly interface
Genomic Sequence Data	Empirical genetic sequences for analysis	SNP data preferred for SNAPP analysis; multiple unlinked loci required
Summary Statistics	Quantifiable features of genetic data	Measures of tree discordance, site patterns, or branch lengths
BEAST 2 Platform	Phylogenetic analysis framework	Provides environment for SNAPP execution and related analyses

Visualization of the Posterior Predictive Checking Workflow

The following diagram illustrates the complete workflow for assessing MSC model fit using posterior predictive simulation:

Diagram 1: Workflow for Posterior Predictive Checking of MSC Models

Advanced Applications in Models with Gene Flow

Addressing Model Misspecification

The standard MSC model becomes misspecified when gene flow occurs between species or populations. Recent methodological advances have developed approaches that explicitly incorporate gene flow into the coalescent framework, including the multispecies coalescent with introgression (MSC-I) and the multispecies coalescent with migration (MSC-M) models [66]. These models differ in their assumptions about the timing and pattern of gene flow:

• MSC-I Model: Models gene flow as short bursts of introgression at specific times, measured using an introgression probability.
• MSC-M Model: Assumes continuous gene flow over extended periods, quantified through a migration rate parameter.

Posterior predictive simulation can help determine which model better fits a given dataset and identify when simpler MSC models are insufficient. Research has shown that incorrectly specifying the mode of gene flow (e.g., assuming continuous migration when pulse introgression occurred) can lead to substantial biases in parameter estimates, including species divergence times and population sizes [66].

Decision Framework for Model Selection

When analyzing empirical data under potential gene flow scenarios, the following decision framework is recommended:

• Initial Assessment: Begin with standard MSC model and perform posterior predictive checks to assess baseline fit.
• Detect Violations: Use discrepancy measures between empirical and posterior predictive distributions to identify potential gene flow.
• Model Expansion: If violations are detected, fit more complex models incorporating gene flow (MSC-I or MSC-M).
• Model Comparison: Use Bayesian model selection techniques, such as Bayes factors, to compare alternative models.
• Sensitivity Analysis: Evaluate the robustness of conclusions to model assumptions and prior specifications.

The pulse introgression model (MSC-I) appears to be more robust to model misspecification and may be preferable in real data analysis unless there is substantive evidence for continuous gene flow [66].

Case Studies and Biological Applications

Empirical Examples from Systematic Biology

Posterior predictive checking has been successfully applied to resolve contentious taxonomic debates where gene flow complicates species delimitation. Key case studies include:

• American Milksnakes (Lampropeltis triangulum complex): Application of model-checking techniques helped distinguish between competing taxonomic hypotheses, supporting a three-species classification rather than excessive splitting into seven species or lumping into a single species [6].

• Longear Sunfish (Lepomis megalotis): Analyses incorporating gene flow models found no support for dividing this widespread species into multiple regional species, demonstrating how proper model assessment can prevent over-splitting of geographically variable populations [6].

These case studies highlight the practical importance of model checking in avoiding both over-splitting of geographically widespread taxa and under-estimation of true species diversity.

Integration with Bayesian Workflow

Posterior predictive checking represents one essential component of a comprehensive Bayesian workflow for phylogenetic analysis [68] [69]. This workflow includes iterative cycles of model building, inference, model checking and improvement, and model comparison. Within this framework, the purpose of model comparison extends beyond simply selecting the "best" model to developing a deeper understanding of how different models represent the underlying biological processes.

Effective Bayesian workflow emphasizes the importance of defining clear research questions at the outset, as these questions determine what data should be collected, what models are appropriate, and how results should be interpreted and communicated [68]. For species delimitation studies, this means framing questions about species boundaries in ways that can be addressed through statistical model comparison and validation.

Technical Considerations and Limitations

The effectiveness of posterior predictive checks depends heavily on the choice of summary statistics. These statistics should capture features of the data that are biologically relevant and sensitive to potential model violations. For MSC models with gene flow, effective summary statistics include:

• Measures of genealogical discordance between loci
• Patterns of shared genetic variation between putative species
• Distributions of branch lengths across gene trees
• Site pattern frequencies that reflect ancestral polymorphism or introgression

The development of novel summary statistics that are specifically powerful for detecting gene flow remains an active research area.

Computational Considerations

Posterior predictive simulation for MSC models is computationally intensive, requiring:

• Substantial memory resources for handling genomic datasets
• Efficient algorithms for simulating coalescent processes with migration
• Parallel computing infrastructure for analyzing multiple predictive replicates

Despite these challenges, the P2C2M.SNAPP package provides a user-friendly implementation that makes these methods accessible to practicing biologists [65].

Posterior predictive simulation provides an essential toolkit for assessing model fit in multispecies coalescent analyses, particularly when gene flow may complicate phylogenetic inference. The P2C2M.SNAPP package offers a practical implementation that enables researchers to detect model violations and make more informed decisions about species boundaries and evolutionary history.

As genomic datasets continue to grow in size and complexity, the importance of rigorous model checking will only increase. Future methodological developments will likely focus on more efficient computational approaches, more powerful summary statistics, and more integrated models that simultaneously account for multiple evolutionary processes. By adopting these model-checking practices, researchers can avoid common pitfalls in species delimitation and produce more reliable estimates of evolutionary relationships.

The multispecies coalescent (MSC) model provides a powerful framework for inferring species phylogenies by integrating the phylogenetic process of species divergences with the population genetic process of coalescence [15]. This model has become a fundamental paradigm in phylogenomics, enabling researchers to estimate species divergence times, population sizes, and species trees while accommodating discordant gene trees [15]. The MSC extends the single-population coalescent model to multiple species, offering a robust statistical approach for analyzing genomic sequence data from multiple species [15].

However, the application of MSC models faces significant challenges, particularly regarding their performance with geographically widespread taxa and their handling of cross-species gene flow. This article explores these challenges through two detailed case studies from vertebrate taxonomy: the American milksnakes (Lampropeltis triangulum complex) and the ocean sunfishes (family Molidae). These cases illustrate both the limitations and appropriate applications of MSC in resolving contentious taxonomic questions, providing valuable insights for researchers working on species delimitation and gene flow analysis.

The Multispecies Coalescent Model: Theoretical Framework

Core Principles and Methodological Approaches

The multispecies coalescent model operates on several key principles that distinguish it from concatenation approaches. It integrates the phylogenetic process of species divergences with the population genetic process of coalescent, providing a mathematical framework for understanding gene tree heterogeneity [15]. The MSC models the genealogical history of genes within the phylogenetic history of species, explicitly accounting for ancestral population sizes and deep coalescence [3].

Methodologically, MSC implementations can be categorized into full-likelihood and heuristic approaches. Full-likelihood methods, such as Bayesian Phylogenetics and Phylogeography (BPP), utilize Markov chain Monte Carlo algorithms to co-estimate species trees and parameters directly from sequence data [15]. Heuristic methods typically operate in two steps: first estimating gene trees from individual loci, then summarizing these trees into a species tree [15]. Recent methodological advances have extended the MSC to incorporate inference of cross-species gene flow, allowing for more realistic modeling of evolutionary history [15].

Limitations and Assumptions

The MSC model carries several critical assumptions that can impact its performance in empirical applications. The model assumes neutral evolution without selection, no gene flow or hybridization after speciation, and free recombination between loci but no recombination within loci [3]. In practice, these assumptions are frequently violated, particularly in complexes of recently diverged or parapatrically distributed species.

As Chambers et al. demonstrated, the MSC model has a particular tendency to over-split species in the case of geographically widespread taxa [70]. This over-splitting occurs because the model often interprets population structure as species-level divergences, especially when sampling is limited or does not adequately represent continuous geographic clines. Additionally, the use of transcriptomic data presents challenges for MSC implementation, as the constituent exons often span large chromosomal regions with substantial recombination [3].

Table 1: Key Features and Challenges of the Multispecies Coalescent Model

Feature	Description	Challenges and Limitations
Theoretical Basis	Extends single-population coalescent to multiple species; integrates phylogenetic and population genetic processes	Assumes neutral evolution, which may not hold for all genomic regions
Gene Tree Handling	Accommodates discordance between gene trees and species tree	Can be misled by other sources of gene tree heterogeneity (e.g., hybridization)
Parameter Estimation	Estimates species divergence times, population sizes, and cross-species gene flow	Requires substantial computational resources for large datasets
Data Requirements	Uses multi-locus genomic sequence data	Performance affected by limited sampling or continuous geographic variation
Species Delimitation	Provides statistical framework for testing species boundaries	Tendency to over-split geographically widespread taxa

Case Study 1: The American Milksnakes (Lampropeltis triangulumComplex)

Taxonomic History and Initial Splitting

The American milksnakes represent a classic example of the taxonomic controversy that can emerge from uncritical application of the MSC model. Historically classified as a single species, Lampropeltis triangulum, with numerous subspecies, this complex was substantially revised by Ruane et al. (2014), who proposed splitting it into seven separate species based on MSC analysis [71] [70]. This revision was quickly adopted by major taxonomic authorities, including the Society for the Study of Amphibians and Reptiles (SSAR) and the Reptile Database, significantly impacting species classifications on citizen science platforms like iNaturalist [71].

The Ruane et al. study utilized the Bayesian Phylogenetics and Phylogeography (BPP) program, a popular implementation of the MSC model, to analyze genetic data from the milksnake complex. Their analysis identified seven putative species within what was previously considered a single biologically variable species [70]. This splitting approach aligned with the trend in modern taxonomy to resolve paraphyly by splitting rather than lumping taxa, though this preference lacks consistent scientific justification [71].

Critical Reassessment and Evidence of Over-Splitting

Subsequent reevaluation of the Ruane et al. study by Chambers et al. revealed significant methodological limitations and sampling deficiencies. Their critical analysis, published in Systematic Biology, demonstrated that the MSC model had likely over-split the milksnake complex due to several factors: inadequate sampling design, failure to account for isolation by distance, and misinterpretation of continuous geographic clines as species boundaries [70].

Chambers et al. found that several of the hypothesized species represented "arbitrary slices of continuous geographic clines" rather than biologically distinct entities [70]. This was particularly evident along the Rio Grande, where snakes on opposite sides of the river were classified as different species despite inhabiting identical habitats and showing evidence of gene flow [71]. Additional re-analyses of the Ruane et al. dataset presented in public talks (though not yet published) reportedly showed substantial gene flow among the proposed species, further undermining the splitting hypothesis [71].

The most robust available evidence supports the recognition of only three species within the complex, rather than the seven proposed by Ruane et al. [70]. This case illustrates the critical importance of examining putative contact zones among delimited species and thoroughly analyzing geographic variation before making taxonomic changes based solely on MSC outputs [70].

Table 2: Comparative Analysis of Milksnake Taxonomic Hypotheses

Taxonomic Concept	Number of Species	Key Evidence	Major Limitations
Traditional View	1 species with multiple subspecies	Morphological consistency, ecological similarities	Does not account for genetic structure
Ruane et al. (2014) MSC Analysis	7 species	Bayesian Phylogenetics and Phylogeography (BPP) analysis	Poor sampling, ignored gene flow, interpreted clinal variation as species boundaries
Chambers et al. (2020) Reassessment	3 species	Analysis of geographic variation, contact zones, and evidence of gene flow	Limited incorporation of morphological and ecological data

Experimental Protocol: MSC Analysis for Snake Taxa

Protocol Title: Multispecies Coalescent Analysis for Delimiting Closely-Related Snake Species

1. Sample Collection and Preparation

• Conduct geographically comprehensive sampling across the entire distribution range, with particular attention to putative contact zones and areas of potential clinal variation.
• Collect tissue samples (scale clips, blood, or muscle) from voucher specimens with precise locality data.
• Extract genomic DNA using standard silica-column or phenol-chloroform methods, with quantification via fluorometry.

2. Sequence Data Generation

• Target approximately 1000-5000 ultraconserved elements (UCEs) or sequence capture of hundreds of nuclear loci via next-generation sequencing.
• Supplement with traditional Sanger sequencing of 5-10 mitochondrial genes (e.g., COI, ND4, Cyt-b).
• Ensure sequence quality through rigorous trimming, alignment validation, and contamination checks.

3. Data Processing and Alignment

• Process raw sequencing data using bioinformatic pipelines such as PHYLUCE for UCE data or customized workflows for targeted sequences.
• Align sequences using MAFFT or MUSCLE, with manual verification of alignments in problematic regions.
• Partition data by locus and apply appropriate nucleotide substitution models determined by ModelTest-NG.

4. Species Tree Estimation

• Implement multispecies coalescent analysis using BPP version 4.0 with guided species delimitation (A11 analysis).
• Run reversible-jump Markov Chain Monte Carlo (rjMCMC) for 100,000 generations, sampling every 100 generations, with at least two independent runs to ensure convergence.
• Apply alternative MSC methods such as StarBEAST2 or SNAPP for comparison with BPP results.

5. Gene Flow Assessment

• Test for cross-species gene flow using the MSC-with-migration model in BPP or specialized tools like D-statistics and TreeMix.
• Analyze geographic clines in allele frequencies and morphological traits using R packages cline and hierfstat.
• Conduct demographic modeling with ∂a∂i or FastSimCoal2 to estimate divergence times and migration rates.

6. Validation and Synthesis

• Compare MSC results with concatenated phylogenetic analyses using IQ-TREE or RAxML.
• Integrate morphological, ecological, and distributional data to assess biological significance of proposed species boundaries.
• Apply model validation techniques including posterior predictive checks and analysis of prior sensitivity.

Case Study 2: Ocean Sunfishes (Family Molidae)

Taxonomic Complexity and Conservation Challenges

The ocean sunfishes (family Molidae) present a contrasting case where taxonomic revision has been more consistently supported by multiple lines of evidence. The Molidae currently contains five recognized species within three genera: Mola mola, Mola alexandrini, Mola tecta, Masturus lanceolatus, and Ranzania laevis [72]. These species include the largest and heaviest bony fishes in the world, distributed from southern Chile to the Arctic Circle [72].

Unlike the milksnake case, the sunfish taxonomic revisions have been driven by a combination of morphological, genetic, and ecological evidence. Recent discoveries include Mola tecta (the hoodwinker sunfish), which was described in 2017 and had previously been mistaken for other Mola species [72]. Additionally, research suggests there may be up to three additional cryptic species of Mola mola that are genetically distinct though morphologically similar [73]. These findings highlight the complex speciation patterns in marine environments and the value of genetic approaches for identifying cryptic diversity.

Conservation Implications and Management Needs

The taxonomic clarity achieved for sunfishes has direct conservation implications. The IUCN Molidae review panel was formed to reassess Red List classifications, with three species previously ranging from Least Concern to Vulnerable, and two species (Mola alexandrini and Mola tecta) requiring initial assessment [72]. This reassessment was urgently needed given serious conservation concerns including vulnerability to mass bycatch, knowledge gaps of basic life history traits, and the lack of any conservation or management plans across their broad ranges [72].

Of particular concern are the unregulated regional fisheries targeting M. mola, M. alexandrini, and Masturus lanceolatus primarily across Japan, Taiwan, and South Korea [72]. These fisheries pose significant threats to sunfish populations, yet landings patterns and population impacts are challenging to assess due to inconsistent catch recording and the difficulty of studying these elusive oceanic species [72]. The IUCN panel emphasized that future data collection, research efforts, and conservation measures should be driven by a species-specific focus, avoiding the common practice of grouping all Molidae species under general labels like 'ocean sunfish' or 'Mola' [72].

Table 3: Conservation Status and Threats to Ocean Sunfish Species

Species	Common Name	IUCN Status	Major Threats	Fisheries Impact
Mola mola	Ocean Sunfish	Previously assessed (Vulnerable to Least Concern)	Bycatch, targeted fisheries, climate change	Targeted in unregulated regional fisheries
Mola alexandrini	Giant Sunfish	Awaiting first assessment	Bycatch, targeted fisheries, climate change	Targeted in unregulated regional fisheries
Mola tecta	Hoodwinker Sunfish	Awaiting first assessment	Bycatch, climate change	Limited information
Masturus lanceolatus	Sharptail Sunfish	Previously assessed (Vulnerable to Least Concern)	Bycatch, targeted fisheries, climate change	Targeted in unregulated regional fisheries
Ranzania laevis	Slender Sunfish	Previously assessed (Vulnerable to Least Concern)	Bycatch, climate change	Limited information

Experimental Protocol: Marine Species Delimitation

Protocol Title: Integrated Taxonomic Protocol for Marine Fishes

1. Sample Acquisition and Preservation

• Obtain specimens from diverse geographic locations across the species' distribution range through scientific surveys, fisheries bycatch, and stranded individuals.
• Collect tissue samples (fin clips, muscle, skin) preserved in 95% ethanol or RNAlater for genetic analyses.
• Document each specimen photographically, record morphometric measurements (standard length, total length, weight), and preserve voucher specimens when possible.

2. Genetic Data Generation

• Extract high-molecular-weight DNA using magnetic bead-based protocols optimized for degraded tissues.
• Prepare whole-genome sequencing libraries with 350-550bp insert sizes and sequence on Illumina platforms to achieve 30x coverage.
• Alternatively, target sequence capture of ultraconserved elements or mitochondrial genome sequencing for broader sampling.

3. Morphological Analysis

• Conduct detailed meristic counts (fin rays, gill rakers) and morphometric measurements using digital calipers.
• Perform geometric morphometrics using landmark-based analysis of body shape.
• Examine osteological features through CT scanning or clearing and staining of representative specimens.

4. Multispecies Coalescent Analysis

• Process raw sequencing data through quality filtering, adapter trimming, and contamination screening.
• Assemble mitochondrial genomes and call nuclear SNPs using reference-based and de novo approaches.
• Implement MSC analysis using SNAPP or StarBEAST2 with appropriate mutation rate priors and chain lengths sufficient for convergence (ESS > 200).

5. Ecological Niche Modeling

• Compile occurrence records from museum databases, literature, and citizen science platforms.
• Extract environmental data (temperature, salinity, chlorophyll-a) from Bio-ORACLE and MARSPEC databases.
• Model species distributions using MaxEnt with proper spatial thinning to reduce sampling bias.

6. Conservation Assessment

• Apply IUCN Red List criteria using the results from integrated taxonomic analysis.
• Estimate population trends from time-series data where available.
• Identify critical habitats and potential threats for conservation planning.

Comparative Analysis and Synthesis

Methodological Lessons and Best Practices

The contrasting outcomes in the milksnake and sunfish case studies provide valuable insights for the application of the MSC model in taxonomic research. The milksnake example illustrates the risks of uncritical application of MSC methods, particularly when sampling is inadequate and population structure is misinterpreted as species boundaries [70]. The sunfish case demonstrates the value of integrating multiple lines of evidence, including morphology, ecology, and distribution data, to validate genetic findings [72].

A critical lesson from these cases is that MSC-based species delimitation should not rely exclusively on genetic data but should incorporate thorough analyses of geographic variation and careful examination of putative contact zones [70]. This is particularly important for widespread taxa with continuous distributions, where the MSC has demonstrated a tendency toward over-splitting [70]. Researchers must be especially cautious when dealing with well-studied, geographically variable, and parapatrically distributed species complexes [70].

Additionally, these cases highlight the importance of considering the practical implications of taxonomic decisions. For sunfishes, species-specific focus is essential for effective conservation management, as different species face varying threats and population trajectories [72]. For milksnakes, excessive splitting may create artificial conservation units that do not reflect biological reality, potentially diverting resources from genuinely threatened populations.

Diagram: MSC-Based Species Delimitation Workflow

Research Reagent Solutions for MSC Studies

Table 4: Essential Research Reagents and Tools for Multispecies Coalescent Studies

Reagent/Tool	Application	Function in MSC Studies
BPP Software	Bayesian species delimitation	Implements MSC model for joint estimation of species trees and delimitation
* Tissue Preservation Buffer*	Field sample collection	Preserves DNA/RNA integrity during transport and storage
Ultraconserved Element (UCE) Probes	Target sequence capture	Provides standardized genomic markers across divergent taxa
Illumina Sequencing Reagents	Whole genome sequencing	Generates high-throughput sequence data for multiple individuals
IQ-TREE Software	Phylogenetic inference	Fast maximum likelihood analysis for gene tree estimation
D-Statistics (ABBA-BABA)	Gene flow detection	Tests for introgression between candidate species
GeoReferencing Tools	Spatial analysis	Documents sampling locations for geographic variation analysis

The case studies of American milksnakes and ocean sunfishes demonstrate both the power and limitations of the multispecies coalescent model in resolving controversial taxonomic questions. When applied critically with adequate sampling and integration of multiple data types, the MSC provides a robust framework for understanding species boundaries and evolutionary history. However, uncritical application without consideration of geographic variation, gene flow, and other biological data can lead to taxonomic over-splitting that does not reflect evolutionary reality.

These insights are particularly valuable for researchers working on species delimitation, conservation prioritization, and understanding cross-species gene flow. As MSC methods continue to evolve, incorporating more realistic models of gene flow and demographic history, their utility in taxonomy and phylogenomics will only increase. Nevertheless, taxonomic decisions with significant conservation implications should always incorporate multiple lines of evidence beyond genetic data alone, including morphology, ecology, behavior, and careful analysis of geographic distributions.

In phylogenomics, choosing an appropriate statistical model is paramount for obtaining an accurate representation of species' evolutionary histories. The multispecies coalescent (MSC) model and concatenation approach represent two fundamentally different frameworks for inferring species trees from genomic data. The concatenation method assumes that all loci share the same underlying genealogy (topologically congruent gene trees), effectively treating the entire multi-locus dataset as a single supergene [74]. In contrast, the MSC model explicitly accounts for the fact that individual gene trees can differ from the species tree due to biological processes such as incomplete lineage sorting (ILS), which is particularly common when species divergences occur in rapid succession [74] [24]. The MSC model naturally accommodates this gene tree heterogeneity by modeling the genealogical relationships of DNA sequences sampled from multiple species within a probabilistic framework that integrates both phylogenetic relationships and population genetic processes [15] [24].

Bayes factors provide a rigorous statistical framework for comparing these competing models by assessing the marginal likelihood of the data under each model, thereby offering a principled approach to model selection in phylogenomics. This model comparison approach is essential because the inconsistency of concatenation methods under certain conditions has been mathematically proven, but these proofs assume that the MSC model is true [74]. Empirical model validation and comparison through Bayes factors thus provide critical evidence for justifying model choice in practical applications. As we will demonstrate, Bayesian model validation and comparison consistently favor the MSC model over concatenation across diverse phylogenomic datasets, revealing that the concatenation assumption of congruent gene trees rarely holds for datasets with more than ten loci [74].

Theoretical Foundation

The Multispecies Coalescent Model

The multispecies coalescent model extends single-population coalescent theory to multiple species, providing a stochastic framework that describes the genealogical relationships of DNA sequences sampled from different species [24]. This model integrates two fundamental evolutionary processes: the phylogenetic process of species divergences and the population genetic process of coalescence, where lineages trace back to their common ancestors [15]. The MSC model accommodates the crucial observation that gene trees from different loci can have topologies that differ from the species tree due to incomplete lineage sorting, particularly when speciation events occur in rapid succession relative to effective population sizes [24].

The probability of gene tree-species tree discordance under the MSC can be quantified mathematically. For a rooted three-taxon tree, the probability that a gene tree is congruent with the species tree is given by:

P(congruence) = 1 - (2/3)exp(-T)

where T represents the branch length in coalescent units, defined as T = t/(2Ne), with t being the number of generations between speciation events and Ne being the effective population size [24]. This mathematical formulation highlights how shorter internal branches (smaller T) and larger effective population sizes increase the probability of gene tree-species tree discordance, creating the "anomaly zone" where the most likely gene tree topology differs from the species tree topology [15].

The MSC model operates under several key assumptions, including complete isolation after species divergence (no gene flow), neutrality of mutations, no recombination within loci, and free recombination between loci [24]. The model parameters typically include species divergence times (τ) and population sizes for both current and ancestral species (θ), which are estimated from the sequence data [24]. The joint probability density of gene trees under the MSC provides the foundation for full-likelihood methods of species tree estimation and allows for inference of species divergence times, population sizes, and, in extended models, cross-species gene flow [15].

The Concatenation Model

The concatenation model, in contrast, combines sequences from multiple loci into a single supermatrix, assuming that all genes share the same underlying tree topology and evolutionary history. This approach effectively treats gene tree heterogeneity as negligible and assumes that any discordance among gene trees results from estimation error rather than biological processes [74]. From a statistical perspective, the concatenation model represents a special case of the MSC model where all gene trees are constrained to be topologically identical [74] [3]. This constraint simplifies the model but introduces potential biases when the assumption of congruent gene trees is violated, which empirical evidence suggests is common in phylogenomic datasets [74].

Bayes Factors for Model Comparison

Bayes factors provide a robust Bayesian approach to model comparison by quantifying the evidence in favor of one model over another, based on their marginal likelihoods. The Bayes factor (BF) comparing Model 1 (MSC) to Model 2 (concatenation) is calculated as:

BF = P(D|MSC) / P(D|Concatenation)

where P(D|M) represents the marginal likelihood of the data D under model M. A Bayes factor greater than 1 indicates support for the MSC model, with values greater than 10 considered strong evidence [74]. The marginal likelihoods for these models are typically estimated using computational approaches such as path sampling or stepping-stone sampling, which are implemented in Bayesian phylogenetic software packages like BEAST 2 [57] [75].

Table 1: Bayes Factor Interpretation Guidelines

Bayes Factor Value	Interpretation	Support for MSC Model
1-3	Anecdotal	Barely worth mentioning
3-10	Substantial	Positive
10-30	Strong	Strong
30-100	Very strong	Very strong
>100	Decisive	Extreme

Quantitative Evidence from Empirical Studies

Comprehensive analyses across diverse phylogenomic datasets provide compelling quantitative evidence favoring the MSC model over concatenation. A landmark study examining 47 phylogenomic datasets collected across the tree of life revealed systematic violations of concatenation assumptions and superior performance of MSC models [74] [76].

Table 2: Model Rejection Rates Across 47 Phylogenomic Datasets

Model Type	Rejection Rate	Primary Causes of Rejection
Substitution Models	44% of loci	Poor fit to sequence evolution patterns, correlated with GC content and proportion of informative sites
Concatenation Model	38% of loci	Violation of topologically congruent gene trees assumption
MSC Model	11% of loci	Gene flow, hybridization, model violations

This study found that the proportion of GC content and informative sites were both negatively correlated with the fit of substitution models across loci, highlighting important factors influencing model adequacy [74]. Furthermore, substantial violations of the concatenation assumption of congruent gene trees were consistently observed across six major taxonomic groups: birds, mammals, fish, insects, reptiles, and other invertebrates [74]. Bayesian model validation strongly favored the MSC over concatenation across all datasets examined, with the concatenation assumption rarely holding for phylogenomic datasets with more than ten loci [74].

Importantly, loci rejecting the MSC model were found to have minimal impact on species tree estimation, suggesting that the MSC framework is robust to minor model violations [74]. This robustness, combined with its superior statistical performance, reinforces the value of the MSC model for phylogenomic inference.

Experimental Protocols

Bayesian Model Comparison Protocol

Objective: To compare the fit of MSC and concatenation models to multi-locus sequence data using Bayes factors.

Materials and Software Requirements:

• BEAST 2 with appropriate model packages [75]
• Sequence alignments for multiple loci
• Computational resources adequate for Bayesian phylogenetic inference

Procedure:

• Data Preparation:

  • Prepare sequence alignments for each locus in NEXUS format
  • Partition data appropriately based on locus boundaries
  • Assess data quality and missing data patterns

• Model Configuration for MSC:

  • Specify the MSC model as the prior for species tree estimation
  • Set up independent substitution models for each locus
  • Configure population size models (constant or variable across species tree)
  • Set prior distributions for divergence times and population parameters

• Model Configuration for Concatenation:

  • Combine all loci into a single alignment
  • Specify a concatenation model that assumes a shared topology across all loci
  • Implement partitioned substitution models to accommodate rate variation among loci
  • Set prior distributions for tree parameters

• Marginal Likelihood Estimation:

  • Perform path sampling or stepping-stone sampling for both models
  • Use at least 50 steps with 100,000 MCMC iterations per step
  • Ensure chain convergence through multiple independent runs
  • Calculate Bayes factors from marginal likelihood estimates

• Model Adequacy Assessment:

  • Perform posterior predictive simulations for both models
  • Compare observed data statistics to predicted distributions
  • Identify specific aspects where models may be inadequate

Troubleshooting Tips:

• If MCMC chains fail to converge, increase chain length and adjust proposal mechanisms
• For large datasets, consider approximate methods for marginal likelihood estimation
• Validate model configurations using simulated datasets with known parameters

MSC Model Fit Testing with P2C2M2

Objective: To assess the fit of the MSC model to a given dataset and identify loci that significantly deviate from model assumptions.

Materials and Software Requirements:

• R package P2C2M2 [57] [75]
• Gene trees and species trees estimated from sequence data
• BEAST 2 for Bayesian phylogenetic inference [75]

Procedure:

• Input Data Preparation:

  • Estimate gene trees for each locus using Bayesian methods in BEAST 2
  • Estimate the species tree under the MSC model
  • Format trees according to P2C2M2 requirements

• Model Fit Test Configuration:

  • Specify the test statistics for evaluating model fit
  • Set the number of posterior predictive simulations (typically 100-1000)
  • Configure computational parameters for the analysis

• Running the Analysis:

  • Execute P2C2M2 with the prepared input files
  • Monitor convergence and completion of simulations
  • Extract test results and summary statistics

• Results Interpretation:

  • Identify loci with significant departures from MSC expectations
  • Investigate potential biological causes for departures (e.g., gene flow, selection)
  • Consider alternative models for problematic loci (e.g., network models)

• Follow-up Analyses:

  • For loci showing significant departures, implement models that account for gene flow
  • Use multispecies network coalescent models when historical gene flow is suspected [57]
  • Compare species tree estimates with and without problematic loci

MSC Model Testing Workflow

The Scientist's Toolkit: Essential Research Reagents and Software

Table 3: Essential Tools for MSC Analysis and Model Testing

Tool/Reagent	Type	Primary Function	Application Notes
BEAST 2	Software Package	Bayesian evolutionary analysis	Implements MSC models, molecular clock dating, and marginal likelihood estimation [75]
P2C2M2	R Package	MSC model fit testing	Assesses adequacy of MSC model, identifies loci with significant departures [57] [75]
PhyloNet	Software Package	Phylogenetic network inference	Infers species networks accounting for both ILS and gene flow [75]
Path Sampling	Algorithm	Marginal likelihood estimation	Estimates model evidence for Bayes factor calculation [75]
Multi-locus Sequence Data	Data	Phylogenomic inference	10-100+ loci recommended for reliable species tree estimation [74]
Posterior Predictive Simulation	Method	Model adequacy assessment	Compares observed data to predictions from fitted model [74]

Application to Complex Evolutionary Scenarios

The MSC framework with Bayes factor model comparison proves particularly valuable when analyzing complex evolutionary histories where multiple biological processes interact. A case study on the Liolaemus wiegmannii lizard species group demonstrated how departures from MSC assumptions can lead to artefactual topologies and biased node height estimates [57]. When strong evidence of model departures was detected in several loci, potentially due to historical gene flow, researchers implemented a multispecies network coalescent model that incorporated reticulation events [57]. This approach revealed previously unreported historical gene flow between the L. wiegmannii and L. montanus groups, highlighting how model testing can uncover complex evolutionary processes [57].

In practice, researchers are encouraged to routinely test the fit of the MSC model when estimating species trees and to follow a phylogenetic network approach when significant departures are detected and historical gene flow is considered a plausible explanation [57]. This model-testing framework ensures that inferences about species relationships accurately reflect the complex interplay of evolutionary processes rather than being artifacts of model misspecification.

Model Selection Decision Pathway

Bayes factor comparison provides a statistically rigorous framework for demonstrating the superiority of the multispecies coalescent model over concatenation in phylogenomic analyses. Empirical evidence from diverse datasets shows that the MSC model consistently outperforms concatenation, with significantly lower rejection rates and better accommodation of biological realities such as incomplete lineage sorting [74]. The protocols outlined here for Bayesian model comparison and MSC model fit testing equip researchers with practical tools for implementing this approach in their phylogenomic studies. As the field progresses, model-aware approaches that test assumptions and select appropriate models will be essential for accurate inference of evolutionary histories, particularly in groups with complex histories of diversification and gene flow [57] [15]. The integration of model testing into standard phylogenomic practice represents a crucial step toward more accurate and reliable reconstruction of the tree of life.

Conclusion

The integration of gene flow into the multispecies coalescent model represents a significant advancement in phylogenomics, moving beyond simplistic tree-like assumptions to embrace the complex, reticulate nature of evolution. The key takeaways are that the MSC framework, particularly the MSC-I and MSC-M models, provides a statistically robust method for inferring species histories that is consistently favored over concatenation, which often fails when gene tree conflict is present. Successful application requires careful attention to model assumptions, the use of validation techniques like posterior predictive simulation, and the integration of complementary data, such as from contact zones, to avoid pitfalls like over-splitting. For biomedical and clinical research, these refined models offer a more accurate picture of population histories and the movement of adaptive alleles, which is crucial for understanding the evolution of pathogens, tracing the history of disease-related genes, and correctly delimiting species that may be reservoirs of infection. Future directions will involve further model refinement to include factors like selection and more complex demographic scenarios, increased computational efficiency for large genomic datasets, and broader application in comparative genomics and association studies.

References

• 1. The performance of coalescent-based species tree estimation methods under models of missing data | BMC Genomics | Full Text [https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-018-4619-8]
• 2. A multispecies coalescent model for quantitative traits [https://elifesciences.org/articles/36482]
• 3. Implementing and testing the multispecies coalescent model [https://www.sciencedirect.com/science/article/pii/S1055790315003309]
• 4. The Multilocus Multispecies Coalescent: A Flexible New Model of Gene Family Evolution - PubMed [https://pubmed.ncbi.nlm.nih.gov/33169795/]
• 5. Computing the probability of gene trees concordant with the species tree in the multispecies coalescent - PubMed [https://pubmed.ncbi.nlm.nih.gov/33333117/]
• 6. Species Delimitation under the Multispecies Coalescent ... [https://sites.utexas.edu/antisense/2025/01/08/species-delimitation-under-the-multispecies-coalescent-with-gene-flow/]
• 7. Simulating gene trees under the multispecies coalescent and ... [https://bmcecolevol.biomedcentral.com/articles/10.1186/1471-2148-13-44]
• 8. Adaptive Introgression as an Evolutionary Force: A Meta ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12181397/]
• 9. How species evolve collectively: implications of gene flow ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC2600545/]
• 10. Gene flow enhances genetic diversity and local adaptation ... [https://www.sciencedirect.com/science/article/pii/S277273512500054X]
• 11. Gene flow | Research Starters [https://www.ebsco.com/research-starters/history/gene-flow]
• 12. A Quantitative Comparison of the Similarity between Genes ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3426559/]
• 13. Multispecies coalescent and its applications to infer ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8692950/]
• 14. Inference of Gene Flow between Species from Genomic Data ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12203007/]
• 15. Multispecies coalescent and its applications to infer ... [https://pubmed.ncbi.nlm.nih.gov/34987842/]
• 16. Review Decoding genomic landscapes of introgression [https://www.sciencedirect.com/science/article/pii/S0168952525001660]
• 17. How to Tackle Phylogenetic Discordance in Recent and ... [https://www.frontiersin.org/journals/plant-science/articles/10.3389/fpls.2021.765719/full]
• 18. Disentangling Sources of Gene Tree Discordance in ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC7875436/]
• 19. An approximate likelihood method reveals ancient gene ... [https://peercommunityjournal.org/articles/10.24072/pcjournal.359/]
• 20. Phylogenomic Discordance is Driven by Wide-Spread ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11906154/]
• 21. Teasing apart the sources of phylogenetic tree discordance ... [https://bmcplantbiol.biomedcentral.com/articles/10.1186/s12870-025-06963-3]
• 22. Resolving phylogenetic conflicts in Pandanales: the dual ... [https://www.frontiersin.org/journals/plant-science/articles/10.3389/fpls.2025.1511582/full]
• 23. Genome-based phylogenetics and species delimitation for ... [https://www.sciencedirect.com/science/article/abs/pii/S1055790325001757]
• 24. Multispecies coalescent process [https://en.wikipedia.org/wiki/Multispecies_coalescent_process]
• 25. Major patterns in the introgression history of Heliconius ... [https://elifesciences.org/articles/90656]
• 26. The Impact of Sequencing and Genotyping Errors on ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12359030/]
• 27. Efficient Bayesian inference under the multispecies coalescent with migration - PubMed [https://pubmed.ncbi.nlm.nih.gov/37871206/]
• 28. 7 Markov Chain Monte Carlo | Doing Bayesian Data ... [https://bookdown.org/ajkurz/DBDA_recoded/markov-chain-monte-carlo.html]
• 29. mstree: A Multispecies Coalescent Approach for Estimating ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC7259675/]
• 30. The BPP program for species tree estimation ... - UCL Discovery [https://discovery.ucl.ac.uk/1473654/]
• 31. A Tutorial on the Use of BPP for Species Tree Estimation ... [https://hal.science/hal-02536475]
• 32. The Spectre of Too Many Species - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC6292489/]
• 33. To split or not to split? Multilocus phylogeny and molecular ... [https://bmcecolevol.biomedcentral.com/articles/10.1186/s12862-019-1422-3]
• 34. The effects of sampling on delimiting species from multi- ... [https://www.sciencedirect.com/science/article/abs/pii/S1055790312002588]
• 35. A Simulator for Network Multispecies Coalescent Models, ... [https://pubmed.ncbi.nlm.nih.gov/37254872/]
• 36. A simulator for network multispecies coalescent models ... [https://datadryad.org/dataset/doi:10.5061/dryad.02v6wwq6x]
• 37. Divergence time estimation using ddRAD data and an ... [https://www.nature.com/articles/s41598-022-07877-y]
• 38. Estimating demographic parameters from large-scale ... [https://bmcgenomdata.biomedcentral.com/articles/10.1186/1471-2156-13-22]
• 39. Recent advances in Bayesian inference of isolation-with- ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6944047/]
• 40. Inference of gene flow in the process of speciation [https://www.sciencedirect.com/science/article/abs/pii/S0040580921000186]
• 41. A New Isolation with Migration Model along Complete ... [https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1003125]
• 42. Inference of Cross-Species Gene Flow Using Genomic ... [https://pubmed.ncbi.nlm.nih.gov/40421982/]
• 43. Challenges in Species Tree Estimation Under the ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC5161269/]
• 44. The impact of taxonomic change on conservation: Does it ... [https://www.sciencedirect.com/science/article/abs/pii/S0006320709003218]
• 45. Disruptions Caused by Changing Taxonomy and Ways of ... [https://forum.inaturalist.org/t/disruptions-caused-by-changing-taxonomy-and-ways-of-reconciling/58796]
• 46. Recombination-aware phylogenomics [https://www.sciencedirect.com/science/article/abs/pii/S016953472500179X]
• 47. Evaluating ARG-estimation methods in the context of ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12005257/]
• 48. Assessing simulation-based supervised machine learning ... [https://www.nature.com/articles/s41437-025-00773-x]
• 49. Hierarchical Heuristic Species Delimitation Under the ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11637770/]
• 50. Species Delimitation Using Genomic Data [https://pmc.ncbi.nlm.nih.gov/articles/PMC11974488/]
• 51. Data: Hierarchical heuristic species delimitation under the ... [https://datadryad.org/dataset/doi:10.5061/dryad.jm63xsjhc]
• 52. Identifying model violations under the multispecies ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6956792/]
• 53. The Importance of Contact Zones for Distinguishing ... [https://pubmed.ncbi.nlm.nih.gov/35993885/]
• 54. Modern phylogenomics: Building phylogenetic trees using ... [https://pure.psu.edu/en/publications/modern-phylogenomics-building-phylogenetic-trees-using-the-multis/]
• 55. Reproductive isolation in a three-way contact zone [https://pubmed.ncbi.nlm.nih.gov/38235507/]
• 56. Species delimitation using genomic data to resolve ... [https://www.sciencedirect.com/science/article/pii/S1055790322002846]
• 57. Revisiting the Multispecies Coalescent Model fit with an ... [https://pubmed.ncbi.nlm.nih.gov/40637759/]
• 58. A structured coalescent model reveals deep ancestral ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11985351/]
• 59. BEAST 2.5: An advanced software platform for Bayesian ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6472827/]
• 60. Unicore Enables Scalable and Accurate Phylogenetic ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12203212/]
• 61. Multiple Sequence Alignment Averaging Improves Phylogeny Reconstruction - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC6657586/]
• 62. LMAP_S: Lightweight Multigene Alignment and Phylogeny eStimation | BMC Bioinformatics | Full Text [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-019-3292-5]
• 63. Evaluating the utility of deep genome skimming for ... [https://www.sciencedirect.com/science/article/pii/S2468265925000836]
• 64. PPC-overview Graphical posterior predictive checking [https://www.rdocumentation.org/packages/bayesplot/versions/1.14.0/topics/PPC-overview]
• 65. Identifying model violations under the multispecies coalescent model ... [https://pubmed.ncbi.nlm.nih.gov/31949994/]
• 66. Examining Gene : Key Flow and Their Insights - Simple Models Science [https://scisimple.com/en/articles/2025-10-05-examining-gene-flow-key-models-and-their-insights--a37llxd]
• 67. Software integration testing for Bayesian models - Adam Howes [https://athowes.github.io/posts/2024-09-15-bayesian-tests/]
• 68. 9. End to End Bayesian Workflows [https://bayesiancomputationbook.com/markdown/chp_09.html]
• 69. 1. Introduction — Bayesian workflow documentation [https://francescacapel.com/BayesianWorkflow/notebooks/introduction.html]
• 70. The Multispecies Coalescent Over-Splits Species in ... [https://pubmed.ncbi.nlm.nih.gov/31180508/]
• 71. Taxonomic Split 18099 (Committed on 2017-03-27) [https://www.inaturalist.org/taxon_changes/18099]
• 72. Recommendations from the IUCN molidae review panel [https://www.sciencedirect.com/science/article/abs/pii/S0308597X23002932]
• 73. Taxonomic Troubles: When is a sunfish not a sunfish? [https://sunfishresearch.wordpress.com/2017/06/19/taxonomic-troubles-when-is-a-sunfish-not-a-sunfish/]
• 74. The Multispecies Coalescent Model Outperforms ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC7302055/]
• 75. Data: Revisiting the multispecies coalescent model fit with ... [https://datadryad.org/dataset/doi:10.5061/dryad.jm63xsjk5]
• 76. The Multispecies Coalescent Model Outperforms ... - PubMed [https://pubmed.ncbi.nlm.nih.gov/32011711/]