High intraspecific trait variation results in a resource allocation spectrum of a subtropical pine across an elevational gradient

Plant functional traits are broadly used to quantify and predict impacts of climate change on vegetation. However, high intraspecific trait variation can bias mean values when few measurements are available. Here, we determine the extent of individual leaf trait variation and covariation across a highly heterogeneous environmental gradient for a widely distributed subtropical pine. We demonstrate the implications of trait variation for characterising species by assessing data availability and variability across the Pinus genus.

Yet, intraspecific trait data availability is limited, with even some of the most data rich species having only moderate coverage in the TRY global plant trait database and uneven data availability between traits, species, communities and regions (Grime, 2006;Niinemets, 2014). Uncritical use of trait data in multi-species analyses risks comparing species with small or widely varying sample sizes that overlook inherent intraspecific variation, with particular care needed when analysing trait data compiled from different sites and time periods (Niinemets, 2014).
Particularly large intraspecific variation is expected in heterogeneous environments, along environmental gradients and in species with large distribution ranges Anderegg et al., 2021;Bussotti et al., 2015;Körner, 2007;Rosas et al., 2019). With marked and rapid declines in temperature with increasing distance above sea level, elevational gradients provide an ideal opportunity to test plant responses to different environmental conditions (Malhi et al., 2010;McGill et al., 2006). Across elevational gradients, after controlling for physiographic factors such as slope and aspect, plant functional traits are expected to vary in response to factors such as declining air temperature and atmospheric pressure, and increasing solar radiation with increasing elevation (Körner, 2007;Körner et al., 1986;Sundqvist et al., 2011). Plants at higher altitudes are typically smaller, with smaller and thicker leaves (Bresson et al., 2011;Körner, 2007) and fewer stomata than their low elevation counterparts, reflecting restricted water availability (Schoettle & Rochelle, 2000) and low temperatures at high elevations (Körner, 2007).
While there are some common trends in plant trait variation across elevational gradients, information is limited on how consistent these fundamental patterns are, restricting our ability to make predictions of how vegetation might change under different environmental conditions. Elevational patterns can be complex and nonlinear, with traits influenced by several interacting drivers such as competition, ecosystem productivity, soil moisture and fertility, clear-sky turbidity, hours of sunshine, wind, season length, geology and human land use (Körner, 2007;Sundqvist et al., 2013). Mountains are highly heterogeneous, with environmental controls differing across local, regional and continental scales (Jobbágy et al., 1996;Midolo et al., 2019;Morley et al., 2018;Sundqvist et al., 2013). Fundamental elevational trends may not hold across latitudes, with information on trait variation in the tropics limited compared to temperate northern hemisphere regions (Wilson et al., 1999, Jetz et al., 2016, but see Chaturvedi et al., 2011Poorter et al., 2008;Wright et al., 2010). Rapid declines in temperature in the tropics by around 5.2-6.5°C every 1000 m a.s.l. with increasing elevation (Colwell et al., 2008), along with high species diversity and ecological trade-offs (Boucher et al., 2013;Wright, 2002), may decouple trait variation from expected trends.
Leaf traits are commonly used to describe physiological differences and ecological strategies among plants (Grime, 2006;Poorter & Bongers, 2006;Westoby, 1998;Westoby et al., 2002;Zanzottera et al., 2020), as they are fundamental for gas and water exchange, carbon assimilation and photosynthesis (Dong et al., 2020;Donnelly et al., 2016;Smith et al., 1997). Leaves can be quick and easy to measure, which is an advantage in remote and challenging areas to access such as tropical mountains (Garnier et al., 2001;Weiher et al., 1999;Westoby, 1998). Leaves can vary in size and in anatomical features K E Y W O R D S environmental gradients, intraspecific variability, leaf economics, Pinaceae, plant functional traits such as stomatal density along environmental gradients, reflecting differences due to genetics and phenotypic plasticity (Donnelly et al., 2016). Plant traits often strongly covary, with suites of traits reflecting different resource-use strategies and providing greater insight on plant strategies than considering traits in isolation. The leaf economic spectrum (LES) is broadly used to describe plant ecological approaches spanning from fast returns on investments to slow growth and persistence (Fajardo & Siefert, 2018;Messier et al., 2017;Niinemets, 2014;Pan et al., 2020;Wright et al., 2004Wright et al., , 2005. The LES is often quantified through the metrics of specific leaf area (SLA) and leaf dry matter content (LDMC), with low SLA and high LDMC associated with persistence and high SLA and low LDMC associated with short-lived and productive strategies (Guo et al., 2018;Poorter & Bongers, 2006;Roche et al., 2004;Wilson et al., 1999). However, large differences between individuals may result in intraspecific differentiation in resource allocation strategies, and hence position along the LES, which may complicate comparison of trait-environment relationships across species (Fajardo & Siefert, 2016Messier et al., 2018;Niinemets, 2014). Determining the extent to which intraspecific variation follows broad LES patterns and identifying whether traits covary predictably will provide a better understanding of the links between vegetation and environment, enabling more accurate predictions to be made of the impacts of environmental change.
Given the need to better understand intraspecific trait variability and its potential drivers across different environmental conditions, we sought to determine the extent to which fundamental patterns in leaf trait variation are consistent across space for a very widely distributed subtropical pine species. By focussing on needles, as fundamental components of conifer physiology, we aimed to capture variation in size dimensions and other characteristics associated with phenotypic plasticity and genetic adaptations (Donnelly et al., 2016).
To better understand implications of trait covariation, we explored the extent to which the LES holds within this species and its predictability with elevation. The Pinus genus is diverse and globally distributed (Ioannou et al., 2014). Therefore, to demonstrate the implications of trait variation for characterisation of the species and the wider genus, we determined how well represented individual Pinus species are by trait data and compared variability within traits between congeneric species. Specifically, we sought to determine: (1) To what extent traits vary predictably across the range of Pinus taiwanensis, (2) whether the measured traits covary, or if lack of covariation leads to changing relationships between traits across the range and, (3) if variability within traits restricts our ability to use species-level trait values to represent a species or differentiate between others within the genus.

| Study system and species
Pinus taiwanensis Hayata (Pinaceae) is widespread across the island of Taiwan, where it is generally considered endemic, although closely related species are found in mainland China and Japan (Fu et al., 1999).
The species spans diverse climatic and habitat conditions, extending over an exceptionally large range in the Central Mountain Range, spanning from around 500 to 3200 m a.s.l. (Taiwan National Forest Inventory, 2013). Taiwan is a subtropical island, experiencing warm and humid conditions, which transition through temperate to alpine with increasing elevation (Li et al., 2013). At low elevations (below 500 m a.s.l.), urban and agricultural land dominates following widespread deforestation, but natural forests are abundant above this height, transitioning from evergreen broadleaved forests into areas of mixed forest, deciduous broadleaved and evergreen broadleaved forest at higher elevations (Li et al., 2013). While P. taiwanensis is found scattered through mixed forest and open habitat at lower elevations, it typically increases in abundance and occurs in monodominant stands at high elevations (Li et al., 2013). As an early successional species, it primarily favours light and humid conditions (Cai & Liu, 2017), regenerating quickly on disturbed land (Chou et al., 2009).

| Site selection and sampling
Mature adult trees >10 cm diameter at breast height (dbh) were sampled between 7 and 30 November 2019 across an elevational gradient in the northern part of Taiwan's Central Mountain Range ( Figure 1, Table S1). Sample sites were identified from local knowledge of the species distribution and selected to achieve even sampling across the elevational gradient. Sites ranged from 495 m to 3106 m a.s.l. and included monodominant stands, mixed forest and scattered individual trees in open habitat, due to changes in forest composition and increasing abundance of P. taiwanensis with elevation preventing sampling from forests with consistent structure.
Coordinates were recorded at each sampling location and elevation, aspect and terrain slope extracted from a high-resolution Digital Elevation Model (NASA JPL shuttle mission, 30 m pixel size) for each sampling location. Several branches were collected per tree, selecting the highest accessible branches (largely 10-30 m above the ground) with fully expanded, mature, sun leaves (Cornelissen et al., 2003). Branches were transported to the laboratory in sealed poly bags containing a small piece of moist filter paper (Kitajima & Poorter, 2010). A total of 92 trees were sampled across 15 sites, with minimum 2 and maximum 12 trees per site.

| Needle trait measurements
Five leaves were randomly selected per tree for functional trait measurements. Eight functional traits were measured: needle thickness (mm), area (mm 2 ), length (mm), fresh and dry mass (mg), specific leaf area (SLA, mm 2 mg −1 ), leaf dry matter content (LDMC, mg g −1 ) and stomatal row density (SD, per mm 2 ). These traits were chosen as they have previously been associated with high intraspecific trait variation (Siefert et al., 2015), ecological changes along elevational gradients (Bresson et al., 2011;Körner, 2007;Schoettle & Rochelle, 2000), sensitivity to environmental variation (Donnelly et al., 2016) and LES strategies of growth and survival (Wright et al., 2004). Needle thickness was measured at the centre of the needle using digital callipers (0.1 mm precision), length using a ruler and area using the Apple iOS iPad application Leafscan (Anderson & Rosas-Anderson, 2017).
Abaxial and adaxial longitudinal stomatal rows were counted using a microscope (Donnelly et al., 2016). Fresh mass was obtained from a further 20 randomly selected needles weighed at 0.1 g precision, as the mass of individual needles was too small to be detected on the balance available to us in the field. Needles were stored in individual tea bags to air dry, further dried in an oven at 60°C for several days then weighed to 0.001 g precision. Individual needle fresh mass was estimated by calculating the difference in fresh and dry mass of the bulk sample, dividing by N = 20 and adding to individual needle dry needle mass measurements.
SLA, the ratio of leaf surface area to dry mass, was calculated: where r = thickness, l = length and m = dry mass, to account for the curved shape of pine needles (Donnelly et al., 2016). LDMC, the ratio of leaf dry to fresh mass, was expressed as dry mass (mg) proportional to fresh mass (g) (Pérez-Harguindeguy et al., 2013;Vaieretti et al., 2007;Wilson et al., 1999). SD per leaf was calculated as:

| Additional trait data
To assess intraspecific needle trait data availability more broadly, the number of data records listed in the TRY global plant database were obtained for the key needle traits of thickness, length, LDMC and SLA across the Pinus genus.

| Analysis
Trait variation was initially explored using descriptive statistics (minimum, maximum, mean, standard deviation, coefficient of variation (CV) and 95% confidence intervals [CIs]). To further explore the variation between needle functional traits and topographic variables of elevation, slope and aspect, we fitted linear mixed effects models (LMMs). All traits met test assumptions according to diagnostic plots on model residuals. Global models included fixed effects of elevation, aspect (as four cardinal directions; north, east, south and west) and terrain slope and random effects of individual trees. Elevation and slope were scaled and aspect was modelled as categorical as it is a circular variable. We incorporated pairwise interactions based on our expectation of interactive effects of topography on trait values. Variance inflation factors of the linear predictors were checked for multi-collinearity, with values between 1 and 2 suggesting no multi-collinearity (Zuur et al., 2010). Guided by our expectation of changing trait values with elevation, all models included elevation, but differed in their inclusion of slope and aspect and their interaction with elevation. Models were ranked by the Akaike information criterion (AICc) from lowest to highest. The most parsimonious models were selected, within two AIC units of the lowest AIC model, and fitted via maximum likelihood (ML). The percentage contribution of fixed and random effects (Pseudo-R 2 conditional ) in explaining functional trait variation was calculated (Nakagawa & Schielzeth, 2013).
A paired samples Wilcoxon test was used to compare abaxial and adaxial SD.
To assess trait covariation in multidimensional space across the full elevation range, we performed a principal component analysis (PCA). The first two components with eigenvalues >1 were retained for further inspection and analysis (Zwick & Velicer, 1986). LMMs were run using component axes as the response and elevation, slope and aspect as predictors, with tree as a random factor. Both component axes met test assumptions according to diagnostic plots on model residuals.
To assess intraspecific needle trait data availability in TRY, we extracted needle length, thickness, LDMC and SLA for individual Pinus species. We quantified how many species had publicly available records for each trait and the mean number of measurements per trait. We calculated how many species had ≥20 measurements per trait and the proportion of species with records with ≥20 measurements, as this measurement intensity has been reported as the minimum acceptable threshold to capture intraspecific variation (Kattge et al., 2011;Kattge et al., 2020). We further investigated a subset of the Pinus species with ≥20 records, selecting the nine most  using the packages 'lme4' (Bates et al., 2015), 'lmerTest' (Kuznetsova et al., 2017) and 'MuMIn' (Barton, 2020).

| RE SULTS
There was substantial intraspecific variation in all traits in P. taiwanensis, with coefficients of variation ranging from 20% to 44% ( Figure 2, Table S2). Needle dry mass was particularly variable across the gradient, with the highest value 18 times larger than the smallest value. SLA increased by a factor of 10 over the gradient, fresh mass and SD increased by a factor of nine and thickness increased by a factor of eight. Needle area and LDMC increased fivefold and length fourfold over the gradient.

| Topographic variation in functional traits
In general, needles became smaller, had higher SD and lower LDMC with increasing elevation, although there was considerable variation around this pattern (Figure 2). Fixed and random (tree) effects (Pseudo-R 2 conditional ) accounted for 47%-95% of the model variance, comprising 0.07%-59% of variance from fixed effects and 24%-57% from random effects (Table 1). There were significant elevation main effects for all traits. Needle size was generally negatively associated with elevation, with significant elevation effects for length (t = −23.89, p < 0.001), area (t = −11.14, p < 0.001), fresh mass (t = −5.00, p = 0.014), dry mass (t = −2.03, p = 0.027) and thickness (t = 0.03, p = 0.006). Models of needle length, area, fresh mass and dry mass were improved with the inclusion of aspect as a fixed effect, with smaller leaves on east facing slopes ( Figure S1, Table S3). Slope gradient, and the interaction between slope and aspect, further explained differences in fresh and dry mass, with marginally lighter needles on slopes with moderate gradients (15-25°) ( Figure S2, Table S3). Overall, fixed (topographic) effects captured substantial variation in needle size, explaining over half of the variance in needle length (59%), over a quarter of the variance in fresh mass (38%) and dry mass (29%), and 22% of the leaf area variance. Variation in needle thickness was less strongly linked to topography, with fixed effects of elevation, aspect, slope and the interaction between aspect and slope explaining 15% of variance.
SD was positively associated with increasing elevation, with a significant elevation effect (t = 0.01, p < 0.001). The SD model was improved by including aspect, with higher SD on east facing slopes ( Figure S1). Fixed effects of elevation and aspect explained 38% of the observed SD variation. SD was significantly higher on the abaxial than the adaxial side of the needle, with a mean of 0.033 and 0.025 per mm 2 respectively (V = 79,307, p < 0.001). Resource-use strategies varied over the elevational gradient, with a significant decline in LDMC (t = −21.41, p = 0.020) with increasing elevation. Including slope and aspect, and their interaction, improved model fit, explaining 8% of the variance in LDMC. However, SLA did not significantly vary with elevation (t = 0.334, p = 0.739) with only 0.07% of the variance explained by elevation. LDMC was marginally lower on east facing slopes ( Figure S1).

| Covariation among functional traits
The first two component axes of the full PCA together explained 65% of the variance ( Figure S3). The first axis (PC1) showed a clear negative association with needle size (length, area and mass) and explained 46% of the total variance. The second axis (PC2) was negatively correlated with needle thickness and SD and accounted for 19% of the total variance. Needle area and length were tightly paired, with close positive associations also present for mass and LDMC ( Figures S3 and S4). PC1 significantly increased with elevation (t = 4.994, p < 0.001), with slope, elevation and aspect fixed effects explaining 42% of the variance and differences between trees attributed to 35% of the variance (Table 1, Figure 2). PC2 significantly declined with increasing elevation (t = −4.878, p < 0.001), with elevation explaining 16% of the variance and differences between trees explaining 51% (Table 1, Figure 2).
The TRY database held information on 123 Pinus species, with 78 species (63%) having at least one measurement for needle length, thickness, SLA and/or LDMC and 45 species (37%) having no data for any of these traits (Tables S4 and S5). The number of species with ≥20 trait measurements varied from 5 to 27 species, comprising 12.5%-53% of species with >0 measurements. The mean number of measurements per species was generally <20, with only needle thickness having a mean > 20. The extent of intraspecific trait variability differed considerably between species. For the subset of Pinus species we investigated further, the difference between the maximum and minimum species trait values ranged from 31 to 266 mm for needle length, 0.19 to 0.75 mm for needle thickness, 244 to 517 mg g −1 for LDMC and 8 to 52 mm 2 mg −1 for SLA (Figure 3).
Most species and traits had much smaller sample sizes than our P. taiwanensis data, suggesting variability is likely to be underestimated in many cases.

| DISCUSS ION
We observed large intraspecific needle trait variation in P. taiwanensis across its broad elevational range. As expected, directional changes linked to elevation were identified, aligning with ecological theory that species traits vary across elevational gradients (e.g. Campetella et al., 2019;Körner, 2007;Sundqvist et al., 2013) and plant function determines species distributions (Hulshof et al., 2013). A spectrum of resource allocation was present across the elevational gradient, suggesting resource-use trade-offs acting across the species range (Fajardo & Siefert, 2018). However, the strength of elevational trends varied between traits, with slope and aspect accounting for a portion of the variation, and often substantial unexplained variance. Our findings agree with suggestions that intraspecific differences contribute an important component of functional trait variability (Ahrens et al., 2021;Albert et al., 2010;Violle et al., 2012).
We provide evidence that pronounced intraspecific variation in resource allocation strategy can occur, likely associated with high environmental heterogeneity across a species range (Figure 2). Such F I G U R E 2 Pinus taiwanensis needle trait values per tree along an elevational gradient from samples taken in the Central Mountain Range of Taiwan. Blue lines represent predicted relationships based on linear mixed effect models, omitting any interaction terms for the purposes of graphical representation. Grey shaded areas show 95% confidence intervals. SLA, Specific leaf area; LDMC, leaf dry matter content; SD, stomatal density variability could potentially impact our ability to forecast responses of species to environmental change, since local-scale drivers might significantly modify LES-climate relationships within and between species. Furthermore, we illustrate the risk of not considering intraspecific variation when using published trait data by highlighting considerable trait variation within and between species across the Pinus genus in conjunction with overall low numbers of reported trait values and variability likely to often be underestimated (Figure 3).
Uncritically deriving and employing mean species trait values to predict community responses to ongoing climate change risks inadequately describing key functional differences between populations and individuals within and between species.

| Variation in functional traits across the species range
Elevation strongly influenced needle size, with trees from low elevations having larger needles (length, area and mass) than those at higher elevations. Needle size is typically influenced by the most limiting environmental factor (Schoettle & Rochelle, 2000), with small leaves expected for plants experiencing extreme drought (Meier & Leuschner, 2008), high or low temperature, nutrient shortages or high radiation (Pérez-Harguindeguy et al., 2013). Our findings follow expectations of patterns across elevational gradients, with needle size likely to be largely temperature driven, since plant productivity is tightly coupled with air temperature and low elevation sites in our study region are characterised by warmer temperatures (Hatfield & Prueger, 2015;Körner et al., 1986;Wright et al., 2017). Large needles at lower elevation suggest water availability is sufficient to prevent desiccation despite high temperatures (Scoffoni et al., 2011), with mean annual rainfall in the lowlands as high as 3756 mm per year (Shiau & Huang, 2015). Trees at low elevations may be maximising resource capture when availability is high.
Unexpectedly, SD increased with elevation. Trees are expected to restrict water loss at high elevations (Schoettle & Rochelle, 2000), where humidity and rainfall are typically low (Körner, 2007;Leuschner, 2000). Declining SD with elevation has been observed in conifers, such as P. flexilis (Schoettle & Rochelle, 2000), P. contorta and Abies lasiocarpa in the dry Mid-Southern Rocky Mountains (Hultine & Marshall, 2000). With high humidity across the elevational gradient in Taiwan, low SD is unlikely to be a drought response per se (Beerling & Chaloner, 1993;Luomala et al., 2005). However, it may be linked to water economy, as stomata may close due to high vapour pressure deficit or temperatures rising above plant photosynthesis thresholds in the humid tropics (Doughty & Goulden, 2009;Duursma et al., 2014;Oren et al., 1999). Changes in CO 2 partial pressure may be influential, with SD increasing as CO 2 partial pressure declines towards higher elevation and plants maximise carbon gain (Mott, 2009;Pato & Obeso, 2012;Woodward & Bazzaz, 1988).
Additionally, a trade-off may exist between SD and stomata size, with small stomata able to respond more rapidly to changing conditions, but present in higher densities (Wang et al., 2014). Our findings TA B L E 1 Summary of the selected models for each trait measured on Pinus taiwanensis needles sampled from the Central Mountain Range of Taiwan. Selected models were the simplest model with the lowest AIC which still included elevation as a predictor. More complex models were chosen only if they reduced AICc by >2 from a simpler model to factors such as changes in temperature, light, humidity and CO 2 concentration with elevation (Gou et al., 2005;Tiwari et al., 2013;Zhou et al., 2012). Trade-offs between water and CO 2 availability across the species range are likely to vary between temperate and tropical systems and could be further explored through common garden and controlled environment experiments to identify the extent of environmental and genetic drivers.
Aspect differences further explained some of the observed needle trait variation for all traits except SLA, while slope influenced needle thickness and mass. Smaller needles and higher SD were observed on east facing slopes, while larger needles were common to moderately steep slopes (~15-25°). High variability in vegetation across Taiwan's mountain forests has been linked to slope and aspect differences, with east and south facing slopes experiencing the greatest gains in forest area between 1963 and 2016 (Morley et al., 2020), and moderately steep, east facing slopes experiencing the highest seedling recruitment (Greenwood et al., 2014(Greenwood et al., , 2015. These differences are likely driven by variation in factors such as microclimate (Lembrechts & Lenoir, 2019) and soil moisture across different topographies (Körner, 2007;Lambrecht & Dawson, 2007), resulting in variations in stand development, composition and interspecific competition which likely further influence trait variation (Vilà-Cabrera et al., 2015). However, uneven sampling across slope and aspect categories across the elevation range, due to local variation in the distribution of the species, must also be considered, with elevational effects potentially masking slope and aspect differences.
Despite significant associations of traits with elevation, slope and aspect, unexplained trait variation ranged from 41% to 92%.
Such variability is unsurprising given the high environmental variability across the elevation range in this highly topographically variable environment. Needle thickness was particularly poorly explained by elevation, slope and aspect (15% of variance). Although the positive correlation between LDMC and needle thickness may suggest an association between thickness and other measures of needle size, needle thickness is generally considered independent from other traits and more complex (Roche et al., 2004). However, it may be linked to drought, nutrient shortages, older leaves (Pérez-Harguindeguy et al., 2013;Roche et al., 2004) or herbivore defence (Hanley et al., 2007). Along with SD, needle thickness was a key contributor to the second component axis (PC2), suggesting a link between water availability and needle thickness, perhaps driven by factors such as soil and microclimate. Identifying the sources of additional key influential factors acting along with elevation related temperature changes to drive trait expression will provide more accurate insights on plant responses to environment.
The high variation among the resource-use traits of LDMC and SLA suggests a spectrum of resource allocation within the species.
LDMC and SLA were related, following the widely accepted expectation that these traits show opposite trends to one another, as observed in the global LES (Wilson et al., 1999). However, overall, LDMC and SLA were poorly explained by differences in topography (fixed effects variance 8% and 0.07% respectively), suggesting a substantial contribution from alternative (unaccounted for) factors. Elevation appeared influential in resource allocation, with significant declines with increasing elevation of LDMC and the PC2 axis of variation, and a significant increase in PC1 with elevation. Differentiation with temperature and water appear likely, due to the close association of PC1 and needle size (length, area and mass) and PC2 with needle thickness and SD. However, observed intraspecific trait patterns may be linked to local-scale biotic interactions varying with elevation, as P. taiwanensis occurs in mixed stands at low elevations and monoculture at higher elevations (Boucher et al., 2013;Cárdenas et al., 2014) and interspecific interactions are considered to be stronger towards lower elevations (Hargreaves et al., 2019). The decline in LDMC with elevation suggests a shift in resource-use strategy from productive at low F I G U R E 3 Minimum, median and maximum needle length, thickness, specific leaf area (SLA) and leaf dry matter content (LDMC) for Pinus species in the TRY global plant trait database. Labelled numbers represent the number of measurements available in TRY for each species and trait. The nine most data rich species are included for each trait except LDMC, where only five species had ≥20 measurements and thus met the threshold outlined in (Kattge et al., 2011;Kattge et al., 2020) elevations to persistence at high elevations.

| Implications of trait variation
Differences in function are relevant in the context of community assembly, population dynamics and ecosystem processes under ongoing environmental change, with populations of the same species capable of responding differently (Bolnick et al., 2011;Kichenin et al., 2013;Rosas et al., 2021;Siefert et al., 2015;Tito et al., 2021;Vilà-Cabrera et al., 2015). For example, intraspecific variation can drive differences in soil stability (Ali et al., 2017), rates of leaf litter decomposition (Lecerf & Chauvet, 2008) and radial growth rates (Laforest-Lapointe et al., 2014). While the extent of intraspecific variation will differ between species, traits and ecological contexts (Kattge et al., 2011;Vilà-Cabrera et al., 2015), our current understanding of these differences is limited. Across species ranges there are many factors that can make sites suit- Although collecting new intraspecific trait data will not always be possible Martin et al., 2017), the TRY global plant trait database provides an opportunity to rapidly quantify intraspecific differences. Yet, many species do not have sufficient data for accurate assessments of intraspecific trait variability, making uncritical comparisons of mean trait values for many species ineffective for addressing differences at the local scale, along heterogeneous environmental gradients and the regional scale for broadly distributed species. Across the Pinus genus, species frequently had insufficient data for intraspecific comparisons and where sufficient data were available, we demonstrate considerable trait variation within and between species. While 20 measurements have been reported as the minimum acceptable threshold to capture intraspecific variation (Kattge et al., 2011;Kattge et al., 2020), larger sample sizes may be required for highly variable species. For example, P. taiwanensis represented by mean leaf length in Huisun Forest would give a value of 178.7 mm, while a mean from the highest sample location on Hehuanshan would be 46% lower at 82.6 mm. Taking population means rather than species means (Albert et al., 2010) or including the standard deviation with mean estimates will improve the species averaging method (Messier et al., 2010). Capturing variability in key areas and across environmental gradients will allow more accurate predictions to be made of changes in vegetation and the ecological implications of key functional differences Martin et al., 2017;Messier et al., 2010;Siefert et al., 2015;Siefert & Ritchie, 2016;Vilà-Cabrera et al., 2015;Yang et al., 2018).
Here, we demonstrate the extent of trait variation across the full elevation range of a particularly widely distributed pine, Pinus taiwanensis. While research on trait variation frequently focusses on differences between species, we provide evidence of large intraspecific needle trait and resource-use strategy variation, with substantial differences between individuals linked to elevational temperature variation and interacting drivers likely further contributing to variation.
We illustrate the risk of underestimating variability, which could have substantial implications for species comparisons, particularly in highly heterogeneous environments. Functional variation is crucial to understand and account for in the context of using mean species traits to predict population resilience, distribution shifts, conservation and management. With intraspecific variation emerging as a fundamental component of functional trait differences, improved understanding of variation across wide environmental gradients will provide critical insight on changes in vegetation over coming decades.

ACK N OWLED G EM ENTS
We thank the forestry students at the Forest Management Laboratory at the National Pingtung University of Science and Technology for their invaluable assistance in the field and laboratory.
Thanks also to Erin Stoll at the University of Stirling for her support with fieldwork data collection, Dr Thiago Silva for suggestions on data collection methods and Dr Daniel Chapman for his comments on initial data analysis. Permits for data collection were provided by