Intraspecific variability of leaf form and function across habitat types

Ecology Letters. 2024;27:e14396. | 1 of 15 https://doi.org/10.1111/ele.14396 wileyonlinelibrary.com/journal/ele LETTER Intraspecific variability of leaf form and function across habitat types GiacomoPuglielli1 | AlessandroBricca2 | StefanoChelli3 | FrancescoPetruzzellis4 | Alicia T. R.Acosta5 | GiovanniBacaro4 | EleonoraBeccari6 | LilianaBernardo7 | GianmariaBonari8,9 | RossanoBolpagni10 | FrancescoBoscutti9,11 | GiacomoCalvia2 | GiandiegoCampetella3 | LauraCancellieri12 | RobertoCanullo3 | MicheleCarbognani9 | MartaCarboni5 | Maria LauraCarranza9,13 | Maria BeatriceCastellani14 | DanielaCiccarelli15 | AndreaCoppi14 | MaurizioCutini5 | AliceDallaVecchia10 | MicheleDalleFratte16 | Maria Carlade Francesco9,13 | PieterDeFrenne17 | MicheleDeSanctis18 | Leopoldode Simone8 | ValterDiCecco13 | GiulianoFanelli18 | EmmanueleFarris19 | AriannaFerrara20 | GiuseppeFenu21 | GoffredoFilibeck12 | CristinaGasperini22 | DomenicoGargano7 | ElisabethKindermann2 | GretaLaBella5 | LorenzoLastrucci23 | LorenzoLazzaro14 | SimonaMaccherini8,9 | MichelaMarignani21 | MicheleMugnai9,14 | LuigiNaselli- Flores24 | Nicodemo GiuseppePassalacqua7 | NicolaPavanetto25 | AlessandroPetraglia10 | FrancescoRota2,26 | Lucia AntoniettaSantoianni13 | AldoSchettino27 | FedericoSelvi22 | AngelaStanisci9,13 | GiacomoTrotta4,11 | PieterVangansbeke17 | MarcoVarricchione9,13 | MarcoVuerich9,11 | CamillaWellstein2 | EnricoTordoni6 Received: 27 November 2023 | Revised: 22 January 2024 | Accepted: 19 February 2024 DOI: 10.1111/ele.14396 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2024 The Authors. Ecology Letters published by John Wiley & Sons Ltd. Giacomo Puglielli and Enrico Tordoni contributed equally to this work. For Affiliation refer page on 12 Correspondence Giacomo Puglielli, Departamento de Biología Vegetal y Ecología, Facultad de Biología, Universidad de Sevilla, Sevilla, Spain. Email: [email protected] Funding information Ente Parco Nazionale del Pollino (Rotonda, Italy) in the frame of the project “Un laboratorio naturale permanente nel Parco Nazionale del Pollino”; National Biodiversity Future Center NBFC, CUP J33C22001190001; European Abstract Traitbased ecology has already revealed main independent axes of trait variation defining trait spaces that summarize plant adaptive strategies, but often ignoring intraspecific trait variability (ITV). By using empirical ITV- level data for two independent dimensions of leaf form and function and 167 species across five habitat types (coastal dunes, forests, grasslands, heathlands, wetlands) in the Italian peninsula, we found that ITV: (i) rotated the axes of trait variation that define the trait space; (ii) increased the variance explained by these axes and (iii) affected the functional structure of the target trait space. However, the magnitude of these effects was rather small and depended on the trait and habitat type. Our results reinforce the idea that ITV is contextdependent, calling for careful extrapolations 2 of 15 | INTRASPECIFIC VARIATION OF LEAF STRATEGIES INTRODUCTION Plant functional traits (Violle etal.,2007) provide a successful tool to understand species sorting along environmental gradients and their assemblage in biological communities (de Bello etal.,2021). The use of plant functional traits across spatial scales and levels of biological organization has allowed identifying specialization axes summarizing the mechanical and energetic tradeoffs constraining plant form and function, for both aboveground and belowground components (Bergmann etal.,2020; Carmona etal.,2021; Díaz etal.,2016; Reich etal.,1997; Weigelt etal.,2021; Westoby,1998; Wright et al., 2004). Aboveground, two largely independent axes of trait variation have emerged across species (Díaz etal.,2016). One axis is defined by traits related to the size of the individual plant and its organs, namely, the maximum height of the plant, the area of the leaf and the dry mass of the seed. This axis describes a tradeoff between the plant and seed size, but also a scaling of organ sizerelated traits with total plant size. A second axis is defined by leaf economic traits, i.e., leaf nitrogen content and leaf mass per unit leaf area, mostly reflecting a tradeoff between carbon and nitrogen allocation between structural vs. metabolic leaf components (the leaf economics spectrum, Wright etal.,2004). Despite the generality of these axes of trait variation (Carmona etal.,2021; Joswig etal.,2022), it remains unclear how intraspecific trait variability (ITV) shapes them (but see De Frenne etal.,2011). ITV is a relevant facet of trait diversity because individuallevel trait variability ultimately represents the raw material for natural selection (Bolnick etal.,2011; Des Roches etal.,2018; Palacio etal.,2024; Westerband etal., 2021). However, ITV has been generally overlooked due to its smaller magnitude compared to betweenspecies trait variation (BTV; Violle etal.,2012) as well as for practical reasons (e.g., sampling effort; Puglielli etal.,2022), especially in studies involving many taxa across large spatial scales. As a result, together with a disproportionate focus on trait patterns at the BTV- level (McGill etal.,2006), there is a huge gap in the availability of ITV- level data for many species across different ecosystems. While large trait databases (e.g., TRY, Kattge etal.,2020) have just started filling this data gap and the importance of ITV is increasingly being acknowledged (Siefert etal.,2015), such an endeavour is still in its infancy. Despite data limitations, we know that the importance of ITV tends to increase with decreasing spatial scale (Albert etal.,2011), for species with broader niches (Albert et al.,2010; Sides et al.,2014) and in communities characterized by lower species richness (Siefert etal.,2015; Thomas etal.,2020). Importantly, a global metaanalysis found that ITV within- and among communities accounts for 25–32% of total trait variation and the magnitude of this effect depends on the considered of ITV patterns across traits and spatial scales. Importantly, our study provides a framework that can be used to start integrating ITV into trait space analyses. KEYWORDS functional diversity, functional traits, intraspecific trait variability, leaf area, leaf mass per unit of the leaf area, plant strategies, trait space Union - NextGenerationEU within the framework of National Biodiversity Future Center (Spoke 4, Activity 4); NBFC to the University of Florence, funded by the Italian Ministry of University and Research, PNRR, Missione 4 Componente 2, “Dalla ricerca all'impresa”, Investimento 1.4, Project CN00000033; NBFC to University of Roma Tre/Department of Science, funded by the Italian Ministry of University and Research, PNRR, Missione 4 Componente 2, “Dalla ricerca all'impresa”, Investimento 1.4, Project CN00000033. Grant of Excellence Departments 2018- 2022, MIUR Italy; NBFC to University of Molise/Department of Bioscience and Territory, funded by the Italian Ministry of University and Research, PNRR, Missione 4 Componente 2, “Dalla ricerca all'impresa”, Investimento 1.4, Project CN00000033, MIUR Italy; MCIN/AEI/ 10.13039/501100011033 and by FEDER “ESF Investing in your future”, Grant/Award Number: PID2021- 122214NA- I00; Grant of Excellence Departments 2018- 2022, MIUR Italy; G. Bo. and SM acknowledge the support of NBFC to University of Siena, funded by the Italian Ministry of University and Research, PNRR, Missione 4 Componente 2, ‘Dalla ricerca all', Investimento 1.4, Project CN00000033 Editor: Roberto Salguero- Gomez 14610248, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ele.14396 by Universidad De Sevilla, Wiley Online Library on [31/05/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 3 of 15 PUGLIELLI etal. trait (Siefert et al., 2015). In fact, leaf economic traits showed a higher degree of ITV compared to sizerelated traits, with the former determined to a greater extent by environmental filtering than by species’ evolutionary history (Capdevila etal.,2023; Flores etal.,2014). Therefore, it is reasonable to hypothesise that traits that are more phylogenetically conserved, such as sizerelated traits (e.g., leaf area, Moles etal.,2005; Ávila- Lovera etal.,2023), could display fewer ITV compared to leaf economic traits and thus be more dissimilar between species than within species. However, so far, to the best of our knowledge, this hypothesis has received little attention (Eichenberg etal.,2015). Other recent largescale analyses focused on demonstrating the adaptive value of trait–trait relationships (Niinemets,2015; Zhou etal.,2022) and traitenvironment relationships (Anderegg, 2023; Anderegg etal., 2018; Griffin- Nolan & Sandel,2023) by comparing them at the BTV- vs. ITV- level. These analyses showed that some trait–trait relationships observed at the BTV- level are not conserved within species, spotlighting the potential adaptive value of alternative trait combinations across environments due to ITV. The analysis of bivariate relationships successfully quantifies the effect of ITV on trait–trait coordination and its potential to expand trait space dimensionality (Griffin- Nolan & Sandel, 2023). However, this approach tests one trait dimension at a time, while strategy schemes are usually defined by a set of independent axes of trait variation at the interspecific level (Laughlin,2023; Westoby,1998). Including ITV in such schemes could reshape their axes, but in a traitdependent manner (Siefert et al., 2015), altering the trait space only in some preferential directions. The magnitude of these effects has never been tested so far; thus, the axes of trait variation that define current strategy schemes may overlook a large part of potentially adaptive trait variability. Additionally, trait space distortions due to the inclusion of ITV have direct consequences on functional diversity metrics derived from trait space analyses by modifying the distributions of individual observations in a trait space (Wong & Carmona,2021). Clarifying the role of ITV on functional diversity is highly important considering the prominent role of individuallevel trait variation and functional diversity in maintaining ecosystem functioning. Here, we explore the effect of ITV on two independent trait dimensions of leaf form and function (Díaz etal.,2016). We selected individual leaf area (LA) as the sizerelated dimension and leaf mass per unit area (LMA) as the leaf economic dimension. LA is a key ecological trait because it determines the amount of lightcapturing surface for photosynthesis and influences leaf thermodynamics, the wateruse efficiency, plant biomechanics and vulnerability to herbivory (Moles,2018; Wright etal.,2017). LMA is a measure of the biomass invested to display a unit of leaf area and leaves with a higher LMA typically have a longer lifespan but a lower photosynthetic capacity (Reich etal.,1997; Westoby etal.,2002; Wright etal.,2004). Given their link to multiple important aspects of leaf form and function, the combinations of these traits provide a robust way to summarize leaflevel trait dimensions related to largely independent axes of ecological specialization (Pierce etal.,2013). In addition, from a practical point of view, these traits are widely measured, guaranteeing information from many individuals across species. We assembled an extensive dataset including 3153 standardized individuallevel measurements of LA and LMA for 167 species (49 families, 132 genera, including both woody and herbaceous plants) throughout the Italian peninsula, a region that is not represented in the global metaanalysis by Siefert etal.(2015), but for which there is now a largely untapped availability of ITV trait data (Chelli etal.,2019). Furthermore, despite its relatively small extension, the Italian peninsula encompasses broad environmental gradients both in latitude and elevation and, therefore, several major habitat types such as coastal dunes, Mediterranean scrub, temperate forests, wetlands and grasslands (Chelli et al., 2019). First, we characterized the two dimensions of the trait space defined by LA and LMA in terms of leading variance directions (eigenvectors) and their associated amount of explained variance (eigenvalues). Then, we quantified the role of ITV in influencing the leaf economics (LMA) and size (LA) dimensions. Specifically, we wanted to answer the following questions: • Q1: Does ITV modify the direction of the leading variance components (i.e., eigenvectors) of a trait space? Given that ITV modifies the trait variance of a population, we hypothesised that comparing the trait space defined at the BTV- level with the trait space accounting for ITV would result in a trait space distortion through rotation of its leading trait dimensions (Figure1a). The magnitude of this rotation could be used to evaluate the overall effect of ITV in modifying the directions along which leading trait dimensions could be identified. • Q2: Does ITV increase the variance explained along each trait space dimension? In addition to altering the direction of the trait dimensions in a trait space, accounting for ITV can modify the amount of variance that is captured along a given trait dimension (Figure1b). We expected that a greater ITV could increase the variance captured along a given trait dimension, even if the magnitude of this effect might be dimensiondependent, by altering the trait space along preferential directions. Additionally, we tested whether the trait space axis with a stronger phylogenetic signal showed a weaker ITV signal (see above). • Q3: Does ITV modify the area and functional structure of the trait space? Including ITV in a trait space ultimately modifies how species/individuals are 14610248, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ele.14396 by Universidad De Sevilla, Wiley Online Library on [31/05/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 4 of 15 | INTRASPECIFIC VARIATION OF LEAF STRATEGIES distributed within a trait space (Figure1a,c). This has consequences on different functional diversity metrics calculated on a given trait space (Wong & Carmona, 2021). ITV may influence (i) functional richness, that is, the area of a trait space, (ii) functional divergence, the average distance of data points from the centre of gravity of the trait space and (iii) probabilistic distributions of trait combinations. We hypothesised that accounting for ITV, due to an increased number of combinations of trait value, would consistently expand the underlying trait space leading to higher functional diversity, that is increased richness, divergence and overlapbased dissimilarity (Carmona etal.,2016), compared to the trait space defined at the BTV- level. Observations in our dataset spanned five habitat types: coastal dunes, forests, grasslands, heathlands and wetlands. Siefert et al. (2015) found that the relative extent of ITV showed little variation with climate across global biomes. However, the relative importance of ITV is expected to increase from coarser to finer spatial scales (Albert etal.,2011). Thus, breaking down our analysis at a smaller spatial scale than the global one, that is, at the level of habitat type, has the potential to shed light on whether different environmental conditions potentially affect the ITV patterns. FIGURE 1 Potential effects of including ITV in trait space. (a) In this example, ITV rotates the trait space. Trait space rotation refers to the angle (α) defined by the eigenvectors (Dim.1 and Dim.2) calculated from the trait's covariance matrix at the interspecific level using species medians (grey eigenvectors) and at the intraspecific level (cyan lines) by replacing species medians with the value of the trait sampled from random individuals. In the trait space defined at the interspecific level (BTV), circles of different colours represent different species, while the trait space including ITV is defined by replacing species with individual observations (see Data analysis). (b) ITV can increase the variance captured along a given trait space dimension. This is tested by evaluating the ratio between eigenvalues associated with either Dim.1 or Dim.2 when including or not including ITV to define the trait space. A ratio >1 implies that including ITV increases the variance captured along a given trait space dimension. (c) ITV might alter the underlying functional structure, that is the shape of the space and density of data points within the trait space. These alterations may be captured using functional diversity metrics such as functional richness, overlap/nestedness and functional divergence. LA, leaf area; LMA, leaf mass per unit area. LA LMA BTV(Species medians) ITV α= trait space rotation EigenvectorBTV EigenvectorITV Ratio=1 VTB/VTIeulavnegiE Dim. 1 Dim. 2 Shifts in densityareas Increasedpoint dispersion Increasedarea Dim.2 Dim.1 DensityBTV DensityITV (a) Q1: Does ITVmodify thedirection of theleading variance components (eigenvectors)of atraitspace? (c) Q3: Does ITV modify the areaand the internal structure of a traitspace? LA LMA (b) Q2:Does ITV increase the variance explained along thetraitspace dimensions? LA LMA Increasing variance explained duetoITV Median Individual Compared to BTV trait space 14610248, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ele.14396 by Universidad De Sevilla, Wiley Online Library on [31/05/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 5 of 15 PUGLIELLI etal. METHODS Dataset We collated multiple unpublished datasets (except for Govaert et al., 2023) including ITV- level observation from 21 research groups across the Italian peninsula. The datasets were provided as plotlevel individual observations for species whose cover summed up to 80% of the relative abundances in their sampling plots (Pakeman & Quested,2007). The species abundances recorded in the plots were also used to classify the community plots according to the first two levels of EUNIS habitat classification (Chytrý etal.,2020; Figure2; TableS1). Hereafter, habitat classification follows EUNIS level I Habitat Types nomenclature and it will be referred to as ‘habitat types’. From this initial collection, we retained all individual measurements of leaf area (LA, mm2) and leaf mass per unit of leaf area (LMA, g m−2) for vascular plants. LA and LMA were always determined using standardized protocols (Pérez- Harguindeguy etal.,2016). Species nomenclature across datasets was homogenized using the R package ‘WorldFlora’ (Kindt,2020, accessed in July 2023) that uses the World Flora Online (https:// www. world flora online. org/ ) as the backbone. From this initial compilation, we retained only seed plants (Spermatophyta), thus removing ferns because we did not want to bias the definition of the functional leaf unit (Vasheka etal.,2019). To identify possible outliers that can strongly affect the degree of ITV, we excluded the values of LA and LMA that exceeded the range reported in global compilations (Carmona et al., 2021; FIGURE 2 Dataset features. (a) Geographic location of the plots in the dataset within the Italian peninsula. Points are coloured according to the EUNIS habitat type classification at level I. (b) Frequency distribution of individual log10- transformed leaf area (log10 LA, mm2) measurements in the whole dataset. (c) Frequency distribution of individual log10- transformed LMA (log10 LMA, g m−2) measurements in the whole dataset. The shaded area in (b) and (c) represents the distribution of log10 LA and log10 LMA in the global dataset of Carmona etal.(2021). The percentage of overlap (Ov.) between trait ranges in our dataset and those in Carmona etal.(2021) is also reported. (d–h) Distribution of observations within the trait space defined by log10- transformed LA and LMA across habitat types. The trait space across habitats was defined using multivariate kernel density. Contour lines show the 0.99 and 0.50 quantiles of the multivariate probability density. Colours follow the same coding as in panel (a). The number of species and observations in the whole dataset (b, c) and per habitat type (d–h) are also shown. 14610248, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ele.14396 by Universidad De Sevilla, Wiley Online Library on [31/05/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 6 of 15 | INTRASPECIFIC VARIATION OF LEAF STRATEGIES Díaz etal.,2016). We decided to keep only those species with at least 4 individual trait measurements to reach a compromise between the number of species retained and the minimum sample size recommended for the selected traits in standardized protocols (n = 5, Pérez- Harguindeguy etal.,2016). Finally, for each species, we kept in the final dataset only unique LA or LMA measurements, that is, LA or LMA values (rounded at the third decimal) that were never duplicated within a single species. This step was needed since some datasets provided trait values at the individual level only for one plot and the same trait value was then kept constant across sampled plots. The final dataset that was used in subsequent analyses included 3153 unique individual measurements of LA and LMA for 167 species (49 families,132 genera) in five habitat types: coastal dunes (16 species), forests (49), grasslands (81), heathlandsscrubs (34) and inland surface waters (15). The EUNIS Level II habitats, together with their proportion of trait observations within each EUNIS Level I category, are listed in TableS1. The geographic coverage of the final dataset and the trait distribution across and within habitat types are shown in Figure2. Data analysis To answer Q1 (Figure1a), we first defined the trait space at the BTV- level (i.e., not accounting for ITV) for all data pooled and for each habitat type using medians of log10- transformed trait values. Since our focus is on variancerelated effects due to the inclusion of ITV, we decided to use medians instead of means because the species had a different number of observations, possibly resulting in a dependence between trait means and variances. Subsequently, we analytically identified the eigenvectors and eigenvalues from the LALMA covariance matrix (Figure1a). This approach uses the principles of major axis regression (Warton etal.,2006) but with a focus on both the eigenvectors and with no attempt to estimate a slope for the scaling between the two traits. This approach proved useful since LA and LMA are largely independent traits (Pearson's r = −0.01, p = 0.47), a condition which is necessary to apply the proposed analytical framework. To account for the overall effect of ITV on the trait space, we built 199 trait spaces per habitat by including at each iteration 1 individual of each species randomly sampled from all the individuals available for that species (Figure1a). Resampling one individual per species at each iteration guarantees that all species have the same weight in defining the trait space at the ITV level, and that the latter is defined on the same number of observations as the BTV trait space. The threshold of 199 repetitions was chosen to reduce the effect of an uneven number of individual trait values per species (mean = 18.8, range = 4–406). That is, we considered only a random subsample of individuals from species with many observations (>199) to balance the possible truncated ITV trait distribution for species with fewer observations (≥4 and ≤199). We, therefore, assessed how ITV modified the major variance axes in the trait space by measuring the rotation of the eigenvectors when including ITV data compared to the BTV- level eigenvectors (Figure1a). The angle between eigenvectors was calculated using the R function provided by Bueno etal.(2023). The significance of the effect of ITV was quantified by testing if the generated angle distribution (n = 199) was greater than zero using the as.randtest function in the ‘ade4’ R package (Dray & Dufour, 2007). Finally, we performed a Procrustes analysis between BTV and each generated ITV trait space across all data and per habitat type to test the correlation between the two spaces. The greater the observed correlation between trait spaces, the smaller the effect of ITV in modifying the position of the species within the trait space. Procrustes analysis was performed using the protest function in the R package ‘vegan’ (Oksanen etal.,2022). To answer Q2 (Figure1b), in each of the previously described iterations, we quantified the ratio between ITV- and BTV- level eigenvalues associated with each eigenvector defining the trait space (Figure1b). In this way, we could quantify whether the ratio between the variance that is captured along an eigenvector at the ITV- vs. BTV- level was significantly greater than a ratio of 1 (using the as.randtest function, as specified above), indicating that accounting for ITV increases the variance captured along a given trait space dimension. Since we hypothesised that the relative effect of ITV would be greater on the trait dimension with the lowest phylogenetic signal, we first retrieved phylogenetic information for the 167 species in the dataset using the phylo.maker function included in the ‘V.phylomaker’ R package (Jin & Qian,2019). The tree was built using the ‘Scenario 3’ option that bounds the missing species at 1/2 point of the family branch with some exceptions (Jin & Qian,2019; Qian & Jin,2016). Polytomies in the phylogenetic tree were randomly resolved generating a set of 100 trees through the ‘ape’ R package (Paradis & Schliep,2019). Then, we calculated the phylogenetic signal associated with the median LA and LMA (dimension 1 and 2 of the trait space, respectively) by computing Pagel's λ (Pagel,1999) for each tree using the ‘phylosignal’ R package (Keck etal.,2016). To answer Q3 (Figure1c), we compared a set of functional diversity metrics, obtained using kernel densitybased methods that rely on trait probability density functions (Carmona et al., 2016), between the BTV- and ITV- trait spaces both for the whole dataset and per each habitat type. Specifically, we assessed functional richness (in our case, the area occupied by the trait space), functional divergence (i.e., the average distance of observations from the centre of gravity of a trait space) and overlapbased dissimilarities calculated 14610248, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ele.14396 by Universidad De Sevilla, Wiley Online Library on [31/05/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 7 of 15 PUGLIELLI etal. as 1—overlap area between trait spaces (Figure1c; Carmona etal.,2016; Mason etal.,2005). The overlapbased dissimilarity allowed us to quantify the differences in the density of trait values between BTV and ITV trait spaces, as well as their degree of nestedness (i.e., the proportion of shared trait space). These calculations were performed at each resampling. However, the resampling procedure inherently modifies the covariance matrix of the sample and, consequently, the trait space rotation. Trait space rotation, while useful to understand the effect of ITV in altering the leading variance directions in the trait space (Figure1a), can alter the interpretation of the overlapbased dissimilarity due to the changing covariance between samples that modifies how the data points occupy the trait space (i.e., their multivariate density). To avoid this artefact, we used previously calculated eigenvectors to reproject LA and LMA data into a trait space with null covariance (i.e., fully orthogonal rotated axes) between the two traits (Figure1c). This approach allowed us to highlight and interpret the differences in functional diversity metrics only in terms of changes determined by the inclusion of ITV. All the analyses were carried out at the 0.99 and the 0.50 quantile of the probability distribution of species in the trait space to better characterize differences when targeting the whole target dataset (0.99 quantile) or the most likely trait combinations (0.50 quantile). All the kerneldensity- based analyses and visualization were implemented by combining the functionalities of the ‘TPD’ (Carmona etal.,2019) and the ‘funspace’ (Carmona etal., 2024) R packages. The differences between the BTV and ITV level metrics were tested using the as.randtest function, as specified above. RESULTS Q1: Does ITV modify the direction of the leading variance components of a trait space? The median rotation angle of the trait space when including ITV was 0.77° for all the data and ranged between 1.05° and 3.80° depending on the habitat (Figure3a–f). In all cases, the distribution of angles FIGURE 3 Trait space rotation by the inclusion of ITV. Trait space rotation refers to the median angle (°) between the eigenvectors calculated from the trait'scovariance matrix at the interspecific (BTV, grey eigenvectors) and at the intraspecific level (ITV, cyan eigenvectors) across 199 resampling for the whole dataset and per habitat type (see Methods). Note that in some cases the cyan lines have an opposite direction compared to the BTV- level eigenvectors. However, it must be noted that eigenvectors are not directional and extend to both directions. Here we just provided a simplified visualization to increase clarity. The 5th and the 95th quantiles of the angle distribution are shown in parenthesis. Number of species and observations per habitat is shown in Figure2 and TableS1. Leaf area (LA, mm2), that is trait dimension 1 (Dim. 1) and the leaf mass per unit area (LMA, g m−2, Dim. 2) were log10- transformed before the analyses. A pvalue <0.05 indicates that the observed angle was statistically greater than zero. 14610248, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ele.14396 by Universidad De Sevilla, Wiley Online Library on [31/05/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 8 of 15 | INTRASPECIFIC VARIATION OF LEAF STRATEGIES between the BTV and ITV eigenvectors across 199 repetitions was statistically greater than zero. These results were robust to differences in the number of species among habitats; a different species data aggregation scenario that simulates data downloading from an online database; and the varying number of individuals in the angle determination (NoteS1 for details, FiguresS1–S5, TableS2). However, the rotation partly depended on the growth form (woody vs. herbaceous, FigureS3). The Procrustes analysis revealed that ITV and BTV trait spaces were always highly and significantly correlated (p = 0.001) when analysing all the data (0.93 ± 0.01) and per each habitat type (0.95 ± 0.01 for heathlandsscrubs, 0.92 ± 0.04 for inland surface waters, 0.92 ± 0.02 for grasslands, 0.90 ± 0.03 for coastal dunes and 0.90 ± 0.02 for forests). Q2: Does ITV increase the variance explained along each trait space dimension? Accounting for ITV increased the amount of variance captured by each trait dimension when analysing the whole dataset, but the magnitude of this effect depended on the considered habitat type (Figure4a,b). As a general pattern across habitats, accounting for ITV mostly increased the variance captured along the second trait dimension (on average 1.62- fold increase, Figure4b), reflecting changes in LMA (see Figure3a–e), compared to the first trait dimension (on average 1.01- fold increase, Figure4a), mostly reflecting changes in LA. However, the observed increase in the ratios was never significant for LA in each habitat type (Figure4a,b). Furthermore, LA showed a strong and significant phylogenetic signal (Pagel's λ = 0.70, p < 0.001) compared to LMA (Pagel's λ = 0.50, p = 0.56). Q3: Does ITV modify the area and the functional structure of the trait space? We found that accounting for ITV overall increased the functional richness of the trait space at each considered quantile compared to the BTV- level trait space (Table1). However, among habitats and at the 0.99 quantile of the probability density function, the functional richness of the trait space that accounted for ITV was significantly greater than that of the trait space defined at the BTV- level only for forests and grasslands (Table1, p < 0.05) and not significant for the remaining habitat types. The same pattern was observed when differences were evaluated at the 0.50 quantile of the probability density function, but the differences were statistically significant only for coastal dunes, forests and grasslands. Accounting for ITV had a small effect on the functional divergence of the trait space compared to the trait space defined at the BTV- level. The only exception was detected for coastal dunes that displayed greater functional divergence at all quantiles in the ITV- level trait space compared to that at the BTV- level. The overlapbased dissimilarity between ITV- level and BTV- level trait space was generally low both for all data (0.20 ± 0.03 and 0.13 ± 0.02 at 0.50 and 0.99 quantiles, respectively, median ± SD) and across habitats irrespective of the type of data aggregation (range = 0.12–0.36 across quantiles, Figure5a,b, TableS2), but we overall observed a higher dissimilarity and a lower nestedness at 0.50 compared to the 0.99 quantile (Figure5a,b, TableS2). DISCUSSION We identified three ways in which ITV affects the trait space defined at the BTV level. In particular, ITV: (Q1) rotates the leading dimensions of a trait space; (Q2) FIGURE 4 The effect of ITV on the variance captured by each trait space dimension. Distribution of the ratio between eigenvalues when including or not including ITV in building the trait space (eigenvalue ITV/BTV, n = 199) for the whole dataset and per each habitat type for (a) dimension 1 (corresponding to leaf area, LA) and (b) dimension 2 (corresponding to leaf mass per unit area, LMA, see Figure3). The dots show the values obtained across resampling (see Data analysis). The horizontal dashed line indicates a ratio of 1. **Indicates a pvalue <0.05 for the test of the observed ratio distribution against a value of 1. 14610248, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ele.14396 by Universidad De Sevilla, Wiley Online Library on [31/05/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 9 of 15 PUGLIELLI etal. TABLE 1 Metrics of functional diversity of trait space at the interspecific (BTV) and intraspecific (ITV) level by habitat type and for all data pooled. Habitat type Quantile Functional richness Functional divergence BTV ITV BTV ITV Coastal dunes 0.99 2.12 2.49 ± 0.47 0.51 0.60 ± 0.03 0.50 0.42 0.64 ± 0.13 0.67 0.71 ± 0.02 Forests 0.99 2.47 2.74 ± 0.28 0.56 0.56 ± 0.01 0.50 0.49 0.57 ± 0.06 0.69 0.69 ± 0.02 Grasslands 0.99 3.42 4.51 ± 0.56 0.54 0.49 ± 0.03 0.50 0.45 0.68 ± 0.07 0.69 0.68 ± 0.01 Heathlandsscrubs 0.99 3.84 4.12 ± 0.46 0.53 0.54 ± 0.03 0.50 0.79 0.86 ± 0.09 0.72 0.71 ± 0.01 Inland surface waters 0.99 2.66 2.62 ± 0.25 0.59 0.64 ± 0.02 0.50 0.70 0.78 ± 0.08 0.70 0.71 ± 0.02 All data 0.99 4.61 5.65 ± 0.38 0.51 0.48 ± 0.02 0.50 0.68 0.86 ± 0.05 0.69 0.70 ± 0.01 Note: Differences in functional richness and functional divergence between a trait space defined at the betweenspecies level (BTV) and one including withinspecies trait variability (ITV) at the 0.99 and at the 0.50 quantile of the multivariate probability density function per each considered habitat type. The values are median ± standard deviations calculated across 199 resamplings (see Data Analysis section). Bold values indicate significant differences (p ≤ 0.05) between the distribution of functional richness and divergence obtained at the ITV level (see the Data analysis section) and its BTV counterpart. FIGURE 5 Overlapbased dissimilarity and nestedness between ITV vs. BTV trait space by habitat type and considering all data pooled. (a) Dissimilarity and (b) nestedness between ITV and BTV trait spaces at the 0.99 and at the 0.50 quantile of the probability density functions (n = 199). 14610248, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ele.14396 by Universidad De Sevilla, Wiley Online Library on [31/05/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License