Imprints of Past Habitat Area Reduction on Extant Taxonomic, Functional, and Phylogenetic Composition

Introduction

Alternating periods of contraction and expansion of suitable environmental conditions should affect demographic and speciation-extinction dynamics in a given habitat over time. Hence, a fundamental goal of biogeography is to better understand how changes in suitable environmental conditions due to past climatic and geomorphologic history have shaped current biodiversity patterns (Bennett, 1990). Much emphasis has been put on the impact of dramatic habitat restriction, whereby species could only persist in reduced areas called refugia (e.g., Weiss and Ferrand, 2007). Most studies so far have focused on how substantial habitat restriction in glacial refugia have shaped the genetic diversity within organisms at present (Avise, 2000). Conversely, we still poorly understand how past restrictions in refugia have influenced extant patterns for other dimensions of diversity in multi-species assemblages, a.k.a., taxonomic, phylogenetic, and functional diversity.

Here, considering the diversity of organisms making a biota in specific environmental conditions (habitat) of a given region, we examine how habitat contraction in refugia and subsequent expansion can affect these three different facets of biodiversity. Firstly, taxonomic diversity, based on the number and relative abundances of species in a local assemblage, depends on species demography (Hubbell, 2001). Changes in the area of available habitat should affect stochastic demographic dynamics, yielding greater abundance variation among species and extinction risk in smaller refugia (Figure 1). Specifically, we expect that the influence of ecological drift be greater in smaller assemblages and thus shape the relative abundances of species (Vellend, 2010; Blanchard et al., 2020). Thus, we expect demographic fluctuations to happen quickly over the timescale of habitat changes, so that extant taxonomic diversity should be more sensitive to recent than to ancient refugia. Secondly, phylogenetic diversity represents how closely related species are in a habitat according to their evolutionary history. We hypothesize that a larger biota in a larger habitat can undergo more or less frequent speciation or extinction events, respectively, and thus greater net diversification (Figure 1). Because diversification occurs at a slower pace than demographic dynamics, we expect phylogenetic diversity at present to keep a signal of ancient refugia. A greater number of old lineages could persist larger ancient refugia and have descendants in current assemblages (Bose et al., 2019). Conversely if the refuge was particularly small relative to the current assemblage size and expansion has only occurred recently, the loss of many lineages in refugia will have caused the most ancient common ancestor to be more recent (Figure 1). Indeed, the probability of sharing a common ancestor should increase in small assemblages (Patterson, 2005). Furthermore, the variation of trait attributes among extant organisms depends on how trait values have evolved across taxa in the past (Pearse et al., 2019). Specifically, assuming that species sharing greater common history have more similar phenotypes, we can account for the tempo of trait evolution to predict the functional composition of assemblages having undergone a common demographic history that shaped long-term diversification events (Revell et al., 2008). The speed of trait evolution could indeed have changed and depended on when and how adaptive radiations happened and shaped trait diversification over time in a given habitat (Bartish et al., 2016; Harmon, 2019). Therefore, extant patterns of trait diversity should differ depending on the pace and the timing of trait evolution. The extent of trait variation in descendants should depend on both trait evolution rate and the size of the biota over time (Figure 1). Specifically, if greater trait evolution occurs during strong habitat restriction, trait diversity among descendant lineages would be reduced.

Figure 1. Expected patterns of taxonomic, phylogenetic, and trait diversity under different scenarios of habitat expansion following a period with restricted refugia. Yellow shaded areas represent the size of suitable habitat conditioning the number of individuals over time (represented as circles). The size of the assemblage instantaneously increases at a given time (parameter t_c) so as to represent habitat expansion (parameter A). Current patterns of diversity will depend on how restrictive past habitat sizes were relative to present (left vs. right panels) and how long assemblages have spent in refugia (top vs. bottom panels). In larger refugia lower extinction rates allow more species in older lineages to persist in extant assemblages, while in smaller refugia more species will be lost to the increased influence of drift. Functional diversity will depend on the pace and timing of trait evolution in relation to diversification events. Trait evolution can deviate from the ancestral trait values (vertical dashed lines), and can experience exponential deceleration (Early Burst—EB represented in green) or acceleration (Late Burst—LB represented in blue), the resulting patterns of functional diversity will depend of whether acceleration or deceleration of trait evolution happens in or out of refugia. For example, if greater trait evolution occurs during strong habitat restriction, trait diversity among descendant lineages would be reduced.

The taxonomic, phylogenetic and functional facets of biodiversity are thus influenced by the interplay of historical variation in habitat area and past evolutionary dynamics. These three facets thereby convey complementary insights into population dynamics, evolutionary diversification, and trait evolution. To understand how past contraction in refugia and subsequent expansion of a habitat can affect taxonomic, phylogenetic and functional diversity, we propose an eco-evolutionary model explicitly representing the basic processes at play. First, we devise a genealogical model of population dynamics depending on the temporal variation in habitat area (Figure 2A). The model reconstructs past co-ancestry in the genealogy of extant individuals, using coalescent theory (Wakeley, 2009). We represent habitat restriction in past refugia with two parameters: one representing the ending time of refugia with subsequent and instantaneous expansion until present, and the other the extent of habitat restriction in the refugia (ratio of refugial area compared to present). Second, we model evolutionary diversification and phylogenetic relationships based on the genealogical structure. We group extant individuals into distinct species under a model of point mutations and protracted speciation of monophyletic genetic groups (Manceau et al., 2015; Figures 2B,C). Third, we model alternative scenarios of trait evolutions in the phylogeny, yielding extant trait composition. To represent different historical trajectories in trait evolution, we define the rate of trait variation in descendants along the phylogeny as either constant (Brownian motion), greater in older time (early burst), or greater in recent time (late burst). Different paces of trait evolution can be considered to reflect evolutive radiations concurring or not with lineages diversification brought about by demographic expansions (Figure 2C).

Figure 2. Eco-evolutionary dynamics of individual organisms in a habitat that experienced reduction in the past. The model relies on three components. (A) The population component reconstructs the effect of habitat expansion on the past co-ancestry in the genealogy of extant individuals, using coalescent theory without assuming to which species they belong to. (B) The evolutionary component represents how divergent genotypes have made new species over time. It is modeled by first sprinkling mutation events over the simulated genealogy conditionally its branch lengths. Then, the phylogenetic relationships amongst individuals and their abundances are obtained by merging paraphyletic clades into single species depending on the topology of the genealogy. (C) The trait evolution component is modeled using Brownian motion (BM), Early Burst (EB) of Late Burst (LB) along the simulated phylogeny and represents how species niches may have changed over time and at which pace. We represent two extreme scenarios: small demographic expansion (relative to refugia) that occurred long ago (top) and important habitat expansion that occurred recently (bottom) and how these have led to different patterns of taxonomic, phylogenetic, and functional diversity.

Here, we design a virtual experiment to examine how increases in favorable habitats following refugia influence taxonomic, phylogenetic and functional diversity. To this end, we explicitly simulate the demography of assemblages, the emergence of new species, their phylogenetic relationships and associated trait values.

Materials and Methods

Simulating Phylogenies and Trait Values in Extant Assemblages

In our model, fluctuations in the relative species abundances in a given habitat are driven by neutral drift dynamics, depending on the global size of the assemblage (Hubbell, 2001). Regardless of which species these individuals belong to, the dynamics of the assemblage they form can be represented by coalescence, i.e., by tracing the shared co-ancestry of extant individuals backwards in time until a single common ancestor is found (Kingman, 1982). Only the ancestry of individuals observed at present is traced back, so that coalescent methods do not require simulating lineages with no extant descendants and are thereby much faster than their forward-in-time alternatives (Munoz et al., 2018). Thus, without assuming to which species each individual belongs, we simulated the genealogies of individuals in an assemblage with constant size, and modeled the effect of habitat expansion on coalescence events following a period with restricted refugia (Figure 2A; Hudson, 2002; Kelleher and Lohse, 2020). We modulated both the amplitude of habitat expansion (A) and the time elapsed since expansion began (t_c), so as to simulate genealogies covering a wide range of refugia age and habitat expansion scenarios (Table 1).

Table 1. Model parameters and their prior distributions used to simulate species assemblages, the phylogenetic relationships, and trait values.

We sprinkled speciation events over the simulated genealogies (Figure 2B). The probability that a mutation occurred on a specific branch, leading to a new variant individual, followed a Poisson distribution with parameter μ^∗B, where μ is the point mutation rate (Table 1) and B the length of the branch. Descendant lineages of mutated individuals formed new clades that, if paraphyletic, were merged to form species (Manceau et al., 2015), thus monophyletic lineages with distinct genotypes and older than two generations were considered a new species (Figure 2C). We derived the number of individuals descending from a speciation event which gave us both the relative abundances of species and their phylogenetic relations to one another. Note that the point mutation rate is fixed for all simulations and thus the time scale (in generations before present) should be considered relative to μ. Additionally, any effect of expansions following periods with restricted refugia will depend on how parameter μ and the date of expansion are related. Indeed, if too few diversification events should happen during later expansion phases, as they may if μ is particularly small, we expect there to be little if not any signature in extant functional and phylogenetic diversity.

We simulated alternative models of trait evolution on the simulated phylogenies (Figure 2C), using the mvMORPH package in R (Clavel et al., 2015). We considered three models of trait evolution from an ancestral value set to 0 (Table 1). In the first model, traits evolved from their ancestral state following a pure Brownian motion (BM), with a variance parameter σ set either to 0.1, 0.2, or 0.4 to represent varying speed of trait evolution. In the second model, trait evolution exponentially decelerated with time (β = -0.01) following an Early Burst (EB). The third model represented Late Burst (LB) trait evolution, with exponential acceleration (β = 0.001).

We simulated 1,000 phylogenies for every past change date (5,000 in total), in which the model parameter values varied uniformly within above-mentioned limits (Table 1). We summarized the amplitude of the demographic expansion by parameter A, which reflects how much smaller the refuge was in the past before the favorable habitat expanded to its current and fixed size. Then, trait evolution was simulated with 100 replicates along each phylogeny according to the 9 parameter combinations given in Table 1.

Characterizing Extant Taxonomic, Functional, and Phylogenetic Diversity

We investigated the influence of the parameters of past expansion (A, t_c) on three different facets of extant biodiversity. We characterized taxonomic diversity by means of (i) species richness, and (ii) species frequencies spectra (SFS). The SFS represented the number of species by octaves of abundance class, and we performed a principal components analysis (PCA) (Lê et al., 2008) to summarize the variation of SFS across simulations. Similarly to using allelic frequency spectra in population genetics to infer population history (Marth et al., 2004), we expected the SFS to retain the signature of past demographic fluctuations (here, periods of habitat restriction) that caused more important drift dynamics affecting species relative abundances in smaller assemblages.

We used complementary metrics of phylogenetic diversity to characterize the evolutionary history of extant species of the assemblage. We quantified (i) Phylogenetic Diversity (PD), (ii) Mean pairwise distance (MPD), and (iii) Mean Nearest Taxon Distance (MNTD) using the picante package (Kembel et al., 2010). MPD is expected to be more sensitive to ancient and MNTD to recent diversification events (Mazel et al., 2016). Additionally, we used the approach described in Lewitus and Morlon (2016) implemented in the RPANDA package (Morlon et al., 2016) to derive the ln-transformed principal eigenvalue (λ^∗) and skewness (ψ) of the Laplacian spectral density profiles of phylogenies. These summary statistics quantified phylogenetic expansion and stem-to-tip distribution of branching events, respectively. Taxonomic and phylogenetic metrics of diversity were related to parameters of demographic expansion (parameters A and t_c) using linear regression models. To evaluate the relative effect of species richness from those of the parameters used to simulated demographic expansions, we used linear models between (i) the residuals of a linear model relating phylogenetic metrics to species richness and (ii) the parameters of demographic expansion. Linear models were fitted using the generalized least squares (GLS) method from the nlme package in R. In the case of PD and the MRCA, allowing the variance to increase with species richness yielded better AIC values when compared to simple linear regressions.

Finally, functional diversity was characterized by the first two moments of the simulated trait distributions (mean and variance) (Fortunel et al., 2014) given each phylogeny and given each trait evolution model parameter combination σ and β.

All simulations were conducted using Python version 3.7.4 (Python Software Foundation) and subsequent analyses performed using the R (3.6.3) statistical platform (R Core Team, 2020).

Results

Taxonomic Diversity

Simulated species richness decreased when the refugia was more restricted compared to current assemblage size, and it did so even faster for more recent changes (Figure 3). We obtained maximum species richness for oldest demographic expansion (i.e., 3,000 generations ago). When the refugia was more ancient, the restrictedness of the refugia did no longer influence extant species richness.

Figure 3. Variation in simulated species richness depending on the amplitude of expansion (A) and the time (in generations) at which it occurred (t_c). Linear regressions are shown in black with the associated R² statistic in the upper-right side of each panel.

The first PCA axis captured 44% of SFS variation and opposed the classes of lowest abundance to classes of intermediate abundance. PCA scores along the first axis were strongly related to species richness (R² = 0.88), and likewise decreased with the amplitude of expansion (Figure 4D). Yet, depending on how ancient the end of refugia was, the amplitude of subsequent expansion did not affect the abundance classes equivalently. Positive PCA scores along the first axis indicated assemblages with many rare species (Figure 4A). Such assemblages were mostly simulated when expansions occurred more than 1,500 generations ago (Figure 4D). Negative PCA scores along the first axis indicated assemblages with greater frequency of common species which, conversely, happened mostly when expansions occurred more recently (less than 500 generations ago). The second axis captured 17.2% of SFS variance and was mostly carried by the largest abundance class, representing assemblages with 0, 1, or 2 hyperabundant species (Figures 4B,C). For changes predating 2,000 generations, the number of abundant species increased with the amplitude of expansion and was maximal for 1,500 generations (Figure 4E).

Figure 4. Principal component analysis (PCA) on the simulated species frequency spectra (SFS). (A) Correlation circle showing the contributions of SFS abundance classes to the two first PCA axes. (B) Individual PCA scores projected onto the two first dimensions of the PCA and clustered according to the number of species in the most abundant abundance class. (C) Contributions of SFS abundances classes to the first (top) and second (bottom) axes of the PCA. (D) Variation in PCA-scores along the first dimension depending on the amplitude of expansion (parameter A) and the time (in generations) at which expansion occurred (parameter t_c). Linear regressions are shown in black with the associated R² statistic in the upper-left side of each panel. (E) Variation in the number of points with, respectively, none, one, or two hyperabundant species depending on the amplitude of expansion (parameter A) and the time (in generations) at which expansion occurred (parameter t_c).

Phylogenetic Diversity

We found that phylogenetic diversity (λ^∗, Faith’s PD, MPD, and the age of the MRCA) decreased overall with the amplitude of expansion—i.e., when the refuge relative to present was smaller, and did so more rapidly for recent change (see Figures 5A–C, 6). The ln-transformed principal eigenvalue λ^∗ correlated positively to PD, MRCA, MPD, and species richness (Supplementary Figure 2), and thus larger λ^∗ values were indicative of lengthy older rooted phylogenies with many lineages and greater overall phylogenetic diversity (Lewitus and Morlon, 2016). We found strong linear relationships between (i) the residuals of the regressions between phylogenetic metrics and species richness and (ii) the parameters of demographic expansions (all R² greater than 0.68; Supplementary Figure 1). This suggests that despite the influence of species richness on phylogenetic metrics of diversity, the latter still carried the signature of past demographic expansions. On the other hand, MNTD did not significantly vary along parameters of demographic change (Figure 5D), even when accounting for the influence of species richness. The density profile skewness ψ was always positive, which indicated stem-to-tip asymmetry favoring terminal branching events. ψ values increased with the amplitude of expansion (i.e, when the size of refugia decreased relative to present) and reached higher values when expansions occurred long ago (Figure 6B).

Figure 5. Variation in simulated metrics of phylogenetic diversity (A) Faith’s phylogenetic diversity (PD), (B) mean pairwise distance (MPD), (C) most recent common ancestor (MRCA), and (D) mean nearest taxon distance (MNTD) depending on the amplitude of expansion (A) and the time (in generations) at which it occurred (t_c). Linear regressions are shown in black with the associated R² statistic in the upper-right side of each panel.

Figure 6. Variation in spectral density profile summary statistics: (A) ln-transformed principal eigenvalue (λ*) and (B) skewness (ψ), of rerooted phylogenies (vis a vis oldest simulated MRCA) depending on the amplitude of expansion and the time (in generations) at which it occurred. Linear regressions are shown in black with the associated R² statistic in the upper-left side of each panel.

Functional Composition

Trait mean—as expected—averaged around the ancestral trait value and trait variance increased with σ (Figures 7, 8). Yet, depending on the trait evolution model (here, EB, BM, or LB, based on the sign of β) simulated trait values deviated more or less importantly from the ancestral value from one simulation to the next (larger IQR) and their variance decreased depending on the amplitude of expansion (Figure 7). More interestingly, the trait IQRs and variance increased depending on λ^∗ (Figure 8A), for the BM or LB trait evolution models. In the case of an EB trait evolution model, the IQR and trait variance did not depend on phylogenies’ spectral density metrics, but instead increased depending on σ. For the BM and LB trait evolution models, the IQR and trait variance increased both depending on σ and λ^∗, more so in the case of a LB. Simulated trait IQRs and variance also increased, with similar trends, for phylogenies with stem-to-tip asymmetry favoring terminal branch lengths.

Figure 7. Variation in median trait moments (mean and variance) as well as their interquartile ranges (IQR) depending on the amplitude of expansion (A). Trait mean and variance did not vary depending on the time at which expansion occurred (not shown). Species richness is fixed for all simulations (S = 31). Brownian Motion is represented in gray, while Early and Late Burst are shown, respectively, in green or blue.

Figure 8. Variation in median trait moments (mean and variance) as well as their interquartile ranges (IQR) along the ln-transformed principal eigenvalue (λ*) (A) and skewness (ψ) (B), summarizing the variation in simulated spectral density profiles according to the parameters used to simulate trait evolution (β in rows and σ in columns in each individual panel). The lines represent the median over 100 replications of trait evolution models and are dashed according to sigma values, while the shaded area represents the replicates’ interquartile ranges (IQR). Species richness is fixed for all simulations (S = 31). Brownian Motion is represented in grey, while Early and Late Burst are shown respectively in green or blue.

Discussion

By simulating the genealogies of individuals in an extant assemblage, for a wide range of scenarios of past restriction in refugia and subsequent expansion we found that the resulting patterns of taxonomic, phylogenetic and functional diversity in extant assemblages, reflect how past refugia have influenced population dynamics, evolutionary diversification, and trait evolution until present.

Relative Species Abundances Keep Track of (More) Recent Habitat Changes

When the size of favorable habitat in refugia was smaller relative to present size (i.e., ampler habitat expansion afterward), our results suggest that more species were lost over time due to greater extinction and ecological drift in reduced assemblages (Vellend, 2010). This phenomenon amplified as the time spent in refugia increased (i.e., when the subsequent habitat expansion was recent), and mostly affected rare species. Reciprocally, species richness and the number of rare species were highest when expansion occurred long enough ago that specific diversity could subsequently reach a new equilibrium, i.e., when it no longer depended on the amplitude of expansion. Hence, our results suggest that extant taxonomic diversity, such as species richness, is influenced by past demographic events having affected whole assemblages. More specifically, for recent habitat changes, taxonomic diversity kept track of the influence of restriction in more recent refugia. Motivated by analogous uses in population genetics, the Species Frequency Spectra (SFS), a multispecies extension of allele frequency spectra used in population genetics (Fisher, 1930), was especially useful to detect the influence of more recent demographic expansions. Because SFS provide more detailed and unabridged information regarding the variation in species abundance in an assemblage, these allowed us to detect the specific signature carried by specific classes of abundance, in particular that of rare species.

Phylogenetic Diversity Keeps Track of (More) Ancient Habitat Changes

Overall, phylogenetic diversity decreased as past assemblages experienced reduced-sized refugia and subsequent habitat expansion. Metrics of phylogenetic diversity (PD, MPD, λ^∗) reflected the higher extinction of lineages in smaller refugia. The effect of past habitat restriction was especially well-captured by metrics of phylogenetic diversity even when expansions were recent. MNTD failed to capture the signature of more recent expansions (Mazel et al., 2016). Nevertheless, the signature of assemblage-wide expansions was still perceptible in other phylogenetic diversity metrics (λ^∗, PD, MPD) despite very ancient and small changes in the past.

There has been a growing interest in how network science can provide novel ways to characterize the shape of phylogenetic trees (Chindelevitch et al., 2019). We found that the shape of a phylogeny, and more specifically its stem-to-tip asymmetry increased—indicating more frequent terminal branching events (Lewitus and Morlon, 2016)—as the size of the refugia decreased compared to present, even more so if expansions occurred long ago. With restricted past refugia, older lineages were lost by extinction, which amounted to very few deep branching events in extant phylogeny. Once the size of favorable habitat increased, so did diversification as new lineages had a greater chance of persisting into extant assemblages. With larger expansions, the greater imbalance between the loss of old lineages in refugia and the diversification of newer lineages after the expansion amounted to more terminal branching events. Thus stem-to-tip asymmetry increased. With recent expansions, there was not as much time for diversification to amount to as many new lineages (for a same rate of speciation for our simulations). Yet, demographic expansions affected both the number of terminal branching events and age of the MRCA. The simulated phylogenies should be all rooted to the oldest MRCA to capture the relative signature of expansions on stem-to-tip asymmetry (Pennell et al., 2012).

Hence, using different metrics of phylogenetic diversity and tree shape provide insights into how past habitat restrictions have impacted the diversification of species. Phylogenetic diversity provided a complementary insight to taxonomic diversity as taxonomic metrics mostly reflected the effect of past assemblage-wide demographic expansions while phylogenetic metrics reflect the effect of long-term speciation/extinction dynamics. Yet, species richness was influenced by global diversification, hence using phylogenetic diversity metrics which did not depend on species richness by capturing tree-shape such as the skewness of the spectral density profiles.

Phylogenetic diversity metrics carried the signature of assemblage-wide past demographic changes. Even if expansions occurred long ago, extant patterns of phylogenetic diversity still carried the signature of the loss of past lineages, unlike metrics of taxonomic diversity.

Functional Composition Reflects the Interplay Between Trait Evolution and Past Diversification in Phylogenies

Simulating trait evolution along phylogenies obtained for different expansion scenarios, trait values deviated from the ancestral value depending on the pace of trait evolution with regards to overall phylogenetic diversity. All replicates remained centered on the ancestral trait value. In the case of a Brownian (BM) or late burst (LB) trait evolution model, we showed that the extent by which traits deviated from the ancestral value depended on phylogenetic diversity (λ^∗), which in turn depended on parameters of past demographic changes. Our results show how patterns of functional composition change with regards to parameters of demographic expansion: past refugia dynamics have driven the levels of diversification of lineages and thereby the opportunity for traits to deviate from an ancestral trait value. In larger past refugia, increased diversification amounting to higher phylogenetic diversity allowed lineages to deviate more from the ancestral trait value. However, in the case of an early burst (EB), most trait evolution happened in smaller assemblages with fewer lineages, thus trait mean and variance remained small and unaffected by subsequent habitat expansion. Thus, the shape of a phylogeny interplays with trait evolution to shape extant functional composition (Tucker et al., 2018; Pearse et al., 2019). Specifically, phylogenetic asymmetry can be related to asymmetry in evolutionary models (Kirxpatrick and Slatkin, 1993; Blum and François, 2006). In this study, using a metric of stem-to-tip asymmetry (Lewitus and Morlon, 2016), we show that a substantial increase in trait variance can be due to a combined effect of increased terminal branching events and the acceleration of trait evolution (in the case of a LB model).

Yet additional metrics capturing other aspects of phylogenetic tree shape may have interesting applications to understand how regional trait diversity emerges (Pavoine et al., 2010; Tucker et al., 2017). For instance, if a phylogeny is particularly unbalanced, i.e., if branch lengths are unevenly distributed in the phylogenetic tree—then trait evolution may introduce multi-modal trait distributions. The emergence of functional outliers should then depend on the pace of trait evolution relative to the depth of diverging branching events (Violle et al., 2017). Reciprocally, we can expect trait distributions to reflect evolutive convergence if the distributions of branch-lengths are especially regular. Thus, specific signatures in other trait distribution moments such as the skewness or kurtosis (Gross et al., 2017) may result from the interplay between trait evolution and phylogenetic tree shape (Tucker et al., 2018).

Past environmental conditions may have constrained assemblages to refugia and removed lineages with ill-adapted sets of traits in a non-random manner (i.e., natural selection) and affect extant patterns of functional diversity (Svenning, 2003; Ordonez and Svenning, 2015). Here, we focus on the patterns of functional diversity emerging from the combination of lineage diversification (under the influence of expansion dynamics) and the pace of trait evolution. These types of approaches can reflect radiations which may or may not be adaptive (Czekanski-Moir and Rundell, 2019; Harmon, 2019). Our approach does not explicitly consider differential trait diversification between species, as may occur in the case of niche differentiation for instance (Dayan and Simberloff, 2005). A promising perspective would be to add a supplementary component to our model to represent more explicitly how natural selection affects the diversification of species niches depending on their adequacy to fluctuating environments (Ewing and Hermisson, 2010; Ovaskainen et al., 2011; Karhunen et al., 2013; Manceau et al., 2016).

Process-Based Eco-Evolutionary Modeling

In this study, we devised an original modeling framework to address how past refugia size and age can affect extant patterns of taxonomic, phylogenetic and functional diversity. Typically, phylogeographic studies focus on how current patterns of genetic diversity within species can retain the influence of past fragmentation and isolation of populations in refugia on intraspecific genetic diversity at present (Avise, 2000). Conversely, comparative phylogeographic approaches aim to grasp congruent or differing influence of refugia on several co-occurring taxa (Arbogast and Kenagy, 2001; Overcast et al., 2019, 2020; Swenson, 2019). Because smaller refugia entail greater extinction risk and fewer opportunities for speciation, the diversification of taxa over time should be shaped by changes in habitat availability. Here we confirm these hypotheses, by proposing an approach in which mechanisms of assemblage-wide dynamics of species abundances and trait values are explicitly simulated depending on changing habitat availability.

By directly representing the consequences of assemblage-wide environmental fluctuations over a wide range of parameter values relating to refugia size and age, we can address the broader question of whether any set of chosen metrics carry the signature of such processes. In fact, this approach echoes a more general challenge in ecology: can a snapshot of diversity patterns unambiguously reflect the signatures of many processes acting simultaneously at different spatial and temporal scales (Levin, 1992; Pearse et al., 2019). Relating extant patterns of diversity to the processes having hypothetically shaped them is an overarching goal in many fields of ecology and evolution. Examining how patterns relate to processes in a virtual setting can allow us to evaluate to what extent we can in turn recover parameter values of refugia dynamics from the observation of extant diversity patterns (Nielsen and Beaumont, 2009; Beaumont, 2010). Approximate Bayesian Computation (ABC) approaches for instance, repeatedly compare large numbers of simulated datasets to observed data, by means of summary statistics (Csilléry et al., 2010; Sunnåker et al., 2013). Thus, it is important that the simulated diversity metrics carry the signature of underlying processes, as we have shown here, because inference is directly based on the information withheld in the derived summary statistics (Grimm et al., 2010; Hartig et al., 2011). From here, we could now compare observed patterns of taxonomic, phylogenetic and functional diversity to those simulated using our process-based eco-evolutionary model in order to infer whether and how past expansion dynamics have shaped actual observations of biodiversity.

Perspectives—From a Regional Pool to Local Communities: Integrating the Influence of Assembly Processes on Local Community Composition

In recent years, there has been a renewed interest in better integrating evolutionary history and community ecology (Overcast et al., 2019, 2020; Swenson, 2019). Yet, studies of the assembly of communities often assume that the regional species pool (and associated phenotypes) represents a snapshot of an evolutionary heritage from which species assemble through local dynamics over recent and short enough time periods that diversification events can be ignored. Such approaches consider community assembly as a top-down process, whereby species present in a regional pool go through a series of filters (dispersal limitation, environmental filtering, biotic interactions) to assemble in local communities (Lortie et al., 2004). In this perspective, the species pool is stable, ignoring how the pool species composition is shaped by biogeographic history, colonization events or the tempo of trait evolution (Tucker et al., 2018; Swenson, 2019). However, as we have highlighted above, it is possible to explicitly represent how such processes together generate species composition and traits at a regional level.

Recent studies have underlined the critical influence of the diversity of the species pool in providing immigrants (Cornell and Harrison, 2014) and called for a more mechanistic modeling approach that directly represents key processes shaping community assembly (Denelle et al., 2019). The difficulty lies in the definition of the species pool as a fixed phenomenological object (Lessard et al., 2012a, b). Depending on how it is delineated, from species’ fundamental niches and their limited dispersal abilities in a broad-scale landscape (Karger et al., 2016) to finer-scale biotic interactions in shaping community assembly (Perronne et al., 2017), our perspective on the community level can vary. Here we show how local demographic processes influence the taxonomic, phylogenetic and functional competition of the species pool under the influence of long-term environmental fluctuations. Yet we could extend our framework to consider how these same large-scale biogeographic historical processes feed back into community assembly at a local scale. This could be done for instance by integrating a subsequent step to our modeling framework in which species can potentially be filtered based on their traits (Munoz et al., 2018).

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

EB and MJ performed the analyses. EB wrote the first draft of the manuscript. All authors contributed to design the study, revisions to the manuscript, and gave final approval for publication.

Funding

MJ was supported by the Institut de Recherche pour le Développement. This work has benefited from an “Investissements d’Avenir” grant managed by Agence Nationale de la Recherche (CEBA, ref. ANR-10-LABX- 25-01).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.634413/full#supplementary-material

References