Data and code from: Diagnosing systematic errors and incomplete lineage sorting in the deep phylogeny of Collembola
Data files
Jul 15, 2026 version files 2.10 GB
-
align_MAFFT.sh
6.18 KB
-
analyses.tar.gz
2.10 GB
-
BUSCO_extraction.sh
7.16 KB
-
calculate_branch_length.py
25.68 KB
-
convert-exchangeabilities.py
3.99 KB
-
convert-site-dists-to-k_eff.py
1.58 KB
-
convert-site-dists.py
2.37 KB
-
extract_position_fasta.py
5.82 KB
-
gene_trees.sh
5.80 KB
-
gene-wise_likelihood.sh
1.94 KB
-
loci_filtering_alignment-based.sh
19.42 KB
-
loci_filtering_tree-based.sh
22.17 KB
-
matrix_generation.sh
7.88 KB
-
README.md
7.74 KB
-
site-wise_likelihood.sh
2.14 KB
-
trimming_alignments.sh
12.75 KB
Abstract
Resolving deep phylogeny is complicated by systematic errors and biological sources of gene tree discordance, a challenge that is clearly presented by Collembola, one of the earliest diverging hexapod lineages. To address the long-standing conflict among the four collembolan orders, we expanded taxon sampling (113 species) and marker representation (4,070 single-copy orthologs), conducting a systematic evaluation of potential error sources. Using both concatenation and multispecies coalescent approaches under site-homogeneous (LG) and site-heterogeneous (CAT-PMSF) models, we recovered two primary competing topologies: a Neelipleona-first and an Entomobryomorpha-first hypothesis. Pronounced branch-length heterogeneity across orders and two exceptionally short internal branches suggested susceptibility to systematic error. Analyses of empirical data, combined with extensive simulations, showed that the Neelipleona-first topology arises under conditions known to induce long-branch attraction, including branch-length imbalance, compositional heterogeneity, and model misspecification, with fast-evolving and compositionally constrained sites further amplifying these artifacts. Coalescent simulations demonstrated that incomplete lineage sorting and gene tree estimation error jointly account for much of the deep gene tree-species tree discordance. In contrast, analyses using site-heterogeneous models and multispecies coalescent approaches, both intended to reduce systematic errors, consistently supported the Entomobryomorpha-first topology, recovering Entomobryomorpha + (Symphypleona + (Poduromorpha + Neelipleona)). Our findings clarify the mechanistic origins of phylogenomic conflict in Collembola and highlight the need to jointly consider systematic error and ILS when resolving ancient radiations. We propose the name 'Brachyantennamorpha' for the clade Poduromorpha + Neelipleona.
Dataset DOI: 10.5061/dryad.dv41ns2cj
This dataset contains all multiple sequence alignments (MSAs), gene trees, species‑tree inferences, locus‑property calculations, systematic‑error diagnostics, and coalescent simulations used in the study Diagnosing systematic errors and incomplete lineage sorting in the deep phylogeny of Collembola.
All analyses were performed using empirical transcriptomic and genomic data from Collembola and outgroup taxa. The dataset enables full reproducibility of alignment processing, phylogenetic inference, support evaluation, and modeling of genealogical heterogeneity.
Description of the data and file structure
All data files are packaged within analyses.tar.gz. Below is an overview of all folders and the relationships among files.
1-individual_MSAs_GeneTrees/
Contains per‑locus multiple sequence alignments (MSAs) in FASTA format and corresponding gene trees in Newick format.
Datasets included:
complete_113spp/— 113 taxa; used for matrix0 and matrix1reduced_64spp/— 64 taxa; used for matrix3
Additional files:
-
metricsXXX.csvSequence‑ and gene tree–based features generated using
extract_msa_tree_features.py(from the MSA-and-tree-metrics-exploration tool). Metrics include alignment length, composition metrics, DVMC, RCFV, treeness, evolutionary rate, RF distance, etc.
2-gene_properties_detection_and_sensitivity_tests/
Contains locus‑property metrics and wASTRAL species trees inferred under four subsampling regimes for each property (DVMC, treeness, occupancy, ABS, RCFV, etc.).
Key file:
wastral.4070.tre— wASTRAL tree inferred from all 4,070 gene trees
Used to evaluate sensitivity of species‑tree inference to locus properties.
3-phylogeny/
Each matrixX/ directory contains:
- Concatenated supermatrix (FASTA or PHYLIP)
- Partition file
- Species‑tree estimates under various models:
- LG+F+R4
- Partitioned LG
- PMSF (GTR+CAT–based)
- wASTRAL
- Intermediate IQ‑TREE and PhyloBayes output:
- log files
.ckp.gzcheckpoint files.sitelh(site‑wise likelihoods).rate(site‑rate estimates)
CF/ subfolder contains:
- gCF gene concordance factors
- sCF site concordance factors
- Supporting plots and tables
4-Phylogenetic_support_distribution/
Contains log‑likelihood–based support summaries:
- SLS: site log-likelihoods
- GLS: gene log-likelihoods
- ΔSLS, ΔGLS between competing topologies
Reference topologies:
- T1 = LG.tre
- T2 = CAT-PMSF.tre
Used to classify genes/sites into pro‑T1, pro‑T2, or ambiguous categories.
5-systematic_errors/
Includes analyses of rate heterogeneity, branch‑length heterogeneity, compositional heterogeneity, substitution‑pattern heterogeneity, and model adequacy.
5.1 Across-site rate heterogeneity
Subfolders include:
modelR4_vs_noR4/— LG+F+R4 vs LG+F comparisonsslow_sites/— removal of fastest‑evolving 25%, 50%, 75% of sites
Outputs:
- filtered MSAs
- inferred trees
- logs
5.2 Branch-length heterogeneity
Subfolders include:
LG+G4/CAT-GTR+G4/
Contain posterior trees, branch lengths, and mapping files, used for branch‑length heterogeneity tests and Figure 2.
5.3 Site-compositional heterogeneity
Subfolders include:
CAT-PMSF_T1/LG_T2/
Contain:
- site‑wise log-likelihoods
- ΔlnL(T2–T1)
- Keff values
- compositional‑constraint metrics
5.4 Substitution‑pattern heterogeneity
Two major analyses:
-
relative_model_fit/: Model comparisons among LG, LG4X variants, C20 variants, ±F, ±G. -
absolute_model_fit/: Posterior predictive simulations (PhyloBayes).Files include:
phylobayes.tar.gz(chains, posterior, tracecomp)simulation_fasta.tar.gz(50 replicates)simulation_trees.tar.gz(trees inferred per replicate)
6-genealogical_heterogeneity/
Includes incomplete lineage sorting (ILS) metrics and coalescent simulations.
6.1 ILS_IH_indices/
Contains:
- ILS-index
- IH-index
- quartet imbalance (Phytop)
6.2 simulation/
Includes:
-
1-brlen_mutation_units/: Species tree with branch lengths in mutation units -
2-theta/: Empirical theta values -
3-ultrametric_tree/: Ultrametric tree (R/ape chronos) -
4-phybase/: 10,000 simulated gene trees; “thetaModifyXXX” indicates modified θ values -
5-MSA_simulation/For each replicate:
- true gene trees (AliSim input)
- simulated MSAs
- estimated gene trees
- logs
Scripts included
A brief overview of all executable scripts:
- align_MAFFT.sh — MAFFT/MAGUS locus alignment
- BUSCO_extraction.sh — extraction of BUSCO orthologs
- calculate_branch_length.py — branch‑length & patristic‑distance computation
- convert-exchangeabilities.py — convert PhyloBayes matrices to PAML/IQ‑TREE format
- convert-site-dists.py — convert CAT site profiles to PMSF format
- convert-site-dists-to-k_eff.py — compute Keff per site
- extract_position_fasta.py — extract selected alignment positions
- gene-wise_likelihood.sh — per‑locus AU tests
- gene_trees.sh — gene tree inference
- matrix_generation.sh — construct supermatrices
- loci_filtering_alignment-based.sh — alignment‑based locus filtering
- loci_filtering_tree-based.sh — tree‑based locus filtering
- site-wise_likelihood.sh — site-wise log-likelihood extraction
- trimming_alignments.sh — alignment trimming via trimAl/BMGE/ClipKIT
All scripts include usage instructions inside the file headers.
Sharing/Access information
Additional public resources
The MSA and gene‑tree metrics used for gene‑property analyses can be reproduced using:
- MSA‑and‑tree‑metrics‑exploration tool
https://github.com/xtmtd/MSA-and-tree-metrics-exploration
Data sources
Data were derived from:
- BUSCO ortholog datasets
- Empirical transcriptome/genome assemblies generated in this study
- Published Collembola transcriptomes (see main article for accession numbers)
- Public taxonomic databases used for backbone tree inference
Code/Software
All scripts and pipelines used in this study are open‑source. Unless otherwise noted, analyses were run on Linux (Ubuntu 20.04/22.04).
Major software used
Sequence & alignment processing
- MAFFT v7.520
- MAGUS v0.1
- TrimAl v1.4.1
- BMGE v1.12
- ClipKIT v1.3.0
- PhyKIT v1.19.4
- SeqKit v2.5.1
- csvtk v0.28.0
Orthology
- BUSCO v3 / v5
- TransDecoder v5.5.0
Gene-tree inference
- IQ‑TREE2 v2.2.5
- TreeShrink v1.3.8
- ASTRAL v5.7.8
Systematic‑error diagnostics
- TAPER v1.0.0
- Python 3.8+ (numpy, pandas, ETE3)
- Julia v1.9
Supermatrix construction
- PhyKIT v1.19.4
- FASconCAT‑G v1.05
Model adequacy & Bayesian inference
- PhyloBayes‑MPI v1.8c
- AliSim (IQ‑TREE package)
- FastTree v2.1 / v2.2
ILS & coalescent simulations
- ape v5.7.1
- phybase v2.0
- Phytop v0.3
General workflow
BUSCO extraction → alignment → trimming → alignment‑based filtering → gene tree inference → tree‑based filtering → matrix construction → species‑tree inference → support analyses → systematic‑error diagnostics → coalescent simulations → branch‑length & genealogical heterogeneity analyses
