Loïc Chalmandrier
email: loic.chalmandrier@ujf-grenoble.fr.

Tamara Münkemüller

Vincent Devictor

Sébastien Lavergne

Wilfried Thuiller

Decomposing changes in phylogenetic and functional diversity over space and time

Keywords: β-diversity, biodiversity, phylogenetic entropy, Shannon entropy, Hill numbers, diversity partitioning, bird assemblages, large-scale monitoring, turn-over

1. The α, β, γ diversity decomposition methodology is commonly used to investigate changes in diversity over space or time but rarely conjointly. However, with the everincreasing availability of large-scale biodiversity monitoring data, there is a need for a sound methodology capable of simultaneously accounting for spatial and temporal changes in diversity.

2.

Using the properties of Chao's index, we adapted Rao's framework of diversity decomposition between orthogonal dimensions to a multiplicative α, β, γ decomposition of functional or phylogenetic diversity over space and time, thereby combining their respective properties. We also developed guidelines for interpreting both temporal and spatial β-diversities and their interaction.

3. We characterised the range of β-diversity estimates and their relationship to the nested decomposition of diversity. Using simulations, we empirically demonstrated that temporal and spatial β-diversities are independent from each other and from α and γdiversities when the study design is balanced, but not otherwise. Furthermore, we showed that the interaction term between the temporal and the spatial β-diversities lacked such properties.

Introduction

Patterns of species diversity, as determined by their functional traits and phylogenetic relationships, have become central to addressing a large range of research questions such as the inference of assembly rules in community ecology [START_REF] Diamond | Assembly of species communities[END_REF][START_REF] Webb | Exploring the phylogenetic structure of ecological communities: an example for rain forest trees[END_REF][START_REF] Mouquet | Phylogenetic ecology: advances and perspectives[END_REF], or the delimitation of biodiversity hotspots in macro-ecology [START_REF] Mazel | Multifaceted diversity-area relationships reveal global hotspots of mammalian species, trait and lineage diversity[END_REF]. Using functional and phylogenetic diversity indices implicitly rejects the assumption that species are equally distinct entities, and instead accounts for their functional similarities and shared evolutionary history [START_REF] Mouquet | Phylogenetic ecology: advances and perspectives[END_REF][START_REF] Violle | Let the concept of trait be functional![END_REF].

Thanks to the extension of large-scale biodiversity monitoring [START_REF] Pereira | Towards the global monitoring of biodiversity change[END_REF] and the development of citizen science [START_REF] Bonney | Citizen science: a developing tool for expanding science knowledge and scientific literacy[END_REF]), large datasets have been made available for investigating the spatial and temporal dynamics of biodiversity [START_REF] Dornelas | Assemblage Time Series Reveal Biodiversity Change but Not Systematic Loss[END_REF]. These are of prime importance in evaluating how species assemblages are responding to on-going changes in climate and land uses. Since these temporal changes are not necessarily homogeneous across space, a depiction of biodiversity changes from both a spatial and temporal perspective [START_REF] Magurran | Long-term datasets in biodiversity research and monitoring: assessing change in ecological communities through time[END_REF]) is required to understand which processes contribute to biodiversity dynamics. An adequate methodology is therefore needed to produce meaningful measures of diversity changes over space and time.

In his seminal paper, [START_REF] Whittaker | Vegetation of the Siskiyou mountains, Oregon and California[END_REF] proposed breaking the regional species diversity (γdiversity) down into the average within-community species diversity (α-diversity) and the between community species diversity (β-diversity). More specifically, Whittaker formulated two laws to link α, β and γ-diversities: an additive law (γ = α+β) and a multiplicative law (γ = α × β). Two decomposition frameworks emerged from these two alternative approaches, each with different properties and drawbacks.

The additive law was adapted by [START_REF] Rao | Rao's axiomatization of diversity measures[END_REF] to the "Quadratic Entropy" index which generalised the Gini-Simpson index to include species dissimilarities such as functional or phylogenetic distances. He further proposed decomposing γ-diversity into several dimensions (e.g. space and time), a procedure called Anodiv [START_REF] Pavoine | Clarifying and developing analyses of biodiversity: towards a generalisation of current approaches[END_REF]. However; the additive decomposition of the γ-diversity, and by extension the Anodiv procedure, has been criticised for its inability to produce β-diversity estimates independent from the γand αdiversities [START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF][START_REF] Baselga | Multiplicative partition of true diversity yields independent alpha and beta components; additive partition does not[END_REF]. This property impedes Anodiv's ability to access temporal or spatial biodiversity changes on large spatial scales. Indeed, large-scale biodiversity monitoring typically covers numerous regions with variable γ and α-diversities (e.g. [START_REF] Devictor | Spatial mismatch and congruence between taxonomic, phylogenetic and functional diversity: the need for integrative conservation strategies in a changing world[END_REF], a consequence of the large-scale environmental filtering and historical contingencies that shape the biogeographic gradients of diversity [START_REF] Hawkins | Productivity and history as predictors of the latitudinal diversity gradient of terrestrial birds[END_REF]. Should β-diversities be compared across regions, it is vital that they only quantify spatial or temporal change within these regions, independently from changes in γ and αdiversities. Otherwise these two effects would become indistinguishable.

The second framework, based on Whittaker's multiplicative law, addresses this issue of independence. When diversity is calculated from an equivalent number [START_REF] Hill | Diversity and evenness: a unifying notation and its consequences[END_REF], it produces estimates of β-diversity which are independent from the α and γ-diversities [START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF][START_REF] Tuomisto | A diversity of beta diversities: straightening up a concept gone awry. Part 1. Defining beta diversity as a function of alpha and gamma diversity[END_REF]. Furthermore, the estimate of the β-diversity is set between 1 and the number of communities in the region studied. This property makes it possible to produce standardised estimates of β-diversity which are not dependent on the study design used in the region [START_REF] Chao | Proposing a resolution to debates on diversity partitioning[END_REF]. This is a particularly important feature since large-scale biodiversity monitoring systems tend to be spatially unbalanced [START_REF] Ficetola | Estimating patterns of reptile biodiversity in remote regions[END_REF]. However, despite these properties, Whittaker's multiplicative law was never adapted to breaking functional or phylogenetic diversity down into different dimensions.

We have built on these two frameworks and their respective advantages to propose a novel methodology for decomposing phylogenetic and functional diversity over space and time, and obtaining measurements of β-diversity which are independent of γ-diversity and αdiversity. This paper first introduces our multiplicative framework for estimating spatial and temporal beta diversities. Secondly, using a simulation-based approach, we demonstrate that, in the case of taxonomic diversity, the estimated β-diversities are pairwise independent from the γ and α-diversity and from each other. Finally, we illustrate its novelty and features in a case study by decomposing the spatio-temporal effects on the functional diversity of the common avifauna in four regions of France over the last decade.

Decomposing diversity over space and time

Definitions

We considered a region containing S sites in which species were recorded at T dates. We defined a community as the species composition of site s at a given date t. We defined a site pool as the pool of all communities at site s pooled for all dates, a time pool as the pool of all communities at a given date t pooled for all sites and the regional pool as the pool of all communities for all sites and dates. The spatio-temporal decomposition of diversity (multiplicative α, β, γ decomposition) will ultimately be calculated based on the ratios of the diversities in these different units (communities and pools).

Diversity index

To calculate the diversity of a given unit, we used Shannon entropy exponential. This index is an "equivalent number", part of the family of Hill numbers [START_REF] Hill | Diversity and evenness: a unifying notation and its consequences[END_REF]). As such, its value ranges from 1 (if one species makes up most of the total abundance in the unit) to the number of species in the unit (if their relative abundances are all equal). It can be interpreted as the number of "equally-abundant virtual species" in the unit [START_REF] Tuomisto | A diversity of beta diversities: straightening up a concept gone awry. Part 1. Defining beta diversity as a function of alpha and gamma diversity[END_REF]):

(Equation 1) With P, the vector {p 1 , p 2 , …, p N } of abundances of the N species present in the unit studied.

To include functional and phylogenetic similarities between species, we used the version of this index formulated by [START_REF] Chao | Phylogenetic diversity measures based on Hill numbers[END_REF] from Allen's phylogenetic entropy (2009).

(Equation 2)

where the summation is made over all branches of an ultrametric phylogenetic or functional tree of tips-to-root distance T, L(b) is the length of branch b, and p B denotes the vector containing for each branch b, the summed relative abundance of its descendent species. We will refer to this index as Chao's index.

Since the index includes species similarities, its absolute value can be interpreted as the number of "equally-abundant and fully distinct virtual species" in the study unit.

Spatio-temporal decomposition

We drew inspiration from the Anodiv procedure [START_REF] Rao | Rao's axiomatization of diversity measures[END_REF][START_REF] Pavoine | Clarifying and developing analyses of biodiversity: towards a generalisation of current approaches[END_REF] to decompose diversity according to two orthogonal factors, here time and space, according to a multiplicative framework using Chao's index. It is expressed as follows:

(Equation 3) P st is the vector of species relative abundance in community at site s and date t. The formulation includes ω st, a weight attributed to a community at site s and time t, that sums to 1 over all s and t based for instance on total abundance or species richness in the community. A perfectly balanced design will involve the absence of missing data and the equal weighting of all communities, i.e. for all s and t, . Unequal weighting is not compulsory but is typically relevant when communities have been sampled with different sampling efforts (for a discussion, see [START_REF] Hardy | Interpreting and estimating measures of community phylogenetic structuring[END_REF]. ω s and ω t are the weights of site pools s and time pools t, respectively, and are calculated as the sum of the weights of their constituent communities. The vector of species relative abundance for the site pool s and time pool t are thus calculated as the weighted mean of the species relative abundances in their constituent communities: and . Finally, the vector of species relative abundances in the species pool is calculated as the weighted mean of species relative abundance across all communities: Equation 3 can be reformulated as:

(Equation 4)

With γ being the γ-diversity of the study region; with α T and α S being respectively the mean α-diversity of time and site pools and α the mean α-diversity of the communities in the study region. Or more simply, (Equation 5) With β T being the temporal β-diversity, β S the spatial β-diversity and β ST the interaction term between the temporal and the spatial β-diversities. If the spatial and temporal structure of the dataset is ignored, the total β-diversity β of the region across time can be expressed as:

(Equation 6)

Properties of the β-diversities- [START_REF] Chiu | Phylogenetic beta diversity, similarity, and differentiation measures based on Hill numbers[END_REF] demonstrated that Chao's index obeyed the "replication principle". This implies that β T and β S have several of the properties enumerated in [START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF] and [START_REF] Tuomisto | A diversity of beta diversities: straightening up a concept gone awry. Part 1. Defining beta diversity as a function of alpha and gamma diversity[END_REF], which facilitate the interpretation of their numerical values:

β T and β S are pairwise independent from γ, and from α T and α S , respectively [START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF][START_REF] Baselga | Multiplicative partition of true diversity yields independent alpha and beta components; additive partition does not[END_REF]. Using simulations, we demonstrate below that β T and β S are pairwise independent from each other and are both pairwise independent from α.

-The values of β S (resp. β T ) are intuitive and can be interpreted as "the number of virtual, fully dissimilar, and equally abundant site pools (resp. time pools)" in the study region.

-The values of β S , β T and β have a range that is only dependent on the weights of respectively the site pools, time pools and the communities:

The minimum possible value of β S , β T and β will be 1 if the site pools, time pools or communities respectively are identical. The maximum possible value of β S , β T and β will be attained if the site pools, time pools or communities respectively do not share species or tree branches. N S , N T and N ST can be interpreted as the equivalent number of sites pools, times pools and communities, respectively. If all communities are weighted equally, they will be equal to S, T and ST, respectively.

In other words, our measurement of change in diversity over space (resp. over time) has a natural minimum and maximum. Thus the absolute values of β S and β T can be standardised by their minimum and maximum value to make their value independent from the number of sites and time periods studied [START_REF] Chao | Proposing a resolution to debates on diversity partitioning[END_REF]: and β ST has a minimum value of 1 and a maximum value constrained by both the value of β S and β T (see Appendix S2 for the demonstration),

Relationship to the nested decomposition of diversity-Another methodological choice that can be made when analysing a spatio-temporal dataset is to consider that space and time are nested (e.g. [START_REF] Sobek | Spatiotemporal changes of beetle communities across a tree diversity gradient[END_REF]). If space is considered as nested within time periods, a new measure β S/T can be formulated to characterise the mean spatial β-diversity of the region within time periods [START_REF] Pavoine | The apportionment of quadratic entropy: a useful alternative for partitioning diversity in ecological data[END_REF][START_REF] Tuomisto | A diversity of beta diversities: straightening up a concept gone awry. Part 1. Defining beta diversity as a function of alpha and gamma diversity[END_REF]).

The decomposition will then be expressed as:

(Equation 7)

meaning that Alternatively, if we consider time periods as nested within space:

(Equation 8) meaning that .

Thus, β S/T and β T/S are not strictly β-diversities because they are the ratio of two mean αdiversities from different hierarchical levels, rather than the ratio of a unit's diversity and the diversity of its subunits. However, we demonstrated (Appendix S1) that like β-diversities, they have fixed minimum and maximum values that are independent of the γ and αdiversities: and .

Interpretation of β-diversities

The different β-diversities can be interpreted on their own and in combination with each other (Figure 1). β T quantifies the change in diversity between time pools, or in other words, the temporal change in regional diversity, irrespective of the spatial patterns of diversity. β S quantifies the change in diversity between site pools, or in other words, the spatial change in diversity after averaging the temporal variability of communities. β ST quantifies the interaction between spatial and temporal turnover, it can be used to quantify finer changes such as a rearrangement of species between sites between two dates which are not quantified by β S and β T (Figure 1). The case where β ST is equal to 1 indicates that there is an identical change of diversity across space and time between communities; for instance, if between two dates a species is introduced in all studied communities at equal relative abundance. On the other hand, a value of β ST over 1 denotes a heterogeneity of change of communities over space and time that is not quantified by β S and β T because it averages out at larger spatial or temporal scales. It is interesting to note that it is possible to have a situation where β S and β T equal one while β ST is higher than one. This is illustrated in Figure 1 (β S = β T =1 and β ST > 1) with an extreme case in which two communities fully inversed their composition between the two dates studied. A concrete example could be the mosaic theory of forest regeneration [START_REF] Remmert | The mosaic-cycle concept of ecosystems-an overview[END_REF]: disturbances in a forested landscape would trigger the same temporal successions but at different times and locations generating a heterogeneous landscape. In this case, there is a spatial and temporal change between communities, but the time pools remain constant (in other words, the diversity of the region changes very little between the two dates) and the site pools also remain constant (in other words, when averaged over time, communities across the landscape have a similar composition).

Test of functional and phylogenetic β-diversities

A common way to test the values of functional or phylogenetic β-diversities is to use a randomisation model to generate a distribution of β-diversities under a certain null hypothesis. Our framework is compatible with any kind of randomisation procedure. In the following case study of French avifauna, we chose the species shuffling procedure that has been shown to be among the most efficient null models in terms of Type I error rate [START_REF] Hardy | Testing the spatial phylogenetic structure of local communities: statistical performances of different null models and test statistics on a locally neutral community[END_REF]. A significant high (resp. low) β-diversity thus indicates that species tend to be replaced by dissimilar (resp. similar) species over time or space. We calculated the effect size of each functional β-diversity, as the observed β-diversity minus the mean of its null distribution divided by the standard deviation of the null distribution. If the β-diversity was higher than expected (positive effect size) then the communities differed more than expected under a random assembly model; if the β-diversity was lower than expected (negative effect size), then the communities differed less than expected under a random assembly model.

Independence properties of the spatio-temporal diversity decomposition

We adapted the simulation procedures developed by [START_REF] Baselga | Multiplicative partition of true diversity yields independent alpha and beta components; additive partition does not[END_REF] to demonstrate the independence properties of our diversity decomposition framework. We used two simulation approaches: 1) a top-down approach where we first chose a γ value then generated a community weight vector and sequentially randomly selected the values of α T , α S and α; 2) a bottom-up approach where we first chose an α value, then generated a community weight vector and sequentially randomly selected the values of α T , α S and γ. Each draw was constrained by minimum and maximal values deduced from the properties of the spatiotemporal decomposition stated above. For a given number of sites (S) and dates (T), we tested 200 initial γ or α values between 1 and 200, each repeated 200 times. Both the topdown and bottom-up approach procedures are detailed in Appendix S3. Both procedures are necessary to demonstrate the pairwise independence of the β-diversities from γ and α (Baselga 2010). In the top-down approach, the β-diversities need to be uncorrelated with γ and in the bottom-up approach, the β-diversities need to be uncorrelated with α.

When no data were missing (i.e. all sites observed at all times) and T =4 and S =10, we found no correlation of γ with Std β T (Figure 2 ). When we further explored alternative T and S parametrizations (for any values of S and T between 2 and 10), we found these results to be robust (Table 1).

We also investigated the specific case where community data were missing in the dataset. For all T and S parameter values, a varying proportion of community weights were set to 0 while maintaining at least one community per site and date. We then studied how the amount of missing data affects the independence properties between γ, β T , β S , β ST and α. We found that in the extreme case of perfect balance (all community weights are equal), all correlations remained close to 0. When we introduced unequal weighting and increased the proportion of missing data, the correlation between β T and β S increased slightly but remained on average close to 0. However there were some extreme correlation values that deviated strongly from 0. We observed the same pattern for the other relationships although they tended to show more robustness (Figure 3). This can be explained by the fact that the sampling design becomes less orthogonal between sites and dates as the amount of missing data increases. The extreme case would be a spatial-temporal dataset where a single community was sampled per date, each time at different locations. Then β T and β S would be equal. However, as indicated by the correlation values, which were on average close to 0, the independence relationship remained on average quite robust and only a few weight vectors resulted in a loss of the independence properties.

Overall we therefore conclude that β T , β S and β ST are pairwise independent from α and γ.

We further conclude that β T and β S are pairwise independent from each other but only when the weighting scheme does not deviate too far from a perfectly balanced sampling design (i.e. no missing data and for all s and t, ).

Our simulation procedure has its limitations. Indeed, it is not clear how the inclusion of a phylogenetic or functional tree between species could further constrain the distribution of the different diversity metrics compared to the maximal values of β T , β S and β ST previously fixed. It is however intuitive that the maximum number of completely distinct communities, site pools or time pools is equal to the number of functional/phylogenetic tree branches that emerge from the root. Therefore although our simulation procedure concords with the current understanding of the diversity decomposition using Chao's index, it should be specified that our approach is more strictly appropriate for borderline cases where species are fully distinct or when the inclusion of a functional or phylogenetic tree does not affect the potential distribution of the diversity measures. In other cases, it is possible that the independence properties we empirically accessed would be altered.

Case study: spatio-temporal changes in functional diversity in French bird assemblages Data

We applied our spatio-temporal framework to the avifauna monitored by the French breeding bird survey programme [START_REF] Julliard | Spatial segregation of specialists and generalists in bird communities[END_REF]). This programme relied on skilled ornithologists to monitor common birds using a standardised protocol from 2001 to 2012. Under this scheme, ornithologists recorded every individual seen or heard during a fiveminute period at 10 count points, evenly distributed within 2×2 km survey sites. The sites were randomly selected around the observer's locality, thus ensuring that a variety of habitats were monitored (including intensive farmlands, forests, suburbs and cities). We selected four regions, each defined as a circular window 100 km in diameter, belonging to two disparate biogeographic regions. Two were situated on the Mediterranean coast (MED1 and MED2) and two on the Atlantic coast (ATL1 and ATL2). Each region included different numbers of survey sites for which temporal trends were available (i.e. with sites monitored twice at least five years apart). "Communities" were defined as the species assemblages recorded at the sites in the different regions. Species similarity was estimated from the ultrametric functional tree taken from [START_REF] Thuiller | The European functional tree of bird life in the face of global change[END_REF] based on body mass, diet and feeding behaviour (see Appendix S4 for details). Each community was given the same weight.

Mathematically, this means that for all s and t, . To further facilitate the interpretation of the result of the spatio-temporal decomposition, we used a double principal component analysis (dpcoa, [START_REF] Pavoine | From dissimilarities among species to dissimilarities among communities: a double principal coordinate analysis[END_REF] to visualise the differences between the sites in each region in terms of functional composition (Appendix S5 and Figure S2). A summary of the information on the four regions is available in Table 2.

Results

Spatial β-diversity

Absolute values of spatial β-diversities differed strongly between the four regions (Figure 4). We found that MED2>ATL2>MED1>ATL1 regarding spatial β-diversity. When focusing on standardised β-diversity, the ranking between regions changed: the spatial β-diversity of the region MED2 (which contained 14 sites) appeared much smaller (decreased rank) such as ATL2>MED1>MED2>ATL1. This was confirmed by the multivariate analysis, where MED2 appeared less spatially structured across the main multivariate axis compared to the other regions (Appendix S5, Figure S2). The region MED1 and the region ATL2 stood out for their high values of standardised spatial β-diversity which indicated that the time pools were more distinct in these regions than in the other two regions. This difference could be due to greater region-wide environmental heterogeneity in MED1 and ALT2 compared to MED2 and ALT1. The null model made it possible to determine that the site pools of the region MED1 were more functionally different (albeit marginally) than expected from their taxonomic composition compared to ATL2 despite having a lower standardised value of spatial β-diversity (MED1, effect size of β S = 1.22; ATL2, effect size of β T = 0.21; Figure 4). The multivariate analysis of MED1 (Figure S2) further suggested that sites at both dates were differentiated mainly according to the body size of their constituent species, suggesting a degree of large-scale environmental filtering acting on the birds' functional traits. On the other hand, the difference between the ATL2 site pools was not significantly different from the random expectation. This suggests that the large functional difference between site pools was due to a large number of different species (hence the high spatial β-diversity) but that between the site pools these species were not particularly distinct in terms of their functional traits.

Temporal β-diversity

Overall, the temporal change in regional diversity was more substantial in the two Mediterranean regions than in the two Atlantic regions (Figure 4). However, the use of the null model showed that the temporal β-diversity in MED1 was much higher than expected from the taxonomic change (effect size = 3.09) while in MED2, the temporal β-diversity was only marginally different from the taxonomic change (effect size = 1.48). This showed that between the two dates studied, the MED1 time pools were significantly different in terms of the functional traits of their constituent species (as suggested by the multivariate analysis, this was most likely due to the relative increase in larger-bodied species in all the sites over time, Figure S2). In contrast, the two Atlantic regions showed low temporal β-diversity, and the null model further showed that it was even lower than expected from the taxonomic change for ATL1 (ATL1, effect size of β T = -1.87; ATL2, effect size of β T = 0.15; Figure 4). This indicated that the sites studied in ATL1remained remarkably constant over time, with species being substituted by other species with similar functional traits. Overall, these results are consistent with previous diachronic analyses, which demonstrated substantial changes in bird communities in inland Mediterranean areas over time due to major changes in land use [START_REF] Preiss | Rural depopulation and recent landscape changes in a Mediterranean region: consequences to the breeding avifauna[END_REF][START_REF] Sirami | Vegetation and songbird response to land abandonment: from landscape to census plot[END_REF].

Interaction term

The interaction term β ST was always higher than 1 across the four regions studied. However, as the value of β ST is dependent on both β T and β S and the study design (number of sites and number of dates), it was difficult to compare it across regions. The null model provides a way of determining whether the interaction term was higher or lower than expected from the taxonomic turnover. We found that the interaction term of β ST was much higher than expected (effect size = 1.99) indicating that individual sites changed more in functional composition than expected from the temporal trend at regional level (β T ). Furthermore, the case of ATL1 is interesting because while β S and β T are both lower than expected, β ST is marginally higher than expected (effect size = 0.26). We interpreted this as a case of a region with communities whose functional composition had changed markedly over time, but in opposite directions. This indicates potential for investigating β ST as an emergent spatio-dynamic component of community changes.

Perspectives and limitations of spatio-temporal decomposition

The novel methodology presented herein has several advantages and strengths. Firstly, it allows a standardised decomposition of diversity across regions with potentially very different γ-diversities. In France, the functional and phylogenetic γ-diversities of avifauna are heterogeneous across space due to different macro-climatic influences [START_REF] Devictor | Spatial mismatch and congruence between taxonomic, phylogenetic and functional diversity: the need for integrative conservation strategies in a changing world[END_REF]. Most notably, the Mediterranean Basin has high taxonomic, functional and phylogenetic γ-diversity, which contrasts with the rest of France. Using our methodology, βdiversities are found to be pairwise independent from the local γ-diversities, thus making it possible to study processes such as landscape heterogeneity or land use change across regions without interference from the biogeographical differences quantified by the γdiversities. In contrast, there was a risk with the original Anodiv procedure [START_REF] Rao | Rao's axiomatization of diversity measures[END_REF] of yielding overestimated values for β-diversities in biogeographical areas with high γdiversities, regardless of the actual spatial or temporal change [START_REF] Baselga | Multiplicative partition of true diversity yields independent alpha and beta components; additive partition does not[END_REF].

Secondly, our methodology accounts for differences in sampling efforts between regions. In our case study, spatial and temporal functional β-diversities were standardised using the maximum and minimum value they could possibly attain in a region (Figure 4), thus producing estimates of β-diversities that were unrelated to the number of sites or dates studied within each region [START_REF] Chao | Proposing a resolution to debates on diversity partitioning[END_REF]. Although our method does not prevent bias if the sampling is not representative of the biodiversity in a region over space and time, the standardisation makes it possible to compensate for an unbalanced study design between regions. This solves a common problem arising from large-scale biodiversity monitoring where remote areas tend to be under-sampled [START_REF] Jiguet | French citizens monitoring ordinary birds provide tools for conservation and ecological sciences[END_REF][START_REF] Ficetola | Estimating patterns of reptile biodiversity in remote regions[END_REF].

Thirdly, our approach also makes it possible to test the space-for-time substitution often used when time-series data are not available. The drivers of diversity change such as climate (e.g. [START_REF] Blois | Space can substitute for time in predicting climate-change effects on biodiversity[END_REF] or land use (e.g. [START_REF] Sirami | Vegetation and songbird response to land abandonment: from landscape to census plot[END_REF]) can be studied both across space and time. Ecologists have thus traditionally used space-for-time substitution as an alternative to expensive and rare long-term studies [START_REF] Pickett | Space-for-time substitution as an alternative to long-term studies. Long-term studies in ecology: approaches and alternatives[END_REF]e.g. Chalmandrier et al. 2013). This substitution assumes that changes in diversity over spatial locations and changes in diversity over time are equivalent and independent under the assumption that they are driven by the same ecological process [START_REF] Fukami | Long-term ecological dynamics: reciprocal insights from natural and anthropogenic gradients[END_REF]. However, this assumption can easily be violated by confounding processes such as dispersal [START_REF] Brotons | Colonization of dynamic Mediterranean landscapes: where do birds come from after fire?[END_REF], biotic interactions [START_REF] Thuiller | Stochastic species turnover and stable coexistence in a species-rich, fire-prone plant community[END_REF]) and delayed responses to changes in the local environment (Devictor et al. 2012). Our methodology provides tools which are adapted to testing the assumption on which space-for-time substitution is based: the pairwise independence of β S and β T allows for the direct comparison of the spatial and temporal components of changes in diversity. However, we showed that this independence property was only maintained if the sampling design is balanced (i.e. there is relatively a low amount of missing data and communities have a similar weight). We therefore recommend testing the sampling design beforehand using the simulation procedures to assess whether β S and β T are theoretically pairwise independent.

Fourthly, this is the first study to adapt the diversity decomposition between multiple factors originally proposed by [START_REF] Rao | Rao's axiomatization of diversity measures[END_REF] to the requirements of β-diversities computations, as recommended by [START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF] and [START_REF] Tuomisto | A diversity of beta diversities: straightening up a concept gone awry. Part 1. Defining beta diversity as a function of alpha and gamma diversity[END_REF]. Furthermore, we generalise this approach using Chao's index (2009) which includes species' functional and phylogenetic distances, thereby combining the advantages of both methods. Shannon entropy exponential and its generalisation are the only equivalent numbers that fully combine the properties of the additive and multiplicative α, β, γ decomposition [START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF]. This opens promising avenues for adapting our framework to methodologies based on the additive decomposition of the Shannon entropy (e.g. Pelissier & Couteron 2007). It is also the only equivalent number where 1) there is more or less a general consensus about the handling of unequal weighting of communities (but see [START_REF] Chiu | Phylogenetic beta diversity, similarity, and differentiation measures based on Hill numbers[END_REF] and 2) the unequal weighting still leads to β-diversities values that are pairwise independent from γ and α-diversities [START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF][START_REF] Tuomisto | A diversity of beta diversities: straightening up a concept gone awry. Part 1. Defining beta diversity as a function of alpha and gamma diversity[END_REF]).

Chao's index belongs to a large family of indices that extend the Hill numbers (1973) to include species' phylogenetic similarities [START_REF] Chao | Phylogenetic diversity measures based on Hill numbers[END_REF]) making it possible to explicitly parametrize the weight given to a rare vs. a dominant species. While our framework is transposable to these indices, some properties (pairwise independence of β-diversities from γ and α-diversities, range of β ST and of nested β-diversities) need to be demonstrated, in particular in the case of missing data and unequal weighting of communities. The Chao's index studied is based on assumptions on how to take into account species abundances, i.e. the contribution of a species to the diversity value is proportional to its relative abundance [START_REF] Chiu | Phylogenetic beta diversity, similarity, and differentiation measures based on Hill numbers[END_REF]. It thus may not be suitable for achieving certain analytical aims: for instance, a conservation approach may want to consider rare and dominant species equally regardless of their relative abundance. On the other hand, a focus on ecosystem functioning may require an emphasis on dominant species as they are expected to be the main contributors to ecosystem functioning [START_REF] Garnier | Plant functional markers capture ecosystem properties during secondary succession[END_REF]; but see [START_REF] Mouillot | Rare Species Support Vulnerable Functions in High-Diversity Ecosystems[END_REF]. Furthermore, recent work has shown the value of analysing diversity patterns with multiple equivalent numbers in order to vary the weighting given to dominant as opposed to rare species, and to disentangle multiple assembly rules [START_REF] Arroyo-Rodriguez | Plant βdiversity in fragmented rain forests: testing floristic homogenization and differentiation hypotheses[END_REF]Chalmandrier et al. in press). We therefore argue for more statistical development to adapt spatio-temporal decomposition to other equivalent numbers, using the generalisation of Hill numbers, and thus adding a supplementary parameter which explicitly examines the impact of species' relative abundances.

Conclusion

Recent years have seen major efforts to unify methodologies for evaluating and decomposing assemblage diversity. We have drawn on these achievements to propose a methodology that overcomes the challenges encountered when studying large-scale diversity datasets which encompass multiple orthogonal dimensions. We have shown that this approach can be used with classical animal survey data (also available for butterflies, fishes and plants) and that it provides clear results. Although more work is required to expand this method to multiple diversity indices, we believe that the properties of our methodology open up promising avenues for evaluating and testing diversity change across multiple dimensions. This will allow thorough analyses of the ever-increasing data produced by biodiversity survey programs worldwide. Functional 1 β-diversity decomposition between spatial (dark-grey), temporal (grey) and spatial-temporal interaction (light-grey) components for each of the four regions. Each bar represents the absolute value of the β-diversity (top) or its value standardised by its minimum and maximal possible value (bottom). Numbers above each bar show the effect size against the null model. A negative value indicates a higher than expected value of βdiversity while a positive value indicates a lower than expected value of β-diversity.

  , r = 0.0039; 95% CI interval: [-0.0059, 0.014]), γ with Std β S (Figure 2, r = -0.00027; 95% CI interval: [-0.010, 0.0095]) or Std β T with Std β S (Figure 2, r = 0.02; 95% CI interval: [-0.013, 0.0069]) in the top-down approach and no correlation of α with Std β T or α with Std β S in the bottom-up approach (Figure S1. α and Std β T : r = 0.0052; 95% CI interval: [-0.0046, 0.015]). In contrast, β ST depended on the previously established β T and β S values, but is independent from both γ and α (Figure 2. β ST and γ: r = 0.00038; 95% CI interval: [-0.0094, 0.010]. Figure S1. β ST and α: r = 0.0037; 95% CI interval: [-0.0061, 0.0135]

Figure 1 .

 1 Figure 1. Graphical interpretation of the eight possible patterns of diversity change over space and/or time and their respective βS, βT, βST value. Each cell of the table represents one of the patterns illustrated by the composition of four communities from two sites (lines) at two dates (column), each represented by a circle. Sites and time pools (see Methods) composed from the pooled composition of communities over sites or dates, respectively, are represented by rectangles. Within communities and pools, symbols represent individuals from a species specified by the geometrical shape.

Figure 2 .

 2 Figure 2.Pairwise scatterplot of the value of γ, Std β T , Std β S , β ST and α obtained from a "top-down" simulation procedure. T was equal to 4 and S to 10. The panels on the diagonal represent the distribution of each diversity estimate over the simulations. The panels from the lower triangle represent the pairwise relationship between two of diversity estimates and the panels from the upper triangle contain the Pearson correlation coefficient between two of diversity estimates.

Figure 3 .

 3 Figure 3.Boxplots of the correlation coefficients obtained from the simulations as a function of the percentage of missing data (e.g. community weights set to 0). The boxplot on the far left (designated by "B") illustrates the case of perfect balance (all community weights are equal). The figure displays the correlations between (a) γ and Std β T given by the top-down procedure, (b) γ and Std β S given by the top-down procedure, (c) γ and β ST given by the topdown procedure, (d) Std β T and Std β S given by the top-down procedure (e) Std β T and α given by the bottom-up procedure, (f) Std β S and α given by the bottom-up procedure and (g) β ST and α given by the bottom-up procedure.

  Figure 4.

Methods Ecol Evol. Author manuscript; available in PMC 2015 February 13.

Acknowledgements

The research leading to these results received funding from the European Research Council under the European Community's Seven Framework Programme FP7/2007-2013 Grant Agreement no. 281422 (TEEMBIO). TM was funded by the ANR-BiodivERsA project CONNECT (ANR-11-EBID-002), as part of the ERA-Net BiodivERsA 2010 call. The LECA is part of Labex OSUG@2020 (ANR10 LABX56).

Supporting informations Additional supporting information may be found in the online version of this article. Appendix S1. Minimum and maximum value of nested β-diversities. Appendix S2. Minimal and maximal value of β ST . Appendix S3. Simulation procedures. Appendix S4. Details of the construction of the bird functional tree. Appendix S5. Comments on the multivariate analysis of the spatio-temporal changes of bird functional assemblages. Table S1. Notations for Appendix S1 and Appendix S2.

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.


Marta Carboni
email: marta.carboni@gmx.net

Tamara Münkemüller

Sébastien Lavergne

Philippe Choler

Benjamin Borgy

Cyrille Violle

Franz Essl

Cristina Roquet

François Munoz

Divgrass Consortium

Wilfried Thuiller

What it takes to invade grassland ecosystems: traits, introduction history and filtering processes

Keywords: Alien species, Functional and phylogenetic similarity, Invasibility, Invasion success, Invasiveness, Rabinowitz

Whether the success of alien species can be explained by their functional or phylogenetic characteristics remains unresolved because of data limitations, scale issues and weak quantifications of success. Using permanent grasslands across France (50,000 vegetation-plots, 2000 species, 130 aliens) and building on the Rabinowitz' classification to quantify spread, we showed that phylogenetic and functional similarities to natives were the most important correlates of invasion success compared to intrinsic functional characteristics and introduction history. Results contrasted between spatial scales and components of invasion success. Widespread and common aliens were similar to co-occurring natives at coarse scales (indicating environmental filtering), but dissimilar at finer scales (indicating local competition). In contrast, regionally widespread but locally rare aliens showed patterns of competitive exclusion already at coarse scale. Quantifying trait differences between aliens and natives and distinguishing the components of invasion success improved our ability to understand and potentially predict alien spread at multiple scales.

Introduction

The repeated introduction of alien species is a major component of ongoing global changes and a major threat to global biodiversity [START_REF] Walther | Alien species in a warmer world: risks and opportunities[END_REF][START_REF] Vila | Invasive species. What everyone needs to know[END_REF]. Understanding and predicting which introduced species will spread, naturalize and finally become invasive requires two major steps: 1) characterizing invasion success and 2) identifying its driving forces. Traditionally, invasion success was studied by focusing either on intrinsic biological characteristics conferring invasiveness, or on characteristics of the resident communities favoring invasion (invasibility; [START_REF] Richardson | Plant invasions: merging the concepts of species invasiveness and community invasibility[END_REF]. Additionally, historical factors, such as time since introduction or introduction pathways, also play a role [START_REF] Wilson | Residence time and potential range: crucial considerations in modelling plant invasions[END_REF][START_REF] Mahoney | Introduction effort, climate matching and species traits as predictors of global establishment success in non-native reptiles[END_REF]. Nevertheless, attempts to link the degree of success across species to intrinsic traits, community invasibility and introduction history in a combined way have been rare and have identified few generalities (e.g. [START_REF] Thuiller | Interactions between environment, species traits, and human uses describe patterns of plant invasions[END_REF]). One reason is that our ability to test the importance of complementary factors on invasion processes hinges on the development of integrative frameworks to better quantify and measure the invasion success of alien species [START_REF] Colautti | Quantifying the invasiveness of species[END_REF].

Studying the success of alien species requires detailed quantification of both their local performance and geographic spread. Most studies attempting to identify alien species of greatest concern have focused on frequencies of species occurrence. However, focusing only on occurrence data tends to over-estimate the success of species that occur in low abundances but are widespread, while under-estimating the success of species that dominate where present but are not widely distributed [START_REF] Colautti | Quantifying the invasiveness of species[END_REF]). Here we propose to account for different components of invasion success by applying for the first time the Rabinowitz' classification of rarity and abundance [START_REF] Rabinowitz | Seven forms of rarity. The Biological Aspects of Rare Plant Conservation[END_REF], Fig. 1) to alien species. Combining frequency of occurrence with geographic range, local abundance and niche breadth gives a more comprehensive measure of invasion success in a given area. In this framework, aliens can be considered most successful when they are regionally widespread, locally abundant and have broad ecological preferences.

Pinpointing intrinsic functional traits associated with invasiveness to profile successful invaders has been a long lasting quest in invasion ecology. For plants, a number of traits have been consistently found to facilitate invasion [START_REF] Richardson | Plant invasions: merging the concepts of species invasiveness and community invasibility[END_REF][START_REF] Van Kleunen | A meta-analysis of trait differences between invasive and noninvasive plant species[END_REF]. Small seed sizes, prolific seed production and short juvenile periods regularly emerge as important determinants of invasiveness, as they correlate with long dispersal distances and high competitive potential at the juvenile stage [START_REF] Richardson | Plant invasions: merging the concepts of species invasiveness and community invasibility[END_REF]. Fast relative growth rates are thought to favor highly invasive species in frequently disturbed or earlysuccessional ecosystems [START_REF] Rejmanek | What attributes make some plant species more invasive?[END_REF]Dawson et al. 2011b;[START_REF] Ordonez | Do alien plant species profit more from high resource supply than natives? A traitbased analysis[END_REF]. Finally, studies highlighted that invasive plant species tend to be larger than noninvasive ones, probably in relation with the ability to capture light [START_REF] Van Kleunen | A meta-analysis of trait differences between invasive and noninvasive plant species[END_REF].

Overall, the finding that successful alien plants share sets of functional characters is also supported by studies showing that certain phylogenetic lineages have higher proportions of invasive aliens than expected by chance [START_REF] Lambdon | Is invasiveness a legacy of evolution? Phylogenetic patterns in the alien flora of Mediterranean islands[END_REF]). Hence, it seems possible to screen successful alien plant species based on trait syndromes, such as the leaf-height-seed (LHS) spectrum [START_REF] Westoby | A leaf-height-seed (LHS) plant ecology strategy scheme[END_REF], or based on their phylogenetic origins, which will capture biogeographic history and potential broad bioclimatic adaptations.

While traits conferring invasiveness can favor invasion, they do not guaranty its success. The abiotic and biotic characteristics of the recipient communities also determine invasion opportunities (invasibility, [START_REF] Thuiller | Resolving Darwin's naturalization conundrum: a quest for evidence[END_REF]. For instance, if environmental filtering is the main process determining alien species spread, the most successful aliens should share similar ecological adaptations to natives [START_REF] Diez | Darwin's naturalization conundrum: dissecting taxonomic patterns of species invasions[END_REF]. If instead certain trait values confer greater competitive ability, the most successful aliens should be the ones that rank higher in the trait hierarchy [START_REF] Kunstler | Competitive interactions between forest trees are driven by species' trait hierarchy, not phylogenetic or functional similarity: implications for forest community assembly[END_REF]Gallien et al. 2014b). Conversely, if we assume more intense competition for resources between ecologically-similar species, local communities are likely to be resistant to aliens that have very similar traits or are closely related to the resident natives [START_REF] Darwin | The Origin of Species[END_REF][START_REF] Elton | The ecology of invasions by plants and animals[END_REF]. Based on these assumptions, specific metrics of functional and phylogenetic similarity have been proposed in an attempt to infer the processes driving invasion [START_REF] Thuiller | Resolving Darwin's naturalization conundrum: a quest for evidence[END_REF]Gallien et al. 2014a). The application of these metrics has highlighted that invasion dynamics can broadly vary depending on local ecosystem properties, spatial scale, and the amount of available resources [START_REF] Diez | Darwin's naturalization conundrum: dissecting taxonomic patterns of species invasions[END_REF][START_REF] Carboni | Darwin's naturalization hypothesis: scale matters in coastal plant communities[END_REF], Gallien et al. 2014b).

However, to date very few studies have explicitly investigated how multiple, interacting drivers of invasion success vary among alien species [START_REF] Pellock | Validity of Darwin's naturalization hypothesis relates to the stages of invasion[END_REF]) and fewer still have considered different components of invasion success at the same time [START_REF] Lemoine | Phylogenetic relatedness and leaf functional traits, not introduced status, influence community assembly[END_REF]. If specific processes are important for regionally widespread, locally abundant and generalist alien species compared to less successful ones, identifying processes linked to invasion success seems of utmost importance. Ultimately, the identification of these crucial processes can help anticipating the most concerning invasions.

Here we tackled the challenge of relating invasion success of alien species to intrinsic traits conferring invasiveness and to processes determining invasibility of native communities, while accounting for species' introduction history. First, we developed a framework for quantifying invasion success based on each alien species' invaded range area, local abundance and realized niche breadth in native communities. Second, we asked whether functional-and phylogenetic-based measures of community invasibility are good predictors of the degree of invasion success of aliens compared to their intrinsic traits and introduction history. While our initial aim was to pinpoint the processes determining invasion success, we also tested our framework in a predictive context. We did so in the context of temperate, C3 dominated permanent grassland communities, based on an extensive dataset spanning the whole of France. Investigating invasion success is of high importance in these grasslands that are already severely threatened by urban development and the abandonment of traditional agricultural practices [START_REF] Wesche | Fifty years of change in Central European grassland vegetation: Large losses in species richness and animal-pollinated plants[END_REF][START_REF] Essl | Diversity of native and alien vascular plant species of dry grasslands in central Europe[END_REF].

Methods

Study system, habitats and vegetation plots

Our study focuses on permanent grasslands (hereafter "grasslands") defined as semi-natural herbaceous ecosystems supporting livestock farming systems. They are usually characterized by low external inputs (e.g. feed, fertilizer, energy). In France the area of grasslands is declining since the 1970s, but still covers ca. 10 million hectares [START_REF] Violle | Vegetation ecology meets ecosystem science: Permanent grasslands as a functional biogeography case study[END_REF].

We used a comprehensive database compiled as part of the DIVGRASS project [START_REF] Violle | Vegetation ecology meets ecosystem science: Permanent grasslands as a functional biogeography case study[END_REF] and reporting the composition of more than 50,000 vegetation plots throughout France, geo-referenced at a precision of around 5 km on average, for a total of 4280 species. The plant communities sampled consist of semi-natural grasslands (mesic, ruderal and trampled grasslands) created and maintained by mowing and grazing, and natural grasslands (mountain and dry calcareous grasslands) covering wide ecological gradients. The data consist of estimated relative cover of all species present within an area ranging from 25 to 100 m 2 (Appendix S1).

Each plot was assigned to one of four broadly defined grassland habitats based on an analysis of modularity of the overall species-site network [START_REF] Carstensen | The functional biogeography of species: biogeographical species roles of birds in Wallacea and the West Indies[END_REF]) and using the algorithm of [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF]. The four habitats were: (i) dry calcareous grasslands, (ii) mountain grasslands, (iii) mesic grasslands and (iv) ruderal and trampled grasslands.

Alien species and introduction history

We focused on 160 herbaceous alien species defined as naturalized in France (Table S2) within the DAISIE database (www.europe-aliens.org), which provides an up-to-date inventory of species alien to European countries [START_REF] Pyšek | Alien vascular plants of Europe. Handbook of alien species in Europe[END_REF]). Both neophytes (new alien species introduced > c. 1500) and few archeophytes (old alien species, introduced < c. 1500) were listed in this database.

Species' introduction history (year of first record and main introduction pathway) in France was taken from DAISIE. If no record was available for France, we used the information from bordering countries as an approximation. We removed species with no information on introduction history leading to a final set of 127 alien naturalized species. To quantify minimum residence time, we determined the year of first record as the earliest year of all the records available for each alien species. Introduction pathways were based on the DAISIE hierarchical classification system (Table S1). We summarized this information by identifying the most common introduction pathway among "Agriculture", "Ornamental and horticulture", "Unintentional through transport" and "Other".

Phylogeny and traits

We reconstructed a molecular and dated genus-level phylogeny for the entire species pool based on DNA sequences available in GenBank [START_REF] Roquet | Building megaphylogenies for macroecology: taking up the challenge[END_REF]; Appendix S2). To calculate distance-based phylogenetic metrics (see below), we extracted the cophenetic distance between all species-pairs from the phylogenetic tree.

We collated information on four functional traits: Specific Leaf Area (SLA; the ratio of leaf area to dry mass), plant maximum height at maturity (Height), seed mass (SM) and growth form. These trait data were extracted from a number of databases and local datasets (see Appendix S3 and Violle et al. 2015 for details). The leaf-height-seed (LHS) scheme captures variation in resource acquisition, establishment and competitive ability among species [START_REF] Westoby | A leaf-height-seed (LHS) plant ecology strategy scheme[END_REF]. SLA reflects species relative growth rate (faster growing species having higher SLA values), nutrient economy and water use efficiency. Seed mass reflects dispersal capability, seed production and longevity on the one hand, and species survival at the seedling stage on the other. Height captures species' ability to dominate the vegetation layer in a given habitat and to intercept light. Finally, we characterized growth form as either graminoid, forb or dwarf shrub. We kept only species for which we had information on at least two traits (2930 species from the initial 4280) and estimated the remaining missing values based on multivariate imputation by chained equations (van Buuren & Groothuis-Oudshoorn 2011; Appendix S3). To calculate distance-based functional metrics (see below), we compiled a pair-wise species dissimilarity matrix for all traits combined using the Gower distance (function "daisy" in the "cluster" package), as a multivariate measurement of functional dissimilarity.

Determining invasion success

We characterized invasion success based on the three components of the Rabinowitz scheme (Fig. 1). First, to quantify regional distribution we calculated the species' frequency of occurrence, a measure that correlated strongly also with geographical spread (Fig. S1). Second, local abundance was taken from the average percentage cover of each alien species within the plots where it occurred. Third, niche breadth was calculated using a hypervolume approach based on three bioclimatic variables and one productivity measure hypothesized to be important for the naturalization and spread of alien species. Climatic variables were derived from the French meteorological model Aurelhy [START_REF] Bènichou | Prise en compte de la topographie pour la cartographie des champs pluviomÈtriques statistiques[END_REF], based on interpolated measurements downscaled at 100 m resolution and summarizing climatic information over the last 30 years . We focused on average annual precipitation, average daily minimum temperature in winter (December-February) and average daily maximum temperature in summer (June-August). We used the Normalized Difference Vegetation Index (NDVI), summed over the whole year as a proxy of gross primary productivity. Data were obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) for the years 2000-2014 and processed as described in [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations of snow cover duration[END_REF]. Based on these variables, we used multivariate kernel density estimation to construct a 4dimensional ecological hypervolume for each alien species [START_REF] Blonder | The n-dimensional hypervolume[END_REF]). All variables were rescaled between 0 and 1 (by subtracting the minimum and dividing by the range) for comparability. Hypervolumes were constructed using a Silverman bandwidth estimator and a 0% quantile threshold [START_REF] Blonder | The n-dimensional hypervolume[END_REF]). We then used the calculated volumes in ecological space as proxies of niche breadth for each alien species.

Based on these three components of invasion success (regional distribution, local abundance and niche breadth), we then used two complementary approaches to characterize alien success: 1) we separated them into eight classes (corresponding to the classical Rabinowitz classes -Fig. 1), and 2) we synthesized the three components into a synthetic continuous index of success. The first approach categorizes aliens into invasion classes that might require specific management strategies, while the second ranks species in terms of their overall invasion success with the aim of identifying the features that can best explain it.

For the first (categorical) approach, we first split species into common vs. rare groups for each Rabinowitz' criterion independently (number of plots, local abundance and niche breadth), based on the sum of squares from species to group centers. Then the species were assigned to one of the eight Rabinowitz classes based on the cross-tabulation of the groups obtained for the three criteria (Figs. 1 and2).

For the second (continuous) approach, we calculated a synthetic index representing the species scores on the first axis of the principal component analysis (PCA) performed on the three components of invasion success. This PCA axis explained roughly 60% of the variability in the data (Figure S1). It was most strongly linked to regional distribution, but it also captured niche breadth and local abundance. This "invasion success" indicator (or PCA1) was then used as a response variable in subsequent analyses to identify the most important drivers of invasion success.

Phylogenetic and functional dissimilarity of alien species

We calculated a set of functional and phylogenetic similarity metrics to native species (invasibility metrics, [START_REF] Thuiller | Resolving Darwin's naturalization conundrum: a quest for evidence[END_REF]) based on pairwise co-occurrences at two scales (Appendix S4). We first identified the native species that most often co-occurred with each alien species by using the V-score as a measure of species co-occurrence [START_REF] Lepš | Multivariate Analysis of Ecological Data using CANOCO[END_REF], modified to account for species abundances). V-scores were calculated at both the plot (local community) and the habitat scale (defined by the set of plots belonging to the same habitat). Then, for each alien species we calculated 1) the (absolute) Mean Distance to the Most often Co-occurring Species (MDMCS), and 2) whether it had higher or lower values than the natives for each trait (i.e. its hierarchical position on each trait gradient, [START_REF] Kunstler | Competitive interactions between forest trees are driven by species' trait hierarchy, not phylogenetic or functional similarity: implications for forest community assembly[END_REF]. MDMCS was based on both phylogenetic and functional distances, whereas the hierarchical index was calculated for each trait independently (Height.hier, SLA.hier and SM.hier). Overall these metrics measure the relative invasibility of native grassland assemblages at multiple scales (Table 1, Appendix S4). In the framework of the hierarchical filtering model of invasion biology (Richardson & Pysek 2006), the habitat-level analyses provide insights into environmental filtering, whilst plot-level analyses allow recovering signals of biotic filtering.

Disentangling the factors explaining invasion success in grasslands

Modeling species independently of each other may bias analyses if closely related species have similar success. We tested for a phylogenetic signal of our continuous indicator of invasion success (i.e. PCA1) on the phylogenetic tree of alien species using Blomberg's K statistic [START_REF] Blomberg | Tempo and mode in evolution: phylogenetic inertia, adaptation and comparative methods[END_REF] and found no significant signal (Appendix S5; Fig. 3; K = 0.074, PIC variance p-value = 0.931). Hence, we considered alien species as phylogenetically independent in further analyses of invasion success.

We used boosted regression trees (BRT) to model invasion success of alien species as a function of introduction history (pathway and residence time), functional traits related to invasiveness (SLA, height, seed mass and growth form) and community invasibility metrics (functional and phylogenetic MDMCS at habitat and plot scales and trait-hierarchy metrics). BRTs allow handling different types of predictor variables, can fit nonlinear relationships, are robust to outliers, automatically handle interaction effects between predictors and are able to extract robust estimations of variable importance [START_REF] Elith | A working guide to boosted regression trees[END_REF]). We estimated variable importance for explaining invasion success based on the reduction of squared error attributable to each variable, normalized to sum up to 100 (relative influence). We estimated the goodness-of-fit of the overall model by assessing the R-square of fitted values against observed values. Models were fitted in R (R-Development-Core-Team 2014), using the 'dismo' package [START_REF] Elith | A working guide to boosted regression trees[END_REF].

Finally, to test whether our approach for understanding invasion success was also able to predict which species from a given pool might become successful, we conducted a set of cross-validations through 1) repeated split-sampling and 2) jackknifing, using several goodness-of-prediction statistics to assess performance (Appendix S6 for details).

Results

Invasion success of alien plant species in French grasslands

Based on our three components of invasion success, we classified alien species into the eight abundance classes of the Rabinowitz framework (Fig. 2, Table S2). In particular, alien species were considered to be rare when: (1) they occurred in less than 13 plots, (2) covered in average less than 2.5% of total cover in the plots where they occurred, and (3) occupied an ecological hypervolume of less than 0.008 sd 4 .

Most alien species in French grasslands were rare, with almost 1/3 of species being in the rarest category (narrow geographic range, sparse abundance and small niche breadth). Fortyfour aliens were classified as abundant for at least two criteria (e.g. Ambrosia artemisiifolia) and 20 alien species fell into the most successful alien species category (class A, see Table S2 for the list of species).

The synthetic index of invasion success gave very consistent results with the categorical classification of invasion success (Fig. 2, right panel). The rarest alien species had the lowest values of the invasion success indicator, species that were common for only one criterion had low to medium values, while species common for at least two criteria had medium to high success indicator values. This result supports the dimension reduction approach through PCA as a mean of synthesizing the Rabinowitz multi-dimensional scheme into a onedimensional continuous success indicator.

Disentangling the factors explaining invasion success

Community invasibility metrics were collectively more important for explaining the invasion success of alien species, than either intrinsic trait values or introduction history (Fig. 4). In particular, the dissimilarity to species co-occurring in the same habitat type was by far the most important predictor of success both from a functional and phylogenetic perspective. Specifically, the most successful alien species were phylogenetically close and functionally similar to the native species co-occurring in the same habitat (Fig. 5A-B), but dissimilar to those co-occurring in the same plots (Fig. 5E). Minimum residence time, a classical feature for explaining invasion success, was also important, but only in third place. Further, having higher or lower seed mass compared to the natives co-occurring within plots (SM.hier.plot) was also determinant. Thereby, the most successful alien species were introduced long ago and had smaller seeds, greater SLA and lower height than the natives co-occurring within plots (Fig. 5C,D,F; Fig. S2). Importantly, while the relative trait hierarchy of alien and native species contributed to invasion success, intrinsic trait values of aliens species did not. Overall, this set of historical, functional and invasibility variables explained well the invasion success of alien species in French grasslands (R 2 : 0.571).

Moreover, contrasting phylogenetic similarity patterns at the habitat scale suggests distinct underlying drivers for the different components of success (Fig. 5A). When considering variation across the Rabinowitz' classes (Fig. S3), we found that different strategies could favor different components of invasion success. In congruence with the BRT results, alien species closely related to natives became both widespread and abundant (classes A and E; Fig. S3, left panel). By contrast, alien species that were phylogenetically most distant from natives became widespread regionally but not locally abundant (class B). In other words, the species in class A were those indicating high invasion success of close relatives to the natives in the BRT response curve (left part of Fig. 5A), while species in class B corresponded to the moderate invasion success obtained by distantly related species (right part of Fig. 5A). Finally the species with very low invasion success (rare aliens) had intermediate phylogenetic distances to the natives (central part of Fig. 5A).

Predicting invasion success

While our framework is mainly designed for for identifying the factors driving the success of invasions in grassland ecosystems, we found that our model also had a reasonable predictive performance to estimate invasion success of newly introduced species (Pearson R = 0.56 and R-squared = 0.33; Tables S3,Appendix S6). In particular, frequent, abundant and generalist species were on average correctly predicted as having greater invasion success than the rare ones (Fig. S11). The highest prediction errors and under-predictions occurred for the invaders of class C (regionally scarce, but locally abundant and generalist) and for some of the species of class A (Fig. S12).

Discussion

By explicitly considering different components of invasion success we were able to unravel the drivers of plant invasions in French permanent grasslands. We found support for many classical hypotheses relating invasion success to species' traits and introduction history. Moreover, we were able to identify signals of competition at large spatial extent and for highly diverse communities. But more importantly, we showed that functional and phylogenetic similarities to natives were by far the most important features, suggesting that they are excellent proxies for processes of environmental filtering and competition driving invasion success. Specifically, environmental conditions filtered aliens similar to the natives at the habitat scale, while a combination of hierarchical competition and limiting similarity led to greater success of shorter aliens with smaller seeds and faster growth rates than the natives at fine scale.

We confirmed that most alien species introduced in French grasslands remained rare and failed to establish large and successful populations [START_REF] Richardson | Plant invasions: merging the concepts of species invasiveness and community invasibility[END_REF]. Interestingly, certain aliens considered a threat in other habitats or countries, (e.g. Cortadeira selloana, one of the 100 worst invasive species in Europe; [START_REF] Pyšek | Alien vascular plants of Europe. Handbook of alien species in Europe[END_REF], appear to still be rare in French grasslands demonstrating the importance of the native community and of the historical context. Still, we identified twenty alien species as highly successful (Fig. 2). Several of them, including Erigeron annuus, Oenothera biennis and Paspalum dilatatum, are already known to be invasive elsewhere, causing negative environmental impacts in a number of different regions and habitats [START_REF] Weber | Assessing the risk of potentially invasive plant species in central Europe[END_REF][START_REF] Pyšek | Alien vascular plants of Europe. Handbook of alien species in Europe[END_REF][START_REF] Kumschick | Comparing impacts of alien plants and animals in Europe using a standard scoring system[END_REF]. A few species, such as Avena sativa and Linum usitatissimum, were instead species that are, or have formerly been, cultivated, suggesting that they profited from high propagule pressure, but are unlikely to cause negative impacts.

In accordance with earlier work we found that invasion success was highest for species with a long residence time, which was the third most important factor predicting alien species success [START_REF] Hamilton | Life-history correlates of plant invasiveness at regional and continental scales[END_REF][START_REF] Wilson | Residence time and potential range: crucial considerations in modelling plant invasions[END_REF][START_REF] Sorte | Extra-regional residence time as a correlate of plant invasiveness: European archaeophytes in North America[END_REF]Dawson et al. 2011a). For example, Veronica persica, an archeophyte in Europe, turned out to be among the most successful alien species. Residence time correlates with a number of factors such as propagule pressure or the overcoming of an initial lag-phase in the process of invasion [START_REF] Richardson | Plant invasions: merging the concepts of species invasiveness and community invasibility[END_REF][START_REF] Wilson | Residence time and potential range: crucial considerations in modelling plant invasions[END_REF], particularly by allowing time for acclimatization and local adaptation [START_REF] Lavergne | Increased genetic variation and evolutionary potential drive the success of an invasive grass[END_REF]. Conversely, we found no strong relationship between invasion success and introduction pathways. This is a surprising result given that most of the naturalized alien plants considered in this study (ca. 50%) have been introduced as ornamental or horticultural plants. These are often suggested to spread faster than other alien species [START_REF] Wilson | Residence time and potential range: crucial considerations in modelling plant invasions[END_REF]) due to selective breeding and hybridization in the cultivation process, high propagule pressure and repeated introduction of new genetic material through the ornamental trade [START_REF] Lavergne | Increased genetic variation and evolutionary potential drive the success of an invasive grass[END_REF]. However, our finding may result from repeated introductions diluting adaptations to new environments.

The most interesting results from our study are those linking invasion success to similarity patterns between plant invaders and native communities. Many studies have already demonstrated the existence of non-random similarity patterns in invasions [START_REF] Diez | Darwin's naturalization conundrum: dissecting taxonomic patterns of species invasions[END_REF][START_REF] Carboni | Darwin's naturalization hypothesis: scale matters in coastal plant communities[END_REF][START_REF] Pellock | Validity of Darwin's naturalization hypothesis relates to the stages of invasion[END_REF]Gallien et al. 2014b;[START_REF] Kraft | Functional trait and phylogenetic tests of community assembly across spatial scales in an Amazonian forest[END_REF].

However, here we showed that even when accounting for traits and introduction history, phylogenetic and functional similarity metrics were by far the most important correlates of plant invasion success in French grasslands.

In accordance with the filtering model of invasion biology (Richardson & Pysek 2006), whether successful aliens co-occurred with similar or dissimilar natives depended on the spatial scale (habitat vs. plot). Specifically, the most successful alien species (both widespread and locally abundant) were those sharing habitat preferences with close relatives and functionally similar species in the native species pool [START_REF] Cadotte | Phylogenetic relatedness and plant invader success across two spatial scales[END_REF][START_REF] Lemoine | Phylogenetic relatedness and leaf functional traits, not introduced status, influence community assembly[END_REF]. This finding suggests that alien species share ecological adaptations with natives when they successfully spread within a particular habitat type. It also implies that widespread, locally abundant and generalist aliens are not competitively limited at the habitat scale. The reason may be either that these very successful species are better competitors than natives, or that they avoid co-occurrence with closely related native competitors at fine scales. The fact that the most successful aliens tended to co-occur with functionally distinct natives at the plot scale (having high MDMCS.plot values, smaller seeds and lower height, Fig. 5) supports the latter hypothesis. These results corroborate previous findings showing that filtering mechanisms differ with the scale of observation [START_REF] Thuiller | Resolving Darwin's naturalization conundrum: a quest for evidence[END_REF][START_REF] Carboni | Darwin's naturalization hypothesis: scale matters in coastal plant communities[END_REF][START_REF] Kraft | Functional trait and phylogenetic tests of community assembly across spatial scales in an Amazonian forest[END_REF]). Here we find evidence of multiple hierarchical filters at an unprecedented spatial extent, thanks to the use of three dimensions of invasion success, rather than just based on dichotomous contrasts of aliens vs. natives or invasive vs. naturalized species.

Interestingly, we did not find strong support for intrinsic functional traits defining an overall successful invader [START_REF] Lemoine | Phylogenetic relatedness and leaf functional traits, not introduced status, influence community assembly[END_REF]. By contrast, all hierarchical indices emphasizing trait differences to the natives at the plot level were much more important than the trait values per se. In particular, successful aliens were shorter, had smaller seeds and greater SLA than natives co-occurring within plots. Small seed mass generally correlates with prolific seed output, wider dispersal (often through wind), long persistence in the seed bank, short juvenile periods and high relative growth rate of seedlings [START_REF] Moles | Seed size and plant strategy across the whole life cycle[END_REF][START_REF] Rejmanek | What attributes make some plant species more invasive?[END_REF][START_REF] Van Kleunen | A meta-analysis of trait differences between invasive and noninvasive plant species[END_REF]. These characteristics are likely to favor rapid colonization by alien species in sites where natives are not as fast-growing, particularly in response to grassland disturbance. Competition for light does not seem to be the main process filtering successful aliens in grasslands, as shorter species were more successful. Shorter species can spread laterally in disturbed habitats, thereby being superior competitors for space and/or soil nutrients, and are also likely to be favored in grazed habitats, due to both grazing avoidance and better ability to deal with higher soil bulk density associated with trampling. Therefore, our data suggest that in grasslands, faster growth rates than the surrounding natives, specifically at the juvenile stage (small seed size), together with lateral spread and good dispersal ability to reach disturbed patches are more important characteristics than high competitive ability for light interception (competitioncolonization trade-off; [START_REF] Kneitel | Trade-offs in community ecology: linking spatial scales and species coexistence[END_REF].

Based on the Rabinowitz framework, we also highlighted distinct ecological processes driving different components of invasion success. In contrast with the most successful alien species of class A (widespread, generalists and abundant), aliens of class B were widespread and generalists but remained locally sparse. Aliens of this latter group were very distantly related to the natives not only at the plot scale but also at the habitat scale. One possible explanation is that these species are poor competitors in most grasslands (e.g. they were mostly species with low height and SLA values). They may only colonize grassland sites when resources become available through resource pulses which weaken competitive exclusion [START_REF] Davis | Fluctuating resources in plant communities: a general theory of invasibility[END_REF] or due to continuous propagule pressure and dispersal from adjacent, more suitable habitats. Finally, the alien species in the rare categories are neither very similar nor dissimilar to co-occurring natives at the habitat scale, suggesting that rather neutral processes determine their presence in these grasslands. We argue that processes driving invasion are not the same across alien species and that specific processes are important for widespread species compared to less successful ones. Future studies should acknowledge that the relative importance of filtering processes affecting the success of alien species is likely to change along the classical introduction-naturalization-invasion continuum (sensu [START_REF] Richardson | Plant invasions: merging the concepts of species invasiveness and community invasibility[END_REF].

Finally, through a cross-validation procedure we highlighted that our model based on similarity to the natives also had the potential to estimate invasion success of newly introduced species. Predictive performance was generally reasonable (Appendix S6), though we found high prediction errors and under-predictions particularly for the invaders of class C. These species are locally abundant and generalist, but they are at the moment restricted to a narrow extent, potentially close to introduction sources (i.e. they may represent recent introductions). The invasion success of these species, which are likely not to be at equilibrium, is the most difficult to predict. However, it is encouraging that, while invasion success of species in class A was not always correctly identified, the most well-known invaders in this class (e.g O. biennis, E. annuus, B. inermis, J. tenuis) were indeed very well predicted by the model (Appendix S6). Such a framework has thus the potential to be used in other types of systems and in other regions.

Conclusions

Using French grasslands as a study system, we showed that complementary measures of invasion success are useful to identify the driving forces of invasion across spatial scales and environmental gradients. The Rabinowitz' components originally created for assessing different forms of species' rarity represent a straightforward framework for categorizing alien species and, in turn, for identifying the most successful and potentially harmful ones. In doing so, this framework provides guidance for differentiating and prioritizing alien species for management. Further, our use of complementary similarity metrics based on cooccurrence patterns allowed the identification of the relevant processes of invasion at two spatial scales. The most successful alien species are specifically adapted to the habitat they invade the most, but they then express different traits than co-occurring native species at fine scale within these habitats. Together with residence time and dispersal ability, these characteristics make up the winning alien species in permanent grasslands. Moreover we showed that these characteristics in a modeling framework are likely to give hints on the invasion success that can be attained by newly introduced species. Taken together, the set of approaches and metrics we have proposed here can be used in other ecosystems and biogeographic regions, not only to anticipate invasions but also to gain a better understanding of the processes shaping invasion success worldwide. We follow the same rationale to argue that alien species that fulfill all of these characteristics have high invasion success. At the opposite extreme, aliens in class H are unsuccessful. In between, alien species common for only one criterion have low success, while species common for at least two criteria have intermediate success. From the three criteria a continuous measure of invasion success can also be derived based on a dimension reduction analysis. Phylogenetic tree of alien species in 50,000 French grassland plots with the corresponding commonness class. Color scheme follows the one in the contingency table in Fig. 1. Variable importance in the Boosted Regression Trees-analysis. Bar colors indicate symmetric invasibility metrics (black), hierarchical competition metrics (dark grey), species traits (light grey) and introduction history (white). MDMCS.Habitat.phylo = Mean phylogenetic distance to the most often co-occurring species at habitat level; MDMCS.Habitat.func = Mean functional distance to the most often co-occurring species at habitat level; MDMCS.plot.func = Mean functional distance to the most often co-occurring species at plot level; MDMCS.plot.phylo = Mean phylogenetic distance to the most often co-occurring species at plot level; SM.hier.plot = mean (signed) difference in seed mass at plot level; Height.hier.plot = mean (signed) difference in height at plot level; SLA.hier.plot = mean (signed) difference in SLA at plot level; SM.hier.Hab = mean (signed) difference in seed mass at habitat level; Height.hier.Hab = mean (signed) difference in height at habitat level; SLA.hier.Hab = mean (signed) difference in SLA at habitat level. Also see Table 1.

Figure 1 .

 1 Figure 1. Conceptual visualization of the Rabinowitz framework. Shades of green indicate different environments, blue dots indicate a species' occurrence, and the size of the dots indicates its local abundance. Depending on their geographical distribution (narrow vs. widespread), local abundance (abundant vs. sparse) and niche breadth (specialist vs. generalist) species can be divided into eight distinct classes (A through H in the figure). Only species in class A, geographically widespread, locally abundant and generalists, can be considered common.We follow the same rationale to argue that alien species that fulfill all of these characteristics have high invasion success. At the opposite extreme, aliens in class H are unsuccessful. In between, alien species common for only one criterion have low success, while species common for at least two criteria have intermediate success. From the three criteria a continuous measure of invasion success can also be derived based on a dimension reduction analysis.

Figure 2 .

 2 Figure 2. Panel on the left: Numbers of alien species in the eight Rabinowitz classes. Red: common based on three criteria; Dark orange: common for two criteria; Light Orange: common for one criterion; Yellow: rare. Panel on the right: Concordance between the continuous invasion success measure (based on the first axis of a PCA) and the Rabinowitz classes for the 127 alien species. The least widespread class is on the left of the graph (class H -small regional distribution, locally not abundant, small niche breadth) and the most successful aliens are on the right (class A -regionally widespread, locally abundant, large niche breadth).

  Figure 3.

  Figure 4.

Figure 5 .

 5 Figure 5. Response curves of invasion success to the five most important variables in the BRT. All variables are rescaled between zero and one (apart for the hierarchical indices to highlight positive vs. negative trait differences). Invasion success is expressed by the first axis of a PCA. Negative values relate to low invasion success while positive values relate to high invasion success. Legend: MDMCS.hab.phylo = Mean phylogenetic distance to the most often co-occurring species within habitats MDMCS.hab.func = Mean functional distance to the most often co-occurring species within habitats; MDMCS.plot.func: Mean functional

Ecol Lett. Author manuscript; available in PMC 2016 August 03.

Acknowledgements

We thank J. Renaud for his help with the preparation of the environmental data, and R. Douzet and three anonymous reviewers for their insightful comments. The research leading to these results received funding from the European Research Council under the European Community's Seven Framework Programme FP7/2007-2013 Grant Agreement no. 281422 (TEEMBIO). MC was supported by a Marie Curie Intra European Fellowship within the European Community's Seventh Framework Program (FP7-PEOPLE-2012-IEF; IASIMOV 327616). MC, WT, SL, TM and F.E. acknowledge support by the ERA-Net BiodivERsA, with the national funders ANR (ANR-13-EBID-0004), DFG and FWF, part of the 2013 BiodivERsA call (project "WhoIsNext"). The LECA is part of Labex OSUG@2020 (ANR10 LABX56). CV was supported by the European Research Council (ERC) Starting Grant Project "Ecophysiological and biophysical constraints on domestication in crop plants" (Grant ERC-StG-2014-639706-CONSTRAINTS). All data were extracted from the DIVGRASS database supported by the French Foundation for Research on Biodiversity (FRB; www.fondationbiodiversite.fr) in the context of the CESAB project "Assembling, analysing and sharing data on plant functional diversity to understand the effects of biodiversity on ecosystem functioning: a case study with French Permanent Grasslands" (DIVGRASS).

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.

Table 1

Summary of distance-based invasibility metrics used in this study, identifying the type of distance (Functional vs. phylogenetic), the spatial scale of reference and the underlying assumptions. The hierarchical indices are reported only for the plot scale for brevity.

Index