Bio Ecological Genomics Species Delimitation

Version Compatibility

Reference examples tested with: ASAP (current CLI), mPTP (current C++), bPTP (Python from github.com/zhangjiajie/PTP), splits 1.0+ for GMYC, BPP 4.7+, SNAPP/BEAST 2.7+, DELINEATE (current Python), Dsuite 0.5+, ape 5.7+, fossil 0.4+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Species Delimitation

"Delineate species boundaries from my DNA barcoding or genomic data" -> Apply distance-based (ASAP) and tree-based (mPTP) methods for primary delimitation, multilocus coalescent (BPP) for genomic confirmation, DELINEATE for speciation-process modeling to address MSC oversplitting, integrative-taxonomy validation across multiple lines of evidence, and Dsuite to test introgression before claiming sister relationships.

• CLI: asap (downloadable C binary) OR web (https://bioinfo.mnhn.fr/abi/public/asap/) for single-locus primary delimitation
• CLI: mptp (https://github.com/Pas-Kapli/mptp) for tree-based multi-rate PTP (faster than bPTP)
• Python: bPTP (https://github.com/zhangjiajie/PTP) for legacy tree-based delimitation — NOTE: Python not R
• R: splits::gmyc() for GMYC on ultrametric trees
• CLI: bpp (https://bpp.github.io) for multilocus MSC delimitation; control-file driven
• CLI: Dsuite for D-statistics / f4-ratio / f-branch introgression testing

The Single Most Important Modern Insight -- MSC Methods Delimit Structure, NOT Species

Sukumaran & Knowles 2017 PNAS 114(7):1607-1612 established that BPP, BFD*, and other multispecies-coalescent methods delimit genetic structure, not species. The mathematical reason: the MSC model treats every panmictic population as a "species" in the parameterization. Applied to data with substructure (isolated demes within a species), MSC methods partition that structure into "species" — leading to systematic oversplitting in published literature, especially in lizards, frogs, geckos, and insects.

Modern best practice (post-2021):

1. Run MSC methods (BPP, BFD*) as PRELIMINARY enumeration of candidate population lineages
2. Apply DELINEATE (Sukumaran 2021 PLoS Comput Biol 17:e1008924) to test which lineages have completed speciation vs are in the speciation process
3. Validate via integrative taxonomy (Padial 2010; Carstens 2013) — morphological, ecological, geographic congruence

A second cornerstone: the barcoding gap is OFTEN absent in real data (Meyer & Paulay 2005 PLoS Biol 3:e422). ABGD and ASAP are theoretically grounded in the gap; when the gap is absent, both fail gracefully but their output is unreliable. Always inspect the pairwise distance histogram.

A third: Dsuite-checked introgression must precede sister-species claims. D-statistics, f4-ratio, and f-branch (Malinsky 2021) distinguish admixture from incomplete lineage sorting; without this check, ILS-driven discordance is misinterpreted as gene flow.

Algorithmic Taxonomy

Decision Tree by Scenario

ASAP — The Modern Successor to ABGD

Goal: Primary species hypothesis from a single-locus alignment via hierarchical clustering of pairwise distances, scored by a new asap-score (Puillandre 2021).

Approach: Run ASAP via the web interface (https://bioinfo.mnhn.fr/abi/public/asap/) OR the downloadable C binary. ASAP ranks candidate partitions by asap-score (lower = better); always inspect the top 5-10 partitions, not just the best, and check for large score gaps that signal robust partitioning. Output: ranked species partitions with asap-score, p-value, and threshold distance per partition.

# ASAP CLI (downloadable C binary)
# Substitution model: K2P standard for COI; p-distance for very closely related; JC for general
./asap -d K2P -o asap_results/ aligned_sequences.fasta

# Inspect output
# - asap_output.csv: ranked partitions by asap-score (lower = better)
# - Always look at top 5-10 partitions; large score gaps = robust partitioning
# - The "best" partition is the top-ranked; report secondary partitions too

# Python ASAP-style analysis for inspection (NOT a substitute for ASAP)
# Useful for understanding the distance landscape
from Bio import AlignIO
from Bio.Phylo.TreeConstruction import DistanceCalculator
from scipy.cluster.hierarchy import fcluster, linkage
from scipy.spatial.distance import squareform
import numpy as np
import matplotlib.pyplot as plt

alignment = AlignIO.read('aligned_sequences.fasta', 'fasta')
calc = DistanceCalculator('identity')  # use K2P-friendly model in practice
dm = calc.get_distance(alignment)

names = [r.id for r in alignment]
n = len(names)
dist_array = np.zeros((n, n))
for i in range(n):
    for j in range(n):
        dist_array[i][j] = dm[names[i], names[j]]

# Inspect distance histogram for barcoding gap (often ABSENT per Meyer & Paulay 2005)
upper = dist_array[np.triu_indices(n, k=1)]
plt.hist(upper, bins=50)
plt.xlabel('Pairwise distance')
plt.ylabel('Frequency')
plt.title('Distance histogram — look for bimodality (gap)')
plt.savefig('barcoding_gap_histogram.pdf')

# Hierarchical clustering at threshold scan
condensed = squareform(dist_array)
Z = linkage(condensed, method='average')
thresholds = np.arange(0.01, 0.10, 0.005)
for t in thresholds:
    clusters = fcluster(Z, t=t, criterion='distance')
    print(f'Threshold {t:.3f}: {len(set(clusters))} groups')

mPTP — Multi-Rate PTP for Heterogeneous Intraspecific Rates

Goal: Delimit species from a phylogeny by detecting the transition from speciation to coalescent branching rates, allowing per-species intraspecific rates (Kapli 2017).

Approach: Run mPTP on a rooted ML tree (e.g., from RAxML or IQ-TREE). mPTP's multi-rate model is more flexible than bPTP's single-rate; for shallow phylogenies with similar rates, bPTP may be competitive. mPTP is C++ (5+ orders of magnitude faster than bPTP), available from https://github.com/Pas-Kapli/mptp.

# mPTP installation: build from C source
# git clone https://github.com/Pas-Kapli/mptp
# cd mptp && ./autogen.sh && ./configure && make

# Run mPTP ML inference
mptp --ml --tree_file rooted_tree.nwk --output_file mptp_ml

# MCMC for posterior support
mptp --mcmc 1000000 --mcmc_sample 1000 --mcmc_burnin 100000 \
     --tree_file rooted_tree.nwk --output_file mptp_mcmc

# Output: species partition with branch annotations
# > 0.95 posterior support: strong; 0.80-0.95: moderate; < 0.80: uncertain

bPTP — Legacy Tree-Based Delimitation (Python, NOT R)

Goal: Bayesian Poisson Tree Processes delimitation with single intraspecific rate (Zhang 2013).

Approach: Run bPTP via the Python package from https://github.com/zhangjiajie/PTP (NOTE: this is a Python tool, NOT an R package; install.packages('PTP') does not exist). For heterogeneous-rate data, prefer mPTP.

# bPTP Python install (NOT R)
pip install git+https://github.com/iTaxoTools/PTP-pyqt5

# Run bPTP MCMC
python -m PTP.PTP -t rooted_tree.nwk -o bptp_results \
    --type bayesian --ngen 100000 --burnin 0.1 --seed 42

# Output: posterior support per species partition (> 0.95 strong)
# bPTP over-splits when populations have strong geographic substructure
# Cross-check with ASAP and/or mPTP for consensus

GMYC with Mandatory Ultrametric Tree

Goal: Detect the threshold on an ultrametric phylogeny where branching transitions from interspecific (Yule) to intraspecific (coalescent).

Approach: Convert ML tree to ultrametric with ape::chronos (penalized likelihood) or BEAST (Bayesian time-calibration), then run GMYC. Single-threshold vs multi-threshold; multi-threshold handles heterogeneous rates across the tree.

library(splits)
library(ape)

tree <- read.tree('rooted_tree.nwk')

# Convert to ultrametric (REQUIRED for GMYC)
# chronos: penalized likelihood; lambda=1 mid-clock
# For rigorous time-calibration, use BEAST2 instead
ultrametric_tree <- chronos(tree, lambda = 1)
class(ultrametric_tree) <- 'phylo'
stopifnot(is.ultrametric(ultrametric_tree))

# Single-threshold GMYC
gmyc_single <- gmyc(ultrametric_tree, method = 'single')
summary(gmyc_single)
cat('GMYC species (single-threshold):', gmyc_single$entity[1], '\n')
cat('LR test p-value:', gmyc_single$p.value[1], '\n')

# Multi-threshold for heterogeneous rates
gmyc_multi <- gmyc(ultrametric_tree, method = 'multiple')
summary(gmyc_multi)

# Compare single vs multi-threshold; report both
species_single <- spec.list(gmyc_single)
species_multi <- spec.list(gmyc_multi)

BPP A10/A11 — Multilocus MSC Delimitation with Calibrated Priors

Goal: Bayesian multilocus species delimitation under the multispecies coalescent (Yang & Rannala 2010; current Flouri 2018 BPP v4).

Approach: Configure BPP A10 (joint delimitation + species tree estimation) with priors on theta and tau ESTIMATED FROM DATA — NOT defaults. Yang 2015 Curr Zool 61:854-865 gives explicit guidance: thetaprior = G(2, 2000) for small/moderate data; tighter for large genomic data. Wrong priors dominate the analysis.

# BPP control file for A10 analysis
# CRITICAL: priors below are CALIBRATED to data, not defaults
# Defaults are tuned for primate-mammal data; insects/fish/plants need adjustment

seed = 12345

seqfile = alignment.phy
Imapfile = imap.txt                   # map individuals -> putative species
outfile = bpp_results.txt
mcmcfile = bpp_mcmc.txt

speciesdelimitation = 1 0 2 0.5       # 1 = delimitation mode; rjMCMC parameters
speciestree = 1                       # 1 = estimate species tree
speciesmodelprior = 1                 # 1 = uniform prior on rooted trees

species&tree = 4  speciesA speciesB speciesC speciesD
                  10 8 12 6
                  ((speciesA, speciesB), (speciesC, speciesD));

# CALIBRATED theta prior (per-locus nucleotide diversity)
# G(2, 2000) -> mean = 0.001 (calibrate from observed per-locus pi)
# Yang 2015 Curr Zool 61:854 has taxon-specific recommendations
thetaprior = 2 2000

# CALIBRATED tau prior (divergence times)
# Calibrate from observed maximum genetic distance / mutation rate
tauprior = 2 2000

# MCMC
nsample = 100000
sampfreq = 2
burnin = 50000

# Run multiple independent chains; check convergence

# Run BPP with multiple independent chains
for i in 1 2 3 4; do
    bpp --cfile bpp_control.bpp --seed $((1000 * i)) --out bpp_run_$i.txt &
done
wait

# Compare chains; posterior probabilities must agree across chains
# PP > 0.95 for a delimitation: strong support
# PP 0.50-0.95: moderate; integrate other evidence

DELINEATE — Addressing the Sukumaran-Knowles Oversplitting

Goal: Test whether candidate lineages from BPP represent fully-formed species OR incomplete-speciation structure (Sukumaran 2021).

Approach: Provide DELINEATE with a guide tree of population lineages and a hypothesis about which are species. The output: probability that each population lineage is a true species vs an incipient lineage. THE post-2021 best practice complement to BPP.

# DELINEATE (https://github.com/jsukumaran/delineate)
# Python tool
pip install delineate

# Configure delineate input
# - guide tree: rooted Newick with population lineages as tips
# - constraint file: which tips are confirmed species, which to test
delineate-estimate partitions \
    --tree-file population_lineage_tree.nwk \
    --config-file delineate_config.json \
    --output-prefix delineate_run

# Output: probability that each lineage is a species
# Use ALONGSIDE BPP; do not claim species from BPP alone post-2021

Dsuite — Introgression vs Incomplete Lineage Sorting

Goal: Test for introgression before claiming sister-species relationships, distinguishing admixture from ILS.

Approach: Compute D-statistic (ABBA-BABA), f4-ratio (admixture fraction), and f-branch (Malinsky 2021 Mol Ecol Resour 21:584-595) on a VCF with population assignments. Non-zero D indicates introgression between H3 and one of H1/H2.

# Dsuite (https://github.com/millanek/Dsuite)
# Sets file: 3 columns (sample_id, population/species, optional outgroup tag)
# Tree file: Newick with population names matching the sets file

# All-trio D-statistic
Dsuite Dtrios -t species_tree.nwk sets.txt input.vcf.gz

# F-branch decomposition (more interpretable than raw D)
Dsuite Fbranch species_tree.nwk Dtrios_output_tree.txt > fbranch.txt

# Plot f-branch matrix
dtools.py fbranch.txt species_tree.nwk

# D > 0 with |Z| > 3: significant ABBA-BABA imbalance -> introgression
# f4-ratio gives the admixture fraction estimate

Per-Method Failure Modes

MSC delimitation oversplits population structure as species

Trigger: Running BPP, BFD*, or other multispecies-coalescent methods on data with within-species substructure (isolated demes), then claiming the delimited lineages are species.

Mechanism: The MSC model treats every panmictic population as a "species" in its parameterization (Sukumaran & Knowles 2017). When applied to data with population substructure, MSC partitions that structure into "species."

Symptom: Many "species" delimited; each with low sample size; lineages defined by geography rather than morphology; BPP posterior support is high but biological reality is unclear.

Fix: Run DELINEATE (Sukumaran 2021) on the BPP output to test which lineages have completed speciation. Validate via integrative taxonomy (Padial 2010): require congruence across genetic + morphological + ecological evidence before publishing species.