Hierarchical heuristic species delimitation under the multispecies coalescent model with migration

The multispecies coalescent (MSC) model accommodates genealogical fluctuations across the genome and provides a natural framework for comparative analysis of genomic sequence data to infer the history of species divergence and gene flow. Given a set of populations, hypotheses of species delimitation (and species phylogeny) may be formulated as instances of MSC models (e.g., MSC for one species versus MSC for two species) and compared using Bayesian model selection. This approach, implemented in the program bpp, has been found to be prone to over-splitting. Alternatively, heuristic criteria based on population parameters under the MSC model (such as population/species divergence times, population sizes, and migration rates) estimated from genomic sequence data may be used to delimit species. Here we extend the approach of species delimitation using the genealogical divergence index (𝑔𝑑𝑖) to develop hierarchical merge and split algorithms for heuristic species delimitation and implement them in a python pipeline called hhsd. Applied to data simulated under a model of isolation by distance, the approach was able to recover the correct species delimitation, whereas model comparison by bpp failed. Analyses of empirical datasets suggest that the procedure may be less prone to over-splitting. We discuss possible strategies for accommodating paraphyletic species in the procedure, as well as the challenges of species delimitation based on heuristic criteria.