Di ﬀ erence in Regeneration Conditions in Pinus ponderosa Dominated Forests in Northern California, USA, over an 83 Year Period

: Forest inventories based on ﬁeld surveys can provide quantitative measures of regeneration such as density and stocking proportion. Understanding regeneration dynamics is a key element that supports silvicultural decision-making processes in sustainable forest management. The objectives of this study were to: (1) describe historical regeneration in ponderosa pine dominated forests by species and height class, (2) ﬁnd associations of regeneration with overstory, soil, and topography variables, (3) describe contemporary regeneration across various management treatments, and (4) compare di ﬀ erences in regeneration between historical and contemporary forests. The study area, a ponderosa pine ( Pinus ponderosae Dougl. ex P. and C. Law) dominated forest, is located within the Blacks Mountain Experimental Forest (BMEF) in northeastern California, United States, which was designated as an experimental forest in 1934. We used 1935 and 2018 ﬁeld surveyed regeneration data containing information about three species—ponderosa pine, incense-cedar ( Calocedrus decurrens (Torr.) Florin) and white ﬁr ( Abies concolor (Grod. and Glend)—and four height classes: class 1: 0–0.31 m, class 2: 0.31–0.91 m, class 3: 0.91–1.83 m, and class 4: > 1.83 m and < 8.9 cm diameter at breast height. We used stocking as proxy for regeneration density in this study. We found that historically, stocking in the BMEF was dominated by shade-intolerant ponderosa pine in height classes 2 and 3. Two variables—overstory basal area per hectare (m 2 ha − 1 ) and available water capacity at 150 cm, which is the amount of water that is available for plants up to a depth of 150 cm from the soil surface—were signiﬁcantly associated with stocking, and a beta regression model ﬁt was found to have a pseudo- R 2 of 0.49. We identiﬁed signiﬁcant di ﬀ erences in contemporary stocking among six management scenarios using a Kruskal–Wallis non-parametric one-way ANOVA. Control compartments had the highest stocking followed by burned compartments. In contemporary forest stands, recent treatments involving a combination of burning and thinning resulted in high stocking in height classes 2 and 3. Overall, the stocking in historical BMEF stands was higher than in contemporary stands and was dominated by ponderosa pine. B.N.I.E.; Resources, B.N.I.E. M.W.R.; Supervision, B.N.I.E.; Visualization, S.N. M.W.R.; Writing—original draft, S.N.; Writing—review & editing, B.N.I.E.


Introduction
Public agencies focus forest management efforts on reducing forest fuels in dry forests of California with thinning and prescribed fire treatments [1]. However, a common secondary objective focuses on promoting more resilient future forests, which are dominated by fire-resistant and shade-intolerant ponderosa pine (Pinus ponderosa Laws.) [2]. The resiliency of dry forests can be characterized in terms of in 1935 and 2018 across BMEF. Yellow pine represented a grouping of Jeffrey pine (Pinus jeffreyi Grev. & Balf) and ponderosa pine, two closely related species. To simplify, we refer to this grouping as ponderosa pine throughout the text. The forest stand conditions observed in 1935 had likely already been changed from forest conditions observed in forests prior to European settlement due to sheep grazing, which may have altered fire frequency [33]. Nevertheless, we will refer to the 1935 stand conditions as historical in the following. Thus, we are not using 'historical' to refer to conditions prior to European settlement, but rather to refer to the 1935 conditions that reflect conditions prior to active management in comparison to contemporary conditions observed in 2018. In this study, we answered the following four questions:

1.
What was the stocking of ponderosa pine, incense-cedar, and white fir for four different height classes in 1935 prior to active forest management across BMEF? 2.
Was historical stocking in BMEF associated with overstory stand structure, soil, and topography? 3.
Are there differences in stocking among contemporary (2018) forests that received several management treatments: control, prescribed fire, combination thinning, combination thinning with prescribed fire, thinning from below, and thinning from below with prescribed fire? If so, which treatments differ? 4.
How does historical stocking (1935) for the observed tree species differ from contemporary stocking (2018) for 12 compartments that received management treatments prior to 2018?
Our study allowed detailed quantitative descriptions of historical (1935) and contemporary (2018) stocking of different species under various stand conditions that were dominated by ponderosa pine. We were also able to quantify differences in contemporary stocking under several management scenarios. Information on historical regeneration by species, height classes, and environmental factors associated with regeneration is critically important to managers for prescribing forest management treatments that will determine sustainability of forest values in the future.

Study Area
This study was conducted at the Blacks Mountain Experimental Forest (BMEF) established in 1934 in northeastern California, United States ( Figure 1). BMEF (40 • 40 N, 121 • 10 W) is managed by the United States Department of Agriculture (USDA) Forest Service, Pacific Southwest Research Station, and ranges between 1700 and 2100 m in elevation [34]. BMEF is dominated by ponderosa pine (Pinus ponderosae Dougl. ex P. and C. Law), with lesser amounts of white fir (Abies concolor (Grod. And Glend) and incense-cedar (Calocedrus decurrens (Torr.) Florin) present at higher elevations [35]. At lower elevations, Jeffrey pine (Pinus jeffreyi (Grev. And Balf)) is found within the stands [36]. Western juniper (Juniperus occidentalis Hook.) is sparsely distributed as individuals or small groves in the lower elevations of the forest [36].

Management History
Due to observed high rates of insect-caused mortality [37,38], between 1937 and 1960, most compartments of the forest were subject to sanitation/salvage cutting that removed many of the decadent and high risk trees in the forest, those of risk class 3 and 4 [39]. These harvests were usually light, removing less than 20% of the standing volume. Remnants of the earlier structure remained in four research natural areas (RNAs; Table 1, Figure 1). The RNAs represent undisturbed and late successional forests [36,40]. During the 20th century, there was no fire at BMEF since 1899 (unpublished fire scar data from Carl Skinner, United States Forest Services, retired). However, after 1997 prescribed fire was applied in two of the four RNAs.
The Bull thinning project (hereafter, thinning from below) was started in 2011 and concluded in 2012. Prescribed burns within this project were performed in fall of 2017. Thinning in a Low Structural Diversity treatment, removing all but the mid-canopy trees (hereafter, combination thinning), was conducted in 1996 and prescribed burns were performed in 1997 [41]. All observations on current conditions were conducted in 2018. Treatments are summarized as follows (see details in Table 1):

Regeneration Data
Data for the analysis came from two sources: (i) 1935 historical regeneration data from 100 compartments ( Figure 1) and (ii) 2018 regeneration data-hereafter referred to as contemporary regeneration data-for 12 compartments spanning a range of current conditions at BMEF.

Historical Regeneration Data
Regeneration data were systematically sampled across the entire BMEF in 1935. BMEF was divided into 4055 strips of approximately 1 ha each. Along the length of each strip, 0.0016 ha (4-milacre) quadrats were installed at every 50 m (164 feet) interval [42]. Within each quadrat, seedling data were tallied by three species-pine, white fir, and incense-cedar. Data were also tallied by height classes of 0-0.31 m (0-1 feet), 0.31-0.91 m (1-3 feet), 0.91-1.83 m (3-6 feet), and >1.83 m (6 feet) for seedlings less than 8.9 cm (3.5 inches) diameter at breast height (dbh) [42]. The pine species in 1935 was a combination of ponderosa pine and Jeffrey pine [42]. The data were aggregated to compartment-level, resulting in seedling density by species and height class for 100 compartments [42].

Contemporary Regeneration Data
In 2018, the 1935 sampling design was imitated within 12 compartments, which spanned across six recent treatments with two replications each across the entire BMEF (Table 1, Figure 1). The 2018 field sampling was specifically designed and conducted in order to compare the contemporary data to the available 1935 data, because little quality data on regeneration were collected since 1935. While in 1935 square 0.0016 ha (4-milacre) quadrats were used, we used circular plots of radius 2.27 m (7.45 feet) split into four 0.000404 ha wedges for ease of implementation in the field. Plots were established at an interval of 50 m spacing among permanently marked grid points at 100 m square spacing. The 2018 seedling data, which were collected within the milacre quadrats, were then aggregated to the compartment level as was done with the 1935 data.

Stocking
Stocking, as defined by Gingrich [43], is a relative term that describes a given stand density based on a management objective or a desirable density. For 1935, only the seedling density per hectare summaries calculated from the 0.0016 ha (4-milacre) plots were available rather than the raw inventory data. Therefore, we were unable to observe stocking, defined as percent of plots with at least one tree seedling or sapling per 0.0016 ha, directly for the 1935 data. Density measurements alone can be misleading because of the known spatial heterogeneity of natural tree regeneration [30,44]. There have been several attempts at relating understory tree density and stocking [31,44,45], because stocking is a more reliable measure in assessing management objectives or desirable tree occupancy [43] and accounts for the spatial heterogeneity of natural regeneration. More recently Ritchie [32], using data from the 2000 Storrie fire on the Lassen National Forest, conducted a validation of the Lynch and Schumacher [31] probit model that predicts stocking as a function of understory tree density, and found that the Storrie fire data matched well with the Lynch and Schumacher [31] model fit. The results from Ritchie [32] were applied here to estimate stocking from the compiled, available seedling density data: where Φ −1 (p) is the inverse of the standard normal (probit) of proportion (p) of stocked plots and x is mean number of trees per plot. To evaluate the fit of the new equation we plotted our observed stocking in 2018 at 0.0016 ha plots over the predicted stocking 95% confidence interval from Ritchie [32].

Explanatory Variables for 1935 Conditions
Three groups of explanatory variables were used: overstory, soil, and topography. We compiled overstory variables-trees per hectare and overstory basal area per hectare (Table 2)-for each of 100 compartments using available census data from 1934. In 1934, all live trees >8.9 cm dbh were tallied by 5.08 cm wide diameter classes within 4055 plots (~1 ha in size) throughout BMEF. If a 1-ha plot straddled two compartments, the overstory variables were weighted by the proportion of the plot that fell into a particular compartment. Two compartments that were primarily composed of meadows were dropped, and we used the remaining 98 compartments for any analyses that included the 1934 overstory data (Table 3).  Topographic variables-aspect, slope, and elevation-were extracted from a digital elevation model (DEM) at a 1-m resolution ( Table 2). The zonal statistic function from ArcMap 10.4.1 summarized the values of the DEM in terms of average elevation, slope, and aspect within the extent specified by the shape file for the 100 compartments [47]. Aspect and slope were incorporated into the model using the transformations suggested by Stage [48]. Vegetation information for BMEF was available from an ecological unit inventory (EUI) shapefile. This EUI map, developed for BMEF in 1994, identified different geological landforms, soils, and potential vegetation types [46]. The compartment map was overlaid with the EUI map and the intersect function in ArcMap 10.4.1 was used to extract the soil variables (Table 2) from the EUI map at compartment level [49]. Multiple categories of the variables could fall in a single compartment. If a compartment straddled multiple soil variable categories, an area-weighted average of the midpoints was assigned to the compartment.
Hydrological processes and the control of topography over such processes can be quantified by the topographic wetness index [50]. For the 100 compartments, topographic wetness index was calculated based on the slope from a 1 m resolution DEM, using the equations provided by Beven and Kirkby [50]. Similarly, solar insolation, which takes into account topographic characteristics such as slope, aspect, and elevation [51], was calculated for each of the 100 compartments using the DEM and Spatial Analyst tool from ArcMap 10.4.1.

Historical Stocking
We conducted exploratory data analyses for the 1935 regeneration data using histograms, with stocking in each compartment and frequency on the xand y-axis, respectively. Similarly, we constructed histograms for stocking by height class and species for the historical data from 1935.

Effect of Overstory, Soil, and Topography on Stocking
Meadows were part of the historical forest landscape. For two out of the 100 compartments, the area of meadows exceeded 70%, which made these compartments different from the other compartments, for which the area of meadows was below 40%, with 77 compartments not having any meadows. Therefore, we dropped the two compartments that were dominated by meadows from our analysis and used the remaining 98 compartments in the model development. Stocking is continuous and restricted to the standard unit interval (0, 1). Therefore, we fit beta regression models that take the distribution of stocking into account [52] using the 'betareg' package in R version 3.6.2 [53] with a logit link function [52]. We developed the models as follows: Our first step was to identify the variables that were significantly associated with stocking. We fit univariable beta regression models for all explanatory variables (Table 2). Then, we selected variables with low p-value (p < 0.05) to be used in multivariable beta regression models across the three groups of explanatory variables (Table 2), keeping variables significant at α-level of 0.05 and dropping variables that had variance inflation factor greater than three [54]. This resulted in a model that contained available water capacity at 150 cm (AWC 150 ) and overstory basal area per hectare as explanatory variables. However, in our exploratory data analysis, the overstory basal area per hectare showed a curvilinear relationship with stocking. In order to account for the trend observed in the data and to mimic this trend, we decided to use a range of exponents between 0.10 to 2 at the interval of 0.10 for the overstory basal area per hectare to identify the value of exponent that resulted in the best model. We selected the model that resulted in high pseudo-R 2 and high precision parameter (phi) values [52] as well as a low Akaike Information Criterion (AIC) value [55]. We also checked the assumption of constant variance and goodness-of-fit using residual vs linear predictor and half normal plots [56], respectively, for all beta regression models [50,51].

Stocking in Contemporary Forests under Different Management Treatments
The effect of six different management treatments on contemporary stocking was investigated using a Kruskal-Wallis non-parametric one-way ANOVA [57]. The treatment factor had six levels ( Table 1). When differences among treatment levels were detected, a pairwise Wilcoxon rank sum test [58] was used to determine which treatment levels had different stocking.

Stocking in Historical vs. Contemporary Forests
Stocking by species and height class was available for 12 compartments in both measurement years. Dot plots were constructed for both measurement years for overall stocking, stocking by three species, and stocking by four height classes to investigate if contemporary stocking was higher than historical stocking.

High Stocking for Ponderosa Pine Species in Historical BMEF
Our results indicated that the stocking observed in 2018 for 0.0016 ha plots corresponded well, for observations below about 1000 seedlings per hectare, with the 95% confidence limits of predicted 0.0016 ha stocking from Ritchie's [32] equation (Figure 2). Although there is evidence of some lack of fit at the highest levels of density where model predictions of stocking may be low by around 10%. However, since for very high levels of density the model approaches one, errors in estimation will actually start to decrease for the extreme density levels observed in 1935. In 1935, the majority of compartments exhibited higher seedling densities than observed in the 12 compartments for 2018 (Table 3). Therefore, it is possible that the stocking predicted from Ritchie's [32] equation may be a bit lower than actual stocking. However, we do not know for sure how the equation performs for high seedling densities over 1500 trees ha −1 , because we did not have any 2018 data in that range for validation ( Figure 2). Stocking by species and height class was available for 12 compartments in both measurement years. Dot plots were constructed for both measurement years for overall stocking, stocking by three species, and stocking by four height classes to investigate if contemporary stocking was higher than historical stocking.

High Stocking for Ponderosa Pine Species in Historical BMEF
Our results indicated that the stocking observed in 2018 for 0.0016 ha plots corresponded well, for observations below about 1000 seedlings per hectare, with the 95% confidence limits of predicted 0.0016 ha stocking from Ritchie's [32] equation (Figure 2). Although there is evidence of some lack of fit at the highest levels of density where model predictions of stocking may be low by around 10%. However, since for very high levels of density the model approaches one, errors in estimation will actually start to decrease for the extreme density levels observed in 1935. In 1935, the majority of compartments exhibited higher seedling densities than observed in the 12 compartments for 2018 (Table 3). Therefore, it is possible that the stocking predicted from Ritchie's [32] equation may be a bit lower than actual stocking. However, we do not know for sure how the equation performs for high seedling densities over 1500 trees ha −1 , because we did not have any 2018 data in that range for validation ( Figure 2). The estimated stocking in 1935 was dominated by ponderosa pine (Figure 3b). Our data indicated that there were no white fir and incense-cedar seedlings in 38 compartments in 1935. Most of the incense-cedar and white fir regeneration was present in the compartments at higher elevations, located in the northeast portion of BMEF. Ponderosa pine had highest stocking (mean proportion = 0.53 and median proportion = 0.54), followed by white fir (mean proportion = 0.10 and median The estimated stocking in 1935 was dominated by ponderosa pine (Figure 3b). Our data indicated that there were no white fir and incense-cedar seedlings in 38 compartments in 1935. Most of the incense-cedar and white fir regeneration was present in the compartments at higher elevations, located in the northeast portion of BMEF. Ponderosa pine had highest stocking (mean proportion = 0.53 and median proportion = 0.54), followed by white fir (mean proportion = 0.10 and median proportion = 0.21) and incense-cedar (mean proportion = 0.10 and median proportion = 0.02) (Figure 3).  The majority of the compartments were dominated by seedlings in height classes 2 and 3 irrespective of species (Figure 4d-i). However, among the three species ponderosa pine had the highest stocking in height classes 2 (mean = 0.33 and median = 33) and 3 (mean = 0.33 and median = 0.34) (Figure 4d,g). Stocking in height classes 1 and 4, were low compared to other height classes for all three species (Figure 4). The majority of the compartments were dominated by seedlings in height classes 2 and 3 irrespective of species (Figure 4d-i). However, among the three species ponderosa pine had the highest stocking in height classes 2 (mean = 0.33 and median = 33) and 3 (mean = 0.33 and median = 0.34) (Figure 4d,g). Stocking in height classes 1 and 4, were low compared to other height classes for all three species (Figure 4).

Overstory and Soil Variables were Associated with Stocking
We chose the model with the lowest AIC value (−191.86) and highest pseudo-R 2 (0.49) and precision parameter (29.27) compared to the other tested models (Table 4). However, given the small differences in AIC values among the models (Table 4), this model is not necessarily superior to the other models in terms of model fit. This final beta regression model showed that basal area per hectare raised to an exponent of 0.10 (p-value < 0.001) and available water capacity at 150 cm (p-value < 0.001) were significantly associated with historical stocking (Table 4).

Overstory and Soil Variables Were Associated with Stocking
We chose the model with the lowest AIC value (−191.86) and highest pseudo-R 2 (0.49) and precision parameter (29.27) compared to the other tested models (Table 4). However, given the small differences in AIC values among the models (Table 4), this model is not necessarily superior to the other models in terms of model fit. This final beta regression model showed that basal area per hectare raised to an exponent of 0.10 (p-value < 0.001) and available water capacity at 150 cm (p-value < 0.001) were significantly associated with historical stocking (Table 4). Stocking appeared to increase with increasing overstory basal area per hectare up to about 12 m 2 ha −1 ; above this level the relationship hits a plateau (Figure 5a). Available water capacity at 150 cm had a positive relationship with stocking. Stocking increased as more water was available to plants at depth of 150 cm ( Table 4). The slope-aspect transformations [48] had p-values between 0.04 and 0.07. Because the inclusion of these model terms did not improve the model substantially but increased model complexity, we decided not to include them. None of the remaining topography variables-elevation, solar insolation, and topographic wetness index-were significant in the model.

Differences in Stocking among Management Treatments in 2018
Stocking in 2018 varied across treatments (Kruskal-Wallis test p-value < 0.0001) (Figure 6a). Control (C) compartments had the highest stocking followed by burned (B) compartments (p-value < 0.0001; Figure 6a). The Wilcoxon rank sum test revealed that stocking in TC-B and TB were not significantly different from each other (p-value = 0.0816). The difference in stocking between TC-B and TB-B (p-value = 0.0028) and TC and TB-B (p-value = 0.0014) was significant. Significant differences in stocking (p < 0.0001) were found between all other treatment levels.

Historical Stocking Higher than Contemporary Stocking
Estimated stocking in 1935 was greater than observed stocking in 2018 across all the compartments (Figure 6a). Control treatments had the highest stocking in 2018 compared to the other treatments. These two compartments also had the highest stocking in 1935 (Figure 6a). The majority of compartments were dominated by stocking of ponderosa pine followed by incense-cedar in both historical 1935 and contemporary 2018 conditions (Figure 6b). Stocking in 1935 was dominated by height classes 2 and 3 ( Figure 6c). In 2018 stocking was dominated by height class 4 in control compartments (Figure 6c), while compartments that received the burn treatment were mostly stocked with seedlings from height classes 1 and 4 (Figure 6c). The compartments that received thinning treatments had very low stocking, with height classes 1 and 2 being a bit more abundant than height classes 3 and 4 (see TC-B, TC and TB in Figure 6c). Forests 2020, 11, x FOR PEER REVIEW 13 of 22

Discussion
Our research was able to meet four objectives, which were to 1) understand historical stocking in terms of species and height class, 2) determine the factors associated with historical stocking, 3) understand the contemporary stocking across various management treatments, and 4) compare historical stocking with contemporary stocking across BMEF. Prior to the onset of any forest management activities, the historical stocking was dominated by ponderosa pine in height classes 2 and 3 throughout BMEF. The major drivers of stocking in the historical forest were overstory basal area and soil moisture. There were differences in stocking among management treatments in contemporary forests, and stocking was higher in historical than in contemporary forests.

Discussion
Our research was able to meet four objectives, which were to (1) understand historical stocking in terms of species and height class, (2) determine the factors associated with historical stocking, (3) understand the contemporary stocking across various management treatments, and (4) compare historical stocking with contemporary stocking across BMEF. Prior to the onset of any forest management activities, the historical stocking was dominated by ponderosa pine in height classes 2 and 3 throughout BMEF. The major drivers of stocking in the historical forest were overstory basal area and soil moisture. There were differences in stocking among management treatments in contemporary forests, and stocking was higher in historical than in contemporary forests.

Historical Stocking Was Dominated by Ponderosa Pine at Low Elevations
In the early 1930s, after sheep grazing [33] and prior to the onset of timber harvesting activities, low elevation sites of BMEF were dominated by ponderosa pine stands [34]. The landscape in BMEF has gentle slopes with the majority of the forest in south or southwestern aspect and with an average slope of less than 15% [34]. Such south-facing gentle slopes have been reported to favor seedling germination of ponderosa pine in northern California [59]. The observed dominance of ponderosa pine stocking in these low elevation sites of BMEF can further be explained by overstory basal area per hectare being a major driver of stocking in the historical BMEF. Ponderosa pine dominated the overstory, and its presence not only facilitates the establishment of regeneration by supplying ample seeds [60], but also ameliorates the harsher environmental conditions caused by high temperatures, solar radiation, and wind [61].
Moderate overstory basal area, ranging from 7 to 23 m 2 ha −1 -96% of BMEF fell within this range-is reported to be conducive for establishment of ponderosa pine seedlings [4,62]. In 1934, BMEF exhibited low to moderate overstory basal area for ponderosa pine (0 to 25 m 2 ha −1 ) and low overstory basal area (0 to 7 m 2 ha −1 ) for incense-cedar and white fir. The observed low to moderate overstory basal area maintained an open park-like structure in the majority of BMEF, which favored germination and growth of shade-intolerant ponderosa pine seedlings.
The observed high stocking for incense-cedar and white fir in some of the compartments in the northeast portion of BMEF may be due to these compartments being located in more mesic and higher elevation sites. Mesic and higher elevation sites have been reported to be more conducive to these species [59]. A mixture of all three species in the overstory may have created a dense overstory, which provided enough shade for the growth of shade tolerant incense-cedar and white fir [63]. At the same time, the presence of mature incense-cedar and white fir trees provided the seeds necessary for germination [63].

Historical Stocking Was Dominated by Height Classes 2 and 3
Land survey records across dry forests in Oregon, California, and parts of northern Arizona in the late 1800s showed that regeneration-seedlings and/or saplings-were observed over 35-57% of the dry-forests [64]. The absence of surface fires after Euro-American settlement failed to control their growth and establishment in these areas [64]. At the time of the 1935 regeneration data collection, BMEF had experienced over 35 fire-free years (unpublished fire scar data from Carl Skinner, United States Forest Services, retired). This fire free period may be attributed to heavy sheep grazing in the area [33]. Hence, seedlings that were established around the early 1900s due to the absence of fire could have stagnated in height growth throughout BMEF. Stand stagnation could explain the observed higher stocking in height classes 2 (0.31-0.91 m) and 3 (0.91-1.83 m) compared to height classes 1 and 4 throughout BMEF in 1935. Lacking age data for the 1935 forest conditions, we are unable to make any specific comments on height and age dynamics of seedlings within the historical ponderosa pine forest in BMEF. Oliver [65] reported dense thickets of saplings that had an average height of 2.7 m at the age of 70 within BMEF as a result of stand stagnation. Similarly, Weaver [9] reported stand stagnation after the establishment of regeneration in ponderosa pine forests in Washington and Oregon, USA, because of fire exclusion. Therefore, the seedlings/saplings that were observed in smaller height classes 1 (0-0.31 m) and 2 could be as old as seedlings/saplings in height classes 3 and 4 (>1.3 m) due to stagnation.
Differences in stocking for various height classes could also be attributed to a lag in regeneration-a period of no regeneration-between two pulses of regeneration [11]. Regeneration lags occur due to the periodicity of ponderosa pine, which is driven by heavy cone production years followed by years of few cones or no cones at all [66]. Ponderosa pine can experience a regeneration lag from 15 to 70 years between the establishments of seedlings [67,68]. During this regeneration lag, ponderosa pine can be predisposed to unusual climatic episodes that can fluctuate between wetter fire-free and drier fire-prone periods [69]. If a long drier and fire-prone period prevails, there will be a longer gap in regeneration before another good seed year [70]. Hence, a gap in regeneration might result in the differentiation of seedlings in various height classes. The differentiation of stocking in various height classes can also be attributed to variability in growth due to site effects [71] and localized competition for nutrients and soil moisture [10]. For example, seedlings that germinated at the same time, but in favorable sites might have gained growth in height compared to seedlings growing in poorer sites, which lagged in height growth. Oliver and Powers [71] in their study in the western Sierra Nevada, southern Cascade Range, and Warner Mountains of northeastern California found that seedlings that germinated at the same time reached diameter at breast height at 13.5 years of age on poor-sites as compared to 4.4 years on better sites.

Available Water Capacity at 150 cm and Overstory Basal Area Were Drivers of Historical Stocking
Interactions of topographic variables-elevation, slope, and aspect-result in heterogeneous forest structure in the dry forests across the Sierra Nevada and Southern Cascades, usually through changing the patterns and amount of regeneration along topographic gradients [72]. However, in our study, topographic variables such as elevation, slope, and aspect were not found to be major drivers of stocking. This is possibly due to the limited range in elevation-1700 to 2100 m-and a gentle slope of less than 15% on average [34].
A good seed crop followed by a favorable moisture regime and receptive seed bed is needed for the establishment of ponderosa pine seedlings [11]. We observed that within BMEF available water capacity at the depth of 150 cm was one driver of stocking in 1935. A large number of seedlings possibly germinated due to ample moisture availability within their rooting zone [59,73,74]. Ponderosa pine seedlings have been shown to survive in areas with low precipitation by maintaining high water use efficiencies and root systems-12 to 50 cm deep-that utilize moisture available in their rooting zone [73,75].
Stocking showed a curvilinear relationship with increasing overstory basal area per hectare (Figure 5a). Overstory basal area was dominated by ponderosa pine basal area and therefore, mainly driven by the presence of ponderosa pine. However, the observed highest stocking and overstory densities in higher-elevation compartments in the northeast portion of BMEF could be due to the presence of all three species [63,76,77]. While the observed relationship between overstory basal area per hectare and stocking was stronger than the relationship between available water capacity at the depth of 150 cm, the overall model fit was not very good, as suggested by the moderately low pseudo-R 2 of 0.49. Stocking tended to be slightly overpredicted for lower overstory basal area per hectare ranges, and variability in the data was quite high.
The fitted beta regression model provides general insights about the variables that explain stocking in the historical BMEF environment. Our model, however, comes with the limitations that the soil variables used in the model were derived from an Ecological Unit Inventory (EUI) map produced in 1994 [46] decades after 1935, when the historical regeneration data were collected. In the absence of field-collected soil data, we have to make an assumption that the information on soil variables derived from the EUI map represents the soil conditions that existed in BMEF around 1935 and has not changed. Additional variability in stocking could possibly be explained by incorporating fire history variables such as time since fire, fire frequency, and fire return interval into the model [8]. However, in this study, availability of fire history variables was limited to 28 of the 100 compartments. The fire history variables were statistically significant when we added them to the model using the 28 compartments. Considering the small sample size of 28, the model was overfit with these additional variables and we did not think it was meaningful to report the results in detail. Yet, this exploratory analysis confirms that fire history variables could explain some of the observed variability in stocking.

Contemporary Stocking Was Low in Thinned Stands and Differed across Treatments
In 2018, compartments that received thinning treatments-thinning from below and combination thinning-with or without prescribed fire had lower stocking than the compartments that received the burn only treatment. The mechanical thinning treatments mostly removed the trees up to a diameter of 40.38 cm including the seedlings and saplings up to a diameter of~7 cm, which can explain the low observed stocking in these compartments [34,37]. In addition, within the compartments that received prescribed fires after the thinning treatments, fire might have killed additional seedlings resulting in lower stocking compared to treatments that received thinning only. Another possible reason for low stocking within the combination thinning and thinning from below could be the lack of seed sources near thinned sites [76,77]. Within the thinned sites, due to reductions in canopy cover, seedlings are likely to get exposed to high levels of direct solar radiation [78,79]. This may lead to seedling desiccation and death within thinned compartments [59], where regeneration was already removed during the thinning treatments.
In contrast, higher stocking in the burn only treatment compared to the treatments that received thinning can possibly be attributed to the fire having reduced heavy accumulation of litter and duff, which allowed seeds to reach mineral soil and extract adequate soil moisture for germination [73]. At the same time, the overstory that survived the prescribed fire, provided favorable environmental conditions at microsites for germination and growth of seedlings [73]. The presented results for the burn only treatment should be taken with caution, as time since management might have played an important role in determining the amount of stocking across various management scenarios at BMEF. Burn only treatments were implemented in 1997 while the compartments that received thinning from below in 2012 were burned in 2017 [41]. We collected the regeneration data in 2018. There was a gap of 20 years between application of treatments and data collection in the burn only treatment. During this 20-year window, seedlings might have established in the burn only treatments and resulted in the observed higher stocking. In the compartments that received thinning from below and combination thinning with or without fire, there was a gap of only one year between application of fire and data collection. This may explain the observed lower stocking in these treatments. Because time since treatment was different for the different treatments, the observed treatment effect may be confounded with time.
The observed lower stocking in burned compartments compared to the control compartments may be explained by the mortality of the seedlings over time due to crown scorch and cambium wounding [17]. In the absence of any management activities and fire for many decades in the control compartments [34], the seedlings could have stagnated within the stands during various episodes of recruitment, which resulted in high stocking [65]. In addition, the control compartments were located at somewhat higher elevations and within mesic sites than the burn only treatments. These conditions are conducive in terms of available microsites for shade-tolerant white fir and incense-cedar [16], which was confirmed by the higher stocking observed for incense-cedar and white-fir in the control compartments. Therefore, within these compartments higher stocking may also be explained by the higher abundance of seedlings of the shade-tolerant species and a more mesic environment, rather than by the fact that they did not receive any management.

Higher Historical Stocking than in Managed, Contemporary Stands
The observed lower stocking in 2018 compared to 1935 may indicate that the contemporary forest structure is not as conducive to regeneration as it used to be prior to 1935. Due to infilling of trees, there may be lack of growing space in contemporary forests [80]. In 1935, however, the forest structure was more open, maintained by frequent surface fires, which reduced shading and resource competition and prepared fine mineral seedbeds for seedling germination [17]. However, the differences in stocking between historical and contemporary forests could also simply be due to the management treatments that were applied in the contemporary forests, where the mechanical thinning treatments removed most of the regeneration and prescribed fires may have killed additional regeneration [34,37].

Conclusions
The trajectories of future forest stand structures in ponderosa pine forests of northern California depend on the regeneration response to management and environmental conditions. Historical regeneration studies allow managers to understand the relationships between regeneration and environmental factors, which in turn shape the overall stand structure.
Our field inventory data driven study showed that historically, the regeneration in BMEF was dense and dominated by shade-intolerant ponderosa pine seedlings in height classes 2 and 3 at the lower elevations of the forest. More mesic sites at higher elevations exhibited high stocking of shade-tolerant incense-cedar and white fir. Seedlings that were established around the early 1900s found the available growing space in open park-like conditions favoring the germination and growth of ponderosa pine seedlings. Additionally, in absence of fire to check the growth, seedlings could have stagnated in height growth throughout BMEF, which resulted in high stocking in height classes 2 and 3.
Overstory basal area per hectare and available water capacity at 150 cm depth were related to historical stocking. Stocking increased with both increasing overstory basal area as well as increasing available water capacity. Topography variables were not significant in the final beta regression model. The amount of regeneration was different across management treatments in contemporary conditions. The time since management played a key role in determining the differences in stocking within contemporary BMEF. Compartments that were recently thinned and burned had lower stocking compared to compartments that were burned and thinned 20 years prior.
The stocking in the contemporary BMEF was lower than the stocking in the historical forests. The primary reason for lower stocking in contemporary conditions was because most of the small trees, seedlings, and saplings were either removed by thinning or killed due to fire by recent management activities.
The results of this study are unique because much of our current knowledge on historical regeneration has been based on dendro-chronological reconstructions, which cannot provide the kind of detailed information about historical regeneration conditions that we have presented in this study.