Overview
We 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 area
Our 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 data
The 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 diameter are only recorded within a small portion of the plot, called the microplot. Our study design prevented the inclusion of the microplot and therefore excludes trees beneath this diameter threshold. Trees that grew into the measurement size class between the previous and most recent measurement are included. The exclusion of small trees and saplings may result in a conservative estimate of the difference between edge and interior BA and BAI, as other studies have found a higher density of small-stemmed woody vegetation along forest edges35. BA is calculated from a single plot measurement, as the summed BA of all live adult trees (>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 <2.5% and >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 growth
Abiotic controls
To 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 composition
Tree 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 selection
All 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 harvest
In 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 < 0.001), increase in harvested BA within both all-edge and anthropogenic edges as compared with the forest interior. We note that the exclusion of small-diameter trees from our study could alter these results if there was differential mortality between edge and interior in smaller tree size classes.
Temperate forests are heavily impacted by forest management49. We tested the robustness of the effect of edges on growth and biomass by withholding all subplots with a record of anthropogenic removals on the full FIA plot (i.e., management; n = 3642). We found no difference in the overall effect of edges nor meaningful differences within forest type groups.
Scaling edge effects on forest growth across the Northeast
Ecoregions are a widely used geographic partitioning of ecosystems into coherent spatial units as defined by abiotic, biotic, and anthropogenic characteristics28. EPA Level IV ecoregions are delineated by differences in environmental characteristics analogous to those that we used to model forest growth and thus are a comparable spatial unit to quantify the aggregated effects of fragmentation.
Quantifying fragmentation
To quantify anthropogenic forest edge area, we identify forest cover within 30 m of a road, development, or agricultural field (sensu Smith et al.6) using a 30 m resolution land-cover product from 2016 of the National Land Cover Database (NLCD)50. Edge forest was defined as all forest pixels adjacent (queen’s rule) to a non-forest cultivated or developed pixel (Supplementary Fig. 6a). Figure 3a shows the percentage of total forest area classified as edge within each ecoregion. We report that 18% of the total forest area in our study domain is adjacent to an anthropogenic edge. Differences from the reported 22% in Smith et al. are likely attributable to the use of a different NLCD product. Note that the definition of forest edge here may differ from that of the FIA analysis, given the constraints on quantification of the distance of edge influence and the spatial resolution of the land cover products.
Ecoregion edge impacts
To scale the effects as illustrated in Fig. 3, we quantified ecoregion forest composition by (1) using a 250 m resolution USFS forest type map27, we aggregated raw forest type values to the aggregated forest type groups included in our regression models (Figs. 3b), (2) we calculated the total area of each forest type group within each ecoregion, then used the average temperature, light, water, and nitrogen deposition in each ecoregion as inputs to our GLM regression models to calculate the BAI of edge and interior forest for each forest type. With the proportional area of each forest type, we calculated an area-weighted mean and then differenced the estimated edge and interior BAI to produce an expected difference of forest growth (BAI m2 ha−1) between edge and interior within each ecoregion (Supplementary Fig. 5). Finally, we combined the proportion of edge forest with the expected growth difference to quantify the estimated difference in percent increases in ecoregion BAI within each ecoregion attributable to increases of forest growth at the edge (Fig. 3c).
Quantifying global forest fragmentation
We quantified the extent of forest fragmentation throughout temperate and tropical forests worldwide at the scale of ecoregions using the Hansen Global Forest Change (v1.7)51 dataset on Google Earth Engine (GEE)52. Tropical and temperate biomes were delineated in a global ecoregion map53, analogous to the more detailed ecoregions described earlier. The tree canopy cover layer from the Hansen dataset provided estimates of percent tree canopy cover for the year 2000 at 30 m resolution globally produced by time series analysis of Landsat images51. To calculate the percentage of edge forest in each ecoregion: (1) a 10% threshold (following the FIA definition of minimum forest cover19) was applied to the tree canopy cover layer to separate forest and non-forest pixels, (2) each forest pixel adjacent (queen’s rule) to a non-forest pixel was classified as edge forest on GEE, and (3) ArcGIS Zonal Statistics Tool was used to calculate the percentage of edge forest in each ecoregion. Definitions of forest cover via % canopy cover vary between studies, therefore we performed a robustness check on our results to the threshold definition of forest cover by re-analyzing with a 30% canopy threshold. While there were differences in the calculated raw area of forest edges, the ratio of area fragmented between temperate and tropical forests did not change meaningfully (Supplementary Fig. 7). We then compared the Hansen-derived forest fragmentation to the 2016 NLCD-derived forest fragmentation used in our previous analysis to assess comparability of the two products. Supplementary Fig. 8 shows the agreement between the percent edge forest values calculated based on the two forest maps for the 247 ecoregions in the Northeast US. The agreement is strong especially in large and more forested ecoregions. The Hansen-derived percent edge forest explained 84.5% of the variance in NLCD-derived percent edge forest with RMSE of 6.1 (%) at ecoregion level. The spatial aggregation to ecoregion level largely reduced the uncertainty in the mapping of forest pixels in both products.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
