Elevated growth and biomass along temperate forest edges
OverviewWe used data from the national forest inventory conducted by the US Department of Agriculture, Forest Service, Forest Inventory and Analysis (FIA) program to quantify tree biomass and growth along forest edges and within the forest interior. We estimated the causal impact of the forest edge environment on patterns of tree biomass and growth, while accounting for potentially confounding variables. We then used the regression models to estimate the aggregate difference in growth attributable to forest edges throughout the northeastern U.S. Finally, to better understand the implications of our findings, we quantified the degree of forest fragmentation throughout temperate and tropical forest biomes world-wide, using a 30 m forest cover map.Study areaOur analyses of edge impacts on forest biomass and growth were conducted throughout twenty-states (1.7 million km2) in the northeastern and upper mid-west of the United States (Supplementary Fig. 1). This region contains 765,000 km2 of forest and encompasses gradients of dominant land-uses, climatic conditions, and forest composition while remaining within deciduous, coniferous, and mixed temperate forest ecosystems.Identifying edges in forest inventory dataThe FIA collects measurements of tree size, growth, and land-use within a nested plot design across the country19. Each FIA plot is composed of four individual subplots; within each subplot, the diameter at breast height (dbh) of every tree >12.7 cm is measured during each measurement period. The re-measurement frequency for FIA plots in our study area is between 5 and 7 years, but this can differ between Forest Service regions. In addition to tree measurements, the database details land-use condition data that includes the proportion of the area that is forested and, on some plots, the land-cover class of the non-forest area (FIA User’s Manual, Condition Table). FIA plots are considered forested if some portion of the plot includes a contiguous forest patch (including potentially outside of the plot area) of greater than 4047 m2 that has more than 10% canopy cover. With a memorandum of understanding between the USFS and Harvard University, we had access to the true, unfuzzed plot coordinates, which are not publicly available. Evaluating >48,000 plots in the USFS Northern Region sampled from 2010 to 2020 and selecting the most recent measurement cycle for each plot, we identified subplots that contained both a forest and a non-forest condition and categorized these as edges (Supplementary Table 1). Only subplots that included a forest condition in both the most recent and previous measurement were included. Subplots where the mapped condition changed from forest to non-forest were excluded. Changes in the amount of mapped forest condition were included and are incorporated into the calculation of response variables using the most recent condition area. We identified FIA plots where all four subplots were fully forested as interior plots to be used for comparison. Subplots located within the same plot as an edge subplot (i.e., edge-proximate subplots) were excluded from this study due to limitations in our ability to quantify their distance from an edge. The spatial configuration of subplots is such that a fully forested subplot may be up to ~65 m away from an identified forest edge within another subplot. Studies suggest that the distance of edge influence in temperate forest does not extend more than 30 m into the forest interior15,33. Since the FIA does not contain information about the geometry of non-forest conditions beyond the subplot boundary, we deemed that the large uncertainty in the relationship between these subplots to a non-forest edge precluded their inclusion in the study. The FIA plot configuration prevented quantification of the distance of edge influence in our analysis; the exclusion of subplots adjacent to edge-subplots may limit direct comparisons with other fragmentation studies.We used the FIA condition data to characterize the non-forest land use in edge subplots. Information on adjacent non-forest land cover is not collected on all FIA plots (4327 of 6607 edge subplots). We aggregated FIA land-cover classification to a binary anthropogenic or unknown edge type designation and present results from all edge subplots and the anthropogenic edge subset (FIA User’s Manual Condition Table, Section 2.4.50).For each subplot (168 m2 in area), we calculated two primary response variables of interest: total live tree BA and BAI. Notably, trees smaller than 12.7 cm dbh) in m2. BAI was calculated on a per-tree basis as the difference in radial growth of live adult trees between the most recent and previous measurements, and then divided by the number of years between measurements (m2 yr−1). In addition, we aggregated individual tree diameter measurements to calculate mean stem density (stems ha−1) and mean tree diameter for each subplot (Fig. 2).We accounted for variable subplot area by normalizing both BA and BAI to a per-hectare of forested area basis, resulting in units of m2 ha−1 and m2 ha−1 yr−1, respectively. To account for potential small-area bias, we performed a sensitivity analysis on the relationship between BA and subplot forested area (Supplementary Fig. 2). We subsequently excluded 1284 subplots under 30 m2 in area as the area to BA relationship asymptotes relationship above this threshold. Finally, we accounted for errors in field dbh measurements, sometimes resulting in negative BAI values, by excluding the 97.5% quantiles of both BA and BAI distributions.Given their spatial configuration, FIA subplots are not fully independent measurements, potentially introducing issues with pseudo-replication and spatial autocorrelation within our dataset. To test for spatial autocorrelation we examined the semivariance of model residuals36, and found that there was high correlation only at distances of less than 1 km. The spatial stratification of the FIA plot design minimizes issues of plot–plot proximity within our study. However, to account for autocorrelation between subplots, we filtered our pre-matched dataset to only including one subplot from each FIA plot. For plots containing multiple edge subplots, we selected the subplot with the largest forested area. For interior plots, we selected the central subplot and excluded all others.Isolating the effect of edges on growthAbiotic controlsTo account for environmental controls on forest growth we included the most critical abiotic predictors of terrestrial vegetation productivity (light, water, temperature, and nitrogen deposition) as covariates in the regression models (Supplementary Fig. 4, Supplementary Table 2). Light, water, and temperature data were drawn from spatial raster maps (0.5° resolution) as unit-less indices of relative limitation on vegetation productivity, ranging from 0 to 13. Nitrogen data were drawn from the 2018 NADP gridded inorganic wet nitrogen deposition product (4 km spatial resolution; kg of N ha−1)37. To interpolate across small gaps in the raster data (usually along water bodies), we used the Nibble tool from ArcGis Pro (ESRI Team). We then used FIA plot locations to extract values from each raster layer for all FIA subplots.Forest compositionTree species may vary in their responses to biogeochemical changes that occur on forest edges. Overall forest community response emerges from complex interactions between species. We used aggregations of tree species, termed forest composition groups (or forest types)38, to assess if species composition influenced the response to altered edge condition. Forest type classifications for each subplot are provided by the FIA (FIA User’s Manual, Condition Table) and are defined in Appendix D therein. We aggregated the FIA forest types into eight broader species groups, following Thompson et al.23, and defined in Supplementary Table 1.Matching, GLM regressions, and model selectionAll statistical analyses and most of the data processing were conducted in R, version 3.439. Using a causal inference framework, we created a quasi-experimental statistical design that included pre-matching followed by a GLM regression analysis40. Matching emulates an experimental design using observational data by identifying control groups of untreated (forest interior) plots that were as similar as possible to treated (forest edge) plots in terms of observable confounders. By capturing key differences in abiotic variables we control for the fundamental drivers of forest productivity, allowing for a direct estimation of the average treatment effect of edges. Similarity was defined by nearest-neighbor covariate matching determined by Malahanobis distance, implemented in the MatchIt library in R41, the simplest and best method when the dataset is robust enough to find a match for every treated plot20. This method excludes forest interior plots that are not matched with an edge plot. Given differences in sample size between the full edge dataset and the subset designated as anthropogenic edges, we performed matching separately on the two datasets. To assess the efficacy of matching on reducing the differences in covariate distributions, we used summary statistics calculated with the MatchIt library and report the pre- and post-matched covariate balance in Supplementary Table 4 and Supplementary Table 5 (sensu Schleicher et al.42). Matching was highly successful, largely eliminating differences in all covariate distributions in both datasets.Our primary response variables of interest, BA and BAI, were right-skewed, non-normally distributed and violated the assumptions of normality necessary for ordinary least squares regression43. We, therefore, used a GLM to better fit the structure of our data. GLMs are an extension of linear regression that allow more freedom in the choice of probability distribution function through the use of a link function to model relationships between predictors and response variables44. The gamma probability distribution is frequently chosen to model BA, given its assumptions of positive, continuous values and flexible model form23,45. We performed a series of GLM regressions on our post-matched datasets, using a gamma probability distribution with an inverse link function to model the relationship of BA and BA with a suite of predictor variables, using the glm function as implemented in the R Core stats package39. Due to differences in sample size between the all-edge dataset and the anthropogenic-edge subset, we modeled these two datasets separately for each of BA and BAI, resulting in four separate regression analyses. We used a model selection framework to identify the most parsimonious model within each of the model sets based on the Akaike Information Criterion (AIC) and residual deviance statistic46,47. We report the model-selection and model-fit results for each of our separate analyses, including model forms, AIC, Nagelkerke Pseudo-R2, and residual deviance in Supplementary Table 2. Across all four regression analyses, the best-performing model was one that included an interaction between the edge-status and forest type categorical variables, as well as the variables of temperature-limitation, light-limitation, water-limitation, and nitrogen deposition.We then used the best performing model from each analysis to compare the differences in BA and BAI between forest edge and interior across each forest type. We estimated the treatment effect of edge-state within each forest type using the ggeffects package48 to calculate marginal effects with the continuous predictors (temperature, light, water, and nitrogen deposition) held at their within-forest type regional means. The results of this analysis are displayed in Fig. 1 and Supplementary Table 3; primary error bars on the interior point show the 95% confidence interval of the marginal effect from the full edge model, while secondary error bars show the CI from the anthropogenic edge model. Due to the smaller sample size in the anthropogenic model, estimates of the mean marginal effect of the interior plots vary slightly (though non-significantly) from those from the full dataset. The main text description reports outputs from both models, calculated from separate interior mean estimates. For visual clarity, we only display one set of interior means in Fig. 1.Mortality and timber harvestIn tropical forests, large reductions in productivity along edges are associated with increased tree mortality.9 To assess differences in tree mortality across our study region, we applied a simplified GLM analysis, including edge-state as our only predictor variable. The FIA differentiates between mortality attributed to timber harvest and that attributed to other, non-harvest causes. The results of this analysis are presented as marginal effects of each edge category in Supplementary Fig. 3. There are no significant differences in biogenic mortality between edge groups and no difference in overall mortality (combined biogenic and anthropogenic); there is a small, but statistically significant (p More
