An evolutionary perspective on leaf economics: phylogenetics of leaf mass per area in vascular plants

In plant leaves, resource use follows a trade-off between rapid resource capture and conservative storage. This “worldwide leaf economics spectrum” consists of a suite of intercorrelated leaf traits, among which leaf mass per area, LMA, is one of the most fundamental as it indicates the cost of leaf construction and light-interception borne by plants. We conducted a broad-scale analysis of the evolutionary history of LMA across a large dataset of 5401 vascular plant species. The phylogenetic signal in LMA displayed low but significant conservatism, that is, leaf economics tended to be more similar among close relatives than expected by chance alone. Models of trait evolution indicated that LMA evolved under weak stabilizing selection. Moreover, results suggest that different optimal phenotypes evolved among large clades within which extremes tended to be selected against. Conservatism in LMA was strongly related to growth form, as were selection intensity and phenotypic evolutionary rates: woody plants showed higher conservatism in relation to stronger stabilizing selection and lower evolutionary rates compared to herbaceous taxa. The evolutionary history of LMA thus paints different evolutionary trajectories of vascular plant species across clades, revealing the coordination of leaf trait evolution with growth forms in response to varying selection regimes.

sites): σ 2 among =16286 vs σ 2 within =1912. In order to limit the inuence of extreme values, we performed comparative analyses on log-transformed LMA values that we hereafter note l.
Information on growth forms and Raunkiaer life forms (Raunkiaer, 1934) was collected from original studies, databases, and electronic and paper ora (see supporting information for references). The main growth forms were graminoid, forb and shrubtree. We distinguished epiphyte, climber, succulent and aquatic (which corresponded to Raunkiaer's helophyte and hydrophyte life forms) growth forms. Ferns and lycophytes (Selaginellaceae and Lycopodiaceae ) formed a distinct group. Leaf habit (deciduous/evergreen) and life-history (annual/perennial) were available for a subset of the species only.
The supertree. Phylogenetic relationships between species were described as an informal supertree at species level based on published phylogenies (Bininda-Emonds, 2004). We  (2006) and Wikström and Kenrick (2001). Within 34 Angiospermae families, we resolved relationships between genera using published family phylogenies (see supporting information for references). In the two large Asteraceae and Poaceae families, we used tribe-level trees to branch genera within tribes as polytomies (nodes with more than two descending daughter clades). Conicting positions of genera across references were resolved by setting a polytomy at higher level. Genera in the unresolved families were branched as polytomies on the family node. Species were branched as polytomies at genus level in all families.

2
The resulting supertree had 1675 ancestral nodes among which 58% were resolved as bifurcations (nodes with two daughter clades). Below the genus level, the percentage of bifurcations among ancestral nodes rose to 75%.
Ages of ancestral nodes and phylogenetic distance. In order to estimate branch lengths, we dated the ancestral nodes of the supertree in a two-step procedure. We rst dated well-identied nodes in the supertree using published ages. A total of 187 (11%) ancestral nodes matched dated nodes from the literature. For clades in Angiospermae, ages were taken from a comprehensive update of divergence times based on the analysis of sequence data (Bell, Soltis, and Soltis, 2010). These estimates were obtained using a relaxed method allowing for heterogeneity in substitution rates along the phylogeny and a number of fossil ages as calibration points. We completed these ages with estimates derived from several sources for non-angiosperm clades (Anderson, Bremer, and Friis, 2005;Bremer, Friis, and Bremer, 2004;Janssen and Bremer, 2004;Wikström, Savolainen, and Chase, 2001). Age consistency was checked across references. In the second step, the tips of the supertree were assigned an age of 0. Remaining unknown ages were calculated as a combination of the two closest known ages, and branch length as the dierence between ages at the two ends. The second step was achieved using the bladj module of the Phylocom software (Webb, Ackerly, and Kembel, 2007). Ages ranged from 535 Myr for the root of the tree, that is the basis of the Tracheophyta clade, to 2.2 Myr for three genera in the Juglandaceae family (Carya, Juglans and Petrocarya ). Node age dates the oldest evolutionary split within the corresponding clade.
We note D the phylogenetic distance matrix across all species, with diagonal elements d ii = 0 and d ij being the estimated age of the most recent common ancestor (MRCA) of species i and j (age of divergence between the lineages i and j).
Models of continuous trait evolution. We present here the models used in the analyses of LMA evolution. These models estimate parameters that describe macro-evolutionary patterns expected from micro-evolutionary processes (Hansen and Martins, 1996). First, we considered the Brownian motion (BM) model, which assumes that the trait evolves by means of genetic drift independently along the branches of the tree and at constant rate σ independently (Freckleton, Harvey, and Pagel, 2002;Hansen and Martins, 1996): during an innitesimal period dt, the variation in the trait value l is dl = σ 2 dt, where the parameter σ measures the magnitude of perturbations or drift (Hansen, 1997). In addition to drift, directional selection produces similar macro-evolutionary patterns when it occurs under environmental conditions uctuating randomly and rapidly compared to evolutionary time (Diniz-Filho, 2001;Hansen and Martins, 1996). In the BM model, the expected variance in the trait value at a given ancestral node k linearly increases with the distance (period) from the root to the node, d rk : V k = σ 2 (d rk ) (Diniz-Filho, 2001;Freckleton, Harvey, and Pagel, 2002). V k also measures the expected covariance between the trait values of any two where l i is the trait value of species i: the older the divergence between those two species, the higher the covariance.
Second, we considered models of LMA evolution by stabilizing selection (Hansen, 1997).
These models conform to evolution following an Orstein-Ulhenbeck (OU) process: the variation of the trait value during an innitesimal period dt sums as the eects of drift, σ 2 dt as in the BM model, and selection towards a phenotypic optimum θ: dl = −α(l − θ) + σ 2 dt. The parameter σ has the same interpretation as in the BM model. The parameter α controls the rate of adaptive evolution to the optimum θ, or the rate at which the past is discounted in the model: a descendant starts out being similar to its ancestor then the similarity exponentially decreases at a rate determined by α (Hansen, 1997). α is analogous to the magnitude of a supposed selective force. If α is large, species tend to adapt rapidly and trait values are weakly scattered around the optima, depending on σ; if α is small, the equilibrium variance among species is large and species tend to drift apart. In model OU, the expected covariance between two species i and j is cov(l i , l j ) = σ 2 2α e −αd ij 1 − e −2α(dr−dij) , which shows the exponential decrease. Note that a null value of α leads to the BM model, which thus appears as 4 a special case of OU models.
The OU model allows dierent selective optima to be specied along one phylogenetic tree (Butler and King, 2004;Hansen, 1997): changes in the selective optimum mimic changes in the selection regime in dierent parts of the phylogeny. After a change in the selecting regime, evolution unfolds independently along each lineage. Here we considered three OU models based on dierent evolutionary scenarii following the approach of Butler & King (2004) (Butler and King, 2004) (see Fig. 5). Model OU1 species one single phenotypic optimum for the whole tree; model OU3 species three dierent optima thanks to a supposed change in the selection that occurred at the divergence at the basis of Spermatophyta and maintained two dierent optima along the Angiospermae and Gymnospermae branches; lastly, model OU5 includes the optima in OU3 and two others along the Monocotyledoneae and Eudicotyledoneae branches. Phenotypic optima were specied so that after divergence, dierent selection regimes occurred along the daughter branches (Fig. 5), instead of supposing that one of them further evolved under the regime before the split as in (Butler and King, 2004).
Models were tted on the complete set of species and on the two sets of herbaceous and woody species separately by pruning the phylogenetic tree accordingly (Fig. S6). We distinguished the major herbaceous (forbs and graminoids) and woody (shrubs and trees) growth forms only in order to ensure sucient numbers in each group. When dealing with OU models, parameter estimation can raise numerical issues (Butler and King, 2004;Hansen, 1997).
We followed Hansen's recommendation (Hansen, 1997) to rst estimate the α parameter and then run the calibration procedure for the other parameters (θ, σ) keeping parameter α xed to its estimate. We compared the models using the Bayesian Information Criterion  We also considered the null model of phylogenetically independent evolution (PI) which assumes that all species evolve independently as if they were placed at equal distance on a star-phylogeny. This model is not properly an evolutionary model in the sense that it ignores phylogenetic relatedness across species. It serves as a null hypothesis to test the existence of phylogenetic structure in the trait data.
Tree-level analysis of phylogenetic signal. We calculated four statistics to quantify the magnitude and test the signicance of the phylogenetic signal at the tree level: Pagel's lambda (Pagel, 1999), the K-statistic (Blomberg, Garland, and Ives, 2003), withâ the estimated trait value at the root of the tree (Garland and Diaz-Uriarte, 1999) and P satisfying the equation PVP = I. A relatively high value of r obs indicates that the phylogenetic tree correctly describes the covariance structure of trait values among the tips. In K's denition, r obs value is divided by the value obtained under the hypothesis of Brownian motion evolution, r BM , in order to allow comparison between studies (Blomberg, Garland, and Ives, 2003). The value of K was tested against the null hypothesis of phylogenetic independence (PI) by permuting the values among the tips a large number of times (Blomberg, Garland, and Ives, 2003). A K less than one indicates lower similarity among close relatives than expected under Brownian motion evolution along the candidate tree.
Finally, we perform an autocorrelation analysis at the tree level using Moran's index, I, which allows to quantify and test similarity in LMA values with respect to the phylogenetic distance between species (Diniz-Filho, 2001;Gittleman and Kot, 1990). Moran's I was calculated as: where the sums are on the n species of the sample, l i is the log-transformed LMA of species i,l is the mean trait value, and the weight ω ij measures the degree of relatedness between species i and j. Weights were proportional to the inverse of the phylogenetic distance d ij between species i and j, namely the age of their MRCA: ω ij = 1 2d ij . The value of I was also tested against the null hypothesis of phylogenetic independence (PI).
We also performed a correlogram analysis calculating I within classes of divergence time chosen to ensure a sucient number of observations within each class. Phylogenetic correlograms indicate how the similarity among species, i.e. interspecic correlation, changes with the period of divergence between them. This temporal similarity pattern theoretically discriminates the BM and OU models (Diniz-Filho, 2001;Hansen and Martins, 1996): similarity is expected to decrease linearly with time since divergence when the trait evolves by random genetic drift only or by drift and rapidly changing directional selection (BM), whereas it decreases exponentially when the trait evolves by drift and stabilizing selection (OU) (Hansen and Martins, 1996). We adjusted linear and exponential ts to the observed patterns in LMA similarity using non-linear least squares and compared the ts using a likelihood ratio test. Similar analyses were performed within growth forms. Phylogenetic correlograms were also calculated using taxonomic levels, which showed consistent patterns compared to the divergence time approach (see Fig. S5).
Clade-level analysis of the phylogenetic signal. We analyzed the phylogenetic signal at clade level using an Analysis Of Traits (AOT Webb, Ackerly, and Kembel, 2007). DW is calculated as the standard deviation of a across daughter clades with respect to the ancestor's value of a: where n i is the number of daughter nodes of ancestral node i, a i its estimated mean and a i,j the estimated mean of daughter node j. The ancestral mean (a) was estimated recursively from the tips' attributes to the root of the tree following Felsenstein's algorithm (Felsenstein, 1985): where b i,j is the branch length between nodes i and j. Signicant divergence widths were detected by testing the statistic against the null hypothesis PI. As for tree-level statistics, the test was conducted by permuting LMA values among the tips (n = 100000 permutations).
P-values were corrected using the Benjamini & Hochberg method (Benjamini and Hochberg, 8 1995) for multiple testing, ensuring a stringent analysis. When the DW of a node was lower than expected at random, the node was characterized as conservative; when the DW value was higher than expected at random, the node was characterized as diversifying. clade size in parentheses).