Western Larch Regeneration Responds More Strongly to Site and Indirect Climate Factors Than to Direct Climate Factors

: Substantial shifts in the distribution of western larch ( Larix occidentalis Nutt.) are predicted during the coming decades in response to changing climatic conditions. However, it is unclear how the interplay between direct climate e ﬀ ects, such as warmer, drier conditions, and indirect climate e ﬀ ects, such as predicted increases in ﬁre disturbance, will impact ﬁre-adapted species such as western larch. The objectives of this study were (1) to compare the relative importance of stand, site, and indirect versus direct climatic factors in determining western larch seedling recruitment; (2) to determine whether seedling recruitment rates have changed in recent years in response to disturbance, post-ﬁre weather, and / or climate; and (3) to determine whether seedlings and mature trees are experiencing niche di ﬀ erentiation based on recent climatic shifts. We addressed these objectives using data collected from 1286 national forest inventory plots in the US states of Idaho and Montana. We used statistical models to determine the relative importance of 35 stand, site, and climatic factors for larch seedling recruitment. Our results suggest that the most important predictors of larch seedling recruitment were indicative of early-seral stand conditions, and were often associated with recent ﬁre disturbance and cutting. Despite indications of climatic niche compression, seedling recruitment rates have increased in recent decades, likely due to increased ﬁre disturbance, and were unrelated to post-ﬁre weather. Compared to sites occupied by mature trees, seedling recruitment was positively associated with cooler, drier climatic conditions, and particularly with cooler summer temperatures, but these climatic factors were generally less important than biotic stand variables such as stand age, basal area, and canopy cover. These results suggest that, for ﬁre-dependent species such as western larch, increased heat and drought stress resulting from climatic change may be o ﬀ set, at least in the near term, by an increase in early-seral stand conditions resulting from increased ﬁre disturbance, although localized range contraction may occur at warm, dry extremes.


Introduction
Recent studies suggest pronounced changes in the distribution of tree species in western North America are likely over the coming decades in response to changes in climate [1][2][3]. Western larch (Larix occidentalis Nutt.), an important tree species in forests of the northern Rocky Mountains USA that is often prioritized in management decisions [4], may be especially impacted. Bioclimate envelope models developed for western larch (Larix occidentalis Nutt.) by Rehfeldt and Jaquish [5] project that much of the future distribution of suitable climatic conditions for this species will expand to locations it does not currently occupy, and these shifts may be evident by as early as 2030. However,

Study Design and Predictor Variables
To assess western larch seedling occurrence and densities, we acquired field data collected from the US Forest Service's Forest Inventory and Analysis (FIA) plots. We defined the study area as the region within the US portion of the Rocky Mountains in the states of Montana and Idaho, where western larch has been observed on FIA plots ( Figure 1). The FIA plot grid consists of approximately one plot per 2248 ha, comprising a probabilistic sample across all land cover types, forest types, and ownership groups [50,51]. Numerous site variables and western larch seedling density were measured at each FIA plot. Each FIA plot consists of four 7.3-m radius subplots, covering roughly 1/15th ha [52]. Nested within each subplot is a 2.1-m radius microplot, wherein field crews record the number of seedlings for each species present [52].
Forests 2020, 11, x FOR PEER REVIEW 4 of 29 decadal climate data at our study sites and plotted seedling and tree presence versus climatic variables and compared climatic tolerances for seedlings versus mature trees. We hypothesized that, in response to recent changes in climate, seedlings would occupy a cooler, moister subset of sites occupied by trees.

Study Design and Predictor Variables
To assess western larch seedling occurrence and densities, we acquired field data collected from the US Forest Service's Forest Inventory and Analysis (FIA) plots. We defined the study area as the region within the US portion of the Rocky Mountains in the states of Montana and Idaho, where western larch has been observed on FIA plots ( Figure 1). The FIA plot grid consists of approximately one plot per 2248 ha, comprising a probabilistic sample across all land cover types, forest types, and ownership groups [50,51]. Numerous site variables and western larch seedling density were measured at each FIA plot. Each FIA plot consists of four 7.3-m radius subplots, covering roughly 1/15th ha [52]. Nested within each subplot is a 2.1-m radius microplot, wherein field crews record the number of seedlings for each species present [52]. Plots with seedlings present are displayed with solid triangles, while the plots with seedlings absent are displayed with circles. Shaded areas represent six ecoregion province/section designations defined in McNabb et al. [53].
The sample consisted of 1286 plots that were measured between 2003 and 2016 containing western larch as a component, which was defined as any plot with any combination of at least one western larch seedling (≥15.2 cm in length), at least one live western larch tree ≥2.5 cm diameter, or at least one dead western larch tree ≥12.7 cm diameter. On each plot, FIA delineates different conditions (i.e., stands) based on attributes such as forest type, stand-size class, and ownership group. Plots with seedlings present are displayed with solid triangles, while the plots with seedlings absent are displayed with circles. Shaded areas represent six ecoregion province/section designations defined in McNabb et al. [53].
The sample consisted of 1286 plots that were measured between 2003 and 2016 containing western larch as a component, which was defined as any plot with any combination of at least one western larch seedling (≥15.2 cm in length), at least one live western larch tree ≥2.5 cm diameter, or at least one dead western larch tree ≥12.7 cm diameter. On each plot, FIA delineates different conditions (i.e., stands) based on attributes such as forest type, stand-size class, and ownership group. Thus, plots may be comprised of a single stand or may contain multiple stands. Because of our desire to retain stand-level variables that cannot be logically aggregated at the plot level (e.g., habitat type), we decided to use data from a single stand to represent each plot containing multiple stands, as opposed to aggregating stand-level data across the entire plot, as has been used in other studies based on FIA data [26,54]. For plots containing multiple stands (179 of 1286 plots), we developed decision criteria to select a 'primary' stand for the plot that favored those with a western larch forest type, containing western larch, and occupying the largest portion of the plot area.
We computed 35 predictor variables, some of which are collinear, for our analyses, and grouped these variables into five broad descriptive categories (Table 1). Most abiotic variables (e.g., aspect) and biotic variables (e.g., live basal area) were measured or computed on FIA plots as described in USDA [52] and O'Connell et al. [51]. For each plot, we determined the number of years since fire disturbance (dating back to 1984) using either FIA data (i.e., DSTRBCD in O'Connell et al. [51]) or data from the Monitoring Trends in Burn Severity (MTBS) database [55]. Plots without fire disturbance during this time period were assigned a time-since-fire of 75 years to facilitate analysis. To reduce the number of ecoregion-province-section designations [53], we grouped all plots within individual sections of the Middle Rocky Mountain Steppe-Coniferous Forest-Alpine Meadow Province into a single designation (M332) and grouped all plots within the Palouse Dry Steppe and Intermountain Semi-Desert Provinces into a single designation (331/342), while leaving plots within the Northern Rocky Mountain Steppe-Coniferous Forest-Alpine Meadow Province in their individual sections (M333A, M333B, M333C, and M333D). We decided to use the forest type determined by field-recorded dominance based on stocking (i.e., FLDTYPCD in O'Connell et al. [51]) rather than type calculated by FIA's stocking-based forest type algorithms (i.e., FORTYPCD in O'Connell et al. [51]) because the latter assign a forest-type of "non-stocked" when stocking is less than 10% [51], which is common following disturbances such as fire that often result in pulses of western larch regeneration. When no live trees occur on a plot, FLDTYPCD is assigned based on seedlings, if any are present, or based on examination of similar, undisturbed stands adjacent to the plot footprint. To reduce the number of forest type designations, we aggregated forest types to the group level [51] and then further aggregated any groups that made up fewer than 2% of observations into a new group designated as 'other'. This resulted in seven forest type groups: Douglas-fir, fir/spruce/mountain hemlock, hemlock/Sitka spruce, lodgepole pine, ponderosa pine, western larch, and others. Folded aspect was calculated from field measurements [56]. Due to the wide area of latitude represented in our study area, we converted actual elevation to equivalent elevation by adding 129.4 m to absolute elevation for every 1-degree difference from the minimum latitude among all plots [57].
Climatic variables (Table 1) were obtained using the ClimateWNA application [58] and are based on PRISM data [59]. We extracted both baseline climate data ) and recent decadal climate data (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). The former likely reflects establishment conditions for older trees, while the latter corresponds closely with the period during which seedlings measured on our FIA plots were establishing and were used in modelling seedling presence/absence and changes in seedling density, as well as for comparing climatic niches of seedlings and trees and for identifying recent relative to baseline climatic changes. Comparisons between baseline and recent decadal climate at our study sites can be found in Appendices A and B. All measures of temperature were higher for the recent decadal period, while all measures of precipitation were lower (Appendix B). We also extracted monthly climate data for the 1984-2016 time period to assess the effects of post-fire weather on seedling recruitment for any plots that burned during this period (Objective 2).

Factors Related to Seedling Presence and Density
To identify significant predictors of western larch seedling presence, we used data from the most recent visit for each FIA plot containing western larch as a component (1286 plots) to construct three statistical models of western larch seedling presence: logistic regression (LR), classification trees (CT), and random forests (RF). We reclassified the FIA seedling count data as a binary (presence/absence) response variable for each plot, where presence of at least one western larch seedling qualified as a presence for the plot. The intent of using three different modelling approaches was to identify common factors that explain seedling presence regardless of model structure [54]. LR was included because it permits assessment of the importance of individual variables by comparing their z-scores and their associated p-values. CT and RF were included because they provide intuitive interpretation of variable importance, are robust to collinearity and interaction among predictor variables, and make no a priori assumptions about the distributions of response or predictor variables [60]. All three models were developed in R [61] with a classification threshold of 0.5.
LR was performed using function glm in baseR as a generalized linear model with a logit link and binomial family [61]. To meet the model's assumption of noncollinearity among predictor variables [62], we performed principal components analysis using the climatic variables and equivalent elevation. The first two principal components resulting from this analysis explained over 90% of the variance. The first principal component was clearly related to temperature and the second was clearly related to moisture. To reduce collinearity, we used these first two principal components in place of the climatic variables and equivalent elevation in logistic regression models for seedling presence. To further reduce collinearity among pairs of stand variables, we computed each variable's point biserial correlation coefficient with seedling presence, and eliminated the variable with the lower value from inclusion in logistic regression models.
Classification trees were built using R package rpart (Recursive Partitioning and Regression Trees; [63]). To determine the appropriate complexity for the classification tree, we identified the maximum complexity parameter (CP) with relative error of less than one standard error [64], and then re-built the classification tree using the chosen CP (CP = 0.027). We developed the RF model in R package randomForests [65] using 500 trees for evaluation (parameter ntree = 500) and considering 7 variables for each split (mtry = 7). We used 10-fold cross-validation and subsequent comparison of confusion matrices to evaluate each model's performance.
We identified important predictors of western larch seedling presence as those that were statistically significant predictors in at least one model. For the LR model, significant predictors are those whose variable coefficients were significantly different than 0 based on z scores (i.e., p (|z|) < 0.01). Significant predictors for the CT are those from the nodes of the final classification tree. For the RF model, significant predictors were identified by their corresponding decreases in overall accuracy and Gini index [60,62].
To complement modelling, chi-square tests of association (for categorical predictor variables) and Wilcoxon rank sum tests (for continuous predictor variables) were performed to test for significant differences between plots with and without seedlings for individual predictor variables (α = 0.05 for each method). The objective of the chi-square analysis was to test whether seedling presence is associated with each categorical variable, based on a null hypothesis that there is no association, where larger chi-square values (and thus smaller p-values) provide evidence for an association [66]. The non-parametric Wilcoxon rank sum analysis [66] evaluates the association between seedling presence and continuous variables by comparing the variables' distributions at sites with seedlings versus at sites with no seedlings. Thus, both of these nonparametric tests indicate whether an association exists between seedling presence and other variables. We calculated Cramer's V scores to determine the strength of association with each predictor variable, where the values of V range from zero (no association) to one (perfect association [67]. Values less than or equal to 0.3 indicate a weak association, those between 0.4 and 0.5 a medium association, and those greater than 0.5 a strong association. Differences in climatic values between plots where seedlings were present and absent were illustrated using histograms. We also assessed the density and distribution of western larch regeneration. We first calculated mean and median seedling density for all plots on which western larch seedlings were present and then performed Kruskal-Wallis (for multiple comparisons) or Wilcoxon rank sum (for single comparison) tests (α = 0.05) [66] to compare densities among levels of categorical predictor variables. The Kruskal-Wallis test corresponds to a nonparametric analog of a one-way analysis of variance, and tests for differences among distributions of continuous variables [66]. We calculated Pearson correlation coefficients to determine correlations with continuous predictor variables (α = 0.05).

Trends in Seedling Presence and Density
To assess recent trends in western larch seedling presence and density, we used data from FIA plots on which measurements have been made at two time periods. In Montana, time 1 data was collected on these plots from 2003-2006 and time 2 data from 2013-2016. In Idaho, time 1 data was collected from 2004-2006 and time 2 data from 2014-2016. Since the number, arrangement, and attributes of stands found on a plot often change over time due to disturbance and stand dynamics, we used only stand-level data from time 2 in our analysis. To determine whether changes in stand attributes that occurred between time 1 and time 2 influenced changes in seedling density, we also calculated change in live basal area of western larch and of all tree species, change in live and total trees per hectare, and change in live canopy cover between the time periods. Selection of the primary stand (described above) for plots containing multiple stands (58 plots) was performed using stand-level data from time 2, and we included only plots with stands that occupied the same footprint on the plot area at time 1 and time 2. These constraints on the initial 1286 plots yielded a total of 406 re-measured plots.
We determined seedling presence and density for each plot at both time 1 and time 2. To compare density among levels of categorical predictor variables, we calculated mean and median change in seedling density and performed Kruskal-Wallis tests (for multiple comparisons) or Wilcoxon rank sum (for single comparison) tests (α = 0.05) [66]. To determine correlations with continuous predictor variables, we calculated Pearson correlation coefficients (α = 0.05). To restrict our assessment to trends in natural regeneration, we did not include any plots that were planted with western larch. We characterized weather conditions during the first three post-fire water years (October 1 to September 30) for each plot that burned between 1984 and 2013. Unfavorable conditions could reduce seedling recruitment during this window of time that is crucial to western larch regeneration [23,28]. We examined three variables that we felt were most indicative of potential heat and drought stress: CMD, GSP, and mean temperature during the warmest month (MWMT). To account for among-site variation in climate settings, we used an approach similar to that employed by Harvey et al. [28] and standardized weather variables for the 3-year post-fire time window to ± SD of the local 1984-2016 average. To test for significant relationships between these post-fire weather variables and seedling presence and density, we used Wilcoxon rank sum tests (for seedling presence) (α = 0.05) and calculated Pearson correlation coefficients for seedling density (α = 0.05).

Climate Niche Differentiation between Seedlings and Trees
We characterized climatic niches for seedlings and mature trees (≥15.2 cm, d.b.h) for both baseline and recent decadal periods. We assumed that the baseline period (1961-1990) represents conditions adequate for survival of mature trees, and the recent decadal period (2001-2010) represents conditions conducive to regeneration. We calculated means and medians for several climate variables, defined niche boundaries using the 5th and 95th percentiles of each climate variable for plots with seedlings and plots with mature trees, and estimated niche tolerances for seedlings vs. mature trees as the difference between the 5th and 95th percentiles [26]. We then used unpaired t-tests to compare mean values for plots with seedlings versus plots with mature trees, and calculated differences between seedling and tree medians, 5th and 95th percentiles, and tolerances. Differences in climatic values between plots with seedlings versus those with mature trees were illustrated using histograms.

Factors Related to Seedling Presence and Density
Collectively, the three models suggest that biotic stand variables related to tree density, basal area, and canopy cover were the most important predictors of seedling presence (Figures 2 and 3). Seedling presence had negative relationships with live basal area (LR, RF), live basal area of western larch (CT), larger stand-size class (RF), and live canopy cover (RF). Furthermore, seedling presence was negatively related to stand age based on the CT ( Figure 3) and RF models, but was positively related based on the LR model. This contradiction can be explained by examination of a histogram of stand age distribution for plots with seedlings present (Figure 4), which shows a large number of these plots had stand ages of less than 20 years. Plot numbers are low in the 41-60 year-old category, before increasing slightly to a second peak in the 101-120 year-old category, after which plot numbers drop off steeply. Inability of the LR model to account for non-linear relationships between continuous predictor variables and seedling presence may also explain the positive relationship detected for live plus dead trees per hectare by this model. The RF model performed slightly better than the other two models, while the CT model had slightly lower performance than the LR model (Table 2).        Numerous climatic variables were also important predictors of seedling presence, but were generally less important than biotic stand variables as indicated by the CT (Figure 3) and the RF model, where climatic variables had much lower values for mean decrease in overall accuracy and Gini index relative to biotic stand variables. In general, seedling presence was negatively related to  Numerous climatic variables were also important predictors of seedling presence, but were generally less important than biotic stand variables as indicated by the CT (Figure 3) and the RF model, where climatic variables had much lower values for mean decrease in overall accuracy and Gini index relative to biotic stand variables. In general, seedling presence was negatively related to climatic variables that indicate moister and warmer site conditions (LR, CT), and positively associated with PRATIO (CT). Lastly, seedling presence was negatively associated with increasing years since fire disturbance (LR) and slope (LR), and positively associated with equivalent elevation (LR).
Our LR model detected significant interactions between the second principle component (indicating site moisture) and both live basal area and live trees per hectare. As site moisture decreased the likelihood of larch regeneration increased with greater live basal area and live trees per hectare. In addition, several significant interactions between stand and climate variables are evident from our CT model (Figure 3). Live basal area of western larch was a significant predictor of seedling presence only on plots where MWMT was ≥18 • C, while PRATIO was a significant predictor only on plots where live basal area of western larch was ≥0.11 m 2 /ha.
Western larch seedlings were detected on a total of 269 of 1286 plots (21%), with mean and median seedling densities of 1739 and 556 seedlings per hectare, respectively, on plots where seedlings were present. The likelihood of seedling presence was significantly associated with each of the 9 categorical predictor variables, although only the association with stand-size class is considered of medium strength (Table 3). Significant differences were detected between plots with and without seedlings for 23 of 26 continuous predictor variables (Table 4). In general, plots with western larch seedlings tended to have (1) experienced more recent fire disturbance, (2) lower live basal area of western larch and of all tree species, (3) lower canopy cover, (4) fewer trees per hectare, (5) a younger stand age, (6) higher densities of seedlings of other tree species, (7) higher equivalent elevation and lower slope, and (8) cooler temperatures, lower amounts of precipitation, and lower reference evapotranspiration levels (Eref). Histograms of sites with and without seedlings, relative to pairs of climatic variables, confirm these patterns and show plots with seedlings clustering toward the cooler, drier portions of climatic space ( Figure 5).    Statistically significant differences in seedling density were detected among levels of forest type, ecoregion province/section, stocking class, and reserved status ( Table 3). The weak nature, although statistically significant in most cases, of correlations between seedling density and the continuous predictor variables (Table 4), suggests that seedling density was less influenced by many of these variables than seedling presence.
Plots with recent fire disturbance (109 plots total) had a much higher likelihood of seedling presence (65% versus 17% for unburned) and tended to occur on drier sites than unburned plots. MAP and WP were significantly lower on burned plots (746 and 206 mm, respectively, for median values), compared to unburned plots (836 and 269 mm, respectively, for median values). Only 7 of 109 (~6%) recently burned plots had mean annual precipitation (MAP) values greater than 1200 mm Figure 5. Histograms of western larch seedling presence/absence for all plots (left column) and for only plots with recent fire disturbance (right column) relative to selected climate variables. Statistically significant differences in seedling density were detected among levels of forest type, ecoregion province/section, stocking class, and reserved status ( Table 3). The weak nature, although statistically significant in most cases, of correlations between seedling density and the continuous predictor variables (Table 4), suggests that seedling density was less influenced by many of these variables than seedling presence.
Plots with recent fire disturbance (109 plots total) had a much higher likelihood of seedling presence (65% versus 17% for unburned) and tended to occur on drier sites than unburned plots. MAP and WP were significantly lower on burned plots (746 and 206 mm, respectively, for median values), compared to unburned plots (836 and 269 mm, respectively, for median values). Only 7 of 109 (~6%) recently burned plots had mean annual precipitation (MAP) values greater than 1200 mm and WP values greater than 400 mm, whereas 148 (~13%) and 203 (~17%) of 1177 unburned plots exceeded those precipitation levels ( Figure 5).

Trends in Seedling Presence and Density
Seedlings presence increased from 67 of 406 plots (17%) at time 1 to 90 of 406 plots (22%) at time 2. Twenty-five plots (6%) had seedlings present at time 1 and not at time 2, 48 plots (12%) had seedlings at time 2 and not at time 1, 42 plots (10%) had seedlings present at both times, and 291 plots (73%) lacked seedlings at both times. Seedling density increased on 61 plots (15%), decreased on 39 plots (10%), and did not change on 306 plots (75%). Predictor variables most strongly associated with the likelihood of an increase in seedling density or changes in seedling density from time 1 to time 2 were basically the same as to those identified for seedling presence and density presented above. None of the 13 climatic variables were significantly associated with the likelihood of an increase in seedling density or were correlated with changes in seedling density. Tables summarizing these results are provided in Appendices C and D.
Seedling presence and density on recently burned plots was unaffected by post-fire weather. We found no significant relationships between seedling presence and any post-fire weather variable, nor any significant correlations between seedling density and any post-fire weather variable.

Climate Niche Differentiation between Seedlings and Trees
Comparison of means, tolerances, and boundaries for climatic variables suggest climatic niches differed between seedlings and mature trees. Seedling means and medians were lower for all climatic variables except DD_0 (where larger values actually represent cooler climatic conditions), median CMD, and PRATIO (Table 5). Tolerances for all precipitation variables, and for CMD and Eref, were narrower for seedlings than for trees, while temperature tolerances were similar or wider for seedlings than for trees (Table 5). Trees had higher values than seedlings for the 95th percentile, and lower values than seedlings for the 5th percentile, for MAP, GSP, and CMD, suggesting seedling niche contraction at both margins for these measures of precipitation and moisture stress. Trees had higher values than seedlings for both the 95th and 5th percentiles for WP and for all for temperature variables, suggesting a shift in seedling niche toward lower WP and cooler climatic conditions, respectively (Table 5). These measures of niche location and breadth indicate contraction of the climatic niche boundary in terms of precipitation and a shift in niche location toward cooler, drier climatic space. Similarly, the spread of plots with trees extends further along the horizontal axes of precipitation and CMD (Figure 6), further illustrating contraction in climatic niche breadth for seedlings. Plots with seedlings also had a higher mean and median PRATIO values than plots with trees (Table 5). Temperature and precipitation from the recent decadal period shifted toward the warmer, drier portions of climatic space, and many plots experienced MWMT (not shown) and mean coldest month temperature (MCMT) values during the recent decadal period that were nearly 2 • C greater than maximum values during the baseline time period (Figure 6).

Non-Climatic and Indirect Climatic Factors Affecting Seedling Recruitment
Changing climatic conditions can impact tree species directly, by increasing heat and drought stress, and indirectly, by altering disturbance regimes [23,34,69]. Our observations that seedling recruitment tended to be most strongly related to stand variables indicative of early-seral conditions suggest that non-climatic factors were of greater importance to western larch regeneration success than direct effects of climate change. However, it is likely that indirect effects of climate change have also strongly impacted regeneration dynamics through the promotion of increased wildfire disturbance. Our observations of greater seedling recruitment in stands with lower live canopy cover and live tree basal area agree with established understanding of regeneration requirements for western larch. Western larch seedlings require the sunlit conditions found in open stands that are often created by recent disturbances such as wildfires or timber harvesting [21,41], which reduce canopy cover and the density of overstory trees. The ability of recent disturbances to create stand conditions favorable to larch regeneration is also likely reflected in our detection of a negative relationship between stand-size class and seedling presence, and of substantially younger mean stand ages on plots where seedlings were present compared to those where they were absent (48 and 88 years, respectively; Table 4).
Recent disturbance history may also explain higher larch regeneration on plots with lodgepole pine and western larch forest types. Following disturbance, these early-seral forest types occupy many sites capable of supporting western larch in the northern US Rocky Mountains [43,68]. We found that plots with lodgepole pine and western larch forest types had younger median stand ages (56 and 71 years, respectively) and were more likely to have experienced recent fire disturbance (26% and 10% of plots, respectively), than plots in other forest types (median stand age 87 years; 5% of plots with recent fire disturbance). We also detected a similar, though weaker, trend for plots with a ponderosa pine forest type. These plots had intermediate levels of seedling presence and density (Table 3), coupled with relatively young median stand ages (51 years) and a relatively high likelihood of recent fire disturbance (16% of plots). The ponderosa pine forest type can be early-seral on sites capable of supporting western larch, especially those at the warmer, drier margins [43,68].
Western larch can successfully regenerate following stand-replacing fires, and also low-intensity, understory burns, which often leave larger, fire-resistant western larch and ponderosa pine (Pinus ponderosa Lawson & C. Lawson) in the overstory [70,71]. We believe that our detection of a second peak in seedling presence for plots with stand ages between 101 and 120 years ( Figure 4) is partially due to recent understory burns in some of these mature stands within our sample. In fact, recent fire disturbance was detected on some stands in our sample with stand ages exceeding 250 years. Our sample also likely contained stands that were treated using seed-tree harvesting methods, which seek to enhance larch regeneration through scarification of the soil surface during seedbed preparation, and retention of scattered, mature larch in the overstory as a seed source [21,72]. Although FIA crews do not determine the specific cutting method when assigning treatment codes, seed-tree cuts are a common silvicultural prescription in western larch stands [72]. Because seedbed preparation activities and retention of mature larch produce conditions similar to those found in mature stands following understory burns, stands treated with seed-tree harvesting methods are also likely to show successful seedling recruitment coupled with older stand ages.
Greater likelihood of seedling presence under decreased site productivity and at higher elevations concurs with our observation that larch seedlings tended to occur on cooler, drier sites (further discussion below in Section 4.2). These cooler, drier sites are typically found at higher elevations and tend to be less productive than other sites capable of supporting western larch [43,68]. However, the strength of the relationship between seedling presence and site productivity is relatively weak (Cramer's V = 0.214), as is the correlation between seedling density and equivalent elevation (r 2 = 0.130), suggesting that climatic or stand factors may better explain seedling presence at cooler, drier sites.
Observed differences in seedling presence and density among different ecological sections and among ownership categories seem to reflect differing amounts of recent fire disturbance and/or cutting. The Flathead Valley (M333B) and Northern Rockies (M333C) sections of province M333, where we detected the highest seedling presence rates and densities (Table 3), have experienced higher rates of recent fire disturbance (12% of plots) relative to other areas (5% of plots). Likewise, lands under 'other federal' ownership, which include National Park Service, Bureau of Land Management, and Department of Defense lands, have experienced higher rates of recent fire disturbance (40% of plots) compared to other ownership categories (4%-9% of plots under National Forest, state, and private/tribal ownership). Lands under state and private/tribal ownership have experienced higher rates of recent cutting (17% and 18% of plots, respectively) compared to those under 'other federal' and National Forest ownership (2% of plots for each). Thus, National Forest lands have experienced the least amount of combined fire disturbance and cutting, which has contributed to lower rates of seedling presence (18% of plots on National Forest lands had seedlings) relative to other ownerships (34% of plots on 'other federal' lands, 22% of plots on State lands, and 27% of plots on private/tribal lands).

Climatic Factors Affecting Seedling Recruitment
Despite the influences of the stand and site factors detailed above, our observation of higher likelihood of seedling recruitment on cooler, drier sites suggests that climatic conditions-particularly seasonal climatic variables-can also limit western larch regeneration. Heat and drought stress are important causes of larch seedling mortality during warm, dry summer conditions that are common on western larch sites [16,21]. Indeed, our classification tree model indicated that seedling recruitment was negatively related to MWMT and positively related to PRATIO. Furthermore, decreases in late-summer precipitation and increases in late-summer temperatures measured across our study areas (Appendix A) suggest late-summer heat and drought-stress have risen in recent years. We suspect this negatively impacted regeneration on warmer sites, while regeneration on cooler sites, such as those occurring at higher elevations, may have been enhanced. Several studies have predicted upslope movement of fire-tolerant montane species including western larch under future climate change [28,73].
Our results highlight the importance of the interplay of climate, stand, and site factors in governing western larch regeneration dynamics. For example, our LR model suggests that the negative relationships between live basal area and/or live tree density and seedling presence found under most conditions is altered as sites become drier, with regeneration enhanced by higher live basal area and tree density under these conditions. Similarly, the results of our CT model suggest that the amount of live basal area of western larch was an important predictor of seedling presence only on sites where MWMT exceeded 18 • C, indicating regeneration success on these warmer sites was more sensitive to stand structure [26] and/or close proximity to a seed source.
Increased wildfire disturbance resulting from the indirect effects of warmer, drier conditions, and subsequent increases in conditions favorable to larch seedling establishment, must be contrasted with the potential negative impact of increased heat and drought stress experienced by seedlings during the one to five-year post-fire window crucial for successful establishment [23,28]. Comparison of seedling presence/absence across all plots and only those with recent fire disturbance (Figure 2) suggest that temperature and moisture stress may limit seedling recruitment in the warmest, driest portion of climate space, while recruitment on wet sites receiving more than 1200 mm of MAP is limited by lack of recent fire disturbance needed to create suitable stand conditions for seedling establishment. Based on our results, it is unclear whether increases in heat and drought stress alone can determine success or failure of larch regeneration at sites that do not fall along the warm, dry margins of climate space. Furthermore, the lack of a relationship between seedling recruitment and post-fire weather provides evidence that the stand conditions created by fire and other disturbances may outweigh the influence of post-fire climatic conditions even on warmer, drier sites. Although several studies have documented recent declines in post-fire conifer regeneration due to increases in post-fire heat and drought stress [15,30], the response of western larch has been mixed. Harvey et al. [28] found no relationship between post-fire drought severity and western larch seedling establishment, while Urza and Sibold [23] found seedling establishment was highly sensitive to wetter post-fire weather conditions.

Trends in Seedling Presence and Density
Despite documentation of warmer, drier climatic conditions across our sites, our results indicate that western larch seedling presence and density actually increased in recent years, leading us to reject our hypothesis that seedling recruitment would decline. In addition, our hypothesis that important predictors of trends in seedling recruitment would be those reflecting heat and moisture stress also proved incorrect. We believe, instead, that trends in seedling recruitment were mostly due to the same site and stand factors that were responsible for overall patterns in seedling recruitment (discussed in Section 4.1). Of primary importance was the presence of early-seral stand conditions resulting from recent fire disturbance and cutting. In contrast, analyses of relationships between trends in seedling recruitment and predictor variables suggest that climatic factors were relatively unimportant in explaining these trends.
As a fire-adapted species, western larch can expect to benefit, at least in terms of regeneration, from increased fire disturbance [23], although regeneration success may be limited at the warm, dry margin of climatic space. Furthermore, our data support predictions that the area occupied by early-seral forest types, such as western larch, will increase in the coming years with increasing fire disturbance [23,74]. This may become particularly evident on moister sites that have experienced relatively little fire disturbance during recent decades.

Climatic Niche Differentiation between Seedlings and Trees
Our results show that western larch seedlings occupy only a subset of climate envelope space compared to that occupied by mature western larch trees. We identify two notable shifts in the distribution of western larch seedlings: first, seedlings occur more often on cooler and drier sites where mature western larch occurs, and second, seedlings rarely occur on sites at the warmest extremes of the distribution of mature western larch ( Figure 6). These results reflect observed shifts in temperature, precipitation, and climatic moisture deficit at sites occupied by western larch (Figure 6). Specifically, MAT, MCMT, and MWMT have increased between the baseline and recent decadal periods (Appendix B), contributing to downward shifts in seedling niche location along these axes ( Figure 6). Little seedling recruitment occurred at sites where MAT exceeded 8 • C, which exceeds the baseline MAT of sites with mature western larch, despite the presence of mature trees at such sites ( Figures 5  and 6). In contrast, decreases in MAP and WP, and increases in CMD, have been less pronounced, resulting in only slight contraction or expansion (WP) of seedling niche boundaries at the dry margins of these precipitation and moisture deficit gradients. The combination of niche differentiation ( Figure 5) and climatic shifts from baseline to recent decadal observations ( Figure 6) suggest that increases in temperature have affected seedling recruitment more than decreases in precipitation or increases in moisture deficits.

Management Implications
Managers seeking to maintain or increase western larch extent or abundance have relied on previous research showing its tendency to regenerate following wildfire [10,21] as well as climate envelope models that make projections about where tree species will occur in the future under changing climatic conditions [5,49]. Many climate envelope models, including those used to develop management plans within the range of western larch, fail to consider whether seedlings occupy different climate envelopes than mature trees of the same species [49]. Our results indicate that sites within the hot and dry margin of western larch's current distribution may no longer be suitable for regeneration of western larch. Therefore, post-disturbance revegetation efforts should account for potentially contracting ranges of western larch to cooler, drier sites.
We found that although western larch regeneration is related to climatic factors, it is more strongly related to non-climatic factors such as disturbance history, stand density, and stand age. However, these non-climatic factors are also affected by increasing temperatures and decreasing precipitation, both of which were observed at our study sites. As temperatures and moisture stress increase, disturbance such as wildfire and insect epidemics are also likely to increase. Following such disturbances, efforts to encourage western larch regeneration may be more successful on cooler microsites such as northern aspects. Additionally, managers should evaluate whether opportunities exist for establishing larch at higher elevation sites that were formerly marginal for this species, keeping in mind that the water-holding capacity of soils on cooler, less productive sites may be less favorable to western larch regeneration.
The importance of climatic and non-climatic factors, and hence the value of projections made by traditional climate envelope models that fail to incorporate site and stand characteristics, depend on the spatial scale in question. This is particularly true in complex mountainous terrain, where slight variations in elevation, slope, aspect, and shading by other vegetation can strongly affect microclimatic factors such as site-specific climatic moisture deficit and potential evaporation. It is assumed that climate governs species' distributions at broad geographic scales, whereas non-climatic factors, such as edaphic characteristics and biotic interactions, determine the specific sites occupied by species at finer spatial scales [7,75]. At broad geographic scales, focus often centers on projecting geographic shifts in suitable climatic conditions, identifying populations at risk, and selecting genetic material suited to future climatic conditions [3,5]. Thus, climate envelope models provide a useful approach to understanding the dynamics that influence species distributions at these broad geographic scales [3,7,75]. At finer spatial scales, the focus turns to predicting the persistence of local populations, and to understanding the factors that explain local shifts in species occurrence, such as subtle shifts in elevation and aspect to track suitable regeneration conditions [76]. Our results primarily address these latter questions concerning that nature of species' responses at finer spatial scales in complex, mountainous terrain. Understanding the dynamics governing western larch occurrence and regeneration at finer spatial scales may be particularly important for forecasting probability and persistence at drier sites along the eastern and southern trailing edge of its range. Here, greater abiotic stress and more intense interspecific competition combine to reduce growth and establishment, while disturbance and growth efficiency-related mortality increases [8,9]. At these sites, managing for a mosaic of seral stages and monitoring regeneration under a variety of site conditions may prove the best way to mitigate the impacts of climate change on western larch over the coming decades.

Limitations of This Study
We were unable to assess the impact of several factors that may affect western larch regeneration but are not measured on FIA plots, including seed availability, burn severity, and local variation in soil and substrate characteristics. Seed availability is closely tied to distance to seed source, which has been identified as an important predictor of regeneration success in western conifers, particularly among species with heavier seeds such as ponderosa pine, Douglas-fir, grand fir, and subalpine fir [21,22,24]. Distance to seed source may be less important for regeneration of western larch [28], whose light, wind-born seeds are dispersed greater distances [21]. Nonetheless, lack of regeneration on some of our sites may reflect absence of a nearby seed source, and not the influence of climatic or stand and site factors. Several studies have noted relationships between the spatial patterns of burn severity, seed dispersal, and patterns of post-fire regeneration [18,22]. Western larch regeneration is favored by higher burn severity [23], as long as patch size is not large enough to hinder seed dispersal [28].
Local variations in soil characteristics likely also influence larch regeneration. Soils with greater moisture-holding capacity and productivity, such as deep, porous soils [21], and volcanic ash-derived soils [77], are known to promote western larch presence and productivity. Soils with greater moisture-holding capacity can moderate the amount of heat and drought-stress experienced by seedlings during late-summer [77]. This ability may play a part in our failure to detect a strong relationship between seedling recruitment and moisture/temperature except at the warm, dry margins of climatic space.
Lastly, despite its importance, regeneration comprises only part of the dynamics determining response of western larch to changing climatic conditions. Growth and mortality rates also affect larch persistence under increasing levels of disturbance and heat and drought stress. In a companion analysis of larch mortality rates on the same plots included in this study (in prep.), we found that rates of natural larch mortality on remeasured plots were lower during the 10-year period between measurements (0.76% mortality/year) than during the 5-year period prior to initial plot establishment (1.17%/year), and were lower compared to co-occurring conifers (1.20%/year for Douglas-fir, 2.73%/year for lodgepole pine, 2.44%/year for subalpine fir, and 1.08%/year for grand fir). This suggests that recent changes in climatic conditions have not resulted in increased mortality of mature larch trees.

Conclusions
Western larch regeneration currently occupies 21% of all sites with western larch present on a probabilistic sample of 1286 sites in the Northern Rockies of USA. We found evidence that moisture availability and precipitation affect seedling distribution, although non-climatic variables (tree density, basal area, canopy cover, stand-size class, and stand age) were stronger predictors of western larch seedling presence than climatic variables. Sites with seedlings present represent only a cool and dry to mesic subset of the climatic space occupied by mature larch trees, supporting our hypothesis that niche differentiation between mature tree versus seedlings of this species may be occurring. Changes in seedling density over time were more strongly related to disturbance history and stand characteristics than to climatic factors, where seedling density was most likely to increase following fire. This result coincides with an increase in fire disturbance and thus suggests that increases in heat and drought may be offset, in the near term, by increased regeneration under early-seral conditions created by disturbance, particularly at cooler sites. Analysis of climate data at western larch sites confirms that recent conditions (2001-2010) are warmer and drier than in past decades , and regeneration has shifted within climatic space as seedlings occupy only a relatively cool and dry subset of the climate space occupied by mature trees. The combination of increased fire disturbance in this region and climatic niche differentiation of seedlings suggests that western larch may retreat from hot, dry sites, including those with recent fire disturbance, and that post-disturbance reforestation efforts should focus on the cooler, drier sites within western larch's range.  Appendix B Table A1. Comparison of mean, median, and 95th and 5th percentile values of climate variables, including P-values of differences in means (α = 0.05) according to paired t-tests for baseline   Appendix C Table A2. Results of analyses of change in western larch seedling density for nine categorical predictor variables. Predictors with significant associations according to Chi-square contingency tests have Cramer's V scores in bold (α = 0.05). Levels of predictors with significant differences in seedling density or change in seedling density according to Kruskal-Wallis (for multiple comparison) or Wilcoxon rank sum tests (for single comparisons) are indicated by different letters next to their median seedling density or mean change in seedling density values. Appendix D Table A3. Summary statistics of 26 continuous predictor variables for plots with and without an increase in western larch seedling density. Predictors with significantly different values between plots with and without an increase in seedling density according to Wilcoxon rank sum tests (α = 0.05) have values in bold and italics. Significant correlation coefficients (α = 0.05) are also indicated in bold and italics. Variables are described in Table 1.