The Evolutionary Legacy of Diversification Predicts Ecosystem Function

Theory suggests that the structure of evolutionary history represented in a species community may affect its functioning, but phylogenetic diversity metrics do not allow for the identification of major differences in this structure. Here we propose a new metric, ELDERness (for Evolutionary Legacy of DivERsity) to estimate evolutionary branching patterns within communities by fitting a polynomial function to lineage-through-time (LTT) plots. We illustrate how real and simulated community branching patterns can be more correctly described by ELDERness and can successfully predict ecosystem functioning. In particular, the evolutionary history of branching patterns can be encapsulated by the parameters of third-order polynomial functions and further measured through only two parameters, the “ELDERness surfaces.” These parameters captured variation in productivity of a grassland community better than existing phylogenetic diversity or diversification metrics and independent of species richness or presence of nitrogen fixers. Specifically, communities with small ELDERness surfaces (constant accumulation of lineages through time in LTT plots) were more productive, consistent with increased productivity resulting from complementary lineages combined with niche filling within lineages. Overall, while existing phylogenetic diversity metrics remain useful in many contexts, we suggest that our ELDERness approach better enables testing hypotheses that relate complex patterns of macroevolutionary history represented in local communities to ecosystem functioning.


Introduction
The phylogenetic structure of ecological communities conflates diversification of lineages, geographic movement of lineages, and their local assembly. For instance, the presence of many ancient lineages in a community may result from particularly strong diversification or geographic expansion of these ancient lineages, which are hence present in the species pool from which the community is sampled. Moreover, the presence of such ancient lineages in a community may result from a high capacity to locally coexist. The composition of ecological communities in turn drives ecosystem functioning (Tilmanet al. 2001). Phylogenetic structure of a community hence integrates the diversification, expansion, and local coexistence of lineages with the functioning of ecosystems. While such a relationship does not provide mechanistic explanation of ecosystem functioning-functional traits remain indispensable for this purpose-it demonstrates that past diversification can potentially drive present ecosystem functioning. The most basic example of such a diversification/functioning relationship is probably the relationship between high species richness and primary productivity. Evolution of numerous species and their local coexistence increase species richness of local communities, and species richness often increases productivity (e.g., Venail et al. 2015). However, increased species richness may represent species stemming from ancient, recent, or continuous diversifications. To date, we do not know how productivity or other ecosystem functions depend on the presence of species representing particularly ancient, recent, or continuous diversifications.
The coexistence of species that diversified anciently, recently, or continuously may differently affect ecosystem productivity. First, many functional traits and biotic interactions show phylogenetic signal: recently diversified species share more functional traits and biotic interactions than do anciently diversified species (Losos 2008, phylogenetic conservatism sensu Wiens et al. 2010;Gomez et al. 2010). A community in which species stemming from ancient diversifications coexist may hence be functionally highly diverse (Cadotte et al. 2008, but see Prinzing et al. 2008), which in turn has been shown to increase ecosystem productivity (Cadotte et al. 2008). Also, phylogenetic signal of traits within plant and herbivore lineages impedes herbivores from switching among major ancient plant lineages (Wiens et al. 2010). Coexistence among such anciently diversified plant species may hence impede exchange of herbivores among plants and reduce herbivore impact on plant productivity. Second, recent diversification within lineages may lead to niche filling-that is, to an increasingly efficient partitioning and exploitation of available resources among species (Ricklefs 2007;McPeek 2008). Niche filling might also decrease competition among species, which may respond by increased growth (Venail et al. 2013). A community in which recently diversified species coexist might hence show increased productivity. Third, trait conservatism and niche filling might act together, so that com-munities composed of species stemming from a constant continuous diversification might be most productive: the species representing the earlier diversifications might occupy distinct major niches, and species representing recent diversification within lineages might efficiently partition these niches. Finally, environments may filter particular traits represented by particular major lineages-that is, lineages originating from an intermediate age of diversification (Grime 2006). Coexisting species resulting from such intermediate diversifications might hence be sufficiently closely related to share traits that permit passage of a local environmental filter and sufficiently distantly related to profit from complementarity effects facilitating coexistence, such as limiting similarity (Silvertown et al. 2006). Both effects would increase productivity. Overall, the above scenarios suggest that primary productivity might increase in the presence of species from most ancient, most recent, continuous, or intermediate diversifications.
Until now, phylogenetic diversification, or more generally the evolutionary history represented within a community, has often been inferred from a simple metric of phylogenetic diversity, such as the mean or sum of phylogenetic branch lengths connecting pairs of species (e.g., Proches et al. 2006;Prinzing et al. 2008;Cadotte et al. 2009; and many other studies). However, these simple metrics often fail in differentiating even very different patterns of accumulation of ancient or recent lineages represented by species in a community ( fig. 1). Accumulated lineages come with different functional traits that ultimately lead to ecosystems with specific functioning. For example, imagine that phylogeny B in figure 1 represents a temperate forest community comprising four broadleaf (angiosperm) trees from different families (e.g., oak, birch, maple, and ash). Phylogeny D, which has identical total branch lengths (i.e., Faith's phylogenetic diversity [PD]), represents a forest comprising two closely related angiosperms (e.g., oaks) and two closely related gymnosperms (e.g., pines). As the host breadth of many herbivorous insect species is confined to a single plant family (e.g., Ødegaard et al. 2005;Weiblen et al. 2006;Pearse and Altermatt 2013), trees in B may not share many herbivores, whereas some tree species in phylogeny D certainly would, potentially increasing phytophagy and reducing productivity (e.g., Yguel et al. 2011). Thus, evolutionary processes (trait evolution of the trees and herbivores and resulting community assembly of trees and herbivores) may influence this specific interactive community and its functioning. Overall, phylogenetic community structure can link these processes of diversification, community assembly, and ecosystem functioning, yet standard descriptors of phylogenetic community structure miss potentially important information, notably on the pattern of lineage accumulation through time.
Evolutionary biologists typically describe diversification patterns of a single lineage by tree imbalance (e.g., Heard 1992) or by lineage-through-time (LTT) plots. However, so far, these plots are mathematically described only as continuous shifts in diversification (Pybus and Harvey 2000;Machac et al. 2013). Moreover, diversification patterns are used to infer continuous shifts in two underlying processes-speciation and extinction-and to reconstruct past diversities (Stadler 2013 and references therein). Community ecologists, however, trim phylogenies to all lineages within a local community, and such community phylogenies can display more complicated relatedness patterns with noncontinuous shifts between recent, ancient, or constant branching. These patterns reflect, besides speciation and extinction, a multitude of processes determining which particular species and lineages actually colonize and survive in local communities (e.g., Mayfield and Levine 2010;Valiente-Banuet and Verdú 2013). An existing application of LTT plots to communities uses only a single axis along which LTT plots are differentiated continuously from ancient to recent diversifying (Martin 2002). More complex LTT patterns have been analyzed only by disintegrating the overall LTT pattern into numerous discrete time steps and analyzing each separately (Stadler 2011). Yet, we are not aware of tools for quantifying integratively the overall shape of a complex LTT pattern with a manageable number of parameters (or without making mechanistic assumptions on speciation and extinction that are unjustified in the case of species communities). Accordingly, it was not possible until now to test how these complex phylogenetic patterns may be related to ecosystem processes.
Here, we report a new approach describing diversification pattern, namely ELDERness (for Evolutionary Legacy of DivERsity), and evaluate its ability to capture the temporal dimension of relatedness patterns within a community on the basis of a comprehensive description of lineage accumulation through time. In brief, ELDERness describes the evolutionary history of lineages within a given commu-nity through the parameters of a polynomial function fitted to the phylogeny of the species in that community ( fig. 2). The polynomial regression curves can then be compared to a perfectly constant diversification represented by a straight line of LTT plots. Specifically, the first and second areas of the curve that are above or below the straight line (hereafter referred to as ELDERness surfaces ES1 and ES2, respectively) can be quantified ( fig. 2). We address two key questions. First, can the ELDERness surfaces differentiate among simulated branching scenarios better than common metrics of phylogenetic diversity, phylogenetic tree imbalance, or diversification metrics? Second, do ELDERness surfaces predict community productivity better than these commonly used metrics? Overall, we link the theoretical scenarios of evolutionary branching of lineages to the ecosystem functioning of local communities ( fig. 2). We compare the observed link to the predictions of the above scenarios.

Simulations and Model Fitting
We simulated four types of phylogenetic trees with differing rates of branch accumulation (constant, recent, ancient, and complex branching) by using a forward-time, purebirth simulation with a length of 30 time units. Four types of phylogenies were simulated using an approach used by Morlon et al. (2011) for their simulation: (1) "constant trees" were simulated with a constant rate of branching (speciation), resulting in a log-linear accumulation of species; (2) "late-branching trees" were simulated using a step function in speciation rates, where the last 25% of the overall time interval of the simulation had a speciation rate that was five times higher than the previous 75%; (3) "ancient-branching trees" were simulated using a step function in speciation rates, where the first 25% of the simulation had a branching rate that was 10 times higher than the subsequent 75%; (4) "complex trees" were simulated with a step function where the first 25% had a branching rate that was 10 times higher than the middle 50% and the last 25% had a branching rate that was five times higher than the middle 50%. Differences in rates of branching were chosen to produce phylogenies while keeping species richness constant. For each type of tree, we simulated 10 trees each of 10, 20, 40, and 80 species to examine the predictive power of ES1 and ES2. Thus, 40 simulated phylogenies of each tree type were created, totaling 160 phylogenies. We explored other types of tree simulations (e.g., more extreme branching rate of different time periods and different mean rates of branching and number of time steps), but all led to the same conclusions as those based on the four tree types presented here. Simulated phylogenies, their LTT plots, and the phylogenetic metrics measured on them are deposited in the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.jd150 (Yguel et al. 2016).
We described the LTT by fitting a curve, because this reduces a complex pattern to a few standardized parameters that can easily be interpreted and compared between branching patterns of communities of different species richness and mean age. Specifically, we fitted to each LTT plot of each simulated phylogeny second-and third-order polynomial regressions that took the form N(t) p 0 1 a 1 t 1 a 2 t 2 1 a 3 t 3 , where N is the number of taxa logarithmically transformed and t is time. Because we always have one lineage at time 0 and always have all species of the community at time t max (i.e., today), we fixed N(0) p ln(1) p 0 and N(t max ) p ln(species richness). Next, we used a change of variables to simplify the equation, where a 1 is fixed at t max and depends on a 2 and a 3 for the thirdorder polynomial: a 1 p ln(species richness)=[t max -a 2 t max -a 3 (t max ) 2 ]. Then, we replaced a 1 by the above expression, and we obtained the following (this could be simplified again but was coded in this form in R for practical reasons): To find a 2 and a 3 , we solved the equation . We used the same method for the second-order polynomial. All these and subsequent analyses were performed in R (R Core Team 2014; see supporting file S1 [supporting files available online] for code to create LTT plots and extract the third-order polynomial parameters). The third-order polynomials generally fit better than the second-order polynomials across all simulations, even penalizing for the number of predictor variables (R 2 and adjusted R 2 between 94.9% and 99.2%; supporting file S2). We note that fits of greater than thirdorder polynomials impeded simple descriptions of curve shape (see below) and thereby straightforward interpretation in terms of ancient versus recent diversification. We also note that second-order polynomials would not produce sufficiently complex curves. Note finally that sigmoidal function might also be adequate but would render the steps of integration below comparatively more difficult.
A given value of a given parameter of a third-order polynomial regression is difficult to interpret. We hence characterized how the polynomial curve deviated from a constant rate of branching (i.e., the ELDERness surfaces; fig. 2 and supporting file S3 for simulated phylogenies). To do so, we measured the surface of the polynomial curve above or below a constant rate of branching using the integral of the polynomial model minus the integral of the constant rate of branching (in the interval [0, t max ]). But in complex LTT patterns, part of the polynomial curve can be above the constant rate of branching (i.e., part of the surface difference is positive) and the other part below the constant rate of branching (i.e., part of the surface difference is negative). To better describe the shape of the polynomial curve, we have thus measured these two surfaces separately. For that, we first identified the intersection between the third-order polynomial curve and the straight line of constant rate of branching. Because our third-order polynomials were fixed at the intercept and at t max , the intersections were found by solving the equation [ln(species richness)=t max ]t p a 1 t 1 a 2 t 2 1 a 3 t 3 , with a 1 expressed as a function of a 2 and a 3 as specified above. This equation can be simplified as t(t max 2 t)[a2 1 a3(t max 1 t)] p 0. It has three solutions: 0, t max , and t int p 2(a 2 =a 3 ) 2 t max . There can thus be an intersection point at t int only if t int lies between t p 0 and t p t max (which are the two other obligatory intersections). If t int is not in the domain [0; t max ], there is only one surface ES1 equal to the integral of the polynomial model minus the integral of the constant rate of branching in the domain [0; t max ]. ES2, the second surface, is equal to 0. If t int is in the domain [0; t max ], the first surface ES1 is measured with the integral of the polynomial curve minus the integral of the constant rate of branching in the domain [0; t int ]. ES2 is measured using the integral of the polynomial curve minus the integral of the constant rate of branching in the domain [t int ; t max ]. All these measures of determinants and surfaces were performed in Excel using the formula indicated.
The ES1 or ES2 can be negative if the surface between the third-order polynomial and the constant rate of branching is below the constant rate of branching, and inversely (i.e., ES1 or ES2 are positive if the surface described previously is above the constant rate of branching; see fig. 2B). A positive ES1 and negative ES2 correspond to a complex pattern with branching rate increasing at the beginning and at the end (see phylogeny A in fig. 2B). A negative ES1 and positive ES2 correspond to a pattern with highest branching rate in the middle (see phylogeny B in fig. 2B). If ES1 and ES2 are close to 0, branching rate will be constant (see phylogeny C in fig. 2B). Our approach also permits identification A B ×40 ( of a range of ancient diversifications: strongly ancient, ancient and very recent, and moderately ancient (D, D 0 , and D 00 , respectively, in fig. 2B). Finally, our approach permits identifying a range of recent diversifications: strongly recent, recent and very ancient, and moderately recent (E, E 0 , and E 00 , respectively, in fig. 2B). We note that, in our study, ES1 varied between very positive and 0. Inversely, ES2 varied between 0 and very negative. In both cases, there hence was no need to take absolute values (except for the ES1 measured on the "Cadotte phylogeny"; see supporting file S4.2). We note also that a small absolute ES1 and a large absolute ES2 correspond to an "ancient" inflexion point of the curve (i.e., the point where the curve intersects the line of constant branching). We therefore explored whether adding the intersection point improved the discriminative and predictive power of the model. We found no such improvement (classification rates and R 2 remained exactly the same). Finally, a model containing only the inflexion point as independent variable yielded an R 2 of only !0.05. We hence do not present results involving intersection points.

Can ELDERness Surfaces Differentiate among Branching
Scenarios Better than Established Metrics Do?
Using ES1 and ES2, calculated from the parameters of the third-order polynomials previously fit to each simulated phylogeny, we asked whether ES1 and ES2 could discriminate between the different phylogeny types simulated with the four different branching patterns. We first conducted a factorial discriminant analysis (FDA) using ELDERness surfaces to predict the four phylogeny types, which included the four levels of species richness and the four phylogeny types (i.e., 160 phylogenies total) using the R packages MASS (v7.3-40) and ade4 (v1.7-2). We then calculated a confusion matrix (cross-validation test) that shows how efficiently FDA can correctly classify simulated phylogenies into their simulated branching pattern type (counting correct vs. incorrect classifications). We compared the rate of good classification with those obtained on the basis of FDAs using either of the commonly used phylogenetic metrics: Faith's PD, mean pairwise distance (MPD), and mean nearest taxon distance (MNTD; see Webb 2002). These metrics were measured using the R package picante (v1.6-2; Kembel et al. 2010). Admittedly, PD, MPD, and MNTD were not intended to measure branching patterns. We hence additionally used Pybus and Harvey's gamma (2000) and Colless's index of phylogenetic tree imbalance (Heard 1992; see also Heard and Cox 2007). Strictly, Colless's index does not measure changes in diversification through time, but the method is nevertheless established in the literature as an approach to describe patterns of diversification (e.g., Heard and Cox 2007). Positive gamma suggests that the branching rate increases through time (i.e., diversification is rather recent); negative gamma suggests that it decreases (i.e., diversification is rather ancient). Colless's index is 0 for a perfectly balanced phylogeny in which all lineages diversify at the same rate, and it is 1 for a maximally imbalanced phylogeny. Gamma and Colless's indices were measured using the R package ape (v3.0-11) and apTreeshape (v1.4-5), respectively. Colless's index was standardized by 2=[(n 2 1)(n 2 2)] as in Heard (1992) with n as the number of terminal taxa. We recall that we created the same branching patterns for each level of species richness, and hence standardization of parameters for species richness was not needed. Finally, we used ANOVA to test the effect of phylogeny shape on ES1 and ES2. A Tukey test was also used to see which phylogeny types had significantly different ES1 and ES2. P values of these ANOVAs, as well as the mean values of ES1 and ES2, are presented in supporting file S5. We also graphically illustrated the discrimination of the different phylogeny types by ES1 and ES2 ( fig. 3) and assessed the covariation between each of the common phylogenetic metrics and the two ELDERness surfaces using correlation tests.

Relationships between ELDERness Surfaces and Productivity of Grassland Communities
Our ultimate aim was to test whether our new descriptors of phylogenetic branching patterns might be able to predict biodiversity-ecosystem functioning relationships. We used a data set of an exceptional long-term experiment on grassland communities and productivity from the Cedar Creek long-term ecological research site (from Tilman and Downing 1994;Cadotte et al. 2009) that uses a given source pool of species to artificially assemble communities. Being designed to test the combined effects of evolutionary production and local coexistence of large or small numbers of species, such experiments also permit testing whether and how ES1 and ES2 of plant communities correlate with their productivity. The data set includes information on annual aboveground productivity from 1996 to 2007 for 151 grassland plots that were maintained with varying community compositions and levels of species richness from 1 to 16 species over that time period. To allow for comparisons with Cadotte et al. (2009) for the same system, we analyzed the data using the phylogeny successfully applied by Cadotte et al. (2009), hereafter referred to as "Cadotte phylogeny." This Cadotte phylogeny was obtained using maximum likelihood analysis of DNA sequences from four genes: matK, rbcl, ITS1, and 5.8s. We transformed it into an ultrametric phylogeny as done in Cadotte et al. (2009) using the R package ape (v3.2). This function estimates the node ages of a tree using a semiparametric method based on penalized likelihood (Sanderson 2002). The branch lengths of this phylogeny from Cadotte et al. (2009) are derived from mean numbers of substitutions (i.e., per site) between 0 (root) and 1 (present) and used as a proxy of time. We also used a more recent phylogeny from Zanne et al. (2014) that is truly dated, referred to as "Zanne phylogeny." This Zanne phylogeny is clearly more informative and advanced than the Cadotte one, and so we leave the analyses based on the Cadotte phylogeny to the supporting files. We calculated ELDERness surfaces for all communities with eight and 16 species, thus resulting in 65 plots; calculating third-order polynomials is nonsensical for the communities with fewer species (here, four species) and may result in bad fit (see supporting file S4.1). With the Cadotte phylogeny, six residual outliers were found in the species communities (see supporting file S4.2 and S4.3). With the Zanne phylogeny, there was no outlier observed. In all cases, we used the age of the root of the phylogeny, including all the eight or 16 species of the "experimental" pool of each of the plots.
We used ELDERness surfaces, species richness (S), Faith's PD, MPD, and MNTD, as well as all measures of functional diversity used in Cadotte et al. (2009), to predict average annual aboveground productivity from 1996 through 2007. We did not try to predict alternative measures that relate the observed productivities to those expected from monocultures, such as net diversity effect size or transgressive overyielding. These alternative measures were not applicable, because monocultures for five of 16 species were missing. The functional explanatory variables used in Cadotte et al. (2009) were the number of functional groups, functional diversity (Petchey and Gaston 2002), functional attribute diversity (Walker et al. 1999), functional diversity from nonmetric multidimensional scaling, presence of C 3 grasses, presence of C 4 grasses, presence of forbs, presence of nitrogen fixers, and the variation coefficients of six individual traits: leaf area, leaf perimeter area ratio, leaf lobiness, specific leaf area, seed weight, and height. We also explained productivity by MPD and MNTD and their standardized version net relatedness index and nearest taxon index using the "taxa shuffle" null model of the package's R picante (Webb 2002) and by Pybus and Harvey's gamma (2000) and the Colless's index (Heard 1992).
We first analyzed the relationship between the independent variables and productivity using the 65 plots, which varied in species richness. We then tested the effect of ES1 and ES2, including as covariables species richness and the best other functional predictor (i.e., presence of nitrogen fixer) to separate the information carried by ES1 and ES2 from the species richness and the best functional predictor. We used the Akaike information criteria corrected for small sample sizes (AICc) and the coefficient of determination (R 2 ) to compare linear models including single and multiple independent variables, ES1 and ES2 (individually or in combination), or either of the above phylogenetic or functional community characteristics. The ES1 and ES2 were often considered combined, because they do not represent two distinct variables measured independently but parameters of a single function. Together, ES1 and ES2 describe the bivariate relationship between time and branching in the same phylogeny of the same species within the same community. We always found the residuals of our simple and multiple regression analyses to approach normality and homoscedasticity and hence used a Gaussian error distribution. All phylogenetic metrics measured on the basis of simulated phylogenies and on the basis of data from Cadotte et al. 2009 are deposited in the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.jd150 (Yguel et al. 2016).

ELDERness Surfaces Correctly Identify Different Branching Types; Established Metrics Do Not
The FDA results suggest a good differentiation of the four different simulated phylogeny types into distinct clusters along the first two axes, which explained 87.63% and 12.37% of the variance, respectively (fig. 3) Colless's index always performed worse than ES1 and ES2. MPD and MNTD together have less discriminative power than ES1 and ES2 (i.e., 78.12% vs. 83.75%, respectively; see tables 1 and S6.8 (tables S2, S5-S8 available online). Overall, ELDERness surfaces outperformed other methods in most conditions. Not surprisingly, species richness, when used alone, consistently failed to predict the different phylogeny types (table 1 and supporting file S6).

Correlation of ELDERness Surfaces with Established Metrics
The two ELDERness surfaces, ES1 and ES2, of simulated phylogenies were poorly correlated with species richness (supporting file S7a). The ES1 was strongly correlated with MPD, Faith's PD, MNTD, and gamma (r p 0:82, 0.55, 0.68, and 20.45; respectively), whereas ES2 was strongly correlated with only MNTD and gamma (r p 0:60 and 20.77; respectively). Overall, ES1 and ES2 correlated positively with Faith's PD, MPD, and MNTD but negatively with gamma (supporting file S7a). We note that ES1 and ES2 were poorly correlated only with each other or with Colless's index (supporting file S7a). Analyzing the data from the Cedar Creek experiments confirms some of these correlations (i.e., ES1 or ES2 vs. gamma, PD, MPD, or MNTD) and identified new ones (i.e., ES1 vs. species richness or Colless's index; supporting files S7b and S7c).

ELDERness Surfaces Predict Grassland Productivity Better than Established Metrics
We compared models predicting productivity of Cedar Creek plots using ES1 and ES2 with different combinations of measures of phylogenetic diversity, species richness, and different types of species functional trait diversity (supporting file S8). Using the Zanne phylogeny (see "Methods"), the two best models included ES1 and the combination of ES1 and ES2 and explained more variation in productivity (14%-16%), compared with the best model with a phylogenetic diversity metric (11%; although AICc did not differ substantially from a model with MNTD: AICc p 758:37 with ES1 or 757.29 with ES1 and ES2 vs. 759.61 with MNTD; see table 2). Gamma and Colless's index explained only 4% of the variance (table 2). A combination of two commonly used parameters, MPD and MNTD, explains less variance than the combination of our ES1 and ES2 (13% vs. 16%; see supporting file S8 for details). The relationship between ES1 and mean productivity was significantly negative ( fig. 4; table 2; see supporting file S8 for details), indicating that the closer the branching rate was to constant, the higher was the productivity. The effect of ES2 alone was not significant, probably because ES2 values were always close to 0 (mean p20:85; SD p 1:21), and hence we do not show the effect of ES2 in figure 4. In- terestingly, the effect of ES1 remained significant even after including the most important trait-and species-based predictors of productivity (i.e., presence of nitrogen fixer and species richness; table S8.3). In fact, both of these predictors explain less variation of productivity than did ES1 or than the combination of ES1 and ES2. Using the Cadotte phylogeny, the results were slightly different, but the conclusion remained the same (tables S8.2, S8.3).

Discussion
Our results from simulated phylogenies and experimental communities show that the ELDERness approach charac-terizes branching patterns in phylogenies better than established measures. Our approach produces ELDERness surfaces, ES1 and ES2, that could distinguish not only simulated phylogenies comprising ancient lineages or recent lineages but also more complicated scenarios with species stemming from both ancient and recent diversifications ( fig. 3; table 1). Most importantly, we found that ELDERness surfaces explained the variation in productivity within a grassland better than current phylogenetic diversity metrics (table 2) and did so independently of and better than species richness and a recognized driver of productivity, the presence of nitrogen fixers (supporting files S7 and S8.3). We showed that the most productive communities had smaller  ELDERness surfaces, which represent more constant branching rate along the phylogeny, than those of less productive communities.

Benefits of ELDERness Surfaces for Characterizing Branching Patterns
Our results indicate that the ELDERness approach may be more powerful than established measures of phylogenetic diversity for describing how branching patterns affect properties of experimental communities. Therefore, a key future direction will be to study whether ELDERness surfaces are stronger correlates of community or ecosystem processes than phylogenetic diversity in other systems as well, notably in naturally assembled communities. ELDERness surfaces also have a useful property of being poorly correlated with species richness (supporting files S7 and S8.3), in contrast to Faith's PD, and hence they do not need to be standardized for species richness (Safi et al. 2011). This means that the correlations between ELDERness surfaces and community properties are likely not due to covariation with species richness and the related sampling of particular species into species-rich communities, contrary to relationships observed for PD (Venail et al. 2015; but see Cadotte 2015). Interestingly, other measures of community-level phylogenetic patterns (e.g., MPD) that are also independent of species richness were not correlated with grassland productivity as strongly as the ELDERness surfaces (table 2 and supporting  file S8). However, while the ELDERness approach in our study often describes the branching patterns in a community better than phylogenetic diversity measures, it also has a clear limitation: LTT parameters cannot be estimated for very species-poor communities. If we have to give a number, we may propose a limit of six to eight species. In our data set with the four-species communities, the fit is less good (supporting file S4.1) and may give unusable ELDERness surfaces. Nonetheless, other phylogenetic diversity measures have also been demonstrated to show problematic behavior in species-poor communities (e.g., Fritz and Rahbek 2012). In addition, results were, to a minor degree, different with a coarsely and more finely dated phylogeny. However, other metrics, such as MPD or MNTD, were even more sensitive to the quality of the phylogeny than ELDERness surface in our case study. In fact, the effects of ES1 or the combination of ES1 and ES2 remained significant whatever the phylogeny, always explaining the highest proportion of variation in productivity (see supporting file S8). Second, we suggest that, with the rapid accumulation of molecular sequence data, detailed and well-dated phylogenies will become increasingly available in the near future, likely resulting in more robust patterns in community phylogenetics.

How Might Different Branching Patterns Explain Variation in Productivity?
In the majority of our analyses, we found that ES1 and ES2 predict productivity better than established metrics of phylogenetic diversity or diversification. We again stress that we do not attempt to explain ecosystem productivity on the basis of only such purely statistical relationships per se, an exercise that ultimately requires knowledge of traits and their precise effects. Instead, the observed relationships permit us to identify links between functioning of ecosystems and the diversification and local assembly of lineages that are consistent or inconsistent with particular mechanisms. Specifically, the observed pattern is not consistent with several of the hypotheses stated in the introduction, namely, increased productivity due to communities composed of species originated from only (i) ancient diversification (e.g., very positive ES1; see fig. 2B) and hence of complementary functional traits or absence of shared enemies, for instance; or (ii) recent diversification (i.e., very negative ES1; see fig. 2B) and hence of niche filling. Rather, the observed pattern is consistent with the hypothesis of increased productivity due to a combination of occupation of distinct niches by older lineages with niche filling by more recent lineages (Ricklefs 2007;McPeek 2008): the most productive communities showed a very steady accumulation of lineages (i.e., ES1 and ES2 were close to 0), rather than a curve with sudden increases and plateaus. We want to stress, however, that this interpretation remains speculative and based on the assumption of niche conservatism and hence the maintenance of traits related to niche occupation (Wiens et al. 2010). To support this interpretation, we would need to reconstruct niches within lineages and, more importantly, to understand the role of evolutionary niche divergence and niche filling for the functioning of present ecosystems (Nuismer and Harmon 2014;Chalmandrier et al. 2015). We admit that any effect of community structure may be due to the effect of particular species represented in a community rather than due to the structure per se. For this study system, however, such sampling effects appear to be unimportant: analyses based on productivity overyielding in Cadotte et al. 2008 are consistent with those based on productivity per se. More importantly, particular species contributing most to productivity will do so due to particular traits, and we have accounted for such traits or trait structures, including nitrogen fixation. ELDERness surfaces performed independently of and even better than these traits in terms of explained variance. This result suggests that the sampling of particular species is unlikely to explain the observed effect of ELDERness surfaces on productivity.
Overall, the ELDERness approach may allow us to link patterns of macroevolutionary branching and local species coexistence to statistically explain ecosystem functioning. First, across entire lineages, recent or ancient branching may correspond to different mechanisms of diversification, such as niche filling (Ricklefs 2007), ecological versus geographic speciation (McPeek 2008), or mass extinctions (Crisp and Cook 2009). Second, branching patterns of lineages may affect the branching patterns maintained in local communities (see also Anacker and Harrison 2012, but more research is needed); and we showed that the latter affect ecosystem productivity. Third, by influencing productivity, local branching patterns influence natural selection, potentially affecting the evolutionary branching of lineages. Such feedbacks between branching patterns and ecosystem functioning would be an essential new insight, and the ELDERness approach might, in the future, help to identify such feedbacks.

Conclusions
The ELDERness approach describes species relatedness patterns within a community in a way that former metrics of phylogenetic diversity do not. We provide for this purpose a ready-and easy-to-use method and script to implement this new approach. ELDERness's in-depth description of the branching patterns represented by a community is important, because it allowed us to identify how the evolutionary history represented in a community affects its productivity. For our study system, ELDERness predicted productivity even independent of species richness or the presence of nitrogen fixers, whereas common phylogenetic diversity measures did not. We stress that particular traits, like nitrogen fixation, remain obviously important to predict productivity, and our aim here was not to replace the use of traits for understanding ecosystem functioning. We also stress that this study is only a first application of the ELDERness approach to one biodiversity experiment. The Cedar Creek experiment seems to be representative of many experiments in showing a distinct effect of species richness on productivity and no independent effect of phylogenetic diversity (Venail et al. 2015; but see Cadotte 2015). Nevertheless, different experiments may give different results. The Cedar Creek experiment is the oldest running experiment, and we do not know whether, in newer experiments that represent more transitive stages, branching patterns already influence productivity. On the other hand, the Cedar Creek experiment and other existing biodiversity/ ecosystem-functioning experiments consider systems under disturbance: mowed grasslands or recently planted forest. We do not know whether, in less disturbed systems, branching patterns still influence productivity. Future applications of the ELDERness approach should hence focus on both newer experiments and more natural undisturbed systems, and they should account also for ecosystem functions other than aboveground productivity of plants.