Keywords: phylogenetic diversity area-relationship, functional diversity area-relationship, species arearelationship, conservation biogeography, diversity indices, Hill's numbers, mammals

Aim-To define biome-scale hotspots of phylogenetic and functional mammalian biodiversity (PD and FD, respectively) and compare them to 'classical' hotspots based on species richness (SR) only.

Location-Global

Methods-SR, PD & FD were computed for 782 terrestrial ecoregions using distribution ranges of 4616 mammalian species. We used a set of comprehensive diversity indices unified by a recent framework that incorporates the species relative coverage in each ecoregion. We build large-scale multifaceted diversity-area relationships to rank ecoregions according to their levels of biodiversity while accounting for the effect of area on each diversity facet. Finally we defined hotspots as the top-ranked ecoregions.

Results-While ignoring species relative coverage led to a relative good congruence between biome top ranked SR, PD and FD hotspots, ecoregions harboring a rich and abundantly represented evolutionary history and functional diversity did not match with top ranked ecoregions defined by species richness. More importantly PD and FD hotspots showed important spatial mismatches. We also found that FD and PD generally reached their maximum values faster than species richness as a function of area.

Main conclusions-The fact that PD/FD reach faster their maximal value than SR may suggest that the two former facets might be less vulnerable to habitat loss than the latter. While this point is expected, it is the first time that it is quantified at global scale and should have important consequences in conservation. Incorporating species relative coverage into the delineation of multifaceted hotspots of diversity lead to weak congruence between SR, PD and FD hotspots. This means that maximizing species number may fail at preserving those nodes (in the phylogenetic or functional tree) that are relatively abundant in the ecoregion. As a consequence it may be of prime importance to adopt a multifaceted biodiversity perspective to inform conservation strategies at global scale.

Introduction

Understanding the ecological and evolutionary processes driving the distribution of life on Earth is essential from both applied and theoretical perspectives. The quantification of biodiversity, central to conservation science, has recently moved from a pure species-count focus (e.g. species richness, SR) to a more integrative approach. Biodiversity assessments now consider the overall evolutionary history embedded within a set of taxa (i.e. phylogenetic diversity, PD) along with the diversity of ecological traits (i.e. functional diversity, FD). In conservation science, this novel approach has redefined the identification of species of conservation interest by taking their high evolutionary or functional distinctiveness into consideration [START_REF] Mouillot | Rare species support vulnerable functions in high-diversity ecosystems[END_REF][START_REF] Isaac | Mammals on the EDGE: conservation priorities based on threat and phylogeny[END_REF] and has also made it possible to detect unique macroecological assemblages [START_REF] Forest | Preserving the evolutionary potential of floras in biodiversity hotspots[END_REF], for example, 'cradles' and 'museums' of life [START_REF] Chown | Areas, cradles and museums: the latitudinal gradient in species richness[END_REF]. Furthermore, the loss of FD or PD per unit of habitat loss is likely to be a better predictor of ecosystem vulnerability than the loss of single species. Indeed, the loss of a given amount of FD or PD, often assumed to be related to particular combinations of functional traits or of a certain lineage respectively, may threaten the functioning of the ecosystem, whereas the loss of a given single species might not be noticeable if redundant species persist [START_REF] Loreau | Biodiversity and ecosystem functioning: synthesis and perspectives[END_REF][START_REF] Srivastava | Phylogenetic diversity and the functioning of ecosystems[END_REF].

This new perspective also provides fundamental insights into community assembly at multiple spatial scales [START_REF] Mouquet | Ecophylogenetics: advances and perspectives[END_REF]. A multifaceted approach may help unravel the different drivers of community structure (e.g. competition or environmental filtering [START_REF] Webb | Phylogenies and Community Ecology[END_REF]) or ecosystem functioning [START_REF] Cadotte | Using phylogenetic, functional and trait diversity to understand patterns of plant community productivity[END_REF][START_REF] Gravel | Phylogenetic constraints on ecosystem functioning[END_REF]. In macroecology, contrasting SR, PD, and FD offers a potential means for disentangling the processes shaping large-scale diversity distribution [START_REF] Safi | Understanding global patterns of mammalian functional and phylogenetic diversity[END_REF][START_REF] Davies | Phylogenetic diversity as a window into the evolutionary and biogeographic histories of present-day richness gradients for mammals[END_REF][START_REF] Huang | Traits, trees and taxa: global dimensions of biodiversity in mammals[END_REF]. For example, the global latitudinal diversity gradient has recently been reinterpreted from a novel evolutionary perspective, merging earth climatic history, phylogenetic diversity and species richness in a unified and testable framework [START_REF] Hawkins | Different evolutionary histories underlie congruent species richness gradients of birds and mammals[END_REF]. A multifaceted perspective thus represents a promising avenue for exploring diversity distribution because it is at the crossroads between ecology, evolution and conservation biology but also paleontology and paleoclimatology [START_REF] Hawkins | Post-Eocene climate change, niche conservatism, and the latitudinal diversity gradient of New World birds[END_REF].

One of the most striking features of biodiversity is the spatial heterogeneity of its distribution, with some regions harboring extraordinary levels of biodiversity: the so-called hotspots of biodiversity [START_REF] Reid | Biodiversity hotspots[END_REF][START_REF] Ceballos | Global mammal distributions, biodiversity hotspots, and conservation[END_REF][START_REF] Guilhaumon | Taxonomic and regional uncertainty in speciesarea relationships and the identification of richness hotspots[END_REF]. They have not only fascinated macroecologists, who try to understand their origins (e.g. historical perspective [START_REF] Wiens | Phylogenetic origins of local-scale diversity patterns and the causes of Amazonian megadiversity[END_REF]), but also conservationists seeking the best opportunities to allocate the limited resources available for global scale conservation. For example the "biodiversity hotspots" concept has been proposed to prevent the extinction of large numbers of endangered species, by protecting places "where exceptional concentrations of endemic species are undergoing exceptional loss of habitat" [START_REF] Myers | Biodiversity hotspots for conservation priorities[END_REF].

The most recent comparisons of the worldwide distribution of hotspots have been limited to different taxonomic groups and components of species richness for a given taxa (e.g. endemic, total, endangered [START_REF] Orme | Global hotspots of species richness are not congruent with endemism or threat[END_REF][START_REF] Ceballos | Global mammal distributions, biodiversity hotspots, and conservation[END_REF][START_REF] Lamoreux | Global tests of biodiversity concordance and the importance of endemism[END_REF]) or when carried out in a multifaceted context, have included only limited functional information (e.g. [START_REF] Huang | Traits, trees and taxa: global dimensions of biodiversity in mammals[END_REF] only used geographic range size and body mass as the sole descriptors of mammals functional diversity to define hotspots). This lack of relevant trait information makes it difficult to adequately represent spatial distribution of functional diversity because geographic range size may not portray properly species niches, rather it is mostly influenced by historical biogeography and macroevolution [START_REF] Gaston | The Structure and Dynamics of Geographic Ranges[END_REF].

Here we conducted an identification of global hotspots of mammal taxonomic, phylogenetic and functional mammalian diversity. We based our analyses on the updated version [START_REF] Fritz | Geographical variation in predictors of mammalian extinction risk: big is bad, but only in the tropics[END_REF] of the dated phylogeny of [START_REF] Bininda-Emonds | The delayed rise of present-day mammals[END_REF] and a set of functional traits encompassing important aspects of mammal resource use, selected to represent independent and informative niche dimensions [START_REF] Safi | Understanding global patterns of mammalian functional and phylogenetic diversity[END_REF]. We used the world's ecoregions [START_REF] Olson | Terrestrial Ecoregions of the World: A New Map of Life on Earth[END_REF] to define geographical units harboring unique species assemblages and ecosystems. Ecoregions have proven valuable for addressing a range of questions in macroecology and more applied conservation issues [START_REF] Guilhaumon | Taxonomic and regional uncertainty in speciesarea relationships and the identification of richness hotspots[END_REF][START_REF] Lamoreux | Global tests of biodiversity concordance and the importance of endemism[END_REF].

To account for expected area-effects on taxonomic [START_REF] Triantis | The island species-area relationship: biology and statistics[END_REF], phylogenetic [START_REF] Morlon | Spatial patterns of phylogenetic diversity[END_REF] and functional diversity [START_REF] Cumming | Contrasting spatial patterns of taxonomic and functional richness offer insights into potential loss of ecosystem services[END_REF] we constructed diversity-area relationships (DARs hereafter) for 13 terrestrial biomes (global scale regions gathering ecoregions experiencing similar environmental conditions such as tundras or mediterranean forests). We used a model-averaging approach that fits 19 mathematical functions to the data [START_REF] Triantis | The island species-area relationship: biology and statistics[END_REF][START_REF] Guilhaumon | Taxonomic and regional uncertainty in speciesarea relationships and the identification of richness hotspots[END_REF] and then computed an AIC-weighted average of the 19 predicted curves. To quantify the different types of diversity, we used a set of unified taxonomic, phylogenetic and functional diversity indices weighting differently species' coverage [START_REF] Chao | Phylogenetic diversity measures based on Hill numbers[END_REF] and corresponding to modified version of Faith [START_REF] Faith | Conservation evaluation and phylogenetic diversity[END_REF]), Allen [START_REF] Allen | A new phylogenetic diversity measure generalizing the shannon index and its application to phyllostomid bats[END_REF] and Rao [START_REF] Rao | Diversity and dissimilarity coefficients: a unified approach[END_REF] quadratic entropy (see Methods). For each diversity indice we identified as hotspots those ecoregions with the largest positive deviations from respectively SARs (species-area relationships), PDARs (phylogenetic diversity-area relationships) and FDARs (functional diversity-area relationships) and investigated their spatial congruences.

Our global exploration of mammals SARs, PDARs and FDARs reveals important mismatches between the spatial scaling and the geographical extremes of SR, PD and FD, calling for integrative approaches. Methods 1. Dataset 1.1. Mammal assemblages-We used the distribution maps provided by the Mammal Red List Assessment (http://www.iucnredlist.org/) for 4616 terrestrial species (for which we have functional traits, see below) to obtain occurrence data for each of the 827 ecoregions defined by [START_REF] Olson | Terrestrial Ecoregions of the World: A New Map of Life on Earth[END_REF]. We retained 782 ecoregions (mean number of ecoregions per biome = 60.1; standard deviation= 53.3; min=17; max=223). Ecoregions are a valuable tool for studying multifaceted hotspots because it also serves as the basis of World Wildlife Fund's conservation planning [START_REF] Olson | The Global 200: a representation approach to conserving the Earth's most biologically valuable ecoregions[END_REF], the Nature Conservancy's international efforts [START_REF] Groves | Drafting a conservation blueprint: a practitioner's guide to planning for biodiversity[END_REF], and the delineation of Conservation International's Biodiversity Hotspots [START_REF] Mittermeier | Hotspots Revisited. Cemex[END_REF] and High Biodiversity Wilderness Areas [START_REF] Mittermeier | Wilderness and biodiversity conservation[END_REF]. Furthermore, it has commonly been used to define taxonomic hotspots [START_REF] Lamoreux | Global tests of biodiversity concordance and the importance of endemism[END_REF][START_REF] Guilhaumon | Taxonomic and regional uncertainty in speciesarea relationships and the identification of richness hotspots[END_REF] because it encompasses relatively homogeneous biological systems. We retained ecoregions harboring more than one mammal species and excluded Mangrove ecoregions and large uninhabited parts of Greenland and Antarctica because of low data reliability or availability for these areas [START_REF] Lamoreux | Global tests of biodiversity concordance and the importance of endemism[END_REF]. Domestic mammals were also excluded from the analysis.

For each ecoregion and species, species coverage (C i ) was calculated as the intersected surface (in km 2 ) between the range of species and the ecoregion. We then computed, for each species i, the following relative coverage (RC i , equation 1).

Equation (1)

Basically a species will have low relative coverage in a given ecoregion if its distribution range is small. The relative coverage was used to calculate diversity indices incorporating relative abundance (see below). By doing so, we were able to differentiate a species poorly represented in an ecoregion, but with a unique evolutionary history (e.g. a monotreme species) or with unique functional traits (e.g a top predator), from species with similar evolutionary history (or functional traits), but with greater occupancy in the ecoregion. This weighting scheme emphasizes species that are well distributed in the ecoregion. Establishing how our measure of species relative coverage is important for conservation and ecosystem functioning is not straightforward. Nevertheless we believe that the evolutionary history/ functional characteristics of a species that show a very small distribution range in a given ecoregion should not have the same theoretical influence on PD/FD than a widespread species in this ecoregion. Although it is unlikely that our measure of relative coverage represents a direct measure of local species abundance, it has been shown that a positive relationship between range size and local abundance is common [START_REF] Gaston | Abundanceoccupancy relationships[END_REF]. Nevertheless departure from this relationship probably exists. First, we did not use the complete range size of the species but only its extent in the ecoregion. Second, we acknowledge that the potential important residual variation that exists around the relationship may depend on species life history traits. For example species with high dispersal abilities (or a species at high level in the trophic hierarchy) may have a large range size but be relatively rare at the local scale. It is also possible that a species with a narrow range may exhibit high local abundances, for example because it uses an abundant resource that is restricted to a small area of the ecoregion. Nevertheless, we believed that our measure of species coverage was a needed first step to incorporate abundances in the definition of PD/FD hotspots.

1.2. Phylogeny and functional traits-We used the calibrated and dated ultrametric phylogenetic tree updated by [START_REF] Fritz | Geographical variation in predictors of mammalian extinction risk: big is bad, but only in the tropics[END_REF] from [START_REF] Bininda-Emonds | The delayed rise of present-day mammals[END_REF].

We computed functional diversity indices using body mass (log-transformed), diet (vertebrates, invertebrates, foliage, stems and bark, grass, fruits, seeds, flowers, nectar and pollen, roots and tubers), habits (aquatic, fossorial, ground dwelling, above ground dwelling), activity period (diurnal, nocturnal, cathemeral, crepuscular), and litter size (data from [START_REF] Safi | Understanding global patterns of mammalian functional and phylogenetic diversity[END_REF]). These traits encompass important aspects of mammal resource use, including the temporal and spatial windows used to get their food. They represent independent and informative niche dimensions for evaluating variability in mammal traits related to important ecosystem processes, such as decomposition and seed dispersal, as well as trophic control (Sekercioglu, 2010).

Diversity indices

A myriad of indices have been proposed in the last years to include species traits in diversity indices [START_REF] Pavoine | Measuring biodiversity to explain community assembly: a unified approach[END_REF]. Here we follow the comprehensive framework from [START_REF] Chao | Phylogenetic diversity measures based on Hill numbers[END_REF], which unifies a set of taxonomic and phylogenetic/functional diversity indices based on Hill numbers. There were three reasons for this: first it unifies most of the TD, PD and FD indices used in the literature (see below and table 1); second it represents equivalent numbers of species to satisfy the replication principle that ensures intuitive results for ecologists and conservation biologists [START_REF] Jost | Entropy and diversity[END_REF][START_REF] Chao | Phylogenetic diversity measures based on Hill numbers[END_REF]. For example if the PD of an assemblage equals d (d being a real positive number), it has the same diversity as a community consisting of d equally abundant and maximally distinct species (i.e. with the maximum distance observed in the phylogeny). Third, we present here one of the only comprehensive and intuitive framework that incorporates species relative coverage (or abundance) into biodiversity indices.

2.1.Phylogenetic diversity-Consider a phylogenetic tree composed of a set B t of i branches. PD can be defined the 'mean diversity of order q over T years' as [START_REF] Chao | Phylogenetic diversity measures based on Hill numbers[END_REF]:

Equation (2)
Where L i is the length of branch i in the set B t , a i is the total abundance descended from branch i (i.e the summed abundance or relative coverage of species descending from this branch), T is the height of the tree. The parameter q impacts the influence of node (or branch segment) abundance on the diversity index: a high q value gives more weight to the nodes with high relative abundances. This general formula encompasses a set of well-known diversity indices. With q=0, Faith cor PD = PDfaith /T, PDfaith being the phylogenetic diversity defined by [START_REF] Faith | Conservation evaluation and phylogenetic diversity[END_REF] and Faith cor PD being the corrected version of Faith PD. With q=1, Allen cor PD = exp(H p /T), H p being the phylogenetic entropy as defined by [START_REF] Allen | A new phylogenetic diversity measure generalizing the shannon index and its application to phyllostomid bats[END_REF] and Allen cor PD being the corrected version of H p .With q=2, Rao cor PD = 1/ (1-QE), QE being the quadratric entropy defined by [START_REF] Rao | Diversity and dissimilarity coefficients: a unified approach[END_REF] and Rao cor PD being the corrected version of QE (see table 1 for details). To summarize, q influences the relative weight of widespread versus rare species into the computation of the diversity index. It gives progressively more weight to widespread species and progressively ignores rare species. This point could be problematic if a species is rare and endemic in this ecoregion because we will progressively ignore this unique species. Nevertheless our study aims at characterize the evolutionary history and the functional characteristics that are widespread in a given ecoregion (and somehow representative), which justifies the use of Chao's framework.

2.2.Functional diversity-

We adapted previous indices for functional diversity. First, we calculated functional distance among pairs of species using Gower distance, which can mix categorical and continuous traits with equal weight and can cope with missing values (some traits were missing for 80 species, representing less than 2% of our dataset). We then applied a hierarchical cluster algorithm to convert the functional distance matrix into a functional dendrogram ensuring the ultrametric property (note that using non ultrametric functional distances did not change our conclusions) using the UPGMA method (function hclust in R (R Development Core Team, 2010)). The corresponding FD indices were named Faith cor FD, Allen cor FD and Rao cor FD (see table 1). Note that Faith cor FD is equivalent to the [START_REF] Petchey | Functional diversity: back to basics and looking forward[END_REF] definition of FD (i.e. 'the total branch length of a functional dendrogram'). Like Faith cor PD, Faith cor FD is intrinsically correlated to SR [START_REF] Huang | Traits, trees and taxa: global dimensions of biodiversity in mammals[END_REF]. This is the case for all the dendrogram-based approach to estimate functional volume. Nevertheless it is interesting to use it here because it is directly comparable to Faith cor PD and represent a diversity volume (or "richness" sensu [START_REF] Pavoine | Measuring biodiversity to explain community assembly: a unified approach[END_REF]). In addition we computed FD using bodymass only to test to which extent the use of multiple traits influence our results. Note also that the expected correlation between SR and FD/PD based on dendrogram lengths become weaker when moving from q=0 to q=2.

2.3. Species diversity-For maximally distinct species (i.e. a star phylogeny or star functional dendrogram), these indices actually constitute species diversity indices, namely species richness when q=0, the exponential of Shannon entropy when q=1, and the inverse of Simpson diversity when q=2 (see table 1, [START_REF] Chao | Phylogenetic diversity measures based on Hill numbers[END_REF]. We used these taxonomic diversity indices to compare appropriately PDAR/FDAR with corresponding SAR (i.e. comparing DARs that are built with diversity indices based on the same q). We, however, always compared PD/FD hotspots with those based on SR (and not Simpson or Shannon indices) since hotspots list based only on taxonomic diversity indices (i.e only quantifying evenness in abundances) might not be appropriate in a conservation context.

Constructing DARs

To account for expected area-effects on SR, PD and FD, we provide a construction of DARs for 13 terrestrial biomes [START_REF] Olson | Terrestrial Ecoregions of the World: A New Map of Life on Earth[END_REF]. Such DARs correspond to a non-overlapping design (i.e. they are built from single data points, which corresponds to a type IV curve in Scheiner's terminology ( 2003)).

Set of models-A wide range of statistical models have been proposed to describe

SARs [START_REF] Tjørve | Shapes and functions of species-area curves (II): a review of new models and parameterizations[END_REF]. Here, 19 models were selected to fit SAR, PDAR, and FDAR following [START_REF] Triantis | The island species-area relationship: biology and statistics[END_REF] see Appendix S1 in Supporting information). Recent attempts to model PDAR only used the power model [START_REF] Morlon | Spatial patterns of phylogenetic diversity[END_REF] but, given the uncertainty regarding the shape of PDAR and FDAR, we tested a large spectrum of models. These models were chosen because they vary in forms (e.g. sigmoid or convex, including asymptotic relationships) and complexity (2 to 4 parameters).

Model fitting-

We constructed 117 datasets (9 indices * 13 biomes) and fitted 19 models to each dataset, for a total of 2223 DARs. We carried out our analyses using another dataset that also adds an 'artificial' point of null diversity and null area (0.001 and 0.001 to avoid computing problems). Models were fitted using nonlinear regression with minimization of the residual sum of squares. Models were further evaluated by examining the normality and homoscedasticity of residuals. To do so, we applied the Lilliefors's test for normality and a Pearson correlation between squared residuals and area for homoscedasticity. Previous studies [START_REF] Guilhaumon | Taxonomic and regional uncertainty in speciesarea relationships and the identification of richness hotspots[END_REF] considered a model valid when the p-value associated with the normality and homoscedasticity tests exceeded the arbitrary threshold of 5%. All DAR analyses were carried out using an updated version of the 'mmSAR' package [START_REF] Guilhaumon | Taxonomic and regional uncertainty in speciesarea relationships and the identification of richness hotspots[END_REF] for R statistical and programming environment (R Development Core Team, 2010).

3.3.

Model averaging-For each dataset, we discriminated between different models using an information-theory framework designed to evaluate multiple working hypotheses [START_REF] Burnham | Model selection and inference: a practical information-theoretic approach[END_REF]. The Akaike Information Criterion (AIC) can be used to evaluate the goodness of fit of different non-nested models on a given dataset. The weights of evidence were then derived from the AIC values to evaluate the relative likelihood of each model given the data and the set of models [START_REF] Burnham | Model selection and inference: a practical information-theoretic approach[END_REF]. Using these weights we derived averaged DARs for each biome and each diversity indice.

Model standardization and comparison-

For each DAR, we divided each predicted diversity value by that of the largest ecoregion in the considered biome in order to report the percentage of maximal diversity reached in this largest ecoregion. The resulting standardized DAR therefore range between 0 and 1, which makes it possible to compare the scaling of different diversity facets with area (see below).

Hotspot lists and spatial congruence

4.1. Hotspot selection-Averaged residuals were calculated from the standardized averaged model (as defined in point 3). A positive residual for a given ecoregion means that observed diversity is higher than expected given its area. Hotspots were defined as the ecoregions with the highest residuals. We ranked ecoregions according to their averaged residuals: the higher the residuals, the higher the concentration of biodiversity in the ecoregion. Note that ranking in terms of original or standardized curve/observed diversity gives exactly the same results because standardization is linear (see Appendix S2). We also derived an averaged rank across SR, PD and FD hotspots to provide an integrative hotspot definition by summing up the ranking for each ecoregion across the biodiversity facets (i.e. SR, PD and FD).

Impact of DAR shape on hotspot lists-We investigated whether PDARs and

FDARs were different enough from corresponding SAR to deeply modify the hotspot rankings. In other words we wanted to test if PDARs/FDARs are needed to define hotspots or if SAR is a good proxy for FDAR/PDAR when defining hotspots. SAR, PDAR and FDAR were directly comparable thanks to the standardization procedure explained above (they are all expressed as a proportion of the maximal diversity predicted for the largest ecoregion and thus vary between 0 and 100%). We computed the difference between the standardized PD/FD in each ecoregion and the proportion of diversity predicted by the area using the SAR (and not PDAR/FDAR as previously done) and ranked these differences to compute hotspot lists. Then, we compared the congruence between PD/FD hotspot lists derived from SAR and the 'natural' PD/FD hotspot lists derived from the PDAR/FDAR (as explained in the point 4.1). If SARs correctly model the scaling of PD/FD with area, the lists of hotspots should be very similar. In this case, SARs would be well suited to directly modeling the spatial scaling of PD/FD to define hotspots and it would not be necessary to construct explicit PDAR/FDAR.

Results

We start first by reporting the general results of the statistical procedures related to the DAR estimations and then by describing the outcomes of this procedure on the hotspot lists.

DAR modeling

1.1. Convergence, homoscedasticity, and normality-One of the 19 models showed unrealistic fits (Epm2, see Appendix S1) and was not considered in the analysis. Of the remaining 2106 fits (13 biomes × 18 models × 9 indices), 1895 (90%) converged. Amongst the different indices, models fitting Rao-based DAR showed the highest nonnormality of residuals (homoscedasticity was not the limiting factor, results not shown). Indeed at the 1% level test of homoscedasticity and normality of the residuals, 53%, 68% and 75% of the models were valid for the Rao, Allen and Faith-based indices, respectively.

Relative model fit-

Model shape-To

illustrate the difference between the rate of increase in SAR and FDAR/PDAR, we plotted the difference between predicted PDAR/FDAR and their corresponding predicted SAR for four biomes that cover the latitudinal gradient (Fig. 1 and Appendix S5). The starting value of the curve was zero in most cases, while differences between PDAR/FDAR tended to zero as area increased. This means that PDAR/FDAR and SAR have the same proportion of diversity when area tends to zero (generally it was 0% of maximum diversity) and also ends at the same point because of the standardization (their respective maximum 100%). At the intermediate area between the two extremes, PDAR and FDAR were in general higher than SAR (i.e. positive difference), indicating that PDAR and FDAR reached their maximum diversity faster than SAR. This difference increased with the q parameter defining the weight given to species coverage in the diversity indices (Faith cor to Rao cor ). PDAR and FDAR had a similar shape in most cases. Results were qualitatively equivalent when fitting DARs without artificial zeros except when area tended to zero: PDAR/FDAR started at a relatively higher percentage of diversity than SARs and thus the difference between PDAR/FDAR and SAR curves tended to start with positive values for some biomes (see Appendix S5).

Functional and phylogenetic hotspots

We extracted residuals (i.e. observed minus predicted diversity) from each averaged DAR and ranked them to define hotspots of diversity. As an example, we mapped diversity ranks for tropical moist forests (Fig. 2 but also see Appendix S6 for all biomes and indices) considering Allen cor PD and FD hotspots as well as the traditional SR hotspots. SR rankings were relatively well distributed in the three continents whereas PD Allen cor hotspots were much more concentrated in Afrotropics (and Central America) or in Afrotropics and Indomalaysian realms for FD Allen cor hotspots. Interestingly, when focusing at the five hottest hotspots for this biome (Table 2), two important results emerged: (1) SR hotspots list shared few ecoregions with the lists of Allen cor PD, FD and integrative hottest hotspots (i.e. two, one and three ecoregions respectively); and (2) the FD and PD hottest hotspots shared only two ecoregions.

The same hotspot mismatches were found across all biomes (Fig. 3). For example, with the cut-off point for defining hotspot set at 5% richest ecoregions we found that congruence ranged from 5% (FD Rao cor hotspots vs SR hotspots) to 74% (PD faith cor hotspots vs. SR hotspots. Interestingly, when compared to hotspots defined with SR, hotspots defined with the Faith cor index strongly matched, while those defined with Rao cor strongly mismatched, those defined with Allen cor falling in-between. In other words, the hotspot rankings were significantly correlated -but not equal-across different indices (see Appendix S7). These differences were robust against the threshold used to define valid DAR models (i.e. the pvalue threshold used to reject a model based on the non normality and/or homoscedasticity of its residuals, see Appendix S8) due to the weak influence of this threshold on the definition of hotspots (see Appendix S9). We also explored to what extent the use of multiple traits influence the definition of hotspots. We show that FD hotspots lists based on body mass only differ from FD hotspots defined with all our complete set of traits (appendix S9).

On average, using SAR instead of PDAR/FDAR to define PD/FD hotspots marginally modified the hotspot list (Appendix S10). However, it turns out that there is still high variability between biomes. For some of these, using SAR instead of PDAR/FDAR dramatically changes the hotspot lists. Note also that the nonlinear fit of the power model alone gave relatively similar results to those obtained when using model averaging to define hotspots (Appendix S11).

Discussion

We found considerable geographical mismatches among global mammal hotspots of SR, PD and FD and, quite importantly, found that the magnitude of the mismatches depends on the index considered, which highlight the importance of considering a variety of indices [START_REF] Huang | Traits, trees and taxa: global dimensions of biodiversity in mammals[END_REF]. Mismatches were higher with Rao-based indices (Rao cor ), lower when using the Faith cor indices and in-between for the Allen cor indices, whatever the facet considered. This is not entirely surprising given the correlation between SR and PD/FD (high with Faith cor , medium with Allen cor and weak with Rao cor ). [START_REF] Rodrigues | Complete, accurate, mammalian phylogenies aid conservation planning, but not much[END_REF] already pointed out a high congruence between Faith cor PD hotspots and species richness hotspots. As a result, they concluded that incorporating phylogenetic information is not a major concern in conservation. Nevertheless, incorporating relative species coverage into the definition of multifaceted hotspots alters this conclusion. Faith cor indices do not incorporate species abundance or coverage and give equal weight to rare and dominant evolutionary history in a given location [START_REF] Chao | Phylogenetic diversity measures based on Hill numbers[END_REF]. However, it seems appropriate to give less weight to particular evolutionary histories (i.e. particular branch paths) that are rare in a given ecoregion because they are less representative than a widespread species in this ecoregion. Allen cor and Rao cor indices give more weight to a given branch if it is long and well represented in an ecoregion. For instance, hotspots defined using the PD Rao cor are mostly concentrated in the Australasian realm because of the presence of marsupials. This group has a unique evolutionary history since they diverged 147 millions years ago from placentals (extant eutherians, grouping the majority of mammals) and are widely distributed (i.e. large coverage) through the Australasian ecoregions (and not in South American ecoregions). These results are congruent with those of a recent study revealing important mismatches between global mammal trait variance hotspots and SR [START_REF] Huang | Traits, trees and taxa: global dimensions of biodiversity in mammals[END_REF]. Although these authors used a different approach (they used grid cells as geographical units and did not use information about species coverage) these close results are probably explained by the fact that Rao cor indices are linked to a measure of variance [START_REF] Pavoine | Measuring biodiversity to explain community assembly: a unified approach[END_REF]. We also showed that the use of body mass alone to define FD hotspots is not sufficient to match the FD hotspots defined with our complete set of traits but that it still represent an acceptable approximation. We also showed that PD and FD hotspots are not always congruent, suggesting that PD is not necessarily a good surrogate for FD (at least for the functional traits selected here).

As well as defining hotspots, DARs have been shown to be useful in both applied and fundamental ecology. We found that Faith cor PDAR and FDAR generally reach their maximum faster than SAR [START_REF] Cumming | Contrasting spatial patterns of taxonomic and functional richness offer insights into potential loss of ecosystem services[END_REF]. This result was expected since Faith cor PD and FD explicitly account for redundancy between species while SR does not. More specifically, it is possible that small sample areas already contain a broad set of phylogenetic history and functional diversity (e.g. a mouse and an elephant), whereas large sample areas contain relatively more redundant species (e.g. several species of mice) and thus PDAR/FDAR reach their maximum faster than SAR. [START_REF] Morlon | Spatial patterns of phylogenetic diversity[END_REF] obtained a similar result for PDAR on nested Mediterranean plant communities ranging from 6.25m 2 to 400m 2 of spatial extent. They used power law (see Appendix S1) to model PDAR and SAR and found that the rate of increase in Faith PD with area (z PDAR ) was slower that in SR (z SAR ). When standardizing DARs, PDAR is above SAR if zPDAR < z SAR . They showed that protected areas in Australian Mediterranean-like regions (representing 13% of the regions) capture 72% of PD, but only 56% of species richness, indicating that PDAR accumulates total diversity faster than SAR.

Our results show that if only a fraction of the total biome area is protected, the percentage of remaining PD (compared to the initial PD) will be higher than the percentage of remaining species. If we consider that PD or FD are better predictors of ecosystem functioning, resistance or resilience [START_REF] Cadotte | Using phylogenetic, functional and trait diversity to understand patterns of plant community productivity[END_REF][START_REF] Gravel | Phylogenetic constraints on ecosystem functioning[END_REF] than SR, it means that ecosystem features might be more robust to species loss than previously predicted (but see [START_REF] Mouillot | Rare species support vulnerable functions in high-diversity ecosystems[END_REF]).

We also found that a key feature of a comprehensive measure of diversity is that when rarely represented evolutionary history is progressively removed (i.e. using different values of q), the differences between PDAR/FDAR and SAR increase. In other words, PD/FD of abundant lineages reaches its maximum faster than when considering all lineages having the same coverage. This result suggests that the evolutionary history or functional traits of wellrepresented taxa are relatively more rapidly sampled when area increases. For example, branches of major mammal lineages (e.g. bats, rodents or carnivores) are probably already well sampled in small ecoregions and thus PD or FD reaches their maximum faster than TD in larger sample areas. It follows that well represented functional/phylogenetic biodiversity might be robust to habitat loss, a point that is not detected when considering all lineages having the same coverage.

Although SARs have been thoroughly investigated [START_REF] Scheiner | Six types of species-area curves[END_REF], we showed that there is not a single best model that fits all the data. Thus the automatic use of a single model (traditionally the linear version of power model) is not justified. Conversely, to date PDAR and FDAR have been subject to very little investigation (but see [START_REF] Cumming | Contrasting spatial patterns of taxonomic and functional richness offer insights into potential loss of ecosystem services[END_REF][START_REF] Morlon | Spatial patterns of phylogenetic diversity[END_REF][START_REF] Helmus | Phylogenetic diversity-area curves[END_REF][START_REF] Wang | Spatial patterns of tree species richness in two temperate forests[END_REF]. Here, given the important variability across biomes and indices, we also show that a single best model does not exist for PDAR and FDAR. Nevertheless the model-averaging framework allows taking into account these uncertainties and we used an averaged prediction to remove the area effect on PD/FD. We also asked whether the averaged SAR could be a good proxy of the averaged PDAR/FADR to remove this area effect to define PD/FD hotspots. We demonstrated that there is a notable difference between PD/FD hotspot lists defined using PDAR/FDAR and those defined using SAR, suggesting that the construction of PDAR/FDAR is required to define functionally-or phylogenetically-based hotspots and that SAR alone can not be used for this purpose.

We constructed DARs using a particular experimental design [START_REF] Scheiner | Six types of species-area curves[END_REF] but we are aware that all methodologies for constructing DARs have their own drawbacks and we suggest that the next challenge in the study of large scale multifaceted DARs is to test different methodological designs. For example the strictly nested design (SNQ) of [START_REF] Storch | Universal species-area and endemics-area relationships at continental scales[END_REF] seems particularly interesting to analyse. Nevertheless, since our work was about delineating hotspots of diversity, we had to construct DARs using a non-overlapping design.

Conclusion

Here we used a unified framework for building large scale DARs for each facet of mammal diversity. The spatial scaling of each facet revealed that PD/FD reach their maximal diversity faster than SAR suggesting that PD/FD might be less vulnerable than SR to habitat loss. In addition, we extracted the area effect on the diversity of individual terrestrial ecoregions to identify multifaceted hotspots of diversity. We showed that multifaceted hotspots are not necessarily congruent and, thus, that SR, PD and FD are not necessarily good surrogates for each other, especially when considering species relative coverage.

Although global hotspot identification is important as an initial coarse-scale assessment of the conservation value of different regions [START_REF] Lamoreux | Global tests of biodiversity concordance and the importance of endemism[END_REF], several challenges would need to be addressed before our results could be directly transferred into conservation planning actions. Rows correspond to different biomes, while columns represent differences between Diversity Area Relationship (DAR) : PDAR-SAR and FDAR-SAR. For each plot, the differences between PDAR/FDAR and SAR are represented for 3 values of q: 0 (Faith cor index), 1 (Allen cor index) and 2 (Rao cor index). Positive differences mean that PDAR or FDAR are higher than SAR. Area is given in km 2 . Trop moist Forest =Tropical moist forest, Medit. F. = Mediterranean forest. Table 1 The set of diversity indices used in the analysis T represents the height of the phylogenetic tree or the functional dendrogram. PD = Phylogenetic diversity; FD = Functional diversity. * The q parameter refers to the Chao's framework (see Equation 2).

. Author manuscript; available in PMC 2014 November 01.

  The variation in diversity indices explained by area was generally high (median R 2 of the best function in each dataset was 0.5, see Appendix S3) but was quite variable. The R 2 of the best model for each dataset ranged from R 2 =0.0001 (asymptotic model fitting SR in Montane grasslands and shrublands) to R 2 =0.95 (P2 function fitting Rao cor PD in Temperate coniferous forest, see Appendix S3). No single best model outperformed across all data sets, with model selection varying markedly across biomes and diversity indices and reveling substantial levels of uncertainty with different models showing equivalent levels of support (see Appendix S4).

Figure 1 .

 1 Figure 1. Differences between predicted Phylogenetic Diversity Area Relationship (PDAR) / Functional Diversity Area Relationship (FDAR) values and corresponding predicted Species Area Relationship (SAR) values

Figure 2 .

 2 Figure 2. Taxonomic, phylogenetic and functional mammal hotspot selection for tropical moist forests For each biodiversity facet [(1) Species Richness, (2) phylogenetic diversity (Allen cor PD) and (3) functional diversity (Allen cor FD)] a map (A) and a diversity area relationship (B) are presented. Graphs B represent species area relationship (SAR), phylogenetic diversity area relationship (PDAR) and functional diversity area relationship (FDAR). Model fits are shown with a colored curve (see legend) and the averaged fit is presented in black. Red circles indicate hotspots, the larger the diameter, the higher the ranking. Maps A represent the derived ranks from the residuals of the averaged model presented in B.

Figure 3 .

 3 Figure 3. Relationships between the threshold used to define hotspots (expressed in % of ecoregions defined as hotspot) and the similarity between corresponding hotspot listsFrom the left to the right column we compared species richness and phylogenetic diversity hotspots (« SR Vs. PD »), species richness and functional diversity hotspots (« SR Vs. FD »), phylogenetic and functional diversity hotspots (« PD Vs. FD »), species richness and integrative hotspots (agreement between the three facets, « SR. Vs. Int. »). From the top to the bottom row we used Faith cor , Allen cor and Rao cor as PD and FD. The dark continuous line represents mean percentage of congruence of hotspot lists averaged across biomes; the shaded polygon is the associated standard error of the mean. The relative congruence among hotspot lists of two biodiversity facets was determined as the number of ecoregions identified as hotspots by both, divided by the total number of ecoregions in a group.

. Author manuscript; available in PMC 2014 November 01.

Acknowledgments

We thank J. Lawler, G. Cumming and an anonymous reviewer, who provided useful comments on an earlier version of this manuscript. 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). MVC, RDL, and JAFDF received productivity research grants from CNPq, Brazil.

Biosketch

F. Mazel is a PhD student mostly interested in macroecology and macroevolution. Specifically, he aims at describing and understanding the distribution of biodiversity in the light of ecology, evolution, paleontology and paleoclimatology.

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.

Table 2

The five hottest hotspots of tropical and subtropical moist broad leaf forest.


Anne-Christine Monnet

Frédéric Jiguet

Christine N Meynard

David Mouillot

Nicolas Mouquet

Wilfried Thuiller

Vincent Devictor
email: vincent.devictor@univ-montp2.fr.

Asynchrony of taxonomic, functional and phylogenetic diversity in birds

Keywords: Beta diversity, breeding bird survey, functional traits, homogenization, Rao, species turnover, temporal dynamics

Aim-We assessed the temporal trends of taxonomic, functional and phylogenetic diversities in the French avifauna over the last two decades. Additionally, we investigated whether and how this multifaceted approach to biodiversity dynamics can reveal an increasing similarity of local assemblages in terms of species, traits and/or lineages.

Location-France.

Methods-We analysed a large-scale dataset that recorded annual changes in the abundance of 116 breeding birds in France between 1989 and 2012. We decomposed and analysed the spatiotemporal dynamics of taxonomic, phylogenetic and functional diversities and each of their α-, βand γ-components. We also calculated the trend in the mean specialization of bird communities to track the relative success of specialist versus generalist species within communities during the same period.

Results-We found large variation within and among the temporal trends of each biodiversity facet. On average, we found a marked increase in species and phylogenetic diversity over the period considered, but no particular trend was found for functional diversity. Conversely, changes

INTRODUCTION

Past as well as recent rapid global change has already triggered community reshuffling and important species range shifts across the globe [START_REF] Parmesan | Ecological and evolutionary responses to recent climate change[END_REF]). An increasing number of studies have also quantified the potential future impacts of global change on the composition of local communities over large areas by sorting particular species, functional traits or lineages [START_REF] Thuiller | Consequences of climate change on the tree of life in Europe[END_REF][START_REF] Viol | More and more generalists: two decades of changes in the European avifauna[END_REF][START_REF] Buisson | Toward a loss of functional diversity in stream fish assemblages under climate change[END_REF]. In this context, it is now widely acknowledged that strategies for conserving biodiversity should no longer focus solely on species richness, but also on the ecological functions performed by species and on the evolutionary history supported by lineages [START_REF] Diniz-Filho | Darwinian shortfalls in biodiversity conservation[END_REF]. Functional diversity allows us to account for the differences among species in their ecological traits, and is a key factor in understanding ecosystem functioning [START_REF] Naeem | The functions of biological diversity in an age of extinction[END_REF]. For instance, the disappearance of species that represent unique functional traits may have drastically different consequences from the disappearance of the same number of species having common functional traits [START_REF] Mouillot | A functional approach reveals community responses to disturbances[END_REF]. Phylogenetic diversity, on the other hand, can represent a proxy for functional diversity (i.e. more lineages represent more functions; [START_REF] Cadotte | Using phylogenetic, functional and trait diversity to understand patterns of plant community productivity[END_REF] or be considered as an important aspect of evolutionary history of conservation interest [START_REF] Winter | Phylogenetic diversity and nature conservation: where are we?[END_REF].

Taken together, the use of a multifaceted approach to describe the spatial distribution of biodiversity has recently proved to be useful for describing how biodiversity covaries in space, unravelling assembly mechanisms [START_REF] Cadotte | The ecology of differences: assessing community assembly with trait and evolutionary distances[END_REF] and understanding the drivers explaining community composition (e.g. [START_REF] Pavoine | Hierarchical partitioning of evolutionary and ecological patterns in the organization of phylogenetically-structured species assemblages: application to rockfish (genus: Sebastes) in the southern California Bight[END_REF] for rockfish; Graham et al., 2012 for hummingbirds; Huang et al., 2012 for mammals;or Bernard-Verdier et al., 2013 for plants).

Future scenarios of change in biodiversity have been recently extended to include phylogenetic diversity [START_REF] Thuiller | Consequences of climate change on the tree of life in Europe[END_REF]; changes in the functional composition of communities have also been documented [START_REF] Villéger | Contrasting changes in taxonomic vs. functional diversity of tropical fish communities after habitat degradation[END_REF]. However, whether these different facets of biodiversity change show similar temporal trends remains largely unexplored [START_REF] Magurran | Long-term datasets in biodiversity research and monitoring: assessing change in ecological communities through time[END_REF]. For example, a net increase in regional species richness may be accompanied by a decrease in functional diversity [START_REF] Villéger | Contrasting changes in taxonomic vs. functional diversity of tropical fish communities after habitat degradation[END_REF][START_REF] Baiser | The relationship between functional and taxonomic homogenization[END_REF]. Metrics other than traditional diversity indices may thus be necessary to picture the actual trends.

The interest in the study of these facets of biodiversity can be fostered by spatially partitioning them into α-, βand γ-diversity (Ricotta, 2005a;[START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF]. The trend in γdiversity provides information about the dynamics of the diversity in a specific area, while the trend in α-diversity reflects the trends of each local community within this area. The complementary use of β-diversity reflects the turnover (taxonomic, functional or phylogenetic) among communities [START_REF] Bernard-Verdier | Partitioning phylogenetic and functional diversity into alpha and beta components along an environmental gradient in a Mediterranean rangeland[END_REF]. This decomposition of biodiversity facets into α-, βand γ-components was shown to be valuable in conservation or macroecological studies [START_REF] Melo | Environmental drivers of beta-diversity patterns in New-World birds and mammals[END_REF][START_REF] Meynard | Beyond taxonomic diversity patterns: how do α, β and γ components of bird functional and phylogenetic diversity respond to environmental gradients across France?[END_REF]. Nevertheless, it remains unclear whether these complementary components have similar temporal dynamics. For instance, after a decrease in taxonomic β-diversity, functional β-diversity may or may not decrease, depending on how functional traits are distributed among communities.

Apart from multifaceted approaches, community-weighted means (CWM) were specifically developed to measure the relative change in abundances (or occurrences) of species with specific traits (or any characteristics) within communities [START_REF] Ricotta | CWM and Rao's quadratic diversity: a unified framework for functional ecology[END_REF]. These traits can be a priori selected to reflect species-specific sensitivities to a given pressure. For instance, following habitat or climate changes, individuals (or species) sensitive to those changes should be replaced locally by other individuals (or species) that benefit from the same changes. This approach has been successfully used to show that, following habitat disturbance, species dependent on a few habitat types (specialist species) tend to be replaced by generalist species [START_REF] Kampichler | Large-scale changes in community composition: determining land use and climate change signals[END_REF][START_REF] Viol | More and more generalists: two decades of changes in the European avifauna[END_REF], a process interpreted as a specific form of biotic homogenization [START_REF] Clavel | Worldwide decline of specialist species: toward a global functional homogenization?[END_REF]. While the complementarity between CWM and other diversity metrics has been proposed to better describe community reshuffling [START_REF] Ricotta | CWM and Rao's quadratic diversity: a unified framework for functional ecology[END_REF], the relevance of these two alternative approaches has never been supported by empirical data at large temporal and spatial scales.

Here, we conducted a temporal analysis of biodiversity dynamics (i.e. taxonomic, functional and phylogenetic diversities and their respective spatial turnover) and of a CWM measuring the average community specialization of bird communities. We analysed data from a highresolution survey of birds conducted over two decades in France (1989-2012), a time period during which large climate [START_REF] Jiguet | Population trends of European common birds are predicted by characteristics of their climatic niche[END_REF] and land-use changes [START_REF] Antrop | Landscape change and the urbanization process in Europe[END_REF] have taken place in Europe. More specifically, we: (1) assessed whether taxonomic, functional and phylogenetic diversities (and their respective spatial turnover) had different temporal trends over the same period; and (2) investigated the complementarity of multifaceted approaches and CWM to reveal potential change in the similarity of local assemblages in terms of species, traits and/or lineages.

METHODS

Bird data

We analysed data from the French breeding bird survey (FBBS), which provided the abundances of French avifauna over a 24-year period on an annual basis (from 1989 to 2012) [START_REF] Jiguet | French citizens monitoring ordinary birds provide tools for conservation and ecological sciences[END_REF]. The monitoring programme used two schemes. These two schemes both used standardized protocols and followed the same basic principles ensuring the possibility of comparing the yearly changes in the relative abundances of species from year to year. During the first scheme (corresponding to the period 1989-2001), survey routes were freely chosen by observers (note, however, that these routes were located in various habitats including farmlands, forests and urban areas) and the same observer monitored the same route from year to year. Along each route, the observer had to monitor from 10 to 15 point counts evenly distributed within the habitats along the route. The same observer monitored exactly the same point counts in the same order over 5 min. Population trends obtained from this first scheme were comparable to those obtained in other European countries with a random sampling of surveyed sites [START_REF] Julliard | Common birds facing global changes: what makes a species at risk?[END_REF]. From 2001, a new scheme was launched, in which the same protocol was used, but instead of being freely chosen by the observers, 2 km × 2 km plots were randomly selected among 80 possible plots around the locality proposed by the observer [START_REF] Jiguet | French citizens monitoring ordinary birds provide tools for conservation and ecological sciences[END_REF]. As in the former protocol, the same observer monitored the same 10 point counts evenly distributed within the plot from year to year.

In both schemes, point counts were monitored on approximately the same day of the year (± 7 days), at the same time of the day (± 15 min, within 1-4 h after sunrise), for 10 min. Sampling sessions were repeated twice a year during the breeding period, once before and once after the pivotal date of 8 May (4-6 weeks apart). This allowed a sampling of early singing species and late migrants. For each point count, in each monitored site (routes or plots), each year, we took the maximum abundance of each species recorded during the two sessions as a measure of the relative abundance of the species. To ensure that each site included the same number of point counts over the 24-year period, we selected 10 point counts per site. Among the 1818 sites retained, each site was monitored at least twice and on average 5.7 ± 2.8 years (mean ± SD; Fig. 1).

The FBBS focused on common birds that regularly breed in France (about 300). Among those, we selected species that have been monitored at least twice at two different sites. To avoid the influence of rare species not correctly monitored by the protocol (e.g. wetland species) we only retained the most common species, i.e. we retained species with individuals representing 99% of all individuals monitored in the database; which correspond to 116 breeding species.

Measuring taxonomic, functional and phylogenetic diversity

We adapted the method developed by [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] to describe the taxonomic, functional and phylogenetic diversity of bird communities. Here, a community was considered as a pool of species that co-occur in a given point count in a given year. In brief, we used the [START_REF] Rao | Diversity and dissimilarity coefficients: a unified approach[END_REF] quadratic entropy index, which offers the same mathematical framework for all three facets of diversity, allowing straightforward comparisons between their trends [START_REF] De Bello | The partitioning of diversity: showing Theseus a way out of the labyrinth[END_REF]. The Rao index is given by the formula , where d ij is the distance between species i and j, p i and p j are their relative abundances and S is the total number of species in the community. The Rao index integrates the abundance of species and a measure of distance between species (d ij ) which can be either functional or phylogenetic. Taxonomic diversity was estimated with the same index, except that d ij = 1 was used between each pair of species (in this case, the Rao index was thus reduced to the Gini-Simpson diversity index; Ricotta, 2005a).

Matrices of pairwise functional and phylogenetic distances were both required to estimate functional and phylogenetic diversity. We obtained pairwise functional distances among the 116 species from a set of 22 functional traits. These traits embraced major life-history traits, feeding habits and morphological characteristics of birds [START_REF] Petchey | Low functional diversity and no redundancy in British avian assemblages[END_REF]Devictor et al., 2010). From these traits (see Appendix S1 in Supporting Information for a detailed description of each trait), we calculated pairwise distances between species using the Gower distance, and then we produced a dendrogram (UPGMA clustering, selected from the consensus method of [START_REF] Mouchet | Towards a consensus for calculating dendrogram-based functional diversity indices[END_REF] to ensure that distances between species pairs were ultrametric. Note that in measuring functional diversity (whatever the metric considered), the included traits and their specific categorization may rely on arbitrary decisions. In our case, we used major ecological traits described in birds as being sensitive to large-scale environmental filters [START_REF] Petchey | Low functional diversity and no redundancy in British avian assemblages[END_REF][START_REF] Meynard | Beyond taxonomic diversity patterns: how do α, β and γ components of bird functional and phylogenetic diversity respond to environmental gradients across France?[END_REF]. Most of these traits are also directly or indirectly involved in ecosystem functioning (Şekercioğlu, 2006). However, other trait selections would also have been possible to test more specific variation in community composition.

Ultrametric phylogenetic distances between the 116 species were directly extracted from a dated, calibrated molecular phylogenetic tree assembled by [START_REF] Thuiller | Consequences of climate change on the tree of life in Europe[END_REF] (see Appendix S2). For this phylogeny, consistent estimates of branch lengths were available. Maximum distances in both the phylogenetic and functional matrices were standardized to a maximum value of 1 before further analysis.

We calculated the three components for each biodiversity facet: γ-diversity (total diversity of a site), β-diversity (turnover between point counts within a site) and α-diversity (diversity in a point count). Because α-diversity can be fully determined by γand β-diversities, we focused only on γand β-diversities which were obtained from the standard decomposition of the Rao quadratic entropy: γ = βα (Ricotta, 2005a). Although the division of biodiversity into α-, βand γ-components has often been applied to a larger spatial scale/ region, this landscape decomposition is meaningful and relevant for bird communities. Indeed, γ-diversity (here considered as representing the bird diversity of a landscape) is composed of heterogeneous habitats, and thus β-diversity (turnover between point counts within a site) represents intra-landscape variability for which birds were shown to be sensitive [START_REF] Devictor | Functional biotic homogenization of bird communities in disturbed landscapes[END_REF]. All indices were then transformed according to [START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF] to normalize the properties of all diversity metrics.

Measuring specialization of communities

To describe the dynamics of the relative proportion of specialized species in communities, we calculated the community specialization index (CSI) for each site in each year. CSI is a community-weighted mean representing the average of the species-specific level of habitat specialization of co-occurring species, weighted by their abundances [START_REF] Julliard | Spatial segregation of specialists and generalists in bird communities[END_REF]. We quantified the degree of habitat specialization for a species (the species specialization index, SSI), following the approach of [START_REF] Julliard | Spatial segregation of specialists and generalists in bird communities[END_REF]. SSI is calculated as the coefficient of variation (SD/mean) of species densities across habitat classes (we used the SSI values available in [START_REF] Devictor | Distribution of specialist and generalist species along spatial gradients of habitat disturbance and fragmentation[END_REF] based on species densities recorded in 2005 by the FBBS). This approach assumes that a given species is more specialized to certain habitat classes if its density is higher there than elsewhere. Note that SSI values were shown to be robust to the change in the habitat classification considered (Le [START_REF] Viol | More and more generalists: two decades of changes in the European avifauna[END_REF]. CSI is then given by the average of each species' SSI, weighted by the species' abundance within each site:

where p ij is the relative abundance of species i in site j and S is the total number of identified species in the site. One expects that the CSI calculated for a given site will decrease following a relative increase in generalists in that site [START_REF] Julliard | Spatial segregation of specialists and generalists in bird communities[END_REF].

To assess the potential contribution of generalist species to the observed changes in functional and phylogenetic diversity, we also tested the linear relationships between species specialization and their functional and phylogenetic originality, using linear regressions. Specialization (SSI) was considered as the response variable and functional (or phylogenetic) originality, measured as the average pairwise functional (or phylogenetic) distance of that species with all others, as the independent variable.

Statistical analyses

To estimate the temporal dynamics of the different community indices, we used a two-step approach. We first modelled the change in the yearly average of each diversity index. To do this we used a generalized additive mixed model (GAMM) designed as follows to meet the general structure of the framework proposed by [START_REF] Dornelas | Quantifying temporal change in biodiversity: challenges and opportunities[END_REF]. The community measures (taxonomic, functional and phylogenetic diversity, or CSI) were considered as the response variables. Sites were considered as a random effect to account for variability in the level of indices among sites. We integrated geographical coordinates in isotropic smooth terms, according to the methods of [START_REF] Wood | Generalized additive models: an introduction with R. Chapman and Hall/CRC[END_REF], to account for structural spatial gradients. In this model, we also accounted for temporal autocorrelation structure by fitting an autoregressive model of order 1. Note that a more formal integration of spatial autocorrelation was not possible in these models that already accounted for temporal autocorrelation. However, the study of the semi-variograms of the residuals of that model confirmed that no spatial structure remained in the residuals after the smoothing of coordinates. In this model, time (in years) was considered as a discrete variable. As we expected that the size of the community could affect diversity indices, taxonomic diversity (TD) was also systematically included as a covariate when analysing functional diversity (FD) and phylogenetic diversity (PD) and CSI. β-TD was also systematically added for the models involving β-FD or β-PD. This first model provided us with the effect of each transition for years t to t + 1 on each index and its corresponding standard error. We arbitrarily fixed all community measures to 100 for the year 2001, considered as the reference year.

Then, a second GAM model was used to summarize the overall changes (and their nonlinear trends, if any) in the yearly estimates of each index over the 24-year period. To do this, we used, for each community measure, its yearly estimates (provided by the first model) as a response variable and we used a smoothing structure with a fixed degree of freedom of 4 for the smoothing function of year, considered as a continuous variable. In this second model, the sampling effort (i.e. the number of sites monitored per year) was also used as a covariate. The variance of the yearly estimates (given by the first GAMM) was also used as a weight in this second model. We tested whether this model could effectively describe a nonlinear trend in each index.

In parallel, we repeated the whole analysis presented above, but with presence/absence data rather than abundances. We conducted this analysis to assess the importance of including abundance in estimating the facets of biodiversity [START_REF] Newbold | Mapping functional traits: comparing abundance and presence-absence estimates at large spatial scales[END_REF].

All calculations of indices, randomizations for the null model and statistical analyses were performed with R, version 3.0.1 (R Core Team, 2013).

RESULTS

Although there is no accepted way to formally estimate model fit for the first GAMMs [START_REF] Wood | Generalized additive models: an introduction with R. Chapman and Hall/CRC[END_REF][START_REF] Zuur | Mixed effects models and extensions in ecology with R[END_REF], each index shows very important variation for a given year and among years (Fig. 2). Beyond these variations, the yearly changes in all diversity metrics could also be described by nonlinear trends from 1989 to 2012 with one or two specific inversions in their direction (Fig. 2). For instance, the general trend in species diversity corresponded to a nonlinear curve, with three successive periods (F 3,24 = 4.8, P = 0.01, R 2 = 0.87): a relatively stable period of 7 years (1989-95), followed by a consistent increase over 13 years (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) and a new stable (or even slightly decreasing) period of 5 years . The trends of local species richness and local abundance were qualitatively similar, showing that the local increase in taxonomic diversity was accompanied by the local increase in both species richness and abundances from approximately 1996-2008 (Appendix S3). Trends in γ-FD and γ-PD (adjusted to variation in taxonomic diversity) had different shapes. While the general trend in γ-FD was nonsignificant (F 3,24 = 2.18, P = 0.12; Fig. 2c), the trend in γ-PD followed a two-step period with a marked and consistent increase from 1995 to 2012 (F 3,24 = 24.9, P < 0.001, R 2 = 0.82; Fig. 2e). Overall, each of these facets (γ-TD, γ-FD and γ-PD) had their own general trajectories.

The corresponding trends in β-diversities (β-TD, β-FD and β-PD) were also different. In particular, β-TD followed a curvilinear trend with only one turning point around 2002 (F 3,24 = 5.51, P = 0.006, R 2 = 0.54; Fig. 2b) and a marked decrease since 1997, while β-FD followed a general trend with three different periods (F 3,24 = 4.9, P = 0.010, R 2 = 0.78; Fig. 2d) and β-PD did not show particular long-term trends (F 3,24 = 2.12, P = 0.12; Fig. 2f). Note that these trends in β-diversities were different from each other and different from the trends in γ-diversities.

Finally, during the same period, CSI, reflecting the relative abundance of specialist species in local sites, followed a two-step period with a marked and consistent decrease from 1999 to 2012 (F 3,24 = 8.14, P < 0.001, R 2 = 0.61; Fig. 2g). In other words, since 1999, local assemblages have tended to include more individuals of generalist species. We further found that, at the species level, the more specialist species are, the higher their values of functional (Student's t-test: t = 3.98, P < 0.001, R 2 = 0.12) and phylogenetic (Student's t-test: t = 2.98, P = 0.003, R 2 = 0.07) originality.

Beyond the change of the average trends over the 24-year period, this temporal analysis showed high yearly variation within and among indices. All trends described above thus mask several short-term increases or decreases of each facet. Noticeably, while some indices increased during a given time window (e.g. β-FD during 2003-12) others decreased during the same period (γ-TD). The γ-diversity of a given facet could also increase for a given time period while the corresponding β-diversity of this same facet was decreasing during the same period (e.g. γ-TD and β-TD, from 1997-2005). Note that the deviance explained by sampling effort ranged from 0.04 to 10.11% for all community measures and that the effect of this variable was never significant (P < 0.05 for all community measures).

We found similar general trends for every diversity index when calculated from presence/ absence data (Appendix S4) showing that beyond the change in the local relative abundances, changes in the identity of species recorded in local sites were also responsible for the observed changes in diversity facets.

DISCUSSION

In this paper we first assessed whether taxonomic, functional and phylogenetic diversities (and their respective spatial turnover) can have different temporal trends over the same period. Among major changes experienced by bird assemblages, a general increase in taxonomic diversity and phylogenetic diversity (since approximately 1995) as well as the overall decrease in CSI, β-TD and β-FD, can be emphasized. However, during this period, the temporal trends of each facet are best described by the succession of mid-term fluctuations than by linear trends.

These fluctuations suggest that each of these facets mirrors different and complementary aspects of community change. For instance, an increase in community taxonomic diversity and the decrease in community specialization occurred simultaneously (Fig. 2a,g). While the populations of some species may have increased recently due to favourable conservation policies [START_REF] Donald | International conservation policy delivers benefits for birds in Europe[END_REF], the general increase in local species diversity coupled with the decrease in community specialization is likely to reflect the biotic homogenization of the French avifauna also described at the European level (Le [START_REF] Viol | More and more generalists: two decades of changes in the European avifauna[END_REF].

We also found that local changes in diversity are scale dependent. In addition, the increase in species richness observed at a given spatial scale is often coupled with a decrease in diversity at a different scale [START_REF] Van Turnhout | Scale-dependent homogenization: change in breeding bird diversity in the Netherlands over a 25-year period[END_REF]. Our results allow us to extend this conclusion to other facets of diversity. For instance, a directional change in β-TD was not necessarily coupled with a similar change in β-FD or β-PD. As also suggested by [START_REF] Baiser | The relationship between functional and taxonomic homogenization[END_REF], these results show that the biotic homogenization of a given component of biodiversity (e.g. taxonomic diversity) is thus not necessarily coupled with the homogenization of other components. Thus, the local increase in traditional diversity metrics (species richness, species diversity measure, total abundance) can mask different trends in other facets [START_REF] Magurran | Long-term datasets in biodiversity research and monitoring: assessing change in ecological communities through time[END_REF]. But our results additionally suggest that this scale dependency is also valid over time. In other words, different diversity indices may have different temporal inertias. To our knowledge, such a potential decoupling of diversity indices at large spatial and temporal scales has hardly been investigated.

Obviously, not a single metric of biological diversity can be considered as a silver bullet. In particular, all metrics of species diversity (e.g. the Gini-Simpson index) derived from information theory have their intrinsic drawbacks (e.g. being bounded within a given range of values, having non-uniform properties across this range, or representing rare or common species differently). The Rao quadratic entropy used in this study also has its strengths and weaknesses (Ricotta, 2005b) but has the great advantage of offering a coherent and consistent framework for measuring several facets of diversity. In particular, we can be confident that the differences in the temporal trends observed among indices (taxonomic, functional and phylogenetic diversity) were not produced by changes in the indices considered (only the pairwise distance between species differed among these indices). However, a similar analysis could benefit from calculating similar temporal trends with complementary approaches (e.g. using the functional diversity indices developed by [START_REF] Villéger | New multidimensional functional diversity indices for a multifaceted framework in functional ecology[END_REF]. Some of the observed variation in yearly indices can be explained by the sampling design itself [START_REF] Magurran | Long-term datasets in biodiversity research and monitoring: assessing change in ecological communities through time[END_REF]. In particular, variation in the sampling effort or in the protocol could induce uncontrolled errors in the fluctuations of species abundances. To draw general trends, we thus systematically included the number of monitored sites as a covariate to account for yearly variation in sampling effort. Note also that the same observer monitored a given site with the same protocol (either from 1989 to 2001 or from 2001 to 2012) ensuring among-year standardization of diversity measures. Moreover, while many unexplained sources of variation of the diversity indices could not be described, this variation does not necessarily represent sources of bias for our general conclusions. Indeed, we were interested in the relative change in several indices measured with the same data. Unless species representing specific lineages, traits or habitat specialization were monitored differently during the period considered, one does not expect to find the asynchrony of indices observed during the same period. In particular, important short-term changes in diversity indices could be observed during periods within which the protocol and the sampling effort were similar (e.g. from 1989 to 2000 and from 2001 to 2012 steep and idiosyncratic changes in some indices can still be observed; Fig. 2).

Even with standardized protocols and constant sampling efforts, abundances monitored from large-scale surveys should still be handled with caution. Many sources of error have been shown to influence the local abundances recorded by a given observer in a given site (e.g. variability in detectability among species, in observers' ability to detect species or even in meteorological conditions, [START_REF] Sauer | Observer differences in the North American Breeding Bird Survey[END_REF]. However, while these sources of error can induce strong bias when true abundances are studied, they should not affect our conclusions derived from the relative changes in abundances from year to year in each local point count [START_REF] Bas | Accounting for weather and time-of-day parameters when analysing count data from monitoring programs[END_REF]. Moreover, the same analysis conducted on presence/absence data (i.e. removing variation in abundances between species) provided similar quantitative results (see Appendix S4). [START_REF] Newbold | Mapping functional traits: comparing abundance and presence-absence estimates at large spatial scales[END_REF] have recently assessed the importance of including abundance information when mapping functional diversity and concluded that presenceabsence data often already yield sensitive indices.

The causes of the observed changes in diversity facets remain to be confirmed. Ideally, to disentangle the deterministic and stochastic processes that are potentially responsible for the local changes in biodiversity, the observed trends (and their variances) should be explicitly coupled with local changes in environmental variables. In practice, however, these analyses are difficult to achieve at a large spatial scale due to the absence of land-cover surveys with sufficiently fine resolutions (e.g. the CORINE land-cover data available for France do not record landscape elements lower than 25 ha and are therefore only able to detect major changes). Using a spatial analysis of how diversity indices were related to habitat structure and compositions, [START_REF] Meynard | Beyond taxonomic diversity patterns: how do α, β and γ components of bird functional and phylogenetic diversity respond to environmental gradients across France?[END_REF] suggested that low levels of human impacts should generally favour all three facets of diversity (taxonomic, functional and phylogenetic) while higher β-diversities should be observed in heterogeneous landscapes. But whether similar relations are observed with temporal dynamics of diversity indices and temporal changes in landscape characteristics has yet to be confirmed. Our results, along with recent theoretical simulations [START_REF] Münkemüller | From diversity indices to community assembly processes: a test with simulated data[END_REF], suggest that, following habitat or climate modifications, some indices could respond faster than others due to the specific delay with which each species can react and the inherent velocity of the dominant assembly rules.

Measuring the consequences of the observed changes is also a major challenge for macroecology and conservation biogeography. Our results reveal rather small relative changes for most facets. However, these changes result from the fluctuation of common species monitored in hundreds of sites distributed all over the country, which represent large variations in the number of species and individuals involved (on average the fluctuations of approximately two species and 40 individuals at local sites was observed during the period considered). Some of the most common species have indeed declined or increased drastically during this period [START_REF] Jiguet | Population trends of European common birds are predicted by characteristics of their climatic niche[END_REF], potentially affecting specific functions and interactions after population depletions [START_REF] Gaston | Commonness, population depletion and conservation biology[END_REF]. We further showed that specialist species tend to have more original functional and phylogenetic characteristics than generalists. A more detailed species-by-species approach targeting the specific contributions of particular species (e.g. specialist versus generalist, protected versus non-protected) should help refine these findings.

Overall, our study presents the first temporal, multifaceted, multiscale analysis of change in community composition over large spatial and temporal trends. Our results emphasize that asynchrony can occur between different diversity facets so that, more than a given trend per se, it is the jointly studied trends of complementary biodiversity facets (including CWM) that provide enlightening information on major changes in the structure and composition of local communities. Beyond its ecological relevance for bird communities, our approach should provide a practical framework for evaluating temporal trends of different biodiversity facets, which can be applied to other taxonomic groups and is a step forward towards adopting more integrative approaches to the study of biodiversity on large spatial and temporal scales.

Figure 1 .

 1 Figure 1. Spatial distribution of sites from the French breeding bird survey. Example of the distribution of 10 point counts (black filled circles) within a plot (a 2 km × 2 km square). Each black circle within the plot corresponds to a local community.

Figure 2 .

 2 Figure 2. Temporal trends of diversity facets over 1989-2012. Taxonomic (a) γ-diversity and (b) βdiversity. Functional (c) γ-diversity and (d) β-diversity. Phylogenetic (e) γ-diversity and (f) β-diversity. (g) Community specialization index (CSI). Yearly changes in diversity indices (and their standard error, in the grey band) were obtained from a model accounting for spatial gradients and temporal autocorrelation. Variations in (γ or β) functional and phylogenetic diversities and in CSI were adjusted to variations in (γ or β) taxonomic diversity. The y-axis represents the relative variation of the facet considered compared with

. Author manuscript; available in PMC 2014 July 25.

Acknowledgments

We sincerely thank the hundreds of volunteers that took part in the French breeding bird survey (STOC EPS programme). This project was supported by the 'Fondation pour la Recherche sur la Biodiversité' (FRB). W.T. received support from the European Research Council under the European Community's Seven Framework Programme FP7/2007-2013 Grant Agreement no. 281422 (TEEMBIO).

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.

BIOSKETCH

The research group combines several skills including theoretical research on community ecology, biogeography and conservation. Anne-Christine Monnet is doing a PhD in the impacts of global change on biodiversity using different approaches from species to communities.