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

Past environmental changes have shaped the evolutionary and ecological diversity of extant organisms. Specifically, climatic fluctuations have made environmental conditions alternatively common or rare over time. Accordingly, most taxa have undergone restriction of their distribution to local refugia during habitat contraction, from which they could expand when suitable habitat became more common. Assessing how past restrictions in refugia have shaped species distributions and genetic diversity has motivated much research in evolutionary biology and biogeography. But there is still lack of clear synthesis on whether and how the taxonomic, functional and phylogenetic composition of extant multispecies assemblages retains the imprint of past restriction in refugia. We devised an original eco-evolutionary model to investigate the temporal dynamics of a regional species pool inhabiting a given habitat today, and which have experienced habitat reduction in the past. The model includes three components: (i) a demographic component driving stochastic changes in population sizes and extinctions due to habitat availability, (ii) a mutation and speciation component representing how divergent genotypes emerge and define new species over time, and (iii) a trait evolution component representing how trait values have changed across descendants over time. We used this model to simulate dynamics of multispecies assemblages that occupied a restricted refugia in the past and could expand their distribution subsequently. We characterized the past restriction in refugia in terms of two parameters representing the ending time of past refugia, and the extent of habitat restriction in the refugia. We characterized extant patterns of taxonomic, functional and phylogenetic diversity depending on these parameters. We found that extant relative abundances reflect the lasting influence of more recent refugia on demographic dynamics, while phylogenetic composition reflects the influence of more ancient habitat change. Extant functional diversity depends on the interplay between diversification dynamics and trait evolution, offering new options to jointly infer current trait adaptation and past trait evolution dynamics.

Past environmental changes have shaped the evolutionary and ecological diversity of extant organisms. Specifically, climatic fluctuations have made environmental conditions alternatively common or rare over time. Accordingly, most taxa have undergone restriction of their distribution to local refugia during habitat contraction, from which they could expand when suitable habitat became more common. Assessing how past restrictions in refugia have shaped species distributions and genetic diversity has motivated much research in evolutionary biology and biogeography. But there is still lack of clear synthesis on whether and how the taxonomic, functional and phylogenetic composition of extant multispecies assemblages retains the imprint of past restriction in refugia. We devised an original eco-evolutionary model to investigate the temporal dynamics of a regional species pool inhabiting a given habitat today, and which have experienced habitat reduction in the past. The model includes three components: (i) a demographic component driving stochastic changes in population sizes and extinctions due to habitat availability, (ii) a mutation and speciation component representing how divergent genotypes emerge and define new species over time, and (iii) a trait evolution component representing how trait values have changed across descendants over time. We used this model to simulate dynamics of multispecies assemblages that occupied a restricted refugia in the past and could expand their distribution subsequently. We characterized the past restriction in refugia in terms of two parameters representing the ending time of past refugia, and the extent of habitat restriction in the refugia. We characterized extant patterns of taxonomic, functional and phylogenetic diversity depending on these parameters. We found that extant relative abundances reflect the lasting influence of more recent refugia on demographic dynamics, while phylogenetic composition reflects the influence of more ancient habitat change. Extant functional diversity depends on the interplay between diversification dynamics and trait evolution, offering new options to jointly infer current trait adaptation and past trait evolution dynamics.

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 longterm 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.
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 ecoevolutionary 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).
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.

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 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.
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 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.
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).
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 Frontiers in Ecology and Evolution | www.frontiersin.org TABLE 1 | Model parameters and their prior distributions used to simulate species assemblages, the phylogenetic relationships, and trait values.

Model components Parameter Prior Description
Demographic fluctuations s 500 (Assemblage) Sample size: number of individuals observed at present in a region of size Ne Ne 10,000 Regional assemblage or pool size at present P (1,000, 5,000) Regional assemblage or pool size in past refugia t c 100, 500, 1500, 2000, 3000 Beginning date of expansion (in generations before present), after which the regional assemblage size changed from P to Ne A = Ne/P (2, 10) Amplitude of expansion 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  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).

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.
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 2 = 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).

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 . 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 2 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 stemto-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).

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.

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 )-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 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.
amounted to more terminal branching events. Thus stem-totip 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 assemblagewide 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 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-totip asymmetry , 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 .
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 cooccurring taxa (Arbogast and Kenagy, 2001;Overcast et al., 2019Overcast et al., , 2020Swenson, 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 assemblagewide 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(Overcast et al., , 2020Swenson, 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 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.