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.

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 intra-specific 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 (Reich et al., 2003; Adler et al., 2014) and parameterisation of dynamic vegetation models (Fyllas et al., 2012; 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 (La) variations reflect some aspects of the whole plant leaf energy and pan class="Chemical">water balances; leaf dry mass per area (pan class="Chemical">LMA) and leaf nutrients concentration variations may reflect contrasting resource allocation strategies; seed mass (Sm) variations reflect seedling survival and colonisation trade-offs; maximum plant height (Hmax) 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 n class="Species">angiosperms under favourable conditions (Reich et al., 2003). 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 (Brodribb et al., 2012). 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 n class="Chemical">pan class="Species">angiospermsn> 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 pan>tterns of trait covariation are similar within PFTs, and (3) interpret the effects of environmental variation on functional trait variation at the PFT and the species level.

TABLE 1

Functional traits and environmental variables abbreviations and units of measurement used in this study.

Functional trait	Abbreviation	Unit
Leaf area	La	cm2
Leaf dry mass per area	LMA	g m–2
Leaf dry matter content	LDMC	g g–1
Leaf thickness	Lt	mm
Leaf C concentration	Cm	mg g–1
Leaf N concentration	Nm	mg g–1
Leaf P concentration	Pm	mg g–1
Leaf Ca concentration	Cam	mg g–1
Leaf Mg concentration	Mgm	mg g–1
Leaf K concentration	Km	mg g–1
Light saturated photosynthetic rate on area basis	Asat,a	mmol m–2 s–1
Dark respiration rate on area basis	Rdark,a	mmol m–2 s–1
Wood density	ρW	g cm–3
Seed mass	Sm	g
Maximum tree height	Hmax	m
Environmental variable
Average annual temperature	TA	°C
Average monthly temperature	Ti	°C
Minimum temperature of the coldest month	Tmin	°C
Annual precipitation	PA	mm
Monthly precipitation	Pi	mm
Driest quarter precipitation	Pdq	–
Soil sand content	Sand	%
Soil clay content	Clay	%
Soil pH	pH
Soil electric conductivity	EC	μs cm–1
Soil organic matter	SOM	%
Soil N concentration	N	%
Soil P concentration	P	%
Soil K concentration	K	cmol kg–1
Soil Ca concentration	Ca	cmol kg–1
Soil Mg concentration	Mg	cmol kg–1
Water holding capacity	WHC	%

Materials and Methods

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). Sites sampled covered an altitudinal range from 374 m above sea level (asl) at Mount Olympos to 1665 m asl at Mount n class="Chemical">pan class="Chemical">Pan>rnassos, with mean annual temperatures (TA) varying from 5.9 (1516 m asl at Mount Rodopi) to 15.0°C (470 m asl at Mount pan>n class="Chemical">Parnitha), and with mean annual precipitations (PA) ranging from 0.37 (470 m asl at Mount Parnitha) to 1.1 m (1358 m asl at Mount Pindos) (Supplementary Table S1). On each of the 15 mountains, we established at least three 30 m × 30 m plots at different elevations: ‘low,’ ‘medium,’ and ‘high.’ In some mountains were many forest and soil types were abundant we established more than three plots, for example at Mounts Rodopi, Pindos, Olympos, and Kissavos. In each plot, diameters at breast (1.3 m) height (D) 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 n class="Chemical">pan class="Species">Abies cephalonican> Loudon, pan>n class="Species">Abies borisii-regis 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 pan class="Chemical">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. Gas-exchange 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 (Asat,a in μmol pan class="Chemical">CO2 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 (Rdark,a in μmol pan class="Chemical">CO2 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 Asat,a and Rdark,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 Asat,a and Rdark,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 pan class="Chemical">water-saturated leaf fresh mass (Ws in g) was measured. Laminar leaf thickness (Lt in mm) was measured with a digital calliper, the leaves were scanned with a portable scanner and the projected leaf area (La in pan class="Disease">cm2) 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 (Wd in g) determined. The leaf dry matter content (LDMC in g g–1) was subsequently determined as the ratio of Ws/Wd. Leaf dry mass per unit area (LMA in g m–2) was estimated as the ratio of Wd to La. The mass-based photosynthesis (Asat,m) and dark respiration (Rdark,m) rates were calculated by dividing the area-based 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 pan class="Chemical">HCl and with Ca and pan class="Chemical">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 Nm) with N and P additionally expressed and used in statistical analyses on an area basis (Na, Pa) 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 (Hmax in m) was approximated using the 0.99 quantile from the tree-by-tree measurements across the plot network (Supplementary Table S2). Seed mass (Sm 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 km2) climate data were extracted from the CHELSA database (Karger et al., 2017), including average monthly (Ti) and annual temperature (TA), total monthly (Pi) and annual precipitation (PA), and total precipitation during the driest quarter of the year (pan class="Disease">Pdq).

In the lab, composite soil samples (from the five pan class="Disease">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 pan class="Chemical">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 pan class="Chemical">HNO3 and pan class="Chemical">H2O2 (Jones et al., 1991).

The exchangeable cations K, Ca, and n class="Chemical">pan class="Chemical">Mgn> were extracted with pan>n class="Chemical">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 CO2 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 (log10 transformed data, apart from ρW), followed by pan class="Disease">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 pan class="Species">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 Hmax and Sm, 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 Hmax and Sm 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 log10 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 (n class="Chemical">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 comn class="Chemical">parison 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 (Tmin) and the total precipitation of the driest quarter (pan class="Disease">Pdq). 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 pan class="Disease">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 (n class="Chemical">pan class="Species">Abies cephalonican>, pan>n class="Species">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).

Results

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 n class="Chemical">pan class="Disease">Cmn> content and ρW and positively with Nm. Needleleaf evergreen and broadleaf deciduous species were perfectly sepan>rated with a small overlap between evergreen and deciduous broadleaves (Figure 1). This analysis provided good support for our a priori PFT definition.

FIGURE 1

Linear discriminant analysis on the full trait dataset indicating an effective separation of PFTs based on their 10 leaf and one wood trait values. Colours indicate different plant functional types (PFTs) (Ne: needleleaf evergreens, Be: broadleaf evergreens, and Bd broadleaf deciduous). Trait abbreviations: La, leaf area; LMA, leaf dry mass per area; LDMC, leaf dry matter content; Lt, leaf thickness; Nm – Pm – Cam – Mgm – Km leaf, N, P, Ca, Mg, and K mass basis concentrations, Asat,a, light saturated photosynthetic rate on area basis; Rdark,a, dark respiration rate on area basis and ρw wood density. See Table 1 for units.

Linear discriminant analysis on the full trait dataset indicating an effective separation of PFTs based on their 10 leaf and one wood trait values. Colours indicate different plant functional types (PFTs) (Ne: needleleaf evergreens, Be: broadleaf evergreens, and Bd broadleaf deciduous). Trait abbreviations: La, leaf area; n class="Chemical">LMA, leaf dry mass per area; LDMC, leaf dry matter content; Lt, leaf thickness; Nm – Pm – Cam – Mgm – Km leaf, N, P, Ca, Mg, and K mass basis concentrations, Asat,a, light saturated photosynthetic rate on area basis; Rdark,a, dark respiration rate on area basis and ρw wood density. See Table 1 for units.

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 La, ρW and Asat,m than either of the broadleaf PFTs, Bd had the highest Nm, Pm, and n class="Chemical">Mgm concentrations and the highest Asat,m and Rdark,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 Nm, Pm, Mgm, Asat,a, and Rdark,a but higher LMA, Lt, Cm and ρW than Bd. Although there was no difference between the three PFTs in terms of their mean Asat,a, when comparing within PFT there were substantial differences between species evident (for example A. borisii-regis vs. P. halepensis and C. avellana vs. C. coggygria).

FIGURE 2

Mean trait values for the 24 dominant species of the MEDIT plot network (see also Supplementary Table S3). See Table 1 for abbreviations and units.

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 n class="Chemical">pan class="Chemical">LMAn>, pan>n class="Chemical">LDMC, and Cm 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 Hmax and Rdark,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).

TABLE 2

Principal components analysis on 15 traits expressing whole-plant economics, aggregated at species level, for the 24 most dominant species.

	PC1	PC2	PC3
Eigenvalue	6.56	2.51	1.87
Portion of variance	43.74	16.76	12.46
La	0.76	0.21	0.46
LMA	−0.95	0.12	0.16
LDMC	−0.91	0.20	0.12
Lt	−0.68	−0.27	0.33
Cm	−0.88	0.00	0.25
Nm	0.92	0.19	0.06
Pm	0.74	0.59	−0.07
Cam	0.50	−0.37	−0.67
Mgm	0.87	0.08	0.10
Km	0.49	−0.45	0.33
Asat,a	0.41	−0.28	0.40
Rdark,a	−0.05	0.53	0.02
ρW	0.03	−0.73	0.38
Sm	0.35	0.31	0.76
Hmax	−0.34	0.80	0.01

FIGURE 3

Principal components analyses (first two axes) on the average species traits dataset (15 traits), for the 24 most dominant species of the MEDIT plot network. Colors indicate different PFTs (Ne: needleleaf evergreens, Be: broadleaf evergreens, and Bd: broadleaf deciduous). See Table 1 for abbreviations and units.

Princin class="Chemical">pal components analysis on 15 traits expressing whole-plant economics, aggregated at species level, for the 24 most dominant species.

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 La, Nm, Pm, and n class="Chemical">Mgm and negatively to LMA, Lt, and Cm (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 Asat,a and Rdark,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.

FIGURE 4

Principal components analyses of 12 foliar and 1 wood traits, across PFTs (A) within needleleaf evergreens (B), within broadleaf deciduous (C), and within broadleaf evergreens (D) (see also Supplementary Table S6). See Table 1 for abbreviations and units.

Princin class="Chemical">pal components analyses of 12 foliar and 1 wood traits, across PFTs (A) within needleleaf evergreens (B), within broadleaf deciduous (C), and within broadleaf evergreens (D) (see also Supplementary Table S6). See Table 1 for abbreviations and units.

Considering then each PFT separately, for Ne the first PC (23%), represented a needle size dimension, with La covarying with pan class="Chemical">LMA, Lt, and Asat,a (Figure 4B). The second conifer axis (17%) was associated positively to pan class="Chemical">LDMC and ρW and negatively to Nm and Pm. 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 Nm and Pm.

For Bd, the first PC (22%) was positively related to n class="Chemical">pan class="Chemical">LMAn>, Lt and Rdark,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 Nm, Pm, and Km) with leaves of high pan>n class="Chemical">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 pan class="Chemical">LMA, pan class="Chemical">LDMC, and ρW and negatively to La and Mgm. 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 Asat,a and Rdark,a and negatively to Cam, 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 Lt – n class="Chemical">pan class="Chemical">LMAn> (Figure 5A), pan>n class="Chemical">LDMC – LMA (Figure 5B) and LDMC – ρW (Figure 5C). Similarly, leaf Nm and Pm 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 Nm and Pm and between Mgm and Cam (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 – La (Figure 5H) and Lt – LDMC (Figure 5I) in the full dataset, negative relationships were revealed when considered within PFTs. Furthermore sign-differences emerged even between PFTs. For example La – LMA scaled positively within Ne, showed no association within the Bd, and had a negative association within Be (Figure 5J). La – Cm 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).

FIGURE 5

Bivariate relationships among functional traits in Mediterranean Forests (A–L). 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.

Although all three PFTs showed similar positive slopes in their Na – pan class="Chemical">LMA and Pa – pan class="Chemical">LMA relationships (Figures 6A,B), at any given LMA, both Na and Pa 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 Asat,a and both LMA and Na 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 Rdark – LMA association (Figure 6G). Also of note in Figure 6 is that, despite significant associations being observed for the Bd Asat,a – Pa, Rdark,a – Na and Rdark,a – Pa 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 Pa – Na 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 Asat,a vs. LMA relationship (Figure 6D) appears negative, even though for both the Ne and Bd groupings the within PFT-association is clearly positive.

FIGURE 6

Bivariate relationships between gas exchange rates (Asat,a and Rdark,a) and leaf dry mass per area (LMA), nitrogen (Na) and phosphorus (Pa) 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.

Bivariate relationships between gas exchange rates (Asat,a and Rdark,a) and leaf dry mass per area (n class="Chemical">pan class="Chemical">LMAn>), pan>n class="Chemical">nitrogen (Na) and pan class="Chemical">phosphorus (Pa) 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.

Trait Variation as Influenced by Plot Location, Plant Functional Type, and Species

Of the 13 studied traits the least variable was leaf n class="Chemical">pan class="Disease">Cmn> content (CV = 0.039) while the most variable was La (CV = 1.219) followed by pan>n class="Chemical">LMA, Lt and Rdark,a (Supplementary Table S8). 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 La (0.85), Nm (0.77), LMA (0.75), Lt (0.74) and ρW (0.60), most of the variation was attributed to the PFT grouping. By contrast, the environmental component was greater for Km (0.50), Ldmc (0.37) and Rdark,a (0.25), suggesting that these traits are also considerably influenced by sampling location. For Cm (0.42 + 0.07), Pm (0.38 + 0.06), Mgm (0.46 + 0.10) and Asat,a (0.09 + 0.20), the Φ + S component was higher than the environmental component, while the variation of Cam 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 Cam, Asat,a, and Rdark,a.

FIGURE 7

Partitioning of the total variance for leaf and wood functional traits. Traits are sorted based on the portion of variance attributed to PFT. See Table 1 for abbreviations and units.

n class="Chemical">Partitioning of the total variance for leaf and wood functional traits. Traits are sorted based on the portion of variance attributed to PFT. See Table 1 for abbreviations and units.

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 class="Chemical">N, P, and Ca concentrations, thus generally reflecting ‘Soil n class="Chemical">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 (Tmin, pan class="Disease">Pdq, 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 pan class="Disease">Pdq) were associated with higher LMA, LDMC, Cm, Km, and ρ and lower La. Higher Tmin was negatively related to Pm. Soil nutrient status seems to be positively associated with higher ρW and lower Pm. Sandier soils were also inferred to lead to increased Cm but lower Km. 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 plot data (Supplementary Figure S3), with the exception of the La – Pdq and the ρW – Pdq and ρW – Soil Nutrient Status association.

TABLE 3

Partial Kendall correlation coefficients between the environmental component of each traits’ variation and the four axes of environmental variation across the MEDIT plot network.

	Tmin	Pdq	Soil nutrient status	Soil texture
La	0.123	0.343	−0.100	−0.078
LMA	0.116	−0.220	0.043	0.076
Lt	0.044	0.073	0.144	−0.025
LDMC	0.145	−0.323	−0.053	0.170
Cm	0.056	−0.258	0.039	0.194
Nm	−0.038	0.179	−0.102	0.110
Pm	−0.395	−0.045	−0.245	−0.013
Cam	0.085	0.167	0.073	−0.034
Mgm	0.017	−0.061	−0.105	0.181
Km	−0.131	−0.231	−0.038	−0.247
Asat,a	0.077	0.010	−0.065	−0.017
Rdark,a	0.026	−0.158	0.079	−0.061
ρW	−0.026	−0.214	0.227	0.060

n class="Chemical">Partial Kendall correlation coefficients between the environmental component of each traits’ variation and the four axes of environmental variation across the MEDIT plot network.

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 (La, n class="Chemical">pan class="Chemical">LMAn>, Lt, pan>n class="Chemical">LDMC, Cam, and Mgm) 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 La responded differently along Tmin variation, LMA, Lt and Cam responded differently along variations in soil texture and LDMC and Mgm varied individualistic along soil the nutrient status gradient. On the other hand Cm, Nm, Pm, Km, Asat,a, Rdark,a, and ρW were best modelled with just random intercepts for plots and species.

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.

In terms of the fixed effects, an increase in LAI positively affected La and negatively n class="Chemical">pan class="Chemical">LMAn> and Rdark,a. An increase in Tmin was positively associated with Km and higher pan>n class="Disease">Pdq had a positive effect on Lt and a negative effect on LDMC and Cm. Soil nutrients availability had a positive effect on Mgm and a negative effect on Km and finally, coarser soils were associated with higher La and Mgm. 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 species-specific trait-environment relationships which may even trend in opposite directions.

Discussion

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 are considered (Figure 1). This reaffirms the a priori Ne, Be, and Bd classification on MMF when ‘approached’ from a trait-based 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 costs (Reich et al., 1997; Kazakou et al., 2007). Interestingly, Asat,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 Hmax 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. 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 pan class="Chemical">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 La and ρW, and Be a relatively high-leaf cost with a generally high ρW and low Hmax (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 n class="Chemical">LMA and lower Nm and Pm than Bd species, yielding a lower mass-based photosynthetic and respiration rate (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, Lt, Cm were higher and Nm, Pm, Mgm, Asat,m, Rdark,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 (n class="Chemical">pan class="Species">Arbutus unedon>, pan>n class="Species">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 Asat and/or Rdark.

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 Rdark,a contributing stronger than Asat,a though (Figure 4C). Interestingly, as also documented in tropical tree species (Patiño et al., 2012), n class="Chemical">Cam 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 La 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 Asat,a and Rdark,a seem to trade-off with Cam. 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 trade-offs 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 n class="Chemical">LMA and La 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 – La 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 Asat,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 Asat,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 La was the most variable trait followed by n class="Chemical">pan class="Chemical">LMAn>, Rdark,a, and Lt. Various evolutionary, ontogenetic and environmental effects control variation in pan>n class="Chemical">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 (La, LMA, Lt) occurred mainly between PFTs, indicating a certain degree of trait conservatism within these groups (Figure 7). Similar findings were reported for LMA and Nm 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, Cm, Nm, and Pm were the least variable nutrients (de la Riva et al., 2018). Higher variability in Km, Cam, and Mgm 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 n class="Chemical">pan class="Chemical">LMAn>, pan>n class="Chemical">LDMC, Cm, Km, 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 La 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 n class="Chemical">pan class="Chemical">LMAn> of all four species, but at the same time variation in soil texture affected interspecific pan>n class="Chemical">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 Cm it seems that a common response along dry season precipitation can be identified, while for others like Asat,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 whole-plant 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 (La, n class="Chemical">pan class="Chemical">LMAn>, Lt, and Rdark,a being the most variable) and some other traits such as Km, pan>n class="Chemical">LDMC, and Rdark,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, Cm, 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 trait-based 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. n class="Chemical">pan class="Disease">CMn> 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. pan>n class="Disease">PD, MA, and JL subsequently contributed to the original manuscript drafted from NF. All authors commented and approved the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.