Mid-Scale Drivers of Variability in Dry Mixed-Conifer Forests of the Mogollon Rim, Arizona

: The structure and composition of southwestern dry mixed-conifer forests have changed significantly, decreasing forest resiliency to uncharacteristic disturbances which also threaten ecosystem services. Restoration of these forests can be informed by historical conditions; however, managers and researchers still lack a full understanding of how environmental factors influence forest conditions. We investigated historical and contemporary variability in dry mixed-conifer forests in northern Arizona and identified important environmental drivers. We utilized forest sample plots and dendrochronological reconstruction modelling to describe forest conditions in 1879 and 2014, respectively. We used correlogram analysis to compare spatial autocorrelation of average diameter, basal area and tree density, and structural equation modeling to partition the causal pathways between forest structure, forest composition, and a suite of environmental factors reflecting climate, topography, and soil. Historical (1879) reconstructed forests had significantly fewer trees, lower basal area, and higher average diameter than contemporarily (2014). Composition has shifted from ponderosa pine dominance towards a more mixed-species composition. Historically, forest structure did not exhibit strong spatial autocorrelation, but contemporary tree density and diameter were strongly autocorrelated. Environmental factors described little variation in historical forest conditions but are more important for contemporary conditions. Managers can utilize this increased understanding of variation to tailor silvicultural prescriptions to environmental templates.


Introduction
Mixed-conifer forests cover approximately 1 M hectares in the southwestern United States [1] and provide critical ecosystem services such as wildlife habitat [2], watershed protection [3], carbon sequestration, and nutrient cycling [4,5]. However, 20th-century land-use practices, specifically unregulated logging, grazing, and active fire suppression, have disrupted natural fire regimes [6][7][8]. This has decreased these forests' resilience [9,10]-the "ability of an ecosystem to regain structural and functional attributes that have suffered harm from stress or disturbance" [11]. Disrupting the historical disturbance regimes of frequent, low-severity fires has altered the structure and composition of these forests [12][13][14] and has decreased their resilience to disturbances such as severe wildfire, insects, disease, and climate change [9,10]. Restoring forest structure and composition can increase ecological resilience, and many restoration prescriptions are guided by historical reference conditions [11].
Reliance on historical reference conditions assumes that these characteristics can provide context for managing contemporary ecosystems, and species conservation can be accomplished by approximating the range of environmental conditions present prior to ecological degradation [15]. While there has been criticism that historical conditions will become irrelevant under future climatic conditions [16], restoring a more characteristic composition and structure would increase ecosystem resiliency [3,9] and allow these ecosystems to resist type conversion due to wildfire and climate warming [9,10,17].
In mixed-conifer forests of the American southwest, species composition varies along a temperature-moisture gradient. At one end of the spectrum, 'warm-dry' mixed-conifer is dominated by ponderosa pine (Pinus ponderosa) and can include other fire-tolerant/shade-intolerant species such as Gambel oak (Quercus gambelii), Douglas-fir (Pseudotsuga menziesii) and southwestern white pine (Pinus strobiformis). 'Cool-moist' mixed-conifer generally lacks ponderosa pine and has a greater composition of fire-intolerant/shade-tolerant species such as quaking aspen (Populus tremuloides), Engelmann spruce (Picea engelmannii), and white fir (Abies concolor) [18]. In the warm-dry type, species composition has shifted towards more shade-tolerant species over the last 100+ years since fire exclusion; ponderosa pine has decreased in dominance and the proportion of southwestern white pine and white fir have increased [12][13][14]19].
Forest density in warm-dry forests has also changed from pre-fire exclusion conditions. Historically, these forests had densities ranging from 51.6 to 245.6 trees ha −1 and basal areas of 7.9 to 28.5 m 2 ha −1 [9]; however, contemporary forest conditions across the southwest are much denser [12,14,[19][20][21][22][23]. These overly dense forests are at an increased risk of burning at high severity in an increasing number of large, uncharacteristic fires across the southwest [24,25]. This contemporary fire regime differs from the natural historical fire regime of frequent, low-to mixed-severity fires [9]. Mean fire return intervals ranged from 2 to 8.5 years on the Mogollon Rim in northern Arizona [13] and to 31.6 years in New Mexico [26] and 32 years in southwestern Colorado [12].
While previous studies have described the natural range of variability, managers still lack full understanding of the relationships that drove historical forest structure. Topography, climate, and soil factors are known to influence forest condition. Topographic models and measures of solar radiation have been used to understand variation in forest structure and composition [27,28]. Climate has strong interactions with fire frequency [12,26,29], and precipitation and climatic water availability drive variation in tree establishment [30][31][32] and density [23,33]. Soil properties, especially parent material, also influence forest structure and composition [23,[34][35][36].
Additionally, a more complete understanding of spatial heterogeneity in dry mixedconifer forests is needed. Variability in forests is spatially structured and hierarchically organized [9,23,37]. Fine-scale (<4 ha) spatial patterns typically describe the arrangement of individual trees within a stand and are nested within mid-scale (4-400 ha) patterns, which describe the variation among stands within a landscape. Mid-scale patterns are further nested within landscape-scale (400+ ha) structure. A review of spatial analyses in firefrequent forests found that most spatial analyses of reference conditions have focused on fine-scale tree patterns of clusters of trees, single trees, and openings that manifest at scales from 0.4 to 4 ha [9]. However, such patterns might be difficult for forest managers to implement in stand-level treatments [37]. Landscape restoration projects that solely focus on fine-scale patterns risk creating landscapes that are heterogeneous at fine scales, but homogeneous over mid-to landscape scales [38]. The lack of mid-scale heterogeneity analyses in southwestern dry mixed-conifer forest limits managers from designing appropriate restoration projects.
Our study addressed these knowledge gaps by investigating drivers of historical variability in dry mixed-conifer forests on the Mogollon Rim in northern Arizona. Specifically, we focused on answering the following research questions: (1) What were the historical structural conditions in these ecosystems? (2) How did structural conditions vary spatially, and has spatial variation changed since fire exclusion? and (3) What were the drivers of variability in these forests, and how have drivers changed in relative importance since fire exclusion? Answers to these questions can help silviculturists design restoration treatments that have appropriate levels of spatial and structural variability and better reflect environmental templates. By analyzing forest variability of historical and contemporary time periods, we can further explore how fire exclusion has impacted these forests.

Study Design
The Mogollon Rim is an escarpment stretching approximately 320 km from northern Arizona to eastern Arizona and forms the boundary between the Colorado Plateau to the north and the Southern Basin and Range Province to the south [39]. Our study area on the Mogollon Rim ranges from 2223 to 2399 m in elevation, has a mean annual temperature of 9.3 °C and a mean annual precipitation of 89 cm (Table A1). The forests fall into the warm/dry classification of mixed-conifer forest [18], with a major component of ponderosa pine, and a mixture of southwestern white pine, Douglas-fir, white fir, Gambel oak, and quaking aspen. Additional tree species found on the Mogollon Rim include New Mexico locust (Robinia neomexicana) and bigtooth maple (Acer grandidentatum) [19]. Lowseverity fires burned frequently on the Mogollon Rim with a mean fire interval of 2 to 8.5 years until 1879 [13]. The forests on the Mogollon Rim are currently the target of restoration efforts to increase ecological resilience [40] and protect important municipal water supplies [41].
This study utilized existing pre-treatment data collected from the long-term ecological assessment and restoration network (LEARN) which was established in 2014 by the Ecological Restoration Institute, in coordination with the U.S. Forest Service Mogollon Rim Ranger District of the Coconino National Forest (see Figure 1). The surveyed area covers approximately 250 ha, broken into six blocks where each block consisted of three ~16 ha treatment units with 15 0.04-ha circular, fixed-radius plots arranged on a 60-m grid within each unit for a total of 270 plots across the study area. Crews completed data collection in 2014, following methods described in detail in Roccaforte et al. [42]. Surveys recorded species, diameter at breast height (DBH: 1.37 m above the ground), diameter at stump height (DSH: 40 cm above the ground), total height, height to the base of the live crown, two crown radii measurements (long and short side) and condition (live tree, snag, log, cut stump, etc.) of all trees taller than breast height and all dead trees that may have predated Euro-American settlement (ca. 1879). Crews collected tree cores on all pre-settlement trees, trees larger than 40 cm DBH, and 10% of all trees smaller than 40 cm DBH. Map of study area. Points represent the location of selected major cities, and colors ranging from brown to blue represent low to high elevations. The Mogollon Rim is located southeast of Flagstaff. The study area on the Mogollon Rim (inset) is arranged into six blocks.

Determining Historical Variability
To determine the historical variability at our site in the dry mixed-conifer forests of the Mogollon Rim, we reconstructed historical forest structure and composition using field plot data and dendroecological reconstruction techniques. These techniques were developed and discussed in detail most recently by Rodman et al. [19] to incorporate additional species. This reconstruction model estimated the diameter of each tree during a set reconstruction year. We used 1879 for the reconstruction year, as a nearby by study [13] reported this to be the date of fire regime disruption. Historical tree diameters were based on dendrochronological data (i.e., cross-dated increment cores collected from trees on field plots-see above) when available, and species-specific "back-growth" regression equations when dendrochronological data was not available. These "back-growth" equations estimated the historical diameter using species-specific growth and bark thickness [43] equations. Historical diameters for dead trees were estimated by using current diameter and decomposition equations based on snag/log condition classes to determine death date [44], and then input into the "back-growth" equations to estimate the diameter during the reconstruction year. Cut stumps were assigned death dates of either 1965 or 1981, corresponding with observed releases in tree growth [45] and harvesting operations on the Mogollon Rim [19]. While this method of reconstructing forest structure may fail to detect some small trees that died and decomposed prior to the contemporary surveys, comparisons to historical surveys indicate that 91 to 94% of pre-settlement trees can be identified by contemporary surveys in areas lacking recent disturbance [46,47].
We used the results of the reconstruction model to summarize average diameter (DBH; cm), tree density (trees ha −1 ), and basal area (m 2 ha −1 ) for each field plot. We quantified composition by calculating the ecological importance value (EIV; Equation (1)) of each species in each plot, as described by Curtis and McIntosh [48]. EIV describes the importance of a given species by accounting for its relative density (abundance) as well as relative basal area (dominance), of the species. This index ranges from 0 to 2 and is calculated using the following equation: where EIVspp is the species-specific ecological importance value; nspp and BAspp are the species-specific tree density and basal area, respectively; and ntotal and BAtotal are the total tree density and basal area, respectively. We summarized these measures of forest structure and composition across all study plots to determine the historical variability, as well as the contemporary conditions at our Mogollon Rim site.

Measuring Spatial Variability
We evaluated spatial variability across the study area by quantifying spatial autocorrelation. We used the 'spdep' package in R [49,50] to calculate Moran's I [51,52]. We chose to use Moran's I to evaluate spatial variability because this metric is useful for describing patterns of continuous phenomena, estimating the size of patches, and is commonly used in research (e.g., [14]). Significant values of Moran's I indicate positive and negative correlations between the attributes of interest for geographically compared plots [53]. We calculated Moran's I for average diameter, tree density, and basal area and then compared the empirical values at specified distances (lags) to simulations of expected values given no significant autocorrelation to evaluate significance.
When calculated globally (i.e., across the entire study area) Moran's I describes the spatial autocorrelation of a variable-that is, whether the variable is distributed independently across a landscape or is dependent upon the value of its neighbors. When calculated locally, Moran's I can also be plotted over increasing lag distances to create a correlogram, describing the lag distance at which observations are highly correlated or exhibit spatial autocorrelation. We calculated a local Moran's I over 17 lags, where each lag generally represents the 60 m distance between survey plots. Comparing the correlograms of the historical and contemporary data, we visually evaluated whether spatial variability had changed since the exclusion of natural frequent fire regimes.
The arrangement of the experimental blocks at the study site posed a challenge to using Moran's I in our analysis. Ideally, we could evaluate a complete range of distances up to the maximum distance between any two points in the study area. However, there are significant gaps between some blocks. Blocks 2 through 5 are contiguous but Blocks 1 and 6 are disjunct from the others (see Figure 1). In Blocks 1 and 6, Moran's I cannot be calculated across distances larger than the maximum distance within the block, so we limited our analysis to 1000 m, approximately the distance across one block. While this hinders our ability to make inferences about landscape spatial patterns, we were still able to evaluate variability across mid-scales.

Identifying Drivers of Variability
We used structural equation modeling (SEM) techniques in Amos software [54] to identify the important relationships between environmental factors and forest structure and composition that drove variation at the study site. SEM is an analytical framework useful for investigating complex ecological systems because it can model multivariate relationships, explicitly partition direct and indirect causal relationships, include unmeasured concepts as latent or composite variables, and provide a test of causal inferences [55,56]. SEM has been successfully used in southwestern forests to understand ponderosa pine regeneration [32] relationships between environmental conditions, fire history, and understory species richness and abundance [57,58] and to explore the drivers of plat-ani-mal interactions [59]. We chose to use SEM in our analysis to evaluate the relative importance of abstract environmental factors, and to explicitly incorporate the relationship between structure and composition as an indirect causal relationship.
To evaluate the relationships between environmental factors, forest composition, and forest structure, prior to and following fire regime disruption, we built two independent models using reconstructed and contemporary field plot conditions. These two models followed the same model building techniques, and started from the same a priori model, but used separate historical and contemporary datasets. While we could have tried multigroup modeling, which would focus on determining which pathways differ most between time periods [60], we were more interested in precise estimates of the effects of the drivers of variability, optimized for each time period.
Our conceptual model describes our hypothesis that environmental factors directly influence forest structure and composition, and indirectly influence forest structure through composition (see Figure 2). Measured variables were grouped into three broad environmental factors: topography, climate, and soils. The measured variables of each environmental factor were all correlated to each other (not shown in Figure 2 for simplicity). Each factor could have been represented by many measures, so we compiled a large pool of potential explanatory variables to select from (see Table A1 for a complete, detailed list of explanatory variables). Conceptual model used at the beginning of structural equation modeling. This model represents our initial hypothesis that topography, climate, and soil directly drive variation in structure and composition, and indirectly drive variation in structure through composition. Dashed boxes represent conceptual groupings of variables, solid boxes represent measured variables, hexagons represent composites of multiple measured variables, and arrows indicate causal relationships in the data. Variables are described in Methods: Identifying drivers of variability and details for all variables are provided in Table A1. Correlation between all environmental variables was included in the model, but for simplicity were not shown in the diagram.
We used composite variables extensively as predictors of forest properties. Composites are an interpretational tool that can be used to compile the influences of multiple conceptually-linked variables into one composite effect. Topographic variables were based on a high-resolution (1 m × 1 m), LiDAR-derived digital elevation model (DEM), and were calculated at 10 m resolution in ArcMap 10 software [61] and in R (version 3.6.1) statistical software using the 'raster' [62] and 'SpatialEco' [63] packages. In the conceptual model we selected from Beer's aspect (or northeastness) [64], the heat load index (HLI) [65], and the solar radiation index (SRI; calculated for the years 1879 and 2014) [66] to represent 'aspect' as a composite variable; we selected from elevation, topographic slope position (TPI), and hierarchical slope position (HSP) to represent 'position' as a composite variable; and we selected from slope, roughness, and the topographic ruggedness index (TRI) [67] to represent 'texture' as a composite variable ( Figure 2).
Climate factors included seven different climate variables: precipitation, mean temperature, minimum and maximum temperatures, mean dewpoint temperature, and maximum and minimum vapor pressure deficits. These data were acquired from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) [68] for each month in 30 year periods from 1895 to 1924 (historical climate) and 1981 to 2010 (contemporary climate). PRISM data were used in ecological analyses where weather data have not been collected on site. We spatially downscaled climate variables from their native 800 m × 800 m resolution to the plot level resolution (60 m × 60 m) using gradient and inverse distancesquared weighting methods as described by Rodman et al. [23]. The downscaled climate variables were summarized as annual averages and seasonal (e.g., summer = 1 June-31 August) averages for the 30-year periods and assigned to each corresponding plot. In the conceptual model we selected from precipitation, dewpoint temperature, minimum and maximum vapor pressure deficit variables into a 'water availability' composite variable; and selected from mean, minimum and maximum temperature variables into a 'temperature' composite variable ( Figure 2).
Soil parent material is known to influence forest conditions [23,35]; however, parent material does not vary substantially across the study area. Data from the Natural Resources Conservation Service's Soil Survey Geographic database (SSURGO) indicated that the entire study area was Kaibab Limestone residuum [69]. In lieu of parent material, we used soil characteristics that do vary across the study area and are still important drivers of forest conditions [32,59]. Soil variables were calculated from the SoilGrids 100 m dataset [70]. Soil variables included six soil properties-percent organic C, total N, bulk density, pH, percent sand, and percent clay-at seven standard soil depths-0, 5, 15, 30, 60, 100, and 200 cm. We selected from all these variables to represent a 'soil' composite variable ( Figure 2).
We used average diameter and tree density from the historical reconstruction and the contemporary survey as indicators of forest structure in the SEMs. To represent forest composition in the SEMs, we used distance-based ordination techniques to reduce the complexity of the community data so that they would be easier to use in variable selection and structural equation modeling. We used the 'vegan' package [71] in R for this analysis. We transformed and standardized plot-level species count data before calculating Bray-Curtis distance to describe the differences between plot overstory communities. We then used nonmetric multi-dimensional scaling (NMDS) to calculate an independent three-axis solution for both the historical and contemporary community data and to calculate species scores within ordination space. Each three-axis solution was rotated to place the maximum variation along the first axis, which we then used to summarize the community composition of each plot (see Appendix B for a more complete description of this approach). Transformed average diameter, transformed density and the first axis scores from the community NMDS, make up the response variables for both our historical and contemporary models.
With more than 100 potential explanatory variables to include in each model, but desiring to keep the models parsimonious, we selected explanatory variables from each category to build each composite variable based on correlation with the response variables. We selected variables with the strongest correlations (i.e., greater than 0.1 or less than −0.1), making sure not to select redundant variables (e.g., the same soil characteristic at multiple depths, or the same climate variable at multiple seasons).
Using these criteria, we identified eleven measures of topography, climate, and soil to include in both the historical and contemporary models. For the topography factor, we selected hierarchical slope position to represent position, the solar radiation index to represent aspect, and the topographic ruggedness index to represent texture. We selected percent organic C at 30 cm, pH at 0 cm, and percent clay at 30 cm to represent the composite soil factor. For the composite climate factor, we selected winter minimum vapor pressure deficit and spring precipitation to represent water, and fall mean temperature, and winter maximum and minimum temperatures to represent temperature. For the contemporary model, we identified nine measures of topography, climate, and soil to include in the model. For topography, the same factors (HSP, SRI, and TRI) were selected. We selected pH at 0 cm, and percent clay at 30 cm to represent the soil factor. For climate, we selected spring precipitation and summer mean dewpoint temperature to represent water, and winter minimum and maximum temperatures to represent temperature.
After variables were selected, we verified the linearity of relationships before making necessary modifications to the model. To simplify the model where possible, we replaced all of the topography composite variables (aspect, position, and texture), which had only one explanatory variable, with direct effects. We found multicollinearity between the explanatory variables of 'temperature' and 'water' which may have caused path coefficient inflation in the model. This issue was resolved by combining these two factors into a single 'climate' composite. While our approach of correlating all environmental variables was conservative, this saturated the model, and model fit statistics could not be calculated with any spare degrees of freedom. To evaluate model fit we removed the least significant correlation from the model, and calculated three model fit statistics (chi-squared goodnessof-fit, Chi 2 ; an adjusted goodness-of-fit index (AGFI); and root mean square error of approximation, RMSEA) over one degree of freedom. Because there were no significant correlations in the paths removed, the final models behaved essentially the same as a saturated model. While saturated models are traditionally not used as the final SEM, this was appropriate in our analysis because our objective was to evaluate the relative importance of the environmental factors rather than test a novel theory of hypothesized relationships in the ecosystem.

Historical and Contemporary Conditions
Historical forest conditions prior to fire regime disruption (1879) differed significantly from contemporary conditions in 2014 ( Figure 3). Mean tree diameter across the study sites averaged 27.5 (95% CI: 13.3-57.0) cm historically and varied widely across sample plots. Contemporary average tree diameter was significantly lower (p < 0.0001) than historical conditions and averaged 20.1 (7.4-39.0) cm.
In 1879, basal area averaged 12.6 (2.0-79.4) m 2 ha −1 (see Figures 3 and 4). Historical basal area was highly variable but significantly lower than contemporary basal area (p < 0.0001). Contemporary basal area averaged 30.8 (12.0-58.3) m 2 ha −1 (see Figures 3 and 4). Historical forests had an average density of 165 (48-352) trees ha −1 . Comparisons with contemporary forest were found to be significantly denser (p < 0.0001), with an average density of 657 (188-2302) trees ha −1 . Historically, forest overstory composition was typical of dry mixed-conifer in the southwest ( Figure 5). Ponderosa pine accounted for approximately half of total EIV (0.934/2.0; see Equation (1)), with minor components of white fir (EIV: 0.368), Douglas-fir (EIV: 0.302), and Gambel oak (EIV: 0.214). Contemporary composition is also characteristic of dry mixed conifer, though ponderosa pine was found to have decreased in importance (EIV: 0.622) and there has been a shift towards more mesic species. For example, white fir now accounts for approximately one-third of total EIV (0.620/2.0), Douglas-fir increased to 0.370 EIV, and southwestern white pine increased from 0.018 in 1979 to 0.102 EIV in 2014. There were also changes in the relative dominance of sprouting deciduous trees: Gambel oak (0.214 to 0.106) and aspen (0.118 to 0.008) decreased; bigtooth maple and New Mexico locust increased (0.042 to 0.146; and 0.004 to 0.024, respectively).

Spatial Variability
Overall, we found little indication of spatial autocorrelation in historical forest structure, suggesting homogeneity of conditions. This homogeneity can be seen in the maps of historical forest structure (e.g., Figure 4). Historically, average diameter had low Moran's I and fell within the envelope of no significant spatial autocorrelation at all distances (Figure 6a). Basal area was significantly autocorrelated at distances of 90 and 210 m but at all other distances was within the envelope of no significant spatial autocorrelation ( Figure  6b). Interestingly, density was significantly negatively autocorrelated at 810 m, but was otherwise not significantly autocorrelated at all other lags. (Figure 6c). Figure 6. Correlograms of historical and contemporary forest structure. Each panel displays the spatial autocorrelation (as measured by Moran's I) for historical and contemporary conditions (differentiated by color). Points connected by solid lines indicate the Moran's I at a given lag distance, and the dotted lines indicate the upper and lower limits of no significant spatial autocorrelation using 95% confidence envelopes. Points that are above the threshold are significantly spatially autocorrelated, points that are below the threshold are significantly negatively autocorrelated. Historically, measures of forest structure were generally not autocorrelated. Contemporarily, average diameter and density are both significantly autocorrelated.
Unlike historical conditions, some contemporary conditions exhibited significant spatial autocorrelation. Average diameter was highly autocorrelated over distances of 0 to 360 m and significantly autocorrelated at distances up to 1000 m ( Figure 6a). Tree density was significantly autocorrelated in contemporary forests at distances up to 360 m and showed slight autocorrelation again at approximately 600 m, but was otherwise insignificant ( Figure 6c). Unlike average diameter and density, contemporary basal area remained largely uncorrelated, only showing slight autocorrelation at distances of 210 and 570 m (Figure 6b).

Drivers of Variability
In both the historical and contemporary datasets, environmental variables exhibited stronger correlations with composition as compared to structure or composition and historical correlation coefficients were generally lower than contemporary values ( Figure 7). Additionally, the relationships between average diameter and most environmental variables switched signs, changing from weakly positive to moderately negative, or weakly negative to moderately positive between historical and contemporary time periods (Figure 7). The Chi 2 , RMSEA, and AGFI evaluations of model fit suggest that the historical SEM adequately fit our data (see Table 1). In these tests the p-value indicates the probability of good fit; commonly, p > 0.05 is considered adequate, though this is only convention. The AGFI has no p-value; commonly, an AGFI > 0.95 is interpreted to represent a good-fitting model, but again, this is only convention. The historical model has good descriptive power for composition (r 2 = 0.483) but has low descriptive power for average diameter (r 2 = 0.088) and basal area (r 2 = 0.101). Topography and climate factors had relatively high importance in driving historical forest structure and composition, while soil factors had a lower impact on forest conditions (see Figure 8). The climate-composition relationship had the highest single path coefficient in the historical model (total absolute value of the path coefficient: 0.66), while topography has a higher overall impact on historical composition than climate, from the cumulative effect of aspect, position, and texture (0.73). Climate was a significant driver of historical density (0.28) and diameter (0.37) but was of approximately equal importance with topography for these structural variables (0.24 and 0.32, respectively). Aspect and texture were not significant drivers of tree density or diameter. Soil had an impact on historical density (0.35) and diameter (0.38) that was similar to that of climate and topography. Soil had a relatively small, though still significant, impact in driving historical forest composition (0.39). Composition did not significantly drive variation in historical tree density but did have an impact on average tree diameter. While statistically significant, this path coefficient (−0.23) was the smallest driver of historical average diameter. See Appendix A Table A2 for full details of the resulting model path components.  Table A2: Historic and contemporary model pathway details. Topography and climate have the greatest influence on forest composition, and climate and soil have the greatest influence on forest structure. For our contemporary model, Chi 2 , RMSEA, and AGFI results also suggest that the contemporary model converged at a solution with good fit ( Table 1). The model had good descriptive power for density (r 2 = 0.296) and diameter (r 2 = 0.317), and even better description of forest composition (r 2 = 0.632). In the contemporary model (Figure 9), climate factors had the strongest relative importance for composition, while topography had the highest relative importance for forest structure. The pathway from climate to composition had the strongest total absolute value of a path coefficient in the model (0.79) while pathways that lead from climate to density (0.21) and average diameter (0.21) had relatively low importance. Topography had the strongest influence on contemporary forest density (−0.36) and average diameter (0.32). Like the historical model, neither aspect nor texture were significant drivers of contemporary forest structure. The cumulative effects of position and aspect on composition (−0.59) was a relatively important pathway. However, texture was not an important driver of contemporary composition. Soil had a relatively weak importance to both contemporary forest density (0.20) and diameter (0.26); this differs from the historical model, where soil was a relatively equal driver of forest structure. The pathway from soil to composition (0.35) was the weakest driver of contemporary composition. The relationship between contemporary forest composition and structure differed from the one suggested by the historical model. Contemporarily, composition was a significant driver of forest density (0.26), but not forest diameter. See Appendix A Table A2 for full details of the resulting model path components. Figure 9. Contemporary structural equation model. The relative importance of each environmental factor (topography, climate, and soil) can be interpreted by the width of the pathways leading to structure and composition, which have been scaled to indicate the magnitude of the path coefficient. Paths with negative coefficients are marked with diagonal stripes. Only paths with coefficients significantly different from 0 are included in the diagram. Coefficients are also reported on the paths, and letters on the path correspond to entries in Table A2: Historic and contemporary model pathway details. Topography and climate have the greatest influence on forest composition, and climate and soil had the greatest influence on forest structure.

Historical Range of Variability
Our results suggest that the abrupt disruption to the historical frequent fire regime drastically altered the forests on the Mogollon Rim, as has been reported in dry mixedconifer forests across the southwest (e.g., [12,[19][20][21][22][23]). Reconstructed plots showed low density, with large and variable tree sizes and low, yet variable, basal area. The historical variability we found was similar to the ranges reported by both Reynolds et al. [9] and Wassermann et al. [72] for dry mixed-conifer forests in the southwest. Additionally, the average basal area and tree density we found were slightly higher than those reconstructed elsewhere on the Mogollon Rim [14,19] and in the San Juan mountains of southwestern Colorado [12]. Historical dry mixed-conifer forest at the Grand Canyon were re-ported to be denser than those we found in our study [20,39] and those reported by Williams and Baker [73] for the Mogollon Rim. When compared to mixed-conifer reference conditions outside the southwest, our results are similar to those found on the Front Range of northern Colorado and southern Wyoming [74,75], and considerably less dense than the results reported from some parts of the Sierra Nevada in California [76]. Historical conditions from other parts of the Sierra Nevada [33,77,78] and mountains in Oregon [79][80][81] had higher basal area and lower tree density than in our study area, suggesting that those forests had fewer and larger trees than those found on the Mogollon Rim.
Contemporary forest conditions diverged significantly from historical in all three measures of forest structure. Average diameter was significantly lower in contemporary conditions, reflecting the legacy of harvesting larger trees for timber and the influx of many small trees due to fire exclusion originating in the mid-1900s [14]. Basal area and tree density increased significantly. These changes are consistent with other studies in dry mixed-conifer forests across the southwest; Fulé et al. [20], Heinlein et al. [22], Rodman et al. [19,23], and Strahan et al. [14] all recorded large increases in tree density and basal area in dry mixed-conifer forests over a similar time frame.
We found a shift in forest composition away from dry mixed-conifer, towards a more wet mixed-conifer composition. For example, ponderosa pine has decreased in importance (EIV), while white fir has increased greatly. Southwestern white pine and Douglas-fir have also showed modest increases in EIV. This trend has been previously reported for forests of the Mogollon Rim [13,19] and in other mixed-conifer forests in the Southwest [14,20,22,26]. Our results parallel those of other studies concluding that frequent fires kept mesic, fire-intolerant species in check, but when released from fire they established and survived in areas where they were previously excluded [9,13].

Changes to Spatial Patterns
Changes to forest structure and composition were accompanied by changes to the spatial pattern to forest conditions, as indicated by our correlogram analysis. Prior to fire regime disruption, the forests on the Mogollon Rim exhibited little significant spatial autocorrelation across scales up to 1000 m (~314 ha). Similarly, Strahan et al. [14] found no significant autocorrelation at similar resolutions, and at distances up to 2500 m in mixedconifer forests near our site. Both random and aggregated fine-scale patterns have been found on the Mogollon Rim [19,23], as well as outside the southwest [76,82]. Random arrangement of forest structural conditions suggests multi-scalar-level structural heterogeneity that has been theorized [9,37], but not well documented.
In our study, contemporary average diameter and tree density were strongly autocorrelated at distances up to 360 m, suggesting that these forests were composed of large homogeneous patches, approximately 40 ha in size. This distance likely reflects the scale of variation for environmental patterns related to geomorphic features (e.g., ridges and valleys) on the study landscape. Autocorrelation at distances up to 1000 m indicated low variation between patches. Interestingly, basal area showed no autocorrelation despite both average diameter and tree density being autocorrelated. Disruption to the natural frequent fire regime, which maintained historical forest patterns, is the likely explanation for this change. Strahan et al. [14] described a similar shift in the mid-scale variability of community traits, and reported significant spatial autocorrelation up to 250 m. Managers seeking to restore historical spatial patterns should break up large and homogeneous stands into smaller, variable patches, and restoration prescriptions should aim to generate random structural variation within landscape patches. Silvicultural treatments such as group selection or variable density thinning may achieve these goals.

Drivers of Variability
Our model of historical drivers of variability had poor explanatory power for describing variation in forest structure, and moderate power for composition. This weak relationship between environment and forest structure was also seen in low bivariate correlations between average diameter, density and almost all environmental predictors. Low variation in historical density may have been due to the natural fire regime of frequent, surface fires. Indeed, fire may have been the main driver of historical variability and diminished importance of environmental factors in determining tree size and stand basal area. However, we were unable to include site-specific measures of fire in our model, and the relative importance of fire versus environmental factors in determining historical structure and composition needs further research.
In our model, climate and topography were the main drivers of historical variation in forest structure and composition. Path coefficients indicated that warm and dry winters, high spring precipitation, and warm fall temperatures correlated with low ponderosa pine dominance. Microsites with higher water availability had higher tree densities, and sites with cold, dry winters were associated with larger trees. Winter temperatures influence snowpack duration into the spring and may in turn determine forest density [33]. Spring precipitation is important for ponderosa pine seedling establishment in some forests [31], while fall temperatures may influence fuel moisture and wildfire behavior. Over longer periods, climate can affect fuel availability and timing of widespread fire years in mixed-conifer forests [12,26,29].
Topography was an important driver in our model of historical variation. Slope position was consistently important, reflecting the effects of microsite variability on forest conditions [27,83]. Aspect also drove variation. Sites with more solar radiation dry faster, and drought-tolerant species such as ponderosa pine and Gambel oak are better able to occupy these sites. Drying could also affect the availability of fuels for combustion and in turn affect fire frequency. Interestingly, Rodman et al. [23] did not find topographic factors to be important in determining historical forest structure at sites near ours, nor did Abella and Covington [34] at other sites in northern Arizona. The systematic sample grid in our study area may have captured a wider range of topographic conditions and greater variation in structure and composition than these previous studies.
To evaluate how the relative importance of drivers of variability may have changed, we compared historical and contemporary models. A major difference between the two models was that they describe forests with and without fire. In addition, other disturbances were present during the contemporary period including livestock grazing and selective logging. We were not able to include these processes in our contemporary model. In addition, we did not have measures of natural disturbances such as insect or disease outbreaks, droughts, or wildfires in either model. Such mortality factors would contribute to variation in forest density and composition. Without plot-level measures of these disturbances, we assumed they affected all plots equally, and were confounded with other changes that occurred between 1879 and 2014, such as fire regime disruption and climate change.
Soil was a relatively unimportant driver of variability in both historical and contemporary models. However, Rodman et al. [23], Abella and Denton [35], Kimsey et al. [36], and Laughlin et al. [58] demonstrated that soil parent material drove significant variation in forest structure and composition at northern Arizona sites. The limited variation in soil parent material across our study site likely explains our findings.
Interestingly, environmental factors exhibited a stronger influence on contemporary forest conditions than on historical structure and composition, possibly due to anthropogenic exclusion of frequent fire. The contemporary model had stronger explanatory power than the historical model, which reflected stronger correlations between contemporary environmental variables and forest measures. Climate was distinctly the most important driver of forest composition on the contemporary landscape, and indicated that sites with high spring precipitation, cold winter low temperatures, and dry summers were associated with low ponderosa pine dominance. Forest composition may have historically been constrained by species adaptedness to fire, but in the absence of fire climate asserted greater control of forest composition. Similar responses to fire exclusion have been de-scribed in other studies [9,13,14,19,20]. We found that climate increased in relative importance contemporarily, and different components of climate became important. Similarly, Mueller et al. [84] described a strengthening fire-climate relationship, and others have reported changes in the timing of widespread fire years relative to periods of drought or above average precipitation since fire regime disruption [85,86].
Topography was less important to contemporary forest composition than in the historical period, perhaps because microsite-fuel relationships [83] were less important. Topography was the strongest driver of both tree density and diameter, but this relationship changed dramatically since fire exclusion. While lower slopes and valleys historically had lower density and larger trees, contemporarily they were dominated by high densities of smaller trees. Historical logging targeted large, commercially valuable trees, which may have been concentrated on lower, more productive microsites. While density has increased significantly across the entire study area, Rodman et al. [23] and Stephens et al. [33] both found that density increases were greatest on mesic sites such as valleys and lower slopes. The increase in density on the Mogollon Rim sites was also related to compositional changes. Composition is now a significant driver of forest density, and intermediate conifers and sprouting species that would have previously been kept in low densities by frequent fire have increased on lower sites.

Conclusions
Overall, our study showed historically variable, heterogeneous, and open dry mixedconifer forests on the Mogollon Rim prior to anthropogenic fire exclusion, with ponderosa pine dominating stands on ridgetops, and more mixed composition persisting in drainages. These changes parallel those seen in dry mixed-conifer forests across the southwest and demonstrate the need for extensive restoration efforts to reduce forest density and increase diversity of structural conditions. Overly dense forests are at a heightened risk of large, severe wildfire, which could cause type change and reduce provision of ecosystem services. Managers restoring dry mixed-conifer forests on the Mogollon Rim can use historical conditions as guides for silvicultural prescriptions, while restoration of dry mixedconifer elsewhere in the southwest may be better served by the more general HRV described in reviews (e.g., [9]).
Based on our findings, we hypothesize that fire was a major driver of historical forest variation, and in the absence of fire, environmental factors assumed greater control of variation. This supposition points to the reintroduction of fire as a restoration tool. Reintroduction of fire would help reduce the density of small diameter, fire-intolerant trees and restore historical forest composition. Small-diameter trees have higher fire-related mortality than large trees, and persistence of fire-intolerant species would have been limited [87]. Low severity fire, even multiple burns at low severity, do not always result in forest density that approximates historical conditions, and it may be desirable to allow some fire to burn at moderate severity on the Mogollon Rim [88,89]. However, recent research suggests that this can also lead to mortality of larger trees [90]. Some loss of large trees may be acceptable in dry mixed-conifer forests when these species are undesirable from a restoration perspective; however, further research and applied management experiments are needed to better understand how mixed-severity fire can be used to achieve multiple restoration objectives.
Contemporary forest conditions are driven by environmental factors, especially those related to climate, but these relationships are stronger and different from those of the past. Historical climate relationships may not serve as an effective guide under novel climate situations, so managers seeking to increase ecological resilience to climate change are advised to use adaptive silvicultural approaches, such as implementing a wide variety of treatments and opportunistically making use of microsite variation. Topographic position may serve as a useful guide for restoring historical forest composition, indicating that forest managers should encourage pine-dominated stands on ridgetops and upper slopes, while allowing a more mixed composition of conifers and shade-tolerant hardwoods in valley bottoms and lower slopes. The topographically complex Mogollon Rim has a diversity of microsites that could serve as refugia for mesic mixed-conifer species as climate change intensifies over the coming decades. The increased understanding of historical ranges, patterns, and environmental drivers of variability presented in our study will be useful for maintaining dry mixed-conifer forests across the Southwest. Acknowledgments: Thank you to C. Ester and the SRP team for making this project possible. Thank you to the staff and field crews at the Ecological Restoration Institute for funding this research, collecting field data, and providing the dataset that is at heart of this research, as well as administrative support.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
Appendix A Table A1. Summary of all environmental variables considered for inclusion in models. Variables are organized by environmental factor, and subgroup if applicable. Summary statistics (mean, minimum, maximum, and standard deviation, brief description, units, and data source) are provided for each variable. Historical and contemporary climate variables are listed separately.  Table A2. Summary of pathways in historic and contemporary models. Letters in the Diagram Key column correspond to the letters on pathways in Figures 8 and 9. Values of NA indicate that this pathway is not included in these model diagrams. Pathway: From and Pathway: To describe the directional relationship between model components that each pathway represents. The Model column indicates whether the following values correspond to the pathway in the historical model or the Contemporary model. The Coefficient column contains the path coefficient, which describes the relative magnitude and direction of each pathway's relationship. Values with * are statistically significant from 0 (p < 0.05). The Components column contains the predictors used to calculate each pathway; if more than one predictor was used, the path coefficients used to calculate the composite are given in parentheses.

Appendix B
To represent forest composition in our models, we used distance-based ordination techniques to reduce the complexity of the community data so that they would be easier to use in variable selection and structural equation modeling. As noted in our methods, we transformed plot-level species count data using Wisconsin double standardization before calculating Bray-Curtis distance to describe the differences between plot overstory communities. We then used nonmetric multi-dimensional scaling to calculate an independent three-axis solution for both the historical and contemporary community data thus resulting in a single species score within our ordination space. Finally, each threeaxis solution was rotated to place the maximum variation along the first axis, which was used to summarize the community composition of each plot ( Figure A1). Figure A1. Ordinations of historical and contemporary forest community data. Points represent the composition of each site; the distance between points represents the similarity between the composition of each site; points that are close together are similar; points that are far apart are dissimilar. The three-axis ordination space is displayed as Axis 1 by Axis 2 (a,b), and Axis 1 by Axis 3 (c,d). Four letter species codes indicate the species scores within the ordination space: ABCO (white fir; Abies concolor), ACGR (bigtooth maple; Acer grandidentatum), PIPO (ponderosa pine; Pinus ponderosa), PIST (southwestern white pine; Pinus strobiformis), POTR (quaking aspen; Populus tremuloides), PSME (Douglas-fir; Pseudotsuga menziesii), QUGA (Gambel oak; Quercus gambelii), and RONE (New Mexico locust; Robinia neomexicana).
For our ordinations of historical and contemporary forest community data, our approach successfully described species composition in three axes, with final stresses of 0.11 and 0.12 (respectively). Each ordination was rotated to orient the most variation along the first axis (Axis 1). Axis 1 explained approximately half of variation in overstory composition (r 2 of 0.425 and 0.511 for historic and contemporary, respectively). In both ordinations, Axis 1 was used to describe a continuum ranging from ponderosa pine dominated sites at the far negative end, to sites dominated by relatively rare sprouting species (bigtooth maple and New Mexico locust) at the positive end, and intermediate values capturing sites with Douglas-fir, white fir, southwestern white pine, Gambel Oak, and aspen ( Figure  A1a,b). In the historical ordination, Gambel oak and aspen were grouped at one end of Axis 2, and Douglas-fir, white fir and southwestern white pine were grouped at the other end. Gambel oak and Douglas-fir were grouped at one end of Axis 3, and southwestern white pine, white fir, and aspen were grouped at the other end ( Figure A1c). In the contemporary ordination, there was less differentiation of species along Axis 2 or Axis 3. Gambel oak is differentiated on Axis 2 ( Figure A1b), and southwestern white pine is differentiated on Axis 3 ( Figure A1d), with all other species clustered together.