Dear Placide, Great analysis.
We replicated your results using a very different maximum likelihood phylodynamics pipeline and different software. We are sharing here, since the replication using very different methods should give more confidence about your main findings.
The full reproducible analysis is here: BDBV 2026 outbreak: phylogenetic dating analysis
Note that this uses a slightly different selection of data than in your analysis.
Here is a brief summary:
The ML analysis is in close agreement on the TMRCA: 18 March (Bootstrap 95%CI: 2 Feb - 17 Apr).
This figure also shows parsimony reconstruction of location at interior nodes. Relatively few of the internal nodes are well-resolved, so I would be cautious about phylogeography for the moment.
Our analysis slightly disagrees on the evolutionary rate: We get 6.3*10^{-4} subst/site/year, but with very wide CI (Bootstrap 95% CI: 3.4-13.4 * 1e-04).
We also estimated the rate using the Bayesian BactDating method, and that rate agrees with our ML analysis. Disagreement on rates is common between different models and software and I find that BEAST estimates tend to be on the higher end. Both of our values are within the range that has previously been reported for EBOV. For now, I suggest checking sensitivity to the tree prior. The discrepancy may also be related to us rooting the tree using outgroups rather than via the clock model.
A profile likelihood is a different way of looking at the CI and gives something a bit tighter and shows how noisy the optimisation is:
We used a more complex model for detecting outliers that adjusts for multiple testing bias (treedater::outlierLineages). Using this criterion, the two samples you flagged don’t meet the threshold for exclusion. I don’t think it hurts your results, but we chose to keep them.
Relatedly, we applied the new DiagnoDating methods for to generate some diagnostic plots. We believe this flags the same samples you identified, but none of them meet threshold if adjusting for multiple testing, and the sample distributions (rates and residuals) look quite close to expected.
There is some clear overdispersion in the QQ plot using the strict clock model; this disappears with the additive relaxed clock, so this would support switching to relaxed clocks already.
We performed a formal test to choose between strict and relaxed clock models (treedater::relaxedClockTest). It’s very marginal, but this still prefers the strict clock. This will likely flip with the collection of more data. In any case, using an additive relaxed clock doesn’t move the dates very much (March 13 TMRCA vs March 18).
We estimated Ne(t) using ML and a parametric bootstrap for confidence intervals. This looks similar to your estimates but we haven’t compared quantitatively. It’s obviously important to know if the saturation in Ne observed in April reflects underlying epidemiology or sampling. I would treat this very cautiously for the moment, since I don’t know how the selection of samples was made. We can say that there isn’t any obvious geographical skew in recent samples that would explain the drop in Ne. Geographic diversity of the samples did not really drop April-> June.
Erik Volz & Vinicius Franceschi, Imperial College London