# Scale-invariant topology and bursty branching of evolutionary trees emerge from niche construction

^{a}Loomis Laboratory of Physics, University of Illinois at Urbana–Champaign, Urbana, IL 61801;^{b}Carl R. Woese Institute for Genomic Biology, University of Illinois at Urbana–Champaign, Urbana, IL 61801;^{c}Institute for Universal Biology, University of Illinois at Urbana–Champaign, Urbana, IL 61801

See allHide authors and affiliations

Contributed by Nigel Goldenfeld, February 6, 2020 (sent for review August 30, 2019; reviewed by Marcus W. Feldman, Joachim Krug, and Kim Sneppen)

## Significance

Phylogenetic trees describe both the evolutionary process and community diversity. Recent work, especially on bacterial sequences, has established that, despite their apparent complexity, they exhibit two unexplained broad structural features which are consistent across evolutionary time. The first is that phylogenetic trees exhibit scale-invariant topology. The second is that the backbones of phylogenetic trees exhibit bursts of diversification on all timescales. Here, we present a coarse-grained model of niche construction coupled to simple models of speciation that recapitulates both the scale-invariant topology and the bursty pattern of diversification in time. These results show, in principle, how dynamical scaling laws of phylogenetic trees on long timescales may emerge from generic aspects of the interplay between ecological and evolutionary processes.

## Abstract

Phylogenetic trees describe both the evolutionary process and community diversity. Recent work has established that they exhibit scale-invariant topology, which quantifies the fact that their branching lies in between the two extreme cases of balanced binary trees and maximally unbalanced ones. In addition, the backbones of phylogenetic trees exhibit bursts of diversification on all timescales. Here, we present a simple, coarse-grained statistical model of niche construction coupled to speciation. Finite-size scaling analysis of the dynamics shows that the resultant phylogenetic tree topology is scale-invariant due to a singularity arising from large niche construction fluctuations that follow extinction events. The same model recapitulates the bursty pattern of diversification in time. These results show how dynamical scaling laws of phylogenetic trees on long timescales can reflect the indelible imprint of the interplay between ecological and evolutionary processes.

Phylogenetic trees represent the evolutionary history of a group of organisms, usually constructed through an appropriate proxy such as the so-called 16SrRNA gene. This gene codes for a part of the translational machinery of the cell and, as such, is presumed to be an essential part of all cellular life. Canonically, this gene (or the 18S variant in Eukaryotes) is used to define operational taxonomic units (OTUs) that correspond to a generalized notion of species. By evaluating the similarity of DNA sequences, the timeline of speciation and subsequent evolution can be estimated. Nodes on a phylogenetic tree represent OTUs, with the external or leaf nodes being extant organisms that are observed and the internal ones being hypothetical organisms that are inferred based on the similarity and the embedded evolutionary process. When a tree is rooted, the top node, or the root, represents the inferred common ancestor of all nodes in the tree. The lengths of the branches of the tree correspond to nucleotide changes, and the timescale is set by a molecular clock assumption.

It is by no means mandatory that evolutionary history should be tree-like in its topology. Prior to the last universal common ancestor, the translational machinery evolved rapidly, and today’s canonical genetic code emerged. This code is not only universal, but also is nearly optimal in minimizing errors, in a sense that has been quantified precisely (1, 2). Simulations of the evolution of the genetic code indicate that it would have been highly unlikely for it to have these properties if the evolution had proceeded in a treelike manner; only with horizontal gene transfer of core translational machinery, and thus a network topology for the phylogeny, can a unique and optimal code evolve rapidly (3). Thus, the evidence suggests that there was a major evolutionary transition that occurred during the emergence of the translational machinery; prior to this transition, the concept of species did not exist in the canonical sense. This transition can be inferred from a topology change in the phylogeny of the 16SrRNA gene. In short, there is much that can be learned about evolution by studying the topology of its timeline.

In this paper, we ask what can be learned from a detailed study of the topology of modern phylogenetic trees, constructed with the benefit of large-scale genomic datasets. In fact, converging lines of evidence strongly suggest the presence of scale invariance in both topological and metric aspects of trees (4⇓⇓⇓⇓–9). These works provide a quantitative and data-rich analysis that is in the same spirit as hypotheses made earlier on the basis of fossil extinction and taxonomic evidence (10, 11) and discussed theoretically with reference to critical models of evolution (12, 13) (see, e.g., refs. 14⇓–16 for detailed reviews of the literature prior to the genomics era). The new element of the recent analyses (4⇓⇓⇓⇓–9) is that they are based on molecular phylogeny, and so contain far more information about the evolutionary process than can be obtained from patterns of extinction.

Nodes of a phylogenetic tree stand for extant organisms (outer leaves) and their hypothesized ancestors (inner leaves). They can be analyzed to quantify the topology of tree branching (4⇓⇓–7). There are several metrics developed in the literature to characterize the topology of phylogenetic trees (4⇓⇓⇓–8). Here, we shall focus on the metrics defined in ref. 5, in which the basic idea is to count the number of subtaxa that diversify from a given node i. We call the first quantity *A*, then it can be shown from the definition of C and A that *B*, would have each node branching into two: one node that simply persists to the edge of the tree without subsequent branching, and a second node that branches just like the parent, into a nonbranching and a branching node. In this case,

How do real data scale? Remarkably, it is found that, over three orders of magnitude of A, there is a power-law scaling of

There have been many theoretical attempts to model the evolution of phylogenetic trees. See refs. 21⇓⇓–24 for comprehensive reviews. The equal-rates-Markov (ERM) model, first developed by Yule in 1924 (25) and later expanded in the literature (26⇓–28), usually serves as a null hypothesis for the evolutionary process of the tree. The ERM assumes that all extant species on the tree have the same speciation rate. The resultant tree is less unbalanced than the observed ones, and, statistically, it should behave on large scales as a balanced tree with the asymptotic scaling

Connecting the nodes, edges of a phylogenetic tree represent the evolution time from a parent node to its daughter(s). They can be analyzed by measuring the edge-length abundance distribution (EAD) (9). The EAD is calculated from the number of tips or leaf nodes k that descend from a given internal node i. For a given node i with k tips,

As reviewed in ref. 9, although the Yule process (25) and the Kingman coalescent (34) produced power-law EAD, they do not generate exponents that are measured with actual phylogenetic trees. Extended Kingman coalescent with time-varying rate (9) and Λ-coalescent (35, 36) can produce the observed value with a tuning parameter. Yet, the specific biological reason of the parameter choice is not clear.

One way to obtain anomalous power-law scaling for tree topology is to directly include power-law aging behavior into the rules for the generation of trees (4, 37⇓–39). For example, by requiring that branching probabilities are a particular power-law function of the branch age (38), it is possible to obtain both logarithmic and power-law scalings for

The structure of phylogenetic trees can be interpreted as arising from the interplay between evolution and ecosystem dynamics. We view the nodes as revealing information primarily about ecological processes that result in the fixation of beneficial mutations, whereas the edges reveal information primarily about evolutionary processes, because their lengths reflect the numbers of DNA mutations. The presence of a power-law aging or long-term memory implies a breakdown of the separation of scales implicit in the identification of nodes with ecological interactions and edges with evolutionary dynamics. This standard identification would be valid if the ratio R of ecological timescales to evolutionary timescales could be assumed to be zero. However, even though

Niche construction (43⇓⇓⇓⇓⇓⇓⇓–51) is a term that describes the fact that organisms modify the environment and, thus, create new ecological niches; in turn, these niches affect the evolutionary trajectory of other organisms that share the environment (44). The resulting dynamics is a coevolution of the coupled dynamical variables for the organisms (52⇓⇓–55) or their genomes (45), as well as the environment itself. This coupled dynamics contains two-way feedbacks between the organisms and the environment, which are local in time. However, phylogenetic trees follow only the dynamics of the organisms themselves. The effective theory for the organismal degrees of freedom can be obtained conceptually by integrating out the environmental variables (e.g., using functional integration), and the resulting description would then contain interactions that are nonlocal in time, leading to an effective long-term memory in the branching process.

We will see that niche construction indeed introduces long-lived memory into the evolutionary process, and even a very simple caricature of niche construction over evolutionary time can capture the power-law scaling of

The niche of a species generally refers to its role or function in an ecosystem and can be thought of as the “variables by which species in a given community are adaptively related” and which control the species population response to each other and their environment (56). The habitat occupied by a given species is distinct from the niche, in this definition, and the two together comprise the ecotope (56). The ecotope involves the environmental factors that the species relies on, including the geographic configuration, the climate, etc., and the interactions with other species in the same ecosystem, represented by the species’ position in the food web and their dynamical history. The phenomenon that organisms modify the environment and thus create new niches is termed "niche construction" (44). In contrast to natural-selection theory, which treats the environment as a static stage on which population dynamics and genetics occur, niche construction theory emphasizes the modification of the environment by the organisms as an explicit process and a key factor in evolution (43⇓⇓⇓⇓–48), although there remain controversies (49⇓–51). Specifically, organisms can shape the environment they live in, change the selection pressure, and, as a result, reroute their own evolutionary path. A similar feedback of organisms on their environment is often referred to in ecology as ecosystem engineering (57⇓⇓–60), which is on a shorter timescale than that of niche construction. Even though it can be argued that niche construction really means ecotype construction, here, we will use the term niche, as is frequently done.

Ecological models, such as the MacArthur–Levins model (61), typically treat a niche as an abstract one-dimensional space which the organisms inhabit (although multidimensional niche models also exist; see, for example, refs. 56 and 62). We will regard a niche as the total available growth space or evolutionary degrees of freedom of the organism. We will call n the available niche value in the analysis below. When we say an organism has a large niche value, we mean that it has a large number of possible ways to adapt to the environment and to eventually survive or to reach genome-type fixation. On the other hand, an organism with a small niche value is very discriminating with regard to its environment, and this will affect its ability to be resilient within a wide range of environmental fluctuations. Our usage of niche means that, following a speciation event, the daughter species have a similar niche value, but with some fluctuation. The niche value can change either by organismal influence on the environment or through the evolution of mutations that enable key innovations to arise.

Our work is, of course, not the first to attempt to model niche construction, but the ingredient here is the invoking of critical scaling theory to explain the observed topological scaling laws of the large-scale evolutionary dynamics. Of earlier work in this area, we specifically draw attention to applied population-dynamics models (52⇓⇓–55) and population-genetics models (45) that study the effect of niche construction or ecosystem engineering on organism populations and evolution.

This article is organized as follows. First, we present a minimal model of the large-scale effects of niche construction, which we call the Niche Inheritance Model. In this model, the descendant species inherit the parent’s niche with fluctuations due to niche construction and evolution. The model is a caricature of the most significant ecological interactions that influence a phylogenetic tree in our assessment, but, based on a huge body of work on scaling laws, we anticipate that such a minimal model will yield nontrivial predictions that can be in agreement with experimental data (40). Next, we show that an apparent power-law regime develops when strong niche construction (destruction) leads nodes to be deactivated due to a lack of niche. The scaling laws are revealed through data-collapse scaling. Finally, we end with some discussion about the use of minimal models in evolutionary ecology at large timescales.

## Results

### Niche Inheritance Model.

We assign each species node three attributes: the amount of available niche n, the speciation rate r, and the extinction probability e. The tree-generation algorithm is as follows. Let the parent node be represented by its parameters **3** and **4**, which will be discussed in the next paragraphs. Then, we test whether the node goes extinct or remains to bifurcate later. The test is done by drawing a uniformly distributed random number in

The speciation rate r is treated as an increasing function of niche to reflect the fact that the more available niche space there is, the more likely it is for a speciation event to be successful. Specifically, with the interpretation of available niche n as the total available growth space or evolutionary degrees of freedom, when an organism has a large niche value, it has a large number of possible ways to adapt to the environment and to eventually survive. In this way, the speciation rate r naturally positively correlates with the niche n. To capture the first-order feature and to develop a minimal description, we set the relation between r and n to be linear, as shown below.

The extinction probability e implicitly incorporates all ecological interactions that can lead to extinction after a species emerges. Despite the fact that there are multiple factors driving species extinction, there is no well-accepted way to quantify the strength of each factor or to braid them together. In this model, we assume that the extinction probability e increases with the speciation rate r, based on the reasoning that a large speciation rate results in a big group of competitors with similar niche and, thus, reduces the survival probability of individual species. In our effort to build a minimal model, we assume a simple functional form of

In a numerical simulation, we start with a root node and evolve the tree based on the above rules until it reaches a certain size. Then, we compute A and C using the following definitions. For an arbitrary node i on the tree, let

### Existence of the Absorbing Boundary.

In the above framework, there exists a boundary case of

Imagine the left node starts by chance with a larger n and, thus, a higher r than the right one. If

### Effect of Niche Construction Strength.

While the absorbing boundary induces imbalance in the tree, the frequency of nodes hitting the boundary also matters. Based on Eq. **2**, we see that

We demonstrate in Fig. 2

In Fig. 2, we average the C values corresponding to the same A and plot the resulting quantity

For zero niche construction, *A*. Notice that the scale is linear-logarithmic, and the *B*. Instead of

Fig. 2*C* shows the averaged

So far, we have shown that niche construction together with the absorbing boundary induces a power-law scaling regime in the

### Singularity Induced by the Absorbing Boundary.

We attempted to describe the scaling behavior using a simplified mean field that ignores the fluctuations of the tree topology due to the random-species birth and death processes (see the nonextinctive mean-field model in *SI Appendix*). Although the mean-field calculation explains the exact scalings for the two extrema of balanced and maximally unbalanced trees, it fails to recapitulate the scaling laws that we observe in the numerical simulations. However, it is well known from the theory of critical phenomena that nontrivial power-law scaling arises from singularities in limit processes (40) that cannot be captured by mean-field theory. Thus, we now focus on the singularity induced by the absorbing boundary.

The imbalance induced by a large niche construction effect is crucially related to the condition

In the above argument, we have implied that nodes are able to reach the

In Fig. 3, we show the dependence of *B*. This is demonstrated as the segment of straight line under the double-logarithmic scale in Fig. 3*A*. As *B*.

### Critical Scaling at the Absorbing Boundary.

Although the behavior of

With this conjecture, we observe that

The data collapse indicates a critical behavior of *B*. From the function

In the above two sections, we have considered a soft absorbing boundary of finite and small

### Power Law in the EAD.

In this section, we show that the niche construction model also reproduces another characteristic of phylogenetic trees: the power law in the EAD, first discovered by O’Dwyer et al. (9). For convenience, we briefly revisit the definition of this distribution. The distribution is denoted by a function

In Fig. 5, EADs of the niche construction model are presented. Each EAD has been scaled differently to reduce overlap in the plot. In the double-logarithmic scale, *k* < 1,000, indicating a power-law behavior. The scaling range is comparable with the observation of real phylogenetic trees (9). The data points on the right side are more scattered because large clades are sampled insufficiently, as in the case of

### The Necessity of Eco-Evolutionary Feedback in a Minimal Model.

In this section, we will connect the scaling behavior of **3**, this coupling represents perhaps the simplest nontrivial feedback between the ecological variable, niche values, and the evolutionary variable, edge lengths. To explore what are the essential ingredients of a minimal model for phylogenetic tree structure, we will demonstrate the effects of modifying our model. We will show that, without this feedback, it is not possible to recapitulate scale invariance and the anomalous scaling laws, and so this is an essential part of a minimal model for phylogenetic tree structure.

To begin with, we recall that the edge length of a particular node is determined by its speciation rate r. More specifically, we required the edge length, or time till speciation, to follow an exponential distribution with parameter r. We could modify this element of the model by requiring all edge lengths to be equal to a constant number instead. We will set *A*, none of the curves have a noticeable range of power-law scaling. It’s worth noting that the finite and constant edge lengths effectively satisfy the infinite time assumption in the mean-field calculation (*SI Appendix*). As a result, when extinction events are not frequent, the analytical result in *SI Appendix* (*SI Appendix*, Eq. **S11**) should apply to this modified model. Indeed, for small fluctuation strengths, the analytical curves agree well with the simulation in Fig. 6*A*. This plot demonstrates that node deactivation alone is not sufficient to produce the power-law behavior in

Eliminating variability of edge lengths might be too drastic a change. Therefore, in the next modification, we allow the edge length to follow an exponential distribution with a constant rate. In contrast to the original model, this modified version does not retain the coupling between the niche value and the speciation rate (Eq. **3**). Fig. 6*B* shows that, even with variable edge lengths, *C*, data for

The above two modifications illustrate that edge lengths need to be not only varying, but also coupled to niche values in order to generate realistic topological structures. We could understand the effect of edge lengths by considering a node with a very small niche value and, hence, a very small speciation rate. Such a node most likely has a very long waiting time before speciation, but because the simulation time in real life is finite, the node cannot speciate. On the other hand, nodes with large niche values will speciate more often and have more child nodes, thus causing the imbalance in the tree. This distinction between large and small niche values is present only because the model couples niche values with speciation rates. Without this coupling, an identical distribution of edge lengths is insufficient to induce enough imbalance in the tree, and *B*. In conclusion, our specific mechanism to assign edge lengths not only reproduces realistic statistics of edge lengths, but also is essential for a realistic topological structure.

All data associated with the manuscript are accessible publicly, via https://github.com/zhiru-liu/niche-inheritance-trees.

## Discussion

We have presented a model to explain the observed universal scaling of phylogenetic trees. We incorporate niche construction as an explicit evolutionary process in the tree growth. By analyzing the Niche Inheritance Model, we make two significant conclusions. First, a large niche construction effect, together with the absorbing boundary, leads to an apparent power-law regime in the tree topology. This is in the same range of A as observed in actual phylogenetic trees (5). The existence of the power-law

Our model has simple rules for the evolution of the tree. The significance is that there is a local in-time interplay between the speciation rate and niche availability and that this can generate a critical behavior in

There are several issues that require further investigation. First, we have predicted a scaling form for the cross-over point

Our results show that niche construction is more than a feedback between evolutionary and ecological processes arising when their timescales are not widely separated. Niche construction not only leads to a perturbation in the evolutionary trajectories of all components of an ecosystem, but also creates an indelible footprint on the evolutionary process that cannot be eliminated, even for very long times. These memory effects manifest themselves through the anomalous scaling laws that characterize observed phylogenetic trees.

## Acknowledgments

We thank James O’Dwyer and Kevin Laland for their critical reading of the manuscript and helpful suggestions. A portion of this paper was adapted from the dissertation of C.X. This work was supported by the NASA Astrobiology Institute under Cooperative Agreement NNA13AA91A issued through the Science Mission Directorate.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: nigel{at}illinois.edu.

Author contributions: C.X. and N.G. designed research; C.X., Z.L., and N.G. performed research; and C.X., Z.L., and N.G. wrote the paper.

Reviewers: M.W.F., Stanford University; J.K., University of Cologne; and K.S., University of Copenhagen.

Data deposition: All data associated with the manuscript are accessible publicly via GitHub (https://github.com/zhiru-liu/niche-inheritance-trees).

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1915088117/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- ↵
- ↵
- K. Vetsigian,
- C. Woese,
- N. Goldenfeld

- ↵
- E. Hernandez-Garcia,
- M. Tuğrul,
- E. Alejandro Herrada,
- V. M. Eguiluz,
- K. Klemm

- ↵
- ↵
- P. J. Maldonado

- ↵
- N. Goldenfeld

- ↵
- C. Colijn,
- G. Plazzotta

- ↵
- J. P. O’Dwyer,
- S. W. Kembel,
- T. J. Sharpton

- ↵
- ↵
- ↵
- ↵
- J. Chu,
- C. Adami

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- A. Masucci

- ↵
- ↵
- A. O. Mooers,
- S. B. Heard

- ↵
- ↵
- ↵
- ↵
- G. U. Yule

- ↵
- D. G. Kendall

- ↵
- ↵
- ↵
- ↵
- ↵
- D. Aldous,
- R. Pemantle

- D. Aldous

- ↵
- ↵
- ↵
- ↵
- ↵
- N. Berestycki

- ↵
- M. Stich,
- S. Manrubia

- ↵
- S. Keller-Schmidt,
- M. Tuğrul,
- V. M. Eguíluz,
- E. Hernández-García,
- K. Klemm

- ↵
- S. Keller-Schmidt,
- K. Klemm

- ↵
- N. D. Goldenfeld

- ↵
- G. I. Barenblatt

- ↵
- L. Y. Chen,
- N. Goldenfeld,
- Y. Oono

- ↵
- D. S. Bendall

- R. C. Lewontin

- ↵
- H. C. Plotkin

- F. J. Odling-Smee

- ↵
- K. N. Laland,
- F. J. Odling-Smee,
- M. W. Feldman

- ↵
- F. J. Odling-Smee,
- K. N. Laland,
- M. W. Feldman

- ↵
- K. Laland,
- B. Matthews,
- M. W. Feldman

- ↵
- ↵
- ↵
- M. Gupta,
- N. Prasad,
- S. Dey,
- A. Joshi,
- T. Vidya

- ↵
- M. W. Feldman,
- J. Odling-Smee,
- K. N. Laland

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- G. Barker,
- E. Desjardins,
- T. Pearce

- G. Barker,
- J. Odling-Smee

- ↵
- ↵
- ↵
- T. Biancalani,
- L. DeVille,
- N. Goldenfeld

- ↵
- S. Bornholdt,
- K. Sneppen,
- H. Westphal

- ↵
- N. Aktekin

- ↵
- S. Nee,
- A. O. Mooers,
- P. H. Harvey

- ↵