Widespread bursts of diversification in microbial phylogenies

TL;DR

This paper introduces a new theoretical framework to analyze microbial diversification, revealing widespread rapid bursts of evolution across diverse habitats and uncovering universal patterns in microbial evolutionary tempos.

Contribution

It presents a novel coarse-grained approach to quantify ecological innovations and diversification bursts in microbial phylogenies, robust to data incompleteness.

Findings

Identification of widespread diversification bursts in microbial evolution

Discovery of universal diversification tempo across habitats

Robust methodology applicable to incomplete microbial data

Abstract

Until recently, much of the microbial world was hidden from view. A global research effort has changed this, unveiling and quantifying microbial diversity across enormous range of critically-important contexts, from the human microbiome, to plant-soil interactions, to marine life. Yet what has remained largely hidden is the interplay of ecological and evolutionary processes that led to the diversity we observe in the present day. We introduce a theoretical framework to quantify the effect of ecological innovations in microbial evolutionary history, using a new, coarse-grained approach that is robust to the incompleteness and ambiguities in microbial community data. Applying this methodology, we identify a balance of gradual, ongoing diversification and rapid bursts across a vast range of microbial habitats. Moreover, we find universal quantitative similarities in the tempo of diversification, independent of habitat type.

Tables

Abbreviation	Short study description	Qiita ID	Accession/study link
AlaFir	Alaskan Fire Chronosequence - Tanana Valley	1030	http://www.ebi.ac.uk/ena/data/view/PRJEB14866
AldFir	Alder/Fir	1031	http://www.ebi.ac.uk/ena/data/view/PRJEB15055
BerMes	Bergen Mesocosm	1222	http://www.ebi.ac.uk/ena/data/view/PRJEB14793
BioFun	Biodiversity And Functional Patterns Of Microbial Assemblages In Postglacial Pond Sediment Profiles	1622	http://www.ebi.ac.uk/ena/data/view/PRJEB14823
BioShr	Bioturbating Shrimp Alter The Structure And Diversity Of Bacterial Communities In Coastal Marine Sediments	678	http://www.ebi.ac.uk/ena/data/view/PRJEB15461
BraAnt	Brazilian Antarctic Cleanup	1033	http://www.ebi.ac.uk/ena/data/view/PRJEB14907
CanSoi	Cannabis Soil Microbiome	1001	http://www.ebi.ac.uk/ena/data/view/PRJEB15461
CatArc	Catlin Arctic Survey 2010	723	http://www.ebi.ac.uk/ena/data/view/PRJEB18098
CatSou	Catchment Sources Of Microbes	894	http://www.ebi.ac.uk/ena/data/view/PRJEB14739
ComMic	Comparison Of Microbial Flora In Ant-Eating Mammals	1056	http://www.ebi.ac.uk/ena/data/view/ERX301013
DevOra	Development Of The Oral Microbiota In Captive Komodo Dragons (Varanus Komodoensis)	1747	http://www.ebi.ac.uk/ena/data/view/PRJEB14602
EffSoi	Effect Of Soil Ph On Soil Metagenome	805	http://www.ebi.ac.uk/ena/data/view/PRJEB18597
EmpNzt	EMP Nztabs	1035	http://www.ebi.ac.uk/ena/data/view/PRJEB3228
EpoSva	EPOCA Svalbard	—	https://doi.pangaea.de/10.1594/PANGAEA.769833
FerSpa	Fermilab Spatial Study	—	—
FriAla	Friedman Alaska Peat Soils	1692	http://www.ebi.ac.uk/ena/data/view/PRJEB15218
GeoLan	Geochemical Landscapes	1036	http://www.ebi.ac.uk/ena/data/view/PRJEB14909
GleCan	Glen Canyon Soils	1526	http://www.ebi.ac.uk/ena/data/view/PRJEB15163
GloPat	Global Patterns Of 16S Rrna Diversity At A Depth Of Millions Of Sequences Per Sample	721	http://www.ebi.ac.uk/ena/data/view/PRJEB18917
GreIgu	Green Iguana Hindgut Microbiome	963	http://www.ebi.ac.uk/ena/data/view/PRJEB15058
GreLak	Great Lake Microbiome	1041	http://www.ebi.ac.uk/ena/data/view/PRJEB14819
GulOil	Gulf Oil Spill Sediment	1197	http://www.ebi.ac.uk/ena/data/view/PRJEB14900
GutMic	Gut Microbiota Of Phyllostomid Bats That Span A Breadth Of Diets	1734	http://www.ebi.ac.uk/ena/data/view/PRJEB14489
HawKoh	Hawaii Kohala Volcanic Soils	1579	http://www.ebi.ac.uk/ena/data/view/PRJEB15174
IceWed	Ice Wedge Polygon	1578	http://www.ebi.ac.uk/ena/data/view/PRJEB9043
IntMic	Intertidal Microbes	662	http://www.ebi.ac.uk/ena/data/view/PRJEB18565
JurAnt	Jurelivicius Antarctic Cleanup	776	http://www.ebi.ac.uk/ena/data/view/PRJEB15611
KniCom	Knight Comp Biogeography	1748	http://www.ebi.ac.uk/ena/data/view/ERP022166
L4Tim	L4 Time Series 2009-2010	1240	http://www.ebi.ac.uk/ena/data/view/PRJEB14864
LatSur	Latitudinal Surveys Of Algal-Associated Microorganisms	933	http://www.ebi.ac.uk/ena/data/view/ERP021699
LdrBio	Ldrd Biological Carbon Sequestration	1043	http://www.ebi.ac.uk/ena/data/view/PRJEB14912
Lts16S	Ltsp	1037	http://www.ebi.ac.uk/ena/data/view/PRJEB14908
MagMon	Magnificent Mongolian Microbes	864	http://www.ebi.ac.uk/ena/data/view/PRJEB15460
MenMic	Mendota Microbial Observatory	1242	http://www.ebi.ac.uk/ena/data/view/PRJEB14911
MetSan	Metcalf Sandiego Zoo Folivorus Primate	2182	http://www.ebi.ac.uk/ena/data/view/PRJEB14631
MicCom	Microbial Community Of The Bulk Soil And Rhizosphere Of Rice Plants Over Its Lifecycle	1642	http://www.ebi.ac.uk/ena/data/view/PRJEB15194
MisBay	Mission Bay Sediment Viromes	1673	http://www.ebi.ac.uk/ena/data/view/PRJEB15214
MonCow	Monensin Cow Hindgut Study Cornell	1621	http://www.ebi.ac.uk/ena/data/view/PRJEB14795
MorDry	Morgankiss Dry Valley Lake Communities Protist Diversity In A Permanently Ice-Covered Antarctic Lake During The Polar Night Transition	638	http://www.ebi.ac.uk/ena/data/view/PRJEB18567
MovPic	Moving Pictures Of The Human Microbiome	550	http://www.ebi.ac.uk/ena/data/view/PRJEB19825
OreTra	Oregon Transect	1038	http://www.ebi.ac.uk/ena/data/view/PRJEB14862
PavLak	Pavilion Lake Research Project	809	http://www.ebi.ac.uk/ena/data/view/PRJEB18097
PerSta	Peralta Starlings	1694	http://www.ebi.ac.uk/ena/data/view/PRJEB14798
PolPol	Polluted Polar Coastal Sediments	1198	http://www.ebi.ac.uk/ena/data/view/PRJEB14880
RioJan	Rio De Janeiro Coastline	1039	http://www.ebi.ac.uk/ena/data/view/PRJEB15044
RusTun	Russian Tundra Samples Cryocarb	1034	http://www.ebi.ac.uk/ena/data/view/PRJEB15045
SalEnv	Saline Environments That May Harbor Novel Lignocellulolytic Activities Tolerant Of Ionic Liquids	1580	http://www.ebi.ac.uk/ena/data/view/PRJEB15178
SoiMet	Soil Metagenome	808	http://www.ebi.ac.uk/ena/data/view/PRJEB18643
Son	Microbiota of freshwater fish slime and gut	940	http://www.ebi.ac.uk/ena/data/view/PRJEB14822
SpaTem	Spatial And Temporal Variation In Nest And Egg Bacteria Of Wild Birds	1098	http://www.ebi.ac.uk/ena/data/view/PRJEB15144
TemBog	Temperate Bog Lakes	1288	http://www.ebi.ac.uk/ena/data/view/PRJEB15148
TemTx	Temple Tx Native Exotic Ppt	1289	http://www.ebi.ac.uk/ena/data/view/PRJEB15146
YelGra	Yellowstone Gradients	925	http://www.ebi.ac.uk/ena/data/view/PRJEB20056

Equations

$$\text{Birth (B):}L\stackrel{b}{⟶}2L,$$

$$\text{Death (D):}L\stackrel{d}{⟶}\varnothing ,$$

$$\text{Innovation (I):}L\stackrel{\rho g_x(k)}{⟶}kL,$$

$$\text{Heterogeneous Innovation (H):}L\stackrel{\eta h_{s,m}(k)}{⟶}kL,$$

$$\text{(B):}{\mathcal{L}}_B(b)=b\left(z-1\right)z{\partial }_z,$$

$$\text{(D):}{\mathcal{L}}_D(d)=d\left(1-z\right){\partial }_z,$$

$$\text{(I):}{\mathcal{L}}_I(\rho ,x)=\rho x\frac{\left(z-1\right)z}{1-xz}{\partial }_z,\text{ and}$$

$$\text{(H):}{\mathcal{L}}_H(\eta ,s,m)=\eta m{}_2{F}_1\left(1,sm+1,s+1;z\right)\left(z-1\right)z{\partial }_z.$$

$$\log \text{Pr}\left[{\mathcal{T}}_K|\theta \right]=\sum _{i\in \mathcal{I}[{\mathcal{T}}_k]}\log {\varphi }^{(k_i)}(t_i,s_i|\theta )$$

$${\theta }^{\ast }=\arg \underset{\theta }{\max }\log \text{Pr}\left[{\mathcal{T}}_K|\theta \right].$$

$${\varphi }^{(k)}(t,s|\theta )=\frac{1}{k!}{{\partial }_y^k{\Phi }_{t,s}(y|\theta )|}_{y=0}.$$

$${\Phi }_{t,s}(y|\theta )=\frac{{\mathcal{U}}_{t-s}\left({\mathcal{U}}_s(0|\theta )+y\left(1-{\mathcal{U}}_s(0|\theta )\right)|\theta \right)-{\mathcal{U}}_t(0|\theta )}{1-{\mathcal{U}}_t(0|\theta )},$$

$${\mathcal{U}}_{\tau }(z|\theta )=\sum _{k=0}^{\infty }P\left(k|\tau ,\theta \right)z^k,$$

$${\partial }_{\tau }{\mathcal{U}}_{\tau }(z|\theta )=\mathcal{L}(\theta ){\mathcal{U}}_{\tau }(z|\theta )$$

$$D=2\left(\log \text{Pr}\left[{\mathcal{T}}_K|{\theta }_{alt}^{\ast }\right]-\log \text{Pr}\left[{\mathcal{T}}_K|{\theta }_0^{\ast }\right]\right)$$

$$G({\mathcal{T}}_K)=2\sum _{\sigma =0}^{K-1}\sum _{k\in {\mathcal{K}}_{\sigma }}{n}_k^{\sigma }\log \frac{n_k^{\sigma }}{N_{\sigma }{\varphi }^{(k)}(t_{\sigma },s_{\sigma })},$$

$${\psi }_t(z)=\sum _{n\ge 0}{p}_n(t)z^n\equiv \sum _{n\ge 0}{p}_n(t)|n⟩=|{\psi }_t⟩$$

$$p_n(t)=\left[z^n\right]{\psi }_t\left(z\right)=\frac{1}{n!}{\frac{{\partial }^n}{\partial z^n}{\psi }_t(z)|}_{z=0}\equiv \frac{1}{n!}⟨n|{\psi }_t⟩.$$

$$⟨m|n⟩=m!{\delta }_{m,n}$$

$$1=\sum _n\frac{|n⟩⟨n|}{n!}.$$

$${\psi }_t(z)=⟨\underline{z}|{\psi }_t⟩=⟨0|e^{z\hat{a}}|{\psi }_t⟩.$$

$$⟨\hat{\mathcal{O}}{⟩}_{\psi }=⟨\underline{1}|\hat{\mathcal{O}}|\psi ⟩=\sum _{mn}{O}_{mn}{p}_n,O_{mn}=\frac{1}{m!}⟨m|\hat{\mathcal{O}}|n⟩.$$

$$\left[{\partial }_z,z\right]=1.$$

$${\partial }_t{\psi }_t(z)=\mathcal{L}\left[z,{\partial }_z\right]{\psi }_t(z)$$

$${\psi }_t(z)=e^{\mathcal{L}[z,{\partial }_z]t}{\psi }_0(z)\equiv ⟨\underline{z}|e^{\mathcal{L}[{\hat{a}}^{†},\hat{a}]t}|{\psi }_0⟩.$$

$$e^{x{\partial }_z}\psi (z,{\partial }_z)e^{-x{\partial }_z}=\psi (z+x,{\partial }_z)$$

$$e^{xz}\psi (z,{\partial }_z)e^{-xz}=\psi (z,{\partial }_z-x),$$

$$e^{xz{\partial }_z}\psi (z,{\partial }_z)e^{-xz\partial z}=\psi (z{e}^x,e^{-x}{\partial }_z).$$

$${\psi }_t(1-f+fz)=\sum _{n\ge 0}{p}_n(t)(1-f+fz{)}^n=\sum _{n\ge 0}\left(\sum _{k\ge 0}{p}_{n+k}(t)\left(\genfrac{}{}{0ex}{}{n+k}{k}\right)(1-f{)}^k\right)f^n{z}^n.$$

$$A\stackrel{b}{⟶}2A$$

Full text

Widespread bursts of diversification in microbial phylogenies

Alice Doucet Beaupré1,2*$\ddagger$, James P. O’Dwyer2,3$\ddagger$,

1 Program in Ecology, Evolution, and Conservation Biology, University of Illinois, Urbana, Illinois, USA

2 Department of Plant Biology, University of Illinois, Urbana, Illinois, USA

3 Carl R. Woese Institute for Genomic Biology, University of Illinois, Urbana, Illinois, USA

$\ddagger$ These authors contributed equally to this work.

• doucetb2@illinois.edu

Abstract

Until recently, much of the microbial world was hidden from view. A global research effort has changed this, unveiling and quantifying microbial diversity across enormous range of critically-important contexts, from the human microbiome, to plant-soil interactions, to marine life. Yet what has remained largely hidden is the interplay of ecological and evolutionary processes that led to the diversity we observe in the present day. We introduce a theoretical framework to quantify the effect of ecological innovations in microbial evolutionary history, using a new, coarse-grained approach that is robust to the incompleteness and ambiguities in microbial community data. Applying this methodology, we identify a balance of gradual, ongoing diversification and rapid bursts across a vast range of microbial habitats. Moreover, we find universal quantitative similarities in the tempo of diversification, independent of habitat type.

Introduction

Large scale microbiome sampling and sequencing [1, 2, 3, 4] have documented global microbial diversity with unprecedented scope and resolution. The tools currently applied to these data allow us to quantify the amount and type of diversity found in microbial communities [5, 6, 7], and yet we know remarkably little about the the underlying community dynamics and tempo of diversification that generated the biodiversity we observe. This gap in our knowledge calls out for robust new ecological and evolutionary theories that will allow us to connect mechanisms to observed patterns [8].

To address this challenge, we introduce a new methodology to bridge the gap between biological process and observed microbial biodiversity. Our approach leverages the inference of dynamical processes from evolutionary trees, previously applied to understand large-scale evolutionary structure[9, 10, 11, 12, 13, 14, 15], and also the dynamics of viral populations on shorter timescales[16, 17, 18, 19]. We also incorporate the recent identification of bursts of diversification in microbial phylogenies [20]. The result is a model which includes traditional, slow processes for gradual speciation (one lineage goes to two lineages, which we call the ‘birth’ of a lineage) and extinctions (one lineage disappears, which we call ‘death’), together with a third set of mechanisms, incorporating the process of ecological innovation potentially followed by radiative diversification.

We apply our framework to data spanning 13,500 individual samples, 56 habitat types, and 29 biomes [1], finding a previously unidentified balance of fast and slow evolutionary processes in these data, and a tendency towards universal behaviour in the quantitative description of diversification. We cannot directly quantify the traits and their changes through time that may have led to a given combination of rapid and gradual processes, but our results are strongly suggestive of a centre-ground in the long-standing debate over phyletic gradualism versus punctuated equilibrium [21].

Materials and methods

0.1 Coarse-grained phylodynamics

Our knowledge of the diversification of a group of organisms is often characterised by the branching of their evolutionary lineages, reconstructed using genetic sequence data sampled in the present day. We can think of an evolutionary tree, also known as a phylogeny, as the input data for inferring theoretical models of diversification, a methodology known as phylodynamics [16, 12, 13, 9, 14, 15, 18, 19, 10, 11]. This approach is usually thought of in terms of speciation and extinction rates. But in recent work [20] we identified patterns of bursts in the branching in microbial phylogenies. Therefore, in addition to the traditional speciation-extinction process, we also parametrize fast (but brief) bursts of diversification. In these burst processes, we propose that a lineage will undergo a much faster rate of diversification, $\sigma$, but for a very short time, $\tau$. The origin of a burst could e.g. be a key functional change that opens the opportunity for a rapid radiation [22, 23], or a disturbance which opens up a new habitat to be colonized. In the following we will refer interchangeably to this third process as bursts or innovations.

There are two problems with inferring the parameters of this generalized innovation process when $\sigma$ is large. First, there will be parts of any reconstructed phylogeny where we may not have enough information in our sequence data to distinguish the true ordering of very fast branching events, as shown in Figure 1. Even if we did have longer sequences, there is always a speed limit on what kinds of process we can accurately infer from these data, making these rates hard to infer. Second, in any realistic evolutionary history we would expect many different rates $\sigma$, corresponding to the idiosyncrasies of individual events. In this new approach, we (partially) bypass both difficulties by applying a method of coarse-graining, where we decompose a sample phylogeny of age $T$ into $K$ slices of width $T/K$ (Fig. 2 Panel A). In a coarse-grained phylogeny we are no longer trying to resolve down to each binary split in the tree—what we have access to are the ‘chunks’ of diversification between time slices (Fig. 2 Panel B), which define equivalence classes of binary trees.

Surprisingly, there is a way to bypass this speed limit, by leveraging the distribution of sizes of these (apparent) bursts of branching. Even though we can’t resolve phylogenies down to the shortest timescales, this distribution still carries information about the parameters of the innovation process. The catch is that we cannot distinguish between different values of $\sigma$ and $\tau$ independently, as the distribution of burst sizes is a geometric distribution which only depends on the product of diversification rate and diversification time, $\sigma\tau$. By looking at the evolutionary history through a blurred lens, we therefore collapse a multi-parameter family of models into a single parameter, reminiscent of the loss of information under coarse-graining in physics—so that at a sufficiently coarse resolution, many different fine-scale models map onto the same effective theory.

The final steps in our pipeline (Fig. 2 Panels C and D) involve computing likelihoods for a given model of innovation using the chunks defined by coarse-graining. We considered three distinct models in this preliminary work: the traditional speciation-extinction process (abbreviated BD, for birth and death); our basic speciation, extinction and individual burst process (BDI), which assumes a single value of the product $\sigma\tau$; and finally what we call speciation, extinction and heterogeneous bursts (BDH), where we compound the innovation processes using a distribution over the product $\sigma\tau$. BDI and BDH offer an agnostic and parsimonious alternative to models with time-varying [9, 10, 11], trait-dependent [12, 13], and diversity-dependent [14, 15] models.

The likelihoods are then calculated via numerical solutions of a master equation, described below, for the probability that an initial lineage (at the beginning of any chunk) will have branched into $k\geq 1$ lineages at the end of the chunk conditional on each lineage having at least one extant descendant. Each chunk represents a properly weighted sum over many histories compatible with uncertain tree structures and therefore helps us bypass the need for an accurate estimate of branch lengths.

Microbial phylogenies usually span enormous amounts of evolutionary time with relatively short sequences and may suffer from a lack of phylogenetic signal sufficient to reconstruct early ancestral states, which in turn imposes an horizon deep in the tree beyond which topological inaccuracies become inevitable [24]. Furthermore, slices closer to the present contribute much more weight in the total likelihood than those from an otherwise noisy and uncertain past because they contain more chunks. For example, in Fig. 3 below, consider the intersection of the cumulative distribution and the y-axis equals the number of chunks in that slice. We see that the slice closest to the present (top-left panel) contributes at least an orders of magnitude more chunks than the earliest ones (panels in bottom two rows) where low phylogenetic signal potentially degrades the quality of the input phylogeny.

What happens if we apply this methodology to a perfect, fully-resolved tree? Reassuringly, we show in the Supplementary Information that in this limit and using a very large number of slices our inference using the BD model exactly recapitulates previous approaches [25, 26, 27, 28, 9]. On the other hand, in cases where we do have limited resolution due to short sequences, our method extends the current applicability of BD models by allowing the inference of rate parameters over incompletely-resolved phylogenies. We can then distinguish between our three nested models using likelihood ratio testing, and we also perform an exact goodness-of-fit test by comparing a given empirical tree to an ensemble of typical trees generated by the model and its parameter estimates constrained to the same empirical size and depth (see Supplementary Information and Fig. S4).

0.2 Data availability

All amplicon sequence data and metadata have been made public through the data portal (qiita.microbio.me/emp) and all accession numbers to the European Nucleotide Archive (www.ebi.ac.uk/ena) can be found in Table S1. All processed individual phylogenetic trees used in this study will be made public prior to publication.

0.3 Code availability

Custom code for the tree coarse-graining, maximum likelihood inference, and Markov Chain Monte Carlo (MCMC) goodness-of-fit test is stored in a private GitHub repository and will be made public prior to publication.

0.4 Dataset and data preparation

As an initial application, we considered the Earth Microbiome Project (EMP) 10K (as of 2015), drawing from 5.59 million representative 16S rRNA sequences (open reference 97% similarity cut-off, V4 region[29]) gathered from 56 different studies and breaking into 13,500 separate samples [1]. Table S1 lists these studies, together with a short description as it appears in the EMP metadata, and a clickable URL to their Qiita entries and accession number or study website. These samples span 29 biomes and 56 habitats, including host-associated habitats from 93 different types of host. The EMP dataset is available on the project website (www.earthmicrobiome.org) as single large phylogenetic tree reconstructed using FastTree [30, 31] and representative sequences for each taxonomic unit from all habitats through the Qiita data portal (qiita.ucsd.edu).

The first critical question to address is to what extent we should aggregate samples when applying our phylodynamic inference. Should all human-associated samples be aggregated, or all soil samples? We know that the ancestors of present-day organisms in each sample were unlikely to be co-located, experiencing the same environmental context for their entire evolutionary history. On the other hand, we don’t want simply to pool all samples together; microbial life from sufficiently different habitats is likely to have different underlying evolutionary and ecological processes. For this study we took the approach of inferring parameters sample by sample, but it would be straightforward to aggregate any particular group of sampled data. For example, questions pertaining to phylogeography may necessitate the use of a spatially hierarchical aggregation scheme. We began with the large EMP phylogeny containing all representative sequences and we transformed it into a chronogram using the method of mean path lengths (MPL) [32, 33] weighted by relative operational taxonomic unit abundances. We chose this method because likelihood and semi-likelihood approaches for ultrametrizing phylogenies do not scale well for trees with more than 1,000 to 10,000 branches or leaves. MPL on the other hand requires only two traversals of the tree and therefore scale linearly with the size of the tree—ultrametrizing the full EMP phylogeny with MPL takes about 5 minutes on a consumer laptop. Alternatively supertree methods[34] could also be used for transforming in a reasonable amount of time very large phylogenies into a chronogram.

This second step is necessary because, like all phylodynamic approaches, our models are based on changes through time, whereas our reconstructed trees are based on sequence divergence. Transforming the tree into a chronogram puts the two on the same footing. We did not need to calibrate the root of the full EMP tree in actual Myr or Byr because only the dimensionless rates and exponents matter for the quantitative comparison and characterization of gradual (slow) vs bursty (fast) evolution. From this large phylogeny we then pruned a smaller subtree for each individual sample. These individual sample phylogenies form the input data for our analyses.

0.5 Model definitions

All models considered in this study (BD, BDI, and BDH) are constructed from four building blocks which create and destroy lineages which we will represent by the letter $L$, namely the stochastic birth process

[TABLE]

with per-lineage birth/speciation rate $b$, the death process

[TABLE]

with per-lineages death/extinction $d$, the innovation process

[TABLE]

where $\rho$ is the per-lineage innovation initiation rate $\rho$ and transition probabilities $g_{x}(k)=(1-x)x^{k-1}$ for $k\geq 1$ and 0 otherwise. The parameter $x=1-e^{-\sigma\tau}$ characterizes the geometric burst-size distribution. Finally we compound the innovation process with a beta distribution over $x$. The beta distribution is parametrized by its precision $s$ and mean $m$ to obtain the heterogeneous innovation process

[TABLE]

where $\eta$ is the per-lineage heterogeneous innovation initiation rate and transition probabilities $h_{s,m}(k)=s(1-m)\frac{\Gamma(s)\Gamma(sm+k-1)}{\Gamma(sm)\Gamma(s+k)}$ for $k\geq 1$ and 0 otherwise. The beta distribution is sometimes parametrized by two shape parameters, $\alpha$ and $\beta$, in our case controlling the behavior of the parameter $x$, now a random variable, around 0 and around 1, respectively. In terms of those shape parameters the precision $s=\alpha+\beta$ and mean $m=\alpha/(\alpha+\beta)$.

For each process there is an associated stochastic generator which encodes the instantaneous transition rates between different number of lineages presented in equations (1) through (4). Generators for the above processes are respectively given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

To combine models one adds generators together. For example the generator of the BDH model is written $\mathcal{L}_{BDH}=\mathcal{L}_{B}+\mathcal{L}_{D}+\mathcal{L}_{H}$ with associated parameter set $\theta_{BDH}=\left\{b,d,\eta,s,m\right\}$.

0.6 Coarse-grained phylogeny likelihood

The coarse-graining/slicing operator $\mathcal{R}_{K}$ decomposes a sample phylogeny $\mathcal{T}$ of age $T$ into slices of width $T/K$. The resulting coarse-grained phylogeny is characterized by a multiset of chunks $\mathcal{T}_{K}=\mathcal{R}_{K}\left[\mathcal{T}\right]=\{(t_{i},s_{i},k_{i})\}_{i\in\mathcal{I}[\mathcal{T}_{k}]}$, with $\mathcal{I}$ the index map given some arbitrary tree traversal. The log-likelihood of observing $\mathcal{T}_{K}$ under a given models with parameter set $\theta$ is written

[TABLE]

and maximum likelihood parameter estimates

[TABLE]

0.7 Chunk likelihoods

The expression for individual chunk likelihoods is given by

[TABLE]

The chunk generating function

[TABLE]

where $\mathcal{U}_{\tau}(z|\theta)$ is the probability generating function

[TABLE]

solution to the master equation

[TABLE]

with initial condition $U_{0}(z|\theta)=z$. The generator $\mathcal{L}(\theta)$ is one of $\mathcal{L}_{BD}(b,d)$, $\mathcal{L}_{BDI}(b,d,\rho,x)$, or $\mathcal{L}_{BDH}(b,d,\eta,s,m)$. The subtraction of $\mathcal{U}_{t}(0|\theta)$ in the numerator and normalization by $1-\mathcal{U}_{t}(0|\theta)$ in the denominator account for conditioning on non-extinction of at least one observed extant lineage. Details on the derivation, generalization to incomplete sampling of lineages, and numerical computation of the above quantities can be found in the Supplementary Information.

0.8 Model comparison and goodness of fit

The BD model is nested into BDI, and BDI into BDH, therefore we can perform a likelihood ratio test by comparing the $D$-statistic

[TABLE]

against a $\chi^{2}_{\text{ddof}}$ distribution with number of degrees of freedom $\text{ddof}=|\theta^{*}_{alt}|-|\theta^{*}_{0}|$. For $(\mathcal{H}_{alt},\mathcal{H}_{0})=(\text{BDI},\text{BD})$ we have $\text{ddof}=2$, and for $(\mathcal{H}_{alt},\mathcal{H}_{0})=(\text{BDH},\text{BDI})$ we have $\text{ddof}=1$.

We perform an exact goodness of fit test by sampling the $G$-statistic over the space of possible input coarse-grained phylogenies constrained by the depth and number of leaves of the empirical sample phylogeny. This statistic is given by

[TABLE]

where the first sum runs over slices and the second one over all chunk sizes found in slice $\sigma$. In the summand $n^{\sigma}_{k}$ represent the number of chunks of size $k$, and $N_{\sigma}=\sum_{k}n^{\sigma}_{k}$ the total number chunks, in slice $\sigma$. The $G$-statistic corresponds to the information divergence between the empirical and the theoretical chunk size distribution. The exact goodness of fit statistic is given by the fraction of coarse-grained phylogeny with $G(\mathcal{T}^{\prime}_{K})$ greater than the empirical value $G(\mathcal{T}_{K})$ where $\mathcal{T}^{\prime}_{K}$’s are sampled using the MCMC Metropolis-Hasting algorithm poised at the maximum likelihood parameter estimates. We describe the algorithm and coarse-grained proposal distribution in the Supplementary Information. This step is necessary because the number of degrees of freedom of trees with fixed depth, fixed number of leaves, and fixed number of slices is unknown.

Results

We rejected the BD hypothesis in favour of innovation (the BDI or BDH model) according to a likelihood ratio test for nested models (BD $\subset$ BDI $\subset$ BDH) at significance level $p~{}=~{}2\times 10^{-11}$, corresponding to a Bonferroni corrected $5\sigma$ family-wise error rate level $\alpha~{}=~{}2.7\times 10^{-7}$. In cases where BD was rejected, we use the maximum likelihood estimates for the parameters of the BDI or BDH process to analyze the balance of slow and fast processes (birth rates vs innovation rates) and the phylogenetic signature of fast processes themselves (the distribution of burst sizes). Figure S1 shows that BD was rejected in favour of BDI in 98% of samples. Subsequently, our basic innovation model was rejected in favour of heterogeneous innovation in 80% of samples. To understand why BDH is clearly selected in many cases it is instructive to look at a typical example. Figure 3 shows how BDH captures the fatter tail of the empirical chunk size distribution across all slices of a coarse-grained phylogeny, while BD (and to an extent BDI) fail. Nonetheless, Fig. S4 shows that for 57% of samples BDI and BDH are sufficient to recapitulate the phylogeny, while for 43% of samples both BDI and BDH fail the goodness-of-fit test, which suggests that we need a more complex model of innovation to account for them.

Beyond looking at single studies, our results can be expressed in terms of three important messages. First, the mixture of bursty (in the form of innovation events) and gradual diversification is clearly preferred over pure gradualism in the vast majority of our samples. To explain the dynamics implied by reconstructed evolutionary trees, we need both ongoing, slow diversification, and bursts of faster diversification that last for a relatively short time. Second, even though fast processes by definition produce more lineages per unit time while they are in play, the initiation of bursts of any size is also more common than slow gradualism in samples that show evidence of heterogeneous innovation. Finally, in the right hand panel of Fig. 4 and in Fig. S2-S3, we document the distribution of parameters controlling the shape of the burst size distribution for the heterogeneous innovation process. The distribution is beta-geometric and at large burst size $k$ this distribution behaves as $\sim k^{-(\beta+1)}$ where $\beta=s(1-m)$. The effective exponent $\beta+1$ of this power law is clustered with median $3.63$ and quartile coefficient of dispersion $(Q_{3}-Q_{1})/(Q_{3}+Q_{1})=0.07$, a surprisingly narrow range of values, independent of the other estimated parameter values, the habitat, or the study. The apparent universality has echoes in recent work [17, 20], and in the long history of studying scaling in evolutionary history, for example in the number of species per genus in a given taxonomic group [36]. Our current analysis goes beyond documenting these patterns, by connecting this universality to a mechanistic interpretation.

Discussion

We have introduced a new methodology to interpret what diversity in environmental sequence data can tell us about the ecological and evolutionary processes that shaped it. The key theoretical step in our new method is to recognize that faster diversification processes which appear intermittently and last only for a short time, still leave a signature in imperfectly reconstructed phylogenies. This signature persists even when the quality and length of our sequence data and consequent resolution of the phylogeny is relatively low compared to the timescale of the processes. Combining this realization with existing methods for inferring slower gradual processes from phylogenies, we were able to quantify the balance of fast and slow processes, and the parameter values that best describe the structure and distribution of burst sizes. Our conclusions in applying this to a large dataset encompassing heterogeneous habitats are stark: we almost always need these heterogeneous faster processes to complement gradual diversification in order to explain these data, and the parameters that best explain observed structure of evolutionary trees are surprisingly universal across studies and environmental context.

Conclusion

These results raise many questions, and open a number of doors for future investigation. Perhaps the primary open question is: what biological changes cause the bursts we observe in empirical trees? Are these genuinely due to innovations, where an adaptation opens the door to many further adaptations [22, 23, 37, 38]? They could also be the result of exploration of new habitats, disturbance opening up niche space to be invaded [39, 40], or something else entirely. Our current analysis cannot answer these questions clearly, but the evidence does clearly show that an explanation is necessary. Second, we have shown that a class of distinct fast processes all map on to the same observable phenomena at coarse temporal resolutions through a combination of their parameters. This is a quantitative example of a long-discussed idea in ecology that only a handful of parameters survive to describe phenomena at larger or longer scale. The assumption is inherent in neutral models, but also in other, simplified models of macroecological patterns [41]. Our approach can form the starting point of quantitatively understanding which parameters and processes ‘upscale’, and which do not. Finally, why do we see such clearly convergent patterns across divergent habitat types? The ecological and evolutionary constraints leading to the patterns we’ve seen deserve a fuller explanation.

Acknowledgments

A.D.B. was supported by a CompGen Fellowship from the University of Illinois and by the Cooperative State Research, Education, and Extension Service, US Department of Agriculture, under project number ILLU 875-952. J.O.D. acknowledges the Simons Foundation Grant #376199, McDonnell Foundation Grant #220020439, NSF #DEB1557192 and Templeton World Charity Foundation Grant #TWCF0079/AB47. Sample processing, sequencing and core amplicon data analysis were performed by the Earth Microbiome Project (www.earthmicrobiome.org), and all amplicon sequence data and metadata have been made public through the data portal (qiita.microbio.me/emp).

Author contributions

A.D.B. and J.O.D. designed the study and wrote the manuscript. A.D.B. developed analytical tools and analysed phylogenetic data.

Corresponding author

Correspondences to Alice Doucet Beaupré (doucetb2@illinois.edu)

Supplementary Information: Widespread bursts of diversification in microbial phylogenies

1 Formalism

1.1 Generating functions (GF) and holomorphic/Fock-space formalism

The equivalence between holomorphic/meromorphic generating functions and Fock-space methods applied to classical objects has a long history in the field of nonequilibrium statistical mechanics of many-body systems, beginning with the seminal papers [42, 43]. Yet it is only relatively recently that those methods have been recognized as potent tools applicable to mathematical biology and ecology[44, 45]. Let the probability generating function (PGF) of the state $\psi$ at time $t$

[TABLE]

where $p_{n}(t)$ represents the probability of having abundance $n$, or more to the point, of having $n$ lineages at time $t$. We adopt the convention that roman letters inside a ket $|n\rangle$ are used to denote a monomial/unit mass/abundance state $|n\rangle\equiv z^{n}$ while Greek letters inside a ket $|\psi\rangle$ denote mixtures like above. Single abundance state probabilities can be extracted by successive derivation

[TABLE]

In the Fock-space formalism we write the explicit scalar product

[TABLE]

and the completeness relation

[TABLE]

The normalization of the probability distribution $p_{n}(t)$ implies that $\psi_{t}(1)=1$. Evaluating a GF at a point $z$ requires the introduction of a left coherent-state, denoted with an underline,

[TABLE]

En passant this innocent looking equation is at the core of the equivalence between holomorphic functions and state vectors in Fock spaces. It allows us to easily translate between the natural, combinatorially intuitive language of generating functions, which we will use profusely in the following to construct our coarse-grained conditioned tree observables (or chunks), and the practical linear algebra methods forming the numerical backbone of this study. Under this equivalence the normalization condition becomes $\langle\underline{1}|\psi_{t}\rangle=1$. We use this particular coherent state to find the expectation value of operators, namely

[TABLE]

Classical stochastic observables are usually diagonal, i.e $O_{mn}=o_{n}\delta_{m,n}$. The holomorphic representation for creation and annihilation operators corresponds to multiplication and derivation by $z$, i.e. $\hat{a}^{\dagger}\equiv z$ and $\hat{a}\equiv\partial_{z}$, and satisfy bosonic commutation relation

[TABLE]

The master equation for a continuous-time stochastic process is written in the language as the partial differential equation (PDE)

[TABLE]

together with initial condition $\psi_{0}(z)$. $\mathcal{L}[z,\partial_{z}]$ is the time evolution generator and encodes all the information about the instantaneous dynamics of the process. Brackets here denote the dependence of the generator on $z$ and $\partial_{z}$, not its multiplication by the bosonic commutator. Eq. (23) admits the formal solution

[TABLE]

This is the starting point for the numerical exponentiation scheme used below and in the main text to construct the likelihood of a coarse-grained (CG) tree.

The holomorphic formalism admits three important similarity transformations[46]: two shifts,

[TABLE]

and

[TABLE]

and a scaling transformation,

[TABLE]

1.2 Incomplete lineage sampling

If we approximate the sampling process by a Bernouilli trial with success probability, or sampling fraction, $f$, then the GF for the joint probability of successfully sampling $n$ individuals out of a population of $n$ or more individuals is given by

[TABLE]

This expression captures the fact that every states with abundance greater than $n$ contribute to the probability of sampling exactly $n$ lineages.

2 Processes

2.1 Birth process

The birth process with per capita birth rate $b$ consists in the transition

[TABLE]

The generator for the birth process is given by

[TABLE]

The instantaneous transition rates between states are given by

[TABLE]

The generator, in matrix form,

[TABLE]

2.2 Death process

The death process with per capita death rate $d$ consists in the transition

[TABLE]

The generator of the death process is given by

[TABLE]

and the instantaneous transition rates

[TABLE]

[TABLE]

2.3 Innovation process

The innovation process consists in initiating at per capita rate $\rho$ a Yule (pure birth) process and identifying all finite-time transition probabilities with infinitesimal transition rates. To obtain the innovation process generator we first need to solve the birth process exactly. We begin by writing its formal solution using Eqs. 29 and 23 as

[TABLE]

Using the change of variable $y=1/z$ we rewrite the evolution equation

[TABLE]

Using Eqs. 25 and 27 we solve

[TABLE]

where $\alpha=1-e^{-\sigma\tau}$. For a given burst time-scale $\tau$, higher fitness processes lead to values of $\alpha$ closer to 1. For a single initial lineage, $\psi_{0}(z)=z$ and we recover the geometric PGF

[TABLE]

The generator of the innovation process

[TABLE]

where $g_{\alpha}(n)=(1-\alpha)\alpha^{n-1}$, is therefore given by

[TABLE]

In the second equality, we made sure to only include off-diagonal transitions $1\rightarrow k\geq 2$ in the positive term by absorbing the false transition $1\rightarrow 1$ in the negative diagonal term. Doing so also highlights the fact that the effective innovation rate is actually equal to $\rho\alpha$ rather than $\rho$. It is tempting to interpret, in line with the Red Queen hypothesis, the rate $\rho(1-\alpha)$ of initiations of “invisible” innovations as capturing perhaps the rate at which lineages must constantly innovate in short bursts of higher fitness just to keep pace with co-occurring lineages and changing environments. The instantaneous transition rates are given by

[TABLE]

and the matrix of the generator

[TABLE]

2.4 Heterogeneous innovation process

To account to the possibility of multiple innovation processes with different values of $\alpha$ we compound the parameter of the geometric part with its conjugate prior the beta distribution. The generator

[TABLE]

In the last two lines we used the reparametrization of the beta distribution in terms of the precision $s=\alpha+\beta$, $s>0$, and the mean $\overline{m}=\alpha/(\alpha+\beta)$, $0<\overline{m}<1$. Similarly as before the instantaneous rates

[TABLE]

For $k=m-n\gg 1$, instantaneous rates follow the power law

[TABLE]

The power law phase is controlled by the $\beta$ parameter of the beta distribution. The parameter $\beta$ in turn controls the shape of the density of the geometric innovation parameters around $x=1-\epsilon$. As $\epsilon\ll 1$ approaches 0 and $x$ approaches 1, the geometric distribution acquires a progressively longer exponential decay that tends towards a uniform improper distribution over all $k\geq 1$. When $\beta<1$, the density diverge algebraically as $\epsilon^{\beta-1}$ and the compounding of wide geometric distributions give rise to power laws with tail exponent between 1 and 2. When $\beta=1$, then the tail exponent becomes exactly 2, and when $\beta>1$, the tail exponent becomes greater than 2. Finally the matrix of the generator

[TABLE]

3 Maximum Likelihood approach

3.1 CG tree observables and likelihood

We extract dynamical information contained in a chronogram (ultrametric phylogenetic tree) by looking at CG transitions, i.e. single lineages going to many lineages at a time closer to the present. Those are called chunks in the main text and to each chunk is associated a tuple $(t,s,k)$ and likelihood chunk $\phi_{f}^{(k)}(t,s)$ where $0<f\leq 1$ is the sampling fraction. A likelihood chunk represents the probability that a unique lineage at time $t$ in the past had exactly $k\geq 1$ descendants at time $s\leq t$ that each had at least 1 extant lineage in the present, while all other lineages above $k$ were unobserved by virtue of having gone extinct or of being missed during sampling of intensity $f$. We can write using the previous formalisms

[TABLE]

where $\mathcal{U}_{t}(z|\theta)$ is the solution of Eq. 23 with initial condition $\mathcal{U}_{0}(z|\theta)=z$ and $\mathcal{U}^{(k)}_{t}(z|\theta)$ is its $k$-th derivative w.r.t. $z$. In terms of probabilities contained in $\mathcal{U}_{t}(z|\theta)=\sum_{k\geq 0}P(k|1,t,\theta)z^{k}$ and $\left(\mathcal{U}_{t}(z|\theta)\right)^{m}=\sum_{k\geq 0}P(k|m,t,\theta)z^{k}$ the expression 46 for chunks is given by

[TABLE]

This is a probability distribution over $k\geq 1$. To show this we sum over $k$

[TABLE]

In the first equality we added and subtracted the extinction term $k=0$. In the second equality we used the semigroup/Chapman-Kolmogorov property. In the third equality we used the shift property Eq. 25. In the fourth equality we cancelled terms inside the argument of $\mathcal{U}_{t-s}$ and finally used the normalization condition $\mathcal{U}_{t-s}(1)=1$. All terms are positive and sum to one and therefore $\phi^{(k)}_{f}(t,s|\theta)$ is a probability distribution over $k\geq 1$.

Combining the method of characteristics and numerical complex derivation, the first (holomorphic) formulation of Eq. 46 allows for the exact numerical computation of chunk likelihoods. Yet we found the second (Fock space) formulation using truncated matrix exponentials and linear algebra to be easier to implement and more stable at large values of $k$. All the information of a given model reside inside the generator $\mathcal{L}(\theta)$. The Greek letter $\theta$ is understood as the set of all parameters of a given model. We will denote the action of CGing/slicing with slice width $\Delta=T/K$ a chronogram $\mathcal{T}$ by

[TABLE]

The index set $I$ runs over chunk indices following some traversal order over the CG chronogram. We can write the likelihood of a given CG chronogram

[TABLE]

The central inference objective of this framework is to seek ML estimates (MLE) $\theta^{*}$ by maximizing the log-likelihood

[TABLE]

3.2 Equivalence with the Morlon et al. likelihood

It is worth mentioning here that Eqs. 46 and 50 recover the likelihood found in [9] in the limit of infinitesimal chunk duration $t-s=\delta t$. This limit is also equivalent to a particular slicing scheme where each branch gets a chunk of size $k=1$ with duration $t-s$ equal the length of the branch, and each node a chunk of size $k=2$ realized with an infinitesimal duration. Let us first show the equivalence between the infinitesimal and branch-node slicing. Consider two consecutive branch segments of length $t-s$ and $s-r$. Together they contribute

[TABLE]

to the full tree likelihood. Recall the Chapman-Kolmogorov semi-group property of generating functions which ensures $\mathcal{U}_{t-s}(\mathcal{U}_{s-r}(z))=\mathcal{U}_{t-r}(z)$. Therefore

[TABLE]

Substituting backward in Eq. 52 it follows that

[TABLE]

This means that the product of consecutive chunks of size $k=1$ compounds into one chunk of size $k=1$ with length equal to the sum of individual chunk lengths. This in turn implies the equivalence between the infinitesimal slicing of a tree and a scheme where each branch $i$ gets its own $k=1$ chunk with the same initial and final times $t_{i}$ and $s_{i}+\delta t$, while its immediate downstream node inherits a chunk of size $k=2$ with initial and final times $s_{i}+\delta t$ and $s_{i}$. Chunks “touching the present” contribute a weight

[TABLE]

Finally node contributions for the BD model in this scheme

[TABLE]

Therefore the contribution of an internal branch starting at $t$ in the past and ending with a node at $s$ has weight

[TABLE]

and the tree likelihood in the limit $K\rightarrow\infty$ for a tree of depth $t$ and $n$ leaves

[TABLE]

This is equivalent to Eq. (1) in [9] for the case of constant birth and death rates modulo a factor $\delta t^{(n-1)}$. Indeed for the BD model

[TABLE]

where the absorbing state/exctinction probability

[TABLE]

and $r=d/b$ the ratio of constant per capita death and birth rates. We already found the Yule limit previously when we constructed the innovation process. In the Yule limit $r\rightarrow 0$, therefore $p_{0}\rightarrow 0$ and $\alpha=p_{0}(t)/r\rightarrow 1-e^{-bt}$. Generalizing to time-varying rates is straightforward and the equivalence remains valid.

3.3 Chunk generating function

Using manipulations similar to those of the previous sections, we find the chunk PGF

[TABLE]

This expression has a very intuitive combinatorial interpretation that could have been guessed from the beginning and used as the starting point for the chunk decomposition. It symbolically stipulates that at $t-s$, for whichever state you are in, substitute each lineage represented by atoms $z$ with a Bernouilli trial of success probability equal to the probability of surviving across time $s$ to the present. In other words, lineages are ”preemptively” split at time $t-s$ into empty atoms (the monomial 1) weighted by the probability of not making it to the present, $\mathcal{U}_{s}(1-f)$, and unit atoms (the monomial $y$) weighted by the probability of making it to the present, $1-\mathcal{U}_{s}(1-f)$. We easily recover the extant chunk PGF when conditioning on survival by omitting the extinction term at $y=0$ and renormalizing. We find

[TABLE]

and finally it is straightforward, if tedious, to verify that

[TABLE]

recovers Eq. 46.

The chunk generating function also satisfy the Chapman-Kolmogorov property

[TABLE]

which we can verify using Eq. 62 and the Chapman-Kolmogorov property for $\mathcal{U}_{t}(z)$. This property of the chunk generating function translates at the level of likelihood chunks thus:

[TABLE]

where $\text{Comp}(k)$ and $\text{Comp}(\pi)$ are respectively the set of compositions of the integer $k$ and by a slight abuse of notation the set of compositions equivalent to an integer partition $\pi$, i.e. a composition is equivalent to a partition iff both multisets of their parts are equal. E.g. $(1,5,4,5)\sim(5,5,4,1)$, but $(1,5,4,5)\not\sim(6,5,4,2)$. $\text{Part}(k)$ represents the set of partitions of the integer $k$. Finally for $\pi$ a partition of $k$, the size of the set of its equivalent compositions is given by the multinomial coefficient

[TABLE]

with $m_{\pi}(i)$ the multiplicity of parts of size $i$ in $\pi$. Eq. 65 describes how to breakup a chunk into a sum-product of chunks of smaller or equal size, effectively enumerating elements of its equivalence class. Applying this decomposition recursively one can see how a coarse-grained phylogeny is equivalently indexed by a set of time slices (the $r$’s so to speak) together with an integer composition for each slice. Moreover it naturally gives rise to the Hastings ratio used in the Metropolis-Hastings MCMC algorithm described in Section 4.3 where coarse-grained phylogenies are sampled using the space of integer multipartitions.

3.4 Combined models, nestedness, and likelihood ratio test (LRT)

The Birth-Death model (BD) is obtained by adding the Birth/Yule (B) and Death (D) processes together to obtain the BD generator

[TABLE]

Similarly for the Birth-Death-Innovation (BDI) model

[TABLE]

and the Birth-Death-Heterogeneous innovation (BDH) model

[TABLE]

All three models are nested. BD is recovered from BDI in the limit $\rho\alpha\rightarrow 0$, and BDI is recovered from BDH in the limit $\overline{m}\rightarrow\alpha$, $s\rightarrow\infty$, and $\eta\rightarrow\rho$. We can therefore perform model selection using the LRT, i.e. comparing the statistic

[TABLE]

against the $\chi^{2}_{\text{ddof}}$ distribution with $\text{ddof}=|\theta_{\text{alternative}}|-|\theta_{\text{null}}|$ degrees of freedom. $\text{ddof}=2$ for BD vs BDI, and $\text{ddof}=1$ for BDI vs BDH.

3.5 Numerical exponentiation and truncation

Optimizing Eq. 51 requires the evaluation of Eq. 46, equivalently Eq. 47, for arbitrary values of $\theta$. To do so we use the completeness relation Eq. 19 to transform Eq. 46 into an explicit matrix-vector multiplication, e.g.

[TABLE]

The vector elements of the left- and right-most bra and kets are

[TABLE]

and

[TABLE]

The matrix elements of the combination operator

[TABLE]

The elements of the two matrix exponentials,

[TABLE]

and

[TABLE]

are obtained using Scipy’s $\mathtt{linalg.expm}$ function[47, 48] and truncated combinations of Eq. 31, 34, 41, and 45. Since those generator matrices are formally of infinite dimension we need an approximation in order to perform any kind of numerics. Similarly to the finite state projection (FSP) algorithm[49], but even more simply, we truncate by taking the $N\times N$ submatrix of $\hat{L}$ which includes transitions amongst states $0\leq m,n\leq N-1$. We omit to correct the diagonal terms and therefore do not compensate for the missing probability flow into states $n\geq N$. This approximation makes $\hat{L}$, and by extension $e^{\hat{L}s}$, non-stochastic and would eventually leak all the probability mass of the initial state. We found this approximation to nonetheless work very well in supercritical cases, e.g. $b>d$ in the BD model, as long as $N$ is chosen at least between $2\times$ or $3\times$ the largest chunk size $\max_{i\in I}k_{i}$. Almost all maximum likelihood (ML) inference fall inside the supercritical region as expected from the simple fact that most phylogenies grow. It also avoids the build up of probability mass into higher abundances state which is caused by the reflecting boundary inherent to a stochastic truncation (i.e. with corrected diagonal terms). In any case we found results to be insensitive to the choice of leaky vs. non-leaky approximations. Finally we put all pieces together and find for a given chunk

[TABLE]

where all sums are now between [math] and $N-1$.

3.6 Full Likelihood construction

In order to accelerate the evaluation of the complete log-likelihood for a CG phylogeny given by Eq. 51 we organize the construction of the likelihood around three observations:

• •

all chunks have identical $t-s=\Delta$ and therefore share the same column-vector elements $[e^{(t-s)\hat{L}(\theta)}]_{m1}$,

• •

all chunks in a given slice share row-vector elements for both $\sum_{m}\left(1-f\right)^{m}[e^{s\hat{L}(\theta)}]_{mn}$ and $\sum_{m}\left(1-f\right)^{m}[e^{t\hat{L}(\theta)}]_{mn}$,

• •

and slices share interfaces such that, going from the past to the present, one slice’s $s_{i}$ is the next slice’s $t_{i}$, leading to cancellation of all first two terms of the form $\log(1-x)$ in Eq. 77 except at the root and leaves.

Thus memoization of select function calls or explicit storage of certain parts of the above expression greatly reduces the number of matrix exponentiations and dot products necessary to complete one evaluation of the likelihood. We used Scipy’s $\mathtt{optimize.minimize}$ function with the L-BFGS-B[50] algorithm to perform the optimization in Eq. 51. We typically used 5 restarts, each one initiated at the best out of 10 random starting points chosen randomly from a centered normal distribution with variance 2 in transformed coordinates $\log b,\log d,\log\rho,\log\eta,\log s$, and $\operatorname{logit}~{}\overline{m}$. We also assume an approximate complete sampling of lineages and set $f=1$.

4 Goodness of fit test

4.1 Exact goodness-of-fit (gof) test

The LRT allows us to compare two models together. We would also like to assess the quality of the MLE without reference to an alternative model. This is the purpose of the gof. For this purpose we choose the $G$-statistic characterizing the information divergence between empirical and theoretical chunk frequencies. For a CG phylogeny $\mathcal{T}_{K}=\mathcal{R}_{K}[\mathcal{T}]$,

[TABLE]

The first sum runs over CG slices $\sigma$ and the second sum over all non-zero chunk sizes within that slice. $n^{\sigma}_{k}$ is the number of chunks of size $k$ in slice $\sigma$, and $N_{\sigma}$ the total number of chunks in slice $\sigma$. Finally $\phi^{(k)}_{f}(t_{\sigma},s_{\sigma})$ is the ML chunk distribution found in Eq. 46. Implementing an exact gof requires that we find the distribution of $G$ over all trees given the MLE of the model parameters and the constraint of a given tree depth with a given number of leaves, namely

[TABLE]

Doing the sum exactly would require exploring the space of all input data (trees) satisfying identical constraints. This would require varying all internal branch lengths, polytomy size/node degrees, and node positions in all possible ways while keeping the tree depth and size constant and is a combinatorially intractable task for even modestly sized trees. To approximate this intractable sum we will instead use a Metropolis-Hastings MCMC algorithm, which we describe next.

4.2 CG proposal distribution

Fortunately the CG representation of trees gives us a way to completely bypass the need to implement multiple types of proposal distributions (topological moves and metric moves) which would have to be followed by an expensive CGing step of the proposed tree. All those types of move can be subsumed into a single proposal distribution: the uniform random sampling of partitions. The intuition behind this proposal is as follows. If a node within a chunk is dragged across a slice interface towards the past, the chunk it was part of gets fragmented and the chunk further in the past into which the node moves increases in size. Moving a node across a slice interface towards the present has the opposite effect, namely it decreases the size of its chunk of origin and coagulates several chunks in its new slice. Topological moves restricted to a given slice simply reshuffle chunk order and has no effect on the chunk frequencies within a slice (for models with independent lineages). Arbitrary topological moves across slices are equivalent to a series of fragmentation and coagulation events and are therefore implicitly realized by the moves described above. Proposing a CG move will therefore consists of four steps:

• •

Choose a random chunk anywhere in the CG phylogeny except for the first slice which contains the leaves. Denote its size $k_{s}$.

• •

Choose (without replacement) $k_{s}$ random chunks in the previous slice (immediately closer to the present). Their sizes are indicated by the partition (multiset) $\pi_{s}=\{\lambda_{j}\}_{j=1}^{k_{s}}$. This is a partition of $k=\sum_{\lambda\in\pi_{s}}\lambda$. Denote the multiplicity of parts in $\pi_{s}$ by $m_{\pi_{s}}(\lambda)$, $\lambda\in\mathbb{N}^{+}$.

• •

Choose uniformly at random a partition $\pi_{t}$ of $k$ and let $k_{t}=|\pi_{t}|$ be the number of parts in $\pi_{t}$. Denote those parts $\lambda_{n}$ with $n$ running from 1 to $k_{t}$ and their multiplicities $m_{\pi_{t}}(\lambda)$, $\lambda\in\mathbb{N}^{+}$.

• •

Replace $k_{s}$ by $k_{t}$, and the chunks with sizes in $\pi_{s}$ by the chunks with sizes in $\pi_{t}$. The number of elements $k_{s}=|\pi_{s}|$ and $k_{t}=|\pi_{t}|$ are not necessarily equal (there is no conservation rule at the interface across iterations), but the sums of their parts are equal because they are both partitions of the same number $k$ (conservation rule within an iteration).

Choosing randomly which $k_{s}$ chunks to repartition followed by the sampling a new partition with an arbitrary number $k_{t}$ of parts seamlessly combines both metric and topological moves, maintains fixed the depth of the tree and the number of leaves, and satisfies the lineage conservation constraint whereby the sum of chunk sizes in a slice is equal to the number of chunks in the adjacent slice closer to the present. We used the excellent and fast algorithm by Arrita and DeSalvo[51, 52] to implement the uniform random generation of partition. We have the following algorithm:

4.3 Metropolis-Hastings MCMC algorithm

• •

Propose a move $\mathcal{T}^{\prime}_{K}=\mathcal{M}[\mathcal{T}_{K}]$,

• •

Accept move with probability $\alpha=\min\left[1,\frac{\text{Pr}[\mathcal{T}^{\prime}_{K}]}{\text{Pr}[\mathcal{T}_{K}]}\frac{g\left(\mathcal{T}_{K}|\mathcal{T}^{\prime}_{K}\right)}{g\left(\mathcal{T}^{\prime}_{K}|\mathcal{T}_{K}\right)}\right]$,

• •

Repeat last two steps and occasionally return a sample tree $\mathcal{T}_{K}$ and statistic $G(\mathcal{T}_{K})$ according to a random sweep schedule.

We burn half the chain and return 100 samples. A sweep is considered complete when the number of proposed moves reaches about twice the number of chunks in the CG tree. We found this rule of thumb sufficient to eliminate autocorrelation between samples. The Metropolis ratio

[TABLE]

While the uniform random sampling of partitions is symmetric in the space of partitions it is not symmetric in the space of trees and we need to introduce the Hastings ratio

[TABLE]

To a given partition $\pi$ corresponds several underlying CG trees and those are enumerated by the set of compositions equivalent to a given partition. For a partition $\pi=\{\lambda_{n}\}_{n}$ with $|\pi|$ parts the number of equivalent compositions is given by the multinomial coefficient $|\text{Comp}(\pi)|=|\pi|!/\prod_{n}m_{\pi}(n)!$. For example, take the partition $\pi=(2,1,1,1)$ of 5. This partition has $k=4$ parts and multiplicities $m(1)=3$, $m(2)=1$, and $m(\lambda)=0$ for all other $\lambda\not\in\{1,2\}$. Therefore $|\text{Comp}(\pi)|=4!/(1!3!)=4$. Indeed the set of equivalent compositions $\text{Comp}(\pi)=\{(2,1,1,1),(1,2,1,1),(1,1,2,1),(1,1,1,2)\}$.

The output of the MCMC is a set of IID trees $\mathcal{T}_{K}^{1}$, $\mathcal{T}_{K}^{2}$, $\mathcal{T}_{K}^{3}$, $\ldots$ distributed according the ML chunk size distribution of a model with all constraints satisfied, together with IID samples $G(\mathcal{T}_{K}^{1})$, $G(\mathcal{T}_{K}^{2})$, $G(\mathcal{T}_{K}^{3})$, $\ldots$ of the $G$-statistic. The Metropolis-Hastings algorithm effectively biases the MCMC walk in the space of CG trees toward trees that are in some sense ”typical” for their size and depth given the model under consideration with parameters at their MLE. The $p$-value of the gof is approximated by the fraction of values of $G$ from the MCMC that are larger than the empirical value $G(\mathcal{T}_{K})$ and quantifies how ”typical” the empirical tree is. An asymptotically ”perfect” tree, i.e. with an empirical chunk distribution that fits the theoretical chunk distribution exactly, would have a $G$-statistic of zero and $p$-value of 1. A ”typical” tree should have a $p$-value fluctuating around 0.5. The $G$-statistic of typical trees will distribute according to a $\chi^{2}$ distribution. We did not find a way to calculate the number of degrees of freedom of this distribution a priori, which would render moot the whole MCMC procedure presented above. Instead we can fit a $\chi^{2}$ distribution to the MCMC output knowing that the number of degree of freedoms $k=\langle G(\mathcal{T})\rangle_{\text{MCMC}}$. This gives a complementary way to find an asymptotic approximation to the $p$-value, e.g. in cases where the MCMC is computationally very expensive and one can only obtain a few samples of the $G$-statistic in reasonable time.

Fig S4 shows the distribution of the gof p-values from a sub-sample of size 366 of all samples. Using the asymptotic p-values, we see that around 57% of model fits pass the exact gof test at significance level 0.05. When we look at the gof in conjunction with the LRT between BDI and BDH, we see that BDH passes the gof and is favorably selected in 39% of samples. BDI is not rejected and passes the gof in 18% of samples. BDI is not rejected and fails the gof in 4% of samples. Finally, and maybe more interestingly, BDH is favorably selected but does not pass the gof in 39% of samples.