Consistent role of Quaternary climate change in shaping current plant functional diversity patterns across European plant orders

Current and historical environmental conditions are known to determine jointly contemporary species distributions and richness patterns. However, whether historical dynamics in species distributions and richness translate to functional diversity patterns remains, for the most part, unknown. The geographic patterns of plant functional space size (richness) and packing (dispersion) for six widely distributed orders of European angiosperms were estimated using atlas distribution data and trait information. Then the relative importance of late-Quaternary glacial-interglacial climate change and contemporary environmental factors (climate, productivity, and topography) as determinants of functional diversity of evaluated orders was assesed. Functional diversity patterns of all evaluated orders exhibited prominent glacial-interglacial climate change imprints, complementing the influence of contemporary environmental conditions. The importance of Quaternary glacial-interglacial climate change factors was comparable to that of contemporary environmental factors across evaluated orders. Therefore, high long-term paleoclimate variability has imposed consistent supplementary constraints on functional diversity of multiple plant groups, a legacy that may permeate to ecosystem functioning and resilience. These findings suggest that strong near-future anthropogenic climate change may elicit long-term functional disequilibria in plant functional diversity.


Supporting Information Overview
This document provides supplementary information not contained in the main text of the article "Consistent role of Quaternary climate change in shaping current plant functional diversity patterns across European plant orders." SI Text-1 provides an extended description the used data (distribution and traits), the trait imputation procedure, functional diversity estimates, and the predictors used in the statistical analyses presented in the main text. SI Text-2 describes the association between multivariate metrics of functional diversity with both species richness and univariate measurements of trait richness and dispersion. It also discusses the difference between observed and expected values of functional diversity based on the local species richness.
SI Text-1. Extended methods-Description of the used traits, the trait imputation procedure, functional diversity estimates, and statistical analyses

Distribution data.
This work uses plant species distribution data from Atlas Florae Europaeae (AFE; ref. 1 ), which maps the distribution of European flora on an equal-area mapping unit of ~50×50km (AFE cells) based on the Universal Transverse Mercator (UTM) projection and a Military Grid Reference System (MGRS). We excluded observations in the former Soviet Union due to uncertainties in these records. This study focuses on angiosperms species, contained in six orders: Caryophyllales -1384spp, Brassicales -973spp, Ranunculales -654spp, Saxifragales -352spp, Rosales -334spp, and Malpighiales -73spp. Plant orders included in this study compose the majority of the AFE and represent a pool of independent evolutionary lineages. For consistency with trait data, occurrences of taxonomic units below the species level (e.g., forms, varieties, and subspecies) were collapsed to those of the corresponding species. The distribution and local species richness of each of the evaluated orders in presented as Fig. S1. Figure S1. Species richness of the six evaluated European Angiosperm orders. Values represent the number of individual species per Atlas Florae Europaeae grid-cells (~50×50 km). Dashed line indicated the northern limit of post-glacial refugia for cold tolerance trees as summarized in ref. 2 . Figure generated using R 3 version 3.0.2. <https://www.r-project.org> based on species counts for each AFEgrid cells.

Used traits.
Five ecomorphological traits were used in the analyses presented in this study: specific leaf area (SLA, cm 2 *g -1 ), seed mass (SWT; mg), maximum stem height (Hmax; m), stem/wood density (WD; kg*m -3 ) and growth form (tree, shrub, herbs/forbs, vines/lianas). These represent some of the best traits for predicting the geographic distribution of vegetation types and assessing the mean plant physiological response in a region. Used traits were also selected due to their use to discriminate between distinctive functional strategies among concurring plants [4][5][6] and the link between these traits and the abiotic environment 7-11 ).
The SLA of a plant indicates how much light-capturing area it produces for each gram of leaf tissue. SLA is negatively correlated with leaf thickness, lifespan, toughness and sensitivity to herbivores and positively associated with mass-based measurements of leaf nitrogen content, photosynthetic capacity, plant relative growth rate, all of which are relevant to plant performance 7 . Therefore, SLA reflects a species position along the leaf economics spectrum 5,7,12 . H max represents the balance between the gains from access to light, the cost of structural support (given the disturbance regime) and water transport, and the sensitivity to biomass loss from mechanical disturbances 4,5,13 . H max is also a proxy for other key traits indicating plant growth 13 . SWT represents the balance between offspring number (seed production) and size (seed size), a trade-off often related to age-dependent survival probability 4,14 . SWT is also negatively correlated with dispersal distance (in the case of winddispersed species) and positively associated with seedling survival probability under light or water limited conditions 14,15 . Lastly, WD represents the trade-offs between growth rates, construction costs, and mortality rates 16 . Within this "wood economics spectrum", a species falls along a continuum going from high volumetric growth rates, low construction costs, and high mortality rates, towards low volumetric growth rates, high construction costs, and low mortality rates 9,16 .

Trait imputation
Mean trait values for each continuous trait were initially determined using the LEDA trait database 17 for as many species as possible. Missing trait information was obtained from multiple data sources: SLA 7,18 ; maximum canopy height 10,[17][18][19] , seed mass 14,20 ; woody density 16 ; and growth form 21 .
To "fill-in" missing values, we used a multiple-imputation approach (MICE; refs. 22,23 ). It is a Markov Chain Monte Carlo (MCMC) in which the missing values are replaced by a small set of modelled alternatives (usually <10). The MICE approach is valid in those situations when a multivariate distribution is a reasonable description of the data, as in the case of trait values. MICE is often preferred over other imputation techniques as it is readily available, provides an approximate solution with good properties, and is not problem-specific. The R package mice 24 was used to generate imputed values.

Supporting Information
Ordonez & Svenning -Functional diversity consistency across orders A predictive mean matching (pmm) approach was used to impute missing values. The pmm is a semi-parametric imputation approach and the default method for continuous variables in MICE. The pmm approach is similar to a regression method, except that the model does not generate the imputed values, but rather serves to construct a metric for matching cases with missing data with similar cases with observed data 25,26 . Imputed values are determined by randomly sampling the observed values based on their match to the results of a simulated regression model. In other words, observations whose regression-predicted values are closest to the regression-predicted value for the missing value are used as imputed values. In our case, the implemented regression model is based on the available trait information, and both evolutionary relatedness (genus) and morphological attributes (growth form). The predictive mean matching method ensures that imputed values are plausible and might be more appropriate than a regression method if the normality assumption is violated 27 . Mean trait values for the four evaluated traits are presented in Fig. S2.

Functional diversity estimation
In this study, the functional composition of a site was described using two complementary functional diversity metrics: functional richness 28 and functional dispersion 29 . Range and variance from each of the evaluated traits were computed for comparisons purposes (see SI Text-2). Functional richness (F Rich ) measures the range of the trait spectrum of a species assemblage. This measurement is based on the multivariate trait space filled by an assemblage as measured by the size of a multivariate convex hull 28 . Functional dispersion (F Disp ) measures the functional packing of the species assemblage. This measurement is based on the distance of individual species to the assemblage mean trait composition 29 . Functional diversity metrics were evaluated using all continuous traits. Before the analyses, all continuous traits were log 10 transformed as they are log-normally distributed. Transformed were also standardized (mean=0 and SD=1) to ensure that all traits contribute equally to the functional diversity estimation and that the units used to measure traits have no influence in the functional diversity estimation 30 .

Environmental information.
A total of seven predictors were selected to summarize historical and contemporary environmental conditions in Europe (Table S1 and Fig. S3). Used predictors describe the ecological mechanisms most often consider as the determinants of diversity patterns at continental and global scales (see refs. [31][32][33][34][35] ).
Historical predictors (Table S1) include climatic stability 34 and accessibility to suitable environments 34,35 at the end of the Last Glacial Maximum (LGM; ~21,000 years ago). Historical climatic stability was measured as the annual temperature and precipitation velocities (measured in kilometers×decade -1 ) from the LGM. Velocities were calculated following ref. 36 as the ratio between temporal anomalies (that is the absolute difference between the LGM and current temperature and precipitation) and the spatial gradients (variability in temperature and precipitation during the LGM on a 3×3 cell neighbourhood). Following the approach in ref. 37 , accessibility to postglacial re-colonization from ice-age refugia (measured in kilometers -1 ) was estimated as the sum of inverse distances between an AFE grid cell and regions considered to be suitable for cool-temperate trees during the LGM. Suitability of cool-temperate trees 21,000 years BP was defined, following ref. 38 , as areas where GDD ≥800°C, mean temperature of the coldest month ≥−15°C, and summer precipitation ≥50mm. Historical climate conditions were determined based on three different

Supporting Information
Ordonez & Svenning -Functional diversity consistency across orders climatic models (CCSM, MIROC, and MPI) for the last glacial maximum from the Paleoclimate Modelling Intercomparison Project Phase II (PMIP 2 ; http://pmip2.lsce.ipsl.fr/). Selected contemporary factors describe water-energy dynamics 31,32 and the importance of spatial heterogeneity 33 . These include (Table S1) mean annual temperature, total annual precipitation (mm×year -1 ), Normalized Difference Vegetation Index (NDVI, unitless), and habitat heterogeneity (variability in elevation measured in m, estimated using ref. 39 elevation model). Table S1. Historical and contemporary predictor variables used to explain spatial variation in European functional richness and dispersion. All variables were summarized for AFE grid cell (~50x50km).

Estimation -Data source
Mean annual temperature velocity -[LGM to present] (km*decade -1 ) Calculated using ref. 36 Figure S3. Contemporary and historical environmental predictors used as explanatory variables of functional richness and dispersion patterns of six orders of European angiosperms. From top to bottom and left to right: mean annual temperature, annual precipitation, NDVI, topographic heterogeneity, temperature and precipitation velocity, and accessibility to glacial refugia. Red polygon shows the maximum extent of the ice sheet 21,000 years BP. Dashed line indicated the northern limit of post-glacial refugia for cold tolerance trees as summarized in ref. 2 . Figure generated using R 3 version 3.0.2 <https://www.r-project.org> using information in the sources described in Table S1.
Spatially corrected pairwise correlations between variables indicated no clear pattern of the correlations between variables (Fig. S4). However, there was a significant negative correlation between historical and contemporary factors. Figure S4. Pearson correlated coefficients between environmental variables used in this study. Significance of the correlation after accounting for spatial autocorrelation determined using Dutilleu-correction 41 . Significance levels: ***=P < 0.001; **=P < 0.01; *=P < 0.05; blank =not significant.

Supporting Information
Ordonez & Svenning -Functional diversity consistency across orders SI Text-2. Extended results-Association between multivariate metrics of functional diversity with both species richness and single trait range and variance.
Across all evaluated orders, as the number of species in grid accumulates the functional richness (F Rich ) increases Fig. S5. This continuous increase is contrasting to similar analyses based on all tree species in Europe 42 , where F Rich decelerated at richness levels below 30 species. By comparison, F Disp association with species richness show multiple shapes across evaluated orders, but the most common shape of these relations is a unimodal trend (Fig. S5). For most of orders, F Disp of those species-rich sites converges to the regional dispersion. Residuals of functional richness/dispersion vs. species richness loess regressions exhibited distinct fine scale variability (Fig. S6). However, regression residuals for both F Rich and F Disp across all orders showed a consistent south to north gradient, and a predominance of negative values (blue to light-green). This indicates a predominance of lower than expected F Rich and F Disp based on the number of species in a location. The strength of the relation of F Rich with trait range and of F Disp with trait variability was similar across evaluated traits. This pattern is consistent across most of the evaluated orders (Fig. S7). Variation in F Rich was strongly and positively associated with the local trait range (richness) of canopy height and seed weight for almost all of the evaluated orders (Fig. S7). Similarly, the variance in wood density, WD and SWT had the strongest effects on F Disp for almost all of the evaluated orders (Fig. S7). The association of F Rich with H max and SWT, and of F Disp with WD and SWT indicate that the trait space and its packing is maximized in areas with a wide array of dispersal strategies and growth forms. Figure S7. Scatter plots of the relation between functional richness and the trait range, and functional dispersion and trait variance for six orders of European angiosperms included in the Atlas Florae Europaeae. Curves are unimodal OLS regressions between species trait range and dispersion and functional diversity.