Data Package Metadata   View Summary

From the Methods section of Crouch 2023 Chapter 2 (dissertation completed in May 2023):

STUDY AREA

Our study area encompassed aspen ecosystems throughout Arizona, USA (Fig. 2.1a, 2.1b) (Little 1971; Perala 1990). In contrast to more northerly latitudes, aspen ecosystems in Arizona are a relatively rare feature on the landscape, occupying less than 2% of forested land (Johnson 1994; Rolf 2001; Gitlin et al. 2006; Zegler et al. 2012). On the southwestern edge of its range, aspen is limited to relatively high elevations, where lower temperatures and higher precipitation allow this drought-intolerant species to survive (Perala 1990; Rehfeldt et al. 2009). Aspen can be found as low as 2,000 m in elevation in the ponderosa pine (Pinus ponderosa var. scopulorum) forest type, where small pockets of aspen occur on north-facing slopes or in drainages with increased water availability (Rasmussen 1941; Covington et al. 1983; Martínez González and González-Villarreal 2005; Fairweather et al. 2008; Zegler et al. 2012). As elevation increases into the mixed-conifer and, in some areas, spruce-fir forest types, the aspen component tends to be more abundant and less aspect-limited (Rasmussen 1941; Merkle 1962; Fairweather et al. 2008; Zegler et al. 2012). In these forest types, aspen occurs not only in pure stands but also in mixed stands with conifers, including ponderosa pine and Douglas-fir (Pseudotsuga menziesii var. glauca) at lower elevations, white pine (Pinus strobiformis or Pinus flexilis var. reflexa) and white fir (Abies concolor) at mid elevations, and subalpine fir (Abies lasiocarpa var. arizonica) and Engelmann spruce (Picea engelmannii) at the highest elevations, where aspen reaches its upper limit above 3000 m.

SITE SELECTION

We sampled 220 aspen plots that represent the range of conditions under which aspen exists in Arizona (Fig. 2.1b). These plots were located across seven major areas: North Kaibab (n = 19), South Kaibab (n = 26), Flagstaff (n = 113), Mogollon Rim (n = 13), White Mountains (n = 25), Prescott (n = 17), and Coronado (n = 7) (Fig. 2.1b). All data were collected during the 2020, 2021, and 2022 growing seasons (June – October), when aspen trees had leaves. Most of our sampling occurred around Flagstaff because of the wide range of sites that aspen occupies in this area (Fig. 2.1c).

To ensure we obtained a representative sample of aspen sites and conditions, we stratified sites across four variables – elevation (≤ 2400 m, > 2400 m); aspect (north/east, south/west); ungulate management (none, fenced exclosure [2 m tall fences built around aspen stands to exclude ungulates] or jackstraw treatment [large piles of woody debris protecting aspen regeneration from ungulate browse]); and fire history (0-2 years post-fire, 2-20 years post-fire, > 20 years post-fire; included wildfire and prescribed fire) – resulting in 24 strata. We first sought to obtain one plot for each stratum, which we accomplished for 21 of the 24 strata, before building out a sample that was proportional to how much aspen occurs in each stratum. We assessed aspen’s actual occurrence in each stratum using a GIS layer of aspen’s observed range on three ranger districts surrounding Flagstaff (Flagstaff and Mogollon Rim Ranger Districts on the Coconino National Forest; Williams Ranger District on the Kaibab National Forest) (DePinte 2018). Although this layer covers only three of the nine ranger districts we sampled, it is the most accurate estimation of where aspen occurs in Arizona because it is a fine-scale layer of aspen’s recent presence based on direct observations from an aircraft (DePinte 2018). We compared the proportion of aspen observed on the landscape, based on area from the GIS layer, to the proportion of aspen plots we sampled, based on the number of plots that fell into each of our strata. We succeeded in obtaining a representative sample across elevation, aspect, and fire history, with proportions of aspen observed in each stratum versus aspen sampled differing by less than 7% for each stratum (Table 2.1). Due to a lack of accurate GIS data documenting where fenced exclosures and jackstraw treatments occur across the three ranger districts, we were not able to assess how much aspen occurs in areas treated for ungulate management. Instead, we sampled these areas evenly across strata, resulting in roughly one third of our plots occurring in ungulate management treatments (Table 2.1).

When possible, we prioritized remeasurement of existing aspen monitoring plots to reduce the number of redundant plots on the landscape and to facilitate research permission on national forest land. We revisited plots previously established by the Coconino National Forest (n = 44), the Apache-Sitgreaves National Forest (n = 5), Zegler et al. (2012) (n = 20), and Northern Arizona University’s Ecological Restoration Institute (n = 12). All four of these networks established plots using stratified or completely random sampling, ensuring the locations of these plots lacked bias. We established the remaining 139 plots by identifying aspen stands that filled target strata, standing on the edge of selected stands, laying out a transect longways through those stands, and establishing plots every 30 m along the transects. The Coconino National Forest, Apache-Sitgreaves National Forest, and Ecological Restoration Institute plots were also established along transects with plot spacings ranging from 100 m to 300 m. In contrast, Zegler et al. (2012) established sites at randomly located points within known aspen stands and sampled plots in each of the four cardinal directions 20 m from those points.

FIELD DATA COLLECTION

Each study plot consisted of two fixed-area, circular plots: an overstory plot (8 m radius) and a nested regeneration plot (4 m radius) sharing the same plot center (Zegler et al. 2012). We collected GPS coordinates at the center of each study plot, recorded whether the plot fell in an area of ungulate management (i.e., fenced exclosure or jackstraw treatment), and noted whether there was evidence of recent conifer removal, as indicated by cut conifer stumps present in or directly adjacent to the plot. For a plot to be included in our study, it had to contain at least five live aspen stems between the 8 m overstory and 4 m regeneration plots combined. In the 8 m overstory plot, all trees with diameter at breast height (dbh; height = 1.37 m) > 12.7 cm were measured. In the 4 m regeneration plot, all trees > 0.02 cm in height and < 12.7 cm dbh were measured. In the regeneration plot, we classified stems into two size classes: regeneration (< 1.37 m tall) and recruitment (> 1.37 m tall and < 12.7 cm dbh). We chose a recruitment threshold height of 1.37 m to be consistent with previous studies of aspen juveniles in Arizona (Binkley et al. 2006; Zegler et al. 2012). For all live aspen, we recorded height and dbh (except for regeneration and recruits that were < 1 cm dbh). For every dead aspen and live tree species other than aspen, we recorded size class and dbh.

For all live aspen, we documented the top three damaging agents present on each tree (Zegler et al. 2012). When more than three damaging agents were present, preference was given to agents with the greatest severity of impact (i.e., most likely to cause dieback and mortality) (Zegler et al. 2012). These damaging agents included insects, diseases, ungulate browse, other animal damage, and abiotic damages. For insects and diseases, we grouped individual species into functional groups to facilitate analysis and because some biotic damages (e.g., defoliating insects) were impossible to identify based solely on the damage they caused. These functional groups included sucking and gall-forming insects (excluding oystershell scale), bark beetles, wood-boring insects, defoliating insects, canker-causing diseases, foliar and shoot diseases, and decay diseases (USDA Forest Service 2013; Steed and Burton 2015). We assessed oystershell scale and certain cankers individually because of their potential to have outsized impacts on aspen tree health compared to native insect species and less pathogenic diseases (Hinds 1985; Zegler et al. 2012; Crouch et al. 2021, 2023). The cankers we assessed individually were Cytospora canker, Hypoxylon canker (caused by Entoleuca mammatum), Ceratocystis canker (caused by Ceratocystis spp.), and sooty bark canker (caused by Encoelia pruinosa). We lumped all abiotic damages together, which included fire scarring of stems, drought scorch on leaves, and chlorosis of leaves. We assessed animal damage to aspen stems, including browse, ungulate barking (i.e., elk chewing aspen bark), and other animal damage. In addition to directly quantifying ungulate impacts on aspen stems, we counted ungulate scat piles within the 8 m overstory plot. We identified scat piles by species (i.e., elk [Cervus canadensis], deer [Odocoileus hemionus or O. virginianus couesi], or cattle [Bos taurus]) and treated piles from the same species as distinct when piles were clearly separated, contained more than three pellets, and differed color or size (Bunnefeld et al. 2006; Rhodes and St. Clair 2018).

DATA CALCULATIONS

Using tree height and diameter data, we calculated our three response variables: density (trees ha-1) of live aspen regeneration, live aspen recruitment, and dead aspen recruitment. We did not use dead aspen regeneration density as a response variable because evidence of dead regenerating stems disappears quickly (Zegler et al. 2012). We also calculated density (trees ha-1) of live overstory aspen, dead overstory aspen, live overstory tree species other than aspen, live overstory conifers, live regeneration of tree species other than aspen, and live conifer regeneration (Table 2.2). We used height and diameter data to calculate basal area of stems > 5.1 cm dbh for live aspen, dead aspen, live tree species other than aspen, and live conifers (Table 2.2). Using the presence/absence data for all damaging agents on each live aspen stem, we calculated the proportion of stems affected by each agent in each plot (Table 2.2).

Using the GPS coordinates we collected at each plot’s center, we calculated elevation, aspect, and slope using a 30 m2 digital elevation model (Table 2.2). We transformed raw aspect into a continuous variable ranging from 0–2 with 0 representing southwest (225°) and 2 representing northeast (45°) (Beers et al. 1966). We also calculated heat load and potential annual direct radiation, two indices that assess site-level temperature based on slope, aspect, and latitude (McCune and Keon 2002). We assessed fire occurrence at each plot for the past 20 years using wildland fire perimeters from the USDA Forest Service Region 3 GIS database (https://www.fs.usda.gov/detail/r3/landmanagement/gis) and prescribed fire perimeters obtained from national forest staff. We assessed fire severity at each plot using data obtained from the Monitoring Trends in Burn Severity program (https://www.mtbs.gov/), which provides fire severity data at 30 m resolution. We created categorical variables to represent both fire occurrence and severity in addition to a binary variable for plots that burned twice in the past 20 years (Table 2.2). Finally, we used GPS coordinates and maps obtained from national forest staff to verify whether plots fell inside areas of ungulate management and conifer removal treatments, and we created binary variables for both ungulate management and conifer removal (Table 2.2).

We obtained soils data from SoilGrids (https://www.isric.org/explore/soilgrids), which provides global soil mapping at 250 m resolution (Poggio et al. 2021). We used 9 of 12 available soil metrics to capture variables that represent soil moisture (e.g., sand content and bulk density), fertility (e.g., cation exchange capacity, nitrogen, and soil organic content), rooting environment (e.g., bulk density, clay content, and coarse fragments), and chemical environment (e.g., soil pH) (Table 2.2). SoilGrids provides data up to 2 m below the surface; however, we aggregated mean values for each variable to a depth of 1 m because most lateral aspen roots occur within the first 1 m of the soil (Jones and DeByle 1985). We obtained climate data for each plot from ClimateNA (https://climatena.ca/), which downscales PRISM data (Daly et al. 2008) at 800 m resolution (Wang et al. 2016). Specifically, we obtained variables representing annual and, when available, seasonal temperature, precipitation, and drought for the five years preceding when we sampled each plot (Table 2.2). We chose five years to be consistent with other studies that have assessed the influence of climate on juvenile aspen (Clement et al. 2019; Reikowski et al. 2022). In addition to climate variables obtained directly from ClimateNA, we calculated monsoon index (summer precipitation ÷ annual precipitation) and annual dryness index (annual degree-days above 5°C ÷ annual precipitation) because of the importance of the monsoon system in Arizona and the important influence of precipitation, in general, on aspen occurrence, growth, and mortality (Rehfeldt et al. 2009; Worrall et al. 2013; Kane et al. 2014; Ireland et al. 2020).

ANALYSIS: SUSTAINABILITY OF REGENERATION AND RECRUITMENT

To determine whether aspen is sustainably regenerating and recruiting, we compared abundance of juvenile aspen to two different thresholds for self-replacement. The first set of thresholds, which we refer to as the WNA (western North America) thresholds, were 2500 stems ha-1 for regeneration and 1250 stems ha-1 for recruits as outlined in the literature (Mueggler 1989; Campbell and Bartos 2001; O’Brien et al. 2010). Because these thresholds were developed for aspen in more northerly parts of its range, we wanted to develop a second set of thresholds specific to aspen in Arizona using size class data from our study plots. These thresholds, which we refer to as the AZ (Arizona) thresholds, are site-specific and based on the overstory aspen present in each plot. We calculated these AZ thresholds based on data from 68 healthy study plots. To be considered healthy, a plot had to contain no oystershell scale and < 20% browse, which is considered the threshold of sustainable browsing (Jones et al. 2005; Rogers and Mittanck 2014). From these 68 plots, we calculated mean density of live overstory (201.9 trees ha-1), recruiting (4411.9 trees ha-1), and regenerating stems (8575.1 trees ha-1), and then we calculated the ratios between overstory trees to regenerating stems (1: 42.5) and overstory trees to recruiting stems (1: 21.9). For each study plot, we then multiplied the density of living and dead overstory aspen by both ratios. For plots with no overstory aspen, we defaulted to the WNA thresholds. We then compared observed densities of aspen regeneration and recruitment across our 220 study plots to both the WNA and AZ thresholds. To facilitate our understanding of where juvenile aspen were observed at sustainable levels, we categorized self-replacing status of regeneration and recruitment across the seven major areas where aspen occurs (Fig. 2.1b).

ANALYSIS: FACTORS INFLUENCING REGENERATION AND RECRUITMENT

We considered 69 variables that could potentially influence aspen regeneration and recruitment, representing eight overarching categories: stand structure, ungulate impacts, damaging agents, fire, management, site factors, soils, and climate (Table 2.2). We conducted two analyses – random forests and structural equation modeling (SEM) – to determine which of these factors significantly influence regeneration and recruitment. We analyzed all data in R version 4.2.1 (R Core Team 2022), using the dplyr package (Wickham et al. 2022) for data manipulation and the ggplot2 package (Wickham 2016) for figure creation. First, we used random forests to help us determine which of the 69 predictor variables had the strongest influence on our three response variables (i.e., density of live regeneration, live recruits, and dead recruits). Random forests are a useful tool for assessing variable importance in regression and classification settings among an array of potential predictors (Breiman 2001). Specifically, we used the VSURF package (Genuer et al. 2015), which used 50 random forest runs, each of which was built using 2000 trees, to rank variable importance for each of our three response variables. VSURF is robust in noisy, high dimensional settings and in the presence of highly correlated predictors (Genuer et al. 2010). VSURF outputs a ranked list of variables based on importance, which is calculated using out-of-box mean square error for each fitted tree, along with a group of variables highly related to the response that is geared towards interpretation (Genuer et al. 2010, 2015). We used both the ranked list of variables and the group of interpretation variables when building SEMs.

Once we obtained a list of the most important variables influencing each response, we used SEMs to assess how those predictor variables and their interactions influence aspen regeneration and recruitment. SEMs are an insightful tool for ecological research because they allow the user to build models based on theoretical understanding of an ecological system, resulting in a network of causal, multivariate relationships with a complete accounting of direct and indirect relationships and the relative strengths of those relationships (Grace 2006; Lefcheck 2016). SEMs are valuable in the specific context of our study because we understand how individual factors influence juvenile aspen (Crouch et al. 2023), but we do not understand how these various factors interact and which are the most important drivers of regeneration and recruitment. Our first step in building SEMs was to construct an a priori model based on our theoretical understanding of how biotic and abiotic factors influence juvenile aspen. This a priori model (Fig. 2.2) applied to all three response variables and accounted for all 69 variables that potentially influence regeneration and recruitment using the eight categories of influencing factors (i.e., climate, fire, site factors, soils, management, stand structure, ungulate impacts, and damaging agents).

For each of the three responses, we built a “full” SEM, which included the highest ranked variable based on random forests from each of the eight categories of influencing factors. We then used a combination of backward and forward selection to optimize model fit (using AIC and Fisher’s C statistic) and explanatory power (using R2 of the response variable). This optimization process included removing variables with low significance in the model and adding in more than one variable per category (e.g., adding a second climate variable) when two variables from one category had high importance values based on random forests. We also tested how swapping in one variable to replace another variable of the same category (i.e., replacing heat load with radiation) affected the model, although only one such swap resulted in both improved model fit and explanatory power (spring climate moisture index [CMI] swapped in to replace winter CMI in the live aspen regeneration SEM). We used the piecewiseSEM package to build our SEMs because this package accommodates use of mixed-effects models (Lefcheck 2016). Prior to fitting individual regressions that underlie the SEMs, we log-transformed the three response variables to satisfy normality assumptions. For the individual regressions that underlie piecewiseSEM, we used the lme4 package (Bates et al. 2015) to fit linear mixed-effects models with the hierarchical, nested structure of our plots (i.e., plots [n=220] within study sites [n=87] within minor areas [n=19] within major areas [n=7]) modeled as random effects. Study site refers to a transect or group of plots that are clustered near each other, whereas minor area refers to a group of such transects or plots in a larger but still confined area (e.g., an individual mountain or fire footprint). Because study site location was accounted for implicitly as a random effect in these mixed-effects models, we did not explicitly include major area, UTM easting, UTM northing, or other spatial variables in SEMs.

Finally, we wanted to explore specific impacts of different ungulate species (i.e., elk, deer, and cattle) on recruitment, and to do so, we fit six simple linear regression models with each of the three species’ scat counts as predictors and density of live and dead recruits as responses. Similar to the linear models that were built for SEM, these linear models were mixed-effects models fit using the lme4 package (Bates et al. 2015).