Data and code from: Estimating divergence times and diversification rates with the unresolved fossilized birth-death process

Lim, Wonseop 1 ; Raskin, Levi 1 ; Li, Jacky1; Huelsenbeck, John 1 ; Nielsen, Rasmus 1

Research facility: University of California, Berkeley

Published Sep 17, 2026 on Dryad. https://doi.org/10.5061/dryad.w6m905r4v

Data files

Sep 17, 2026 version files 18.45 GB

Examples.zip

18.45 GB
README.md

14.55 KB

Abstract

Tip-dating using Fossilized Birth-Death (FBD) models offers a favorable alternative to conventional node-dating analyses of the evolutionary timescale, bypassing various conceptual and technical difficulties. However, despite their increasing popularity in phylogenetic studies, their use has remained rather limited in taxonomic scope because of the difficulty in collecting morphological data, which is often required for tip-dating analyses, and computational intractability when a large number of fossils are present. Recent studies, however, suggest that FBD models can also be used in a principled way as priors for molecular clock analyses, even in the absence of morphological data, relying only on temporal information and given clade assignments. Heath et al. (2014) introduced an "unresolved" FBD model in which uncertainty in fossil placements is analytically integrated out when morphological data are unavailable, providing a potential solution to this problem. Here, we show that the topological enumeration scheme used by Heath et al. (2014) to marginalize over unresolved fossil attachments overcounts configurations in which multiple fossil taxa can be resolved as a clade, and we provide a correction to this factor, yielding a valid marginalized unresolved likelihood. We apply our implementation to empirical datasets of Cetacea, Emydidae, and Palaeodictyopterida, and show that our method converges and yields reasonable estimates within a reasonable amount of time, even when hundreds of fossils are present.

Dataset DOI: 10.5061/dryad.w6m905r4v

Description of the data and file structure

Datasets and outputs from the empirical analyses, as well as the source code for the program UFBD (also available at: https://github.com/wnsplim/UFBD), accompanying the paper: "Estimating Divergence Times and Diversification Rates with the Unresolved Fossilized Birth-Death Process" (Lim et al., 2026).

Files and variables

File: UFBD-main.zip

Hosted on Zenodo

Description: Source code for the program UFBD used to analyze the empirical data. Its README documents the input formats, every configuration option, and the output columns. Build with CMake (see README), or use the static Linux binary in releases/.

File: Examples.zip

Description: Input data, configuration files, and MCMC output of the three empirical analyses of Lim et al. (2026). Running ufbd -config <config file> inside a folder reproduces its analysis.

Cetacea/
    config_cetacea.cfg
    cetacea.treefile
    cetacea.clades
    cetacea.fossils
    in.BV
    output/
Emydidae/
    config_emydidae.cfg
    emydidae.treefile
    emydidae.clades
    emydidae.fossils
    emydidae.fasta
    emydidae_partition.nex
    output/
Palaeodictyopterida/
    config_palaeodictyopterida.cfg
    palaeodictyopterida.fossils
    output/

Cetacea (whales, dolphins and porpoises): 72 extant species in the backbone tree; the sequence information enters as a Hessian file (approximate-likelihood analysis); 535 fossil species and 22 unsampled extant species; nine time bins for the rates.

Emydidae (pond turtles): 50 extant species, analyzed directly from a 15-locus alignment (full-likelihood analysis); 21 fossil species and 7 unsampled extant species; time-homogeneous rates.

Palaeodictyopterida (an extinct Paleozoic insect group): no backbone tree; 446 fossil species; seven time bins.

The sequence file and backbone tree for Cetacea were obtained from McGowen et al. (2020), after pruning all non-cetacean outgroups. The sequence file and backbone tree for Emydidae were obtained from the global turtle phylogeny of Thomson et al. (2021), retaining only emydid species. Fossil records for all three clades were obtained from the Paleobiology Database (https://paleobiodb.org; accessed July 13-26, 2026). Details of the MCMC settings and paleontological data-cleaning procedures can be found in the main text and Supplementary Appendix 3 of Lim et al. (2026).

Input files

Tabular input files are tab-separated text with a header row. <clade> stands for the folder name in lowercase.

config_<clade>.cfg — The UFBD configuration of the analysis: input file names, output prefix, MCMC settings, the conditioning of the FBD process and the prior on the origin age, the time bins and each rate, and the substitution and clock models. See UFBD README for details.

<clade>.fossils — One row per fossil species or unsampled extant species.

Column	Meaning
`taxon`	Species name, written `Genus_species`.
`min_age`, `max_age`	Youngest and oldest possible age of the fossil (Ma). The true age is treated as unknown within this interval. A row with both ages 0 is an extant species that is not in the backbone tree (an unsampled extant).
`clade`	Clade the species belongs to, as defined in `<clade>.clades`. `whole` is the entire study group.
`assignment`	Where the species can attach to the backbone tree, relative to its clade. `CROWN`: among the crown group of the clade. `STEM`: on the stem branch of the clade, rootward of that ancestor. `TOTAL`: either.

<clade>.clades — Clade definitions. clade is the name used in the fossils file. taxa is a comma-separated list of backbone tip names; the clade is every tip descended from the most recent common ancestor of the listed tips, so two tips suffice.

<clade>.treefile — The backbone tree of the extant species in Newick format, rooted and bifurcating, without branch lengths. The topology is held fixed, and the node ages are estimated.

in.BV (Cetacea) — Maximum-likelihood branch lengths, gradient and Hessian of the log-likelihood for each of the six sequence partitions of McGowen et al. (2020), in the in.BV format of PAML (Yang, 2007). It replaces the alignment in the approximate-likelihood analysis. Each block starts with the number of taxa (72), then the tree with branch lengths, the gradient vector, and the Hessian matrix. Blocks are numbered 0–5 in file order, and this is the numbering of the clock partitions in the output.

emydidae.fasta (Emydidae) — Aligned DNA sequences of the 50 backbone species, 13,659 sites, FASTA format, - for gaps and missing data.

emydidae_partition.nex (Emydidae) — NEXUS charset block giving the site range of each of the 15 loci, named by gene. Loci are numbered 0–14 in this order in the output.

Output files (`output/`)

Each analysis ran four independent MCMC chains numbered 0–3. Files named <prefix>_chain<N>.* belong to chain N and include the burn-in. Files named <prefix>.* pool the four chains after discarding the first 25% of each as burn-in. <prefix> is the clade name in lowercase.

The Emydidae folder holds two analyses. emydidae uses the non-centered parameterization of the clock-rate variance and is the analysis reported in the paper. emydidae_centered is the same model under the centered parameterization, run to compare the mixing of the two. The centered run was stopped after about 191 million generations without reaching the convergence targets, so it has per-chain files only, and its console file covers the part of the run after a restart from a checkpoint at about generation 47.4 million.

<prefix>.console.txt — Everything UFBD printed to the screen: the full command line, the time interval of every rate column (under Rate intervals:), the progress of each chain, the acceptance rates of every proposal, and the generation at which the convergence targets were reached.

<prefix>_chain<N>.log, <prefix>.log — The MCMC samples, one row per sampled generation. Column names are built from the model, so the three analyses have different columns.

Column	Meaning
`n`	MCMC generation of the sample. In the pooled `<prefix>.log`, a running counter.
`posterior`, `likelihood`, `prior`	Natural logarithm of the posterior density, the likelihood, and the prior density.
`x0`	Age of the origin of the clade (Ma), the start of the stem branch leading to the root.
`x1` … `x<K>`	Ages of the K internal nodes of the backbone tree (Ma). `x1` is the root. The same labels mark the nodes in the `.trees` and `.tree` files.
`originAge`	Origin age (Ma) in the Palaeodictyopterida analysis, which has no backbone tree.
`lambda`, `mu`, `psi*`	Speciation, extinction and fossil sampling rates (per lineage per Myr). With one rate for all time (Emydidae), the columns are `lambda`, `mu`, `psi`. With time bins, the suffix names the bin: `lambda0` … `lambda8` in Cetacea for the bins [0,5), [5,10), …, [35,40), [40,∞) Ma; `lambda_b0` … `lambda_b6` in Palaeodictyopterida for [250,260), [260,270), …, [300,310), [310,∞) Ma. Likewise for `mu` and `psi`. The `Rate intervals:` block of the console file lists this map.
`ou_lambda_theta`, `ou_lambda_sdEq`, `ou_lambda_nu`, and the same for `mu`	Hyperparameters of the Ornstein–Uhlenbeck process that smooths the log-rates across bins: long-term mean of the log-rate, stationary standard deviation of the log-rate, and reversion rate (per Myr).
`nSA`	Number of fossils that are sampled ancestors in the sample, that is, placed on a branch as a direct ancestor of a later sampled species instead of on a side branch of their own.
`clockMean_<k>`, `clockSigma2_<k>`	Mean substitution rate (substitutions per site per Myr) and variance of the log branch rates in clock partition k. Cetacea: k = 0–5, the six blocks of `in.BV`, geometric Brownian motion clock. Emydidae: k = 0–14, the 15 loci, uncorrelated lognormal clock.
`exch<k>_0` … `exch<k>_5`, `freq<k>_0` … `freq<k>_3`, `alpha<k>` (Emydidae)	GTR substitution model of locus k. The six exchangeability rates in the order A–C, A–G, A–T, C–G, C–T, G–T (sum to 1); the equilibrium frequencies of A, C, G, T (sum to 1); and the shape of the gamma distribution of rate variation among sites, with four rate categories.

<prefix>_chain<N>.trees, <prefix>.trees — The sampled backbone trees in NEXUS format, one per row of the corresponding .log; TREE STATE_<n> matches column n. Branch lengths are in Myr and internal nodes carry the labels x1 … x<K>. Fossils are not drawn; their positions are in the latent log. Absent for Palaeodictyopterida.

<prefix>.tree — The posterior summary tree in NEXUS format. Each internal node sits at its posterior mean age, and the annotation [&height_95%_HPD={a,b},age_95%_HPD={a,b}] gives the 95% highest posterior density interval of that age (Ma). height is measured from the youngest tip and age from the present; with every tip extant, they coincide. Absent for Palaeodictyopterida.

<prefix>_chain<N>_latent.log — The latent variables of each sample, one row per row of the chain's .log. Not pooled.

Column	Meaning
`n`	MCMC generation, as in the `.log`.
`y_<taxon>`	Age of the fossil (Ma), sampled within its `min_age`–`max_age` interval.
`z_<taxon>`	Attachment age (Ma): the age at which the lineage of that fossil or unsampled extant species branches off the backbone tree.
`sa_<taxon>`	1 if the fossil is a sampled ancestor in that sample, 0 otherwise.
`az_<taxon>`	Attachment zone of the fossil, an integer decoded by `<prefix>_chain<N>_zones.tsv`. Written only for fossils that can attach in more than one zone.
`rate_<node>_p<k>`	Substitution rate (substitutions per site per Myr) of the branch leading to `<node>` in clock partition k, where `<node>` is a tip name or an internal-node label `x<i>`.

<prefix>_chain<N>_zones.tsv — Legend for the az_ columns. Each zone is one clade–assignment pair that occurs in the fossils file: zone_index is the integer written in the latent log, clade is the clade name, and type (crown, stem or total) is the assignment. Identical across chains. Absent for Palaeodictyopterida.

<prefix>_chain<N>.log.ckp — Checkpoint of the full state of the chain (generation, random-number-generator state and all parameter values), used to restart an interrupted run with ufbd ... -resume. Not meant to be read directly.

References

Lim, W., Raskin, L. Y., Li, J. K., Huelsenbeck, J. P., & Nielsen, R. (2026). Estimating divergence times and diversification rates with the unresolved fossilized birth-death process. bioRxiv. https://doi.org/10.64898/2026.09.10.749957

McGowen, M. R., Tsagkogeorga, G., Álvarez-Carretero, S., dos Reis, M., Struebig, M., Deaville, R., Jepson, P. D., Jarman, S., Polanowski, A., Morin, P. A., & Rossiter, S. J. (2020). Phylogenomic resolution of the cetacean tree of life using target sequence capture. Systematic Biology, 69(3), 479–501. https://doi.org/10.1093/sysbio/syz068

Thomson, R. C., Spinks, P. Q., & Shaffer, H. B. (2021). A global phylogeny of turtles reveals a burst of climate-associated diversification on continental margins. Proceedings of the National Academy of Sciences of the United States of America, 118(7), e2012215118. https://doi.org/10.1073/pnas.2012215118

Yang, Z. (2007). PAML 4: Phylogenetic analysis by maximum likelihood. Molecular Biology and Evolution, 24(8), 1586–1591. https://doi.org/10.1093/molbev/msm088