Intraspecific trait variation and species turnover in successional tropical forests: assessing trait imputation for community-weighted means

Accounting for intraspecific trait variation (ITV) is crucial to plant ecology for vegetation modeling efforts. ITV can be substantial; however, it remains unclear how ITV influences community-weighted mean (CWM) trait estimates. We use leaf and root trait data from 423 trees of 72 species from 15 Angiosperm families in combination with community data from 164 25 × 25 m plots comprising 580 species to evaluate the contributions of ITV and compositional turnover to CWMs, comparing unlogged primary tropical forest to selectively logged and clear-cut secondary tropical forest. We also examine the effect of imputing missing trait values using phylogenetic generalized linear modeling (PhyloPars) on CWMs. For six of the seven traits, ITV negatively covaried with community compositional turnover to generate larger CWM differences between forest types than observed if ITV was not integrated. For example, plot average-weighted mean specific leaf area was 10.7 and 10.4 m2 kg−1 for primary and secondary forests, not accounting for ITV, but shifted to 9.8 and 11.1 m2 kg−1 after doing so. Our results from 72 species assemblages were largely consistent with results using phylogenetically imputed traits for the entire community. The contribution of ITV to CWMs ranged from 25 to 75%, with ITV, not species turnover, driving CWM trait variation among successional forest types. CWM trait estimates became more conservative with forest age, whereas ITV for many traits showed an opposing acquisitive shift (i.e., increasing in leaf area or root length) and because of negative covariation between ITV and species turnover, weighted mean trait differences between successional forest types increased.


Introduction
Plant functional traits broadly scale along a fast-slow axis related to plant life-history variation and demographic rates (Wright et al. 2010;Reich 2014;Díaz et al. 2016;Ruger et al. 2020). For example, the most broadly measured plant functional traits distill to two main axes of functional variation-one related to leaf economics (e.g., specific leaf area-SLA, leaf nitrogen) and another related to whole plant size (Díaz et al. 2016). Despite this, there is substantial functional redundancy among plant taxa (i.e., within communities, Pillar et al. 2013) and considerable Communicated by Simon Pierce. intraspecific trait variation (ITV, Siefert et al. 2015) that muddles trends, for example, weakening trait-demography relationships (Poorter et al. 2018a;Yang et al. 2018). Functional trait variation also has implications for understanding plant niche partitioning of environmental resources (Laughlin and Messier 2015;Kunstler et al. 2016;Laughlin et al. 2020). For example, at the community level, functional traits have popularly been used to infer community shifts in plant resource-acquisition strategies across environmental gradients (e.g., environmental filtering or other mechanisms of community assembly, Escudero and Valladares 2016).
In tropical forests, functional diversity accounts for differences in forest productivity more so than taxonomic diversity (van der Sande et al. 2018), richness (Rosenfield and Müller 2020), or forest dynamics in space and time, especially regarding forest succession (Lohbeck et al. 2012;Letcher et al. 2015;Muscarella et al. 2016a;Hogan et al. 2018). It is often more straightforward to infer plant community dynamics or plant-environment relationships from a trait perspective than through changes in individual species densities or community demographics (Poorter et al. 2018a(Poorter et al. , 2018b, which is especially true as community diversity increases (Baraloto et al. 2010). When diversity is high, using a trait perspective can uncover complex associations between the environment (e.g., soil or light habitat) and plant strategies (Kraft et al. 2008), which may be obscured by traditional taxonomy. Moreover, when a functional perspective couples phylogenetics, inference about environmental-or plant lineage-structured trait and compositional variation can be strengthened (Baraloto et al. 2012;Draper et al. 2021); however, most studies use species-level trait data, ignoring ITV. Hence, there is motivation to understand whether incorporating phylogenetically structured ITV strengthens our ability to infer plant community-environment relationships (Violle et al. 2012;Siefert et al. 2015). This is of particular importance at the community scale because functional traits are used to inform the plant functional types or benchmark aggregate ecosystem functions in forest demography or ecosystem models, which predict tropical forest responses to global change (Warren et al. 2015;Fisher et al. 2018;Powell et al. 2018;Sakschewski et al. 2016;Xu and Trugman 2021).
Firstly, ITV arises because of genetic variation among individuals, the environment, and the interaction between the two (Bolnick et al. 2011;Siefert et al. 2015). A portion of ITV is due to the environment; thus, properly accounting for ITV across quantified environmental gradients can reveal trait-environment relationships (Bruelheide et al. 2018;Craven et al. 2018;Hofhansl et al. 2021). Moreover, slight differences in functional trait values have been reported for tree ontogenetic stages (i.e., canopy-dwelling adult trees > 10 cm dbh vs. smaller trees < 10 cm dbh, Spasojevic et al. 2014). Thus, ITV has consequences for estimating communityweighted means (CWM); however, accounting for ITV is labor intensive, requiring extensive trait sampling of individuals. Yet conceptually, ITV is central to the ideological reconciliation of ecological and evolutionary frameworks. On the one hand, ecological frameworks seek to account for ITV in predicting community attributes or ecosystem properties. On the other hand, the evolutionary perspective aims to understand how trait selection has shaped communities, ecosystems, and their constituent biota (Bolnick et al. 2011;Violle et al. 2012;Turcotte and Levine 2016). Studies seeking to account for ITV in plant communities have traditionally looked at ITV between populations, communities, and ecosystems (e.g., Albert et al. 2010). Several studies have looked at ITV across successional or environmental gradients in tropical forests (Husholf and Swenson 2010, Messier et al. 2010Paine et al. 2011;Auger and Shipley 2013); yet only a few still have sought to assess how ITV influences CWMs (but see Muscarella et al. 2016a, b;Katabuchi et al. 2017;Derroire et al. 2018;Subedi et al. 2019).
Additionally, forest communities typically move from acquisitive to conservative functional composition as tree composition changes with succession (Christensen and Peet 1984;Garnier et al. 2004;Swenson et al. 2012;Lohbeck et al. 2013;Letcher et al. 2015;Muscarella et al. 2016a). Thus, CWM functional trait values shift toward traits associated with species conservativeness (or slow life-history strategy, e.g., lower SLA and higher wood density-WD), partly because of species turnover and interspecific functional differences (Chazdon et al. 2010). The contribution of ITV in a successional context is less understood; however, these studies indicate that contributions can be substantial. For example, ITV of key leaf traits (leaf area, SLA, leaf nitrogen content, and leaf phosphorus content) accounted for 25-75% of CWM leaf trait variation in a Chilean temperate rainforest (Fajardo and Siefert 2018). Root traits are interesting because they are abstruse and ecologically complex (Bardgett et al. 2014;Laliberté, 2017), rarely studied for ITV, and relate to ecosystem function via the flow of nutrients and water into trees (Freschet et al. 2020). Thus, understanding how CWM root traits vary with environment and forest age is potentially key to modeling the functioning of tropical forests in future (Bardgett et al. 2014;Warren et al. 2015;Sakschewski et al. 2016;McCormack et al. 2017).
The lack of studies aimed at quantifying the contribution of ITV to CWMs is likely due to the laborious undertaking of measuring individual-level functional traits and community composition across environmental or successional gradients (Paine et al. 2011). This is particularly difficult in tropical forests, where large trees and high diversity complicate trait sampling (Baraloto et al. 2010). Methodological approaches that have sought to remedy this issue by imputing missing trait values are a potential solution. Still, each imputation approach has its assumptions, and some are more useful than others (reviewed by Swenson 2014). One of the first proposed methods was multivariate imputation chained equations (MICE, van Buuren and Groothuis-Oudshoorn 2010), which assumed randomness in missing values and uses the existing relationships and uncertainties of nonmissing data to yield statistically consistent imputations; however, this method has since been shown to be flawed for several reasons, perhaps the most consequential of which is not considering phylogenetic relatedness (Penone et al. 2014). Two more recent methods include Bayesian hierarchal probabilistic matrix factorization (BHPMF, Schrodt et al. 2015), which uses a hierarchical taxonomic structure (genera-family-order, etc.) to infer phylogenetic structure when imputing missing trait values, and PhyloPars (Bruggeman et al. 2009;Goolsby et al. 2017), which uses phylogenetic generalized linear models (pGLM) to impute missing trait values using branch lengths. A recent study compared the three methods, concluding that PhyloPars outperforms MICE and BHPMF for trait imputation, producing the most accurate estimates of missing values (Johnson et al. 2021). However, despite the choice of the imputation method, estimates of missing data can still be inaccurate, even with only a small portion (5-10%) of missing data, and can be sensitive to taxonomic biases (Johnson et al. 2021).
Few studies have examined how CWM patterns in root traits change with forest succession and how they may be different from leaf traits. A study from a temperate forest that used BHMPF to impute missing root trait data (Caplan et al. 2019) found that CWMs shifted toward conservative root strategies (thicker diameters, lower specific root length, higher root tissue density, and lower root nitrogen content). A different study sampled a strong edaphic chronosequence in temperate rainforest in New Zealand and found that CWMs became more conservative with increasing forest age (Holdaway et al. 2011). The difference between these two studies highlights the need for more studies that quantify the strength and direction of soil environment-structured root trait variation, forest succession-driven variation, and other ITV sources on CWMs. Here, we use root and leaf trait data collected in a paired sampling design along a 6.6-km transect in tropical forest of Hainan, China to scale up trait measurements to the community level using a complementary network of small vegetation plots (Xu et al. 2015a, b).
The motivation of this study is twofold-(1) to apply the best available methodology to quantify the contribution of ITV to weighted mean trait values using an a priori systematic trait sampling that measured forest age-structured ITV for a phylogenetically diverse subsample of trees, while (2) validating if the best available methods for phylogenetic trait imputing drastically affect inference at the whole community level. Specifically, we address two research questions-(1) to what extent do species turnover (i.e., changes in community composition because of forest age and environmental filtering effects), ITV, and their covariation influence CWM functional trait estimates? (2) How sensitive are these patterns to phylogenetic trait imputation for missing taxa? Concerning question 1, we expected that species turnover with forest age would influence estimates of community-weighted traits to a greater degree than ITV, because interspecific trait differences are typically larger than ITV at the landscape scale (Auger andShipley 2013, Derroire et al. 2018). For question 2, the testable null hypothesis is no difference in the relative contributions of species turnover and ITV when using directly measured versus imputed functional trait data. Accepting the null hypothesis would validate trait imputation methods and potentially guide functional trait measurement strategies in diverse communities.

Study site: Jianfengling, Hainan Island, China
The Jianfengling forest reserve (JFL) in Hainan, China (18°23′-18°15′N, 108°36′-109°05′E, Fig. 1) is a 47,200-ha montane tropical forest with a history of logging that dates to 1957, when about two-thirds of the area was either clear cut or selectively logged (Zhou 1995;Xu et al. 2015b). All logging ceased in 1994 under a state-wide logging ban, followed by a national logging ban in 1998 (Wenhua 2004;Zhou 1995). The reserve encompasses several vegetation life zones, from tropical semi-deciduous monsoon forest at the lower elevations to mossy high-elevation elfin forest, with evergreen-monsoon forest dominated by Podocarpaceae intermixed throughout at elevations < 1000 m (Huang et al. 1995). The most common vegetation life zone is tropical montane rain forest, which occurs at elevations between 600 and 1100 m. It is characterized by a mix of palms and broadleaf evergreen trees that reach an average canopy height of 18 m (Jin et al. 2013).
JFL has a tropical monsoon climate with seasonal rainfall, where the wet season occurs between May and October (Zeng 1995). Cumulative annual rainfall over the last 30 years has averaged 2700 mm (Wu 1995). The soils are derived from porphyritic granite and are classified as lateritic and humic yellow soils (Wu 1995). They are infertile, leached tropical soils with an exchangeable base content of 30 mL per kg and some aluminum accumulation. Nutrient supplies in these soils are principally derived from the surface accumulation of organic matter, with rates of mineral and organic matter cycling in the soils being slower than in other tropical soils (e.g., latisols, Wu 1995).

Root and leaf functional trait data
Functional traits were collected from May to July 2017. Four hundred and twenty-three saplings and small, 1 3 non-canopy-dwelling trees (i.e., trees > 1 m tall, but less than < 10 cm diameter at 1.3 m, for ease of sampling) of 72 species from 15 selected Angiosperm families were measured for leaf and root functional traits. Trait collection was done systematically along a 6.6-km transect that spanned primary and secondary forest with logging history (Fig. 1, Hogan et al. 2020). This sampling design targeted families across the Angiosperm phylogeny representative of the functional variation in root morphologies across the angiosperm phylogeny and sampled many individuals within each species (average of 6) to assess successionally structured ITV (Hogan et al. 2020). Three fully sun-exposed leaves and five fine root systems containing the most distal 3-4 root orders (McCormack et al. 2015) per tree were collected, washed, scanned, dried, and weighed following standard trait measurement protocols. Morphological traits were measured from scanned image analysis using WinFolia 2007b and WinRhizo 2016 (Regent Instruments, Quebec, Canada). For each tree, measured trait values were averaged across the three leaves and five root samples.
These individual-level trait data can be condensed to a few axes of key morphological variation among roots and leaves. A principal component analysis showed that three leaf traits-leaf area, specific leaf area (SLA), and leaf thickness and four root traits-root system average diameter, specific root length (SRL), root tissue density (RTD), and root branching intensity (RBI) accounted for most of the above-and belowground trait variation among Fig. 1 The study site is the Jianfengling Forest Reserve (JFL) of Hainan Island, China. Hainan Island is a small continental island of China's southeastern coast, shown in the red box of the inset map. The 47,200-ha JFL boundary is shown with red boundary. The 6.6-km transect of where functional traits of saplings and small trees were sampled is shown in black. The first half of the transect is in secondary forest with a history of logging, and the second half of the transect is in unlogged primary forest. Small 25 × 25-m plots (Xu et al. 2015a, b) in unlogged, primary, and either clear-cut or logged secondary forest are shown as pink and green points, respectively. Logging and forest disturbance in the secondary plots have happened in the last 20-55 years individuals (Hogan et al. 2020). Therefore, we focus on those seven traits in these analyses.

The network of small plots
We integrate data from a network of small plots previously established in the JFL forest reserve (Xu et al. 2015a, b). Local soil environments between the transect and the network of small plots were systematically sampled and found to be similar (see Online Resource 1). Between August 2007 and June 2009, a network of 164 1/16th-ha (25 m × 25 m) plots were established in the JFL reserve (Xu et al. 2015b, Fig. 1). The plots span the local variability in tropical montane rain forest range in altitude that was captured along the 6.6-km transect. Roughly one-third (52 of 164) of the small plots are in old-growth forest, with no record or visible evidence of logging in the last 200 years (Xu et al. 2015b). The remainder of the plots (112) are in forest areas that were either selectively logged or clear cut, with human land use ceasing between 20 and 56 years ago (Xu et al. 2015b). These secondary forests vary in past disturbance intensity, forest age, and structure, but may all be classified as secondary tropical forests (Xu et al. 2015b). Selective logging practices in these forests resulted in the removal of 30-40% of the mature stems ≥ 40 cm in diameter of commercially valuable timber species (Xu et al. 2015b;Zhou 1995). In each of the 164 small plots, all free-standing stems ≥ 2.5 cm in diameter at 1.3 m height were mapped, measured, and identified to species (see Fig. 2).

Data analyses: Lepš's intraspecific variability effects and assemblage-weighted means
We use Lepš et al.'s (2011) method for partitioning weighted mean trait values into components of compositional turnover, ITV (due to forest successional status), and their covariance. In the interest of scaling up individual trait data to the scope of JFL using tree community data from the network of small plots, we used an assemblage (i.e., partial community)weighted mean approach, with basal area species weights. We call this an assemblage-weighted mean (AWM) instead of the classical community-weighted mean (CWM), because we could not measure traits for all species found in the plot communities. This created the motivation to evaluate trait imputation, despite that CWMs typically include data for all species. Species for which we sampled functional traits accounted for 12.4% of species (72 of 580) and 37% of stems (20,673 of 65,144) in the network of small plots (Fig. 2). Basal area in the plots ranges from 16.7 to 77.8 m 2 ha −1 , averaging 44.2 m 2 ha −1 . The basal area of species in those plots for which functional traits were sampled along the 6.6-km gradient was between < 1 and 42.3, averaging 18.4 m 2 ha −1 . Thus, the proportion of plot basal area for which there was trait coverage ranged from < 1 to 82.3%, averaging 43.6%.
For each of the seven traits of interest, using species mean trait values from all individuals (n = 423) along the transect (x i ), we calculated AWMs for each plot. We term these the fixed averages (after Lepš et al. 2011), thus  (Hogan et al. 2020). Fifty-two of those plots are in well-conserved primary forest, while 112 are in previously logged or clear-cut secondary forest. Relative stem densities of these 72 species sum to 44% of stems in primary forest plots and 33% of stems in secondary forest plots. Colors of species names correspond to plant family in order from phylogenetically oldest (i.e., more basal) to youngest (i.e., more derived) from left to right where p i is the basal area fraction of that species within the plot (i.e., the community weight) and S is the species richness. For the AWMs, species weights were calculated after removing all species for which trait data were not sampled; this was necessary as weights must sum to one. We then calculated AWM trait values for each plot using habitat-specific mean trait values for each species (x i_habitat ). For example, for plots in secondary forests, trait means were calculated from individuals sampled within secondary forest (n = 225). We term these the specific averages (again after Lepš et al. 2011), as such Thus, the fixed average is the weighed mean irrespective of ITV due to forest successional status, which the specific average accounts for. The intraspecific variability effects were calculated as the difference between fixed and specific averages. Two-way analyses of variance (ANOVAs) assuming homoscedasticity were conducted for each of the seven functional traits chosen. Then, using the sum of squares of these tests, the variance in trait values was decomposed into intraspecific and community (i.e., species turnover) portions and their covariation (i.e., interaction). This was done using the 'trait flex ANOVA' function (see supplemental material of Lepš et al. 2011) in R v.3.5.0 (R Core Team 2018).

Data analyses: trait imputation and comparison to assemblage-weighted means
To validate the results regarding intraspecific variability effects using AWMs, we employed an approach to impute trait data for unsampled taxa. As we have explained, Phylo-Pars (Bruggeman et al. 2009) appears to be one of the best tools available for imputing functional traits (Penone et al. 2014;Johnson et al. 2021). Additionally, due to the high degree of phylogenetic conservatism in certain root traits (e.g., root diameter), PhyloPars is particularly appropriate. It uses a pGLM that calculates phylogenetic covariance of taxa from a Brownian motion model of trait evolution. Parameter estimates from the phylogenetic and functional trait variance-covariance matrix are then used to impute missing data values.
A coarse phylogeny for all 580 taxa in the network of small plots was obtained from Phylomatic v3 (www. phylo diver sity. net/ phylo matic, Webb and Donoghue 2005), using the 'slik2015' base tree, which contains the most comprehensive phylogenetic skeleton for tropical trees to date (Slik et al. 2018). Wood density (WD) was used as an additional trait in the pGLM, because we wanted to have one trait with data for all taxa. WD is a highly phylogenetically conserved trait (Swenson and Enquist 2007;Kraft et al. 2010), which can be used in the pGLM to constrain modeled trait values for taxa and clades with no measured trait data. WD values were obtained from the global database of WD values where available Zanne et al. 2009), and the CTFS-ForestGEO WD database (http:// ctfs. si. edu/ Public/ Datas ets/ CTFSW oodDe nsity/) via the 'BIOMASS' package ('getWoodDensity' function, Réjou-Méchain et al. 2017). Where species WD values were unavailable, mean values from a higher taxonomic classification (genus or family) were used (see Online Resource 3).
WD data were combined with three separate incomplete trait matrices, one using fixed trait values (x i , i.e., species mean values for all individuals measured along the transect), and two habitat-specific matrices (x i_habitat ), with the habitats being secondary and primary forest, as delineated along the sampling transect (see Hogan et al. 2020). pGLMs were fit separately to each incomplete trait matrix and predicted to impute missing data using PhyloPars (http:// zeus. few. vu. nl/ progr ams/ phylo pars, Bruggeman et al. 2009) and the Phylomatic-generated phylogeny. For the model settings in PhyloPars, we allowed for correlated evolution of the different traits but did not allow for intraspecific variation in the models, because, in our case, the separate fitting of the pGLMs accounts for this. We then used the imputed trait matrices and basal area species weights (p i ) of species in the small plots to recalculate fixed and specific CWMs for each plot. We re-ran the trait flex ANOVAs, comparing the AWM analysis, using only the species for which we measured functional traits, to the CWM analysis, which used the PhyloPars-imputed trait data for all 580 species found in the plots.

Estimates of assemblage-weighted trait valuesfixed vs. specific means
In one-third of cases, fixed AWM values for leaf functional traits differed from specific AWMs (Table 1, Fig. 3 left panels). For example, plot average fixed AWMs for leaf area were 77.5 cm 2 and 78.2 cm 2 for secondary and primary forest, respectively (a 0.7 cm 2 difference). Those estimates shifted to 70.5 cm 2 and 81.5 cm 2 for specific AWMs (a difference of 11 cm 2 , Fig. 3A). Fixed AWM estimates for SLA were 10.7 and 10.4 m 2 kg −1 for small plots in primary and secondary forest, respectively. In contrast, specific AWM estimates for SLA showed a greater difference between forest types, being 9.3 and 10.8 m 2 kg −1 for small plots in primary and secondary forest, respectively, (Fig. 3B). The trend for leaf thickness was similar to the trend for leaf area and SLA. AWM differences in leaf thickness between forest successional types were not significantly different for  2 and 3). AWM analyses use small plot tree community data and functional trait data for 72 species (i.e., only those for which trait data were directly measured along the 6.6-km transect in JFL), whereas CWM analyses use small plot tree community data and functional trait data for 580 species (i.e., the entire community for which we applied pGLM data imputation). For the AWM analysis, fixed mean values use all functional trait data from individuals sampled along the transect (n = 423), whereas specific assemblage-weighted means match functional trait data from plant individuals to plots with corresponding forest-type classification (n = 198 for primary forest and n = 225 for secondary forest). For the CWM analysis, fixed values use pGLM-imputed functional trait data for a matrix containing species means irrespective of forest type (i.e., all trait data collected along the 6.6 km transect), whereas specific values use pGLM-imputed functional traits data corresponding to either secondary or primary forest portions of the transect. Intraspecific variability = fixed mean-specific mean df degrees of freedom, SS sum of squares, F F-statistic, p probability for F-statistic Probabilities are denoted as *p < .05, **p < .01, and ***p < .001. Note that in this case, the mean-squared error is equal to the sum of squares because there is 1 degree of freedom; therefore, mean-squared errors are not shown. (See Online Resource 2 for AWMs and Online Resource 3 for CWMs)  fixed averages but were for specific averages (Table 1), being about 200 µm greater in secondary forest than in primary forest (Fig. 3C). Similarly, for root functional traits, fixed AWM values were similar and statistically indistinguishable (in all cases) between primary and secondary forest plots. At the same time, specific AWMs showed statistically significant differences in three of eight cases (Fig. 4). For example, fixed and specific AWMs for average root system diameter (Fig. 4A) showed opposing trends, with fixed estimates being slightly (< 0.01 mm) narrower in primary than secondary forest, but with specific estimates being slightly thicker (0.03 mm) in primary than secondary forest. For both the fixed and specific AWMs for root diameter, differences among forest types were not significant (Table 1). Fixed AWMs for SRL, the belowground analog to SLA, were 4.85 m g −1 in primary forest plots and 5.00 m g −1 in secondary forest plots, while the specific CWM estimates for SRL were 3.91 m g −1 in primary forest plots and 5.11 m g −1 in secondary forest plots (Fig. 4B). The difference of specific AWMs for SRL between secondary and primary forest plots was 1.11 m g −1 instead of a difference of 0.15 m g −1 for fixed AWMs. These differences went from being statistically non-significant for the fixed AWM comparison to significant for the specific AWMs (Table 1). The results for RTD were similar to those for average root system diameter (Fig. 4C). Fixed AWMs for RTD decreased from primary to secondary forest, being about 0.02 g cm −3 less dense in secondary than primary forest plots. Specific AWMs for RTD showed the opposite trend, with AWMs for RTD being about 0.3 g cm −3 denser in secondary than in primary forest. Differences between fixed AWMs were not statistically different, but the difference between specific AWMs were statistically significant (Table 1). Lastly, AWM RBI showed no significant differences when specific means were compared across forest types, with values averaging about 1.97 root tips cm −1 (Fig. 4D, Table 1). However, fixed AWM RBI values were statistically different, averaging 1.73 tips cm −1 in primary forest plots and 2.13 tips cm −1 in secondary forest plots (Fig. 4D), a weighted mean difference of 0.4 tips cm −1 .
In summary, the only difference in fixed AWM values by forest successional type for any of the seven functional traits considered was for RBI (Table 1, Figs. 3 and 4), but statistical differences in habitat-specific AWM trait estimates occurred for all traits except for root diameter, which had a marginally significant difference (p = 0.08, Table 1). Even when differences between specific AWMs and forest types were not statistically significant, they were always greater for 164 1/16th-ha plots in Jianfengling, Hainan, China. Points are mean CWM values (± 95% confidence intervals) for plots by forest type (52 in primary and 112 in secondary). Fixed CWMs (in black) use trait values irrespective of forest type, while specific CWMs (in gray) use trait data that account for ITV between forest types. The left panel in each figure shows assemblage-weighted mean (AWM) values, which use only measured data from 72 species, while the right panel shows PhyloPars pGLM-imputed CWM functional trait values using the entire community (580 species) than differences in fixed AWMs. In some cases, they were also different in sign (i.e., for SLA, leaf thickness, root diameter, and SRL, Figs. 3 and 4).

Comparing assemblage-weighted means to community-weighted means-evaluating pGLM-imputed trait values
Trends in CWM means confirmed the results of the AWM analyses. Like the AWM analysis, the pattern of having larger differences in habitat-specific means than for fixed means between forest types was consistent across all leaf (Fig. 3, Online Resource 3: Fig. S3.1) and root (Fig. 4, Online Resource 3: Fig. S3.2) CWMs using pGLM-imputed trait matrices. When pGLM-imputed functional trait matrices for the entire community were used, the variability among weighted trait estimates for plots decreased (i.e., 95% confidence interval breadth contracted). This decrease in estimate variability led to more statistically significant differences for specific CWMs relative to specific AWMs (Table 1). For leaf traits (leaf area, SLA, and leaf thickness), no qualitative difference emerged from CWMs using pGLMimputed trait data for taxa, versus only measured trait data (Fig. 3, Table 1, Online Resource 3: Fig. S3.1). The main difference between the AWM and CWM values in Fig. 3 is that pGLM-imputed CWM values for SLA were translated about 1 m 2 kg −1 greater across the fixed and specific means for both forest types (Online Resource 3: Fig. S3.1). Additionally, the difference between fixed CWMs for RBI was not statistically significant, whereas it was for the fixed AWMs for RBI (Table 1). Also, the difference in specific CWMs for root diameter was statistically significant (Table 1), but the difference between AWMs was marginally significant (Table 1).
Comparing AMWs to CWMs for root traits, we see similar patterns (Online Resource 3: Fig. S3.2). Using AWMs resulted in no differences between fixed and specific estimates for root diameter (Fig. 4A). However, when the pGLM-imputed data were used, the specific estimates for root diameter increased slightly for primary forest and decreased for secondary forest, creating a statistically significant difference between forest types (Table 1). For the three other root traits (SRL, RTD, and RBI), differences were accentuated slightly (Fig. 4B, C, D). This discrepancy in forest-type-weighted trait estimates resulted more from ITV than from species turnover (Table 1, Fig. 5). In many cases, for example, for SLA, root diameter, SRL, and RTD, ITV and species turnover had negative covariance values (Fig. 5), meaning that beta diversity in species composition drove weighted trait values one way (in our case, usually toward the conservative end of the plant economics spectrum, but not always), while ITV with environment caused shifts in the opposite direction (in our case, root morphologies became more conservative, while leaf morphologies were more acquisitive from primary to secondary forest, Hogan et al. 2020).

Interspecific variability effects-the relative effect of ITV and turnover on weighted mean trait estimates
Results from the interspecific variability effects from AWMs and CWMs (Table 1) can be used to infer the relative contribution of ITV and species turnover by forest type to weighted mean trait estimates. They can also be Fig. 5 Decomposition of the total variability in root and leaf traits for the tree communities in 164 1/16th-ha small plots in Jianfengling, Hainan China. The blue portion of the bars corresponds to the percent of variance explained by species turnover, while the red part shows the intraspecific variability effect-or the percent of variance explained by ITV. Black bars denote total variation (i.e., the variation or total sum or squares-SS, in specific weighted mean averages). The space between the top part of the column and bar corresponds to ITV and species turnover covariation; if the bar is above the column, the covariation is positive, and if below, the covariation is negative. The values are standardized by the total variation in specific CWM averages. The figure follows Fig. 3 in Lepš et al. (2011) compared to understand the effect of pGLM imputation on the relative contribution of ITV and species turnover by forest type to weighted mean trait estimates. Between the AWM and CWM analyses, ITV and species turnover had similar contributions to total functional trait variance for the seven traits considered (Fig. 5, left panel). The relative contribution of ITV in explaining the trait variation due to forest type increased markedly from the AWMs to CWMs (Fig. 5, middle panel). This decrease came as the contribution of ITV to trait flex ANOVA error decreased from the AWM analysis to the CWM analysis (Fig. 5, right panel). For example, when using only measured data (the AWM analysis), the total relative contribution of forest type for leaf area, SLA, and leaf thickness was 5, 11, and 15 percent, respectively (Fig. 5). The relative contribution of forest type for those three leaf traits increased to 26, 24, and 19 percent when pGLMimputed trait matrices were used in the CWM analysis (Fig. 5). Nearly, all this variation was due to ITV and not to species turnover (relative contribution < 1%) among plots. Moreover, the relative contribution of forest type to functional trait variation in root diameter, SRL, RTD, and RBI was 2, 19, 4, and 49 percent (Fig. 5), respectively, using only directly measured functional traits (AWMs). In contrast, the relative contribution to forest type increased to 12, 24, 9, and 77 percent when using pGLM-imputed data (CWMs).
Overall, the trait flex ANOVA models were strikingly similar between the AWM analysis and the CWM analysis (Fig. 5). For each of the seven traits considered, the trait flex ANOVA models had comparable total explanatory power (Fig. 5, left panel). Considering total weighted mean trait variation, five traits (SRL, SLA, root system diameter, leaf thickness, and leaf area) showed negative covariation for the AWM trait flex ANOVAs. For the CWM analysis, only three traits (SRL, root system diameter, and leaf area) showed negative covariation in relation to total trait variation; however, the magnitude of their negative covariation increased. For weighted mean trait variation due to forest type, both the AWM and CWM analyses were mostly attributable to ITV and not to compositional turnover, although the magnitude differed slightly, with slight differences in ITV and species turnover covariation (Fig. 5, middle panel). Lastly, trait flex ANOVA model errors were comparable between AWM and CWM analyses. The general pattern was that the error associated with ITV decreased when using pGLMimputed complete functional trait matrices as opposed to strictly using measured trait data. The magnitude of error attributable to species turnover and the covariation between ITV and species turnover was similar between AWM and CWM models (Fig. 5, right panel).

Discussion
In response to our first research question on whether species turnover is a more significant contributor to variation in CWM trait values than ITV-for six of seven traits (i.e., all traits except RBI), species turnover outweighed the contribution of ITV in explaining functional trait variation. However, trait variation due to forest successional status was nearly entirely explained by ITV, which was consistent across all root and leaf functional traits. This likely reflects both functional trait variation maximized by the sampling design because of forest successional status and actual ITV among sampled trees. Concerning our second research question, on the sensitivity of patterns of CWMs to phylogenetic trait imputation-we find some support for the PhyloPars method, consistent with recent work (Johnson et al. 2021). The inference was not entirely unchanged quantitatively by trait imputation which would support the null hypothesis of no difference in the relative contributions of species turnover and ITV when using directly measured versus imputed functional trait data; however, patterns were strikingly similar and changed very little qualitatively. This result speaks to the need to create ecologically relevant designs for the sampling of plant functional traits (Baraloto et al. 2010;Lepš et al. 2011) that capture ITV efficiently relative to measurement labor. The approach we demonstrate here using a plant lineage-based sampling design accounts for species turnover among communities and quantifies ITV using a priori knowledge of habitat variation across the landscape (i.e., soil fertility, texture, and forest age). Its application captured important differences in CWM trait estimates due to past logging history and current forest structure (i.e., light and soil environmental differences, Hogan et al. 2020)., which would have been overlooked using traditional analyses with species mean trait values.

The utility of habitat-structured functional trait measurements for assessing ITV in species-rich tropical forests
Functional perspectives (Keddy 1992;Westoby and Wright 2006) help frame complex and subtle changes in species and functional composition with the environment into ecologically meaningful information. This has been demonstrated regarding patterns of community assembly (Kraft et al. 2008;Baraloto et al. 2012;Spasojevic and Suding 2012;Laughlin 2014) and their underlying variation in plant traits (Letcher et al. 2015;Muscarella et al. 2016b). The data used in this study confirm that ITV in tropical tree communities is substantial (Paine et al. 2011;Poorter et al. 2018a;Hogan et al. 2020) and affects CWM trait estimates (Baraloto et al. 2010;Muscarella and Uriarte 2016). ITV may translate to vital physiological or ecological differences among individual trees, such as rates of photosynthesis and growth (Wright et al. 2010;Shipley et al. 2016), intrinsic rates of population growth (Adler et al. 2014), or variation in ecosystem productivity (Lohbeck et al. 2013;Poorter et al. 2017;van der Sande et al. 2018). Consistent with a burgeoning interest in ecologically accounting for species trait variance rather than the mean trait values, (Hulshof and Swenson 2010;Violle et al. 2012;Siefert et al. 2015), we find that ITV due to environmental variation in resource availability (i.e., soil nutrient content, see Hogan et al. 2020), as affected by forest successional state (Odum 1969;Christensen and Peet 1984;Guariguata and Ostertag 2001;Chazdon et al. 2010) influences CMW trait estimates.
Undoubtedly, some amount of ITV is related to environmental variation, but it remains difficult to quantify dimensions of niche space or specific ecosystem properties that systematically correlate with this variation (Lavorel and Garnier 2002;Shipley et al. 2016). Furthermore, there is the possibility of genetic variation and genetic × environment-interactive effects (Soliveres et al. 2014;Kunstler et al. 2016;Yang et al. 2018). Whether trait measuring campaigns focus on accurately measuring interspecific or intraspecific trait variation, the two approaches employ drastically different sampling designs (Baraloto et al. 2010). On the one hand, diversity differences are favored by a focus on trait variation among taxa. On the other hand, when ITV is the focus, sampling effort is maximized for one or a few species potentially at the expense of diversity-driven differences. In the tropics, the costs of this trade-off are particularly acute, as species diversity and environmental variation are both extensive. A second practical difficulty arises in that species compositional differences are often confounded with environmental variation (i.e., environmental filtering), making it challenging to assess ITV for all taxa, especially less common ones. Despite this, Lepš et al. (2011) propose an appealing compromise, grouping trait sampling by a priori ecological classification (i.e., measure habitat or experimental-treatment-specific trait values). This approach proved useful in assessing how fertilization and mowing shape the functional traits of European grasslands (Lepš 1999;Lepš et al. 2011).
In our case, the principal axis of environmental variation, irrespective of compositional differences among forest types, was a gradient in soil fertility and texture, with soil texture and fertility decreasing along the transect from more secondary to primary forest (Online Resource 1, Hogan et al. 2020). ITV along this gradient was consistent and measured 7 cm 2 for leaf area, 0.8 m 2 kg −1 for SLA, 0.4 mm for root diameter, 0.35 m g −1 for SRL, and 0.3 tips cm −1 , on average (Hogan et al. 2020). Additionally, there was considerable interfamilial variation among the taxa, which were preselected to capture a range in root functional strategies. We carefully sampled leaves and roots from individual juvenile trees that spanned environmental variation in soil environment and forest type. These considerations led to the trait data being well suited to scale to the community level using the forest demography data from the network of small plots (Online Resource 1, Xu et al. 2015a, b). Thus, the middle-of-the-road solution proposed by Lepš et al. (2011) of structuring functional trait sampling across known environmental gradients (e.g., light availability or soil nutrient gradients or gradients in land use or disturbance history) is also useful when carefully applied to forests.

Consequences of ITV for CWM estimation with tropical forest successional status
CWM patterns of trait variation can help to interpret ecological processes shaping community assembly from regional species pools (Kraft et al. 2008;Muscarella and Uriarte 2016;Spasojevic et al. 2016;Katabuchi et al. 2017). However, CWMs are gross generalizations of the aggregate function of species assemblages that may or may not reflect broader forest level or ecosystem functioning (e.g., net ecosystem exchange, nutrient uptake, or biomass production capacity of a stand; Warren et al. 2015;Fisher et al. 2018). Their accuracy depends on proper estimates of species traits in relation to the environment, which must account for both among and within species variation (Violle et al. 2012;Shipley et al. 2016;Yang et al. 2018). Measuring the morphology of all the plants within a community may one day be feasible but is currently not possible. Yet, the ability to accurately estimate trait values (e.g., leaf area, root-tip abundance) is critical to accurately model how vegetation will respond to future changes in the biosphere (Warren et al. 2015;McCormack et al. 2017;Fisher et al. 2018).
Scaling-up trait measurements from a 6.6-km transect to a whole forest ecosystem revealed that ITV compounded with compositional differences to create discrepancies in community-weighted trait values concerning forest successional status. That is, primary and secondary habitatspecific weighted means were much more divergent than weighted means not considering successionally structured ITV, which was consistent for both the AWM and CWM analyses. These differences would not have been detected had we not sampled functional traits along an environmental gradient consistent with the environmental variation found in the network of small plots (Figs. 3 and 4). This is due to the negative covariation of ITV and compositional turnover for most (six of seven) of the traits examined (for either the AWM or CWM analysis, Fig. 5). For example, along the JFL gradient, we measured a decrease of about 0.35 m g −1 in root system SRL from secondary to primary forest, which was consistent across 14 Angiosperm families, (Hogan et al. 2020). Accounting for this ITV, the specific mean difference in SRL between primary and secondary plots was 1.11 m g −1 for the AWMs and about 0.65 m g −1 when pGLM-imputed data were used (Fig. 4B). Covariation between ITV and species turnover for SRL was positive for the AWM estimate but was negative for the more complete pGLM-trait CWM analysis (Fig. 5). Regardless, like SRL, the difference in CWM trait estimates (this study) exceeded the magnitude of ITV measured across the JFL gradient (Hogan et al. 2020). This is because ITV usually compounds with β-diversity to shift CWM trait values toward the conservative end of the plant economics spectrum (Tables 1, Figs. 3 and 4) as species composition shifts toward more conservative species in older forest (Lohbeck et al. 2013;Letcher et al. 2015;Muscarella and Uriarte 2016;Spasojevic et al. 2016;Caplan et al. 2019).
Several studies have reported negative covariation between species turnover and ITV for a variety of functional traits (grass leaf traits in response to mowing and leaf nitrogen in response to fertilization: Lepš et al. 2011, leaf traits: Fajardo and Siefert 2018, leaf traits: Kichenin et al. 2013, SLA, leaf phosphorus, and leaf thickness: Derroire et al. 2018, plant height and leaf traits excepting SLA: Zuo et al. 2017). Negative trait covariances between drivers of trait variation (e.g., ITV because of environment or disturbance and species turnover in our case) confer resistance in the change of CWM trait values to a single driver, because the effect of one driver limits the trait response to the other driver. This dampening can be understood in the context of the portfolio effect (Schindler et al. 2015), where negative covariance between system components (e.g., drivers of community trait variation) leads to dampened dynamics at higher-order levels (e.g., the community) (Bolnick et al. 2011). The magnitude of the negative covariation in species turnover and ITV depends on the strength of species sorting (i.e., relative abundance) along environmental gradients versus the degree to which the environmental gradient affects ITV (Muscarella and Uriarte 2016;Denelle et al. 2019). On the other hand, positive covariation between species turnover and ITV has been found in some other studies from an equal variety of traits (leaf nitrogen: Kichenin et al. 2013, leaf area, leaf dry matter content, leaf nitrogen: Derroire et al. 2018, SLA: Zuo et al. 2017. The idiosyncrasy of these results points to the context dependency in the covariation of ITV in decomposing the functional trait variation of plant communities (Kichenin et al. 2013). In the AWM analysis we conducted, we found slight negative covariation for SLA, RTD, root system diameter, leaf thickness, and leaf area (Fig. 5). Extending the analysis to the entire community using pGLM-imputed functional trait matrices created some minor variation in these patterns. The sign of the covariation because of forest type changed for three of the seven traits (leaf area, leaf thickness, and SRL). Thus, although these changes were small, trait covariation seems sensitive to the assemblage of species for which functional traits are measured, potentially pointing to the role of species-specific responses in CWM trait patterns (Muscarella and Uriarte 2016;Katabuchi et al. 2017).
With respect to tropical forests, differing patterns in the contribution of ITV versus species turnover and their covariation may emerge for wet and dry tropical forests (Poorter et al. 2019). We report more negative than positive covariation in ITV and species turnover in CWM trait estimates. However, positive covariation of ITV and turnover was reported for most leaf traits in the Guanacaste dry forest of Costa Rica (Derroire et al. 2018). Patterns of β-diversity among sampled communities in relation to the functional uniqueness of the species contributing to those patterns (i.e., those lost or gained across assemblages) (Spasojevic et al. 2016) drives the covariation of the drivers of trait variation within and across plant communities. In communities with high levels of functional redundancy, such as tropical forests, responses may be muted as individual species play a less impactful role in driving overall patterns (Baraloto et al. 2010;Laughlin 2014;Umaña et al. 2017).
Author contribution JAH and CB designed the study, JAH and HX collected the data, and JAH analyzed the data and wrote the manuscript. All authors contributed intellectually and aided in the editing of manuscript drafts.
Funding Partial financial support was received from CTFS-ForestGEO at the Smithsonian. We are grateful for many small personal donations that funded the soil analyses (http:// www. exper iment. com/ china roots).

Data availability
The functional trait data used in this study are archived on figshare: https:// doi. org/ 10. 6084/ m9. figsh are. 79963 28. v4. The tree demographic data are available from Dr. Han Xu at the Chinese Academy of Forestry upon request.
Code availability All code used is freely available as part of R. The code for the trait flex ANOVA function is available in the supplemental material of Lepš et al. (2011). Example code for the analyses and creation of figures is available from JAH upon request.

Conflict of interest
The authors have no competing financial or nonfinancial interests to disclose.