Functional Trait Variation Among and Within Species and Plant Functional Types in Mountainous Mediterranean Forests

Plant structural and biochemical traits are frequently used to characterise the life history of plants. Although some common patterns of trait covariation have been identified, recent studies suggest these patterns of covariation may differ with growing location and/or plant functional type (PFT). Mediterranean forest tree/shrub species are often divided into three PFTs based on their leaf habit and form, being classified as either needleleaf evergreen (Ne), broadleaf evergreen (Be), or broadleaf deciduous (Bd). Working across 61 mountainous Mediterranean forest sites of contrasting climate and soil type, we sampled and analysed 626 individuals in order to evaluate differences in key foliage trait covariation as modulated by growing conditions both within and between the Ne, Be, and Bd functional types. We found significant differences between PFTs for most traits. When considered across PFTs and by ignoring intraspecific variation, three independent functional dimensions supporting the Leaf-Height-Seed framework were identified. Some traits illustrated a common scaling relationship across and within PFTs, but others scaled differently when considered across PFTs or even within PFTs. For most traits much of the observed variation was attributable to PFT identity and not to growing location, although for some traits there was a strong environmental component and considerable intraspecific and residual variation. Nevertheless, environmental conditions as related to water availability during the dry season and to a smaller extend to soil nutrient status and soil texture, clearly influenced trait values. When compared across species, about half of the trait-environment relationships were species-specific. Our study highlights the importance of the ecological scale within which trait covariation is considered and suggests that at regional to local scales, common trait-by-trait scaling relationships should be treated with caution. PFT definitions by themselves can potentially be an important predictor variable when inferring one trait from another. These findings have important implications for local scale dynamic vegetation models.

Plant structural and biochemical traits are frequently used to characterise the life history of plants. Although some common patterns of trait covariation have been identified, recent studies suggest these patterns of covariation may differ with growing location and/or plant functional type (PFT). Mediterranean forest tree/shrub species are often divided into three PFTs based on their leaf habit and form, being classified as either needleleaf evergreen (Ne), broadleaf evergreen (Be), or broadleaf deciduous (Bd). Working across 61 mountainous Mediterranean forest sites of contrasting climate and soil type, we sampled and analysed 626 individuals in order to evaluate differences in key foliage trait covariation as modulated by growing conditions both within and between the Ne, Be, and Bd functional types. We found significant differences between PFTs for most traits. When considered across PFTs and by ignoring intraspecific variation, three independent functional dimensions supporting the Leaf-Height-Seed framework were identified. Some traits illustrated a common scaling relationship across and within PFTs, but others scaled differently when considered across PFTs or even within PFTs. For most traits much of the observed variation was attributable to PFT identity and not to growing location, although for some traits there was a strong environmental component and considerable intraspecific and residual variation. Nevertheless, environmental conditions as related to water availability during the dry season and to a smaller extend to soil nutrient status and soil texture, clearly influenced trait values. When compared across species, about half of the trait-environment relationships were species-specific. Our study highlights the importance of the ecological scale within which trait covariation is considered and suggests that at regional to local scales, common trait-by-trait scaling relationships should be treated with caution. PFT definitions by themselves can potentially be an important predictor variable when inferring one trait from another. These findings have important implications for local scale dynamic vegetation models.

INTRODUCTION
In recent years the study of plant functional traits has gained particular interest in ecological and ecophysiological research. This interest arises from the idea that functional traits can provide a stable basis for re-expressing fundamental ecological processes from first principles (McGill et al., 2006): this even leading to the suggestion that there may be some sort of ecological equivalence to the elemental periodic table (Winemiller et al., 2015). Recently developed global plant traits databases (Kattge et al., 2011), and analyses of trait co-variation across wide geographical scales (Reich et al., 1997;Díaz et al., 2016), have provided valuable inroads toward that objective. Nevertheless, considerable intraspecific variation at local scales and environmentally induced variation at larger scales have been repeatedly observed, with their implications still not fully explored (Siefert et al., 2015).
Comparative studies of functional trait variation across species, or plant functional types (PFTs), provide one basis for the identification of life history strategies Adler et al., 2014) and parameterisation of dynamic vegetation models Atkin et al., 2015). Some common leaf, wood, and seed traits are assumed to reveal the way plants acquire resources, reproduce, and compete with other plants (Westoby et al., 2002). For example, leaf area (L a ) variations reflect some aspects of the whole plant leaf energy and water balances; leaf dry mass per area (LMA) and leaf nutrients concentration variations may reflect contrasting resource allocation strategies; seed mass (S m ) variations reflect seedling survival and colonisation tradeoffs; maximum plant height (H max ) variations are indicative of a plant's ability to capture light and disperse seeds; and wood density (ρ W ) variations broadly reflect a trade-off between growth and mortality (Díaz et al., 2016). Although these general dimensions of trait variation that identify fundamental plant strategies have been observed globally, recent studies suggest that at local scales these relationships may not be robust (Messier et al., 2017) and with trait-by-trait scaling relationships differing between sites characterised by different species and/or growing conditions (Schrodt et al., 2015;Lira-Martins et al., 2019) or within species (Anderegg et al., 2018). This means that generic equations that predict one trait from another may not be possible across a broad spectrum of scales. Functional trait variation can be related to species taxonomy as well as to the environmental conditions a particular individual is growing under (Fyllas et al., 2009;Garnier et al., 2016). Geographic gradients, where a set of species is repeatedly found under different conditions, provide natural laboratories for exploring the relative effects of taxonomy and environmental plasticity on trait variation (Asner et al., 2014;Turnbull et al., 2016).
Mountainous Mediterranean forests (MMF) occur across extended mountain ranges, such as those of the Grecian Peninsula, where forest type frequently covaries with elevation from relatively dry Mediterranean to more temperate profiles (Box and Fujiwara, 2015). Three basic woody PFTs dominate these mountains: needleleaf evergreen (Ne), broadleaf deciduous (Bd) and sclerophyllous broadleaf evergreen (Be). These PFTs are usually a priori defined based on their leaf form and habit with the Ne and Be species typically having higher leaf longevity than Bd. Nevertheless an evaluation of other functional traits is necessary to better understand and model the underlying basis of their differential distributions across the Mediterranean region (Carnicer et al., 2013), or their contribution across communities of different successional stage and/or disturbance history (José Vidal-Macua et al., 2017). Understanding how these PFTs interact with each other and why a particular PFT may dominate under particular environmental conditions is important in order to predict vegetation dynamics under global change at both the regional and planetary scale.
Needleleaf evergreen species are mainly found in relatively adverse environments and disturbed habitats and are considered to be able to survive under extreme conditions due to their relatively high cavitation resistance and nutrient use efficiency (Bond, 1989;Brodribb et al., 2012). The high cavitation resistance of Ne species is related to the wide hydraulic safety margin under which they operate (Choat et al., 2012) and thus a generally lower hydraulic conductivity, which may, in turn, reduce their competitive ability against angiosperms under favourable conditions . However, recent studies suggest that single hydraulic traits should not be used to explain differences in whole plant hydraulic strategy (Anderegg, 2015;Gleason et al., 2016), that PFT classification might not adequately capture the impacts of drought on tree mortality (Anderegg, 2015) and that the severity of drought might be a stronger predictor of tree mortality compared to PFT grouping (Greenwood et al., 2017). Ne species tend to be superior colonisers in disturbed sites but they are also able to tolerate low disturbance regimes . Under conditions of high and stable water and nutrient availabilities, broadleaf deciduous species are considered to generally outcompete Ne (Berendse and Scheffer, 2009). This is thought to be due to their higher hydraulic conductivity, lower LMA, and higher photosynthetic capacity which place them toward the 'acquisitive' part of the leaf economic spectrum (Wright et al., 2004). Broadleaf evergreen species of MMF on the other hand, prevail at drier conditions (Box and Fujiwara, 2015) by deploying long-lived schlerophyllous leaves, with high construction costs (Williams et al., 1989) and high LMA that can survive extended dry periods and maintain photosynthesis through and beyond periods of extended soil water deficits (Givnish, 2002;De Micco and Aronne, 2012).
Recent ecophysiological studies suggest that the contrasting trait syndromes between angiosperms and gymnosperms, can lead to different responses to increased temperature and drought (Carnicer et al., 2013). It is thus important to understand how functional traits express 'plant strategies, ' and if trait covariation differentiates between PFTs. In this study we therefore systematically measured 12 leaf morphological and biochemical traits plus wood density across a range of species and environmental conditions all along a forest plot network on Mediterranean mountains in Greece (Table 1). In addition, we estimated species-specific maximum height from tree-by-tree biometric measurements made within each plot, and we extended the trait database with seed mass information from the literature. We aimed to: (1) test whether PFT definition in MMF correspond to different suites of functional traits, (2) test if the patterns of trait

Study Plots and Species
The study extended across 61 plots of the MEDIT network (Fyllas et al., 2017b) and covered the most important mountains of continental Greece in terms of both species diversity and ecosystem productivity (Supplementary Figure S1) were determined for all individuals with D ≥ 10 mm. We also estimated the height (H) of more than 50% of the individuals in each plot, ensuring to cover the full range of observed D. Five soil pits to 0.3 m depth were also dug in each plot and all soil, below the litter horizon, was extracted for subsequent analysis of physical and chemical properties. Plot-level leaf area index (LAI) was estimated using an ACCUPAR 2000 as the average values of 20 measurements made at 1 m above the forest floor. Within our plots we sampled individuals from 39 tree and shrub species. Individuals from the 24 dominant species (defined as those species that contributed at least 5% to the total stand basal area) were used for further analysis. Each species was assigned to one of the three PFTs already discussed, viz needleleaf evergreens including Abies cephalonica Loudon, Abies borisiiregis Mattf., Picea abies Karst., Pinus halepensis Miller, Pinus nigra Arnott and Pinus sylvestris L., broadleaf evergreens including Arbutus unedo L., Arbutus andrachnae L., Quercus coccifera L., Quercus ilex L., and Phillyrea latifolia L. and broadleaf deciduous including Acer campestre L., Betula pendula Roth, Carpinus orientalis Mill., Castanea sativa Mill., Corylus avellana L., Cotynus coggygria Scop., Fagus sylvatica L., Fraxinus ornus L., Ostrya carpinifolia Scop., Pistacia terebinthus L., Quercus cerris L., Quercus frainetto Ten., and Quercus pubescens Willd.

Functional Trait Measurements
Within each plot (apart from the monospecific stands) at least 10 individuals were selected for functional trait measurements (Supplementary Table S12). The number of individuals per species was selected based on the relative contribution of each species to the stand's basal area. One, fully sunlit branch from mature individuals was cut by climbing on the tree and/or using telescopic scissors, and immediately placed in a water bucket where it was recut prior to leaf gas-exchange measurements. For species found only in the understory, sunlit individuals outside the plot were sampled. Whilst still attached to the recut branch, healthy fully expanded leaves were selected and placed within the gasket of a LICOR-6400 infrared gas analyser (LI-COR, Lincoln, NE, United States). Only current year's leaves were used. Gasexchange was monitored and when leaves had reached a stable photosynthetic rate, with a stomatal conductance higher than 0.05 µmol s −1 m −2 at an incident photon irradiance (I) of 1500 µmol quanta m −2 s −1 , the area-based light saturated net photosynthetic rate (A sat,a in µmol CO 2 s −1 m −2 ) was recorded as the mean value of five measurements per leaf, made across 3 s intervals. The area-based leaf respiration rate (R dark,a in µmol CO 2 s −1 m −2 ) was estimated from the average value of five (with 3 s intervals) measurements made on leaves that were placed for at least 5 min in the dark. All measurements were made with a chamber temperature near 25 • C and a relative humidity between 50 and 70%. The average temperature of the chamber during the measurements was recorded, and subsequently both A sat,a and R dark,a were re-expressed at a common temperature of 25 • C, using the equations from Tjoelker et al. (2001) and Higgins et al. (2016) respectively.
In addition to the leaves used in the A sat,a and R dark,a measurements, at least two, fully developed leaves of the same branch were placed in sealed bags with moist tissue paper and left in dark conditions for 24 h before their water-saturated leaf fresh mass (W s in g) was measured. Laminar leaf thickness (L t in mm) was measured with a digital calliper, the leaves were scanned with a portable scanner and the projected leaf area (L a in cm 2 ) was estimated using the image analysis software Image-J (NHI, version 1.47). Once back in the laboratory, leaves were dried at 80 • C for 48 h and their dry weight (W d in g) determined. The leaf dry matter content (LDMC in g g −1 ) was subsequently determined as the ratio of W s /W d . Leaf dry mass per unit area (LMA in g m −2 ) was estimated as the ratio of W d to L a . The mass-based photosynthesis (A sat,m ) and dark respiration (R dark,m ) rates were calculated by dividing the areabased rates with LMA. For the determination of wood density (ρ W in g cm −3 ), a piece of each cut branch was transferred to the lab, where its dry weight (at 70 • C for 48 h) and volume, estimated via the water displacement method, were measured (Williamson and Wiemann, 2010).
For determinations of leaf cations and P composition, 0.5 g of ground leaf material was heated at 450 • C for 5 h, in 1N HCl and with Ca and Mg concentrations subsequently determined using an atomic absorption spectrophotometer, K with a Corning 410 flame photometer, and P by the vanade-molybdate method (Jones et al., 1991). Plant N content (1 g samples) was determined by the Kjeldahl wet-oxidation method (Bremner and Mulvaney, 1982). Total C and S were determined by a LECO CNS 2000 analyser (TruSpec Micro, St. Joseph, MI, United States). All leaf nutrient concentrations are expressed on a per mass basis (in mg g −1 , denoted hereafter with a 'm' subscript, for example N m ) with N and P additionally expressed and used in statistical analyses on an area basis (N a , P a ) due to the obvious relevance of area based metrics when looking at photosynthesis-nutrient associations (Lloyd et al., 2013).
For some analyses, the functional trait dataset of leaf and wood traits, was complemented by estimates of maximum species height and seed mass. Maximum species height (H max in m) was approximated using the 0.99 quantile from the tree-bytree measurements across the plot network (Supplementary Table S2). Seed mass (S m in g) data for all species were extracted from the Seed Information Database (Royal Botanic Gardens Kew, 2018).

Climatic Data and Edaphic Properties
For each plot, long-term high resolution (∼1 km 2 ) climate data were extracted from the CHELSA database (Karger et al., 2017), including average monthly (T i ) and annual temperature (T A ), total monthly (P i ) and annual precipitation (P A ), and total precipitation during the driest quarter of the year (P dq ).
In the lab, composite soil samples (from the five pits) were air dried, crushed and 2 mm sieved, prior to determination of soil particle size by the hydrometer method (Bouyoucos, 1962). Soil pH and electrical conductivity were also estimated in a suspension of 1:1 water:soil (Doran et al., 1996), organic C determined by the Walkley-Black wet oxidation method (Nelson and Sommers, 1982), and total N by Kjeldahl wet-oxidation (Bremner and Mulvaney, 1982). Soil total P was determined by wet-acid digestion with HNO 3 and H 2 O 2 (Jones et al., 1991).
The exchangeable cations K, Ca, and Mg were extracted with 1N ammonium acetate at pH 7 (Thomas, 1982), with K concentrations subsequently measured by Corning 410 flame photometer, and Ca and Mg, by Varian AA400 Plus atomic absorption. Carbonate content was determined using the Bernard method by measuring the evolved CO 2 after addition of HCl (Nelson, 1982). Maximum water holding capacity (WHC) was measured for each soil sample (Gardner, 1986): each soil sample was saturated with water in a cylinder, and WHC was calculated based on the weight of the water held in the sample vs. the sample dry mass (dried at 105 • C for 24 h). By measuring WHC without taking into account stone and/or rock content, we refer to the intrinsic ability of mineral soils to hold water regulated mainly by pedogenetic factors that determines soil type and soil organic matter dynamics.

Statistical Analysis
All statistical analyses and figures were made with R (R Core Team, 2019). Initially, a linear discriminant analysis (LDA, package MASS) was performed on the full trait dataset (including intraspecific variation) to verify that functional traits (predictor variables) can be efficiently used to define PFTs (response variable). Analyses of variance (log 10 transformed data, apart from ρ W ), followed by Tukey HSD post hoc tests were used to explore for differences among mean trait values between the three PFTs, with species treated as a random effect (packages lme4, lsmeans) and including intraspecific trait variation by using the full dataset. Within angiosperms we additionally applied a phylogenetic ANOVA (package phytools; Revell, 2012) to test for trait differences between Bd and Be species, using in this case the across site mean trait values per species. The latest GBOTB tree was used to take into account the phylogenetic history of the study species (package V.PhyloMaker; Jin and Qian, 2019). The correlation matrix of the per species average trait dataset (no intraspecific variation), extended with the species-specific mean H max and S m , was analysed with a Principal Components Analysis (PCA, package FactoMineR), to identify the major functional dimensions. Further PCAs were performed, at the full dataset (including intraspecific variation) as well as after aggregating the full dataset into PFTs (Ne: n = 128, Bd: n = 181, Be: n = 35), so as to explore whether the major functional dimensions differed between PFTs. In these analyses H max and S m were excluded as data was not available for each measured individual. For all PCAs we estimated the radius of the equilibrium circle of descriptors (Legendre and Legendre, 1998) to assess the contribution of each trait to each principal component.
For both the full-dataset, as well as testing separately within each PFT, Pearson's correlation coefficients were estimated for all trait pairs. Analyses were performed on log 10 transformed traits values, with the exception of ρ W . For significantly correlated trait pairs, standardised major axis regressions (SMA, package smart) were fitted in order to test whether the scaling relationships between traits were similar across PFTs.
Each trait's variability was estimated by the coefficient of variation (CV). In addition multilevel linear models (package lme4) were used to quantify the sources of trait variation (Fyllas et al., 2009) that account for: (a) between PFT variation ( ), (b) interspecific variation (S), (c) variation among regions (R), i.e., between mountains (this reflecting plastic and/or filtering responses to the wide environmental gradient (climatic and edaphic) of our plot network), and (d) between plot variation (P), this reflecting natural environmental variability of plant growing conditions. Our multilevel model can be written as: with µ, the overall mean value for each trait (T), and ε the residual term, which includes both within-species variability, as well as any measurement error. As only sun-leaves were collected, we expect micro-environmental effects to be minor in comparison to taxonomic and environmental effects.
After fitting the multilevel model for each trait, the derived components were extracted for further analysis. In particular the environmental effect on each trait's variation was estimated by adding the regional and plot (E = R + P) components. The E component expresses the value a trait would take in each plot after removing the effects of PFTs and species, revealing the 'true' effect of environmental variability (Fyllas et al., 2009). Climatic variability between plots was expressed by variation in the minimum temperature of the coldest month (T min ) and the total precipitation of the driest quarter (P dq ). These two variables were used to express the two main climate factors limiting plant growth in Mediterranean plants during the year, i.e., winter cold and summer drought (Mitrakos, 1980). In order to identify the key axes of edaphic variability across our plot network a PCA on the correlation matrix of the soil variables was performed, and the scores of the plots on the first two axes were used as edaphic predictors for subsequent analyses. We note that our soil measurements were made at the top 30 cm and might not be appropriate for species with deep roots, but can in general be considered proxies for soil fertility and water retention ability (Ordoñez et al., 2009). The effects of local growing conditions (climatic and edaphic) on each trait's environmental component were tested using Kendall partial correlation analysis (package ppcorr).
Finally, for a subset of four species (Abies cephalonica, Pinus nigra, Fagus sylvatica, and Quercus frainetto), with functional trait measurements in at least five individuals of these species across a minimum of five plots, linear mixed effect models were used to explore whether the effects of environmental variability on trait variation were independent of species identity. Trait values were z-score standardised, based on the trait mean value and standard deviation, making effect sizes comparable across traits. In this analysis, LAI (index of light availability) and the two climatic and two edaphic gradients were used as fixed effects, with species and plot used as random effects. For model selection and validation we started with models that used all fixed effect terms, and searched for the optimal random structure, by systematically allowing both intercepts and slopes to vary for each species across the climatic and edaphic gradients of the plot network (Zuur et al., 2009). All models were fitted using REML and the model with the lowest AIC was selected. Subsequently we searched for the optimal fixed components by sequentially removing the non-significant fixed-effect terms (package lmerTest).

Plant Functional Types in MMF
Linear discriminant analysis found the first axis (LD1) to explain 95% of the total variance, with LD1 being negatively correlated with leaf C m content and ρ W and positively with N m . Needleleaf evergreen and broadleaf deciduous species were perfectly separated with a small overlap between evergreen and deciduous broadleaves (Figure 1). This analysis provided good support for our a priori PFT definition.
Differences both within and between PFT for most of the studied traits were identified (Figure 2 and Supplementary Tables S3, S4). For example, Ne had lower mean L a , ρ W and A sat,m than either of the broadleaf PFTs, Bd had the highest N m , P m , and Mg m concentrations and the highest A sat,m and R dark,m and Be had the highest ρ W . The differences between the two broadleaf PFTs were maintained even when their phylogenetic history was considered (Supplementary Table S5), with Be having a lower N m , P m , Mg m , A sat,a , and R dark,a but higher LMA, L t , C m and ρ W than Bd. Although there was no difference between the three PFTs in terms of their mean A sat,a , when Colours indicate different plant functional types (PFTs) (Ne: needleleaf evergreens, Be: broadleaf evergreens, and Bd broadleaf deciduous). Trait abbreviations: L a , leaf area; LMA, leaf dry mass per area; LDMC, leaf dry matter content; L t , leaf thickness; N m -P m -Ca m -Mg m -K m leaf, N, P, Ca, Mg, and K mass basis concentrations, A sat,a , light saturated photosynthetic rate on area basis; R dark,a , dark respiration rate on area basis and ρ w wood density. See Table 1 for units.
Frontiers in Plant Science | www.frontiersin.org comparing within PFT there were substantial differences between species evident (for example A. borisii-regis vs. P. halepensis and C. avellana vs. C. coggygria).
A PCA on the species means dataset, identified three functional dimensions ( Table 2). The first dimension (PC1 -44% of total trait variance) can be considered to describe a leaf dimension that contrasts thick leaves with high LMA, LDMC, and C m with large nutrient rich (N, P, Mg) leaves: separating species across a leaf construction cost dimension. The second dimension (PC2 -17%) is positively related with H max and R dark,a , and negatively with ρ W , and could be considered to reflect a trade-off between height gain and persistent life strategy. The third axis (12% of total variance) was mainly related to seed mass. The 24 studied species occupied distinct areas of the multidimensional trait space (Figure 3): Ne occupying the high leaf construction cost (high PC1 scores), Bd having low construction cost/high nutrient (dry weight basis) leaves (low PC1 scores), while Be seem to adopt an overall conservative tissue construction strategy (low PC1) coupled with small adult stature and high ρ W (low PC2 scores).
Further probing for PFTs differences through an assessment of the trait inter-relationships within each PFT (i.e., treating each tree as an separate observation rather than using species means as above), the first two axes of the full dataset PCA explained 64% of the total variance (Supplementary Table S6). PC1 (45%) was strongly positively related to L a , N m , P m , and Mg m and negatively to LMA, L t , and C m (Figure 4A). This is similar to the first dimension identified in the species level analysis above. On the other hand PC2 (12%) was in this case mainly related to A sat,a and R dark,a , suggesting that leaf gas exchange, expressed on an area basis, is largely independent to leaf resource allocation. PC3 (10%) was associated to LDMC and ρ W , representing a tissue toughness dimension. Considering then each PFT separately, for Ne the first PC (23%), represented a needle size dimension, with L a covarying with LMA, L t , and A sat,a ( Figure 4B). The second conifer axis (17%) was associated positively to LDMC and ρ W and negatively to N m and P m . Thus in Ne species leaf construction and gas exchange seems to be independent from N and P concentration, at least when the nutrients are expressed on a dry-weight basis. On the other hand, tissue density (expressed by LDMC and ρ W ) are both negatively associated with N m and P m .
For Bd, the first PC (22%) was positively related to LMA, L t and R dark,a , indicating that thicker leaves have a higher per area maintenance cost ( Figure 4C) and similar to Ne, the second Bd dimension (21%) contrasts leaves of high nutrient investment (mainly N m , P m , and K m ) with leaves of high LDMC and ρ W .
Although the PCA results for the Ne and Bd were broadly similar, for Be a different pattern was observed with the first dimension (30%) positively related to LMA, LDMC, and ρ W and negatively to L a and Mg m . This indicates that within this PFT leaf and wood construction traits integrate along a common axis ( Figure 4D). The second axis (15%) is positively related to A sat,a and R dark,a and negatively to Ca m , suggesting for this PFT that variations in photosynthetic capacity and respiration are independent of leaf construction costs.

Differences in Bivariate Trait Relationships Between PFTs
Numerous significant bivariate trait relationships were identified both within and across PFTs (Supplementary Figure S2 and Supplementary Table S7). In many cases the sign of the association was similar between PFTs, for example the positive relationships between L t -LMA (Figure 5A), LDMC -LMA ( Figure 5B) and LDMC -ρ W (Figure 5C). Similarly, leaf N m and P m concentration scaled negatively with LMA for all groups (Figures 5D,E). Only in few cases a common slope relationship was identified, for example between leaf N m and P m and between Mg m and Ca m (Figures 5F,G). However, in most cases (78%) the common slope test indicated a significant difference in the scaling exponent between PFTs, suggesting that scaling relationships depend on PFT for most of the bivariate trait associations examined. Interestingly, when data across PFTs were pooled, some of the relationships become of an opposite sign than when the different PFTs are considered separately. For example, although a significant positive association was identified between ρ W -L a ( Figure 5H) and L t -LDMC ( Figure 5I) in the full dataset, negative relationships were revealed when considered within PFTs. Furthermore sign-differences emerged even between PFTs. For example L a -LMA scaled positively within Ne, showed no association within the Bd, and had a negative association within Be (Figure 5J). L a -C m scaled positively within Bd and negatively within Be species (Figure 5K). LMA scaled negatively with ρ W within Ne, and positively within Be and Bd (Figure 5L).
Although all three PFTs showed similar positive slopes in their N a -LMA and P a -LMA relationships (Figures 6A,B), at any given LMA, both N a and P a were higher for Bd than either Ne or Be (Supplementary Table S7). A contrast between Ne and Bd was also evident for the relationships between A sat,a and both LMA and N a where the bivariate associations differed in terms of slope and intercept respectively (Figures 6D,E), although in both cases no significant relationship was observed for Be. This was also the case for the R dark -LMA association ( Figure 6G). Also of note in Figure 6 is that, despite significant associations being observed for the Bd A sat,a -P a , R dark,a -N a and R dark,a -P a associations, no significant relationships were found for these three bivariate associations for either Ne or Be. Also, despite most of the area-based bivariate associations of Figure 6 being PFT dependent, for the P a -N a association all three PFTs essentially fall along the same line. Finally we note that, as for some of the relationships in Figure 6, when the data are pooled (without consideration of PFT) the slope of the A sat,a vs. LMA relationship ( Figure 6D) appears negative, even though for both the Ne and Bd groupings the within PFT-association is clearly positive.

Trait Variation as Influenced by Plot Location, Plant Functional Type, and Species
Of the 13 studied traits the least variable was leaf C m content (CV = 0.039) while the most variable was L a (CV = 1.219) followed by LMA, L t and R dark,a (Supplementary Table S8).  Table 1 for abbreviations and units.
Partitioning this variation according to Eq. 1, for most traits the proportion of the variance attributable to PFT ( ) and species (S) components, surpassed that attributed to environmental conditions (plot and region effects) (Figure 7). For example, for L a (0.85), N m (0.77), LMA (0.75), L t (0.74) and ρ W (0.60), most of the variation was attributed to the PFT grouping. By contrast, the environmental component was greater for K m (0.50), L dmc (0.37) and R dark,a (0.25), suggesting that these traits are also considerably influenced by sampling location. For C m (0.42 + 0.07), P m (0.38 + 0.06), Mg m (0.46 + 0.10) and A sat,a (0.09 + 0.20), the + S component was higher than the environmental component, while the variation of Ca m was equally attributable between environment (0.24) and taxonomy (0.29). We note that for most traits there was a significant within species variation and error term, particularly high for Ca m , A sat,a , and R dark,a .

Environmental Effects
Two main axes of edaphic variation across our plot network were identified. The first PCA axis (48%) was positively correlated with WHC, pH, SOM, N, P, and Ca concentrations, thus generally reflecting 'Soil Nutrient Status' variations, while the second was mainly related to 'Soil Texture' (20%) being positively associated with coarser textured soils (Supplementary Table S9).
We then estimated the effects of the two climatic and the two edaphic (T min , P dq , Soil Nutrient Status, and Soil Texture) gradients on the environmental component of each measured trait using partial Kendall's τ ( Table 3). This shows that drier conditions (lower P dq ) were associated with higher LMA, LDMC, C m , K m , and ρ W and lower L a . Higher T min was negatively related to P m . Soil nutrient status seems to be positively associated with higher ρ W and lower P m . Sandier soils were also inferred to lead to increased C m but lower K m . We note that these trait-environment relationships express the cross-species, 'pure' environmental-driven trait variation, since the effect of leaf habit and taxonomy were removed from this analysis. The same partial correlation analysis was performed using the raw plot level average trait values (Supplementary Table S10). Most of the significant partial correlations were common between the environmental component and the average  Table S6). See Table 1 for abbreviations and units. Figure S3), with the exception of the L a -P dq and the ρ W -P dq and ρ W -Soil Nutrient Status association.

plot data (Supplementary
A summary of the linear mixed effect model analysis for the four most widely measured species in our study is presented in Supplementary Table S11, with the inferred relationships shown in Figure 8. For almost half of the traits (L a , LMA, L t , LDMC, Ca m , and Mg m ) the lowest AIC model included random intercepts associated with plot identity (meaning that even after accounting for the independent covariate that there were significant systematic effects of sampling plot on mean values of the trait under investigation) with random intercepts and slopes for species along one of the environmental gradient (suggesting that different species may respond to environmental variations in fundamentally different ways). In particular L a responded differently along T min variation, LMA, L t and Ca m responded differently along variations in soil texture and LDMC and Mg m varied individualistic along soil the nutrient status gradient. On the other hand C m , N m , P m , K m , A sat,a , R dark,a , and ρ W were best modelled with just random intercepts for plots and species.
In terms of the fixed effects, an increase in LAI positively affected L a and negatively LMA and R dark,a . An increase in T min was positively associated with K m and higher P dq had a positive effect on L t and a negative effect on LDMC and C m . Soil nutrients availability had a positive effect on Mg m and a negative effect on K m and finally, coarser soils were associated with higher L a and Mg m . Slope estimates among the four species along the main climatic and edaphic gradients differed, in some cases even in sign (Figure 8), suggesting that although a common trait response to environmental variability can in some cases be identified (see previous paragraph), that there are also speciesspecific trait-environment relationships which may even trend in opposite directions.

Plant Functional Types in Mountainous Mediterranean Forests
Our results demonstrate a clear grouping of the studied species to PFTs based on leaf habit, when species average trait values . Colours indicate individuals' PFT. When a significant relationship was identified a SMA fit is shown in the respective colour, broken black lines indicate significant relationships in the full dataset (see Supplementary Table S7 for coefficient estimates). The LR tests indicate significant differences between the slope of the PFT specific SMA lines. See Table 1 for abbreviations and units.
are considered (Figure 1). This reaffirms the a priori Ne, Be, and Bd classification on MMF when 'approached' from a traitbased perspective. Our multivariate analysis of species average trait values further identified three functional dimensions in support of the Leaf-Height-Seed (LHS) framework (Westoby, 1998). The first axis (Figure 3 and Table 2) expresses a leaf economic spectrum and contrasts species with cheap leaf construction against species with expensive leaf construction FIGURE 6 | Bivariate relationships between gas exchange rates (A sat,a and R dark,a ) and leaf dry mass per area (LMA), nitrogen (N a ) and phosphorus (P a ) area content (A-I). Colours indicate individuals' PFT. When a significant relationship was identified a SMA fit is shown in the respective colour, broken black lines indicate significant relationships in the full dataset (see Supplementary Table S7 for coefficient estimates). The LR tests indicate significant differences between the slope of the PFT specific SMA lines. See Table 1 for abbreviations and units.
costs (Reich et al., 1997;Kazakou et al., 2007). Interestingly, A sat,a is not associated with cheaper construction-cost leaves in agreement to the notion that variations in leaf photosynthetic capacity may occur more or less independently of variations in construction cost/leaf longevity (Lloyd et al., 2013), as is also suggested by subsequent studies that have found leaf structure, gas exchange and hydraulic traits to be effectively decoupled (Li et al., 2015;Costa-Saura et al., 2016;Rosas et al., 2019).
The second functional dimension expresses a maximum adult stature vs. wood density trade-off, that contrasts tall species (that also have higher leaf maintenance cost) against shorter species with a persistent life style. Taller species are better competitors for light (Westoby et al., 2002) and H max has been associated with higher growth rates in both tropical (Poorter et al., 2008) and Mediterranean (Martínez-Vilalta et al., 2010) forests, at least in terms of wood volume increment.  Table 1 for abbreviations and units. Bold values indicate statistically significant associations (p < 0.05), after controlling for the effect of all other environmental dimensions, and italic indicate marginally significant associations (p < 0.1). Trait abbreviations: L a , leaf area; LMA, leaf dry mass per area; LDMC, leaf dry matter content; L t , leaf thickness; N m -P m -Ca m -Mg m -K m , leaf N, P, Ca, Mg, and K mass basis concentrations; A sat,a , light saturated photosynthetic rate on area basis; R dark,a , dark respiration rate on area basis; ρ w , wood density.
On the other hand, higher wood density is related to lower growth rates (Enquist et al., 1999), higher survival/endurance (Martínez-Vilalta et al., 2010) and resistance to drought-induced xylem cavitation (Hacke et al., 2001). The third dimension mainly relates to seed mass, expressing a trade-off between seed production and establishment rate (Westoby et al., 2002), and interestingly with leaf calcium concentrations showing a negative association with this third trait dimension, as was also found for tropical trees in the Amazon Basin study of Patiño et al. (2012). In this trait space, Ne display high leaf-construction costs and adult stature, Bd low leaf-construction costs along a range of L a and ρ W , and Be a relatively high-leaf cost with a generally high ρ W and low H max (Figure 3). Overall our analysis supports the existence of the LHS framework in MMF. Nevertheless, this conclusion contrasts with the study of de la Riva et al. (2016), who, working across 38 Mediterranean woody species found no evidence for the LHS orthogonal dimensions. This discrepancy may arise from the taxonomically wider set of species measured in our study (including needleleaf evergreens) allowing for more pronounced inter-species contrasts.

Differences Between PFTs
For most of the studied foliar properties, Ne were found at the conservative region of the trait spectrum (Figure 3 and Supplementary Table S4), in agreement to the 'slow seedling' hypothesis that has been used to explain Ne exclusion from faster growing broadleaved species in productive habitats (Bond, 1989;Brodribb et al., 2012). For example Ne species had higher mean LMA and lower N m and P m than Bd species, yielding a lower mass-based photosynthetic and respiration rate FIGURE 8 | Linear mixed effect models for the measured functional traits, across the major axes of environmental variation identified along the MEDIT forest plot network for the four best-studied species. For traits that regression lines are presented, the analysis suggested that along the respective axis of environmental variation the optimum random structure required different slope for each species (Supplementary Table S11). In cases with no regression lines the optimum random structure required only varying intercepts.
(with differences eliminated when expressed on an area basis) as also reported in other studies (Lusk et al., 2003). These photosynthetic differences might be attributable to the higher hydraulic capacity and stomatal conductance of angiosperms that enables them to sustain higher transpiration rates (Lusk et al., 2003;Brodribb et al., 2005). As expected, Ne were characterised by the lowest ρ W between the three PFTs due to their simpler wood structure consisting mainly from tracheids in contrast to the more complex angiosperm wood structure (Zhang et al., 2017). The wood anatomical differences are considered part of a suite of traits that form two generic hydraulic strategies with gymnosperms operating on safer hydraulic margins and having higher cavitation resistance and lower xylem recovery capacity in contrast to angiosperms (Choat et al., 2012;Carnicer et al., 2013). However recent studies show PFT classifications may not capture hydraulic differences (Anderegg, 2015) with the severity of drought being a stronger predictor of tree mortality than PFT grouping (Greenwood et al., 2017). Among the studied angiosperms a continuum of 'fast vs. slow' plant strategies (Reich, 2014) seems to emerge, with more conservative trait values observed for Be compared to Bd species. For example, LMA, L t , C m were higher and N m , P m , Mg m , A sat,m , R dark,m , and ρ W were significantly lower in the studied Be species compared to their Bd counterparts. In tropical and temperate forests a potential coordination of gas exchange and hydraulic architecture has been reported (Brodribb and Feild, 2000;Maherali et al., 2006;Zhu et al., 2013), with ρ W found to negatively correlate with hydraulic conductivity and photosynthetic rates (Santiago et al., 2004;Hoeber et al., 2014).
In addition ρ W has been proposed as a proxy for cavitation vulnerability with denser wood species showing higher cavitation resistance (Hacke et al., 2001;Rosner, 2017). The above seems to agree with the more conservative and drought resistant strategy followed by Mediterranean Be species that also seems to explain their occurrence at drier environmental conditions (Costa-Saura et al., 2016).
Most of the broadleaf evergreen species (Arbutus unedo, Arbutus andrachnae, Quercus coccifera, Quercus ilex, and Phillyrea latifolia) and some of the broadleaf deciduous species (Fraxinus ornus, Ostrya carpinifolia, and Pistacia terebinthus) in this study are commonly found in the understory of MMF. Functional traits such as LMA are known to be sensitive to variation in light availability, with higher irradiance leading to higher LMA (Poorter et al., 2009). Intraspecific LMA variability between sunlit and shaded leaves has been shown for all three PFTs (Wyka et al., 2012;Gratani, 2014) with lower values recorded in shaded leaves and potentially leading to downregulation of gas exchange rates (Chen et al., 2014). In our study, all understory species that had no individual at the sunlit part of the canopy were sampled for trait measurements outside the plot, in sunnier places in order to reduce such environmental effects and enable between species and sites comparisons. However such LMA values will not be representative of the conditions experienced by individuals in the understory of our stands, and could for example lead to overestimation of A sat and/or R dark .
Our analysis found evidence for differences and similarities in multi-trait coordination when the three PFTs were considered separately (Figure 4 and Supplementary Table S6). In Ne the first functional dimension ( Figure 4B) reveals a needle economic spectrum, with bigger needles characterised by a higher dry mass per area, thickness and light saturated photosynthetic rate (Oren et al., 1986). The first Bd dimension is similar to that of Ne, with R dark,a contributing stronger than A sat,a though ( Figure 4C). Interestingly, as also documented in tropical tree species (Patiño et al., 2012), Ca m is negatively associated with PC1 for broadleaf deciduous species, suggesting a leaf construction cost dimension. In particular leaves with low LMA and high mineral content might emerge from thinner, less lignified cell walls and potentially associated with higher levels of organic acid (Poorter and de Jong, 1999). The second Ne and Bd dimensions integrate leaf and wood traits, with species of higher leaf and wood tissue density (LDMC and ρ W ) characterised by lower leaf nutrient concentrations, and highlight LDMC as a better indicator of resource capture and use strategy than LMA (Wilson et al., 1999;Hodgson et al., 2011). This resource use dimension seems to explain community dynamics with high (leaf and wood) tissue density (conservative) species exhibiting higher resistance to physical damage, higher drought tolerance and survival compared to low tissue density species (Markesteijn et al., 2011;Lasky et al., 2014). In Be species on the other hand, an integrated leaf and wood dimension emerged in PC1 ( Figure 4D) indicating trait converge in a previously documented plant economic spectrum (de la Riva et al., 2016). At the 'fast turnover' end of the spectrum there are plants with relatively low LMA, LDMC, ρ W , and high L a and nutrient/cations contents in contrast to the 'slow turnover' end. However as in the across PFTs analysis, the fast-slow turnover dimension identified in the broadleaf evergreen type seems to be independent to the gas exchange axis identified in PC2, where A sat,a and R dark,a seem to trade-off with Ca m . The different pathways leaf and wood economics are realised between PFTs, could potentially relate to the different leaf habit/longevity and/or environmental conditions under which they grow. For example the spatial scale of the analysis might be important, with de la Riva et al. (2016) showing that coordination between functional traits becomes weaker or disappears when considering species belonging to environmentally similar conditions. Thus our results suggest that trait covariance patterns depend strongly on the unit of organisation probed (Anderegg et al., 2018).
Trait inter-correlations are useful to identify functional tradeoffs and plant strategies, while at the same time are frequently used in dynamic vegetation models to infer one functional trait from another (Sakschewski et al., 2015). The next generation of dynamic vegetation models represent functional diversity within traits' spectra, rather than mean species or mean PFT values (Scheiter et al., 2013;Fyllas et al., 2014), and thus universal scaling relationships are particularly useful. Nevertheless, recent studies suggest that many trait covariances may not hold at local spatial scales (Messier et al., 2017;Anderegg et al., 2018). We found differences in some key trait-pair scaling relationships between PFTs in MMF (Supplementary Table S7). For example we found a negative relationship among LMA and L a in the full dataset and within Be, but with positive relationship found for Ne ( Figure 5J). In the Ne case this may be attributable to bigger needles being also thicker, with previous studies showing that the strength of the LMA -L a association was lost when considered within gymnosperms (Ackerly and Reich, 1999). In our study ρ W scaled positively with LMA across angiosperms (both Bd and Be), suggesting a potential coordination of leaf and wood traits (Santiago et al., 2004;Chave et al., 2009;Patiño et al., 2012), in accordance with the coordination of these two traits within Fagaceae (Vilà-Cabrera et al., 2015). Across PFTs though a negative association was found (Figure 5L), in agreement with the findings of Vilà-Cabrera et al. (2015) that included both members from Pinaceae and Fagaceae. However, a consistent positive LDMC -ρ W relationship was observed across all species and PFTs (Figure 5B), highlighting the use of LDMC as a potentially better surrogate for component leaf traits than LMA (Richardson et al., 2013). Interestingly, the negative relationship among A sat,a and LMA across PFTs is actually positive when considered within Bd and Ne species separately ( Figure 6D). This shift could be explained by the contrasting sources of LMA variation that could lead to different relationships between A sat,a (Osnas et al., 2018). Thus across PFTs, where LMA variation relates to structural toughness and leaf longevity, a negative association is expected. These differences highlight that the use of global scaling relationships could be problematic when parameterising dynamic vegetation models, particularly at regional and local scales. We therefore suggest that PFT specific parameterisations need to be developed so as to better represent trait covariation relationships that are usually embedded in such models (Fyllas et al., 2017a).

Functional Trait Variation Along Environmental Gradients
Across our dataset L a was the most variable trait followed by environmental effects control variation in LMA (Poorter et al., 2009), which may be best considered as a part of an integrated trait complex (Poorter et al., 2014). In this study, variation in the three leaf structural traits (L a , LMA, L t ) occurred mainly between PFTs, indicating a certain degree of trait conservatism within these groups (Figure 7). Similar findings were reported for LMA and N m between angiosperms and gymnosperms in Spain (Vilà-Cabrera et al., 2015). Wood density, on the other hand, was rather 'stable, ' although again most of the variation was found at the PFT level (Vilà-Cabrera et al., 2015). In terms of leaf elemental concentrations, C m , N m , and P m were the least variable nutrients (de la Riva et al., 2018). Higher variability in K m , Ca m , and Mg m suggests that these macronutrients are probably scarce across our MMF plots network, in accordance to the 'hypothesis of the stability of the limiting elements' that postulates that limiting nutrients are less variable than more abundant ones (Han et al., 2011). Additionally variation in some nutrients such as Ca and Mg, could also be explained by the high variation in pH and amount of calcareous rocks across our sites, that is known to have a strong influence on the ability of plants to absorb some cations (Goulding, 2016). At the same time, different patterns of variation among scales were observed between the measured macronutrients, indicating different processes of nutrient regulation across our plots network.
For most of the measured traits climatic effects were stronger predictors of environmental variation than edaphic effects ( Table 2). Drier summer conditions lead to increased LMA, LDMC, C m , K m , and ρ W (Supplementary Figure S3), i.e., more conservative leaf and wood deployment. One plant adaption to drier conditions is to reduce the surface area from which they lose water (Givnish, 1987), as also shown in the current and previous studies (Costa-Saura et al., 2016;Rosas et al., 2019). According to Nardini et al. (2014), the reduction in L a could lead to higher vein density and thus higher LMA, offering drought tolerance in Mediterranean plants. It should be noted thought that trait-specific shifts with environmental variables might be better viewed as a coordinated response, especially for traits that have a high number of linkages with other traits, such as LMA (Nardini et al., 2014;Costa-Saura et al., 2016). At the same time, ρ W is known to increase with water deficit across and within species (Preston et al., 2008;Patiño et al., 2012) and even intra-annually within the same individual (Bouriaud et al., 2005). Wood density is related to minimum leaf water potential (Ackerly, 2004;Bucci et al., 2004), cavitation resistance (Hacke et al., 2001) and lower mortality rates (Martínez-Vilalta et al., 2010) and thus plants with high ρ W should, generally speaking, found toward the conservative end of the fast-slow spectrum (Reich, 2014). We also note that all the environmental -trait associations reported here refer to the true environmental effect on trait variability, after removing variation attributed to the taxonomic and/or PFT classification, and thus refer to adaptive trait response to changes in growing conditions.
The trait analysis of the four most common species in our plots network revealed that different environmental variables control intraspecific trait variation, leading to trait-specific and idiosyncratic species responses. For example denser stands (higher LAI) had a common positive effect on the LMA of all four species, but at the same time variation in soil texture affected interspecific LMA variation of Abies cephalonica in a stronger way (Figure 8). In a similar way although an increase in dry season precipitation lead to lower LDMC for all species the effect of soil nutrient availability was much stronger for Quercus frainetto individuals. However, for some other traits like C m it seems that a common response along dry season precipitation can be identified, while for others like A sat,a neither PFT classification nor any of the environmental predictors used in this study could adequately capture intraspecific variation. All the above illustrate that traits may respond individualistically within species across some key environmental gradients, sometimes even when comparing within the one PFT.

CONCLUSION
In this study we verify differences in key functional traits, between the three most abundant PFTs on Mediterranean Mountain Forests. Multivariate analysis of a set of traits, expressing wholeplant economics, support the existence of the three independent axes as suggested by the LHS framework (Westoby, 1998). However, when foliar and wood trait covariation was examined within each PFT, we found different bivariate associations and different functional dimensions, suggesting that trees within each PFT might optimise their coordinated trait responses in alternative ways (Wyka et al., 2012). Also, some traits showed a greater taxonomic variability than others (L a , LMA, L t , and R dark,a being the most variable) and some other traits such as K m , LDMC, and R dark,a were more responsive to environmental variation. Our analysis also shows that drier conditions may lead to more conservative trait syndromes as exemplified by increased LMA, LDMC, C m , and ρ W . However, when explored within populations of the same species, environmental gradients may drive trait variations in different directions for different species. Our findings highlight the effects of source of variation and local environmental conditions on trait values and trait covariation. We thus suggest the use of regional and local data wherever possible when modelling forest function with traitbased approaches.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
NF and JL conceived the study. NF made the gas exchange measurements and undertook all statistical analyses. CM and AG undertook plant sampling and functional trait measurements. EE and CT undertook soil physical and chemical properties measurements. JZ-C helped with gas exchange measurement protocols. PD, MA, and JL subsequently contributed to the original manuscript drafted from NF. All authors commented and approved the manuscript.