Multi-scale spatial controls of understory vegetation in Douglas-fir–western hemlock forests of western Oregon, USA

Forest understory vegetation is influenced by broad-scale variation in climate, intermediatescale variation in topography, disturbance and neighborhood interactions. However, little is known about how these multi-scale controls interact to influence observed spatial patterns. We examined relationships between the aggregated cover of understory plant species (%CU) and multi-scale controls using a largescale experiment including treatments of low (LS), moderate (MS) and variable (VS) disturbance severity replicated in second-growth Douglas-fir (Pseudotsuga meziesii )–western hemlock (Tsuga heterophylla) forests spanning climate and topographic gradients. We developed hierarchical models using a multi-step selection process, assessing changes residual spatial autocorrelation associated with progressively broader spatial scales of influence and interaction. To examine the role of plant traits in mediating multi-scale controls, we contrasted effects for early(%CES) and late-seral (%CLS) species cover. At neighborhood scales, decreases in %CU with overstory density were accelerated with increases in the relative importance of hemlock in the overstory in the in all but the LS treatment. At intermediate scales, %CU was lower in areas with higher potential radiation in undisturbed control treatments but that trend was reversed in the harvested, disturbed areas. When separated, effects of multi-scale controls differed between %CES and %CLS. Rates of increases in %CES with reductions in density increased with disturbance severity and decreased with increases in %CLS. At broader scales, %CES increased with climatic moisture deficit where potential radiation was high, and %CLS low. Similarly to %CU, %CLS was related to a threeway interaction between overstory density, disturbance and hemlock abundance. %CLS declined with increases in climatic moisture deficit where overstory density and the relative abundance of hemlock was high, and decreased with local increases in %CES. Multi-scale controls explained a portion of the observed spatial autocorrelation for %CES but not %CLS, suggesting the spatial patterning of %CLS is related primarily to unmeasured processes. Results show how understory responses to overstory density and disturbance severity vary across the landscape with moisture and potential radiation, at fine scales with neighborhood structure, and with species traits. Hence, understory responses to climate change likely depend on overstory composition and structure, disturbance and species traits.


INTRODUCTION
Forest understory vegetation, including all shrub and herbaceous species, exhibits spatial variation in relation to multi-scale controls including site characteristics, forest structural characteristics, disturbance and the traits of individual species (Miller et al. 2002, Scheller and Mladenoff 2002, Burton et al. 2011. As most studies have focused on a single factor and/or scale, less is known about how factors interact across scales to influence the spatial patterning of understory vegetation. Yet, because most plant species are limited to the understory, which functions to provide a range of ecosystem services (e.g., Gilliam 2007, Neill and Puettmann 2013, Sabatini et al. 2014, understanding the processes underlying spatial patterns carries important implications for forest management and conservation (Levin 2000, Rietkerk et al. 2004, Peters et al. 2007, Puettmann et al. 2009).
There is general agreement that factors operating over multiple spatial scales are important in controlling understory vegetation. Broad-scale spatial variation in regional weather patterns can influence understory vegetation directly by affecting physiological processes and populations, as well as indirectly by affecting overstory vegetation Dyrness 1988, Ohmann andSpies 1998). For instance, increases in climatic moisture deficit (i.e., the difference between potential evapotranspiration and precipitation) can decrease overstory productivity (Grier and Running 1977). This lower productivity in turn can lead to higher resource availability for the understory thereby increasing the productivity of understory vegetation (Campbell and Norman 1998, Reich et al. 2012, Yelenik et al. 2013. At intermediate spatial scales, variation in topography affects potential direct incident radiation (PDIR) in mountainous regions with resulting consequences for understory plant communities (Whittaker 1960, Hemstrom and Logan 1986, McCune 2007. Over finer spatial scales, disturbance can affect understory vegetation by reducing overstory densities and by directly reducing the cover of established vegetation and exposing mineral soil (Roberts 2007). At local neighborhood scales (e.g., 0-100s of square meters; Puettmann and D'Amato 2002), plants influence environmental conditions and interact directly via competition, facilitation, and allelopathy. The relative abundances of different overstory tree species can influence understory vegetation when there are strong differences in species effects on the local environment and resources (Canham et al. 1994, Berger and Puettmann 2000, Burton et al. 2011, as well as the density of tree seedlings and saplings (Frelich et al. 1998, Keeton and Franklin 2005, Dodson et al. 2014. Additionally, canopy gaps and decreases in overstory density can increase resource availability locally (Canham et al. 1990, Gray et al. 2002, and thereby increase the productivity of understory vegetation (Dyer et al. 2010). Furthermore, with a relatively taller stature, tree saplings (i.e., trees !1.37 m tall and ,5 cm diameter at breast height; dbh) can usurp resources and suppress herbaceous and shrub species (Givnish 1995). Reductions in the abundance of early-seral vegetation and associated organisms are often attributed to competition from saplings within high-density forest plantings Berger 2006, Swanson et al. 2011, but see Halpern and Lutz 2013). Similarly, spatial patterns of ground-layer vegetation tend to be more fine-grained in forests with high densities of saplings in the understory (Miller et al. 2002).
Within a forest type, local variation in PDIR, disturbance severity and neighborhood structure is often presumed to be relatively more important for controlling understory vegetation than broad-scale variation in climate. However, the role of broad-scale variation in climate may depend on local conditions; similarly, the role of local variation may vary with climatic conditions. For instance, the overstory canopy may moderate the effects of climate on understory vegetation by buffering climate extremes (De Frenne et al. 2013) or by exacerbating them (Smith et al. 2008, Neill andPuettmann 2013). On the other hand, understory vegetation at xeric sites may be more sensitive, or responsive, to intermediate to finescale variation in PDIR and neighborhood structure than at mesic sites (Whittaker 1960, Hemstrom and Logan 1986, Franklin and Dyrness 1988, Ohmann and Spies 1998. Furthermore, understory responses to multi-scale controls may be contingent upon ground disturbance, which can expose mineral soil and provide opportunities for the colonization of v www.esajournals.org new species (Fahey and Puettmann 2007, Roberts 2007, Neill and Puettmann 2013. Therefore, the roles of factors that vary at one scale may depend strongly on the local and/or broader context. Plant species exhibit variation in numerous traits associated with adaptations to disturbance and stress (Givnish 1988, Tilman 1990, Grime 2001, Reich et al. 2003. Early-seral species are characterized by high rates of colonization, fast growth rates, short lifespans, tall stature and an intolerance of stressful, low-resource conditions (e.g., shade, drought;Halpern 1989). In contrast, late-seral species are generally short in stature with long pre-reproductive periods, clonal reproduction, short-distance dispersal, slow rates of growth and tolerate low resource conditions in the understory (Givnish 1982, Lezberg et al. 1999, Nelson et al. 2007). Thus, early-seral species may be more responsive to variation in resources and disturbance, and more competitive in highresource environments than late-seral species, which may persist, or thrive, in low-resource conditions.
The objectives of this study were to (1) evaluate the roles of factors that vary over multiple scales and their interactions in controlling understory vegetation, (2) examine the relationship between multi-scale controls and the observed pattern of spatial autocorrelation, and (3) assess the role of plant traits in mediating relationships, contrasting early-and late-seral species as an example. We hypothesized that disturbance severity and overstory composition and structure mediate responses of understory vegetation to broadscale variation in climatic moisture deficit, and intermediate-scale variation in PDIR, or conversely, that the relationships between understory vegetation, disturbance severity, and overstory composition and structure vary at intermediate scales with PDIR and over broader scales with climatic moisture deficit (Fig. 1A). Secondly, we postulated that these cross-scale interactions explain similarities among neighboring plots in understory cover (i.e., patterns of spatial autocorrelation within sites; Fig. 1B). Finally, as a result of key differences in traits affecting colonization, competition and stress tolerance, we hypothesized that early-seral species would be more responsive to multi-scale environmental controls, which would explain a larger extent of the observed range of spatial autocorrelation than for late-seral species (Fig. 1A). Fig. 1. Conceptual graphic depicting the factors controlling the cover of understory vegetation over multiple spatial scales (A), and relationships with spatial patterns quantified using semivariogram parameters (B). In panel (A), maps show the differences in the spatial scales over which factors vary. We hypothesize accounting for multi-scale controls and cross-scale interactions in explanatory models of understory cover decreases the range of spatial autocorrelation in model residuals i.e., represented in an idealized semivariogram graphs in panel (B). The model intercept, or nugget parameter, is average semivariance at the smallest lag (closest distance between pairs of points), the sill is the plateau in semivariance reached at the distance at which observations are no longer spatially correlated (i.e., the range parameter).

Study area
We utilized manipulative field experiments replicated seven times across a latitudinal gradient stretching along the Coast Range and in the Western Cascades of western Oregon, USA (Fig.  2). The Mediterranean climate is characterized by maritime and topographic influences resulting in mild and wet winters contrasted with climatic moisture deficits developing during warm and dry summers (Hemstrom andLogan 1986, Franklin andDyrness 1988). Mean annual precipitation ranges from 1222 to 2559 mm  average) with little precipitation falling from June to August (Wang et al. 2012). Temperature and precipitation vary with longitude, latitude and elevation. Soils range from moderately deep to deep and well-to poorly drained Ultisols and Inceptisols. Topography is heterogeneous both within and among sites with slopes ranging from 0 to 508 (0-120%), elevations from 162 to 798 m above sea level, and aspects ranging from 0 to 3608, resulting in variation in potential direct incident radiation from broad to intermediate scales within and among sites (Appendix A).
Located within western hemlock zone of western Oregon (Franklin and Dyrness 1988), the overstory is dominated by 44-66-year-old Douglas-fir trees. Western hemlock varies within and among sites and is an important canopy species at the two sites in the western Cascades (i.e., Keel Mountain and Delph Creek; Fig. 2). Other conifer species, such as western red cedar (Thuja plicata), and hardwood species including bigleaf maple (Acer macrophyllum), red alder (Alnus rubra), Pacific dogwood (Cornus nutalli ), Pacific madrone (Arbutus menziesii ), and golden chinquapin (Chrysolepis chrysophylla) are minor components of the overstory. Forest understory communities are characterized by dominant species including sword fern (Polystichium munitum) and salal (Gaultheria shallon) that vary in relative abundance in association with Oregon oxalis (Oxalis oregonum) at mesic sites and oceanspray (Holodiscus discolor), common whipplea (Whipplea modesta) and (Hieracium albiflorum) at the driest sites (Hemstrom and Logan 1986, Franklin and Dyrness 1988, Ares et al. 2009). Aspects of known disturbance history vary among stands. Stands regenerated naturally following clearcutting or seed-tree harvests with the exception of portions of two sites (Keel Mt. and OM Hubbard), in which Douglas-fir seedlings were planted. Slash management following harvest operations varied among sites with some evidence of burning and in some cases is unknown. There is evidence of historical grazing prior to canopy closure at one site (Green Peak) and three sites were thinned precommercially. For a more detailed descriptions of the study site and management history, see Cissel et al. (2006).

Experimental design
Three disturbance treatment units (14-69 ha) and an untreated control unit (16-24 ha) were delineated at each of seven 44-66 year old stands in a randomized complete block, with the exception of two sites where a random assignment was not possible due to operational constraints (Cissel et al. 2006). The disturbance treatments were implemented based on (1) high uniform residual tree density of 300 trees per ha with undisturbed 'leave islands' (low severity), (2) moderate residual tree density of 200 trees per ha with leave islands and canopy gaps (moderate severity) and (3) variable residual densities of 300 trees per ha, 200 trees per ha and 100 trees per ha with leave islands and canopy gaps (variable severity). Live tree density in the untreated control areas (CON) ranges from 387 to 972 trees per ha. Three to eleven percent of the area was reserved from thinning operations as circular v www.esajournals.org 'leave islands' in all treatments. Circular gap openings were created in 3-10% and 8-17% of the area of the treatment units in the moderate severity (MS) and variable severity (VS) treatments, respectively. Canopy gaps and leave islands varied among 0.1, 0.2 and 0.4 ha in area, were well spaced throughout treatment units without bias with respect to overstory or understory conditions to accommodate harvesting operations. Disturbance treatments were applied in 1997 at three sites, in 1998 at one site and in 2000 at the remaining three sites. Nearterm findings and detailed protocols and have been reported elsewhere (Cissel et al. 2006, Fahey and Puettmann 2007, Ares et al. 2009, Burton et al. 2013, Neill and Puettmann 2013.

Field sampling
Vegetation was sampled in 75-82 permanent 0.1-ha circular overstory plots within treatment units at each site during the sixth year after thinning. Fourteen plots were randomly installed in the control areas and 21 plots in the disturbance treatments, with the stipulation that plots were .15 m from another plot or from treatment boundaries. To ensure sampling of the full range of environmental conditions, additional plots were located in gap centers. Slope and aspect were measured for each plot in the field using a clinometer and hand compass. Within each overstory plot, four 0.002-ha circular understory vegetation subplots were installed and centered at 9 m in each cardinal direction from plot centers, resulting in ;300 subplots distributed across each site. Geographic coordinates for subplots were calculated based on distance from plot center. Potential direct incident solar radiation was estimated for each plot from relationships between simulations of solar trajectories and latitude, slope and aspect using non-parametric multiplicative regression (McCune 2007). Following the omission of plots for which geographic coordinates were incorrect and a single outlying subplot, the data set consisted of 2175 observations at the subplot scale, 544 plots, 28 treatment units, and 7 sites.
In each overstory plot, all trees !5 cm in diameter at 1.37 m aboveground (dbh) were tagged, identified by species, and measured for dbh using a diameter tape. For each tree, the distance and azimuth from plot center was recorded. Cover of individual vascular plant species ,6 m tall in the understory was estimated in all subplots using cover classes: 1%, 5%, 10% and then in 10% increments to 100%, total cover of all species for a given subplot can therefore exceed 100%.
We calculated total basal area of all live trees !5 cm dbh, and the basal area of live western hemlock trees in 908 wedges extending from plot centers to the plot boundary around each subplot. The percent cover of all understory taxa was summed separately for early-and late-seral species and for all species combined. Classification of seral status (i.e., early vs. late) was based on whether a species achieves peak abundance in disturbed or mature undisturbed forest, for early-and late-seral species, respectively (Halpern 1989, Burton et al. 2013).

Statistical analysis
To examine relationships among the aggregated cover of all understory vascular plant species (i.e., herbs and shrubs; %C U ), fine-scale spatial variation in overstory composition and structure, and broad-and intermediate-scale variation in climate and potential direct incident radiation (PDIR) we used mixed effect hierarchical linear models (Littell et al. 1996) in a multi-step modeling process described below (Table 1).
Step 1: Neighborhood-scale effects.-In step 1, we developed 13 alternative general linear mixed models representing alternative hypotheses for the role of fine-scale effects in overstory density (basal area, m 2 /ha), species composition (%Hemlock) and sapling density (stems/ha) and plausible two-way interactions thereof (e.g., Appendix B). Based on our scientific understanding, all potential two way interactions among overstory density, overstory composition and sapling density are theoretically possible. Although considering all potential two-way interactions can be undesirable in circumstances where a narrower range of alternative hypotheses exist (Burnham and Anderson 2002), we did not rule them out in our list of alternative models because we considered all effects ecologically plausible. To limit the number of candidate models, only models including basal area (BA) were considered (Burton et al. 2013, Neill andPuettmann 2013), and two-way interactions were only included in a model if the main effects were v www.esajournals.org also included.
Step 2: Disturbance severity.-In step 2, the bestsupported model selected in step 1 was augmented to include indicator variables to represent effects of disturbance treatments. Because disturbance effects can interact with overstory density, overstory composition and saplings to affect understory vegetation (Fahey andPuettmann 2007, Neill andPuettmann 2013), alternative models in step 2 included plausible two-and three-way interactions of those variables with the treatment variables (Appendix B).
Step 3: Intermediate-scale variation in potential incident solar radiation.-In step 3, effects of potential direct incident solar radiation (PDIR) were added to the model selected in step 2. Similarly to step 2, effects of PDIR may depend upon variables selected in the previous steps. Thus, in addition to the main effect, the alternative models included plausible scenarios with two-way interactions with other continuous variables, and three-way interactions that include categorical treatment effects (Appendix B).
Step 4: Broad-scale climate effects.-The effects of climate were added to the best supported model from step 3. To understand the integrated effects of temperature and precipitation on understory vegetation, we used annual climatic moisture deficit (CMD) calculated as the sum of the monthly differences between potential evapotranspiration and precipitation (Stephenson 1998, Wang et al. 2012. Seasonal and annual climate averages  were obtained for each site based on geographic coordinates and elevation from down-scaled spatial interpolations of monthly data, accounting for effects of local topography, rain shadows, coastal influences and temperature inversions using ClimateWNA (Wang et al. 2012). Similar to step 3, alternative models considered main effects and plausible two-way interactions between continuous variables and three-way interactions with disturbance (i.e., categorical treatment variable).
Random effects.-All models included random effects for site, treatment units within sites, and plots within treatment units. Additionally, understory vegetation exhibits spatial autocorrelation associated with disturbance history, overstory composition, soil properties, and species traits (Miller et al. 2002, Scheller and Mladenoff 2002, Burton et al. 2011). Thus, for each set of fixed effects, we compared models with and without spatial covariance. The spatial covariance model accounted for a spatially autocorrelated error structure at the subplot scale with a range, sill and nugget parameter (Fig. 1B). Exponential models were used for the results reported here because initial exploratory analyses indicated they provided superior fit (i.e., resulting in lower minimized weighted sums of squares) compared to spherical models. Models were fitted using the mixed and variogram procedures in SAS version 9.3 (SAS Institute 2013). To account for non-constant variance in the response and non-linear relationships between explanatory variables, all cover values were transformed using a logarithm transformation after adding one to percent cover values (i.e., for zeros). Model performance was assessed using AICc, a bias corrected version of Akaike's information criteria (AIC) for small samples (Burnham and Anderson 2002). Table 1. Multi-step modelling modeling process used to select the best model of fixed effects, that vary and interact across scales in relation to understory cover (%C U ), the cover of early-seral species (%C ES ) and the cover of late-seral species (%C LS ).
Step Process Variables Scale of measurement Comparisons among scales.-To assess whether the addition of explanatory variables that vary over increasingly broader scales led to improvements in model performance, we compared the top models selected in each of the steps (1-4) to each other, and to an intercept only model. Models were compared using AICc following the same procedure applied to select models within steps. To assess the contribution of each additional scale to improving model performance we calculated AIC weights (w) and weight ratios (i.e., the ratio of the weight from the best model with the lowest AIC to the weight of the model under consideration). Larger weights indicate more support for a model, whereas larger weight ratios indicate less support (Burnham and Anderson 2002). In addition to changes in AICc, the extent to which weights increase, and weight ratios decrease, with additional scales can provide information about the contribution of that scale to improving model performance.
Assessing the variability explained by mixed models is not straightforward, and most methods suffer from known disadvantages and/or a lack of validation. However, we also computed a conditional and marginal ''pseudo'' R 2 values (R 2 c and R 2 m ) to summarize the role of fixed effects and the entire model (fixed plus random effects), respectively, in explaining observed variation by dividing the variation in the predicted values by the sum of all variance components (Edwards et al. 2008, Nakagawa andSchielzeth 2013).
Analysis of residual spatial autocorrelation.-The pattern of residual variation can also provide information related to the scales at which unexplained variation occurs (e.g., McIntire and Fajardo 2009). To examine how cross-scale interactions accounted for by fixed effects explain the observed extent of spatial autocorrelation of percent cover within sites, we summarized changes in the range parameter associated with increasing variables from an intercept-only model with random effects only to models of fixed effects selected in steps 1-4. Confidence intervals (95%) were calculated to determine whether changes were statistically significant.
Species traits.-Finally, we developed and compared separate models for early-and lateseral species cover (%C ES and %C LS , respectively) following the four-step modeling process out-lined for total cover (%C U ) above to assess whether the broad traits associated with these categories (e.g., differences in colonization abilities and competition in high-resource environments versus competition/persistence in stressful environments) mediate responses of understory cover to multi-scale controls. In contrast to the models for %C U , when analyzing %C LS and %C ES , we included the cover of %C ES and %C LS , respectively, as explanatory variables in the model for fine scale effects to account for potential interactions among early-and late-seral species (e.g., competition, facilitation, interference/priority effects). Alternative models accounting for these effects and plausible interactions were assessed at the end of step 1, in step 1b before step 2. The model selected in step 1b was then carried to steps 2, 3 and 4. We also analyzed the effects of model parameters in steps 1-4 on residual spatial autocorrelation for %C ES and %C LS following the procedure outlined for %C U above.

Total cover
Effects of cross-scale interactions.-Total cover (%C U ) was related primarily to cross-scale interactions among intermediate-scale variation in potential direct incident radiation (PDIR), disturbance and local neighborhood-scale variation in live overstory composition and structure. Because residuals from all models were spatially correlated, model comparisons were based on models including spatial covariance (i.e., for all alternatives, models with spatial covariance outperformed models that did not account for spatial correlation of residuals appreciably; Appendix B). There was strong evidence for twoand three-way interactions among factors (Table 2, Appendix B). First, effects of BA varied with disturbance severity (i.e., experimental treatments) and with overstory composition (i.e., %Hemlock). In the untreated control, moderateseverity (MS) and variable-severity (VS) treatments, increases in %C U associated with reductions in BA decreased with increases in %Hemlock. However, %C U was higher and increased at a greater rate with reductions in BA in MS and VS treatments than in the control and low-severity (LS) disturbance treatment (Fig.  3). Additionally, relationships between %C U and PDIR depended on disturbance (Fig. 4). In controls, %C U decreased as PDIR increased while in disturbance treatments, %C U increased with PDIR (Fig. 4). Although the model excluding both PDIR and climatic moisture deficit (CMD) had the lowest AICc, there was no strong evidence against the model selected in step 3 with an effect of PDIR (DAIC ¼ 1.9; Table 2). Fixed effects in the final model (selected in step 3) explained a low proportion of the total variability (R 2 m ¼ 0.18 compared to R 2 c ¼ 0.31; Table 3). Spatial autocorrelation.-Similarities in %C U among neighboring subplots resulted in spatial autocorrelation. In the intercept only model, the range of spatial correlation for residuals was 113 m (95% CI ¼ 99-126). Model variables (i.e., fixed effects) added in step 2 (i.e., disturbance effects and interactions) led to reductions in the range of spatial autocorrelation (i.e., non-overlapping confidence intervals with the range parameter for residuals from a null, intercept-only model). However the range of residual spatial autocorrelation increased again following the addition of effects of PDIR and did not differ from the null, intercept only model (Fig. 10). Table 2. Final models by step for total cover (%C U ), early-seral cover (%C ES ), and late-seral cover (%C LS ). DAIC is the difference in AIC between the model and the model with the lowest AIC among steps 1-4. Support for alternatives is indicated by low DAIC (zero ¼ best), higher weights (w) and lower weight ratios. All alternative models assessed within each step, comparisons of DAIC, weights and weight ratios are given in Appendix A: Table A1.

Early-and late-seral species
Effects of cross-scale interactions.-Similarly to total cover, models for early-and late-seral cover (%C ES and %C LS , respectively) accounting for residual spatial covariance outperformed those that did not. Therefore all comparisons included an exponential model of spatially correlated residuals. For both %C ES and %C LS , the model selected in step 4 outperformed all other models, but not appreciably so for %C LS (Table 2). Additionally, the roles of specific scales, variables and interactions differed between early-and late-seral species. At fine, neighborhood scales, increases in %C ES associated with reductions in overstory basal area (BA) depended on the disturbance. Whereas increases in %C ES with reductions in BA were modest in CON, increases in %C ES were pronounced (i.e., six-to sevenfold compared to controls at the lowest levels of BA) in MS and VS treatments and intermediate where disturbance severity was low (Fig. 5A). Earlyseral cover also decreased, albeit relatively modestly, with increases in %Hemlock (Fig. 5B). Basal area, CMD and sapling density interacted with the cover of late-seral species to affect %C ES . Early-seral cover decreased with BA at a greater rate at high relative to low cover of late-seral species (Fig. 6A). Additionally, %C ES increased with broad-scale variation in CMD where %C LS was low, but did not vary with CMD where %C LS was high (Fig. 6B). Early-seral cover also increased with increases in sapling density, particularly where %C LS was also high (Fig.  6C). Finally, %C ES decreased with CMD at low PDIR, but increased with CMD at high PDIR ( Fig. 6D). Random effects explained only a small proportion of variability above fixed effects in the final model for Table 3). Similarly to total cover, %C LS increased with decreases in BA and %Hemlock. However this effect was less pronounced in the LS thinning treatment (Fig. 7). Late-seral cover was also related to two way interactions of CMD with BA and %Hemlock, with %C LS declining with CMD only where BA and %Hemlock were relatively high (Fig. 8). Finally, %C LS declined Fig. 4. Interactive relationships of total cover of understory species (%C U ) with potential direct incident radiation PDIR and disturbance treatment: control ¼ CON, LS ¼ low severity, MS ¼ moderate severity, and VS ¼ variable severity disturbance treatment. Lines show model estimates; random effects are not shown. v www.esajournals.org with increases in early-seral cover, while effects of disturbance resulted in separate intercepts showing %C LS in moderate-severity (MS) treatment . CON . LS . VS (Fig. 9). However, these fixed effects explained a very small proportion of variability (R 2 m ¼ 0.04) relative to random effects (R 2 c ¼ 0.22; Table 3). Spatial autocorrelation.-Residuals from the intercept-only model were autocorrelated to a range of 123 m (95% CI ¼ 115-132 m) for earlyseral cover and to 134 m (95% CI ¼ 108-160) for late-seral cover. For early-seral cover, the range of residual spatial autocorrelation decreased after accounting for fine-scale, neighborhood interactions, and disturbance ( Fig. 10). Accounting for variation in cover with these effects (Table 2) reduced the estimated range parameter of spatial autocorrelation to 45 m (95% CI ¼ 36-54 m) in step 4 ( Fig. 10). In contrast, although there was a decreasing trend, confidence intervals overlapped across all steps for late-seral species (%C LS ), suggesting that spatial variation in fixed effects did not account for the observed range of spatial autocorrelation.

Total cover
The total cover of understory species reflects understory biomass and productivity, which can function as wildlife habitat, and to moderate nutrient cycling and tree seedling recruitment  Table 3. Pseudo-R 2 measures for mixed models selected at each step/scale (boldface indicates values for final, best performing models). Conditional and marginal R 2 values (R 2 c and R 2 m , respectively) measure the proportion of total variance explained by the entire model (i.e., fixed plus random effects) and the fixed effects alone, respectively.

Neighborhood interactions
Disturbance severity PDIR Climate Step 1 Step 1b Step 2 Step 3 Step 4 v www.esajournals.org processes in forests (e.g., Gilliam 2007, Neill and Puettmann 2013, Dodson et al. 2014). Our results show that the relationship between the aggregated cover of understory species and overstory density varies with fine-scale neighborhood variability in overstory composition, intermediate-scale variation in potential direct incident radiation and disturbance severity. In western Oregon, dense forests comprised primarily of Douglas-fir regenerated following logging with the relative abundance of western hemlock varying locally within stands and across the landscape (Keeton and Franklin 2005). Reduc-tions in overstory density associated with canopy gaps lead to increases in resources available for understory plants (e.g., light, moisture and nutrients), in addition to lower minimum and higher maximum soil and air temperatures (Barg and Edmonds 1999, Gray et al. 2002, Schatz et al. 2012. Meanwhile, with a relatively shorter stature and horizontal foliage orientation across a larger surface and volume than Douglas-fir, western hemlock casts deeper shade and lowers temperatures locally (Thomas and Winner 2000). Additionally, relatively lower leaf nitrogen and higher leaf lignin content of western hemlock   v www.esajournals.org has been observed in our previous work (Fahey andPuettmann 2007, Neill andPuettmann 2013) and was associated with increases in the richness and cover of early-seral species (Burton et al. 2013, Neill andPuettmann 2013). Our results here are consistent with previous results, showing greater responses to density reductions in disturbed than control treatments. Furthermore, these relationships are influenced by other, intermediate-and neighborhood-scale factors. For example, the more regular spatial pattern of the harvesting disturbance may be responsible for the small influence of overstory species composition (i.e., %Hemlock) on the response of total understory cover to overstory density in low-severity (high residual live tree density) stands. The amount of harvesting disturbance also mediated plant cover responses to PDIR as influenced by intermediate-scale variation in topography. The relationship between understory cover and PDIR was likely influenced by a denser overstory taking advantage of higher light availability thus shading out understory plants.
Where disturbances lowered the overstory density or exposed mineral soil (Roberts 2007), understory vegetation cover directly reflected the PDIR trends as influenced by intermediate scale factors, such as topography.

Early-versus late-seral species
The final models selected for early-and lateseral species show how species traits can moderate relationships between percent cover and environmental effects that interact across spatial scales. At local neighborhood scales, although overstory density interacted with disturbance to affect both early-and late-seral cover, these effects also interacted with overstory composition to affect only late-seral species. Cover increased with reductions in overstory density to a greater extent in moderate-(MS) and variable-severity (VS) treatments compared to low-severity (LS) treatment and controls (CON) for early-but not late-seral species. Furthermore, late-seral cover was less responsive to decreases in density where the relative abundance of hemlock was low in CON, MS and VS. In contrast, negative effects of hemlock do not interact with density or disturbance for earlyseral species. Local increases in the relative abundance of hemlock trees in the overstory may therefore favor the persistence late-seral species over early-seral species at low overstory densities.
The negative relation of early-seral cover to late-seral cover at low overstory density and the negative relation of late-seral cover to early-seral cover suggest neighborhood-scale interactions between early-and late-seral species are largely competitive in nature, although the mechanisms may differ. For example, late-seral species may compete via clonal spread of relatively shadetolerant late-seral species while early-seral species tend to colonize via seed and allocate Fig. 10. Range parameters associated with residual spatial autocorrelation for selected models from steps 1-4 for total cover of all species (%C U ), early-seral species (%C ES ) and late-seral species (%C LS ). Bars show means with 95% confidence intervals. Reference line shows estimate (solid) and 95% confidence intervals (dashed) from an intercept only model with no fixed effects. Refer to Table 1 for a summary of the step-wise model building process and Table 2 for model variables.
v www.esajournals.org proportionally more resources towards height growth (Givnish 1982, Tappeiner et al. 1991, Tappeiner and Zasada 1993, Givnish 1995, Lezberg et al. 1999, Nelson et al. 2007). Hence, the establishment of early-seral species at low overstory densities may depend in part on the cover of late-seral species, and the extent to which disturbance reduces the cover of late-seral species. Our results indicate that early-seral species clearly benefit from increases in disturbance severity (Fig. 5A), while late-seral species cover decreases at high disturbance severity (Fig.  9). However, it is possible that such declines were partly compensated for by increased growth of residual plants (Halpern 1989, Puettmann andBerger 2006).
The cover of early-seral species was also related to increases in sapling density, particularly where late-seral cover was also high. The correspondence of early-and late-seral species with high sapling densities may result from the omission of an important environmental variable in the model that jointly affects early-seral cover, late-seral cover and sapling density (McIntire and Fajardo 2009). For example, local variation in unmeasured below-ground resources, e.g., local nutrient pools, can be influenced by past tree species cover in the region (Burton et al. 2011, Cross andPerakis 2011). In highly productive microsites, more structural components of the understory (early-and late-seral herbs and shrubs) may be abundant until the saplings grow large enough to close the canopy. Early-seral, late-seral and tree species may co-exist by partitioning resources across the soil profile (e.g., saplings and late-seral species may use deeper soil horizons compared to early-seral species; Lezberg et al. 1999). At the moment, however, there is no evidence to suggest that natural regeneration processes in variable-density Douglas-fir-western hemlock stands result in sapling densities that are sufficiently high to preclude the development of early-seral vegetation locally (Puettmann and Berger 2006, Swanson et al. 2011, Dodson et al. 2014. Complex interactions of maritime and continental air masses with topography result in climatic variation within and between the Coast Range and western Cascades (Waring et al. 1975, Daly et al. 2002. Climate effects on vegetation are generally pronounced on stressful xeric sites relative to mesic sites Dyrness 1988, Ohmann andSpies 1998). Our results show that at broad to intermediate scales, late-seral cover declines with increasing moisture deficits at high overstory densities and where the relative abundance of hemlock is high. Likewise, reductions in density and hemlock abundance have the greatest effect on late-seral cover where moisture deficits are high, extending this idea to stressful low resource conditions in dense forest understories (e.g., Smith et al. 2008, Neill andPuettmann 2013). In contrast, early-seral cover increases with CMD primarily where late-seral cover is low and PDIR is high. Increasing summer moisture deficits are known to decrease overstory leaf area index (LAI; Grier and Running 1977). Decreased overstory LAI can be associated with increased resources (e.g., light) available to understory plants during seasons when moisture deficits are not the prevailing limitation to growth, thereby explaining the greater increases in early-seral cover under greater moisture deficits (e.g., Campbell and Norman 1998, Reich et al. 2012, Yelenik et al. 2013. Differences between early-and late-seral species may be related to differences in physiological thresholds in resources for germination, photosynthesis, growth and survival. Increases in radiation in the understory (e.g., where CMD is high) are more important for early-than lateseral species, the latter having lower physiological compensation points (e.g., Givnish 1988, Neufeld andYoung 2003). These results suggest that the consequences of increased moisture limitations for early-seral cover vary as a function of other environmental conditions (e.g., PDIR), and as a function of overstory composition and structure for late-seral species.
Moisture limitation is not only determined by rainfall and temperature patterns (i.e., Climatic Moisture Deficits as examined here) but are also influenced by soil properties such as soil depth and texture (Landsberg andWaring 1997, Stephenson 1998,). The lack of soil data, in conjunction with potential errors in climate data due to the interpolations, may partially explain the variability not accounted for in the models. In the case of two ''similar'' sites (OM Hubbard and N. Soup; Fig. 2), this variability may have resulted in a ranking of moisture conditions, which is inconsistent with other indicators of site quality (e.g., site index, the average height of dominant and co-dominant Douglas-fir trees at 50 years in an even-aged forest stand; Cissel et al. 2006). To assess our calculation methods, we also calculated CMD using the Penman-Montieth equation (Landsberg andWaring 1997, Coops et al. 2000) for these two sites. The results of these calculations matched the data derived from ClimateWNA, which uses the Hargreave's equation (Wang et al. 2012) (r 2 ¼ 0.99, data not shown).

Residual spatial autocorrelation
The large amount of unexplained variation and persistence of spatial autocorrelation in model residuals is striking. Unexplained variability within and among sites may be related to a variety of factors not included in our model including differences in age, soil properties and disturbance history. For example, variability in slash treatment among sites at the time of stand initiation and/or in association with the experimental treatments, and associated spatial patterns of slash and ground disturbance, may have effects on understory plant communities that persist through the stem exclusion phase of stand development (e.g., Halpern 1989, Nelson et al. 2007, Halpern and Lutz 2013, Harrington et al. 2013. The fact that models explained less of the variability and spatial autocorrelation in lateseral cover than early-seral cover is consistent with differences in traits between the two groups. Residual spatial autocorrelation may be caused by important spatially correlated abiotic and/or biotic factor(s) not measured or included in our models (McIntire and Fajardo 2009). For example, spatial patterns in total direct solar radiation are disjunct from canopy openings in the northern hemisphere due to the lower zenith angle and may result in high levels of transmittance under high levels of basal area (Canham et al. 1990, Schatz et al. 2012. Solar radiation may be more important for explaining residual patterns in early-seral than late-seral cover because early-seral species are typically better colonizers of disturbed sites and more responsive to increases in resource levels. Additionally, colonization processes can result in spatially aggregated unoccupied locations in otherwise suitable environments. These recruitment limitations, or 'colonization deficits', weaken relation-ships between vegetation and environmental variables (Vellend et al. 2007, Ozinga et al. 2009, Burton et al. 2011) and may be relatively more important for late-seral species, which are tolerant of stressful low-resource conditions and slower colonizers. Residual spatial autocorrelation in late-seral cover may also be related to autogenic processes leading to self-organization (Frelich et al. 1998, Rietkerk et al. 2004. For example, clonal integration among interconnected ramets can lead to a decoupling of the spatial distribution of abundance from the environment (Lezberg et al. 1999). Similarly, positive and negative feedbacks between plants and soils arising through resource use and mycorrhizal mutualisms can lead to self-organized patterns in vegetation (Mallik 2003, Ehrenfeld et al. 2005. However, given the broad range of traits and strategies within these two groups (i.e., earlyand late-seral species), we cannot attribute patterns to a single specific mechanism. Future studies accounting for continuous variation in a range of above-and belowground traits (e.g., Lezberg et al. 1999) can use similar methods however to test specific hypotheses.

Summary and conclusions
This study expands our understanding of how factors interact across scales to influence spatial patterns of understory vegetation, improving our ability to predict the consequences of environmental changes over space and time for understory structure and function. Our results show how broad-scale variability in climate, intermediate-scale variability in topography and potential direct incident solar radiation (PDIR), disturbance severity and fine-scale neighborhood variation in forest structure interact to affect the cover of understory species.
Understory responses to variation in overstory density and PDIR depend on disturbance, which may permit the establishment of early-seral species by exposing mineral soil and reducing the abundance of late-seral species locally. Understory cover and the cover of late seral species decreases with increased proportion shade tolerant hemlocks in the overstory at high overstory densities in moderate-and variable-severity treatments as well as un-treated controls, but to a lesser extent in the uniform low-severity treatment. In contrast, negative effects of hemlock on early-seral cover were independent of treatment and density. Effects of broad-scale variation in climatic moisture deficit depend on intermediatescale variation in topography and potential direct incident radiation, as well as variation in neighborhood structure (i.e., overstory composition and density, associated species in the understory). Conversely, effects of density management depend on the local climate and composition of the overstory and understory. Hence, the effects of climate change may depend not only on environmental context (e.g., topography), but may be mediated by species composition and density.
These trends appear to be moderated by species traits as the specific interactions as effects differed between early-and late-seral species. However, these multi-scale controls do not explain similarities in cover between neighboring plots (i.e., spatial autocorrelation) alone. Accounting for continuous variation in specific traits (e.g., clonality, plant-soil feedback, colonization rate) may lead to mechanistic explanations for high levels of residual spatial autocorrelation and unexplained variability, particularly among late-seral species.   Note: The cover of early-seral (%C ES ) and late-seral (%C LS ) species were analyzed as either an independent and dependent variable depending on the model. CMD ¼ climatic moisture deficit. à PDIR ¼ potential direct incident radiation.

Summary of variables analyzed by site
§ BA ¼ basal area of overstory trees ! 5 cm dbh.

APPENDIX B
Comparison of alternative models developed in steps 1-4 for total cover, early-seral cover and late-seral cover Table B1. Models developed in steps 1-4 for total cover. Information criteria for models with and without spatial covariance are shown. DAICc is given for models with spatial covariance only, since AICc was always lower for these models than models without spatial covariance. Model weights and weight ratios based DAICc show relative support for candidate models, where a higher weight and lower weight ratio indicate greater support for a candidate model.     v www.esajournals.org Table B2. Models developed in steps 1-4 for early-seral cover. Information criteria for models with and without spatial covariance are shown. DAICc is given for models with spatial covariance only, since AICc was always lower for these models than models without spatial covariance. Model weights and weight ratios based DAICc show relative support for candidate models, where a higher weight and lower weight ratio indicate greater support for a candidate model.
v www.esajournals.org Table B3. Models developed in steps 1-4 for late-seral cover. Information criteria for models with and without spatial covariance are shown. DAICc is given for models with spatial covariance only, since AICc was always lower for these models than models without spatial covariance. Model weights and weight rations based DAICc show relative support for candidate models, where a higher weight and lower weight ratio indicate greater support for a candidate model. v www.esajournals.org