Intraspecic trait variation and species functional turnover in successional tropical forests: assessing gap-lling for community weighted means

Accounting for intraspecic trait variation (ITV) is central to plant ecology and crucial for vegetation modeling efforts. ITV can be substantial; however, it remains unclear how ITV inuences 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 small plots comprising 582 species to evaluate the contribution of ITV to CWMs, comparing unlogged, primary forest to selectively-logged and clear-cut secondary forest. We examine the effect of gap-lling missing trait values via phylogenetic generalized linear modeling (PhyloPars) on CWMs. For six of seven traits, ITV negatively covaried with species turnover to generate larger CWM differences than observed if ITV was not integrated. For example, plot average CWM specic leaf area was 10.7 and 10.4 m 2 kg − 1 for primary and secondary forest, not accounting for ITV, but shifted to 9.8 and 11.1 m 2 kg − 1 after doing so. Specic root length showed a similar trend. Our results from 72-species assemblages were supported by the results from the gap-lled analysis using the entire community, where the contribution of ITV to CWMs ranged from 25 to 75%, with nearly all trait variation due to forest type attributable to ITV. Therefore, CWM trait estimates became more-conservative with forest age, whereas ITV for many traits showed an acquisitive shift, and because of negative covariation between ITV and species turnover, forest age-related CWM differences increased. Differences were unaffected, if not strengthened, by gap-lling incomplete functional trait matrices.


Introduction
Plant functional traits broadly scale along a fast-slow axis related to plant life-history variation and demographic rates (Díaz et al. 2016, Reich, 2014, Wright et al. 2010. For example, the most broadly measured plant functional traits distill to two main axes of functional variation -one related to leaf economics (e.g., speci c 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 functional trait variation within taxa (i.e., intraspeci c trait variation, hereafter ITV, Siefert et al. 2015) that muddles trends, for example weakening trait-demography relationships (Poorter et al. 2018a, Yang et al. 2018. Additionally, a portion of ITV is due to the environment; thus, properly accounting for ITV across quanti ed environmental gradients can reveal trait-environment relationships (Bruelheide et al. 2018, Craven et al. 2018) and drivers of variance among them (Hofhansl et al. 2021).
Functional trait variation also has implications for plant tness and how species partition environmental resources within niche space (Laughlin et al. 2020, Laughlin and Messier 2015, Kunstler et al. 2016. Thus, at the community level, functional traits have popularly been used to infer community shifts in plant resourceacquisition strategies across environmental gradients (e.g., environmental ltering, community assembly, Escudero and Valladares, 2016). For example, in tropical forests, functional diversity accounts for differences in forest productivity more so than taxonomic diversity (van der Sande et al. 2018), richness (Rosen eld and Müller, 2020), or forest dynamics in space and time, especially regarding forest succession (Letcher et al. 2015, Lohbeck et al. 2012, Muscarella et al. 2016a, Hogan et al. 2018. Moreover, it is far easier to infer plant community dynamics or plant-environment relationships from a functional trait perspective than through changes in individual species or community (i.e., multivariate) demographics (Poorter et al. 2018a, Poorter et al. 2018b). This trend is exacerbated with increasing plant community richness (Baraloto et al. 2010). In hyper diverse tropical forests, where there can be tens or hundreds of species within genera, using a functional trait perspective often uncovers environmental associations between functional plant strategies (e.g., Kraft et al. 2008), trends which would be obscured by a traditional taxonomic approach. Moreover, when functional trait studies in hyper diverse tropical forests couple a functional perspective to phylogenetic data, patterns of inference can be strengthened (Baraloto et al. 2012). However, most of these studies use species-level trait data, ignoring ITV. Hence, there is motivation to understand whether incorporating 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, Xu and Trugman, 2021. Conceptually, ITV arises because of genetic variation among individuals, the environment, and the interaction between the two (Siefert et al. 2015, Bolnick et al. 2011. Accounting for ITV is labor-intensive, requiring extensive trait sampling of individuals; however, 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, yet on the other, the evolutionary perspective seeks to better understand how natural selection shapes communities, ecosystems, and their constituent biota (Bolnick et al. 2011, Turcotte and Levine, 2016, Violle et al. 2012. 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). Few studies have looked at ITV across successional or environmental gradients in tropical forests; moreover, none have sought to assess how such variation scales up to in uence estimates of community-weighted mean (CWM) trait patterns.
The lack of studies aimed at quantifying the contribution of ITV to CWM trait patterns is likely due to the laborious undertaking of measuring the functional traits for individuals and the community composition across environmental or successional gradients (Paine et al. 2011). This is particularly di cult in tropical forests, where large trees and high levels of diversity complicate trait sampling (Baraloto et al. 2010). Methodological approaches that have looked for ways around this issue by gap-lling or imputing trait values for missing taxa 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 rst proposed methods was multiple imputation chained equations (MICE, van Buuren and Groothuis-Oudshoorn 2010); however, this method has since been shown to be awed for several reasons, perhaps the most-supreme of which, in the context of plant functional traits, is a lack of phylogenetic consideration (Penone et al. 2014). Two of the next most-desirable methods include Bayesian hierarchal probabilistic matrix factorization (BHPMF, Schrodt et al. 2015), which can use a hierarchical taxonomic structure (e.g., genera within families within orders) to infer phylogenetic structure in the imputation of 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, thereby directly accounting for species relatedness via branch lengths. A recent study compared the three methods, concluding that PhyloPars outperforms MICE and BHPMF for trait imputation (Johnson et al. 2021).
Forest communities typically move from acquisitive to conservative functional composition as species composition changes with forest succession (Christensen and Peet 1984, Letcher et al. 2015, Lohbeck et al. 2013, Muscarella et al. 2016a, Swenson et al. 2012, Garnier et al. 2004). Thus, CWM functional trait values shift toward traits associated with species conservativeness (or slow life-history strategy, e.g., SLA and higher wood density -WD), partly because of turnover in community composition alone (Chazdon et al. 2010). The contribution of ITV in a successional context is less well understood; however, contributions can be substantial. For example, ITV of key leaf traits (i.e., leaf area, SLA, leaf nitrogen content, and leaf phosphorus content) accounted for >25% to about 75% of the variation in CWM leaf trait values in a temperate rainforest in Chile (Fajardo and Siefert 2018). Root functional traits are of interest because they are abstruse and ecologically complex (Bardgett et al. 2014, Laliberté, 2017, rarely studied in the context of ITV, and control ecosystem function via the ow of nutrients and water through trees (Freschet et al. 2020). Thus, understanding how CWM root traits vary in the context of environmental variation and forest age is potentially key to modeling the functioning of tropical forests into the future (Bardgett et al. 2014, McCormack et al. 2017, Warren et al. 2015, Few studies have examined how CWM patterns in root traits change with forest succession and how they may be different or similar to leaf traits. A recent study from a temperate forest that used BHMPF to gap-ll missing root trait data (Caplan et al. 2019) found that CWMs shifted toward conservative root strategies (thicker diameters, lower speci c root length, higher root tissue density, and lower root nitrogen content) with increasing forest age and that acquisitive (e.g., pioneer) species and invasive species with acquisitive root strategies were able to persist in disturbed communities. A different study sampled a strong edaphic chronosequence in temperate rainforest in New Zealand and found that ITV in root traits showed similar patterns to CWM patterns in root functional traits (Holdaway et al. 2011). The differences in these two studies highlight the need for more studies that quantify the degree to which soil environment-structured variation, forest succession-driven variation, and other ITV sources have on in uencing CWM estimates of root functional traits. Here we use root and leaf trait data collected in a paired sampling design along a 6.6-km transect in a tropical forest in Hainan, China, to scale up trait measurements to the community level using a complementary network of small vegetation plots (Xu et al. 2015a, Xu et al. 2015b Speci cally, we address two research questions: 1) to what extent do species turnover (i.e., changes in community composition because of forest age and environmental ltering effects), ITV, and their covariation in uence community-weighted mean functional trait estimates? 2) how sensitive are these patterns to phylogenetic trait imputation for missing taxa (i.e., gap-lling of missing trait data). Concerning question 1, we expected that changes in plant community composition and species relative abundances (i.e., species turnover) with forest age would in uence estimates of community-weighted traits to a greater degree than ITV because the magnitude of interspeci c trait variation would be more important if ITV were substantial. For question 2, the testable null hypothesis is no difference in the relative contributions of species turnover and ITV between using directly measured and imputed functional trait data. Accepting the null hypothesis here would validate trait imputation methods and provide a way forward for future functional trait measurement strategies. two-thirds of the area was either clear-cut or selectively-logged, meaning 30-40% of large timber-valuable trees were extracted (Xu et al. 2015b, Zhou 1995. All logging ceased in 1994 under a state-wide (Hainan Island only) logging ban, followed by a Chinese national logging ban in 1998 (Wenhua 2004, Zhou 1995. The reserve encompasses several ecological life zones of vegetation, from tropical semi-deciduous monsoon forest at the lower elevations to mossy high elevation 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 (per the Holdridge life zone classi cation), 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).

Materials And Methods
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 at JFL has averaged 2700 mm (Wu 1995). The soils are derived from porphyritic granite and are classi ed 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 other tropical soils (e.g., latisols, Wu 1995).

Root and leaf functional trait data
Functional traits were collected from May to July 2017. Saplings (i.e., trees >1m tall, but less than <10 cm diameter at 1.3m) of 72 species from 15 selected Angiosperm families, comprising 423 saplings (roughly 6 per species), were measured for leaf and root functional traits. Trait collection from saplings was done systematically along a 6.6 km transect that spanned primary and secondary forest with logging history ( These individual-level trait data can be condensed to a few axes of key morphological variation among roots and leaves. A principal components analysis (Hogan et al. 2020) showed that three leaf traits -leaf area, speci c leaf area (SLA), and leaf thickness, and four root traits -root system average diameter, speci c root length (SRL), root tissue density (RTD), and root branching intensity (RBI) accounted for most of the two to three dimensional above and belowground trait variation among individuals. Therefore, we also 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, Xu et al. 2015b). Between August 2007 and June 2009, a network of 164 1/16 th -ha (25m × 25m) plots were established in the JFL reserve (Xu et al. 2015b, Fig. 1). The plots span the local variability in vegetative life-zones and range in altitude from 259 to 1135 m. 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 that was either selectively logged or clear cut, with human land-use ceasing between 20 and 56 years ago (see Table 1 in Xu et al. 2015b). Selective-logging practices in the area typically result 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 identi ed to species.

Data analyses: the network of small plots and Lepš' Intraspeci c variability effects
In the interest of scaling up our functional trait measurements from individual trees along the gradient to the scope of the JFL forest reserve 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) in lieu of the classical community-weighted mean (CWM) because we could not measure traits for all species found in the plant communities. However, we did employ a sampling design that controlled for plant lineage, targeting species from 15 families across the Angiosperm phylogeny. This sampling design targeted families representing the functional variation in root morphologies across the angiosperm phylogeny and sampled many individuals within each species (average of 6) to assess ITV, where ITV was found to be structured across forest types (Hogan et al. 2020). 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.
Basal area (measured in August 2007 to June 2009) in the array of small plots ranged 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 mean trait values from all individuals for each species sampled along the transect (x i ), we calculated AWMs for each plot. We term these the xed averages (after Lepš et al.

2011), thus:
where p i is the basal area fraction of that species within the plot (i.e., the community weight). 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-speci c trait values for each species (x i_habitat ). For example, for plots in secondary forests, trait values from saplings sampled in secondary forest. We term these the speci c averages (again after Lepš et al. 2011), as such The intraspeci c variability effects for the assemblages were calculated as the difference between the xed and speci c averages. Two-way analyses of variance (ANOVAs) assuming homoscedasticity were conducted with respect to each of the seven functional traits chosen. Then using the sum of squares of these tests, the variance in trait values was decomposed into intraspeci c and community (i.e., species turnover) portions and their covariation (i.e., interaction). This was done using the 'trait ex ANOVA' function (see supplemental material of Lepš et al. 2011) in R v.3.5.0 (R Core Team, 2018)

Data analyses: gap-lling traits and validation of assemblage-weighted means
To validate the results regarding intraspeci c variability effects using AWMs, we employed an approach to gap-ll trait data for unsampled taxa. As we have explained, PhyloPars (Bruggeman et al. 2009) appears to be one of the best tools available for gap lling 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 estimates missing data values.
A rudimentary phylogeny for all 582 taxa in the network of small plots was obtained from Phylomatic v3 (www.phylodiversity.net/phylomatic, 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, in that we wanted to have one trait that was imputed for all taxa. WD is a highly phylogenetically-conserved trait (Kraft et al. 2010, Swenson andEnquist, 2007), which can be used in the pGLM to constrain modeled trait values for taxa and clades with no trait data. WD values were obtained from the global database of WD values ) where available, and the CTFS Forest-GEO WD database (http://ctfs.si.edu/Public/Datasets/CTFSWoodDensity/) via the 'BIOMASS' package ('getWoodDensity' function, Réjou-Méchain et al. 2017). Where species WD values were unavailable, values from a higher taxonomic classi cation (e.g. genus or family) were used (see Online Resource 3).
WD data was combined with three separate incomplete trait matrices, one using xed trait values (x i , i.e., species mean values for all individuals measured along the transect), and two habitat-speci c 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 t separately to each incomplete trait matrix and predicted to complete them for all 582 taxa using PhyloPars (http://zeus.few.vu.nl/programs/phylopars, Bruggeman et al. 2009) and the Phylomatic-generated phylogeny. For the model setting in PhyloPars, we allowed for correlated evolution of the different traits, but did not allow for intraspeci c variation in the models because, in our case, the separate tting of the pGLMs accounts for this. We then used the gap-lled trait matrices and basal area species weights (p i ) of species in the small plots to recalculate xed and speci c averages (i.e., weighted means) for each plot. We re-ran the trait ex ANOVAs, comparing the AWM analysis, using only the species for which we measured functional traits (i.e., 72 species), to the CWM analysis, which used the PhyloPars-imputed trait data for all 580 species in the communities of the small plots.

Results
Linking the transect to the network of small plots -species forest-type preferences and soils comparison Of the 72 species that were sampled for functional traits along the 6.6 km transect, plot community composition showed that 65% of them (47 species) preferred primary forest. In comparison, the remainder (25 species) favored plots in secondary forest (Fig. 2). The basal area in primary and secondary forest plots was roughly equal, averaging 40.5 (± 0.2, standard error) m 2 ha -1 for small plots in primary forest and 42.7 (± 0.2) m 2 ha -1 for those in secondary forest. Soils in secondary forest were signi cantly ner in texture than soils from the primary forest t (% sand: F 1,298 = 63.9, p < .001, % silt: F 1,298 = 53.8, p < .001, % clay: F 1,298 = 38.5, p < .001) and signi cantly less acidic (average pH of 2.5:1 soil in water of 4.39 vs. 4.27, F 1,298 = 14.1, p < .001). Soils in the secondary forest had signi cantly greater total exchangeable base content (F 1,298 = 35.7, p < .001), effective cation exchange capacity (F 1,298 = 11.62, p < .001), and base saturation (F (1,298) = 22.5, p < .001) than the soils in the primary forest. These trends from the soils of the transect mirrored trends in soil properties between small plots in secondary and primary forests (Online Resource 1). Thus, the range of variability in soil conditions between the transect and the network of small plots was similar, making comparisons between functional trait measurements from the transect to the community composition of the network of small plots unbiased relative to the soil environment.
Estimates of assemblage-weighted trait values -xed vs. speci c means Our analytical approach involved calculating AWM estimates of functional traits using all data irrespective of forest type (i.e., xed means) and using only data from each forest type (i.e., speci c means, see complete analyses in Online Resource 2). Fixed AWM values for leaf functional traits differed drastically from speci c AWMs (Table 1, Fig. 3). For example, plot average xed AWMs for leaf area were 77.5 cm 2 and 78.2 cm 2 (a 0.7 cm 2 difference) for secondary and primary forest, respectively. Those estimates shifted to 70.5 cm 2 and 81.5 cm 2 (a difference of 11 cm 2 ) for speci c AWMs (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, speci c 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 those for leaf area and SLA. AWM differences were not signi cantly different for xed averages but were for speci c averages (Table 1), being about 0.02 mm greater in secondary forest than in primary forest (Fig. 3C).
Similarly, for root functional traits, xed AWM values were similar and statistically indistinguishable (in all cases) between primary and secondary forest plots. At the same time, speci c AWMs showed apparent and often statistically signi cant differences (Fig. 4). For example, xed and speci c AWMs for average root system diameter ( Fig. 4A) showed an opposite trend, with xed estimates being slightly (<0.01 mm) more narrow in primary than secondary forest, but with speci c estimates being marginally thicker (0.03 mm) in primary than secondary forest. For both the xed and speci c AWMs for root diameter, these differences were not signi cant (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 speci c 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 speci c 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 xed AWMs for SRL. These differences went from being statistically non-signi cant for the xed AWM comparison by forest type to signi cant for the speci c AWMs (Table 1). The results for RTD were similar to those for average root system diameter (Fig. 4C). Fixed AWMs for RTD showed a decreasing trend from primary to secondary forest, with roots being about 0.02 g cm -3 less dense in secondary than primary forest. Speci c AWMs for RTD showed the opposing trend, with AWMs for RTD being about 0.3 g cm -3 denser in secondary than in primary forest.
Differences between xed AWMs were not statistically different, but the difference between speci c AWMs were statistically signi cant (Table 1). Lastly, AWM RBI showed no signi cant differences when speci c means were compared across forest types, with values averaging about 1.97 root tips cm -1 ( Figure 4D, Table 1). However, xed 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 xed AWM values by forest type for or any of the seven functional traits considered was for RBI (Table 1, Figs. 3 and 4), but statistical differences in speci c AWM trait estimates occurred for all traits except for root diameter, which had a marginally signi cant difference (p = .08, Table 1). Even when differences between speci c AWMs between forest types were not statistically signi cant, they were always greater than differences in xed 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 con rmed the results of the AWM analyses. Like the AWM analysis, the pattern of having larger differences in habitat-speci c means than for xed means between forest types was consistent across all leaf (Fig. 3) and root (Fig. 4) 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% con dence interval breadth contracted). This decrease in estimate variability led to more statistically signi cant differences for speci c CWMs relative to speci c AWMs (Table 2). For leaf traits (leaf area, SLA, and leaf thickness), no qualitative difference emerged from CWMs using pGLM-imputed trait data for taxa, versus only measured trait data (Fig. 3, Table 2). The only noticeable 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 xed and speci c means for both forest types. Additionally, the difference between xed CWMs for RBI was not statistically signi cant (Table 2), where it was for the xed AWMs for RBI (Table 1). Also, the difference in speci c CWMs for root diameter was statistically signi cant (Table 2), where the difference between AWMs was marginally signi cant (Table 1).
Comparing AMWs to CWMs for root traits, we see similar patterns. Using AWMs resulted in no differences between xed and speci c estimates for root diameter (Fig 4A). However, when the pGLM-imputed data were used, the speci c estimates for root diameter increased slightly for primary forest and decreased for secondary forest, creating a statistically signi cant difference between forest-types (Table 2). For the three other root traits (SRL, RTD, and RBI), differences were accentuated slightly, as pGLM-imputed functional trait matrices were used (Fig. 4 B, C, D). This discrepancy in forest-type weighted-trait estimates resulted more from ITV with environment than from species turnover among plots (Table 2, Fig. 5). In many cases, for example for SLA, root diameter, SRL, and RTD, ITV and species turnover had negative covariance values (Table 4, Table S4), meaning that beta diversity in species composition drives weighted-trait values one way (in our case, usually slightly toward the conservative end of the plant economics spectrum, but not always), while ITV with environment causes shifts in the opposite direction (in our case, roots morphologies became more conservative while leaf morphologies were more acquisitive from primary to secondary forest, Hogan et al. 2020).
Interspeci c variability effects -the relative effect of ITV and turnover on weighted-mean trait estimates Results from the interspeci c variability effects from AWMs (Table 1) and CWMs (Table 2) can be used to infer the relative contribution of ITV and species turnover by forest type to functional trait variability. They can also be compared to understand the effect of pGLM gap-lling on the relative contribution of ITV and species turnover by forest type to functional trait variability. 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-ex ANOVA error decreased from the AWM analysis to the CWM analysis (Fig. 5, right panel). For example, when using only measured data (i.e., 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, respectively, when pGLM-imputed trait matrices were used in the CWM analysis (Fig. 5). Nearly all this variation was due to ITV and not to species turnover (i.e., β-diversity, relative contributions < 1%) among the small vegetation 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 (Table 4), respectively, using only directly-measured functional traits (i.e., AWMs). In contrast, the relative contribution to forest-type increased to 12, 24, 9, and 77 percent, respectively, when using pGLM-imputed data (i.e., CWMs).
Overall, the trait-ex ANOVA model were strikingly similar between the AWM analysis and the CWM analysis (Fig.  5). For each of the seven traits considered, the trait ex ANOVA models had comparable total explanatory power (Fig. 5, left panel). Considering total weighted-mean trait variation, ve traits (SRL, SLA, root system diameter, leaf thickness, and leaf area) showed negative covariation for the AWM trait-ex 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 analysis were mostly attributable to ITV and not to turnover, although the magnitude differed slightly, with small differences in ITV and species turnover covariation (Fig. 5, middle panel).
Lastly, trait ex ANOVA model errors were comparable between AWM and CWM analyses. The general pattern is that the error associate with ITV decreased when using data from the entire plant community (i.e., pGLM-imputed complete functional trait matrices as opposed to the 72 species assemblages). 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
The utility of functional perspectives to species-rich tropical forests The application of functional perspectives (Keddy 1992, Westoby andWright 2006) to understand the ecology of plants in relation to environmental variation is a promising way to reduce complex and subtle changes in species composition and variation in functional composition of communities into ecologically meaningful information. This has been demonstrated regarding patterns of community assembly (Baraloto et al. 2012, Kraft et al. 2008, Laughlin, 2014, Spasojevic and Suding, 2012 and their underlying variability in plant traits (Letcher et al. 2015, Muscarella et al. 2016b. The data used in this study con rm that ITV in tropical tree communities is substantial (Hogan et al. 2020), which has a consequence for CWM trait estimates (Baraloto et al. 2010, Muscarella andUriarte, 2016). ITV has been shown to account for vital ecological differences among species, such as rates of photosynthesis and growth (Shipley et al. 2006, Wright et al. 2010, intrinsic rates of population growth (Adler et al. 2014), or variation in ecosystem productivity (Lohbeck et al. 2013, van der Sande et al. 2018 Yet, there has been a burgeoning in examining the role of ITV in functional ecology (i.e., ecologically accounting for species trait variance rather than the mean trait values; Siefert et al. 2015, Violle et al. 2012, Hulshof and Swenson 2010. Undoubtedly, some amount of ITV is related to environmental variation, but di culty exists in quantifying dimensions of niche space or speci c ecosystem properties that systematically relate to this variation (Lavorel and Garnier 2002, Shipley et al. 2016. Furthermore, there is the possibility of genetic variation and genetic environment interactive effects (Kunstler et al. 2016, Soliveres et al. 2014, Yang et al. 2018). Whether to focus on interspeci c or intraspeci c 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 in among taxa trait variation. On the other, when ITV is the focus, sampling effort is maximized for one to 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 di culty arises in that species compositional differences are often confounded with environmental variation (i.e., environmental ltering), making it challenging to assess ITV for all taxa, especially less-common ones. Despite this, Lepš et al. (2011) propose a middle-of-the-road solution that appeals, group trait sampling by a priori ecological classi cation (i.e., measure habitat or experimental-treatment speci c 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. We employed a similar approach to assess how environmental variation in resource availability (i.e., soil nutrient content and light availability, see Hogan et al. 2020), roughly summarized by the successional progression of a forest (Odum 1969, Christensen and Peet 1984, Guariguata and Ostertag 2001, Chazdon et al. 2010 leads to functional and compositional community change. 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 collect 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, Xu et al. 2015b). Thus, the middle-of-the-road solution proposed by Lepš et al. (2011) of structuring functional trait sampling across known environmental variability is useful when carefully applied to forests, as well.

Consequences of ITV for CWM estimation with tropical forest successional status
On one-hand interpreting, CWM patterns of trait variation can help interpret ecological processes shaping community assembly from regional species pools (i.e., environmental ltering of species or even phenotypes) (Kraft et al. 2008, Muscarella and Uriarte, 2016, Spasojevic et al. 2016). On the other, they are gross generalizations of the aggregate function of species assemblages that may or may not re ect broader forest-level or ecosystem functioning (e.g., total photosynthetic, nutrient uptake, or biomass production capacity of a stand; Warren et al. 2015. 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, or more so even, belowground traits like root-tip abundance) is critical to accurately modeling how vegetation will respond to future changes in the biosphere (Warren et al. 2015, McCormack et al. 2017).
Scaling up trait measurements from along the JFL gradient to the forest using data from a network of small plots revealed that ITV compounded with compositional differences (i.e., species basal area weights) to create discrepancies in community-weighted trait values concerning forest successional status. That is, speci c weighted-means were much more different than xed weighted-means, 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 plant lineage-consistent decrease of about 0.35 m g -1 in root-system SRL from secondary to primary forest (Hogan et al. 2020). Accounting for such ITV, the speci c 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 was positive for the assemblage-weighted mean estimate but was negative for the more-complete analysis (Fig. 5). Regardless, like SRL, the difference in CWM trait estimates exceeded the magnitude of ITV measured across the JFL gradient (Hogan et al. 2020). This is because ITV compounds with β-diversity to move CWM trait values toward the conservative end of the plant economics spectrum (Tables 1 and 2, 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. 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, leaf thickness: Derroire et al. 2018, plant height and leaf traits excepting SLA: Zuo et al. 2017). Negative trait covariances between drivers of trait variability (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). 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-speci c responses in CWM trait patterns (Muscarella and Uriarte 2016).
With respect to tropical forests, differing patterns in the contribution of ITV versus species turnover and their covariation and may emerge for wet and dry tropical forests (Poorter et al. 2019). We report more negative than positive covariation in ITV and specie 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).
Additionally, Zuo et al. (2017) found positive covariation between trait variation because of species turnover and ITV for several leaf traits, excepting SLA, in dune grassland communities in Mongolia.
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 variability 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 signi cant role in driving overall patterns (Baraloto et al. 2010, Laughlin, 2014, Umaña et al. 2017).

Conclusions
In response to our rst research question on whether species turnover is a more signi cant contributor to variation 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 type was nearly entirely explained by ITV, which was consistent across all root and leaf functional traits. This likely re ects both functional trait variation maximized by the sampling design because of forest type and actual ITV among tropical trees. Concerning our second research question, on the sensitivity of patterns of CWMs to gap-lling via phylogenetic trait imputation, we nd support for the PhyloPars method, which is in congruence with some recent work (Johnson et al. 2021). Inference was not unchanged quantitively by trait imputation (i.e., accepting the null hypothesis as outlined in the introduction), however patters were strikingly similar (i.e., 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 e ciently relative to measurement labor.
The approach we demonstrate here accounts for trait variation across the community using phylogenetic sampling, as well as for ITV across gradients of environmental variation, like soil fertility, texture, and forest age.

Declarations
Authors' contributions statement: JAH and CB designed the study, JAH and HX collected the data, JAH analyzed the data and wrote the manuscript. All authors contributed intellectually and aided in the editing of manuscript drafts.
following people who helped in the eld and with the processing of plant samples: Shaojun Ling, Yaxin Xie, Jaming Wang, Siqi Yang, Wenguang Tang, Shitaing Ma, Qiqi Zhang, and Jiazhu Shi. Professor Yu from the Jianfengling Forest Bureau provided taxonomic eld identi cation assistance.
Funding: JAH received support via a short-term fellowship from CTFS-ForestGEO at the Smithsonian. We are grateful for many small personal donations that helped fund the soil analyses (www.experiment.com/chinaroots).
Con icts of interest: The authors no competing nancial or non-nancial interests to disclose. values use all functional trait data from individuals sampled along the transect (n = 423), whereas speci c assemblage-weighted means match functional trait data from plant individuals to plots with corresponding foresttype classi cation (n = 198 for primary forest and n = 225 for secondary forest). Intraspeci c variability = xedspeci c. Abbreviations are as follows: df = degrees of freedom, SS = sum of squares, F = F-statistic, p = probability of obtaining the F-statistic. Probabilities are denoted as: *p < .05, **p < .01, ***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,  Statistics on residual variance are not shown. These ANOVAs were done using small plot tree community data and functional trait data for 582 species (i.e., the entire community for which we applied pGLM data imputation to). 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 speci c values use pGLM-imputed functional traits data corresponding to either secondary or primary forest portions of the transect.
Intraspeci c variability = xed -speci c. Abbreviations are as follows: df = degrees of freedom, SS = sum of squares, F = F-statistic, p = probability of obtaining the F-statistic. Probabilities are denoted as: *p < .05, **p < .01, ***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.  Figure 1 The study site: the Jianfengling Forest Reserve (JFL) of Hainan Island, China. Hainan Island is a small continental island off 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 were sampled is shown in black. The rst 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 1/16th-ha plots (Xu et al. 2015a, Xu et al. 2015b) 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 has happened in the last 20-55 years. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Figure 2
Relative stem densities (% of stems ha-1) for the 72 species in 164 1/16th-ha plots in Jianfengling, Hainan Island, China ( Fig. 1, Xu et al. 2015b), for which trait data were sampled (Hogan et al. 2020). Fifty-two of those plots are in well-conserved, primary forest, while 112 are in secondary forest (i.e., previously logged or clear-cut). 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): Lauraceae, Magnoliaceae, Annonaceae, Moraceae, Fagaceae, Juglandaceae, Burseraceae, Anacardiaceae, Sapindaceae, Rutaceae, Pentaphylacaceae, Sapotaceae, Ebenaceae, Theaceae, Styracaceae. 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 intraspeci c 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 speci c 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 speci c CWM averages. The gure follows Figure 3 in Lepš et al. (2011).

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. S1.pdf S2.pdf S3.pdf