Modelling species interactions in macroevolution and macroecology Xu, Liang

University of Groningen

Modelling species interactions in macroevolution and macroecology Xu, Liang

DOI:

10.33612/diss.125954510

IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from it. Please check the document version below.

Document Version

Publisher's PDF, also known as Version of record

Publication date:

2020

Link to publication in University of Groningen/UMCG research database

Citation for published version (APA):

Xu, L. (2020). Modelling species interactions in macroevolution and macroecology. University of Groningen.

https://doi.org/10.33612/diss.125954510

Copyright

Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).

Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.

Download date: 24-06-2021

Modelling species interactions in

macroevolution and macroecology

ii

Printed by: Ridderprint www.ridderprint.nl

Front & Back: Beautiful cover art is designed by Li Huang.

Copyright 2 2 by Liang u

ISBN

- 4- 4-2 -1 printed version

- 4- 4-2 4- electronic version

iii

Modelling species interactions in macroevolution and macroecology

PhD thesis

to obtain the degree of PhD at the University of Groningen

on the authority of the Rector Magnificus Prof. C. ijmenga

and in accordance with the decision by the College of Deans.

This thesis will be defended in public on Friday June, 2 2 at 12:4 hours

by

Liang Xu

born on 1 April 1

in Chong ing, P.R.China

Supervisor Prof. R.S. Etienne Co-supervisor

Prof. G.S. van Doorn Assessment Committee

Prof. L. Harmon

Prof. F. Bokma

Prof. I.R. Pen

To my family

Contents

1. Introduction . . . . 1

2. Detecting local diversity-dependence in diversification 11 3. Inferring the effect of species interactions on trait evolution . . . 29

4. The spatially explicit phylogenetic Janzen-Connell effect predicts realistic macroecological and macroevolutionary patterns . . . 63

5. Synthesis . . . 83

6. Bibliography . . . 95

7. Summary . . . 115

Samenvatting . . . 119

8. Acknowledgements . . . 123

Curriculum Vitae . . . 129

vii

1

Introduction

1

1

2 Introduction

A biological interaction is the effect of the activities of one organism on another organism. It has been widely recognized as the fundamental basis of a vari- ety of ecological and evolutionary phenomena. On the species scale, a biological interaction among species is known as interspecific interaction, and includes mutu- alism, commensalism, parasitism, amensalism, neutralism and competition, which are classified by distinct degrees of benefit or harm species cause to others. The interspecific interactions often have a long-term effect, affecting, for instance, bio- diversity and morphological trait evolution. Thus, studying the underlying mecha- nisms how species interact helps us to understand species assemblages. Biological interactions also happen on the individual scale within species, which is normally known as intraspecific interaction. This type of interaction is mainly due to mem- bers of the same species competing for limited resources such as food, water, living space and mates etc. The conse uence of this competition is the slowdown of population growth that ultimately leads to an e uilibrium population size if the environment does not change, although competition may also result in population cycles or chaos. Biologists use carrying capacity to express the number of in- dividuals or species that a stable environment can maintain on the corresponding scales. The interactions of both scales can operate and interplay on the same time scale, resulting in a complex process influencing biodiversity.

Ecologists and evolutionary biologists have made a great effort trying to re- veal the underlying mechanisms of biological interaction that explains biodiversity.

However, while significant progress in understanding biological processes has been achieved, our knowledge is still limited by the large spatial and temporal scales involved in those processes. For example, ecological process occurs anywhere and anytime such that it is usually difficult to fully track. Macro-evolutionary changes that are significant enough for detection often take a large amount of time, for instance, in unit of million years for speciation and extinction. Thus, a practical method that is able to tackle uestions on large spatial and temporal scales is imperative. Theoretical approaches along with the development of computer sci- ence make this possible. Theoreticians exploit mathematical modelling to study the fundamental interaction patterns hidden in empirical data, which is collected from the extant species and from the fossil record that stores information of the extinct species in the past.

Modelling macroevolution

In 1 4 , based on Corbet’s work on the distribution of butterflies Corbet and illiams’ data of numbers of moths of different species caught in a light-trap

Fisher et al. Part 2 , Fisher Fisher et al. developed a mathematical the-

ory to predict the relative numbers of animals of different species obtained when

sampling at random from a ecological community. Later on, illiams illiams

21 , 21 showed that the logarithmic series used to address the problem of the

fre uency of occurrence of species in a random sample from a mixed population

can also be applied to a great variety of other biological distributions. These suc-

cesses in phenomenological descriptions by mathematics drew people’s interest to

Introduction

1

3

developing theoretical models to explain the biological distributions. Soon, in 1 4 , endall endall tried a number of discontinuous Markov processes that lead to distributions of negative binomial and logarithmic series forms and hoped future researches may find some biological meaning in them. One of the stochastic pro- cesses is a birth-death immigration process, which became the basis of theoretical models of species diversification.

In 1 4, inspired by endall’s work, Nee et al . Nee et al. 1 described a probabilistic null model appropriate for molecular phylogenies, which records lin- eages that have given rise to at least one contemporary descendant. For this model, which is also called the reconstructed process , Nee et al . derived a geometric distribution for the number of lineages existing at any particular time in the recon- structed process, and the distribution of waiting times between birth events, and generalized the results to time-dependent birth and death rates. A likelihood func- tion for a reconstructed phylogeny is derived to estimate the speciation rate and extinction rate. However, the diversification rate was not yet related to other biolog- ical processes, such as trait evolution. Almost at the same time, Pagel Pagel 144 noticed the importance of trait evolution and presented a new statistical method for analysing the relationship between two discrete characters that are measured across a group of hierarchically evolved species or populations. Although likeli- hood calculations for models involving speciation and extinction rates Moore et al.

12 , Nee et al. 1 and rates of character state change have been described separately Pagel 144 , they have not yet been fully integrated. But biologists’

attention started to be drawn towards this direction. In the next decade, based on Nee et al .’s work some research on the state-dependent speciation rates was done by Pagel Pagel 14 , Paradis Paradis 14 and Ree Ree 1 . Pagel’s model allows different rates of speciation to different character states, but it assumes no extinction and states change only at speciation events. Paradis Paradis 14 and Ree Ree 1 presented likelihood-based methods that use reconstructed ances- tral states to compare speciation rates between states but also excluded extinction.

Shortly after, the gap was filled. In 2 , Maddison et al. Maddison et al.

11 developed a binary-state speciation and extinction model BiSSE that allows speciation and extinction rates to depend on the character state of a lineage at each point in time and allows the character state to change over time. A method to calculate the likelihood with two states was presented. They showed how this calculation leads to new methods for parameter estimation and hypothesis testing of a binary character’s effect on diversification along a full phylogenetic tree that contains extinct species. The BiSSE model soon attracted biologists’ interest and inspired many other trait-dependent diversification models, becoming the basis of these models.

However, the BiSSE model assumes that the phylogenetic tree is complete and

fully resolved, and that all character state information is known. These assump-

tions restrict its applicability as few published phylogenies are both complete to the

species level and large enough to detect differential diversification. So, FitzJohn

Fitzjohn et al. improved the model to make it applicable to incomplete phylo-

genies. He developed likelihood calculations that compensate for incomplete phy-

1

4 Introduction

logenetic knowledge in two cases: 1 incomplete random sampling of all extant species within a group and 2 unresolved clades instead of tips in phylogenies. In 2 1 , FitzJohn FItzJohn 4 described a new comparative phylogenetic method, the uantitative state speciation and extinction QuaSSE model, which is based on the BiSSE model, for inferring the effect of uantitative traits on speciation and extinction rates. Also, the likelihood of a phylogenetic tree and the distribution of character states among species are derived. FitzJohn’s work promoted the devel- opment of macroevolution models in trait-dependent diversification rates, opening a door for biologists to further focus on a variety of trait types.

Among the traits of interest studied by the invoked trait-dependent diversifica- tion research trend, geographic traits that crucially influence diversification rates stand out. Several researches emerged to study geographic range as traits that are likely to influence speciation and extinction Cardillo et al. 2 , Jablonski and Roy 2 , Mckinney 122 , Pigot et al. 1 , Ribera et al. 1 . Mc inney Mck- inney 122 found that species with larger ranges are commonly considered less prone to extinction than those with smaller ranges. Range contraction events will therefore increase the rate of extinction, and they may cause anagenetic speciation stochastic extirpation of local populations or cladogenetic speciation speciation fragmenting an ancestral range into smaller descendant ranges . These effects are counterbalanced by range expansion, achieved through dispersal events that estab- lish new populations. Rosenzweig, Chown , Gaston and Chown Chown 2 , Gaston and Chown , Rosenzweig 1 showed that changes in range size may have a variety of effects on speciation rates.

There are many uestions related to geography and diversification. Are regional differences in species richness and endemism driven by spatial asymmetries in spe- ciation and extinction rates, or by asymmetries in direction of lineage dispersal Chown and Gaston 2 , Mora and Chittaro 12 , Stebbins 1 4 ? How do range size and range evolution affect speciation rate Gaston and Chown , Jablonski and Roy 2 , Pigot et al. 1 , agner and Erwin 21 ? hat are the most common modes by which speciation divides ancestral ranges in terms of range size or degree of sympatry Anderson , Barraclough and ogler 11 , Gaston , Phillimore et al. 1 , Pigot et al. 1 ? How do speciation rates compare on the mainland, within an island, and directly after dispersal to the island Gillespie and Roderick , Ricklefs and Bermingham 1 ? Are specialists more likely to arise from generalist ancestors or directly from other specialists Barnett and Simp- son 1 , Nosil et al. 141 ? The need for phylogenetic models to infer these interactions has been raised previously in the context of biogeographic frameworks that assume constant speciation and extinction rates and hence allow consideration of geographic characters only on a static tree Lamm and Redelings 1 , Ree and Smith 1 , Ree and Sanmart n 1 1 .

In 2 11, Goldberg et al. Goldberg et al. developed a model for geo- graphic traits that encapsulates the basic concepts introduced above. The model is also based on the mathematical framework of the BiSSE model Maddison et al.

11 , modified to include biogeographic parameters and to allow state change at

speciation and through local extinction. The geographic state speciation and ex-

Introduction

1

5

tinction GeoSSE model specifies the states as the geographic traits and uses the new likelihood-based approach to estimate region-dependent rates of speciation, extinction, and range evolution from a phylogeny. Also, Goldberg et al. applied it to the evolution of habitat occupancy in Californian plant communities, where they found higher rates of speciation in chaparral than in forests and evidence for expanding habitat tolerances.

However, in this series of speciation and extinction models, diversity-dependent diversification rates are rarely considered. Empirically, palaeontological evidence for diversity-dependence has been found from the relatively constant diversity within numerous higher taxa over millions of years Alroy , Ezard et al. , Gould et al.

, Rabosky and Sorhannus 1 4 , Raup et al. 1 , Stanley 1 . If diversity- dependence and non-zero extinction are true, models of diversification should in- clude both processes. Several existing models incorporate diversity-dependence but without extinction. This is in part because simulations reveal that extinc- tion erases the signature of a reduction in the speciation rate through time Liow et al. 1 , Quental and Marshall 1 , Rabosky and Lovette 1 2 , leading to the suggestion that extinction rates have probably been low in clades that show slowdowns in lineage accumulation Rabosky and Lovette 1 2, 1 . Another explanation is just for convenience to compute likelihood Bokma 1 , Rabosky and Lovette 1 2 . In 2 12, Etienne et al . Etienne et al. 1 used a hidden Markov model HMM considering diversity-dependence in speciation rate and con- stant non-zero extinction rate to numerically compute the likelihood of a phylogeny.

They considered incomplete sampling of species and presence of other species that have gone extinct but affected diversification rates in the past. The model was applied to the phylogenies of clades with an existing fossil record Cetacea and Cenozoic of macroperforate planktonic foraminifera and clades without fossil ev- idence Dendroica , Plethodon and Heliconius . They concluded that the method performs uite well in estimating a non-zero extinction rate from molecular phylo- genies with and without fossil record in some case studies. They suggested that the diversity-dependent diversification model with extinction should be preferred over the constant-rate birth-death model or pure birth model as a more biologically realistic model for macroevolution.

Nevertheless, the diversity-dependence diversification model is designed as a simple model of diversity-dependence, serving more like a null model. It assumes a global scenario, thus, ignores spatially distributed clades that do not interact but may have different influence on clade-specific speciation. In addition, the spatial structure of the community is widely recognized as an essential driver of biodiversity Chown 2 , Gaston and Chown , Rosenzweig 1 . Thus, this simplification may reduce the power of the likelihood formulation for parameter inference.

Modelling trait evolution

Although evolutionary biologists have long been fascinated by paleontology, sys-

tematics, morphology, and genetics, the interactive dynamics of ecological and evo-

lutionary processes received little attention Schoener 1 . Only recently people

1

6 Introduction

started to realize the substantial effect of both directions between ecology and evolution. In fact, the influence of ecological processes on evolution has been doc- umented by many empirical phenomenons. One of the best known instances is the evolution of Gal pagos ground finches, which are also called Darwin’s finches to memorize Darwin’s study by Lack in his book Darwin’s Finches Lack 1 2 in 1 4 . Darwin’s Finches comprise a group of about 1 species. Their beak sizes are evidenced to associate with the size of seeds that are available in environment, showing apparent evolutionary evidence responding to ecological changes. In the species Geospiza fortis , people found that the beak size of the birds is larger when large seeds are more available Grant and Grant . Conversely, when small seeds become more available for some years, small beaks are pervasive in the pop- ulation. The study of the other direction, i.e. of evolution affecting ecology, is tricky because such processes were believed to be too slow to observe. However, recent studies have looked at ecological and evolutionary processes at similar time scales Hairston et al. 2 , Reznick et al. 1 , Reznick and Ghalambor 1 4 , showing that evolution can be fast, thus, plays a role in ecological processes. One empirical example is the evolution of Caribbean lizards Strauss et al. 1 . Before preda- tors invaded the community, the small lizards Anolis sagrei were happy living on the ground and tree trunks. The invasion of the predator Leiocephalus carinatus - a large lizard with curly tail that feeds on Anolis sagrei - results in reduction of the population size of Anolis sagrei . Few small lizards with shorter limbs that better fit the higher branches of the trees where the predator cannot reach, survive. After a few generations, this selection rapidly caused the shift of the limb size distribution towards shorter limbs and increased the abundance of Anolis sagrei .

The first theoretical attempts to model trait evolution occur in the 1 s. In 1 2, right right 21 developed a model to illustrate a basic principle of the evolution of gene fre uencies. The model describes that in a static environment with random mating, selection causes the gene fre uency at a locus to change in order to maximize the mean fitness of individuals in the population until e uilibrium where the fitness reaches a maximum. The fitness here means that how individuals are adapted to the environment thus have a higher chance to survive and give birth to offspring.

In 1 , Lande Lande 1 4 clarified the concepts and constructed a mathe- matical model to describe an adaptive topography for the average phenotype in a population:

Δ ̄𝑧(𝑡) = ℎ 𝜎 𝜕𝑙𝑛 ̄ 𝑊

𝜕 ̄𝑧(𝑡) 1.1

where ̄𝑧(𝑡) denotes the mean phenotypic trait value of the population at time 𝑡.

The change of the mean trait Δ ̄𝑧(𝑡) is affected by the change of the mean fitness

with respect to the mean trait. The formula uantitatively measures how the mean

trait of a population changes to better fit the environment. In addition, it consid-

ers population genetics. The heritability ℎ is determined by the genetic system,

the breeding structure of the population and the environment. It measures the

proportion of the variation in a given trait within a population that is explained

by inheritance from the parent generation instead of the environment or random

Introduction

1

7

chance. The variance of the trait is denoted by 𝜎 , which is normally assumed to be constant and independent from the mean trait Falconer . However, Simpson Barnett and Simpson 1 and an alen an alen 2 in their case studies also argued that this may not be the case. A changing variance in the cause of evolution is more realistic but also leads to a much more complex system.

Nevertheless, Lande’s formula provides a uantitative tool to model trait evolution.

In fact, at the meantime, there are two alternative models dominating the trait evolution studies even till now, i.e. Brownian motion process and Ornstein Uhlenbeck O-U process. hile the Brownian motion process mimics random walks of traits in evolution Felsenstein and the O-U process describes a selective force pulling traits towards an optimum, they assume independent evolution for each species and thus do not account for species interactions Pennell and Harmon 14 . Nevertheless, because of analytical tractability, they are majorly favored to study trait evolution. Apart from these two models, Lande’s model is well biolog- ically interpretable. It assumes a fitness landscape such that species struggle to climb up the hills of the fitness landscape via evolving their traits to better fit the environment.

However, early studies of trait evolution ignored phylogenetic information. Only in 1 Felsenstein Felsenstein 1 introduced a phylogenetic independence con- trasts method to test hypothesis of modes of trait evolution. Felsenstein pointed out that given the fact of species being part of a hierarchically structured phylogeny, if the traits are drawn independently from the same distribution they cannot be used for statistical tests among models. This can be circumvented by considering the phylogenetic relatedness of the species in traits. The novel idea of integrating phylogenetic information opened a new era and spawned a family of uantitative methods to apply models of population and uantitative genetics, paleobiology, and ecology to data including a phylogenetic tree, the family of which is now well known as phylogenetic comparative methods PCMs . In fact, the models mentioned in Section Modelling macroevolution all fall in this family, normally regarded as model-based approaches to investigate rates of diversification and the influence of characters on diversification.

Nowadays, the progress and directions of PCMs start to move towards a more integrative comparative biology. Phylogenetic community ecology is one of the ar- eas in which PCMs have played an increasingly important role Pennell and Harmon 14 . The study of phylogenetic community ecology is to exploit the phyloge- netic relatedness among species in traits to explore how the community structure and patterns are formed ebb et al. 211 . However, species differences in rela- tive abundance, competition and community functioning are rarely incorporated in phylogenetic comparative methods.

Fortunately, the mismatch between community ecology and phylogenetic com- parative methods has started to attract theoreticians’ attention. In 2 1 , Nuismer

& Harmon 14 presented a model that allows the traits of species to evolve along

a phylogeny in an interactive manner to infer the rate of trait evolution and inter-

action. Later on, Drury et al . Drury et al. 4 , 44 derived a likelihood function

to the spatial extension of the Nuismer & Harmon’s model for estimating parame-

1

8 Introduction

ters of interest and tested the statistical power among several relevant models of detecting the impact of species interactions on trait evolution. Their work filled the gap of trait evolution among interacting species. However, the influence of distinct species abundances is still not accounted for.

Modelling individual interactions in community ecology

Ecologists have spent decades on studying species abundance distributions while phylogenetic comparative methods usually assume an e ual population size for all species under consideration Pennell and Harmon 14 . To study the influence of differences in species abundances on community ecology, individual-based simula- tion models that make use of individual properties are thus a powerful approach.

Through simulating ecological processes such as birth, death, migration and specia- tion on an individual basis, one can track the population size of species and therefore assess its impact on interspecific and intraspecific interactions, and further study how community patterns are formed. One of the studies of individual-based com- munity ecology traces back to Hubbell’s neutral theory in 2 1 Hubbell , which borrowed the idea of the neutralist theory of evolutionary genetics developed by imura, Crow, Ewens Ewens , imura , imura and Crow and built upon MacArthur and ilson’s theory of island biogeography MacArthur and ilson 114 . Neutrality means that all individuals in the community share identical rates of birth, death and migration regardless of their species identities. The neutral the- ory has achieved a great success in explaining species abundance distributions in some empirical cases. hilst it opened a door to modelling the change of species abundances along with producing a phylogeny. However, there are mismatches of neutral theory’s predictions and observed empirical patterns. For example, it cannot produce highly diverse communities while at the same time having very dominant species, and it produces phylogenies that show a strong pull-of-the-present. hile neutralists have sought the explanation of these mismatch in simplifications of the model with respect to speciation mode and the way of modelling space, many ecol- ogists argue that this happens because neutral theory ignores differences between species and among individuals.

Speaking of species differences, the negative density effect proposed by Janzen

Janzen and Connell Connell 4 independently is a well-known explanation

of the maintenance of high biodiversity and species heterogeneity, which is also

known as the Janzen-Connell hypothesis. The theory states that the most abundant

species are at a disadvantage because they are more prone to be explored and

attacked by natural enemies. Thus, rare species have more chance to grow in the

vacant sites, promoting species richness. One might think that combining of the

Janzen-Connell hypothesis and the neutral theory may resolve the insufficient power

of the neutral theory to explain high biodiversity in communities. In 2 1 , Levi et

al . Levi et al. 1 presented a model with similar initiatives to study whether

Janzen-Connell effect can maintain a practical species diversity and to answer how

strong Janzen-Connell effects must be to maintain observed levels of tree species

richness.

Introduction

1

9

However, pathogens and herbivores are rarely monophagous and usually at- tack several closely related hosts Agrawal and Fishbein 2 , Gilbert and ebb 4 , Novotny et al. 142 . Thus, closely related species should have a greater negative effect on the survival rate of each other than distantly related species. This mechanism has been demonstrated by several observations in tropical rain forest Gilbert and ebb 4 , Novotny et al. 142 , ebb et al. 211 . For example, from the perspective of the host, ebb et al . ebb et al. 212 used community-wide seedling mortality data to analyze the relationship between seedling survival and phylogenetic neighborhood effects in a tropical rain forest. They found that seedling survival was enhanced in a more heterogeneous neighborhood where phylodiver- sity is high. They further suggested that interactions with pathogens may explain positive effects on seedling survival. From the perspective of the natural enemy, by experimental inoculation of plant leaves with fungal pathogens in a tropical rain for- est, Gilbert & ebb Gilbert and ebb 4 found that the chance that a pathogen is able to infect two different plant species increases with increasing phylogenetic relatedness.

Furthermore, the spatial structure of the community and dispersal ability of in- dividuals play an important role in biological interactions. Not all individuals of a species can invade suitable habitat patches within a community. Many species are strongly dispersal-limited Clark and Clark 2 , Seidler and Plotkin 1 . The geographical distribution of a species depends on its dispersal ability. Dispersal limitation determines patterns of diversity and can be used to explain species co- existence in neutral theory Etienne et al. , Macdougall and Turkington 11 . However, the conse uences of dispersal limitation for species diversity depend on the processes that determine the diversity and abundance of species within sites and exhibit a complex effect on species diversity Cadotte 22 , Schuler et al. 1 . Tilman applied experimental introductions to relax dispersal limitation and demon- strated that species diversity increases with the introduction of species in a grass- land community Tilman 2 4 . It might be because that the increase of dispersal helps species invade empty sites. In contrast, Mou uet and Loreau Mou uet and Loreau 1 found that species diversity declines as dispersal increases because of increased homogenization of the metacommunity. That is to say, the influence of dispersal on diversity patterns may differ for different spatial scales. In addi- tion, the density of a species is not constant in areas with different distances to an empty site. Thus, the J-C effect should show spatial pattern Hubbell , Levi et al. 1 .

Both the negative density effect Janzen-Connell hypothesis and the phyloge-

netic relatedness effect contribute to biodiversity where dispersal ability and spatial

pattern of the effect play a role. hile empirical evidence of the combining effects

has been found recently Liu et al. 11 , theoretical analysis is still lacking. Fur-

thermore, to what extent the phylogenetic J-C effect could make realistic predictions

on macroevolutionary and macroecological patterns is still unclear.

1

10 Introduction

Thesis outline

Closely bonded to the three sections of the introduction aforementioned, the fol- lowing chapters study the mechanism of biological interactions from the perspective of macroevolutionary diversification on species level, the perspective of the inter- and intra-specific eco-evolutionary interactions on trait evolution and the perspec- tive of individual’s birth, death and migration processes. All results and conclusions are summarized in Synthesis.

The slowdown of species accumulation is a widely observed phenomenon in reality. The hypothesis of diversity-dependence diversification is one of the expla- nations that achieve a great success in interpreting such pattern. However, current diversity-dependence diversification models ignore the geographic spatial structure, thus, can only assess the diversity-dependence signal on a global scale. In Chapter 2, I develop a spatial extension to the global diversity-dependence diversification model to mimic macroevolutionary events such as speciation, extinction, migration and contraction on species scale. I generate evolutionary history stored in a phy- logenetic form and apply the global analytic likelihood formulation to estimate the macroevolutionary rates, i.e. speciation rate and extinction rate. Furthermore, I exploit bootstrapping techni ue to examine whether the global approach is capable to detect the diversity-dependence signal on a region of different spatial levels.

As species abundance plays a role in species interaction, thus, further helps form trait evolution, which is usually ignored, in Chapter 3, I introduce a trait evolution model that takes population dynamics into account during the cause of the trait evolution. The model allows traits to evolve along a phylogeny and un- der environmental attraction towards an optimum trait and competitive repulsion.

An approximate Bayesian computation evolutionary algorithm is applied to inves- tigate the strength of stabilizing selection and competition. I, later on, illustrate our method on the body size of the baleen whales and compare the results with a trait evolution model without population dynamics and a trait evolution model with metabolic dynamics. Our work induces such a family of trait evolution models that can accommodate any type of species interactions in the trait evolution process.

In Chapter 4, I follow the intuition of the mechanism of the phylogenetic

Janzen-Connell hypothesis in an explicitly spatial model and develop a mathematical

formulation to calculate the colonization probability of individuals that reflects the

essence of the phylogenetic Janzen-Connell mechanism in a spatial manner. I con-

duct substantial simulations under different parameter combinations that account

for different strengths and interaction distances of negative density effect, phylo-

genetic relatedness effect and individual dispersal ability. I am interested in the

resulting species-abundance distribution, species distribution, characters of phylo-

genetic trees and species richness. The comparison of the results among a variety

of parameter explorations reveal the underlying mechanism of how those effects

affect species assemblage.

2

Detecting local diversity-dependence in diversification

Liang Xu, Rampal S. Etienne

This chapter has been published in Evolution 72, 2 1 22 .

11

2

12 Inferring local diversity-dependence

Abstract

Whether there are ecological limits to species diversification is a hotly debated topic. Molecular phylogenies show slowdowns in lineage accumulation, sug- gesting that speciation rates decline with increasing diversity. A maximum likelihood method to detect diversity-dependent diversification from phyloge- netic branching times exists, but it assumes that diversity-dependence is a global phenomenon and therefore ignores that the underlying species inter- actions are mostly local, and not all species in the phylogeny co-occur locally.

Here, we explore whether this maximum likelihood method based on the non- spatial diversity-dependence model can detect local diversity-dependence, by applying it to phylogenies, simulated with a spatial stochastic model of local-diversity-dependent speciation, extinction and dispersal between two local communities. We find that type I errors (falsely detecting diversity- dependence) are low, and the power to detect diversity-dependence is high when dispersal rates are not too low. Interestingly, when dispersal is high the power to detect diversity-dependence is even higher than in the non- spatial model. Moreover, estimates of intrinsic speciation rate, extinction rate and ecological limit strongly depend on dispersal rate. We conclude that the non-spatial diversity-dependent approach can be used to detect diversity- dependence in clades of species that live in not too disconnected areas, but parameter estimates must be interpreted cautiously.

Key words: Macroevolution, diversity dependence, parametric bootstrap,

simulations, phylogeny.

Inferring local diversity-dependence

2

13

Introduction

Understanding the potential ecological limits to species diversification remains a hotly debated topic Harmon and Harrison , ozak and iens 1 , Rabosky and Hurlbert 1 1 . The rising availability of molecular data to create phylogenies has motivated the development of a variety of methods to interpret lineage diversi- fication and better understand its mechanisms. Such methods include the lineages- through-time plot LTT - a semi-logarithmic plot that tracks the number of species that have descendants at the present through time. LTT plots indicate that species accumulation slows through evolutionary time Moen and Morlon 12 . This de- creasing rate of diversification has often been interpreted as a sign of diversity- dependence Phillimore and Price 1 4 , Pybus and Harvey 1 , Rabosky and Lovette 1 2, 1 , eir 214 , resulting in the absence of a correlation between the crown age of phylogenies and current-day diversity. Nevertheless, other expla- nations also exist including time-dependent speciation and/or extinction rates, or the protracted nature of speciation Etienne and Rosindell , Moen and Morlon

12 .

To infer the presence of diversity-dependent diversification from molecular phy- logenies containing only extant taxa, the standard procedure is to compare the fit of a diversity-dependent DD model Sepkoski 1 , alentine 2 to a model with no disversity-dependence, which is commonly known as the constant-rates CR birth-death model Raup et al. 1 . Diversity-dependent models assume that evolutionary radiations are facilitated by ecological opportunity Schluter 1 , and that speciation is more likely to happen when diversity is low. Importantly, while extinct species leave no descendants at present, they may have affected di- versification and hence also the phylogenetic patterns that are observed at present.

An algorithm to compute the likelihood of a model based on this idea from a species- level molecular phylogeny of present-day species which may be incomplete as long as the number of species not represented in the tree is specified was developed a few years ago Etienne et al. 1 . This likelihood not only allows for estimation of lineage diversification rates but can be used in likelihood-based tests to com- pare the model to other, diversity-independent, models. Standard tests based on the likelihood ratio and the corrected Akaike Information Criterion have recently been reported to be inade uate for the comparison of DD vs CR models because of violation of some of the assumptions leading to the chis uared distribution used in these tests, but a a bootstrap likelihood ratio test is available as an alternative Etienne et al. 2 . In summary, we currently have the tools to check whether and when diversity-dependence can be detected.

However, current models used to detect diversity-dependent diversification on

molecular phylogenies assume that the global species richness of a clade deter-

mines its rate of diversification, even if the species belonging to the clade do not

interact, for example because of disjunct spatial distributions. Hence, the uestion

arises how we can detect diversity-dependence in such occasions. The ideal solu-

tion would be a test with a spatial model that incorporates diversity-dependence. In

2 11, Goldberg et al. constructed a spatial model, the geographic state speciation

2

14 Inferring local diversity-dependence

and extinction model GeoSSE Goldberg et al. , which includes biogeographic states and allows state changes at speciation and through local extinction. How- ever, it is built on the mathematical framework of the binary state speciation and extinction model BiSSE Maddison et al. 11 and thus inherits the assumption from the BiSSE model that all the evolutionary parameters are constant or time- dependent Rabosky and Glor 1 , but not strictly diversity-dependent. Com- puting the likelihood for a spatial diversity-dependence model remains a challenge, however, because it needs to keep track of all species, even currently extinct ones, in all spatial locations. An alternative solution is to test whether the above men- tioned bootstrap likelihood ratio test based on the non-spatial diversity-dependence model can detect local diversity-dependence. In this paper we explore this option.

e extend the diversity-dependent diversification model to two locations con- nected by dispersal, where both speciation and dispersal are diversity-dependent.

In this spatial diversity-dependence model, we incorporate both allopatric speci- ation and sympatric speciation and assume constant extinction because diversity- dependent extinction seems at odds with empirical phylogenies Etienne et al. 1 . e simulate phylogenetic trees following this model using various values for its parameters, to subse uently estimate parameters using a non-spatial diversity- dependent model Etienne et al. 1 . e employ the bootstrap likelihood test to explore whether we can detect diversity-dependencewhen data are simulated under the spatial diversity-dependence model.

Materials & Methods Model

e introduce the simplest spatial diversity-dependent diversification model by assuming two regions, denoted 1 and 2. e call this model the spatial model. It is an extension of the diversity-dependent diversification model of Etienne et al. 2 12, which has no spatial structure, and hence will be called the non-spatial model. Our spatial model considers local macro-evolutionary processes sympatric speciation and local extinction as well as species interactions between locations through dispersal and allopatric speciation . Our aim is to explore whether the simpler non- spatial model can detect diversity-dependence from simulations under the more complicated spatial model, and whether parameters estimated using the non-spatial model relate in an informative way to the true parameters of the generating spatial model.

e assume that sympatric speciation rates are linear functions of the number of species present on the locations. e denote the number of species on locations 1 and 2 by 𝑛 and 𝑛 , respectively. Synpatric speciation rates 𝜆 (𝑛 ) and 𝜆 (𝑛 for both locations are defined as follows:

𝜆 (𝑛 ) = max(0, 𝜆

− (𝜆

− 𝜇) 𝑛

𝐾 ) 2.1

𝜆 (𝑛 ) = max(0, 𝜆

− (𝜆

− 𝜇) 𝑛

𝐾 ). 2.2

Inferring local diversity-dependence

2

15

Here, 𝜆

and 𝜆

are the intrinsic speciation rates of the two locations; these are the rates when diversity is 0. Furthermore, 𝐾 and 𝐾 can be interpreted as the carrying capacities for the two locations. e can rewrite these expressions as

𝜆 (𝑛 ) = max(0, 𝜆

(1 − 𝑛

𝐾 )) 2.

𝜆 (𝑛 ) = max(0, 𝜆

(1 − 𝑛

𝐾 )). 2.4

where we have defined

𝐾 = 𝜆

𝐾 /(𝜆

− 𝜇). 2.

The parameter 𝐾 can be interpreted as the maximum number of niches that the species in the clade can occupy Etienne et al. 1 , and hence it is an ecological limit to diversity.

Dispersal between the two regions is also assumed to be diversity-dependent :

𝑀

(𝑛 ) = max(0, 𝑀 (1 − 𝑛

𝐾 )) 2.

𝑀

(𝑛 ) = max(0, 𝑀 (1 − 𝑛

𝐾 )) 2.

where 𝑀 is the intrinsic dispersal rate when diversity is 0 in the receiving region, and the notation 𝑎 → 𝑏 stands for dispersal from location 𝑎 to location 𝑏. E . 2. and E . 2. show that dispersal rates are dependent on the diversity of the location species are dispersing to. Diversity-dependence is often based on a niche- filling argument: as diversity increases, it is increasingly harder for a new species to enter the community and find its own niche to establish in the community. Entering the community can occur either through speciation or through immigration. Hence, the rate of sympatric speciation and of dispersal both depend on the diversity in the location that the new species enters.

The conse uence of dispersal is that some species inhabit both regions at the same time; we will refer to these as widespread species . In contrast, we will call species residing on a single location endemic species . In our model we incorporate allopatric speciation, i.e. the split of a species that is present on both locations into two species, each present on one location. The allopatric speciation rate is assumed to be negatively related to the intrinsic dispersal rate

𝜆 = 𝜆

𝑀 2.

where 𝜆

is the allopatric speciation rate when the dispersal rate e uals unity.

E . 2. shows that as species dispersal between locations increases, allopatric

speciation becomes less likely. Finally, we consider local extinction rates to be

constant, because empirical phylogenies suggest they do not increase with diversity,

and we consider them e ual for the two locations 𝜇

= 𝜇

= 𝜇 for simplicity.

2

16 Inferring local diversity-dependence

hen the widespread species goes extinct on one location, it becomes an en- demic species. e call this evolutionary process range contraction . For widespread species, complete extinction can only occur by two consecutive local extinction events without species dispersal between these events, i.e. contraction followed by local extinction. Thus we do not allow global extinction, i.e. immediate com- plete extinction for widespread species which is in line with the geographic state speciation and extinction model Goldberg et al. .

Theoretically, it is possible to compute the likelihood of our model given a phy- logeny using the hidden Markov approach of Etienne et al. 2 12. However, be- cause we have to consider all the possible combinations of endemic and widespread species richness i.e. (𝑎, 𝑏, 𝑐) with 𝑎 endemic species on location A, 𝑏 endemic species on location B, and 𝑐 widespread species , not only for the lineages in the phylogeny, but also for now-extinct species, the state space of the model is huge leading to severe computational and numerical problems. Hence, our aim here is to explore whether the computationally manageable non-spatial model Etienne et al.

1 can be used for inferring diversity-dependence from phylogenies simulated under the spatial model.

Simulation

e simulated trees starting with two ancestral species, one in each region. e used the Gillespie algorithm Gillespie to calculate the waiting time between two evolutionary events; this time is exponentially distributed with the sum of all rates as parameter. The probability of each event occurring is proportional to its rate relative to the sum of rates. A speciation event produces a new species while an extinction event eliminates one existing species. Species dispersal and contraction do not change the number of species but alter the character of species, switching between endemic and widespread. The simulation is performed for a given amount of time the crown age and conditional on survival of the crown lineages i.e.

the simulation is restarted if one or both become extinct to guarantee that both ancestors have descendants at present after which the phylogenetic tree of the extant species is constructed from the history of events. Here we show a series of trees see Fig. 2.1 and Fig.S1-S2 in the Supporting Information for trees under various scenarios to be discussed below to demonstrate how trees are shaped under different parameter combinations.

e simulated the phylogenies under a variety of parameter values. To explore how the ecological limit to diversity affects the detection of the diversity-dependent signal, we designed three spatial scenarios differing in ecological limits: two scenar- ios with identical limits on each location Scenario 1: 𝐾 = 20, Scenario 2: 𝐾 = 40 , and one scenario with different ecological limits Scenario :𝐾 = 20, 𝐾 = 40 . For comparison with the non-spatial model, we additionally simulated two non-spatial scenarios differing in ecological limit Scenario 4: 𝐾 = 20 and Scenario : 𝐾 = 40 . e assumed a crown age of 15 time units, which can be interpreted as 15 million years. e fixed the values for the intrinsic speciation rates:

𝜆

= 𝜆

= 0.8, 𝜆

= 0.2.

Inferring local diversity-dependence

2

17

11 11 11 11 11 11 11 11 11 11 22 22 22 22 22 22 22 22 22 22

11 11 11 11 12 11 11 12 11 22 11 11 22 22 22 22 22 22 22 22

12 12 11 21 22 21 22 22 11 11 21 22 11 21 11 23 12 21 21 2

33 12 33 12 11 23 12 33 12 31 23 12 23 33

33 3 33 3 33 3 33 3 33 3 33 3 33

11 11 11 11 11 11 11 11 12 22 22 22 22 22 22 22 22 22

22 11 11 12 11 11 11 11 13 22 21 22 21 12 22 22 2

32 12 22 12 21 11 12 22 21 13 32 31 21 11 33 13

21 31 23 12 12 33 11 13 12 23 31 21 32 23 22

33 3 33 3 33 3 33 3 33 3 33 3 33

11 11 11 11 11 11 11 12 22 22 22 22 22 22

22 21 11 11 11 12 22 11 21 11 11 22 22 22 12 2

21 32 11 11 12 22 22 31 33 22 11 12 22 22 11 13

13 31 32 33 33 13 12 31 33 21 31 22

33 3 33 33 3 33 3 13 33 3 33 3 32

Endemic species Endemic species Widespread species

µ = 0 µ = 0.1 µ = 0.2 µ = 0.4

M

= 1000 M

= 5 M

= 1 M

= 0.15 M

= 0

Figure 2.1: Examples of phylogenetic trees produced in Scenario 1. Because the trees for migration

rates between and 1 are very similar, we only display values of extinction , . , , , .

The branches are colored by the location of species. Sympatric speciation and allopatric speciation are

also distinguishable by the color of the nodes and the daughter species.

2

18 Inferring local diversity-dependence

e looked at the same set of extinction rates as in Etienne et al. 1, 2 : 0, 0.1, 0.2, 0.4. Finally, we studied the behavior of the model and the inference under a gradient of intrinsic dispersal rates: 𝑀 = 0, 0.05, 0.1, 0.15, 0.3, 0.5, 1, 5, 1000. The case 𝑀 = 0 corresponds to a birth-death process occurring on two independent locations. As 𝑀 increases, the model tends towards the non-spatial model with one important difference, see Results and species at the tips become increasingly widespread species. In all, we simulated 36 parameter sets for each scenario. For each parameter set, we generated 100 phylogenetic trees.

Inference

e applied a bootstrap likelihood ratio test Etienne et al. 2 , Gudicha et al.

1 , Tekle et al. 2 1 to the simulated data to determine the power of the non- spatial model to detect diversity-dependence in the spatial model. The 𝜒 likelihood ratio test cannot be used due to the mismatch between type I error rate and the significance level used as reported in Etienne et al. 2 . The bootstrap likelihood ratio test Etienne et al. 2 proceeds as follows:

1. Collect an empirical data set of phylogenetic branching times. One can also simulate data under another model for a specific parameter set which was the case for our study, where we simulated under the spatial model . 2. Estimate from these data the maximum-likelihood ML parameters under the

constant-rates CR model and the diversity dependent DD model the non- spatial model . Then calculate the likelihood ratio which is denoted by 𝐿𝑅 . . Generate a bootstrap sample by simulating 𝑋 data sets under the CR model

using the parameter estimates obtained for the CR model in step 2.

4. For each of these 𝑋 simulated CR data sets, estimate the parameters under the CR model as well as the DD model and compute the likelihood ratio 𝐿𝑅 for data set 𝑖 .

. Compare the observed 𝐿𝑅 with the distribution of 𝐿𝑅 -values (𝑖 = 1..𝑋 from the bootstrap simulations. Count the number of simulations with 𝐿𝑅 larger than 𝐿𝑅 and denote the number by 𝑅 . The p- value of the test is defined as (𝑅 + 1)/(𝑋 + 1).

. A significance level 𝛼 e.g. 0.05 is set to accept or reject the CR model by comparison with the p- value. Record the 𝐿𝑅 associated with this 𝛼, 𝐿𝑅 . . To assess the power of the test, simulate 𝑋 times under the DD model with

the ML parameters estimated under the DD model in step 2.

. For these 𝑋 data sets simulated in step , estimate parameters under both CR and DD model and compute the 𝐿𝑅 for each data set.

. The larger the number of the likelihood ratios exceeding 𝐿𝑅 , the clearer is the

signal of diversity-dependence. Denote the number of the 𝑋 simulations

Inferring local diversity-dependence

2

19

where the 𝐿𝑅 is larger than 𝐿𝑅 by 𝑅 . Define the power of the test by 𝑅 /(𝑋 + 1).

e performed this method for all the parameter sets. e thus have 36 param- eter sets of 100 simulations each with 2000 bootstrap samples, totaling .2 million simulations and parameter estimations for each scenario. Given that each parame- ter estimation takes a few minutes, the total computation time for scenarios was to 1 million minutes, roughly, 1 - 2 years on a single computer. Hence, we performed these calculations on a high-performance computing cluster, but even then computational time was substantial. e therefore provide all simulations and data as supplementary material.

Results Model Behavior

To study how the model behaves under different dispersal and extinction rates, we plotted the species-through-time STT plots that include both extant and ex- tinct species under different 𝐾 settings see Fig.2.2 for Scenario 1 and Fig.S -S , Supporting Information, for other scenarios . The STT plots show how the to- tal number of species changes due to macro-evolutionary events. The STT plots that we show here are for a single location because in our model the diversity- dependence is defined as local dynamics. e also plotted the non-spatial STT plots tracking the total number of species in the system, i.e. for both locations together as supplementary results Supporting Information, see Fig.S1 -S12 . As expected, from the local STT plots we observed a positive correlation between species disper- sal and species richness and a negative correlation between extinction and species richness. However, in the non-spatial STT plots dispersal seems to have a complex influence on the global species richness. Although the effect of dispersal is small, it gets larger with increasing extinction rate. e will discuss it later in the section of parameter estimation. To test the model behavior under high species dispersal rate, we additionally explored an extreme case where dispersal rate is extremely large 𝑀 = 1000 . In this case, all parameter settings varying only in extinc- tion rates lead a similar increasing pattern in species richness and the diversity in both locations reach the ecological limit rapidly for example, 𝐾 = 20 for Scenario 1, see Fig.2.2 and Fig.S -S , Supporting Information, for other scenarios . This phenomenon is similar to a pure birth process due to the extremely high dispersal rate. The biological explanation is that once an endemic species is produced, it spreads out to the other location immediately which makes it almost impossible to go globally extinct. Therefore, the system is filled with widespread species and a few endemic species at the e uilibrium level, which is identical to the ecological limit.

Furthermore, we studied lineages-through-time LTT plots for extant species

for both locations together, which allows comparison with LTT plots from the non-

spatial model. e observed a pattern of an early burst and the pull of the present

ubo and Iwasa 1 1 , Nee et al. 1 Fig.2.4 for Scenario 1 and Fig.S1 -S14

2

20 Inferring local diversity-dependence

µ = 0 µ = 0.1 µ = 0.2 µ = 0.4

Number of lineages

2 105 2040 80

2 105 2040 80

2 5 1020 40 80

2 5 1020 4080

2 105 20 4080

2 105 2040 80

2 105 2040 80

2 5 1020 40 80

2 105 2040 80

−15 −10 −5 0−15 −10 −5 0−15 −10 −5 0−15 −10 −5 0

M

= 0

M

= 0.05

M

= 0.1

M

= 0.15

M

= 0.3

M

= 0.5

M

= 1

M

= 5

M

= 1000

Time

Figure 2.2: Species-through-time STT plots that include extinct species for one location across all

parameter settings of Scenario 1. Lower extinction accelerates species accumulation. Species dispersal

increases the number of species at e uilibrium. The dashed line at value shows the input value of

. The black line denotes the median STT plot, the gray shading represents the uantiles minimum,

2. th percentile, 2 th percentile, th percentile, . th percentiles, maximum .

Inferring local diversity-dependence

2

21

p−valuePower

0.000.05 0.25 0.50 0.75 1.00

0 0.05 0.1 0.15 0.3 0.5 1 5 1000 0 0

µ = 0

0.00 0.25 0.50 0.75 1.00

0 0.05 0.1 0.15 0.3 0.5 1 5 1000 0 0 0.000.05 0.25 0.50 0.75 1.00

0 0.05 0.1 0.15 0.3 0.5 1 5 1000 0 0

µ = 0.1

0.00 0.25 0.50 0.75 1.00

0 0.05 0.1 0.15 0.3 0.5 1 5 1000 0 0 0.000.05 0.25 0.50 0.75 1.00

0 0.05 0.1 0.15 0.3 0.5 1 5 1000 0 0

µ = 0.2

0.00 0.25 0.50 0.75 1.00

0 0.05 0.1 0.15 0.3 0.5 1 5 1000 0 0 0.000.05 0.25 0.50 0.75 1.00

0 0.05 0.1 0.15 0.3 0.5 1 5 1000 0 0

µ = 0.4

0.00 0.25 0.50 0.75 1.00

0 0.05 0.1 0.15 0.3 0.5 1 5 1000 0 0 Scenarios

S1 S4 S5

Dispersal rate

Figure 2. : -values and powers of the test of spatial scenario 1 and non-spatial scenarios 4 and : as the dispersal rate increases, the -value declines approaching while the power of the test rises up to . The signal of diversity-dependence tends to be detected with high dispersal and low extinction.

Especially, in the case of the pure birth process, all the scenarios show such a strong signal that the distribution bars of -values and powers are compressed to thick black lines. hen extinction rate is . , diversity dependence is not detected statistically until dispersal rate reaches . In the box plots,

thick solid lines, boxes and whiskers denote the percentiles of , and , respectively.

for other scenarios , except for the highest extinction rate 𝜇 = 0.4 and lowest dispersal rate 𝑀 = 0 , for which the shape of the LTT plot approaches a straight line.

Detecting diversity-dependence

Diversity-dependence can be detected with high power except when extinc- tion is high larger than 0.4 and species dispersal is low smaller than 1 at the significance level 𝛼 = 0.05 see Fig.2. for Scenario 1 and Fig.S -S , Support- ing Information, for other scenarios . This suggests that extinction tends to erase the signature of diversity-dependence while species dispersal strengthens the sig- nal. hen relating this to the STT and LTT plots, we observe that weak signals of diversity-dependence are accompanied with a low rate of species accumulation. In contrast, strong evidence for diversity dependence often occurs for low extinction and high dispersal. Both of these situations lead to intense species interactions. e also observe substantial early bursts for LTT plots whenever diversity dependence is detected.

To explore whether the diversity-dependent signal would be stronger in the

2

22 Inferring local diversity-dependence

scenario that has a higher ecological limit to diversity, we studied the power of the test for different scenarios with different ecological limits. Table 2. shows power to detect diversity-dependence under different parameter combinations of three spatial scenarios. e observe that systems with a higher ecological limit to diversity show a broader range of high detection power in parameter space. In particular, the scenario with distinct limits 𝐾 = 20, 𝐾 = 40 on two locations shows an intermediate strength of diversity-dependence between two scenarios with identical limits, stronger than Scenario 1 𝐾 = 20 and weaker than Scenario 2 𝐾 = 40 .

e next explored whether the partition of the community into two locations would weaken the strength of the diversity-dependent signal. The non-spatial sce- narios 4 and have the same value of ecological limits as the spatial scenarios 1 and 2, respectively, but constrain the species diversification to only one single location. The spatial structure indeed affects the diversity-dependence detection but in a complex manner see Fig.2. and Fig.S -S , Supporting Information, for other spatial scenarios . hen the locations are more isolated, i.e. they have little species interaction between them, the non-spatial scenarios show stronger diversity-dependence than the spatial scenarios. hen dispersal rate increases, this pattern is reversed: because species dispersal reduces extinction thus leading to a high rate of species accumulation.

Parameter estimate accuracy and precision

The performance of parameter estimation depends strongly on the extinction and dispersal rates. Accurate parameter estimations are obtained for low extinction and dispersal. The median estimates for the ecological limit are around the sum of the local limits Fig.2. , Fig.S -S4, Supporting Information, for other scenarios when both extinction and dispersal rates are low. But bias in parameter estimates increases for larger dispersal and extinction rates. This is due to the fact that both dispersal and extinction strongly control the species richness of the system.

Extinction has a negative effect on diversity so we find that our estimate of the ecological limit decreases with increasing extinction. The influence of dispersal on species richness is more complex. On the one hand, dispersal promotes the conversion of endemic species to widespread species thereby decreasing species richness. On the other hand, dispersal reduces extinction and thereby increases species richness. e observe this phenomenon in our simulation study, especially for high extinction. In all scenarios, the estimates of the ecological limit increase at first but then drop with the dispersal rate increases. This also explains the pattern that the e uilibrium of species richness in the non-spatial STT plots first increases and then declines with increasing dispersal rate.

Speciation and extinction estimates are robust when both extinction and dis-

persal rates are low. However, when species dispersal increases the speciation

estimates are biased upward while extinction is biased downward. Interestingly,

the speciation estimates are biased up to a value e ual to the sum of the local

speciation rates of the two locations. The extinction estimates are biased down to

zero which agrees with the explanation that dispersal reduces extinction.