How Ecology and Landscape Dynamics Shape ...

Viewer
Transcript

Systematic Biology Advance Access published May 12, 2015 Syst. Biol. 0(0):1–18, 2015 © The Author(s) 2015. Published by Oxford University Press, on behalf of the Society of Systematic Biologists. All rights reserved. For Permissions, please email: [email protected] DOI:10.1093/sysbio/syv014

How Ecology and Landscape Dynamics Shape Phylogenetic Trees FANNY GASCUEL1,2,3,∗ , R ÉGIS FERRIÈRE3,4 , R OBIN AGUILÉE5 , AND AMAURY LAMBERT1,6 1 Center

for Interdisciplinary Research in Biology, CNRS UMR 7241, Collège de France, Paris, France; 2 Sorbonne Universités, UPMC Univ Paris 06, CNRS UMR 7625, Laboratoire Écologie et Évolution, Paris, France; 3 Institut de Biologie de l’Ecole Normale Supérieure, CNRS UMR 8197, École Normale Supérieure, Paris, France; 4 Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, AZ USA; 5 Laboratoire Évolution et Diversité Biologique, CNRS UMR 5174, Université Paul Sabatier, Toulouse, France; and 6 Sorbonne Universités, UPMC Univ Paris 06, CNRS UMR 7599, Laboratoire Probabilités et Modèles Aléatoires, Paris, France ∗ Correspondence to be sent to: Center for Interdisciplinary Research in Biology, Equipe Stochastic Models for the Inference of Life Evolution, Collège de France, 11 place Marcelin Berthelot, 75005 Paris, France; E-mail: [email protected]. Received 23 July 2014; reviews returned 23 September 2014; accepted 9 March 2015 Associate Editor: Susanne Renner

How species diversity originates and how it changes across multiple scales of time and space is a perennial question in evolutionary biology. To infer the history of species diversiﬁcation, reconstructed phylogenies can be ﬁtted by mathematical models of tree growth, such as the pure-birth process (also called Yule process; Yule 1925), which assumes that the species tree is only generated by constant-rate speciation; and the birth–death process (Kendall 1948; Raup et al. 1973), which assumes that the species tree is generated by speciation and extinction events occurring at constant rates per lineage. Because the Yule and birth–death processes assume constant evolutionary rates, they are referred to as null models of phylogenetic tree growth. To evaluate the hypothesis that a given clade diversiﬁed at constant speciation and extinction rates, the shape of the clade’s reconstructed phylogeny (i.e., tree of extant species) can be compared to that of trees generated under null models (Mooers and Heard 1997). There are two main properties that characterize phylogenetic tree shape: branching tempo and (im)balance. Branching tempo refers to the distribution of relative branching times from the root

to the tips of the tree; it thus reﬂects variation in the rate of lineage accumulation through time. Imbalance refers to the distribution of branchings among subclades: diversiﬁcation rates that differ between subclades can lead to unbalanced trees, while homogeneous rates lead to balanced trees. Empirical phylogenies have been shown to often differ from trees generated by the null models both in terms of branching tempo and balance. Real phylogenies tend to be unbalanced (e.g., Guyer and Slowinski 1991; Heard 1992; Guyer and Slowinski 1993; Slowinski and Guyer 1993; Mooers 1995; Purvis 1996; Mooers and Heard 1997; Blum and François 2006), and branching tends to slow down (e.g., Nee et al. 1992; Zink and Slowinski 1995; Lovette and Bermingham 1999; Pybus and Harvey 2000; Rüber and Zardoya 2005; Kozak et al. 2006; Seehausen 2006; Weir 2006; McPeek 2008; Phillimore and Price 2008; Rabosky and Lovette 2008a; Jønsson et al. 2012). Statistical models of phylogenetic trees have been developed to characterize variations in evolutionary rates of speciation and extinction, through time and/or among subclades (Sanderson and Donoghue 1996; Ricklefs 2007; Stadler 2011; Pyron and Burbrink

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

Abstract.—Whether biotic or abiotic factors are the dominant drivers of clade diversiﬁcation is a long-standing question in evolutionary biology. The ubiquitous patterns of phylogenetic imbalance and branching slowdown have been taken as supporting the role of ecological niche ﬁlling and spatial heterogeneity in ecological features, and thus of biotic processes, in diversiﬁcation. However, a proper theoretical assessment of the relative roles of biotic and abiotic factors in macroevolution requires models that integrate both types of factors, and such models have been lacking. In this study, we use an individualbased model to investigate the temporal patterns of diversiﬁcation driven by ecological speciation in a stochastically ﬂuctuating geographic landscape. The model generates phylogenies whose shape evolves as the clade ages. Stabilization of tree shape often occurs after ecological saturation, revealing species turnover caused by competition and demographic stochasticity. In the initial phase of diversiﬁcation (allopatric radiation into an empty landscape), trees tend to be unbalanced and branching slows down. As diversiﬁcation proceeds further due to landscape dynamics, balance and branching tempo may increase and become positive. Three main conclusions follow. First, the phylogenies of ecologically saturated clades do not always exhibit branching slowdown. Branching slowdown requires that competition be wide or heterogeneous across the landscape, or that the characteristics of landscape dynamics vary geographically. Conversely, branching acceleration is predicted under narrow competition or frequent local catastrophes. Second, ecological heterogeneity does not necessarily cause phylogenies to be unbalanced—short time in geographical isolation or frequent local catastrophes may lead to balanced trees despite spatial heterogeneity. Conversely, unbalanced trees can emerge without spatial heterogeneity, notably if competition is wide. Third, short isolation time causes a radically different and quite robust pattern of phylogenies that are balanced and yet exhibit branching slowdown. In conclusion, biotic factors have a strong and diverse inﬂuence on the shape of phylogenies of ecologically saturating clades and create the evolutionary template in which branching slowdown and tree imbalance may occur. However, the contingency of landscape dynamics and resource distribution can cause wide variation in branching tempo and tree balance. Finally, considerable variation in tree shape among simulation replicates calls for caution when interpreting variation in the shape of real phylogenies. [adaptive radiation; allopatric speciation; competition; eco-evolutionary feedbacks; ecological speciation; geographic isolation; Macroevolution; phylogeny]

1

[18:29 7/5/2015 Sysbio-syv014.tex]

Page: 1

1–18

2

SYSTEMATIC BIOLOGY

[18:29 7/5/2015 Sysbio-syv014.tex]

(e.g., MacArthur 1972; Rosenzweig 1978; Pimm 1979; Walker and Valentine 1984; Nee et al. 1992; Lovette and Bermingham 1999; Weir 2006; Phillimore and Price 2008; Jønsson et al. 2012). Alternative non-ecological hypotheses (reviewed in Moen and Morlon 2014) include an increase in genetic constraints over time (Gavrilets and Losos 2009) and a nonrandom acceleration of extinction (Zink and Slowinski 1995, although the latter may have an ecological mechanistic underpinning). Ecological niche ﬁlling is supposed to slow down diversiﬁcation through an increase in the number of species competing for limited resources (Mayr 1942, 1947; Simpson 1953; Schluter 2000); by contrast, in the early stage of diversiﬁcation the limited number of competitors is expected to strengthen directional and disruptive selection (Schluter 2000, 2001). Furthermore, spatial heterogeneity in the availability and width of ecological niches is expected to impact species diversity in subclades evolving locally in space, and thus to explain the imbalance of the species trees (e.g., Simpson 1944; Rosenzweig 1995; De Queiroz 2002; Purvis et al. 2011). Alternative hypotheses for tree imbalance include the effect of key innovations (e.g., Simpson 1953; Heard and Hauser 1995; Phillimore et al. 2006; Jønsson et al. 2012; Magnuson-Ford and Otto 2012), sexual selection (Darwin 1871; Barraclough et al. 1995; Mitra et al. 1996), and demographic asymmetries between mother and daughter lineages at speciation (Davies et al. 2011). Does the ubiquitous pattern of phylogenetic branching slowdown and imbalance call for ecological niche ﬁlling and spatial heterogeneity in ecological features as the main explanation, and does this lend weight to the biotic (Red Queen) view of species diversiﬁcation? Here we aim at testing in silico: (i) whether ecological niche ﬁlling is the general explanation for the tendency of branchings to decelerate in phylogenetic trees; (ii) whether wide competition limiting the number of available ecological niches accentuates this slowdown; and (iii) whether spatial heterogeneity in ecological features is the general explanation for the tendency of clade diversiﬁcation to proceed at different rates in different subclades, thus generating unbalanced phylogenetic trees. We use an individual-based model of adaptive radiation driven by ecological speciation. Individuals are characterized by their ecological phenotype, which drives exploitative competition; heritable variation in the ecological phenotype fuels the processes of withinspecies local adaptation (anagenesis) and speciation (cladogenesis) in a dynamic landscape where local communities can split and fuse repeatedly through evolutionary time, thus allowing not only allopatric speciation, secondary contact, character displacement, reinforcement, but also hybridization. This dynamic landscape allows us to investigate how the effects of biotic factors on macroevolutionary patterns are altered by changes in the abiotic environment that occur on longer timescales (e.g., climatic changes or tectonic activity). The latter are indeed also recognized to inﬂuence diversiﬁcation (Keast 1961; Cracraft 1986;

Page: 2

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

2013). These models, which consider either diversitydependent (Rabosky and Lovette 2008b; Etienne et al. 2012), trait-dependent (Maddison et al. 2007; FitzJohn 2010), or time-dependent evolutionary rates (Alfaro et al. 2009; Morlon et al. 2011; Stadler 2011), ﬁt empirical patterns of diversiﬁcation better than the null, constant-rates models. However, the question of the mechanisms causing changes in evolutionary rates through time (possibly mediated by change in diversity) or among clades (possibly mediated by character differentiation), remains out of reach of such models that are phenomenological in essence (Nee 2006; Pyron and Burbrink 2013). A major challenge is now to identify what mechanisms are responsible for variations in evolutionary rates through time and across clades. More mechanistic models of evolutionary diversiﬁcation, either species- or individual-based, have recently been developed (e.g., Doebeli and Dieckmann 2003; Gavrilets and Vose 2005; McPeek 2008; Aguilée et al. 2011; Pontarp et al. 2012; Aguilée et al. 2013) with the goal of investigating how speciﬁc components of the diversiﬁcation process inﬂuence phylogenetic tree shape. Geographical structure involved in allopatric speciation (Pigot et al. 2010), demographic asymmetry at speciation (comparing point mutation, where new species arise as singleton lineages, and random ﬁssion, where ancestral species are randomly cleaved into two daughter species; Davies et al. 2011), and ecological differentiation at speciation (McPeek 2008), have been shown to inﬂuence the shape of reconstructed species trees; conditions were found for these processes that can lead to branching slowdown and/or unbalanced topology. However, these models describe the process of speciation only phenomenologically, as opposed to speciation being an emergent phenomenon driven by underlying genetic and ecological mechanisms at the individual level. Thus, in these models the dynamics of speciation (rate, location) were assumed, rather than being an output of the model. A critical role for individual-level ecological processes in species diversiﬁcation has long been envisioned and discussed ever since Darwin (Darwin 1859; Simpson 1953; Stanley 1979; Grant 1986; Schluter 2000): interactions between organisms and their environment may drive directional selection (e.g., Grant 2003), and interactions among organisms may generate disruptive selection (through inter- or intraspeciﬁc competition, or through predation; e.g., Dieckmann and Doebeli 1999; Bailey et al. 2013). According to the Red Queen view of evolution (Van Valen 1973), ecological processes may even drive diversiﬁcation. Thus, there is a need to develop macroevolutionary models that explicitly take into account ecological processes, and this has been identiﬁed as a major challenge at the interface of evolutionary and ecological theory (Pigot et al. 2010; McInnes et al. 2011; Fritz et al. 2013; Pyron and Burbrink 2013; Wang et al. 2013; Morlon 2014). The conventional view is indeed that ecological niche ﬁlling is responsible for branching slowdown

1–18

2015

GASCUEL ET AL.—ECOLOGY, LANDSCAPE DYNAMICS AND SPECIES TREE SHAPE

METERIALS AND METHODS Model Overview We use Aguilée et al.’s (2013) model of diversiﬁcation. Hereafter, we highlight key features of the model; for a fuller description, see Aguilée et al. (2013). The model describes the dynamics of a community of competitively interacting species that inhabit a patchy landscape. Organisms are diploid and reproduce sexually. The model tracks the population dynamics of each species and their spatial distribution in a patchy landscape that varies in time and space. The model integrates biotic processes that operate at the level of individuals: individual birth and death, local interactions between individuals, dispersal; and abiotic processes that operate at the level of patches: fragmentation or fusion of extant patches that drive the landscape dynamics; heterogeneity in resource distribution that translates into among-patch variation in carrying capacity. The community history of species origination and species loss emerges from the joint operation of the individual-level biotic and patch-level abiotic processes. Local adaptation drives anagenesis, and cladogenesis is driven by the evolution of reproductive isolation. The community can lose species because of extinctions, which can be due to demographic stochasticity, competitive exclusion, hybridization, or environmental catastrophes. Thus, speciation and extinction rates are

[18:29 7/5/2015 Sysbio-syv014.tex]

not predeﬁned; they are model outputs. From one ancestral species the model predicts the growth of the whole species tree, and how the tree is shaped by the underlying processes and their parameters. Individual-level Processes We assume that resources vary along gradients on two spatial dimensions, and that the resource utilization strategy of an individual is described by two ecological traits, x1 and x2 , one for each spatial dimension. The pair (x1 , x2 ) is the individual’s “ecological phenotype”. Individuals also express a ‘choosiness trait’ a, which measures the degree to which ecological phenotypic similarity inﬂuences the probability of mating; choosiness can favor mating between either more similar or more dissimilar ecological phenotypes. All three traits are genetically variable. They are determined by Lk diploid, additive quantitative loci on autosomal chromosomes, where k denotes x1 , x2 , or a. We use Lk = 16, which is large enough to make sympatric speciation unlikely (Gavrilets 2003; Coyne and Orr 2004; Waxman and Gavrilets 2005). Local (i.e., within each patch) population growth is assumed to be logistic, with competition for resources being both locally density-dependent (Roughgarden 1972) and frequency-dependent, that is, stronger between individuals with more similar ecological traits (Dieckmann and Doebeli 1999; Doebeli and Dieckmann 2003). The birth rate of each individual is assumed to be constant. The death rate depends on the local carrying capacity and on the strength of local competition. In each geographical site, the carrying capacity function, which models a continuous distribution of resources, is a Gaussian function with maximum K at phenotype (x1∗ , x2∗ ), called the ecological optimum of the site, and standard deviation K . The strength of local competition is the sum of the competition strengths with all other individuals of the same patch (or set of connected patches - see below); competition strength between each pair of individuals is also given by a Gaussian function of the Euclidian distance between their ecological phenotypes, with maximum value one when the distance is zero, and standard deviation C . The standard deviation of the carrying capacity function, K , is chosen to be larger than that of the competition kernel, C , to ensure that ecological optima would be “evolutionary branching points” in the corresponding trait-substitution sequence model for an asexual population (Geritz et al. 1998; Dieckmann and Doebeli 1999; Vukics et al. 2003; Doebeli and Ispolatov 2010). Gaussian kernels have been widely used to model competition, but also criticized for introducing structural instability in deterministic models of resource use evolution (Barton and Polechová 2005; Gyllenberg and Meszéna 2005; Leimar et al. 2008). By including multiple sources of stochasticity (in individual life history, mutation, dispersal, and landscape dynamics), our results should be immune to the spurious effects of strictly Gaussian competition kernels.

Page: 3

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

Barnosky 1999; Ezard et al. 2011), by striking lineages at random and potentially causing speciation or extinction (“Court Jester” hypothesis, Barnosky 1999). We consider three aspects of environmental changes formulated in terms of landscape dynamics. First, we explore the role of the pace of landscape dynamics, involving the recurrent fusion and fragmentation of geographical sites over time. Landscape dynamics are known to play a major role in diversiﬁcation (Mayr 1942; Greenwood 1965; Rossiter 1995; Aguilée et al. 2013), due to the recurrence of ecological divergence in allopatry followed by reinforcement or character displacement at secondary contact. Second, we investigate the effect of the time in isolation, measured as the expected proportion of time that two adjacent sites spend in isolation of one another. Third, we address the inﬂuence of the frequency of local catastrophes that strike geographical sites at random and cause the extinction of all populations living in that site. Local catastrophes tend to keep metapopulation size below carrying capacity, hence potentially limiting the ecological and evolutionary effects of competition. These analyses provide new insights into the causal relationship between ecological niche ﬁlling and spatial environmental heterogeneity on the one hand, and the branching tempo and balance of phylogenetic trees on the other. By comparing our theoretical results to existing empirical data, we discuss how biotic and abiotic factors may interact and together shape the branching pattern of real phylogenetic trees.

3

1–18

4

SYSTEMATIC BIOLOGY

with close ecological traits (x1 , x2 ), as in Aguilée et al. (2013). Then, we cluster these populations into sets of independently evolving metapopulations by calculating the probability of reproductive isolation between each pair of populations (Mayr 1963): each set of populations that can interbreed (probabilities of reproductive isolation below a threshold thrri = 99%) but are reproductively isolated from other populations (probabilities of reproductive isolation above thrri ), is considered to be a species. We compute probabilities of reproductive isolation using the same function as that used for mating between individuals, but using the average values of traits in each population, x1 , x2 , and a. This approximation reduces computation time signiﬁcantly while having no major effect on our results (online Appendix 2 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51). To follow the evolutionary history of species, at each time step t (every 100 generations, with similar results as 10 generations - see online Appendix 3 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51) we determine (i)

Landscape Dynamics The metacommunity evolves in a structured landscape, made of n×n sites distributed over a two-dimensional grid. Each site is characterized by a different ecological optimum (x1∗ , x2∗ ), with x1∗ and x2∗ increasing by a constant amount x between adjacent sites in either direction. The landscape is dynamic in the sense that sites undergo fragmentation and fusion, causing the repeated alternation of allopatry and sympatry. Fragmentation is due to barriers that appear between adjacent sites at rate f ; fusion is due to extant barriers disappearing at rate c. When a set of sites is connected, migration occurs at birth as in the island model (Hanski and Gaggiotti 2004), with equal probability for the offspring to reach any of the connected sites. We, however, also investigated results under another version of the model, where migration is more likely toward the sites where the offspring’s parents are located than toward other sites, but this led to similar phylogenetic tree shapes (see online Appendix 1 available as Supplementary Material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51). Finally, our model incorporates local catastrophes that strike sites at random and cause the extinctions of all populations present in the site. The rate of local catastrophe is denoted by cr. Delineating Species and Following Their Evolutionary History We delineate species using De Queiroz’s (2007) uniﬁed concept, considering species as “separately evolving metapopulation lineages”—metapopulations (sets of connected populations) that evolve through time independently from other populations. We ﬁrst identify populations by grouping all individuals

[18:29 7/5/2015 Sysbio-syv014.tex]

the Nt existing species [St ]i∈{1,...,Nt } as explained above, and then determine how they descend from the species (j)

existing at the previous time step, [St−1 ]j∈{1,...,Nt−1 } . Descent relationships are established by comparing the phenotypic traits of species and populations at time t to those at time t−1, and by minimizing Euclidian distances between these average traits. Our method ensures that: (i) each species keeps the same identity as long as it remains reproductively isolated from other species, even if its phenotype evolves; (ii) speciation can happen through differentiation of a species into various reproductively isolated populations; (iii) extinction can happen when a population reaches a very small size or through species hybridization; and (iv) descent relationships during hybridization are consistent with probabilities for each species to transmit its genotypes to the future lineage. Further details on how we delineate species and follow their evolutionary history can be found in the online Appendix 4 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51.

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

Reproduction is sexual. We use the same mating probability function as in Aguilée et al. (2013): a Gaussian function with the minimum biological realism required, where mating probability depends on the similarity in the ecological phenotypes of the two partners and on the choosiness of the more choosy one (either toward similarity or dissimilarity). If mating is rejected, another potential mate is drawn at random and the process is repeated until mating succeeds, or until 50 potential mates have been tried. The cost of choosiness is thus assumed to be small (Schneider and Bürger 2006; Kopp and Hermisson 2008), facilitating the evolution of assortative mating. When mating occurs, the offspring sex is determined randomly assuming a balanced sex ratio; the offspring ecological and choosiness traits are determined from the random independent segregation of the Lk parental loci that code for each trait k (k = x1 , x2 , or a). Mutation occurs at each locus with probability k and the mutant allelic value is drawn from a normal distribution with mean equal to the parental allelic value and with standard deviation k 2Lk .

Simulating Phylogenies and Analysing Their Shape Diversiﬁcation was initiated from one ancestral species, spread over all geographical sites (150 individuals per site). The ecological and choosiness traits (x1 , x2 , and a) of this species were determined by Lk , k ∈ {x1 ,x2 ,a}, diploid loci randomly chosen in a Gaussian distribution centered at zero. The ancestral species was, therefore, adapted to the geographical site with intermediate ecological optimum (x1∗ , x2∗ ) in the two dimensions of the environmental gradient. This prevents any border effects in the initiation of diversiﬁcation, and represents a wide range of empirical situations leading to eco-evolutionary diversiﬁcation.

Page: 4

1–18

2015 TABLE 1.

5

GASCUEL ET AL.—ECOLOGY, LANDSCAPE DYNAMICS AND SPECIES TREE SHAPE Parameter values and symbols used in this study, except if a different value is given in the legends

Parameter

Deﬁnition

Default value

r Lk k 2 x 2 a

Per-capita birth rate Number of loci coding each trait k Mutation probability at each locus Expected phenotypic variance of trait x Expected phenotypic variance of trait a

K K C

Carrying capacity at the ecological optimum in each site Standard deviation of the maximum carrying capacity per site Standard deviation of the competition kernel

150 (individuals) 1.0 0.4

n2 x f c cr

Number of sites in the landscape Difference in optimal trait values between adjacent sites Rate of border appearance Rate of border disappearance Rate of catastrophic local community extinction

9 1.0 10−3 (birth rate) 5.10−5 (birth rate) 0 (birth rate)

[18:29 7/5/2015 Sysbio-syv014.tex]

differentiated species never hybridize (if they do so, rapid hybrid collapse always occurs as a result). Steady state in number of individuals and number of species was typically reached after 30,000 to 50,000 generations. This is much more rapid than typically observed in adaptive radiations (Schluter 2000), and reﬂects high mutation rates (see Table 1; similar as in previous individual-based models of diversiﬁcation, Dieckmann and Doebeli 1999; Gavrilets and Vose 2005; Pontarp et al. 2012; Aguilée et al. 2013). Such high mutation rates make it possible to observe species diversiﬁcation in reasonable computational time, and are not expected to affect the overall species tree shapes, at least in this model. On the one hand, the speed of phenotypic trait evolution depends on the product of the mutation rate with the phenotypic variance, and thus similar diversiﬁcation patterns can be observed with smaller mutation rates and higher expected phenotypic variance (e.g., with 2 = 2 = 0.01, 2 = 0.04, and x1 = x2 = a = 10−5 , x 2 a x1 Lx1 = Lx2 = La = 6; Aguilée et al. 2013). On the other hand, patterns of species diversiﬁcation depend on how the rate of phenotypic trait evolution scales with the rate of landscape dynamics (Claessen et al. 2007; Aguilée et al. 2011): similar patterns of diversiﬁcation can also be observed on longer temporal scales with smaller mutation rates given smaller rates of landscape dynamics. At steady state, most simulated metacommunities contained 12 to 70 different species for a total of 1800 to 9000 individuals. For each full-run simulation, we extracted the tree of all extant species every 100 generations and characterized its shape both in terms of branching tempo and balance. We analysed the branching tempo of phylogenetic trees by computing the statistic, which tests for departures from a time-homogeneous, lineagebased pure-birth process of diversiﬁcation (Cox and Lewis 1966; Pybus and Harvey 2000). In “tippy” trees, in which nodes are concentrated close to the present (due to an increase in the rate of lineage accumulation toward the present), is positive; in “stemmy” trees, in

Page: 5

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

Indeed, once a key innovation (novel phenotypic trait that brings the ability to invade new ecological niches) arose, all the surrounding landscape may potentially offer space and ecological niches for new species to evolve, and thus many species may originate from the ancestral species that acquired the innovation (Simpson 1944, 1953; MacArthur 1972; Walker and Valentine 1984; Schluter 2000). This sequence of events may also occur following the colonization of new geographical areas with few ecological competitors. We nonetheless studied the inﬂuence of border effects, by considering an ancestral species adapted to marginal ecological conditions: results show that border effects also inﬂuence phylogenetic tree shape (see online Appendix 5 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51). For each set of parameters, we simulated the process of diversiﬁcation for 100,000 generations, running 50 replicates. Not all simulations resulted in successful adaptive radiations: in up to 60% of them, diversiﬁcation was initiated but failed at some point in time due to the hybridization of all species. This diversity collapse results from the hybridization of two ecologically close species, which increases their phenotypic variance and decreases their level of assortative mating, triggering a cascade of hybridization events throughout the metacommunity. Hybrid collapses have been observed in previous models (Gilman and Behm 2011; Aguilée et al. 2013), and are suspected to occur in nature (Seehausen 2004; Vonlanthen et al. 2012). Here, they are expected to occur at some (possibly very large) time in every single simulation, because of the absence of any mechanism that would prevent species to hybridize. In nature, species are expected to evolve irreversible genetic incompatibilities, which prevent them from hybridizing. However, accounting for such long-term reproductive isolation mechanism does not have any noticeable effect on our results (see online Appendix 6 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51). This is because in simulations without hybrid collapse,

1 16 (loci) 0.001 0.0016 0.01

1–18

6

−2 −4

VT= 0.9 VS= 1.01

50000

80000

20000

Time (generations)

80000

Time (generations)

2

(5)

VT= 2.24 VS= 4.24

20000

50000

80000

Time (generations) (4)−(5)

Phylogenetic trees

(2)−(3)

(4)

FIGURE 1. Species diversity and phylogenetic tree shape (branching tempo and balance ) as a function of time. Time is measured from the common ancestor’s introduction in the landscape. Simulations were replicated 50 times. Top row: Black lines give median values; 95% conﬁdence interval shown in gray. Dashed horizontal lines indicate = 0 and = 0. VT and VS give the average variance, respectively over time (corrected according to the median among simulation replicates) and among simulation replicates. Vertical gray lines and associated labels highlight ﬁve different stages of and variation across time. Bottom row: Examples of species trees corresponding to four stages of diversiﬁcation. Competition parameters: C = 0.4, K = 1.1. See Table 1 for other parameter values.

which nodes are concentrated close to the base (due to a slowdown in the rate of lineage accumulation toward the present), is negative (Pybus and Harvey 2000). The statistic has been used extensively in macroevolutionary studies (e.g., Ricklefs 2007; McPeek 2008; Phillimore and Price 2008; Rabosky and Lovette 2008b; Liow et al. 2010; Quental and Marshall 2010; Davies et al. 2011); we were thus able to compare our results to those of previous theoretical studies and empirical estimates. Although values, as well as their variability, were shown to decrease with clade size (Pybus and Harvey 2000; McPeek 2008; Boettiger et al. 2012), it should not affect the interpretation of our results, as the latter reveal on the contrary an increase in values with age and with species number at steady state (Figs. 1–3). We analysed the balance of phylogenetic trees by calculating the statistic, which is the value that maximizes the likelihood in the -splitting model of clade growth (Aldous 1996, 2001; Blum and François 2006). This statistic allows one to compare the topology of trees with different numbers of tips (Aldous 1996), and has been widely used to describe tree balance (Blum and François 2006; Pigot et al. 2010; Davies et al. 2011).

[18:29 7/5/2015 Sysbio-syv014.tex]

Yule trees (pure-birth process) have = 0, trees more unbalanced than Yule trees have < 0 and trees less unbalanced than Yule trees have > 0 (Yule 1925). Finally, for both statistics and , we computed a measure of variance over time, VT, and a measure of variance among simulation replicates, VS. VT is the mean among simulation replicates of the variance of the difference between the statistic and its median value (calculated on all simulation replicates at the same time step) taken over all time steps. When VT is high, in each simulation the statistic tends to exhibit high temporal variation around its median value among simulation replicates. VS is the mean over time of the variance of the statistic among simulation replicates.

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

(1)−(2)

50000

(1) (2) (3)

−2

(5)

−6

(4)

0

60 40 20 20000

D

C (1) (2) (3)

Balance, log(β +2)

B

−10

(5)

2

(4)

Branching tempo, γ

(1) (2) (3)

0

Number of species

A

80

SYSTEMATIC BIOLOGY

Testing the Effect of Competition and Interaction with Abiotic Factors Default parameter values are given in Table 1. They were chosen to enable diversiﬁcation and span a broad range of possible effects on phylogenetic tree shapes, while keeping computational time under control.

Page: 6

1–18

2015

7

GASCUEL ET AL.—ECOLOGY, LANDSCAPE DYNAMICS AND SPECIES TREE SHAPE

Branching tempo, γ

Balance, log(β +2)

−2

*

−10

−4

20

* ●

2 2

−2 −6

20

−2 −4

*

2 −2 −10

−4

*

2 2

−2 −10

−4

−2

*

2 −2 −10

−4

*

2

20000

50000

80000

Time (generations)

*

*

−10

−4

−6

−2

2

●

* ●

−2

*

0

60 40 20 0

D Local catastrophes

*

●

*

−6

2

●

* ●

−2

*

0

60 40 20 0

●

*

−6

0

60 40 20 0

●

80

C2 Long isolation time

*

* ●

80

C1 Short isolation time

*

●

*

−6

2

●

* ●

−2

0

60 40 20 0

*

80

●

*

−10

0

60 40

●

80

0

B1 Slow dynamics B2 Fast dynamics

*

* ●

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

* ●

*

−6

2 0

●

* ●

−2

40

*

80

0

A Default

60

2

80

Number of species

20000

50000

80000

Time (generations)

20000

50000

80000

Time (generations)

FIGURE 2. Effect of the scaled width of competition and interaction with abiotic factors on the evolution of species diversity, branching tempo and balance of phylogenetic trees as a function of time. Time is measured from the common ancestor’s introduction in the landscape. Simulations were replicated 50 times. Black lines give the median (with 95% conﬁdence interval in dark gray) under wide competition (C /K = 0.75); white lines give the median (with 95% conﬁdence interval in light gray) under narrow competition (C /K = 0.25). Stars indicate statistically signiﬁcant differences (p-value < 0.05) between wide and narrow competition (stars near black lines for differences in median values, t-test; stars on top of dark gray areas, for differences in variance, F-test), either at the beginning or at the end of the simulations (respectively ﬁrst and last 20,000 generations). Dashed horizontal lines indicate = 0 and = 0. Rows B1 and B2 test for the effect of the pace of landscape dynamics (default A ≈ 4.76.10−5 , B1 ≈ 9.9.10−6 , and B2 ≈ 9.8.10−5 ). Rows C1 and C2 test for the effect of time in isolation (default A =≈ 0.95, C1 ≈ 0.52, and C2 ≈ 0.99). Row D tests for the effect of local catastrophes (default A = 0 and D = 2.10−4 ). See Table 1 for other parameter values.

[18:29 7/5/2015 Sysbio-syv014.tex]

Page: 7

1–18

8

SYSTEMATIC BIOLOGY Branching tempo, γ

Balance, log(β +2)

* −6 −10

−4

20 ●

2 2

−2

*

−10

−4

2 −2 −6 −10

2 0 −2

2 2

−2 −10

−4

−2

*

2 ●

*

−6

−2

* −10

−4

●

* ●

−2

2 0

60 40 20 0

*

*

−6 −10

−4

●

−2

2 0

●

*

●

−2

60 40 0

20

*

*

2

*

20000

50000

80000

Time (generations)

−6 −10

−4

●

*

●

−2

2

●

−2

*

0

60 40

*

20 0

E Heterogeneous landscape dynamics

80

D Local catastrophes

2

80

C2 Long isolation time

*

*

−6

0

60 40 20 0

●

* ●

80

C1 Short isolation time

*

*

*

−4

60 40 20

●

* ●

80

0

B2 Fast dynamics

*

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

●

*

−6

●

−2

0

60 40 20 ●

* ●

80

0

B1 Slow dynamics

*

*

−2

2 0

●

* ●

−2

40

*

*

80

0

A Default

60

2

80

Number of species

20000

50000

80000

Time (generations)

20000

50000

80000

Time (generations)

FIGURE 3. Effect of spatial heterogeneity in the scaled width of competition and interaction with abiotic factors on the evolution of species diversity, branching tempo and balance of phylogenetic trees as a function of time. Time is measured from the common ancestor’s introduction in the landscape. Simulations were replicated 50 times. Black lines give the median (with 95% conﬁdence intervals in dark gray) with heterogeneity; white lines give the median (with 95% conﬁdence intervals in light gray) without heterogeneity. Stars indicate statistically signiﬁcant differences (p-value < 0.05) with versus without heterogeneity (stars near black lines for differences in median, t-test; stars on top of dark gray areas for differences in variance, F-test), either at the beginning or at the end of the simulations (respectively ﬁrst and last 20,000 generations). Dashed horizontal lines indicate = 0 and = 0. Rows B1 and B2 test for the effect of the pace of landscape dynamics (default A ≈ 4.76.10−5 , B1 ≈ 9.9.10−6 , and B2 ≈ 9.8.10−5 ). Rows C1 and C2 test for the effect of time in isolation (default A = ≈ 0.95, C1 ≈ 0.52, and C2 ≈ 0.99). Row D tests for the effect of local catastrophes (default A = 0 and D = 2.10−4 ). Row E shows results when spatial heterogeneity impacts landscape dynamics rather than competition width (heterogeneity intensity ≈ 0.46). See Table 1 for other parameter values.

[18:29 7/5/2015 Sysbio-syv014.tex]

Page: 8

1–18

2015

GASCUEL ET AL.—ECOLOGY, LANDSCAPE DYNAMICS AND SPECIES TREE SHAPE

[18:29 7/5/2015 Sysbio-syv014.tex]

appearance f , picked in {4.10−4 , 2.10−3 , 4.10−3 }, and a rate of disappearance c, picked in {2.2.10−5 , 5.10−5 , 8.10−5 }. We compared the resulting tree shapes to those generated without heterogeneity (all borders having f = 10−3 and c = 5.10−5 ). These values ensure the same average pace of landscape dynamics, 4.76.10−5 , and the same average time in isolation, 0.95, in both sets of simulations. They also enable us to keep a similar level of spatial heterogeneity in scaled competition width and in landscape dynamics. Indeed, the ratio var(paceh )/paced of the standard deviation of all the equiprobable values paceh of pace of landscape dynamics with heterogeneity over paced , the default pace without heterogeneity, is equal to the ratio var((K /C )h )/(K /C )d of the standard deviation of all the equiprobable values (K /C )h with heterogeneity over (K /C )d , the default value without heterogeneity. The model was coded in C language using the GNU Scientiﬁc Library (Galassi et al. 2009) for random number generation. Computation and analyses of phylogenetic trees were performed using R (R Development Core Team 2012) and the packages ape (Paradis et al. 2004) and apTreeshape (Bortolussi et al. 2006). Codes are available as a ready-to-use R function, available in the Dryad data repository doi:10.5061/dryad.3bp51.

RESULTS Phylogenetic Tree Shape Through Time - Example First, we describe some of the salient features of the phylogenetic trees that come out when the model is run for a speciﬁc set of parameter values (Fig. 1). Simulations were replicated over 105 generations. Figure 1A shows that species diversity increases smoothly up to a stationary phase, with only small variation among replicates. Rather unexpectedly, the shape of the phylogenetic trees shows more complicated patterns, both across time and among replicates (Fig. 1B,C). Trees of different age, that is, time elapsed between the root and the time at which the ‘leaves’ are observed, can differ markedly in branching tempo and imbalance. After an initial phase where is highly unpredictable and rapidly peaks then decreases (age class labeled (1) in Fig. 1B,C), trees tend to be stemmy and unbalanced (median < 0, median < 0; age class (2)). However, after this initial decrease, starts increasing, followed by from age class (3). In older trees, the topology stabilizes around a stationary state, which under this parameter set tends to be overbalanced (median > 0 from age class (4)), with variation among replicates decreasing over time. Meanwhile, continues to increase and trees become tippy (median > 0). In the class of oldest trees (age class (5)), the branching tempo of phylogenetic trees also reaches a stationary state. Five remarks can thus be made on the dynamics of diversiﬁcation and tree shape for this particular set of parameter values: (i) Combining the ecological

Page: 9

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

We studied how the shape of phylogenetic trees is affected by the scaled width of competition, C /K (Doebeli and Dieckmann 2003). The ratio of the standard deviation of the competition kernel, C , over the standard deviation of the maximum carrying capacity per site, K , measures the range of resources utilized by each species relatively to the total range of available resources in each site. It is inversely proportional to the number of available local niches. Under wide (scaled) competition (low K , high C ), resources are concentrated close to the ecological optimum while the intensity of competition declines slowly with phenotypic distance, which broadens the range of resources utilized by each species, resulting in fewer available local niches. In contrast, under narrow (scaled) competition (high K , low C ), resources are more widely distributed while the intensity of competition declines rapidly with phenotypic distance, which narrows the range of resources utilized by each species, resulting in a higher number of available local niches. We tested values of C /K in {0.25, 0.375, 0.5, 0.75}, using values of C in {0.3, 0.4, 0.5, 0.6} and K in {0.8, 1.0, 1.1, 1.2}. To investigate how the effect of the scaled width of competition on the shape of phylogenetic trees interacts with abiotic factors, ﬁrst we looked at the effect of the pace of landscape dynamics (cycles of fusion and fragmentation of adjacent sites), as measured by 1/(1/f + 1/c) (with f the rate of border appearance and c the rate of border disappearance). The default pace is ≈ 4.76.10−5 (f = 10−3 , c = 5.10−5 ; Table 1). Low pace was set to 9.9.10−6 (f = 10−3 , c = 10−5 ), and fast pace was set to 9.8.10−5 (f = 5.10−3 , c = 10−4 ). Second, we looked at the effect of the time in isolation, measured by the fraction of time that a border is closed on average, f /(c+f ). When the time in isolation is long, the species in adjacent geographical sites tend to spend more time in allopatry, and less time in sympatry. Default, low, and high values of isolation time are 0.95, 0.52, and 0.99, respectively. Third, we looked at the effect of the rate of local catastrophe, cr. We compared the nolocal catastrophe scenario (cr = 0) to a case where the local catastrophe rate is high enough to decrease the species richness at steady state, but low enough to allow diversiﬁcation (cr = 2.10−4 ). We further investigated the inﬂuence of spatial heterogeneity by assuming that each geographical site could have its own values of C and K . In this case, C and K were drawn randomly from {0.3, 0.6} and {0.8, 1.2}, respectively. This allowed us to keep the same average ratio K /C = 2.5 over all sites in both the with- and without-heterogeneity simulations (ensuring a similar number of available local niches in the landscape). Here again, we explored how this effect is modulated by abiotic factors (pace of landscape dynamics, time in isolation and local catastrophes), using the same parameter values as above. Finally, we addressed the effect of spatial heterogeneity in landscape dynamics by randomly assigning to each geographical border a rate of

9

1–18

10

SYSTEMATIC BIOLOGY

[18:29 7/5/2015 Sysbio-syv014.tex]

time the early imbalance and stemminess of the trees have been overwhelmed by subsequent diversiﬁcation due to landscape dynamics and by turnover due to competition and demographic stochasticity. Moreover, the proportion of species that originated from the early allopatric radiation has been strongly reduced, and so has their inﬂuence on the shape of phylogenetic trees. The above processes of diversiﬁcation are also responsible for the increase in variance of and the decrease in variance of among simulation replicates over time. First, young phylogenetic trees have lower variance in than older ones because allopatric speciation events generate less variability in branching tempo than diversiﬁcation due to landscape dynamics: The early allopatric radiation develops similarly in all simulation replicates, whereas landscape dynamics depend on random variation in the number and distribution of removal of geographical barriers. Second, young phylogenetic trees have higher variance in than older ones because they contain fewer species (online Appendix 7 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51, showing the relationship between species richness and variance in in phylogenetic trees simulated under the birth–death model). Effect of Ecological Processes - Competition Width Figure 2A illustrates the effect of the scaled width of competition on the shape of phylogenetic trees of different age. The scaled width of competition is measured by the ratio C /K of the width of the competition kernel over that of the carrying capacity function. It is negatively related to the number of local ecological niches in the stationary state of the metacommunity composition. Patterns obtained under narrow versus wide competition are shown in light versus dark gray. Increasing the scaled width of competition results in lowering both and values and keeping them negative over a wider range of tree age; and may even remain negative at all tree ages, as in the wide competition regime shown in Figure 2A. Under narrow competition, the number of local ecological niches is high, and thus many speciation events occur due to landscape dynamics, driving and toward positive values. This pattern is reinforced by the large number of species and their ecological similarity at steady state, which increase the actual effects of competition and raise species turnover. Under wide competition, the number of local ecological niches is low, diversity is more limited, and thus few speciation events occur due to landscape dynamics. As a consequence, and tend to retain the footprint of the initial allopatric radiation and to remain negative. This pattern is reinforced by the small number of species at steady state and their large ecological differentiation, which reduces the intensity of competition and hence reduces the species turnover.

Page: 10

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

processes of intra- and interspeciﬁc local competition with large-scale, stochastic landscape factors does not prevent a smooth buildup and stabilization of species diversity from emerging. (ii) Trees of different ages have different shapes; this is true even among older communities that have reached the diversity plateau. (iii) Deceleration in branching ( < 0), which conventional thinking associates with communities that are old enough to be ecologically saturated, is in fact only observed among the young trees (age class (2) in Fig. 1’s example, while species diversity is not even close to half the stationary level). (iv) Phylogenetic trees tend to shift from stemmy and unbalanced ( < 0, < 0) to tippy and over-balanced ( > 0, > 0). (v) Young phylogenetic trees exhibit more among-replicates variability in balance , but less among-replicates variability in branching tempo than old phylogenetic trees. We can relate these observations to the underlying biotic and abiotic factors of speciation and extinction. Initially, only one ancestral species is present in the landscape, but rapidly an allopatric speciation event occurs at each geographical site across the landscape, due to directional selection toward the ecological optimum (local adaptation). Because most of these speciation events directly originate by sequential differentiation from the ancestral species, young trees tend to be unbalanced. Within this class of young trees, tends to be negative because all allopatric speciation events occur within a short-time window (they occur independently and from the same process of directional selection), and much earlier than the next speciation events driven by landscape dynamics. Indeed, landscape dynamics cause speciation only if diverging populations are sufﬁciently ecologically different and sufﬁciently choosy for reinforcement at secondary contact to lead to stable sympatric coexistence. Choosiness only begins to evolve after populations have adapted to the local ecological optimum of their geographical site, due to disruptive selection, and independently of landscape dynamics. Therefore, landscape dynamics may cause speciation only after the initial allopatric radiation and evolution of assortative mating. As trees grow into age class (2), the species formed in allopatry in each site start to split due to landscape dynamics. Because all geographical borders in the landscape have the same rates of appearance and disappearance, these speciation events are homogeneously distributed among lineages, and starts to increase. The concurrent fact that continues to decrease indicates that the most recent speciation events occuring due to landscape dynamics do not occur fast enough to compensate for the increase in length of the branches formed by the allopatric speciations during age class (1). This compensation starts as trees grow into age class (3), when then starts to increase. This reﬂects an increase in species diversity due to landscape dynamics, associated with species turnover. When age class (4) is reached, tree balance plateaus. In contrast, continues to increase until age class (5), and only then stabilizes around a stationary value. At this point in

1–18

2015

GASCUEL ET AL.—ECOLOGY, LANDSCAPE DYNAMICS AND SPECIES TREE SHAPE

How Abiotic Factors Interact with Competition Pace of landscape dynamics.—The effect of competition on the shape of phylogenetic trees is not altered by the pace of landscape dynamics. Figure 2A,B1,B2 show that at either slow or fast pace, wide competition always decreases and by a similar amplitude. Under both wide and narrow competition, the pace of landscape dynamics essentially scales the effect of tree age on and : fast landscape dynamics promote a rapid increase in and following the allopatric radiation, and increase their stationary values. This is because landscape dynamics drive diversiﬁcation after the initial allopatric radiation (Aguilée et al. 2013); when fast, they obliterate the inﬂuence of the allopatric radiation (negative and negative ) on phylogenetic tree shape. However, the pace of landscape dynamics inﬂuences how the scaled width of competition affects the variance in and . First, slow landscape dynamics amplify the increase in the variance of under wide competition. This effect is likely related to the relationship between species richness and variance in (online Appendix 7 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51). Second, both slow and fast dynamics amplify the difference in the variance of under narrow and wide competition, but in different ways: slow landscape dynamics reduce the variance in under wide competition, whereas fast landscape dynamics reduce the variance in under narrow competition. This results from three different effects acting on the variance in : the importance of landscape dynamics compared to that of the initial allopatric radiation (which increases the variance in under narrow competition and fast dynamics), the stochasticity in the location of the branching events within small phylogenetic trees (which increases the variance in under wide competition and slow dynamics), and the stochasticity of landscape dynamics (which increases the variance in under slow dynamics).

[18:29 7/5/2015 Sysbio-syv014.tex]

Time in isolation.—Inspecting Fig. 2A,C1,C2 shows that isolation time strongly interacts with competition to shape phylogenetic trees. Short isolation time reduces the effect of variation in scaled competition width. Indeed, short isolation time causes populations to spend less time in allopatry, which limits their potential to differentiate, hence decreases the likelihood of speciation and increases the chance of hybridization and the risk of extinction in sympatry. As a consequence, both the initial allopatric radiation and the subsequent diversiﬁcation due to landscape dynamics occur more slowly, and species diversity at steady state is lower. This has two consequences for phylogenetic tree shape. First, species take more time to evolve assortative mating and thus many initial species can originate from one another sequentially (rather than from the ancestral species). This yields phylogenetic trees that are balanced under both narrow and wide competition. Second, the initial allopatric radiation is now temporally overlapped by the subsequent diversiﬁcation due to landscape dynamics, and the species diversity resulting from landscape dynamics is lower at steady state. Consequently, temporal variation in branching tempo is much less pronounced and phylogenetic trees are buffered against the strong inﬂuence of competition width. Overall, these results underscore the importance of geographical isolation for the ecological effects of competition to manifest themselves. Lasting isolation favors diversiﬁcation while allowing ecological interactions (competition) to exert a strong, negative inﬂuence over the number of species at stationary metacommunity composition. Local catastrophes.—Frequent local catastrophes decrease the effect of the scaled width of competition on both branching tempo and balance (Fig. 2A,D). When the rate of local catastrophes is high, branching tempo under wide competition becomes close to its value under narrow competition. Wide competition yields a large proportion of old lineages (formed during the early allopatric radiation), of which many are removed by local extinctions; as this happens, ecological niches open for new species to arise. The average age of the nodes in the phylogenies decreases as a consequence, which translates into an increasing branching tempo . These effects do not manifest under narrow competition because most species are recent, and thus most species impacted by local extinctions are also recent. Local catastrophes also increase tree balance under wide competition, reducing the difference with its value obtained under narrow competition (Fig. 2A,D). Indeed, local extinctions cause the loss of old lineages that originated from the same ancestral species, which in turn permits the evolution of more recent species more homogeneously spread among the different subclades. Moreover, the decrease in species richness (e.g., median value going from about 12 to 9 under wide competition, and from 62 to 51 under narrow competition, in Fig. 2A,D) may increase median , as

Page: 11

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

Furthermore, increasing the scaled width of competition affects the variability of branching tempo and balance . On the one hand, wide competition increases the variance in among simulation replicates (Fig. 2A, with VT=2.30, VS=3.96 under narrow competition, and VT=2.54, VS=8.24 under wide competition). Wide competition indeed limits species diversity, increasing variance in due to stochasticity (Fig. 1, online Appendix 7 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51). On the other hand, wide competition decreases the variance in , both through time and among replicates (Fig. 2A, with VT=1.21, VS=1.30 under narrow competition, and VT=0.46, VS=0.58 under wide competition). Under wide competition, most branching events result from the initial allopatric radiation, which generates less variability in than diversiﬁcation due to lansdcape dynamics (Fig. 1).

11

1–18

12

SYSTEMATIC BIOLOGY

expected under the birth–death model (online Appendix 7 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51). Eventually, local catastrophes increase variability. On the one hand, their random occurrence in time and space generates stochastic effects on tree shape. On the other hand, they cause an overall decline in species richness - an effect that in itself causes more variability in under the birth–death model (online Appendix 7 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51).

[18:29 7/5/2015 Sysbio-syv014.tex]

DISCUSSION We have reported on the shapes of phylogenetic trees generated by an individual-based model of adaptive radiation driven by ecological speciation in a dynamic landscape. Our model is likely to apply to a broad range of real diversiﬁcation histories. Macroevolution often proceeds in bouts of radiations, and landscape dynamics are likely to have played a role in many (Schluter 2000). Celebrated examples include the radiation of the cichlid ﬁshes in the Great African Lakes (Sturmbauer 1998; Schwarzer et al. 2012; Aguilée et al. 2013), groups of plants and insects in the Macaronesian islands (Rijsdijk et al. 2014), and certain clades of birds in Amazonian forests (Haffer 1969; Terborgh 1992; Haffer 1997; Hill and Hill 2001; Weir 2006; Sedano and Burns 2010). However, phylogenetic tree shape is expected to depend on key features of the diversiﬁcation process. Pigot et al. (2010) and Davies et al. (2011) showed that vicariance associated with high range expansion, or peripatric speciation producing many species of lower abundance, can generate unbalanced phylogenetic trees by causing strong asymmetry in population size at speciation. Pontarp et al. (2012) highlighted that when within-site disruptive selection occurs before between-site directional selection, phylogenetic trees are more balanced and the deceleration in branching is less pronounced. Similarly, frequent sympatric speciation (Dieckmann and Doebeli 1999), recurrent allopatric speciation following successive landscape fragmentation events (Pigot et al. 2010), or sexual selection (e.g., Darwin 1871; Barraclough et al. 1995; Mitra et al. 1996) may all inﬂuence phylogenetic tree shape. Dependence of Phylogenetic Tree Shape on Clade Age Our results show that the tree shape is expected to change as a clade ages. This temporal pattern of variation results from the combination of two different mechanisms of diversiﬁcation. First, species differentiate under directional selection as they adapt locally to different ecological optima. This is the classical scenario of adaptive radiation (Schluter 2000). Because these speciation events occur at fairly similar rates across lineages and before the evolution of assortative mating, and because most of the derived species originate from the same ancestral species, phylogenetic trees initially exhibit decelerating branching and imbalance ( < 0 and < 0). Second, these species may further differentiate due to cycles of directional selection and reinforcement caused by the dynamics of the abiotic environment. This process is germane to the “species pump” mechanism (Greenwood 1965; Rossiter 1995; Aguilée et al. 2013). The new speciation events are homogenously

Page: 12

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

Spatial Heterogeneity Overall, spatial heterogeneity in scaled width of competition (i.e., different competition widths among geographical sites) has little qualitative effect on the patterns of branching tempo and imbalance - even though quantitative differences tested in the oldest trees are all signiﬁcant (Fig. 3A). Although heterogeneity in competition width has virtually no effect on the dynamics of species richness, it tends to delay the stabilization of and, to a lower extent, . As a consequence, the value of may switch sign from positive to negative in old trees. These results cannot be explained by differences in the total number of available local niches in the landscape with and without heterogeneity because our simulations control for a constant K /C = 2.5 across the landscape in both scenarios. However, plotting the stationary values of as a function of the scaled width of competition C /K without heterogeneity helps understand this effect (online Appendix 8 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51): the rate of decrease in increases non linearly with competition width, which suggests that in the presence of spatial heterogeneity, may be driven downward by geographical sites in which competition is wider. The inﬂuence of spatial heterogeneity on the branching tempo of phylogenetic trees is ampliﬁed by fast landscape dynamics and long isolation time, but prevented by slow landscape dynamics, short isolation time, or the frequent occurrence of local catastrophes (Fig. 3A–D). The latter are abiotic conditions that indeed limit species diversiﬁcation or species diversity at steady state, reducing the inﬂuence of competition width on phylogenetic tree shape, and thus that of spatial heterogeneity in competition width. We also examined the effect of spatial heterogeneity in landscape dynamics (Fig. 3E), which also causes both and to decrease, but to a lower extent than spatial heterogeneity in scaled competition width. Species close to slowly changing geographical borders tend to diversify less rapidly than species close to borders which undergo fast dynamics, reducing the balance and branching tempo in phylogenetic trees. However, this effect on the process of diversiﬁcation is limited, and is

smaller than that of spatial heterogeneity in competition width, which determines the distribution of species in geographical sites at steady state.

1–18

2015

GASCUEL ET AL.—ECOLOGY, LANDSCAPE DYNAMICS AND SPECIES TREE SHAPE

Effects of Competition Width Similarly to McPeek (2008), we found that taking the ecological context of speciation and extinction into account is critical to understand patterns of phylogenetic

[18:29 7/5/2015 Sysbio-syv014.tex]

tree shape. A major difference of our model from McPeek (2008) is that we not only considered the intensity of ecological differentiation at speciation events, but also accounted for sustained competition among organisms. This allowed us to study not only the effect of ecological niche ﬁlling (as in McPeek 2008, but with different results - see previous paragraph), but also the effect of competition all along the metacommunity’s history. Our results show that under wide competition (wide competition kernel and narrow distribution of resources), the branching tempo and balance of phylogenetic trees are determined primarily by the initial allopatric radiation: in the long run, they keep a decelerating branching tempo and remain unbalanced ( < 0 and < 0). In contrast, under narrow competition (narrow competition kernel and wide distribution of resources), the shape of phylogenetic trees evolves rapidly: at steady state, trees are overbalanced, with positive branching tempo ( > 0 and > 0). Thus, no simple causal relationship should be expected between ecological niche ﬁlling and branching deceleration. These results generate some predictions similar to McPeek (2008). Clades with negative (young clades, and old clades evolving under wide competition) should harbor species that are ecologically well differentiated and exhibit little overlap of their geographic range. Clades with positive (old clades evolving under narrow competition) should contain species that are poorly differentiated, with substantial overlap in their geographic range. Another implication of these results is that negative values are not necessarily associated with zeroextinction rates. Some earlier theoretical studies argued that negative could only be generated under zeroextinction rates (e.g., Weir 2006; Rabosky and Lovette 2008b), a conclusion that was difﬁcult to reconcile with data from the fossil record (e.g., Stanley 1979; Raup 1991; Gilinsky 1994; Alroy 1996; Foote 2000). Our results yield a more nuanced prediction: negative occurs under biotic or abiotic conditions which are conducive to extinction rates at steady state that are low but non-zero (e.g. under wide competition - Fig. 1, short isolation time - Figs. 2C1 and 3C1, or spatial heterogeneity - Fig. 3A,E). Finally, our results uncover the important inﬂuence of the width of competition on tree balance, (Fig. 2A). By controlling the magnitude of diversiﬁcation due to landscape dynamics, competition width determines to what extent phylogenetic tree shape depends on early history (dominated by the allopatric radiation; wide competition) versus late history (driven by landscape dynamics; narrow competition) of diversiﬁcation.

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

distributed among subclades and occur sequentially, thus generating an increase in branching tempo and balance , potentially reaching positive values. In our model, the early phase of diversiﬁcation, and corresponding tree shape, is inﬂuenced by the site of introduction of the ancestral species. Introduction in a peripheral site causes border effect, resulting in less temporal variation in species tree shape, and a tendency for trees to be more balanced (online Appendix 5 available as Supplementary material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51). Whereas the clade age-dependence of may not have been, to our knowledge, recognized or even hypothesized in earlier studies, previous studies have already considered or built hypotheses around the variation in branching tempo with clade age (e.g., MacArthur 1972; Rosenzweig 1978; Pimm 1979; Walker and Valentine 1984; Nee et al. 1992; Lovette and Bermingham 1999; Weir 2006; McPeek 2008; Phillimore and Price 2008; Liow et al. 2010; Jønsson et al. 2012). The conventional view is that young and still diversifying clades should exhibit phylogenetic trees with positive branching tempo, whereas ecological saturation due to competition should result in negative branching tempo. Accordingly, McPeek (2008) found that slowtrait evolution opposing ecological niche ﬁlling results in positive values, whereas fast-trait evolution leading to ecological niche ﬁlling results in negative values. This view was challenged by Liow et al. (2010): their diversity-dependent lineage-based model with signiﬁcant extinction predicted negative in young clades and positive in older ones. This pattern is similar to our ﬁndings (Fig. 1B), and results from diversity dependence associated with species turnover, which erodes the footprint of initial radiation and increases the number of recent nodes in older clades. However, unlike Liow et al. (2010), our model uncovers a diverse array of quantitative effects of competition and abiotic factors on , which may eventually result in either positive or negative in old clades (Figs. 2 and 3). Our results conﬁrm that the speed of convergence toward the stationary metacommunity size differs from that toward the stationary branching tempo of phylogenetic trees, and from that toward their stationary balance (e.g., Figs. 2A,B2 and 3A,B2,C1, Liow et al. 2010). Thus, negative or positive values should not be taken as revealing an actual decrease or increase in the diversiﬁcation rate through time. They may be strongly related to early variation in the diversiﬁcation rate as well as to temporal changes in the rate of species turnover. This clearly illustrates (McPeek’s 2008) cautionary note that only quantiﬁes the overall variation (since the root) in the rate of lineage accumulation in reconstructed phylogenies.

13

Effects of spatial Heterogeneity Spatial heterogeneity, either in scaled width of competition or in rates of landscape dynamics, can decrease the balance of phylogenetic trees (Fig. 3A–D). Spatial heterogeneity in competition width generates differences among geographical sites in the

Page: 13

1–18

14

SYSTEMATIC BIOLOGY

Modulation by Abiotic Factors Our results show that the inﬂuence of ecology on macroevolutionary patterns is modulated by abiotic factors. Yet, different abiotic factors interact differently with ecological factors to inﬂuence phylogenetic tree shape. The effect of the scaled width of competition on tree shape is almost unaltered by the pace of landscape dynamics. Even though the latter controls the potential rate of diversiﬁcation after the initial allopatric radiation, and therefore inﬂuences tree shape, its effect is largely independent of that of competition width. Landscape dynamics present the metacommunity with opportunities for diversiﬁcation; ecological interactions (competition) determine the successes and failures of the diversiﬁcation process. Both the time spent in isolation and the occurrence of local catastrophic extinctions can almost obliterate the relationship between competition and phylogenetic tree shape. Short isolation time decreases the rate of lineage accumulation and limits species richness in metacommunity stationary state, thus hampering ecological niche ﬁlling. Local catastrophic extinctions remove species produced by ancient branching events, thus promoting new speciation events driven by

[18:29 7/5/2015 Sysbio-syv014.tex]

landscape dynamics and again limiting ecological niche ﬁlling. Both processes thus erase the footprint of early history and strengthen the recent history in phylogenetic trees, as anticipated by Haydon et al. (1994). Ecology and the Shape of Real Reconstructed Phylogenies To generate hypotheses on the mechanisms that could be responsible for branching slowdown (e.g., Nee et al. 1992; Zink and Slowinski 1995; Lovette and Bermingham 1999; Rüber and Zardoya 2005; Seehausen 2006; Weir 2006; McPeek 2008; Rabosky and Lovette 2008a; Jønsson et al. 2012) and imbalance (e.g., Guyer and Slowinski 1991, 1993; Mooers 1995; Purvis 1996; Mooers and Heard 1997; Blum and François 2006) in most real phylogenies (see also online Appendix 9 available as Supplementary Material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51), we examined the conditions, in terms of clade age, competition width, spatial heterogeneity, and abiotic conditions, under which our model predicts negative versus positive and (Figs. 2 and 3). It turns out that young phylogenetic trees all tend to exhibit < 0 and < 0, except those generated under short isolation time. In old clades, at least two mechanisms may keep and at negative values: wide competition preventing further diversiﬁcation, and spatial heterogeneity. These negative and values in old trees are, however, only obtained under low rates of local catastrophes, and a sufﬁciently long time in isolation. These results provide new insights on how to interpret and from real phylogenies. First, negative should not necessarily be interpreted as evidence for ecological niche ﬁlling or wide competition: is negative in young trees under narrow competition, and in trees of any age under any competition width provided the time in isolation is short. Similarly, positive is not necessarily related to early steps of diversiﬁcation or to narrow competition: is positive in trees of old clades evolving under wide competition if frequent local catastrophes occur. Second, negative does not necessarily reﬂect spatial heterogeneity: is negative without spatial heterogeneity if trees are young, or if competition is wide. Similarly, positive does not necessarily depend on the landscape homogeneity: is positive in the presence of spatial heterogeneity if the time in isolation is short. Pigot et al. (2010) argued that their geographic model of species’ birth and death predicts a pattern of phylogenetic imbalance correlating with decelerating branching. We also found a tendency toward a correlation between the signs of median and median (Figs. 2 and 3). This is due to three different processes affecting and in the same direction: the initial allopatric radiation (which decreases both and ), diversiﬁcation driven by landscape dynamics (which increases both and ), and spatial heterogeneity (which decreases both and ). However, as highlighted above, the correlation does not always hold: it varies with clade age and may be prevented by abiotic

Page: 14

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

number of species that can coexist at steady state. Spatial heterogeneity in landscape dynamics generates differences among subclades in the propensity to diversify. These effects, however, are fairly small, and do not generate the degree of imbalance seen in real reconstructed phylogenies (e.g., Guyer and Slowinski 1991; Heard 1992; Guyer and Slowinski 1993; Slowinski and Guyer 1993; Mooers 1995; Purvis 1996; Mooers and Heard 1997; Blum and François 2006). This may be due to the limits that our model sets on spatial heterogeneity (i.e., low variance of C /K and pace of border dynamics across the landscape), and space subdivision (only 9 distinct sites). In the wild, spatial heterogeneity is likely to be the rule at all spatial scales, both in ecological interactions and environmental dynamics. Therefore, we expect spatial heterogeneity in both to contribute signiﬁcantly to real phylogenetic tree imbalance, in line with previous claims (e.g., Simpson 1944; Rosenzweig 1995; De Queiroz 2002; Purvis et al. 2011). For example, in birds Owens et al. (1999) showed that clade size is associated with feeding habits and dispersal abilities: at higher taxonomic level, heterogeneity in feeding habits and dispersal abilities among clades is, therefore, expected to generate imbalanced phylogenetic trees. Such effects might be enhanced further by temporal heterogeneity in biotic and abiotic features (e.g., Grant 2003). Interestingly, spatial heterogeneity in competition and landscape dynamics also slows down branching, by preventing diversiﬁcation in sites that experience wide competition or have boundaries that rarely open. Ubiquitous natural heterogeneity of ecological systems is, therefore, expected to result in branching slowdown in real phylogenetic trees.

1–18

2015

GASCUEL ET AL.—ECOLOGY, LANDSCAPE DYNAMICS AND SPECIES TREE SHAPE

Generalizing the Use of Individual-Based Models of Species Diversiﬁcation Our model (in open access in the Dryad data repository at http://dx.doi.org/10.5061/dryad.3bp51) is easy to modify and run to test for the effect of other factors or processes on phylogenetic tree shape, such as the mating system (role of sexual selection), the landscape structure and dynamics, ecological optima that ﬂuctuate over time, key innovations, the structure of correlations among ecological traits, and interspeciﬁc ecological processes beyond symmetrical competition (e.g., asymmetric competition, consumer–resource interactions, mutualism, or parasitism). Our modeling framework indeed mostly applies to phylogenetic trees documented at the scale of genera, at which all species are more likely to share the same trophic level, and to interact through interspeciﬁc competition for resources. As more phylogenies become available at higher

[18:29 7/5/2015 Sysbio-syv014.tex]

taxonomic levels and among ecologically connected groups, such extensions will be particularly relevant. Individual-based models of species diversiﬁcation should allow investigation of a wide range of questions in macroevolutionary biology. They allow researchers to integrate biotic and abiotic processes across multiple scales, from the genetic determination of individuals’ traits to large-scale environmental changes, and to study the complete history of diversiﬁcation by speciation and extinction. This should help build the theory needed to ﬁll the gap between paleoenvironmental data, the fossil record, and phylogenetic trees reconstructed from extant molecular sequences.

SUPPLEMENTARY MATERIAL Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.3bp51.

ACKNOWLEDGMENTS We are grateful to the Institut National de la Recherche Agronomique (INRA) for giving us access to the Migale bioinformatic plateform and its high-performance computing cluster. We thank Guillaume Achaz, Florence Débarre, and Mike Barker, Katrina Dlugosch, Mike Sanderson and their laboratory members at the University of Arizona for discussions and suggestions, as well as Susanne S. Renner, Ally Phillimore and two anonymous reviewers, who greatly helped to improve a previous version of this manuscript.

FUNDING This work was funded by grants MANEGE (ANR-BLAN-0215), EVORANGE (ANR-09-PEXT011) from France Agence Nationale de la Recherche, France Laboratory of Excellence TULIP (ANR-10LABX-41 ; ANR-11-IDEX-0002-02), France Centre National de la Recherche Scientiﬁque (Mission pour l’Interdisciplinarité, «Pépinière Eco-Evo-Devo»), and the Partner University Fund (Ecole Normale Supérieure — University of Arizona partnership).

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

factors. Interestingly, short isolation time yields a unique combination of median < 0 and median > 0, which is not very sensitive to clade age, competition width, or spatial heterogeneity. Unlike real metacommunities, our simulations include species with very small population sizes (on the order of 100 individuals per species). In larger populations, the rate of speciation is expected to be higher during the phase of diversiﬁcation (i.e., before the steady state in species richness), due to faster phenotypic divergence. However, at steady state the species turnover should be lower, due to reduced demographic stochasticity. Across wider geographical space, in which the number of available local ecological niches is similar but populations are larger, we therefore expect fewer recent species, especially under narrow competition, and therefore both and might be lower than in our simulations. However, the interaction of competition width and landscape dynamics will still determine the dynamics of diversiﬁcation. Given the underlying mechanisms as we have understood them from this study, we do not expect qualitative change in the effects of competition width and landscape dynamics on phylogenetic tree shape. Eventually, our simulations show high variability in and through time and among replicates, in agreement with the analysis of real phylogenies (e.g., Purvis 1995; Pybus and Harvey 2000; Barraclough and Vogler 2002; Linder et al. 2003; Turgeon et al. 2005, and online Appendix 9 available as Supplementary Material on Dryad at http://dx.doi.org/10.5061/dryad.3bp51). This variability originates from the intrinsic stochasticity of the diversiﬁcation process, and is ampliﬁed by small species numbers. As a result, global patterns of variation in and may hint at important processes, but individual values for particular clades may be difﬁcult to interpret: two clades may exhibit quite different values of and even though they evolved under similar biotic and abiotic conditions.

15

REFERENCES Aguilée R., Claessen D., Lambert A. 2013. Adaptive radiation driven by the interplay of eco-evolutionary and landscape dynamics. Evolution 67:1291–1306. Aguilée R., Lambert A., Claessen D. 2011. Ecological speciation in dynamic landscapes. J. Evol. Biol. 24:2663–2677. Aldous D. 1996. Probability distributions on cladograms. In: Aldous D., Pemantle R., editors. Random discrete structures. New York: Springer. P. 1–18. Aldous D. 2001. Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today. Stat. Sci. 16:23–34. Alfaro M., Santini F., Brock C., Alamillo H., Dornburg A., Rabosky D., Carnevale G., Harmon L. 2009. Nine exceptional radiations plus

Page: 15

1–18

16

SYSTEMATIC BIOLOGY

[18:29 7/5/2015 Sysbio-syv014.tex]

Gavrilets S. 2003. Perspective: models of speciation: what have we learned in 40 years? Evolution 57:2197–2215. Gavrilets S., Losos J.B. 2009. Adaptive radiation: contrasting theory with data. Science 323:732–737. Gavrilets S., Vose A. 2005. Dynamic patterns of adaptive radiation. Proc. Natl. Acad. Sci. USA 102:18040–18045. Geritz S., Kisdi E., Meszéna G., Metz J. 1998. Evolutionary singular strategies and the adaptive growth and branching of the evolutionary tree. Evol. Ecol. 12:35–57. Gilinsky N. 1994. Volability and the Phanerozoic decline of background extinction intensity. Paleobiology 20:445–458. Gilman R.T., Behm J.E. 2011. Hybridization, species collapse, and species reemergence after disturbance to premating mechanisms of reproductive isolation. Evolution 65:2592–2605. Grant P.R. 1986. Ecology and evolution of darwin’s ﬁnches. Princeton (NJ): Princeton University Press. Grant P.R. 2003. What Darwin’s ﬁnches can teach us about the evolutionary origin and regulation of biodiversity. BioScience 53:965–975. Greenwood P.H. 1965. The cichlid ﬁshes of Lake Nabugabo, Uganda. Bull. Br. Mus. Nat. Hist. Zool. 25:139–242. Guyer C., Slowinski J.B. 1991. Comparison of observed phylogenetic topologies with null expectations among three monophyletic lineages. Evolution 45:340–350. Guyer C., Slowinski J.B. 1993. Adaptive radiation and the topology of large phylogenies. Evolution 47:253–263. Gyllenberg M., Meszéna G. 2005. On the impossibility of coexistence of inﬁnitely many strategies. J. Math. Biol. 50:133–160. Haffer J. 1969. Speciation in Amazonian forest birds. Science 165: 131–137. Haffer J. 1997. Alternative models of vertebrate speciation in Amazonia: an overview. Biodivers. Conserv. 6:451–476. Hanski I., Gaggiotti O.E. 2004. Ecology, genetics and Evolution of Metapopulations. London: Elsevier Academic Press. Haydon D., Radtkey R.R., Pianka E.R. 1994. Experimental biogeography: interactions between stochastic, historical, and ecological processes in a model archipelago. In: Ricklefs R.E., Schluter D., editors, Species diversity in ecological communities: historical and geographical perspectives. Chicago (IL): University of Chicago Press. P. 117–130. Heard S. 1992. Patterns in tree balance among cladistic, phenetic, and randomly generated phylogenetic trees. Evolution 46:1818–1826. Heard S., Hauser D.L. 1995. Key evolutionary innovations and their ecological mechanisms. Historical Biol. 10:151–173. Hill J.L., Hill R.A. 2001. Why are tropical rain forests so species rich? Classifying, reviewing and evaluating theories. Prog. Phys. Geog. 25:326–354. Jønsson K.A., Fabre P.-h., Fritz S.A., Etienne R.S., Ricklefs R.E., Jørgensen T.B. 2012. Ecological and evolutionary determinants for the adaptive radiation of the Madagascan vangas. Proc. Natl. Acad. Sci. USA 109:6620–6625. Keast A. 1961. Bird speciation on the Australian continent. Bull. Mus. Comp. Zool. 123:305–495. Kendall D. 1948. On some modes of population growth leading to R. A. Fisher’s logarithmic series distribution. Biometrika 35:6–15. Kopp M., Hermisson J. 2008. Competitive speciation and costs of choosiness. J. Evol. Biol. 21:1005–1023. Kozak K.H., Weisrock D.W., Larson A. 2006. Rapid lineage accumulation in a non-adaptive radiation: phylogenetic analysis of diversiﬁcation rates in eastern North American woodland salamanders (Plethodontidae: Plethodon). Proc. R. Soc. Lond. B. 273:539–546. Leimar O., Doebeli M., Dieckmann U. 2008. Evolution of phenotypic clusters through competition and local adaptation along an environmental gradient. Evolution 62:807–822. Linder H.P., Eldenäs P., Briggs B.G. 2003. Contrasting patterns of radiation in African and Australian Restionaceae. Evolution 57:2688–2702. Liow L.H., Quental T.B., Marshall C.R. 2010. When can decreasing diversiﬁcation rates be detected with molecular phylogenies and the fossil record? Syst. Biol. 59:646–659. Lovette I.J., Bermingham E. 1999. Explosive speciation in the New World Dendroica warblers. Proc. R. Soc. Lond. B. 266:1629–1636.

Page: 16

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

high turnover explain species diversity in jawed vertebrates. Proc. Natl. Acad. Sci. USA 106:13410–13414. Alroy J. 1996. Constant extinction, constrained diversiﬁcation, and uncoordinated stasis in North American mammals. Palaeogeogr. Palaeocl. 127:285–311. Bailey S.F., Dettman J.R., Rainey P.B., Kassen R. 2013. Competition both drives and impedes diversiﬁcation in a model adaptive radiation. Proc. R. Soc. Lond. B. 280:1–8. Barnosky A.D. 1999. Does Evolution Dance To The Red Queen Or The Court Jester? Journal of Vertebrate Paleontology 19(Supplement to no. 3):31A.1998 Barraclough T., Harvey P., Nee S. 1995. Sexual selection and taxonomic diversity in Passerine birds. Proc. R. Soc. Lond. B. 259:211–215. Barraclough T.G., Vogler A.P. 2002. Recent diversiﬁcation rates in North American tiger beetles estimated from a dated mtDNA phylogenetic tree. Mol. Biol. Evol. 19:1706–1716. Barton N.H., Polechová J. 2005. The limitations of adaptive dynamics as a model of evolution. J. Evol. Biol. 18:1186–1190. Blum M., François O. 2006. Which random processes describe the tree of life? A large-scale study of phylogenetic tree imbalance. Syst. Biol. 55:685–691. Boettiger C., Coop G., Ralph P. 2012. Is your phylogeny informative? Measuring the power of comparative methods. Evolution 66: 2240–2251. Bortolussi N., Durand E., Blum M., François O. 2006. apTreeshape: statistical analysis of phylogenetic tree shape. Bioinformatics 22:363–364. Claessen D., Andersson J., Persson L., de Roos A.M. 2007. Delayed evolutionary branching in small populations. Evol. Ecol. Res. 9: 51–69. Cox D., Lewis P. 1966. The statistical analysis of series of events. London: Methuen. Coyne J., Orr H. 2004. Speciation. Sunderland (MA): Sinauer Associates. Cracraft J. 1986. Origin and evolution of continental biotas: speciation and historical congruence within the Australian avifauna. Evolution 40:977–996. Darwin C. 1859. On the origin of species by means of natural selection. London: John Murray. Darwin C. 1871. The Descent of Man and Selection in Relation to Sex. London: John Murray. Davies T.J., Allen A.P., Borda-de Água L., Regetz J., Melián C.J. 2011. Neutral biodiversity theory can explain the imbalance of phylogenetic trees but not the tempo of their diversiﬁcation. Evolution 65:1841–1850. De Queiroz A. 2002. Contingent predictability in evolution: key traits and diversiﬁcation. Syst. Biol. 51:917–929. De Queiroz K. 2007. Species concepts and species delimitation. Syst. Biol. 56:879–886. Dieckmann U., Doebeli M. 1999. On the origin of species by sympatric speciation. Nature 400:354–357. Doebeli M., Dieckmann U. 2003. Speciation along environmental gradients. Nature 421:259–264. Doebeli M., Ispolatov I. 2010. Complexity and diversity. Science 328:494–497. Etienne R.S., Haegeman B., Stadler T., Aze T., Pearson P.N., Purvis A., Phillimore A.B. 2012. Diversity-dependence brings molecular phylogenies closer to agreement with the fossil record. Proc. R. Soc. Lond. B. 279:1300–1309. Ezard T.H.G., Aze T., Pearson P.N., Purvis A. 2011. Interplay between changing climate and species’ ecology drives macroevolutionary dynamics. Science 332:349–351. FitzJohn R.G. 2010. Quantitative traits and diversiﬁcation. Syst. Biol. 59:619–633. Foote M. 2000. Origination and extinction components of taxonomic diversity: paleozoic and post-Paleozoic dynamics. Paleobiology 26:578–605. Fritz S.A., Schnitzler J., Eronen J.T., Hof C., Böhning-Gaese K., Graham C.H. 2013. Diversity in time and space: wanted dead and alive. Trends Ecol. Evol. 28:509–516. Galassi M., Davies J., Theiler J., Gough B., Jungman G., Alken P., Booth M., Rossi F. 2009. GNU scientiﬁc library reference manual. 3rd ed. New York: Network Theory Limited.

1–18

2015

GASCUEL ET AL.—ECOLOGY, LANDSCAPE DYNAMICS AND SPECIES TREE SHAPE

[18:29 7/5/2015 Sysbio-syv014.tex]

Rabosky D.L., Lovette I.J. 2008a. Density-dependent diversiﬁcation in North American wood warblers. Proc. R. Soc. Lond. B. 275: 2363–2371. Rabosky D.L., Lovette I.J. 2008b. Explosive evolutionary radiations: decreasing speciation or increasing extinction through time? Evolution 62:1866–1875. Raup D. 1991. Extinction: bad genes or bad luck? New York: Norton. Raup D.M., Gould S.J., Schopf T.J.M., Simberloff D.S. 1973. Stochastic models of phylogeny and the evolution of diversity. J. Geol. 81: 525–542. Ricklefs R.E. 2007. Estimating diversiﬁcation rates from phylogenetic information. Trends Ecol. Evol. 22:601–610. Rijsdijk K.F., Hengl T., Norder S.J., Otto R., Emerson B.C., Ávila S.P., López H., van Loon E.E., Tjørve E., Fernández-Palacios J.M. 2014. Quantifying surface-area changes of volcanic islands driven by Pleistocene sea-level cycles: biogeographical implications for the Macaronesian archipelagos. J. Biogeogr. 41:1242–1254. Rosenzweig M.L. 1995. Species diversity in space and time. Cambridge (MA): Cambridge University Press. Rosenzweig M. 1978. Competitive speciation. Biol. J. Linn. Soc. 10: 275–289. Rossiter A. 1995. The cichlid ﬁsh assemblages of Lake Tanganyika: ecology, behaviour and evolution of its species ﬂocks. Adv. Ecol. Res. 26:187–252. Roughgarden J. 1972. Evolution of niche width. Am. Nat. 106:683–718. Rüber L., Zardoya R. 2005. Rapid cladogenesis in marine ﬁshes revisited. Evolution 59:1119–1127. Sanderson M.J., Donoghue M.J. 1996. Reconstructing shifts in diversiﬁcation rates on phylogenetic trees. Trends Ecol. Evol. 11: 15–20. Schluter D. 2000. The ecology of adaptive radiation. Oxford: Oxford University Press. Schluter D. 2001. Ecology and the origin of species. Trends Ecol. Evol. 16:372–380. Schneider K., Bürger R. 2006. Does competitive divergence occur if assortative mating is costly? J. Evol. Biol. 19:570–588. Schwarzer J., Swartz E.R., Vreven E., Snoeks J., Cotterill F.P.D., Misof B., Schliewen U.K. 2012. Repeated trans-watershed hybridization among haplochromine cichlids (Cichlidae) was triggered by Neogene landscape evolution. Proc. R. Soc. Lond. B. 279:4389–4398. Sedano R.E., Burns K.J. 2010. Are the Northern Andes a species pump for Neotropical birds? Phylogenetics and biogeography of a clade of Neotropical tanagers (Aves: Thraupini). J. Biogeogr. 37:325–343. Seehausen O. 2004. Hybridization and adaptive radiation. Trends Ecol. Evol. 19:198–207. Seehausen O. 2006. African cichlid ﬁsh: a model system in adaptive radiation research. Proc. R. Soc. Lond. B. 273:1987–1998. Simpson G. 1944. Tempo and mode in evolution. New York: Columbia University Press. Simpson G. 1953. The major features of evolution. New York: Columbia University Press. Slowinski J., Guyer C. 1993. Testing whether certain traits have caused ampliﬁed diversiﬁcation - an improved method based on a model of random speciation and extinction. Am. Nat. 142:1019–1024. Stadler T. 2011. Mammalian phylogeny reveals recent diversiﬁcation rate shifts. Proc. Natl. Acad. Sci. USA 108:6187–6192. Stanley S.M. 1979. Macroevolution: pattern and process. San Francisco (CA): Freeman. Sturmbauer C. 1998. Explosive speciation in cichlid ﬁshes of the African Great Lakes: a dynamic model of adaptive radiation. J. Fish. Biol. 53:A18–A36. Terborgh J. 1992. Diversity and the tropical rain forest. New York: Scientiﬁc American Library, W.H. Freeman. Turgeon J., Stoks R., Thum R.A., Brown J.M., McPeek M.A. 2005. Simultaneous Quaternary radiations of three damselﬂy clades across the Holarctic. Am. Nat. 165:E78–E107. Van Valen L. 1973. A new evolutionary law. Evol. Theor. 1:1–30. Vonlanthen P., Bittner D., Hudson A.G., Young K.A., Müller R., Lundsgaard-Hansen B., Roy D., Di Piazza S., Largiader C.R., Seehausen O. 2012. Eutrophication causes speciation reversal in whiteﬁsh adaptive radiations. Nature 482:357–362. Vukics A., Asbóth J., Meszéna G. 2003. Speciation in multidimensional evolutionary space. Phys. Rev. E. 68:041903.

Page: 17

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

MacArthur R. 1972. Geographical ecology: patterns in the distribution of species. Princeton (NJ): Princeton University Press. Maddison W.P., Midford P.E., Otto S.P. 2007. Estimating a binary character’s effect on speciation and extinction. Syst. Biol. 56:701–710. Magnuson-Ford K., Otto S.P. 2012. Linking the investigations of character evolution and species diversiﬁcation. Am. Nat. 180: 225–245. Mayr E. 1942. Systematics and the Origin of Species. New York: Columbia University Press. Mayr E. 1947. Ecological factors in speciation. Evolution 1:263–288. Mayr E. 1963. Animal species and evolution. Cambridge (MA): Harvard University Press. McInnes L., Baker W., Barraclough K., Dasmahapatra T.G., Goswami A., Harmon L. 2011. Integrating ecology into macroevolutionary research. Biol. Lett. 7:644–646. McPeek M. 2008. The ecological dynamics of clade diversiﬁcation and community assembly. Am. Nat. 172:E270–E284. Mitra S., Landel H., Pruett-Jones S. 1996. Species richness covaries with mating system in birds. Auk 113:544–551. Moen D., Morlon H. 2014. Why does diversiﬁcation slow down? Trends Ecol. Evol. 29:190–197. Mooers A., Heard S. 1997. Inferring evolutionary process from phylogenetic tree shape. Q. Rev. Biol. 72:31–54. Mooers A.O. 1995. Tree balance and tree completeness. Evolution 49:379–384. Morlon H. 2014. Phylogenetic approaches for studying diversiﬁcation. Ecol. Lett. 17:508–525. Morlon H., Parsons T.L., Plotkin J.B. 2011. Reconciling molecular phylogenies with the fossil record. Proc. Natl. Acad. Sci. USA 108:16327–16332. Nee S. 2006. Birth-Death models in macroevolution. Annu. Rev. Ecol. Evol. Syst. 37:1–17. Nee S., Mooers A.O., Harvey P.H. 1992. Tempo and mode of evolution revealed from molecular phylogenies. Proc. Natl. Acad. Sci. USA 89:8322–8326. Owens I.P.F., Bennett P.M., Harvey P.H. 1999. Species richness among birds: body size, life history, sexual selection or ecology? Proc. R. Soc. Lond. B. 266:933–939. Paradis E., Claude J., Strimmer K. 2004. APE: analyses of phylogenetics and evolution in R language. Bioinformatics 20:289–290. Phillimore A., Freckleton R., Orme C., Owens I. 2006. Ecology predicts large-scale patterns of phylogenetic diversiﬁcation in birds. Am. Nat. 168:220–229. Phillimore A.B., Price T.D. 2008. Density-dependent cladogenesis in birds. PLoS Biol. 6:e71. Pigot A.L., Phillimore A.B., Owens I.P.F., Orme C.D.L. 2010. The shape and temporal dynamics of phylogenetic trees arising from geographic speciation. Syst. Biol. 59:660–673. Pimm S. 1979. Sympatric speciation: a simulation model. Biol. J. Linn. Soc. 11:131–139. Pontarp M., Ripa J., Lundberg P. 2012. On the origin of phylogenetic structure in competitive metacommunities. Evol. Ecol. Res. 14: 269–284. Purvis A. 1995. A composite estimate of primate phylogeny. Philos. T. Roy. Soc. B. 348:405–421. Purvis A. 1996. Using interspecies phylogenies to test macroevolutionary hypotheses. In: Harvey P., Leigh Brown A., Maynard Smith J., Nee S., editors, New Uses for New Phylogenies. Oxford: Oxford University Press. P. 153–168. Purvis A., Fritz S.A., Rodríguez J., Harvey P.H., Grenyer R. 2011. The shape of mammalian phylogeny: patterns, processes and scales. Philos. T. Roy. Soc. B. 366:2462–2477. Pybus O.G., Harvey P.H. 2000. Testing macro-evolutionary models using incomplete molecular phylogenies. Proc. R. Soc. Lond. B. 267:2267–2272. Pyron R.A., Burbrink F.T. 2013. Phylogenetic estimates of speciation and extinction rates for testing ecological and evolutionary hypotheses. Trends Ecol. Evol. 28:729–736. Quental T.B., Marshall C.R. 2010. Diversity dynamics: molecular phylogenies need the fossil record. Trends Ecol. Evol. 25:434–441. R Development Core Team. 2012. R: a Language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing.

17

1–18

18

SYSTEMATIC BIOLOGY

Walker T.D., Valentine J.W. 1984. Equilibrium models of evolutionary species diversity and the number of empty niches. Am. Nat. 124: 887–899. Wang S., Chen A., Fang J., Pacala S.W. 2013. Speciation rates decline through time in individual-based models of speciation and extinction. Am. Nat. 182:E83–E93. Waxman D., Gavrilets S. 2005. 20 Questions on adaptive dynamics. J. Evol. Biol. 18:1139–1154.

Weir J.T. 2006. Divergent timing and patterns of species accumulation in lowland and highland neotropical birds. Evolution 60:842–855. Yule G. 1925. A mathematical theory of evolution, based on the conclusions of Dr. J.C. Willis, F.R.S. Philos. T. Roy. Soc. B. 213: 402–410. Zink R., Slowinski J. 1995. Evidence from molecular systematics for decreased avian diversiﬁcation in the Pleistocene Epoch. Proc. Natl. Acad. Sci. USA 92:5832–5835.

Downloaded from http://sysbio.oxfordjournals.org/ at INIST-CNRS on June 17, 2015

[18:29 7/5/2015 Sysbio-syv014.tex]

Page: 18

1–18

Version dated : 3 mars 2015 ECOLOGY, LANDSCAPE DYNAMICS AND SPECIES TREE SHAPE

How Ecology and Landscape Dynamics Shape Phylogenetic Trees 1

´gis Ferrie `re3,4,+ , Robin Aguile ´e5 , and Amaury Fanny Gascuel1,2,3,∗ , Re Lambert1,6,+ Center for Interdisciplinary Research in Biology, CNRS UMR 7241, Coll`ege de France, Paris, France ; 2 Sorbonne Universit´ ´ es, UPMC Univ Paris 06, CNRS UMR 7625, Laboratoire Ecologie et ´ Evolution, Paris, France ; 3 Institut de Biologie de l’Ecole Normale Sup´ ´ erieure, CNRS UMR 8197, Ecole Normale Sup´erieure, Paris, France ; 4 Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, USA ; 5 Laboratoire Evolution ´ et Diversit´e Biologique, CNRS UMR 5174, Universit´e Paul Sabatier, Toulouse, France ; 6 Sorbonne Universit´ es, UPMC Univ Paris 06, CNRS UMR 7599, Laboratoire Probabilit´es et Mod`eles Al´eatoires, Paris, France

∗

Corresponding author : Fanny Gascuel, Center for Interdisciplinary Research in Biology, Equipe Stochastic Models for the Inference of Life Evolution, Coll`ege de France, 11 place Marcelin Berthelot, 75005 Paris, France ; E-mail : [email protected].

APPENDIX

Appendix 1 : Sensitivity of phylogenetic tree shape to migration rate across connected sites In our main model, migration is uniform among connected geographical sites, i.e. connected sites are panmictic and species therein evolve in sympatry. We tested for the importance of this assumption by considering another version of the model where, within sets of connected sites, migration is more likely towards the sites where the offspring’s parents are located than towards other sites (parapatry). Figure A1 shows that this assumption has no major influence on the shape of phylogenetic trees.

Branching tempo, γ

Balance, log(β +2) 2 −2

●

−6

−4

VT= 0.86 VS= 0.99

VT= 2.79 VS= 6.98

* −6

−2

VT= 0.88 VS= 1

−4

20

−2

*

●

*

VT= 3.01 VS= 8.83

−2

●

20000

50000

80000

Time (generations)

−6 20000

50000

80000

Time (generations)

VT= 3.72 VS= 7.36

−10

VT= 0.9 VS= 0.97

−4

●

*

●

−2

*

0

40

60

2

2

●

*

●

−10

40

0

60

2

2

−10

2 0

●

−2

60 40 20 ●

20 0

Parapatry Low migration

80 0

Parapatry High migration

80 0

Sympatry

80

Number of species

20000

50000

80000

Time (generations)

Figure A1 – Effect of the probability of migration at birth among connected sites on the evolution of species diversity, branching tempo γ and balance β of phylogenetic trees through time. Time is measured from the introduction of the ancestral species in the landscape. Probability of migration at birth, p, decreases from top to bottom : p=1.0 (sympatry), p=0.8 (parapatry with high migration probability), and p=0.6 (parapatry with low migration probability). Data was computed over 50 simulation replicates. Black lines show the median values and grey areas give the 95% confidence interval. Dashed horizontal lines show γ=0 and β=0. VT and VS give the average variance, respectively over time (corrected according to the median among simulation replicates) and among simulation replicates. Stars show the statistical differences (p − value <0.05, t-test) with results with sympatry, either at the beginning or at the end (respectively first and last 20,000 generations) of the simulations. See Table 1 for other parameter values.

Appendix 2 : Sensitivity of phylogenetic tree shape to the approximation on the reproductive isolation between populations We verified that computing the probability of reproductive isolation between pairs of populations using their average trait values x1 , x2 and a, which reduces computation time significantly, does not affect our results. Indeed, exact probabilities of reproductive isolation between populations should be calculated as the average of the probability of reproductive isolation between each pair of individuals from both populations. We therefore compared the temporal patterns of variation in phylogenetic tree shape obtained using the approximation to those expected without approximation : Figure A2 shows that considering the mean trait values of populations does not strongly modify the shape of reconstructed species trees.

Branching tempo, γ

Balance, log(β +2)

−2

●

−6

−4

VT= 0.86 VS= 0.99

2 2

20000

50000

80000

Time (generations)

*

−6

−2

* VT= 0.85 VS= 0.96

20000

50000

80000

Time (generations)

VT= 2.54 VS= 4.98

−10

0

●

−4

●

* ●

−2

40

*

VT= 2.79 VS= 6.98

−10

2

●

−2

0

2

80 60 40 20 60

80 0

●

20 0

Without approximation

With approximation

Number of species

20000

50000

80000

Time (generations)

Figure A2 – Effect of approximating the probability of reproductive isolation between pairs of populations on the evolution of species diversity, branching tempo γ and balance β of phylogenetic trees through time. Time is measured from the introduction of the ancestral species in the landscape. In the top row, probabilities of reproductive isolation between pairs of populations were approximated using the mean trait values of populations, whereas in the bottom row they were calculated exactly as the mean over the probabilities of reproductive isolation between each pair of individuals of both populations. Data was computed over 50 simulation replicates. Black lines show the median values and grey areas give the 95% confidence intervals. Dashed horizontal lines show γ=0 and β=0. VT and VS give the average variance, respectively over time (corrected according to the median among simulation replicates) and among simulation replicates. Stars show the statistical differences (p − value <0.05, t-test) between both sets of results, either at the beginning or at the end (respectively first and last 20,000 generations) of the simulations. See Table 1 for other parameter values.

Appendix 3 : Sensitivity of phylogenetic tree shape to the time interval between each species identification We verified that the frequency at which we compute the existing species and identify their origin during the course of diversification was high enough to enable a precise determination of the relationships between species. Figure A3 shows that the time step of 100 generations that we chose leads to quite similar results as a time step of 10 generations : it is short enough to accurately capture real phylogenetic tree shapes. Branching tempo, γ

Balance, log(β +2)

−2

●

−6

−4

VT= 0.86 VS= 0.99

2 2

20000

50000

80000

Time (generations)

*

−2 −6

VT= 0.94 VS= 1.02

20000

50000

80000

Time (generations)

VT= 2 VS= 5.57

−10

0

*

●

−4

●

* ●

−2

*

VT= 2.79 VS= 6.98

−10

2

●

−2

0

2

80 60 40 20 40

60

80 0

●

20 0

Every 10 generations

Every 100 generations

Number of species

20000

50000

80000

Time (generations)

Figure A3 – Effect of the time interval between each species identification on the evolution of species diversity, branching tempo γ and balance β of phylogenetic trees through time. Time is measured from the introduction of the ancestral species in the landscape. In the top row, species were identified every 100 generations (default value), and in the bottom row, every 10 generations. Data was computed over 50 simulation replicates. Black lines show the median values and grey areas give the 95% confidence intervals. Dashed horizontal lines show γ=0 and β=0. VT and VS give the average variance, respectively over time (corrected according to the median among simulation replicates) and among simulation replicates. Stars show the statistical differences (p−value <0.05, t-test) between both sets of results, either at the beginning or at the end (respectively first and last 20,000 generations) of the simulations. See Table 1 for other parameter values.

Appendix 4 : Delineating species and following their evolutionary history To follow the evolutionary history of species and determine the phylogenetic tree generated during diversification, at each time step t (typically every 100 generations), we (i) first delineate the Nt existing species [St ]i∈{1,...,Nt } (Fig. A4), and then determine their (j) ancestry relatively to the species [St−1 ]j∈{1,...,Nt−1 } existing at the previous time step (Fig. A5).

a)

Individuals

b)

Populations

c)

Species

Figure A4 – To delineate species, we (a) extract the distribution of all individuals in the phenotypic space (x1 , x2 ), according to cells of width σµx1 × σµx2 ; (b) determine high density phenotypic cells, and group individuals within adjacent high density cells into populations ; (c) group populations that can interbreed (probability of reproductive isolation below the threshold thrri =99%) into species, using their average ecological and choosiness traits (x1 , x2 , a ; calculated over all their individuals) and the same function as for mating between individuals. In this second step (Fig. A5), our algorithm proceeds as follows. Because each (i) species St is a set of populations, we start by determining the identity of all these populations. This identity is the one of the most ecologically similar population at time t − 1 (based on mean ecological traits x1 and x2 ). If all the populations included into (i) (j) (i) species St have the same species identity St−1 , St takes this identity. However, if species (i) (j) St includes populations with different species identities [St−1 ]j∈{1,...,Nt−1 } , hybridization (i) has occurred, and we fix the identity of St as the one of the ecologically most similar of (j) these previous species [St−1 ]j∈{1,...,Nt−1 } . Those containing more individuals are more likely (i) to be closer to St and thus to transmit their identity (consistently with their higher (i) probability to transmit their genotype to the lineage descending from St ). Once identities (i) have been attributed to each species [St ]i∈{1,...,Nt } , we check that two or more different

a)

Time t-1

b)

Time t

c)

Phylogenetic tree

t-1

t

Time

(i)

Figure A5 – To determine the ancestry of the species [St ]i∈{1,...,Nt } delineated at time step t (as presented in Fig. A4), and thus the changes in the phylogenetic tree between time step t−1 and time step t (c), we compare the phenotypic traits of species and populations at time t (b) to those at time t − 1 (a). Descent relationships are determined by minimum Euclidian distances between average traits of entities at time t and at time t − 1. Speciation occurs if (i) (j) different species St descend from the same species St−1 ; in that case (e.g. the red and green (j) (i) (j) species at time t), the species St most similar to St−1 is St−1 (here the red species), whereas (i) (j) other species St (here the green species) are new ones, descending from St−1 . Hybridization (i) occurs if a species St includes populations descending from populations of several different (j) (i) species St−1 ; in that case St (e.g. the red species at time t), has evolved from the most (j) similar of its parental species St−1 (here the red species), partly due to hybridization, whereas (j) other species St−1 (here the yellow species) might go extinct. (i)

species do not have the same identity. If several species [St ]i∈{1,...,Nt } have the same (j) identity St−1 , speciation has occurred between time t − 1 and time t. We therefore update (j) their identity by calculating their phenotypic distance to St−1 : only the most similar (j) (j) species keeps the identity of St−1 , all other species being new ones descending from St−1 .

Appendix 5 : Competition width and phylogenetic tree shape in the presence of border effects In the default conditions, the ecological and assortative mating traits (x1 , x2 and a) of the ancestral species were determined by Lk , k ∈ {x1 , x2 , a}, diploid loci randomly chosen in a Gaussian distribution centered at zero. The ecological optima (x∗1 , x∗2 ) of the 9 geographical sites were taken in {-1,0,1}×{-1,0,1}. Thus, the ancestral species was rather well adapted to the central geographical site. We explored the influence of border effects acting during early diversification, by considering an ancestral species adapted to marginal ecological conditions. In that case, the ecological and assortative mating traits (x1 , x2 and a) of the ancestral species were determined by Lk , k ∈ {x1 , x2 , a}, diploid loci randomly chosen in a Gaussian distribution centered at one. The ancestral species was thus rather well adapted to the geographical site with ecological optimum (x∗1 =1, x∗2 =1). Results show that border effects influence the shape of phylogenetic trees, both initially and at steady state. First, γ exhibit less temporal variation, because the initial allopatric radiation is more gradual : all allopatric speciations do not occur so closely in time. The initial allopatric radiation is therefore not so decoupled from diversification due to landscape dynamics. Second, β is more variable but on average less negative initially, and therefore does not stay negative at steady state in clades undergoing wide competition. Indeed, during the initial allopatric radiation, species originating in geographical sites located on the opposite side (e.g., with ecological optimum x∗1 =-1, x∗2 =-1) of the ancestral species do not descend from the latter, but from recent species with intermediate phenotypes, thus generating more balanced trees.

Branching tempo, γ

Balance, log(β +2) 2

50000

80000

Time (generations)

−6

*

−2 −10

−6

2 0

20000

*

−4

20 ●

* ●

−2

60 40

2

−10

−4

*

●

*

*

−2

2 0

*

●

*

0

B Border effects

80 0

* ●

* ●

−2

60 40 20

A Default

80

Number of species

20000

50000

80000

Time (generations)

20000

50000

80000

Time (generations)

Figure A6 – How border effects alter the effect of the scaled width of competition on the evolution of species diversity, branching tempo γ and balance β of phylogenetic trees through time. Time is measured from the introduction of the ancestral species in the landscape. In row B, unlike in row A (default, from Fig.2.A), the early diversification was constrained by border effects (the ancestral species being adapted to marginal ecological conditions). Data was computed over 50 simulation replicates. Black lines surrounded by dark grey areas give the median and 95% confidence intervals under wide competition (σC /σK =0.75), while white lines and light grey areas give them under narrow competition (σC /σK =0.25). Stars show the statistical differences (p − value <0.05, t-test) between results under wide and narrow competition, either at the beginning or at the end (respectively first and last 20,000 generations) of the simulations. Dashed horizontal lines show γ=0 and β=0. See Table 1 for other parameter values.

Appendix 6 : Hybridization and influence of irreversible reproductive isolation on phylogenetic tree shape There is empirical evidence for the existence of hybrid collapse and hybrid speciation (e.g., Ungerer et al. 1998; Taylor et al. 2006; Behm et al. 2010; Elowsky et al. 2013; Hasselman et al. 2014; Kleindorfer et al. 2014), but mostly among species pairs. After some time, species are likely to evolve irreversible Dobzhansky-Muller genetic incompatibilities, which result from negative epistatic interactions between alleles that arose in independent genetic backgrounds, and prevent them from hybridizing or giving birth to a viable and fertile descent. In our model, we did not incorporate such long-term isolation mechanism. We indeed considered ”reproductive isolation” between species to be based only on the divergence of ecological phenotypes and evolution of assortative mating traits. Because traits evolve, reproductive isolation is thus reversible. Hybridization can happen between species even if the latter diverged a long time before. To explore the influence of hybridization on tree shape and species turnover at steady state (relatively to that of competition and demographic stochasticity) and to investigate tree dynamics uninfluenced by hybrid collapses, we developed a version of the model which accounts for irreversible reproductive isolation between species. We ran simulations under different competition width, σC /σK , because competition width have a strong effect on the proportion of simulations resulting in a hybrid collapse (the latter decreases with competition width), thus potential on the rate of species hybridization. We incorporated this isolation mechanism by preventing reproduction between individuals (and thus hybridization between species) with high number of loci harbouring incompatible alleles. New born individuals inherit alleles from their parents following a random independent segregation of parental loci, and may acquire new alleles following an infinite site model, with mutation probability µ. According to previous mathematical analyses (Orr 1995; Orr and Turelli 2001) and empirical results (e.g., Matute et al. 2010; Moyle and Nakazato 2010), the number of incompatible genetic loci between two individuals, J, increases rapidly (”snow-ball” effect) with their genetic distance d (i.e. the number of different mutations that they carry) : d J =q , 2 d where q is the probability of incompatibility between two loci, and is the number of 2 pairs of loci among the d different loci which separate both individuals. Following Gavrilets (2003), we computed the probability of reproductive isolation due to genetic incompatibilities between pairs of individuals, Q(C, J), as a normalized incomplete gamma function of the expected number of incompatible genetic loci, J, with parameter C quantifying the threshold in number of incompatible loci needed to generate reproductive isolation (Fig. A7) : R∞ Q(C, J) = 1/Γ(C) J tC−1 exp(−t)dt.

1.0 0.8 0.6 0.4

● ● ●

0.2

● ● ● ●

0.0

Probability of genetic incompatibility

●

0

10

20

30

Threshold C 0.1 1 3 5 10 20 40

40

Figure A7 – Probability of genetic incompatibility Q(C, J) as a (normalized incomplete Gamma function) of the genetic distance (number of different mutations) between pairs of individuals. Here, q = 0.1 µ = 8.10−5 , d varies between 0 and 50, and C varies between 0.1 and 40.

50

Genetic distance

We simulated diversification with the accumulation of genetic incompatibilities between lineages as explained above, using parameter values µ = 8.10−5 , q = 0.1 and C varying between 0.1 and 40 (Fig. A7). These parameter values ensure that fully differentiated species cannot hybridize : the parameter C determines the genetic distance over which species are reproductively isolated and, for the smaller value we used (C = 0.1), only two mutations are likely to generate reproductive isolation (Fig. A7). These parameter values also ensure that reproductive isolation following from genetic incompatibilities does not contribute to speciation. Previous studies indeed showed that speciation rates are often decoupled from the rate of evolution of reproductive isolation (Rabosky and Matute 2013), and we are here interested in diversification driven by ecological speciation. Results show that accounting for reproductive isolation due to the accumulation of genetic incompatibilities favors the persistence of diversity over long time scales, by preventing the generation of hybrid collapse (Fig. A8). However, accounting for genetic incompatibilities does not influence the effects of the scaled width of competition on phylogenetic tree shape (predicted in scenarios without hybrid collapse in the absence of genetic incompatibilities ; Fig. A9). Indeed, neither speciation nor extinction rates are influenced by the evolution of genetic incompatibilities (Fig. A10). On the one hand, even at low threshold C, genetic incompatibilities accumulate too slowly to increase speciation rate. On the other hand, if rather distant species hybridize, it always results in a hybrid collapse. Indeed, if distant species could hybridize without generating a hybrid collapse, one would expect that preventing hybridization between species that diverged not long ago (i.e. low threshold in the number of mutations needed to generate irreversible reproductive isolation C) would decrease the species turnover at steady state. As a result, we would expect an increase in speciation and extinction rate with the threshold C. On the contrary, we observe on Figure A10 that speciation and extinction rates are unaffected by the value

1.0 0.8 0.2

0.4

0.6

Figure A8 – Probability of stable persistence (no hybrid collapse in the 100,000 first generations) as a function of the threshold C of the number of incompatible loci needed to generate reproductive isolation. C=∞ corresponds to simulations without genetic incompatibility. µ = 8.10−5 , q = 0.1, σC = 0.3, and σK = 1.2. See Table 1 for other parameter values.

0.0

Probability of stable persistence

of C. This suggest that hybridization only occurs when diversification collapses into a hybrid swarm. Two conclusions follow. First, considering only scenarios in which hybrid collapse does not occur (as we did) is akin to implementing long-term irreversible reproductive isolation between species. Second, in the absence of genetic incompatibilities, species turnover at equilibrium is not due to hybridization between species, but only to competition and demographic stochasticity.

0.1 10 20 40 Threshold C

∞

Branching tempo, γ

Balance, log(β +2) 2

0

−6 2

−10

−4

20000

50000

80000

Time (generations)

−2 −10

●

*

−4

20

*

*

−6

2 0

●

* ●

−2

60 40

80 0

*

*

0

With genetic incompatibilities

B

* ●

*

−2

2

40

●

* ●

−2

60

*

20

Default

A

80

Number of species

20000

50000

80000

Time (generations)

20000

50000

80000

Time (generations)

Figure A9 – Influence of long-term irreversible reproductive isolation due to genetic incompatibilities on the effect of the scaled width of competition on the evolution of species diversity, branching tempo γ and balance β of phylogenetic trees through time. Time is measured from the introduction of the ancestral species in the landscape. In row B, unlike in row A (default, from Fig.2.A), genetic incompatibilities generate long-term reproductive isolation between species. Data was computed over 50 simulation replicates. Black lines surrounded by dark grey areas give the median and 95% confidence intervals under wide competition (σC /σK =0.75), while white lines and light grey areas give them under narrow competition (σC /σK =0.25). Stars show the statistical differences (p−value <0.05, t-test) between results under wide and narrow competition, either at the beginning or at the end (respectively first and last 20,000 generations) of the simulations. Dashed horizontal lines show γ=0 and β=0. These simulations were performed with µ = 8.10−5 , q = 0.1, and C = 0.1. See Table 1 for other parameter values.

40

Threshold C

∞

5e−04 3e−04 1e−04

Extinction rate, µ

6e−04 4e−04 2e−04 0e+00

Speciation rate, λ

0.1 10 20

0.1 10 20

40

∞

Threshold C

Figure A10 – (Per species) speciation (λ) and extinction (µ) rates at steady state as a function of the threshold C of the number of incompatible loci needed to generate reproductive isolation. Boxplots show the distribution of the median values during the last 10,000 generations of the simulations, calculated over 50 simulation replicates ; vertical boxes represent the first, second and third quartiles, and whiskers give the upper and lower values (of maximum 1.5 times the interquartile distance). C=∞ corresponds to simulations without genetic incompatibility. µ = 8.10−5 , q = 0.1, σC = 0.3, and σK = 1.2 (parameter values associated with high probability of hybrid collapse). See Table 1 for other parameter values.

Appendix 7 : Variability in tree balance β and clade size

−2

●● ●●●●● ●● ● ●●●●● ●● ●●●●●●● ●●● ●●●● ●●● ● ● ● ●●● ● ●●● ● ●●●● ● ● ● ●●●● ● ●●●●●● ●● ●● ●●● ●● ● ● ● ●● ●●●●● ●● ●● ● ● ●●●●●●● ●●●●●● ●● ●●●●●●●● ●●● ● ● ● ●● ● ●●●●● ●●● ● ●● ●● ● ● ● ● ● ●●● ●● ●● ●● ●● ●● ●● ●●●●●●●● ● ●●●●● ● ●● ●● ● ● ●● ● ● ●●● ● ●●●● ● ●●● ●● ● ●● ●●● ●●●●●● ●●● ●● ● ●● ● ●● ● ● ●●●●●● ●● ●● ● ●● ●● ●●●● ●●● ●● ●●● ●● ●● ●●●●● ●●● ● ● ●● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ●●● ●●●● ●● ●●●●●●●●●●● ●●●● ●● ● ●● ● ●●● ●●● ● ● ●● ●●●●● ●● ● ● ●●●● ● ●●●● ●●● ●● ●●●● ●●● ●● ●●● ●●● ● ●● ●● ●●● ● ● ●● ●● ●● ● ●●● ● ●● ●● ●●● ●●●●● ●●●●●●● ● ● ●●● ● ●● ● ● ●● ●●● ●● ●●● ●●● ●● ●● ●●● ●● ● ●●● ●● ●●●● ●● ●● ● ●● ● ● ●●● ●● ●●● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●● ●●●●●● ●●● ●●● ●● ●● ●● ●● ●●● ●●●●●●●●●● ● ●●● ●● ●● ● ●●●●● ●●● ●●●●●●● ● ●● ● ●● ●●● ●●●●●● ●● ●● ●●● ●●● ●● ●●● ● ● ● ●● ●●● ●●● ●● ●●● ●●● ●●● ● ●● ●●●● ● ●● ●● ●● ●● ● ● ●● ●●● ●●● ● ●● ● ●● ●● ● ●●● ●● ●●●●●●● ●● ● ● ●● ●● ●● ●●●● ●● ●● ●● ●● ● ●● ●● ● ●● ●● ●● ●● ● ● ● ●●● ● ●●● ● ●● ●● ●● ●● ●● ●● ● ●● ●● ● ● ● ●● ● ●● ● ●●●●● ● ●● ●● ● ●●● ●● ●● ● ● ●●●● ● ●● ● ●●● ●●● ● ●●●●● ● ●● ●● ●● ●● ●●● ●● ● ●● ●●● ● ●● ● ●●● ●● ● ●● ●● ●●● ●● ●●●● ●● ● ●● ●● ●●●● ● ●● ●● ●● ●● ● ●●● ● ●●● ●● ●●●● ●● ● ● ● ●● ●● ● ●●● ●● ● ● ●● ●●● ●●● ●●●● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ●●●●●● ●● ●●● ●● ●● ●● ●● ●● ● ●● ● ● ●●●● ●● ● ● ●● ●●●●● ● ●● ●● ●● ●●● ●● ● ● ● ● ●● ● ●●●● ●● ● ●● ● ●● ● ● ●●● ● ●

●

●●

●

−6 −10

Balance, log(β+2)

2

Figure A11 shows that the variability in the balance β of phylogenetic trees depends on clade size : under the birth-death process (Yule 1925; Kendall 1948; Raup et al. 1973), phylogenetic trees which contain few species exhibit higher variability in β values (Fig. A11). Moreover, very small clades are expected to be characterized by a positive median β value.

●●●●●●●

0

50

100

150

Number of species

Figure A11 – Effect of clade size (number of species) on variability in the balance β of phylogenetic trees generated under the birth-death process. The dashed horizontal line shows β=0, and the red line shows the median value of β over time (computed on intervals of 10 units in species numbers). Trees were simulated on 50,000 time units, with birth rate = 0.0002 and death rate = 0.00018, using the R package (Harmon et al. 2008).

Appendix 8 : Non-linear relationship between competition width and phylogenetic tree shape at steady state

● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

0.25

● ● ● ● ● ● ● ● ●

●

0

● ●

● ● ●

● ● ●

●

−10

−2

0

2

● ● ● ● ● ● ● ●

Balance, log(β +2)

2

● ● ● ● ●

−6 −4 −2

● ●

−4

Branching tempo, γ

4

Branching tempo γ and the balance β of phylogenetic trees vary non linearly with the scaled width of competition : γ and β strongly decrease under the largest competition width (Fig. A12). Therefore, under spatial heterogeneity in competition width, γ and β may be driven downward by geographical sites in which competition is wider (i.e. having high σC /σK ).

0.5

0.75

Competition width

0.25

0.5

0.75

Competition width

Figure A12 – Branching tempo γ and balance β of phylogenetic trees at steady state as a function of the scaled width of competition (i.e. σC /σK ). Boxplots show the distribution of the median γ and β values during the last 20,000 generations of the simulations, calculated over 50 simulation replicates ; vertical boxes represent the first, second and third quartiles, and whiskers give the upper and lower values (of maximum 1.5 times the interquartile distance). Data were obtained from simulations without spatial heterogeneity, with values of σC in {0.3, 0.6} and σK in {0.8, 1.2} (similar to those used with spatial heterogeneity in scaled competition width). See Table 1 for other parameter values.

Appendix 9 : Shape of real reconstructed species trees at the genus level Here we search for patterns in the shape of real phylogenetic trees at the genus level, at which all species are assumed to share the same trophic level, and competition is potentially significant. We extracted recently published phylogenies of genera and characterized their shape in terms of branching tempo γ and balance β. We used phylogenies published by Phillimore and Price (2008), McPeek (2008), Jetz et al. (2012) and Bininda-Emonds et al. (2007), accounting for different groups of living organisms. We excluded all phylogenies of genera which contained less than 5 species or had less than two thirds of their nodes resolved. This led us to use 621 phylogenetic trees. We found that, as for trees not limited to genera (Nee et al. 1992; Pybus and Harvey 2000; McPeek 2008; Guyer and Slowinski 1991, 1993; Mooers 1995; Blum and Fran¸cois 2006), most genus-level phylogenies have negative γ and negative β (Fig.A13). However, γ and β show large variation (Fig.A13). A substantial fraction of trees have negative γ and positive β (including many estimates of β at its upper limit, 10), consistent with a pattern of evolution under short isolation time, or recent diversification in the presence of border effects.

●● ●

2

●

● ●

●

0 −2

●● ● ● ●● ● ● ●● ● ●● ●●● ● ● ●●● ●● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ●● ●● ● ● ●●● ● ● ●● ● ● ● ●● ● ● ●●● ●● ●● ● ● ●● ● ● ● ● ●

●

●

●

●

● ●

●

● ●●

● ● ● ● ●●● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ●● ● ●●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●●●●● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ●● ● ●● ● ●● ●●● ● ● ● ● ●● ● ● ●●●●● ●● ● ● ●● ● ● ● ● ● ●● ● ●● ● ●● ● ● ● ●●● ● ●●● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ●● ● ●● ● ●● ● ●● ●●● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ●●● ●● ● ● ●● ●● ● ● ●● ● ●●● ● ● ● ●● ● ●●● ● ●● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ●● ● ●● ● ● ● ●● ● ● ● ● ●● ●●●● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ●

●

●

●

●

Figure A13 – Branching tempo γ and balance β for phylogenies of genera, compiled from published data. Vertical and horizontal dashed grey lines respectively show γ=0 and β=0. Colors are according to the taxonomic group and data set from which we extracted each phylogeny (Phillimore and Price 2008; McPeek 2008; Jetz et al. 2012; Bininda-Emonds et al. 2007).

●

● ●

●

●

Bird genera (Phillimore and Price − 2008)

● ●

−4

●

Arthropoda genera (McPeek − 2008)

● ●

Chordata genera (McPeek − 2008)

−6

● ●

Magnoliophyta genera (McPeek − 2008)

●

−8

●

Mollusca genera (McPeek − 2008)

● ●

−10

Balance, log(β+2)

●

Bird genera (Jetz et al. − 2012)

● ●

●

● ● ●●

●● ●● ●●●● ●● ● ● ● ●● ●● ●● ● ● ● ●●

Mammal genera (Bininda−Edmonds et al. − 2007)

−10

−5

0

Branching tempo, γ

5

*

References Behm, J. E., A. R. Ives, and J. W. Boughman. 2010. Breakdown in postmating isolation and the collapse of a species pair through hybridization. Am. Nat. 175 :11–26. Bininda-Emonds, O. R. P., M. Cardillo, K. E. Jones, R. D. E. MacPhee, R. M. D. Beck, R. Grenyer, S. A. Price, R. A. Vos, J. L. Gittleman, and A. Purvis. 2007. The delayed rise of present-day mammals. Nature 446 :507–512. Blum, M. and O. Fran¸cois. 2006. Which random processes describe the tree of life ? A large-scale study of phylogenetic tree imbalance. Syst. Biol. 55 :685–691. Elowsky, C. G., I. E. Jordon-Thaden, and R. B. Kaul. 2013. A morphological analysis of a hybrid swarm of native Ulmus rubra and introduced U. pumila (Ulmaceae) in southeastern Nebraska. Phytoneuron 44 :1–23. Gavrilets, S. 2003. Perspective : models of speciation : what have we learned in 40 years ? Evolution 57 :2197–2215. Guyer, C. and J. B. Slowinski. 1991. Comparison of observed phylogenetic topologies with null expectations among three monophyletic lineages. Evolution 45 :340–350. Guyer, C. and J. B. Slowinski. 1993. Adaptive radiation and the topology of large phylogenies. Evolution 47 :253–263. Harmon, L. J., J. T. Weir, C. D. Brock, R. E. Glor, and W. Challenger. 2008. GEIGER : investigating evolutionary radiations. Bioinformatics 24 :129–131. Hasselman, D. J., E. E. Argo, M. C. McBride, P. Bentzen, T. F. Schultz, A. A. Perez-Umphrey, and E. P. Palkovacs. 2014. Human disturbance causes the formation of a hybrid swarm between two naturally sympatric fish species. Mol. Ecol. 23 :1137–1152. Jetz, W., G. Thomas, J. Joy, K. Hartmann, and A. Mooers. 2012. The global diversity of birds in space and time. Nature 491 :444–448. Kendall, D. 1948. On some modes of population growth leading to R. A. Fisher’s logarithmic series distribution. Biometrika 35 :6–15. Kleindorfer, S., J. A. O’Connor, R. Y. Dudaniec, S. A. Myers, J. Robertson, and F. J. Sulloway. 2014. Species collapse via hybridization in Darwin’s tree finches. Am. Nat. 183 :325–341.

Matute, D. R., I. A. Butler, D. A. Turissini, and J. A. Coyne. 2010. A test of the snowball theory for the rate of evolution of hybrid incompatibilities. Science 329 :1518–1521. McPeek, M. 2008. The ecological dynamics of clade diversification and community assembly. Am. Nat. 172 :E270–E284. Mooers, A. O. 1995. Tree balance and tree completeness. Evolution 49 :379–384. Moyle, L. C. and T. Nakazato. 2010. Hybrid incompatibility ”snowballs” between Solanum species. Science 329 :1521–1523. Nee, S., A. O. Mooers, and P. H. Harvey. 1992. Tempo and mode of evolution revealed from molecular phylogenies. Proc. Natl. Acad. Sci. USA 89 :8322–8326. Orr, H. A. 1995. The population genetics of speciation : the evolution of hybrid incompatibilities. Genetics 139 :1805–1813. Orr, H. A. and M. Turelli. 2001. The evolution of postzygotic isolation : accumulating Dobzhansky-Muller incompatibilities. Evolution 55 :1085–1094. Phillimore, A. B. and T. D. Price. 2008. Density-dependent cladogenesis in birds. PLoS Biol. 6 :e71. Pybus, O. G. and P. H. Harvey. 2000. Testing macro-evolutionary models using incomplete molecular phylogenies. Proc. R. Soc. Lond. B. 267 :2267–2272. Rabosky, D. L. and D. R. Matute. 2013. Macroevolutionary speciation rates are decoupled from the evolution of intrinsic reproductive isolation in Drosophila and birds. Proc. Natl. Acad. Sci. USA 110 :15354–15359. Raup, D. M., S. J. Gould, T. J. M. Schopf, and D. S. Simberloff. 1973. Stochastic models of phylogeny and the evolution of diversity. J. Geol. 81 :525–542. Taylor, E. B., J. W. Boughman, M. Groenenboom, M. Sniatynski, D. Schluter, and J. L. Gow. 2006. Speciation in reverse : morphological and genetic evidence of the collapse of a three-spined stickleback (Gasterosteus aculeatus) species pair. Mol. Ecol. 15 :343–355. Ungerer, M. C., S. J. Baird, J. Pan, and L. H. Rieseberg. 1998. Rapid hybrid speciation in wild sunflowers. Proc. Natl. Acad. Sci. USA 95 :11757–11762. Yule, G. 1925. A mathematical theory of evolution, based on the conclusions of Dr. J. C. Willis, F.R.S. Philos. T. Roy. Soc. B. 213 :402–410.

Version dated: June 12, 2015 ECOLOGY, LANDSCAPE DYNAMICS AND SPECIES TREE SHAPE

How Ecology and Landscape Dynamics Shape Phylogenetic Trees 1

´gis Ferrie `re3,4,+ , Robin Aguile ´e5 , and Amaury Fanny Gascuel1,2,3,∗ , Re Lambert1,6,+ Center for Interdisciplinary Research in Biology, CNRS UMR 7241, Coll`ege de France, Paris, France; 2 Sorbonne Universit´ ´ es, UPMC Univ Paris 06, CNRS UMR 7625, Laboratoire Ecologie et ´ Evolution, Paris, France; 3 Institut de Biologie de l’Ecole Normale Sup´ ´ erieure, CNRS UMR 8197, Ecole Normale Sup´erieure, Paris, France; 4 Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, USA; 5 Laboratoire Evolution ´ et Diversit´e Biologique, CNRS UMR 5174, Universit´e Paul Sabatier, Toulouse, France; 6 Sorbonne Universit´ es, UPMC Univ Paris 06, CNRS UMR 7599, Laboratoire Probabilit´es et Mod`eles Al´eatoires, Paris, France

∗

Corresponding author: Fanny Gascuel, Center for Interdisciplinary Research in Biology, Equipe Stochastic Models for the Inference of Life Evolution, Coll`ege de France, 11 place Marcelin Berthelot, 75005 Paris, France; E-mail: [email protected].

MODEL VERSION 2: Refining the criterion for species delineation In Gascuel et al. (2015), our method to delineate species and follow their evolutionary history was based on the probabilities of reproductive isolation between individuals, calculated from their phenotypic and assortative mating traits. This method did not account for geography: populations living in different geographical sites were grouped into the same species if their individuals had similar phenotypes. In an additional version of our model (presented in the online Appendix 6 of Gascuel et al. 2015, available at http://dx.doi.org/10.5061/dryad.3bp51), we wanted to account for the additional assumption that if species have diverged for long time in separated geographical sites, the accumulation of genetic incompatibilities may prevent them from hybridizing, even if their phenotypic traits are identical. We therefore introduced a mechanism of long-term

irreversible isolation resulting from the accumulation of genetic incompatibilities. This mechanism prevented mating between individuals that diverged a long time before. In the new version on the model (Version 2, introduced herein), we also account for genetic incompatibilities to delineate species and follow their evolutionary history. When genetic incompatibilities influence mating between individuals but not species delineation (as in Appendix 6 in the Dryad repository 10.5061/dryad.3bp51), the model would assume that two genetically incompatible species (either in the same or in different sets of connected sites) that converge phenotypically would eventually hybridize, with one of them going extinct. It would also assume that if the species’ traits diverged again after some time, a new speciation event would occur on the phylogeny. But in fact neither hybridization leading to extinction nor speciation actually occurred. This may generate discrepancies between predicted reproductive isolation between populations and their actual occurrence. In conclusion, accounting for genetic incompatibilities only at the individual level might lead to overestimate species turnover, underestimate the number of species, and distort descent relationships between species. In Gascuel et al. (2015), Appendix 6 was designed to investigate the dynamics of species trees when hybrid collapses are prevented (by the accumulation of genetic incompatibilities at the individual level). With the new version of the model, we account for the effect of genetic incompatibilities on species delineation to revisit how hybridization influences species turnover at steady state. We further investigate how this refined criterion for species delineation impacts the structure and dynamics of species trees.

How the model can be modified to account for genetic incompatibilities in species delineation In this new version of the model, we delineate species and follow their evolutionary history during diversification as in Gascuel et al. (2015, with methodological details in the online Appendix 4, at http://dx.doi.org/10.5061/dryad.3bp51), but with the following refinements: - We extract the distribution of individuals in the phenotypic space (x1 , x2 ) and group individuals within adjacent high density cells into populations independently for each set of connected sites; - Populations are grouped into the same species if (i) they can interbreed (as defined Gascuel et al. 2015, if the probability of reproductive isolation is below the threshold thrri =99%), and (ii) if the two populations are also genetically compatible. We assess genetic incompatibility between populations using the same function of genetic distance as for individuals (see Appendix 6). The genetic distance between two populations is the number of different genetic loci harbouring incompatible alleles carried by their individuals, weighted by their frequency in each population. - We determine descent relationships between species at time t and species at time t − 1 based on the minimum number of genetic incompatibilities between their populations. We only use minimum Euclidian distances between their average phenotypic traits to

discriminate the line of descent if populations of a species at time t have an equal (minimum) number of genetic incompatibilities with several populations from different species at time t − 1. Figures 1 and 2 illustrate this new method to delineate species and follow their evolutionary history (update of Figures A4 and A5 in Appendix 4). a)

Individuals

b)

Populations

c)

Species

Set of connected sites A

Set of connected sites B

Figure 1: How to delineate species when accounting for genetic incompatibilities. In this example, we consider two sets of connected geographical sites (A and B), and represent individuals (black dots) in the phenotypic space. We first group individuals into populations by (a) extracting, for each set of connected sites, the distribution of individuals in the phenotypic space (x1 , x2 ), according to cells of width σµx1 × σµx2 ; and (b) determining high density phenotypic cells (in grey), and group individuals within adjacent high density cells into populations. Then, (c) we consider populations of all sets of connected sites (A, B,...) and group into species (same color) those that (i) can interbreed (probability of reproductive isolation below the threshold thrri =99%) and (ii) are not genetically incompatible. Here, the green populations in the sets of connected sites A and B can interbreed and are not genetically incompatible, the orange and pink populations could interbreed based on their phenotypic and assortative mating traits but are genetically incompatible, and the yellow population albeit not genetically incompatible with the green ones (due to recent divergence) cannot interbreed with them due to quite different phenotypic traits.

a) Set of connected sites A

Time t-1

b)

Time t

c)

Phylogenetic tree

t-1

t Set of connected sites B

Time

Figure 2: How to determine species ancestry when accounting for genetic in(i) compatibilities. To determine the ancestry of the species [St ]i∈{1,...,Nt } delineated at time step t (as presented in Fig. 1), and thus the changes in the phylogenetic tree between time step t − 1 and time step t (c), we compare the loci carrying genetic incompatibilities and the phenotypic traits of species and populations at time t (b) to those at time t − 1 (a). Descent relationships are determined by minimum number of genetic incompatibilities between entities at time t and at time t − 1 and, if the latter are equal for multiple species at time t − 1, by minimum Euclidian distances between average phenotypic traits. Speciation (i) (j) occurs if different species St descend from the same species St−1 ; in that case (e.g. the (i) (j) (j) red and yellow species at time t), the species St most similar to St−1 is St−1 (here the red (i) species), whereas other species St (here the yellow species) are new ones, descending from (j) (i) St−1 . Hybridization occurs if a species St includes populations descending from popula(j) (i) tions of several different species St−1 ; in that case St (e.g. the red species at time t), has (j) evolved from the most similar of its parental species St−1 (here the red species), partly due (j) to hybridization, whereas other species St−1 (here the orange species) might go extinct.

New results

●

Turnover at steady state 0e+00 1e−04 2e−04 3e−04

Species richness at steady state 35 40 45 50 55

In the following section, we first explore the influence of the rate of accumulation of genetic incompatibilities on species turnover at steady state. By doing so, we want to revisit our previous conclusions (presented in Appendix 6 of Gascuel et al. 2015) on the importance of hybridization in the species turnover in Gascuel et al. (2015). Then, we use these results to calibrate the rate of accumulation of genetic incompatibilities in order to prevent hybridization between distantly related species. This allows us to investigate how accounting for genetic incompatibilities to delineate species impacts the structure and dynamics of species trees, by comparing the shapes of species trees generated by this new version of model to those presented in Gascuel et al. (2015). Figure 3 shows that about fifty percent of the species turnover observed at steady state in Gascuel et al. (2015) was due to hybridization events that are predicted but do not actually occur (i.e. due to trait convergence of distantly related species in different sets of connected sites).

●

0e+00

2e−03 4e−03 Mutation rate

● ● ●●

● ● ● ● ● ● ●

0e+00

2e−03 4e−03 Mutation rate

Figure 3: Species richness and turnover at steady state as a function of the mutation rate of loci harbouring incompatible alleles. Simulations were replicated 5 times. Parameter values: C=0.1, q=0.1, µ = 1.10−3 , σC =0.3, and σK =0.8. See Table 1 in Gascuel et al. (2015) for other parameter values. To address the effect of these observed hybridization events on phylogenetic tree shape, we ran again all simulations from Gascuel et al. (2015) with this new version of the model. Below, we provide the equivalent of figures 1 to 3 in Gascuel et al. (2015, using the same parameter values). Parameter C quantifies the threshold in number of incompatible loci needed to generate reproductive isolation; it was set to C=0.1. Parameter q gives the probability of incompatibility between two genetic loci; it was set to q=0.1. The mutation rate for new alleles was set to µ = 1.10−3 . The probability of genetic incompatibilities was

VT= 1.65 VS= 1.85

(5)

2

(4)

−2

●

−6

Balance, log(β +2)

●

(1) (2) (3)

VT= 2.53 VS= 4.57

−10

0

2

(5)

−2

40

(4)

20000 50000 80000

20000 50000 80000

20000 50000 80000

Time (generations)

Time (generations)

Time (generations)

(1)−(2)

(2)−(3)

(3)−(4)

(4)−(5)

Phylogenetic trees

D

●

C (1) (2) (3)

−4

(5)

60

(4)

Branching tempo, γ

B (1) (2) (3)

20 0

Number of species

A

80

computed as explained in the online Appendix 6. The above parameter values ensure that genetic incompatibilities accumulate fast enough to prevent hybridization between distantly related species (reaching a plateau in species turnover, see Fig.3), but slow enough to remain computationally tractable and not be the cause of speciation (so that the drivers of speciation are kept ecological, as opposed to genetic).

Figure 4: Species diversity and phylogenetic tree shape (branching tempo γ and balance β) as a function of time. Time is measured from the common ancestor’s introduction in the landscape. Simulations were replicated 50 times. Top row: Black lines give median values; 95% confidence interval shown in grey. Dashed horizontal lines indicate γ=0 and β=0. VT and VS give the average variance, respectively over time (corrected according to the median among simulation replicates) and among simulation replicates. Vertical grey lines and associated labels highlight five different stages of γ and β variation across time. Bottom row: examples of species trees built from the end of each of the four stages of diversification. Competition parameters: σC =0.4, σK =1.1. See Table 1 in Gascuel et al. (2015) for other parameter values.

Figure 5: Effect of the scaled width of competition and interaction with abiotic factors on the evolution of species diversity, branching tempo γ and balance β of phylogenetic trees as a function of time. Time is measured from the common ancestor’s introduction in the landscape. Simulations were replicated 50 times. Black lines give the median (with 95% confidence interval in dark grey) under wide competition (σC /σK =0.75); white lines give the median (with 95% confidence interval in light grey) under narrow competition (σC /σK =0.25). Stars indicate statistically significant differences (p − value <0.05) between wide and narrow competition (stars near black lines for differences in median values, t-test; stars on top of dark grey areas, for differences in variance, F-test), either at the beginning or at the end of the simulations (respectively first and last 20,000 generations). Dashed horizontal lines indicate γ=0 and β=0. Rows B1 and B2 test for the effect of the pace of landscape dynamics (default A ≈ 4.76.10−5 , B1 ≈ 9.9.10−6 and B2 ≈ 9.8.10−5 ). Rows C1 and C2 test for the effect of time in isolation (default A = ≈ 0.95, C1 ≈ 0.52, and C2 ≈ 0.99). Row D tests for the effect of local catastrophes (default A = 0, and D = 2.10−4 ). See Table 1 in Gascuel et al. (2015) for other parameter values.

0

2

0

2

20

*

●

●

*

20000

Time (generations)

50000

80000 2

−6

−10

40

−2

20

−4

*

●

−2

●

80

2

60

0 ●

−2

2

0

20

−10

−4

−2

40

2

0 ●

−2

2

0

−4

*

*

●

*

20000

Time (generations)

50000

80000 −10

*

−6

80

60

20

−6

−2

40

80

2

60

0

−2

●

−6

2

0

0

*

−10

−2

80

60 ●

−4

40

B1 Slow dynamics

* 2

0

−10

−4

20

−6

−2

40

80

2

60

0

A Default

−2

●

−2

0

B2 Fast dynamics ●

−6

80

60

C1 Short isolation time ●

−10

−2

40

C2 Long isolation time

*

−4

20

D Local catastrophes

2

Number of species Branching tempo, γ Balance, log(β +2)

*

*

* ●

*

*

* ●

*

* ●

*

*

●

*

*

*

●

*

*

●

*

20000

Time (generations)

50000

80000

Figure 6: Effect of spatial heterogeneity in the scaled width of competition and interaction with abiotic factors on the evolution of species diversity, branching tempo γ and balance β of phylogenetic trees as a function of time. Time is measured from the common ancestor’s introduction in the landscape. Simulations were replicated 50 times. Black lines give the median (with 95% confidence intervals in dark grey) with heterogeneity; white lines give the median (with 95% confidence intervals in light grey) without heterogeneity. Stars indicate statistically significant differences (p − value <0.05) with versus without heterogeneity (stars near black lines for differences in median, t-test; stars on top of dark grey areas for differences in variance, F-test), either at the beginning or at the end of the simulations (respectively first and last 20,000 generations). Dashed horizontal lines indicate γ=0 and β=0. Rows B1 and B2 test for the effect of the pace of landscape dynamics (default A ≈ 4.76.10−5 , B1 ≈ 9.9.10−6 and B2 ≈ 9.8.10−5 ). Rows C1 and C2 test for the effect of time in isolation (default A = ≈ 0.95, C1 ≈ 0.52, and C2 ≈ 0.99). Row D tests for the effect of local catastrophes (default A = 0, and D = 2.10−4 ). Row E shows results when spatial heterogeneity impacts landscape dynamics rather than competition width (heterogeneity intensity ≈ 0.46). See Table 1 in Gascuel et al. (2015) for other parameter values.

0

●

0

60

2

0

20

*

●

20000

Time (generations)

50000

80000

−6

2

−10

−2

40

●

●

−2

2

80

2

0

−10

−4

20

−2

40

2

0 ●

−2

2

0

−10

−4

20

−6

−2

40

0

60

●

−2

2

80

2

0

−10

−4

20

−6

−2

40

80

2

60

0

−2

●

*

*

*

*

*

●

20000

Time (generations)

50000

80000

−6

80

60

2

0

−10

−4

20

40

−2

−6

80

2

60

0

A Default

*

−6

0

60

B1 Slow dynamics

−2

●

−2

●

−4

20 ●

−10

2

0

0 ●

−2

80

*

60

B2 Fast dynamics ●

−4

40

C1 Short isolation time

*

−2

2

80

C2 Long isolation time ●

−6

−2

40

D Local catastrophes

*

−10

−4

20

E Heterogeneous landscape dynamics

2

Number of species Branching tempo, γ Balance, log(β +2)

*

* ●

* ●

*

* * ●

*

* ●

*

* *

* ●

*

●

*

*

*

●

*

*

20000

Time (generations)

50000

80000

Conclusions Accounting for long-term irreversible reproductive isolation between species increases the species richness and decreases the species turnover at steady state, but does not change phylogenetic tree shape under any of all biotic or abiotic conditions we tested (compare above Figures 4, 5 and 6 to Figures 1, 2 and 3 in Gascuel et al. 2015). Therefore, all results on the role of ecology and landscape dynamics on phylogenetic tree imbalance and branching tempo from Gascuel et al. (2015) extend to this refined version of the model. Future work will use this new version of the model, in which species delineation and evolutionary history are determined by accounting for both populations’ evolutionary history and individuals’ reproductive isolation. This new version allows trait convergence to occur without hybridization between species evolving in disconnected geographical sites.

* References Gascuel, F., R. Ferriere, R. Aguilee, and A. Lambert. 2015. How ecology and landscape dynamics shape phylogenetic trees. Syst. Biol. 0:1–18.

Wildlife-And-Landscape-Ecology-Effects-Of-Pattern-And-Scale.pdf ...

Ph.D. position in landscape ecology and species distribution modeling ...

Position â Postdoc Position in the 'Ecology/Landscape of Fear ...