Biodiversity and productivity in eastern US forests

In all,188 tree species were recorded in the eastern US inventory plots (mean = 7.1, SD = 2.6, range = 3 to 20), comprising 157 angiosperm and 31 gymnosperm species. Trait values varied widely and were only weakly correlated with each other (unweighted mean trait correlations across species ranged from −0.01 to 0.56, SI Appendix, Fig. S1A; basal area-weighted trait correlations across plots ranged from −0.01 to 0.55; SI Appendix, Fig. S1B). For the species included in our analyses, gymnosperms had lower specific leaf area (SLA), rooting depth, and leaf nitrogen (LN) content than angiosperms; wood density (WD) and maximum height did not differ between the two groups (SI Appendix, Figs. S1A, S2, and S3).

To quantify functional trait and phylogenetic diversity, we considered metrics that are intrinsically independent of species richness (SP)—mean pairwise trait distance weighted by basal area (MTD) and mean pairwise phylogenetic divergence (branch length) weighted by basal area (MBL). We also considered two commonly used metrics that mathematically depend strongly on species richness—functional attribute diversity (FAD) (30) and Faith’s phylogenetic diversity (PD) (19). Across plots, species richness was positively correlated with FAD (r = 0.98) and PD (r = 0.81) but negatively correlated with MTD (r = −0.69) and MBL (r = −0.61) (Fig. 2 and SI Appendix, Fig. S4), which was also true within each ecoprovince (SI Appendix, Fig. S5). The strong positive relationship of FAD and PD with species richness is expected from the mathematical formulation of these indices, which are intrinsically correlated with species richness, and thus each other (Materials and Methods, Measuring Diversity).

To understand the negative correlation between species richness and MTD—and MBL to the extent that phylogenetic diversity reflects functional trait diversity (r = 0.77)—we examined how species richness related to species packing. Functional convex hull volume [functional richness, FRic (31)] measures the total volume that the species assemblages occupy in functional trait space. As expected, hull volume increased as the number of species increased (SI Appendix, Fig. S6A). However, as hull volume increased, MTD decreased (SI Appendix, Fig. S6B). Together, these patterns indicate that although increasing species richness does tend to add some additional species beyond the edges of the trait space, most species additions occur in the interior of the trait space. That is, as species richness increases, species tended to be more tightly packed in functional trait space.

Productivity (aboveground biomass growth rate, G) increased with species richness, FAD, and PD but decreased with MTD and MBL (Fig. 2 and SI Appendix, Figs. S4, S5, and S7). Mixed-effects regression models explained ~21% of the variation in productivity, with diversity metrics explaining a highly significant, albeit small (up to ~1.7%), percent of the variation that could not be explained by ecoprovince, stand age, and initial biomass (partial R² for diversity effects, which is measured as the differences in R² between mixed-effects regression models including vs. not including diversity metrics; Table 1 and SI Appendix, Tables S1 and S3). Although each diversity measure had statistically significant independent effects on productivity (Fig. 3, Table 2, and SI Appendix, Tables S2 and S4), any single measure explained nearly as much variation as when measures were combined with each other (Table 1 and SI Appendix, Tables S1 and S2). The low model R² values are characteristic of the large range in plot-level biomass growth in the inventory data (32, 33) and likely reflect multiple factors. The small size of the inventory plots (four subplots, each ~168 m², distributed over ~0.4 ha) leads to large sample variability. Furthermore, there is considerable spatial and temporal environmental heterogeneity within eastern US ecoprovinces that cannot effectively be captured by available broad-scale soil and meteorological datasets (34, 35), which constrains the explained variation in our analysis.

While productivity always had a positive relationship with species richness, the relationship of productivity to functional trait and phylogenetic diversity depended on the metric and the model (Fig. 3 and Table 2). For phylogenetic diversity, PD always had a positive relationship, and MBL always had a negative relationship regardless of the model. For functional trait diversity, MTD always had a negative relationship. In contrast, FAD had a positive relationship when considered singly or combined with PD, but a negative relationship when combined with species richness. This discrepancy highlights the importance of considering mathematical relationships between diversity metrics and species richness, which can influence the magnitude and sign of correlations between a given biodiversity measure and ecosystem functioning.

Stand age, initial plot biomass (at the beginning of the interval during which productivity was measured), and ecoprovince explained a substantial portion of the variation in productivity (Table 1). This is unsurprising given the large range in plot stand ages and initial biomass (SI Appendix, Fig. S8) and the broad environmental differences among eastern US ecoprovinces (29). It also raises the issue of whether the general trends concerning the various biodiversity metrics interacted with these factors. Therefore, we conducted additional analyses that grouped plots by stand age and initial biomass bins and by ecoprovince.

Within stand age bins [subsets of FIA (Forest Inventory and Analysis) plots with similar stand age], the models explained from 17 to 29% of the variation in productivity (G). Species richness was negatively related to productivity in young stands (<30 y) and positively related in older stands (>50 y, SI Appendix, Fig. S8A and Table S6). The overall positive relationship for the entire dataset (Fig. 3) likely occurs because most of the stands (67.0%) are in the older group. In contrast, the effect of MTD was always negative or not significantly different from zero. For MBL, the effect was not different from zero for all but one age class (80 to 90 y), where it was negative.

Within initial biomass bins (subsets of FIA plots with similar initial biomass), the models explained 14 to 21% of the variation in productivity (G). The effect of biodiversity was the strongest in plots with the lowest initial biomass (<50 Mg, partial R² = 6.8%, SI Appendix, Fig. S8B and Table S7). Within each initial biomass bin, the sign of biodiversity effects was consistent with the results based on all plots (Fig. 3). The effects of species richness were always positive or not different from zero. In contrast, the effects of MTD and MBL were always negative or did not differ from zero. There were no obvious trends in the sign or magnitude of diversity effects across initial biomass bins (SI Appendix, Fig. S8B).

When the plots were divided by ecoprovince, the direction of effects of biodiversity on productivity (G) varied (SI Appendix, Fig. S9A and Table S8). Overall, the models explained 9 to 27% of the variation in productivity (G). While the results were generally similar to the broad-scale analysis that spanned the entire eastern United States (Fig. 3), some differences were found. Results for species richness and MBL were mostly consistent between the broad-scale and ecoprovince-scale analyses. Within ecoprovinces, species richness and MBL effects were typically opposite one another, being positive for richness and negative for MBL (as in the broad-scale analysis) or not different from zero. The one exception was the Southern Mixed Forest, where the relationships were reversed—negative for species richness and positive for MBL. Results for MTD were less consistent between the broad-scale and ecoprovince-scale analyses. Effects of MTD on productivity (G) were significantly negative in three ecoprovinces and significantly positive in two, which contrasts with the negative MTD effects in the broad-scale analysis (Fig. 3B). The inconsistent signs of MTD effects across ecoprovinces may reflect differences in community assembly processes, which may warrant additional exploration.

Finally, although SP has no intrinsic correlation with MTD and MBL, the variances of MTD and MBL tend to decrease as SP increases [see Materials and Methods and an analogous result in Fig. 2E of Laliberté and Legendre (36)]. To evaluate the potential effects of this heteroskedasticity on our results, we repeated our regression analysis within groups of FIA plots separated into SP bins, which constrains the SP range and thus the degree of heteroskedasticity (each bin was roughly three species wide; see SI Appendix, Table S9). The positive dependence of productivity on SP and the negative dependence of productivity on MTD and MBL, as observed in the analysis of all FIA plots combined (Fig. 3B), were also observed within most SP bins (SI Appendix, Fig. S9B and Table S9). Eight of the 15 diversity effects (three metrics in five SP bins) were nonsignificant (SI Appendix, Table S9), which may be due to the limited sample size and limited range of diversity (especially for SP) within each bin. To further evaluate the potential effects of heteroskedasticity, we translated the observed MTD and MBL in each FIA plot to standardized effect size (SES) values relative to a null model, and we repeated our regression analysis with these SES values substituted for the original MTD and MBL values. In this null model–adjusted analysis, productivity had a significant positive dependence on SP and a nonsignificant negative dependence on MTD and MBL (SI Appendix, Fig. S10). Together, these analyses suggest limited issues with heteroskedasticity, a robust positive dependence of productivity on SP, and a negative or nonsignificant dependence of productivity on MTD and MBL (Fig. 3B and SI Appendix, Figs. S9B and S10).

Species with missing trait data generally accounted for a small proportion within any given FIA plot (less than 5% of the species or less than 6% of the total basal area in a plot, SI Appendix, Table S10). Even for this well-studied flora and set of ecologically important traits, we were missing about 20% of all species-trait combinations, albeit mostly for rare species. The effect of missing data on functional trait diversity was further reduced for the metrics weighted by basal area as missing values were mostly associated with less abundant species (Eq. 2). The exception was rooting depth (R_max), which was missing for approximately 29% of all species or 23% of the total basal area. Due to this large percentage of missing data, we did not include R_max when calculating trait diversity (MTD and FAD) in our primary analyses presented above. To evaluate the potential effects of excluding an important trait, we repeated the analyses with R_max included (with weights in Eq. 1 set to 0 for missing data). The results with R_max included (SI Appendix, Tables S1 and S2) were generally the same as the results without R_max (Tables 1 and 2). Finally, to evaluate the effect of variation in the number of traits used for distance measures of different species on estimates of species distances in trait space, we imputed missing values (Materials and Methods, Missing Data). The results with imputation (SI Appendix, Tables S3 and S4) were generally the same as the results without imputation (Tables 1 and 2). Because both sets of additional analyses were nearly identical, we conclude that our results are robust to missing data.

We quantified relationships among forest productivity, functional trait diversity, PD, and species richness in nonplantation forests of the eastern United States. Our analysis is based on 1,821,107 individual tree measurements in 23,145 systematically sampled forest inventory plots that were measured and remeasured from ~2000 to 2020 (Fig. 1). These forest inventory data were combined with publicly available trait data and a phylogeny derived from the Open Tree of Life. As an index of productivity, we calculated the rate of aboveground biomass production (Mg/ha/y) from plot remeasurements; this biomass growth rate (G) comprises ~40% of total net primary production in temperate forests (63) and accounts for ~70% of changes in total ecosystem carbon stocks in US forests over recent decades (64). We detail the assembly of these data products below.

Our analyses used national-scale, systematically sampled forest inventory data from the FIA program of the United States Department of Agriculture Forest Service. Because eastern US forests are considerably more diverse both functionally and phylogenetically than western US forests (65), and because FIA remeasurement intervals are longer in the western United States (~10 y) than in the eastern United States (~5 y), we restricted our analyses to plots from the eastern United States (Fig. 1). Each plot consists of four 7.32-m-radius subplots, one centrally located and the other three spaced 36.6 m apart in a triangular arrangement. Within each plot, trees with a diameter at breast height (DBH) greater than 12.7 cm are inventoried. Trees with a DBH 2.54 to 12.7 cm are inventoried in 2.07-m-radius “microplots” (one per subplot). Data reported for each inventoried individual include DBH, height, species identity, aboveground biomass, and other tree-level variables (66).

We used FIA database version 9.0, downloaded in November 2022 from https://www.fia.fs.usda.gov/. We restricted our analyses to plot locations that were measured at least twice (i.e., at least one remeasurement) using the current national standardized annual sampling design (67), which was implemented in most US states beginning around 2000. We further restricted our analyses to plot measurements on accessible forest lands that were classified into a single condition (e.g., stand age and soil type) and not classified as artificially regenerated plantations. We did not include plots incapable of supporting at least 1.40 m³/ha/y of wood volume growth (which excludes plots assigned to the lowest of FIA’s seven site productivity classes) or plots that lost more than 20% of their standing biomass stock to natural mortality or timber harvest between measurements (see SI Appendix, Table S5 for details). The dataset we analyzed included 62,698 surveys at 23,145 plot locations. Each plot location was assigned to an ecoprovince by aggregating the ecological subsection codes reported by FIA to the province level (29). Each ecoprovince (Fig. 1) represents an area of similar geology, soil type, and climate that supports a similar potential natural vegetation type.

For each plot, we calculated the plot-level aboveground biomass growth rate (G; Mg/ha/y) of trees with DBH ≥ 2.54 using a mass balance approach, described in detail below. In our dataset, 11,114 plots were measured twice (resulting in one G estimate), 8,143 plots were measured three times (two G estimates), 3,399 plots were measured four times (three G estimates), and 489 plots were measured five times (four G estimates). For plots with multiple G estimates, we used the mean value in subsequent analyses so that each plot was represented just once in each analysis.

To estimate G, we first calculated the plot-level live aboveground biomass stock (B_t; Mg/ha) of trees with DBH ≥ 2.54 for each plot at each time t by summing the product of individual tree aboveground biomass estimates (b_i for tree i) and their expansion factors (TPHA_i, the number of trees per hectare that a tallied tree represents) for all trees (i = 1 to n) alive at time t: $B_t={\sum }_{i=1}^n\left(b_i\times {\text{TPHA}}_i\right)$   . Data on b_i and TPHA_i are in the FIA Tree table for all trees with DBH ≥ 2.54 (66), where b_i (DRYBIO_AG from the FIA Tree table) is estimated by combining DBH measurements with allometries and the component-ratio-method (66, 68), and TPHA_i is calculated as the inverse of the tree’s sampling area. In our analysis, b_i is converted from pounds to metric tons (Mg), and the TPHA values are a constant (2.47 acre/ha) times the “trees per acre unadjusted” values reported in the FIA Tree table. We then estimated G based on the biomass dynamics of a plot between times t and t+Δt: B_t+Δt − B_t = Δt·G − M_t − C_t, where M_t and C_t, respectively, are the time-t biomass stocks of trees that died or were cut (harvested) between times t and t+Δt (69). Rearranging yields the estimate of G over the interval t to t+Δt: G =(B_t+Δt − B_t + M_t +C_t)/Δt.

We selected five functional traits for analysis that were available for most eastern US tree species. Four of the traits—WD (g/cm³), LN (mg/g), plant maximum height (H_max, m), and SLA (mm²/mg)—are included in the global spectrum of plant form and function (39). For supplementary analyses, we also included a fifth trait, maximum rooting depth (R_max, m); greater variability in rooting depth among species is expected to increase the efficiency of belowground resource uptake [e.g., shallow vs. deep-soil water and nutrients (70)]. WD, LN, and SLA values were obtained from the TRY database (27) (accessed in May 2022). R_max values also were obtained from TRY and augmented using the latest data from Tumber‐Dávila et al. (71) and Guerrero-Ramírez et al. (72). WD values were augmented using data reported by FIA (WOOD_SPGR_GREENVOL_DRYWT in the FIA REF_SPECIES table). H_max was estimated for each species using individual tree heights reported by FIA. Specifically, we defined H_max as the 95th and 99th percentile of heights [HT in the FIA Tree table (66)] reported for a given species. The H_max values based on the 95th and 99th percentiles were highly correlated with each other (r = 0.98); thus, only the 99th percentile of heights was used in measuring H_max in further analyses. Because not all species had measurements of all traits, we evaluated the impact of missing data by performing additional analyses (see below).

The trait values of a species vary in space, which may bias trait distance measurements in each FIA plot. To account for variation in traits within species among plots, we used a geographical buffer prior to calculating trait distances among species in each FIA plot, making use of the fact that multiple measures of georeferenced trait values existed in the data sources that we used. Specifically, we set a circular spatial buffer (moving window) for each FIA plot location, using the “st_buffer” function in the sf R package (73), and then intersected the FIA plot locations (using the approximate coordinates reported in the public FIA data) with the georeferenced trait records using the “st_intersection” function. The geographical coordinates of trait records were obtained from FIA plots (for H_max) and the TRY database (for the other four traits). The trait value of a species in a plot was assigned the mean value of all trait records located within the buffer. The diameter of the buffer ranged from 100 to 900 km. The buffer size varied by trait and was determined by maximizing trait coverage (i.e., at least 1 georeferenced trait record located within the buffer of a plot) for all FIA plots at the finest possible geographic scale. Specifically, WD of a species in a given plot was assigned the mean value of WD of this species in all plots within a 900-km diameter buffer; the buffers for LN, SLA, R_max, and H_max were 600 km, 400 km, 500 km, and 100 km, respectively. After the buffer process, some trait data of a certain species may be missing in a FIA plot because the trait data are outside the buffer of the FIA plot. In these instances, we replaced missing values with the mean values of all available trait data for the species. Distances in trait space among species coexisting in the same plot were calculated following Gower (74):

where D_ij is the trait distance between species i and j, and x_ki and x_kj are the values of trait k of species i and j, respectively. x_k.max and x_k.min are the maximum and minimum values of trait k, respectively, and w_kij is the weight of trait k for the species i-j pair, which is 1 if both species have values for trait k and 0 if trait k is missing in either species i or j.

Phylogenetic trees for the species in FIA plots (SI Appendix, Fig. S3) were extracted from the ALLMB tree presented by Smith & Brown (28), which uses GenBank data and the Open Tree of Life (52, 53) with a phylogenetic backbone provided by Magallón et al. (75). Prior to calculating diversity metrics, species names from FIA plots, trait databases, and the phylogeny were standardized using the Taxonomic Name Resolution Service V5.0 (https://tnrs.biendata.org/, accessed: 6 April 2022), following Tropicos (https://www.tropicos.org/) and World Flora Online (http://www.worldfloraonline.org). Pairwise branch distances were then calculated.

Because site productivity is the sum of the productivity of each individual (1, 76), functional trait diversity and phylogenetic diversity were both calculated using a basal area-weighted approach, where the weight for each tree was its DBH² multiplied by its TPHA value (see above); this is equivalent to weighting each tree by its basal area. MTD and MBL were calculated as:

where X is either MTD or MBL, d_ij is either the pairwise phylogenetic distance (branch length) or trait distance (Eq. 1) between species i and j, n_i is the total basal area of species i, N is the total basal area of all species, and S is species richness.

We also calculated two commonly used metrics, Faith’s PD (19), the sum of all branch lengths [in million years (Ma)]:

where L_ib is the length of the bth branch segment of the ith species on the phylogenetic tree, and FAD (30), the sum of the pairwise trait distances:

Because both PD and FAD are sums of branch lengths or trait distances, they inherently contain species richness as a component.

Unless stated otherwise, the sampling unit in all analyses was an inventory plot location. All variables were calculated for each plot remeasurement and then averaged across all remeasurements at each location, as described above for G. For example, for a plot measured three times (t₁, t₂, and t₃), all diversity metrics and other covariates were calculated at t₂ and t₃ and then averaged. Thus, each plot location occurred once in each analysis. First, we calculated pairwise correlations among growth and the diversity metrics and among the SD of the five traits weighted by basal area. Then, we explored the association between growth and the diversity metrics using linear mixed-effects models. To control for stand history and broad edaphic-climatic factors, each model included a random intercept term for ecoprovince (ECO) and fixed effects for stand age (STD), initial biomass (IniBio), and between one and three diversity components: species richness (SP), phylogenetic diversity (PD or MBL), and functional trait diversity (MTD or FAD), or the basal area-weighted SD of a single trait. For the single-trait models, the basal area-weighted SD was determined for each trait in each plot as the square root of the basal area-weighted variance, which was calculated using the "wtd.var" function in the Hmisc package version 5.1-0 (77) in R, with basal area weights as described above. All numeric variables (i.e., all variables except ecoprovince) were standardized to zero mean and unit variance prior to fitting the models so that the model coefficients (slopes) are standardized to a common, unitless scale. We compared model fits using R² and AIC. We recognize that biodiversity–productivity relationships often are not linear over their entire range (78). We used linear model forms, however, because they can easily accommodate random effects and are straightforward to interpret in terms of the sign and magnitude of effects and because preliminary analyses using generalized additive (nonlinear) models yielded qualitatively similar inferences. The robustness of our results to linear vs. nonlinear forms likely indicates that most FIA plots fall within the linear portion of the biodiversity–productivity relationships.

For additional exploration of the effects of stand age, the initial biomass of the plots, and environmental heterogeneity as captured by ecoprovince, we divided the plots into groups of different stand ages (10-y intervals) or initial biomass (50-Mg intervals). We then repeated linear mixed-effects models for each group. For each group of stand ages, biodiversity metrics (i.e., SP, MBL, and MTD) and initial biomass are fixed factors, and ecoprovinces are random factors; for each group of initial biomass, biodiversity metrics and stand age are fixed factors, and ecoprovinces are random factors. We also explored biodiversity–productivity relationships within each ecoprovince. For each ecoprovince, we conducted ordinary least squares regressions, with species richness (SP), phylogenetic diversity (MBL), trait diversity (MTD), stand age (STD), and IniBio included as independent variables.

To account for missing data, we did three sets of analyses: 1) excluding the trait with substantial missing information, R_max, 2) including R_max because it is known to be important for productivity in trees, and 3) excluding R_max but imputing missing values for the other four traits. We imputed the missing values for LN and SLA using the “phylopars” function in Rphylopars version 0.3.9 (79). The missing values were replaced by restricted maximum likelihood values that were estimated based on the trait values of other species and the phylogeny under a Brownian evolution model. The results of the first set are reported in the main text (Tables 1 and 2). The other two sets (SI Appendix, Tables S1–S4) provide tests of the robustness of our conclusions by either including a potentially important functional trait or by estimating missing values so that all pairwise measures of trait distances are based on the complete set of traits. The number of species with missing values (no data for a given trait), out of the 188 species, was 0 for WD, 54 for LN, 21 for SLA, 108 for maximum rooting depth (R_max), and 0 for maximum height (H_max). For each trait, we measured the average percent of species without data and the average percent of basal area (m²/ha) of these species across FIA plots. The average percentage across plots for each trait was calculated as the average percentage among species within each plot and then averaged across all plot locations.

The form of Eq. 2, a mean value, ensures that the MBL and MTD metrics that we used for phylogenetic and functional trait diversity, respectively, have no intrinsic correlation with species richness (SP). Thus, any correlations that arise in any given empirical dataset will be due to processes affecting community assembly.

An additional issue is heteroskedasticity in the relationships between species SP and other diversity metrics. If all species in the regional species pool are present at a single site, SP takes on its maximum possible value, and both MBL and MTD metrics become highly constrained as the only variance is due to differences in species abundances. Conversely, as the number of species in a site declines, the range of possible values increases for both MBL and MTD [see results for a related metric, FDis, in Fig. 2E of Laliberté and Legendre (36)]. This heteroskedasticity might affect our conclusions, depending on the amount of heteroskedasticity over the observed range of SP. In our data, SP values for single inventories of FIA plots ranged from 3 to 21, out of a total of 188 species in the eastern US FIA dataset. Because all variables, including SP, were averaged over multiple inventories of each FIA plot to calculate means for subsequent statistical analyses, the SP values used in our analysis were often not integers and ranged from 3.0 to 19.5.

To explore the potential effects of heteroskedasticity on our results, we conducted two analyses. The first analysis involved implementing the mixed-effects model reported in Fig. 3B for groups of FIA plots that were separated into the following SP bins: (3, 4.5), (4.5, 7.5), (7.5,10.5), (10.5,13.5), and (13.5,19.5), with the last bin being wider due to the small number of FIA plots with SP > 13.5 (see sample sizes in SI Appendix, Table S9). The mixed-effects model reported in Fig. 3B was fit separately to the FIA plots in each bin so that the SP range and potential heteroskedasticity were constrained. The second analysis involved translating the observed MBL and MTD values into SES values using a null model and then using the SES values in place of the observed MBL and MTD as explanatory variables in the mixed-effects model reported in Fig. 3B. The null model was based on simulated communities following the methods of Laliberté and Legendre (36). Specifically, we created 10,000 simulated communities for each observed SP value from 3 to 21 as follows: In each simulated community, each species identity was randomly assigned from the pool of 188 species, and each abundance was randomly drawn from a lognormal distribution; the log-scale parameters were mean 0 and SD 1, but the results are insensitive to the lognormal parameter values because numeric variables were standardized prior to fitting the mixed-effects model. We then paired each FIA plot with the set of simulated communities whose SP level matched the rounded (nearest integer) SP value of the FIA plot, and we calculated the FIA plot’s SES for MTD as (observed MTD value − mean of simulated MTD values)/(SD of simulated MTD values). We calculated SES values for MBL in the analogous way. Finally, we fit the same mixed-effects model reported in Fig. 3B, except using the SES values for each FIA plot as explanatory variables instead of the observed MTD and MBL values.

Y.L., J.A.H., J.W.L., R.P.G., D.E.S., P.S.S., and S.M.S. designed research; Y.L. and J.A.H. performed research; Y.L., J.W.L., R.P.G., and S.M.S. analyzed data; and Y.L., J.A.H., J.W.L., R.P.G., D.E.S., P.S.S., and S.M.S. wrote the paper.

The authors declare no competing interest.