Over half of western United States' most abundant tree species in decline

Field observations

Since 1999, the FIA program has operated an extensive, nationally consistent forest inventory designed to monitor changes in forests across all lands in the US⁶¹. We used FIA data from 10 states in the continental western US (Washington, Oregon, California, Idaho, Montana, Utah, Nevada, Colorado, Arizona, and New Mexico) to quantify shifts in relative live tree density, excluding Wyoming due to a lack of repeated censuses (Fig. 1). This region spans a wide variety of climatic regimes and forest types, ranging from temperate rain forests of the coastal Pacific Northwest to pinyon-juniper woodlands of the interior southwest⁶². Although the spatial extent of the FIA plot network represents a large portion of the current range of all species examined in this study (Table 1), substantial portions of some species ranges (e.g., Douglas-fir) extend beyond the study region into Canada and/or Mexico and therefore were not fully addressed here.

The FIA program measures forest attributes on a network of permanent ground plots that are systematically distributed at a rate of ~1 plot per 2428 hectares across the US⁶¹. For trees, 12.7 cm DBH and larger, attributes (e.g., species, DBH, live/dead) are measured on a cluster of four 168 m² subplots⁶¹. Trees 2.54–12.7 cm DBH are measured on a microplot (13.5 m²) contained within each subplot, and rare events such as very large trees are measured on an optional macroplot (1012 m²) surrounding each subplot⁶¹. In the event a major disturbance (i.e., >1 acre in size, resulting in mortality or damage to >25% of trees) has occurred between measurements on a plot, FIA field crews record the primary disturbance agent (e.g., fire) and estimated year of the event. In the western US, one-tenth of ground plots are measured each year, with remeasurements first occurring in 2011. Please see Data Availability for more information on forest inventory data accessibility.

Forest stability index

Allometric relationships between size and density of live trees make it difficult to interpret many indices of forest change¹⁹. Live tree density is expected to decline as trees grow in size, owing to increased individual demand for resources and growing space (i.e., competition)^16,23. The expected magnitude of change in tree density, given some change in average tree size, varies considerably across forest communities⁶³, site conditions⁶⁴, and stand age classes²³. Thus, we posit it is useful to contextualize observed changes in live tree density relative to those expected given shifts in average tree size within a stand. To this end, we developed the FSI, a measure of change in relative live tree density that can be applied in stands of any forest community and/or structural type.

To compute the FSI, we first develop a model of maximum size-density relationships for tree populations in our study system (Fig. 6). This model describes the theoretical maximum live tree density (({N}_{max }); in terms of tree number per unit area) attainable in stands as a function of their average tree size ((overline{S})) and will be used as a reference curve to determine the proportionate live tree density of observed stands (i.e., observed density with respect to theoretical maximum density). We use average tree basal area as an index of tree size (one, however, could also use biomass, volume, or other indices of tree size). For stand-type i, the general form of the maximum tree size-density relationship is given by

where a is a scaling factor that describes the maximum tree density at (bar{S}=1) and r is a negative exponent controlling the decay in maximum tree density with increasing average tree size. Such allometric size-density relationships (i.e., power functions) are widely accepted as quantitative law describing the behavior of even-aged plant populations under self-thinning conditions^16,17, and have been used extensively to describe relative stand density in forests^23,24. As detailed below, we allow both a and r to vary with stand-type i, as maximum size-density relationships have been shown to vary across forest communities and ecological settings^63,65. Allowing a and r to vary by forest community type, for example, allows us to acknowledge that the maximum tree density attainable in a lodgepole pine stand is likely to differ from that of a pinyon-juniper stand with the same average tree size.

We next define an index of the relative density of a population of trees j (e.g., species, Pinus edulis) within a stand of type i (e.g., forest community type, pinyon/juniper woodland)

where N is the density represented by tree h (in terms of tree number per unit area), and S is an index of individual-tree size (e.g., basal area, as used here). The denominator of Equation (2) represents the maximum tree density attainable in a stand of type i with average tree size equal to the size of tree h. We therefore express the relative density of a population j within stand-type i as a sum of the relative densities represented by individual trees within the stand. RD can be interpreted as the proportionate density, or stocking, of a population of trees within stand, where values range from 0 (population j is not present within a stand) to 1 (population j constitutes 100% of a stand and the stand is at maximum theoretical density given its size distribution). As we do in this study, one may apply any range of estimators to summarize the expected relative density of a population of trees j across a range of different stand-types (e.g., estimate the mean and variance of RD_j across a broad region containing many different stand-types).

It is important to note that Equation (2) is approximately equal to a simpler method using aggregate indices (i.e., (frac{{sum }_{h = 1}{N}_{hij}}{{N}_{max }(overline{{S}_{i}})})) when tree size-distributions are normally distributed (even age-structures). However, the use of aggregate indices introduces class aggregation bias that results in overestimation of relative density in stands with non-normal size distributions (i.e., uneven age-structures), consistent with other indices of relative tree density⁶⁶. In contrast, summing tree-level relative densities eliminates such bias and allows RD to accurately compare density conditions across stands in very different structural settings (e.g., even-aged plantation vs. irregularly structured old forest). Furthermore, the partitioning of relative density into tree-level densities allows RD to be accurately summarized within tree size-classes⁶⁶. That is, it is possible to explicitly estimate the contribution of tree size-classes to overall stand density using RD.

For a given population j within stand-type i, we define the FSI as the average annual change in relative tree density observed between successive measurements of a stand

where Δt is the number of years between successive measurement times t₁ and t₂ and ΔRD is the change in RD over Δt (i.e., ({{rm{RD}}}_{{t}_{2}}-{{rm{RD}}}_{{t}_{1}})). The FSI may also be expressed in units of percent change (%FSI), where average annual change in relative tree density is standardized by previous relative density

Here, stability is defined by zero net change in relative tree density over time (i.e., FSI equal to zero), but does not imply zero change in absolute tree density or tree size distributions. For example, a population exhibiting a decrease in absolute tree density (e.g., trees per unit area) may be considered stable if such decline is offset by a compensatory increase in average tree size. However, populations exhibiting expansion (i.e., ({{rm{RD}}}_{{t}_{1}}<{,}{{rm{RD}}}_{{t}_{2}})) or decline (i.e., ({{rm{RD}}}_{{t}_{1}}> {,}{{rm{RD}}}_{{t}_{2}})) in relative tree density will be characterized by positive and negative FSI values, respectively.

Statistical analysis

We computed the FSI for all remeasured FIA plots in the western US (N = 24,229). We included plots on both public and private lands and considered all live stems (DBH ≥2.54 cm) in our analysis. As forest management can effect regional shifts in tree density, we excluded plots with evidence of recent (i.e., within 5 years of initial measurement) silvicultural treatment (e.g., harvesting, artificial regeneration, site preparation). All plot measurements occurred from 2001 to 2018, with an average remeasurement interval of 9.78 years (±0.005 years). For brevity, we restricted our analysis to consider the eight most abundant tree species in the western US. We identified the most abundant tree species using the rFIA R package⁶⁰, defining abundance in terms of estimated total number of trees (DBH ≥ 2.54 cm) in the year 2018. We excluded species that exhibit non-tree growth habits (i.e., shrub, subshrub) across portions of the study region. All statistical analysis was conducted in Program R (4.0.0)⁶⁷.

We developed a Bayesian quantile regression model to estimate maximum size-density relationships for stand-types observed within our study area. Here, we use TPH as an index of absolute tree density, average tree basal area ((overline{{rm{BA}}}); equivalent to tree basal area per hectare divided by TPH) as an index of average tree size, and forest community type to describe stand-types. We produced stand-level estimates of TPH and (overline{{rm{BA}}}) from the most recent measurements of FIA plots that (1) lack evidence of recent (within remeasurement period or preceding 5 years) disturbance and/or silvicultural treatment and (2) exhibit approximately normal tree diameter distributions (i.e., even-aged). Here we define an approximately normal tree diameter distribution as exhibiting Pearson’s moment coefficient of skewness between −1 and 1.

We transform the nonlinear size-density relationship to a linear function by taking the natural logarithm of TPH and (overline{{rm{BA}}}), and use a linear quantile mixed-effects model to estimate the 99th percentile of TPH conditional on (overline{{rm{BA}}}) (i.e., in log-log space) for all observed forest community types. We allowed both the model intercept and coefficient to vary across observed forest community types (i.e., random slope/intercept model), thereby acknowledging variation in the scaling factor (a) and exponent (r) of the maximum tree size-density relationship across stand-types. We place informative normal priors on the model intercept (μ = 7, σ = 1) and coefficient (μ = 0.8025, σ = 0.1) following the results of decades of previous research in maximum tree size-density relationships^16,23,63,65.

The FIA program uses post-stratification to improve precision and reduce non-response bias in estimates of forest variables⁶⁸, and we used these standard post-stratified estimators to estimate the mean and variance of the FSI for each species across their respective ranges within the study area (see Code Availability for all relevant code). Further, the FIA program uses an annual panel system to estimate current inventories and change, where inventory cycles consist of multiple panels, and individual panels are comprised of mutually exclusive subsets of ground plots measured in the same year within a region. Precision of point and change estimates can often be improved by combining annual panels within an inventory cycle (i.e., by augmenting current data with data collected previously). We used the simple moving average estimator implemented in the rFIA R package⁶⁰ to compute estimates from a series of eight annual panels (i.e., sets of plots remeasured in the same year) ranging from 2011 to 2018. The simple moving average estimator combines information from annual panels with equal weight (i.e., irrespective of time since remeasurement), thereby allowing us to characterize long-term patterns in relative density shifts. We determine populations to be stable if the 95% confidence intervals for range-averaged FSI included zero. Alternatively, if confidence intervals of range-averaged FSI do not include zero, we determine the population to be expanding when the estimate is positive and declining when the estimate is negative.

To identify changes in species-size distributions, we used the simple moving average estimator to estimate the mean and variance of the FSI by species and size class across the range of each species within our study area. We assign individual trees to size-classes representing 10% quantiles of observed diameter distributions (i.e., diameter at 1.37 m above ground) of each species growing on one of seven site productivity classes (i.e., inherent capacity of a site to grow crops of industrial wood). That is, we allow size class definitions to vary among species and along a gradient of site productivity, thereby acknowledging intra-specific variation in diameter distributions arising from differences in growing conditions. The use of quantiles effectively standardizes absolute size distributions, simplifying both intra-specific and inter-specific comparison of trends in relative density shifts along species-size distributions.

We assessed geographic variation in species relative density shifts at two scales: ecoregion divisions and subsections⁶⁹. Ecoregion divisions (shown for our study area in Fig. 1) are large geographic units that represent broad-scale patterns in precipitation and temperature across continents. Ecoregion subsections are subclasses of ecoregion divisions, differentiated by variation in climate, vegetation, terrain, and soils at much finer spatial scales than those represented by divisions. We again used the simple moving average estimator to estimate the mean and variance of the FSI by species within each areal unit (i.e., drawing from FIA plots within each areal unit to estimate mean and variance of the FSI). As a direct measure of changes in relative tree density, spatial variation in the FSI is indicative of spatial shifts in species distributions during the remeasurement interval (i.e., range expansion/contraction and/or within-range relative density shifts). That is, the distribution of populations shift toward regions increasing in relative density and away from regions decreasing in relative density during the temporal frame of sampling. We map estimates of the FSI for each areal unit to assess spatial patterns of changes in relative density and identify regions where widespread geographic shifts in species distributions may be underway.

We sought to quantify the average effect of forest disturbance processes on changes in the relative density of top tree species in the western US over the interval 2001–2018. To this end, we developed a hierarchical Bayesian model to determine the average severity and annual probability of disturbances (i.e., wildfire, insect outbreak, and disease) on sites where each species occurs. Average severity was modeled as

where y_jk is the FSI of species j on plot k, α_j is a species-specific intercept, β_jl is a species-specific coefficient corresponding to the binary variable x_lk that takes the value of 1 if disturbance l occurred within plot k measurement interval and 0 otherwise. The intercept and regression coefficients each received an uninformative normal prior distribution. The species-specific residual standard deviation ς_j received a uninformative uniform prior distribution⁷⁰.

On average, disturbance will occur at the midpoint of plot remeasurement periods, assuming temporal stationarity in disturbance probability over the study period. As plots in this study are remeasured on 10-year intervals, we assume that tree populations have, on average, 5 years to respond to any disturbance event. Hence, our definition of disturbance severity, β_jl’s, cannot be interpreted as the immediate change in relative tree density resulting from disturbance. Rather, disturbance severity (as defined here) includes the immediate effects of disturbance, as well as 5 years of change in relative tree density prior to and following disturbance (where disturbance is assumed to be functionally instantaneous).

Annual probability of disturbance l on plot k was modeled as

where Δt_k is the number of years between successive measurements of plot k, viewed here as the number of binomial “trials,” and ψ_jl is the species-specific probability for disturbance which was assigned a beta(1,1) prior distribution. Hence, annual probability of disturbance is assumed to vary by species j and by disturbance type l.

We estimate the mean effect of forest disturbance processes on changes in species-specific relative tree density by multiplying the posterior distributions of β_jl and ψ_jl. That is, we multiply species-specific disturbance severity by disturbance probability to yield an estimate of the mean change in relative density caused by disturbance over the study period. We then standardize these values across species by dividing by the average relative density of each species at the beginning of the study period. Thus, standardized values can be interpreted as the annual proportionate change in the relative tree density of each species resulting from disturbance over the period 2001–2018.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com