Reproductive maturity and cone abundance vary with tree size and stand basal area for two widely distributed conifers

. Understanding potential limitations to tree regeneration is essential as rates of tree mortality increase in response to direct (extreme drought) and indirect (bark beetle outbreaks, wild ﬁ re) effects of a warming climate. Seed availability is increasingly recognized as an important limitation for tree regeneration. High variability in seed cone production is a trait common among many northern temperate conifers, but few studies examine the determinants of individual tree cone production and how they vary with stand structure. In subalpine forests in the southern Rocky Mountains, USA, we monitored >1600 Picea engelmannii (Engelmann spruce) and Abies lasiocarpa (subalpine ﬁ r) trees for cone presence (an indicator of reproductive maturity) and a subset of those trees for cone abundance (an indicator of seed production) from 2016 to 2018. We constructed mixed models to test how individual tree cone presence and cone abundance were affected by tree size and age as well as forest attributes at the neighborhood- and stand-scales. The probability of cone presence and cone abundance increased with tree size and age for A. lasiocarpa and P. engelmannii . The youngest ages of trees with cones present were more than 100 yr later for individuals in high basal area (BA) stands (>65 m 2 /ha) relative to low BA stands (<25 m 2 /ha). P. engelmannii produced many more cones than A. lasiocarpa at similar sizes, especially in young, low BA stands. Our ﬁ ndings reveal how differences in tree sizes and stand structures typically associated with time since last disturbance can affect seed production patterns for decades to well over a century. The consistent regional pattern of earlier and more abundant post ﬁ re establishment of P. engelmannnii vs. the delayed post ﬁ re establishment by A. lasiocarpa may be partially explained by species ’ differences in cone abundance by stand structure. The increasing loss of large, dominant cone-producing trees will signi ﬁ cantly reduce seed production to support future tree regeneration and maintain forest cover. However, seed availability and resilience following disturbances may be less limiting than expected for species like P. engelmannii that have the capacity to produce more cones in open-canopy forests, such as recently disturbed areas.


INTRODUCTION
Tree regeneration is required for the persistence of current forest cover and to permit trees to colonize new habitats that become available as climate changes (Walck et al. 2011, Corlett andWestcott 2013) or following coarse-scale disturbance (Pickett et al. 1987). Identifying limitations for tree regeneration may help forecast shifts in forest cover (Enright et al. 2015). For trees that primarily or exclusively reproduce via seed (i.e., obligate seeders), seed availability is an important limitation for tree regeneration (Grubb 1977, Grime 2001. In the context of increasing rates of tree mortality from the direct (e.g., extreme drought) and indirect (e.g., increased fire and insect outbreaks; van Mantgem et al. 2009, Berner et al. 2017, Senf et al. 2018) effects of climate warming, understanding spatial and temporal dynamics of seed production is critical for assessing the potential limitations for tree regeneration.
Studies of seed production for many tree species across the globe indicate that seed production is highly variable from year to year, but spatially and temporally synchronous in quantity among individuals within populations (Kelly 1994, Pearse et al. 2016). This widely studied phenomenon, called mast seeding, requires the alignment of many reproductive processes from flower initiation to seed maturity (Owens 1995) and is hypothesized to be driven by the complex interplay between environmental cues and internal plant resource dynamics Rapp 2014, Pearse et al. 2016). Years of high seed production are important for satiating appetites of animal populations, which allows some seeds to escape predation (Janzen 1971). The surplus often results in widespread seedling establishment under suitable environmental conditions such as the case of the 1918-20 pulse of Pinus ponderosa establishment in the southwestern United States (Pearson 1923).
Studies of northern temperate conifers highlight several-fold differences in seed cone production (hereafter, cone production) among individual trees of the same size and across tree sizes within a year (e.g., Lamontagne andBoutin 2007, Davi et al. 2016). Such variability may be caused by a tree's access to, and allocation of, resources for reproductive processes (Greene et al. 2002, Sala et al. 2012. Indeed, cone production is influenced by characteristics of an individual tree (tree-scale), competition from a focal trees' immediate neighbors (meters, neighborhood-scale), and variability in structure among stands of trees (hectares, stand-scale; Calama andMontero 2006, Nygren et al. 2017). Quantifying variability within and among individuals in a population is necessary to understand the effects of climate warming on tree demographic processes (Enright et al. 2015, Davis et al. 2018.
To initiate cone production (i.e., reproductive maturity), trees must reach a minimum age and/ or size (Owens 1995). Once that threshold size or age is attained, cone production generally increases with tree size (diameter at breast height [dbh] or crown size; Greene et al. 2002, Davi et al. 2016) and age (Viglas et al. 2013). Larger trees (and generally older trees) support greater crown area for male and female flower buds as well as greater access to above (e.g., light; Shea 1987, Greene et al. 2002 and belowground resources (e.g., water, nitrogen and phosphorus; Sala et al. 2012). For some conifers, tree size has been shown to be a better indicator of cone production than tree age (Linhart andMitton 1985, Viglas et al. 2013), possibly because tree size more strongly reflects a tree's access to resources than tree age. Neighborhood-and stand-scale variation in tree densities and sizes, often created by disturbances (e.g., treefalls, wildfire), can affect cone production by altering competition for light and nutrients (Canham et al. 2006, Haymes andFox 2012) or influencing a focal tree's characteristics, such as crown size (Greene et al. 2002). For example, Picea abies individuals in low basal area (BA) stands can produce many more cones than comparably sized individuals in stands of higher BA (Nygren et al. 2017). Differences among species in the allocation of resources to growth, reproduction, and defense and species' response to environmental conditions (e.g., shade tolerance) further complicate the detection of patterns in reproductive parameters like cone production (Koenig and Knops 1998). Improved understanding of the factors limiting cone production is needed for assessing forest resilience and possible shifts in species composition in the context of increased rates of tree mortality (Redmond and Barger 2013, Temperli et al. 2015.
Picea engelmannii Parry ex Engelm. (Engelmann spruce) and Abies lasiocarpa (Hook.) Nutt. (subalpine fir) are two obligate seeding species for which little is known about the tree-to stand-scale determinants of individual tree cone production, despite their widespread distribution across high-elevation forests of western North America. At broad spatial scales, synchronous production of abundant quantities of seed has been linked to climatic cues for A. lasiocarpa in the Pacific Northwest (Woodward et al. 1994) and P. engelmannii in the southern Rocky Mountains (Buechling et al. 2016). During intervening years, low to intermediate levels of seed production (including some years of no cone production) are common for both species (Alexander 1987). In subalpine forests of the southern Rocky Mountains, variability in tree sizes and stand structures related to coarse-scale disturbances (e.g., bark beetles, fire, blowdown) as well as sitespecific tree growth responses to climatic variability (Veblen 1986, Veblen et al. 1994, Buechling et al. 2017) are likely to influence resource availability for cone production. Silvicultural studies have established the minimum size and age of reproductive maturity for spruce and fir (Alexander 1987), but they do not account for the spatial variability in reproductive maturity or cone production among individuals created by differences in stand structure. Better knowledge of the factors affecting cone production and the spatial variability in cone production is needed for more realistic parameterization of demographic and landscape models that examine climate-caused shifts in forest cover and distribution (Schumacher et al. 2006, Temperli et al. 2015, Conlisk et al. 2017, Davis et al. 2018. We examined >1600 P. engelmannii and A. lasiocarpa trees for cone presence (an indicator of reproductive maturity) and a subset of those trees for cone abundance (an indicator of seed production) for 3 yr (2016)(2017)(2018) in the Colorado Front Range (CFR), Colorado, USA. Our objective was to investigate how cone presence and cone abundance are affected by individual tree characteristics (tree-scale), competition from neighboring trees (neighborhood-scale), and variability among stands in structure (standscale). Based on findings for other northern temperate conifers, including other species of Picea and Abies, we expect the following for P. engelmannii and A. lasiocarpa: (1) The more shade-tolerant A. lasiocarpa will reach reproductive maturity at an earlier age and/or smaller size than P. engelmanni, (2) tree size will be a stronger predictor of cone presence and cone abundance (for trees producing cones) than tree age, and (3) greater competition at a neighborhood-scale (e.g., 0-6 m from the focal tree) and greater total BA at a stand-scale will increase the minimum tree age and size threshold required for cone presence and decrease cones per tree.

Study area
The study was conducted in forests dominated by P. engelmannii and A. lasiocarpa on the eastern slope of the CFR (Fig. 1). Picea-Abies forests in the CFR extend from the lower ecotone of the subalpine zone (~2600 m) to alpine treeline (~3500 m). Most stands included in the study were composed exclusively of P. engelmannii and A. lasiocarpa, but some stands included Pinus contorta Douglas ex Loudon (lodgepole pine), Populus tremuloides Michx. (aspen), and Pinus flexilis James (limber pine). Both Pinus species and P. tremuiloides are considered to be seral to P. engelmannii and A. lasiocarpa (Peet 1981, Veblen 1986). The subalpine zone is characterized by a continental climate with long, cold winters and a short, dry growing season. The majority of mean annual precipitation ( precipitation generally falling as rain during convective storms (Kittel et al. 2015). No trends in precipitation have been detected from 1978 to 2010 (Kittel et al. 2015), but temperature increases reflect recent warming trends across western North America. From 1953 to 2016, the average annual temperature was 1.8°C with temperatures reaching an average monthly minimum in January of À1.9°C and an average monthly maximum in July of 19.8°C (NWT LTER 2016). Mean (0.2°C per decade) and maximum (0.44°C per decade) annual average temperatures increased in the subalpine zone in the CFR from 1953 to 2008(McGuire et al. 2012.

Field data collection
Data on cone presence were collected in 15 plots, and data on cone abundance were collected in 10 of the 15 plots located in three general sampling areas: (1) Rollins Pass and (2) Niwot Ridge Long Term Ecological Research site at University of Colorado's Mountain Research Station, and (3) Rocky Mountain National Park (Fig. 1, Table 1). Seven permanent forest plots (hereafter, permanent plots) were located at Niwot Ridge, and eight temporary plots (hereafter, temporary plots) were located at Rollins Pass (four plots), Rocky Mountain National Park (three plots), and Niwot Ridge (one plot). Temporary plots were selected to broaden the range of stand structure characteristics (e.g., stand ages since stand-replacing fire, stand density, and stand BA) represented by the permanent plots to allow us to test stand-scale effects on cone presence/absence and cone abundance. Collectively, plots were located in stands of postfire origin that ranged in age from~120 yr (young postfire) to~250-375 yr (old seral postfire) or stands without evidence of fire for >450 yr (old Picea-Abies stands; Veblen 1986, Kienast and Schweingruber  Rebertus et al. 1991). Stand densities in the more recently burned stands are overlapping with the range of densities in much older (>350 yr) stands, but forest canopies are more open and stand BA is considerably lower in these young, postfire stands (Table 1).
In the permanent plots, we randomly selected and observed >25 P. engelmannii and >25 A. lasiocarpa trees ranging in size from 4 cm dbh to the maximum tree size in each stand. Five permanent plots were observed for cone presence in a single year (2015) and the other two plots were observed for cone presence and cone abundance from 2016 to 2018 (Table 1). All observations of the cone-bearing area (top one-third of tree) were made using binoculars (Eagle Optics Triumph, 8 9 25 power). Cone presence for each tree was defined as the presence of any of the following structures: woody cone axis (with all seeds abscised), fully developed cones from previous years, or new, developing (but immature) cones. If none of these structures were present, the tree was marked as absent of cones. We included all cone structures because our goal was to quantify the capacity to produce cones rather than the interannual variability in cone presence (see Appendix S1 for sensitivity analysis of sampling cone presence/absence in a single year vs. multiple years).
In the temporary plots, we randomly located a plot center within a stand (i.e., an area~1 ha or larger with relatively uniform canopy closure, species composition, and BA) and selected a restricted random sample of more than six but less than ten trees (>4 cm dbh) in three height classes (codominant, intermediate, subcanopy) for P. engelmannii (>24 total) and A. lasiocarpa (>24 total) in summer 2016. An additional 2-4 juveniles (1-4 cm dbh) per plot were randomly selected for cone observations in 2018 after preliminary logistic models using Vuong's test (Vuong 1989) demonstrated that dbh was a stronger predictor of cone presence than tree age for both P. engelmannii and A. lasiocarpa. We measured dbh (cm), height (m), height class (codominant, intermediate, subcanopy), canopy dimensions (shortest and longest canopy radii at canopy base, canopy base height), and number of trees in three size classes (canopy, intermediate, and subcanopy) within 1, 3, and 6 m of the focal tree to calculate a competition index. To estimate tree age of all trees >4 cm dbh (trees < 4 cm dbh excluded) sampled for cone presence and abundance, we extracted tree cores at~20 cm above the ground. All age estimates were corrected for missing rings between innermost ring on tree core and the chronological center of the tree following Duncan (1989) and age at coring height following Villalba and Veblen (1997). We conducted a survey of cone presence (defined above) and cone abundance (see next section) on individual trees in the temporary plots from 2016 to 2018 using the same methods as above. To quantify stand structure in the temporary plots, we measured dbh (cm) and recorded species of all trees (>4 cm dbh) in a 10 9 40 m area centered on plot center.

Estimating cone abundance
We estimated cone abundance in late August (prior to the annual period of seed release in September) using two methods to accommodate the wide range in P. engelmannii cone abundance (0 to >1000 cones). For trees with <100 cones, cone abundance was estimated as the sum of cones counted on two sides of the tree (cones on side 2 visible from side 1 not counted) by one observer (always R. Andrus) located at a distance equal to the height of the tree (max. tree height, 20 m). For trees with >100 cones, we estimated cone abundance from photographs taken of the cone production area (see Appendix S1 for details; Lamontagne and Boutin 2007).

Statistical analysis
Statistical models of the probability of cone presence and total cone abundance were created to test the effect of tree-, neighborhood-, and stand-scale predictor variables. Specifically, we selected tree dbh (cm), tree age, and crown volume (m 3 , calculated using the geometric formula for a cone with an elliptical base) as well as species (P. engelmannii and A. lasiocarpa) as tree-scale predictors. Tree crown measurements to compute crown volume were only available in a subset of the trees in the temporary plots. As such, we computed pairwise spearman correlations among crown volume, dbh, and tree age instead of including crown volume in the models. At the neighborhood-scale, we calculated a neighborhood competition index for each focal tree at two spatial scales (see below). At the stand-scale, we quantified stand structure with the following three measures: (1) age of oldest tree as an estimate of time since last fire, (2) total live stem density (stems/ha), and (3) total live BA (m 2 /ha). In the temporary plots, stand structure variables (BA and stem density) were derived from the 10 9 40 m plot. In the permanent plots, tree-, neighborhood-, or stand-scale variables were available or derived from a 2016 census (dbh [cm], BA, stem/ha; Andrus et al. 2018b) or the original plot installation (tree ages, age of oldest tree, and competition index from stem maps; Veblen 1986). Trees with evidence of stress (e.g. herbivory from insects, needle dieback), crown damage (e.g., broken top), or other factors affecting cone production were excluded as candidate focal trees in this study (Mueller et al. 2005).
To account for localized tree neighborhood competition, we computed a simplified version of the Neighborhood Competition Index, a density-dependent estimate of competition (Canham et al. 2004). This competition index (CI) for each tree was computed as follows, where dbh i is the dbh of the neighboring tree and dist ji is the distance between the target tree (j) and the midpoint of the concentric circle in which the neighboring tree (i) was rooted (0.5 m for 0-1 m, 1.5 m for 1-3 m, and 4.5 m for 3-6 m). The CI was not scaled by focal tree dbh to reduce collinearity between CI and tree-scale variables. The CI was computed at two spatial scales 0-3 m and 0-6 m to test for the spatial scale that had the greater effect on cone presence/absence and cone abundance.
To model the probability of cone presence (defined above), we constructed binomial generalized linear mixed models (using the function glmer in R, Bates et al. 2015) using all data from temporary and permanent plots for cone presence/absence as a function of tree-scale (dbh, age, species), neighborhood-scale (CI from 0-3 and 0-6 m), and stand-scale (age of oldest tree, BA, and density) predictor variables. We selected plot as a random intercept to account for the hierarchical structure of the data (i.e., trees within stands/ plots). To account for collinearity in predictor variables (Appendix S2: Fig. S1), we compared nested models for each neighborhood-and stand-scale variable and selected the variable at each spatial scale that resulted in the lowest AIC (Appendix S2: Table S1; Bozdogan 1987). Interaction terms between species and the best predictor at each spatial scale were tested and retained only if the interaction term reduced AIC by <2 (Bozdogan 1987). We present two models with either dbh or age, BA (best stand-scale predictor), CI 0-6 m (best neighborhood-scale predictor), and species (P. engelmannii and A. lasiocarpa) to compare the relative importance of predictor variables at multiple spatial scales. All predictor variables were standardized to allow comparison of model coefficients by subtracting the mean and dividing by the standard deviation. Because dbh and age models were not constructed from the same set of observations, we used pseudo r 2 (MuMIn package in R; Barto n 2018) to compare model fit and evaluate expectations.
To assess how cone abundance (estimate of total annual cones produced) is affected by tree-, neighborhood-and stand-scale variability, we constructed zero-inflated generalized linear mixed models (glmmTMB in R; Brooks et al. 2017) with a negative binomial error structure to account for significant zero-inflation (score test for zero-inflation, P < 0.001; van den Broek 1995) and overdispersion found in Poisson models (DHARMa package in R; Hartig 2018). Models were constructed with trees from all temporary plots and a subset of the permanent plots with cone abundance data, because cone abundance data were not available for all trees in the permanent plots (Table 1). In the zero model and count model, we included tree-(age or dbh and species), neighborhood-(CI 0-6 m), and stand-scale (young, low BA or old, high BA; see Table 1) predictor variables as well as year to account for interannual variability in cone abundance. We selected a two-level categorical variable to represent stand structure, because sample size of trees producing >10 cones per plot was low (~10 trees per plot) and stand structure predictor variables within each stand type exhibited similar responses to cone abundance. We included a nested random effect of tree within plot to account for repeated measures and a random effect for year to account for interannual variability in cone abundance. All predictor variables were standardized to allow comparison of model ❖ www.esajournals.org coefficients by subtracting the mean and dividing by the standard deviation. Model comparison of the two colinear tree-scale variables (dbh and age) with the same set of observations was tested using AIC (Appendix S4: Table S1) in models also containing species, CI (0-6 m), and stand type. Interaction terms between species and predictors at each spatial scale were tested using a 2 pt. AIC threshold (Bozdogan 1987). Because our models of cone abundance were created using data from a subset of the plots included in the models of cone presence/absence, we present the zero-inflation model results (i.e. model of cone presence/absence) in Appendix S4. Models constructed using the high and low estimates of cone abundance (see Estimating cone abundance) resulted in similar estimates of model coefficients (AE2% for all predictor variables).

Effect of tree-, neighborhood-, and stand-level factors on probability of cone presence
Tree-and stand-scale variability had the greatest effect on the probability of cone presence, whereas neighborhood-scale competition and species had a minimal effect ( Fig. 2A,B). The probability of cone presence for A. lasiocarpa and P. engelmannii increased with tree size (dbh) and age and decreased with higher competition from neighboring trees (CI) and higher stand BA ( Fig. 2A,B; Appendix S2). P. engelmannii had a slightly higher probability of cone presence than A. lasiocarpa but adding species to the model did not improve model fit (DAIC < 2; Fig. 2; Appendix S2: Table S3), illustrating minimal differences in the age and size of cone presence for P. engelmannii and A. lasiocarpa. Both tree age and size were strong predictors of the probability of cone presence at the tree-scale, but the model including tree size explained more variability in cone presence (74%) than the model including tree age (33%). Though crown volume was not included in the model, the correlation between tree size and crown volume (r s = 0.85, P < 0.01) was stronger than the correlation between tree age and crown volume (r s = 0.46, P < 0.01; Appendix S2: Fig. S1). Competition from more neighbors and/or larger neighbors within 6 m of the focal tree had a greater effect on cone presence than competition within 3 m of the focal tree (DAIC > 2; Appendix S2: Table S1). Comparing stand-scale predictors, stand BA was a stronger predictor of cone presence than stem density or age of oldest tree in the stand (BA lowest AIC; Appendix S3: Table S2). The mean tree size (age in parentheses) at which probability of cone presence was >0.50 ranged from 5 to 7 cm dbh (25-30 yr) in young postfire stands with BA <25 m 2 /ha to 16-18 cm dbh (175-300 yr) in old Picea-Abies stands with BA >65 m 2 /ha (Fig. 2) when the CI was held at the median value of the plot (Fig. 2C,D).

Interannual variability in cone abundance by species and size class
Median cone abundance was greater in main canopy trees (>20 cm dbh) than smaller size classes (4-20 cm dbh), greater for P. engelmannii than A. lasiocarpa, and greater in 2016 and 2017 than 2018 for both species (Fig. 3). During the year of highest cone abundances for both species (2016), median cone abundance for A. lasiocarpa (P. engelmannii in parentheses) was 11 (113) Fig. 3C, F). The percentage of cone-producing trees (at least 1 cone produced during the 2016-2018 study period) with zero cones was on average 21% (mean, range 14-25%) for P. engelmannii and 33% (mean, range 14-59%) for A. lasiocarpa in each year with the highest percentages of zero cone abundance occurring in 2018 (Fig. 3). Over the three-year period of study, >96% of main canopy and intermediate trees for both species produced their maximum cone abundance in 2016 or 2017, and this trend was consistent across the three sampling areas (Fig. 3). However, different individuals contributed to these two consecutive years of higher cone abundance (Appendix S4: Fig. S1), and both A. lasiocarpa and P. engelmannii individuals produced on average <70% fewer cones in the year following maximum cone abundance (2017 or 2018).

Effect of tree-, neighborhood-, and stand-level factors on cone abundance
After initial cone presence (see Results above and zero-model results in Appendix S4), cone abundance was best predicted by tree-and standscale variables (P. engelmannii), whereas ❖ www.esajournals.org neighborhood-scale variables had little effect (Fig. 4). Cone abundance increased with tree size (dbh) and age across both stand types (Fig. 4A-F), but tree size was a stronger predictor of cone abundance than tree age for both species (DAIC = 336; Appendix S4: Table S1). P. engelmannii produced more cones than A. lasiocarpa at similar sizes (species variable, Fig. 4A, C, E). At the neighborhood-scale, cone abundance decreased slightly as the quantity and size of competitors increased (competition index [CI] variable, Fig. 4A, B). At the stand-scale, P. engelmannii cone abundance was considerably lower in older, higher BA stands for trees of similar size or age (species 9 stand type interaction variable, Fig. 4A-F). In contrast, A. lasiocarpa cone abundance was relatively unchanged by stand type (stand type variable; Fig. 4A-F) after initial cone presence, but A. lasiocarpa in young, low BA stands did attain the same size as trees in older, Fig. 2. Probability of cone presence modeling results (i.e., coefficients) from logistic generalized linear mixed models (A, tree size model; B, tree age model) and predictions of cone presence for tree size model (C) and tree age model by stand basal area for Picea engelmannii (PIEN) and Abies lasiocarpa. Models of the probability of cone presence included either tree size (A) or tree age (B), tree species (P. engelmannii and A. lasiocarpa), competition index (CI), and stand basal area. Predictor variables were standardized by standard deviation (SD) and centered on zero to allow comparison among predictor variables. Error bars are 95% parametric bootstrapped confidence intervals (A, B). See Appendix S2 for numerical model results. In (C) and (D), each line represents a stand and lines are shaded by stand basal area (lightest gray is 12 m 2 /ha, and darkest gray is 88 m 2 /ha). Models assume median plot CI and no difference in the threshold dbh or age of initial cone presence for P. engelmannii and A. lasiocarpa.

DISCUSSION
Our study demonstrates that cone presence and cone abundance for P. engelmannii and A. lasiocarpa are affected by tree characteristics (treescale), competition from neighboring trees (neighborhood-scale), and variability in stand structure (stand-scale) in the southern Rocky Mountains. Whereas interannual seed variability within a population is linked to climate cues (Woodward et al. 1994, Buechling et al. 2016, the quantity of cones produced by individuals varies primarily with tree size (for both species) and stand age and BA (P. engelmannii only). By linking variability in cone abundance with tree-to standscale factors, we interpret how seed availability may be affected by tree and stand attributes associated with time since different kinds of disturbances and how seed supply may differ between species in post-disturbance landscapes.
Effect of tree-, neighborhood-, and stand-scale factors on the age/size of reproductive maturity Our finding that tree age and tree size were strong predictors of the timing and size of reproductive maturity (as estimated by the >0.50 probability of cone presence) for P. engelmannii and A. lasiocarpa is consistent with findings in other conifers (Viglas et al. 2013, Davi et al. 2016. In our study, however, tree size was a stronger predictor of reproductive maturity than Fig. 4. Cone abundance model results (i.e., coefficients) for (A) dbh and (B) age zero-inflated, negative binomial generalized linear mixed models for Picea engelmannii and A. lasiocarpa, and predictions of cone abundance for P. engelmannii (C, dbh model; D, age model) and Abies lasiocarpa (E, dbh model; F, age model). Models of cone abundance include dbh or age (tree-scale), species (P. engelmannii, PIEN; A. lasiocarpa, ABLA; tree-scale), competition index (CI; neighborhood-scale), and stand type (young, low basal area or old, high basal area; stand-scale) predictor variables. Predictor variables were standardized to allow comparison, and X indicates an interaction term between predictor variables. Error bars are 95% confidence intervals on coefficient estimates (A, B). Modeled cone abundance (solid line; C, D, E, and F) and prediction error (shading; 95% confidence interval) are marginal effects of dbh or age and stand type, assuming median CI at the stand type level and mean values of random effects in 2016, a high cone abundance year for both species. Observed and predicted cones per tree are plotted on a log scale for interpretability in the graph only (i.e., models were not constructed with log-transformed response variable). Models account for zero-inflation by multiplying the predicted abundance by the predicted probability of presence. Y-axis is different for (C, D) Picea engelmannii and (E, F) A. lasiocarpa. Points represent observed cone abundance values in 2016 (C, D, E, and F). tree age. As cone production requires a significant expenditure of resources (Sala et al. 2012), the transition to reproductive maturity appears to occur when trees attain a size at which they have access to sufficient resources for reproduction, and the age that corresponds to this size can vary considerably for both species. The timing of suppression-release cycles of juveniles in closedcanopy forests can lead to differences of over a century in the age at which P. engelmannii and A. lasiocarpa, and probably many other northern temperate conifers, reach the size of reproductive maturity (Veblen 1986, Wright et al. 2000. Differences in structure among stands (standscale) had a much greater effect on the age and size of reproductive maturity than the spatial clustering of trees within a stand (neighborhoodscale competition) in our study. However, both greater competition and higher stand BA appear to affect the age and size of reproductive maturity. For example, the timing of reproductive maturity occurred >150 yr later and the size of reproductive maturity increased by >10 cm dbh across the range of stand BAs (~20-80 m 2 /ha), which are representative of the stand structural conditions of the focal species (Peet 1981, Aplet et al. 1988. At the stand-scale, we found that stand BA was a stronger predictor than stand density or age of oldest tree (estimate of years since stand initiation). In Picea-Abies forests in the southern Rocky Mountains, stand density can be relatively similar in young, postfire stands and older stands, but the effect of moisture availability on radial growth can result in twofold differences in stand BA with time since stand initiation (e.g., stands BW2 and BW3; Chai et al. 2019). Stand BA may more accurately represent the availability of light and competition for resources than stand density or time since stand initiation for trees approaching the age and size threshold for reproductive maturity (Greene et al. 2002).
P. engelmannii and A. lasiocarpa attained reproductive maturity at roughly the same tree size and age after accounting for variability at the tree-to stand-scale. Abies balsamea (L.) Mill. (balsam fir) can produce cones at smaller sizes than Picea glauca (Moench) Voss (white spruce; Greene et al. 2002), and at the generic level, Abies spp. are considered to be more shade-tolerant than Picea spp (Messier et al. 1999). Therefore, we expected A. lasiocarpa to produce cones at younger ages and at smaller sizes than P. engelmannii, especially in closed-canopy forest. The similarity in the age and size of P. engelmannii and A. lasiocarpa reproductive maturity may have important implications for the durations of their respective periods of cone production. P. engelmannii commonly lives for >600 yr, whereas A. lasiocarpa typically only lives for 350 yr (Oosting andReed 1952, Bigler et al. 2007). In stands with high BA, P. engelmannii and A. lasiocarpa may not begin producing cones until trees arẽ 200 years old, because seedlings can persist in the seedling bank for decades before dying or releasing into the canopy. This could result in significantly longer periods of cone production for the average P. engelmannii than A. lasiocarpa individuals. Given that years of low cone abundance are more common for A. lasiocarpa than P. engelmannii in the southern Rocky Mountains (Alexander 1987), A. lasiocarpa at the individual tree level are unlikely to make up for the shorter period of cone production.

Interannual variability in cone abundance
We observed several-fold differences in cone abundance among individuals within a size class and zero cone abundance for nearly 20% of reproductively mature individuals each year. High variability within and among tree sizes appears to be a common trait among many northern temperate conifers, such as Abies alba (Davi et al. 2016), Abies balsamea (Greene et al. 2002), Picea glauca (Lamontagne and Boutin 2007), and Picea abies (Nygren et al. 2017). Despite high variability in cone abundance, we observed consistency in years of high and low cone abundance across~50 km north to south. The two consecutive years of higher cone abundance (2016, 2017) across all sample areas suggest climatic conditions and pollen supply, two factors known to limit cone production (Pearse et al. 2016), were suitable for cone production during 2016 and 2017 for both P. engelmannii and A. lasiocarpa. However, different individuals contributed to these two consecutive years of higher cone abundance. Resource depletion within trees that produced a high abundance of cones in 2016 may have limited P. engelmannii and A. lasiocarpa individuals from producing abundant cones again in 2017 Rapp 2014, Pearse et al. 2016).

Effect of tree-, neighborhood-, and stand-scale factors on cone abundance
Once P. engelmannii and A. lasiocarpa reached reproductive maturity, cone abundance increased with tree size and age, which is commonly observed for conifers (e.g., Turner et al. 2007, Viglas et al. 2013. However, P. engelmannii and A. lasiocarpa exhibited key differences in the response of cone abundance to stand structure. P. engelmannii in open-canopy forests (young, low BA stands) produced many more cones at smaller tree sizes and younger ages than A. lasiocarpa grown under the same conditions or P. engelmannii in closed-canopy forest. Differences in the effect of stand structure on P. engelmannii and A. lasiocarpa cone abundance may reflect A. lasiocarpa's greater shade tolerance than P. engelmannii (Carter et al. 1988, Brodersen et al. 2006 or the higher outcrossing rates (i.e., lower rates of cone abortion from inbreeding) in P. engelmannii than A. lasiocarpa seed in younger trees (Shea 1987). In the context of expected declines in subalpine forest density from increasing mortality of large trees (van Mantgem et al. 2009, Stovall et al. 2019) and increasingly unsuitable climate conditions for establishment (Andrus et al. 2018a), the capacity for P. engelmannii to produce higher abundances of cones at smaller sizes and younger ages than A. lasiocarpa in open-canopy forests implies that P. engelmannii seed availability may be less limiting than A. lasiocarpa for future seedling establishment.

Potential implications of effects of coarse-scale disturbances on cone production
Results from our study offer insight into initial effects of different types of disturbances on cone production in subalpine forests and may partially explain observed differences in timing and abundance of P. engelmannii and A. lasiocarpa seedling establishment following disturbances. Bark beetle outbreaks and blowdown events in Picea-Abies forests primarily affect large-diameter trees (> 20 cm dbh) within a stand (Veblen et al. 2001, McMillin et al. 2003, Jenkins et al. 2014, implying significant loss of cone production for one or both species following either disturbance. For example, mortality of all large-diameter P. engelmannii (>20 cm; Hart et al. 2014) and many intermediatesized trees (10-20 cm dbh) in old-growth Picea-Abies forests (approximately > 200 yr ;Kulakowski et al. 2016) is a common scenario following spruce beetle (Dendroctonus rufipennis) outbreaks. Such an event would drastically reduce P. engelmannii cone production and the chances for recovery to a species composition similar to the pre-disturbance forest (i.e., compositional resilience), especially as severity and spatial extent of an outbreak increases (Jenkins et al. 2014). Whereas spruce beetles only attack a host species, severe and widespread blowdown events-such as the~10,000-ha Routt-Divide blowdown in northern Colorado-can result in mortality of high percentages of large (>20 cm dbh), dominant seed-producing P. engelmannii and A. lasiocarpa across broad spatial extents (Veblen et al. 2001, Baker et al. 2002. Height growth release of unaffected juveniles, rather than new seedling establishment, is the predominant initial recovery mechanism following bark beetle outbreaks ) and blowdown (Kulakowski and Veblen 2003). Thus, the loss of the dominant cone producers is unlikely to affect the initial recovery, but low abundances of cone-producing trees may significantly limit recruitment of new seedlings for decades, reducing resilience to subsequent disturbance by another type of lethal insect (Temperli et al. 2014). However, species like P. engelmannii may respond to the increased light availability by producing more cones at smaller sizes.
Large, severe (i.e., killing most or all trees in patches >100s of ha) wildfires are characteristic of subalpine forests dominated by P. engelmannii and A. lasiocarpa in the Southern Rocky Mountains , thus resulting in a scarcity of cone-producing trees over large areas. Colonization of P. engelmannii and A. lasiocarpa in stand-replacing patches occurs via wind dispersal from the forest edge and seed availability is a common limitation for postfire regeneration, particularly at distances >100 m from unburned forest (Alexander 1987, Coop et al. 2010. In studies of postfire regeneration, seed availability is commonly assessed as distance to nearest seed source, but our results show high variability in individual tree cone production, suggesting all individuals should not be treated equally in their potential to supply seed. Different responses of P. engelmannii and A. lasiocarpa cone production to stand structure may partially explain consistently greater relative abundance of P. engelmannii than A. lasiocarpa in ~100-yr-old forests recovering from stand-replacing fires and delayed postfire A. lasiocarpa establishment in the southern Rocky Mountains (Veblen 1986, Aplet et al. 1988, Rebertus et al. 1991, Donnegan and Rebertus 1999. In younger, low BA stands (e.g., along fire perimeters or partially burned stands), our study shows that P. engelmannii produce much higher abundances of cones at younger ages than A. lasiocarpa. This would presumably lead to greater P. engelmannii than A. lasiocarpa seed availability along the unburned forest edge to support new seedling establishment within the fire perimeter, and this is supported by higher abundances of spruce than fir seed along the edge of clear-cuts (Noble and Ronco 1978). P. engelmannii seeds likely disperse further into large, high severity patches, because P. engelmannii are generally taller and have lighter weight seeds (Alexander 1987, Shea 1987. Higher probability of P. engelmannii than A. lasiocarpa seedling establishment is further supported by higher outcrossing rates of P. engelmannii (Shea 1987), greater quantities of filled seeds for P. engelmannii (Noble and Ronco 1978), higher germination rates of P. engelmannii (Shea 1987), and the greater suitability of postfire substrates (e.g., bare mineral soil) for P. engelmannii establishment (Knapp and Smith 1982). The capacity for P. engelmannii establishing after fires to produce greater quantities of cones at smaller sizes and younger ages than A. lasiocarpa would also help explain the pattern of earlier and more abundant colonization of severely burned sites by P. engelmannii.

CONCLUSION
Our study contributes to a broader understanding of factors influencing conifer reproduction by explaining variability in two reproductive parameters, cone presence and cone abundance, for P. engelmannii and A. lasiocarpa. Although the greater shade tolerance of A. lasiocarpa compared to P. engelmannii might be expected to yield greater access to resources for reproduction and therefore attainment of reproductive maturity at an earlier age or smaller size, the threshold age and size of reproductive maturity were similar for both species across a large gradient in stand structure. On the other hand, cone abundances in our study were highly variable within and among tree size classes each year for both species, which appears to be a common pattern among many northern temperate conifers. We found that variability in cone abundance was most strongly affected by tree size and that cones per tree for P. engelmannii were considerably greater than A. lasiocarpa in open-canopy forests. The increasing loss of large trees (>20 cm dbh in our study), which are the dominant cone producers, to drought and disturbances (e.g., bark beetle outbreaks) implies that seed availability to support future tree regeneration may be increasingly limited. However, the capacity of some species (e.g., P. engelmannii) to produce higher abundances of cones at smaller sizes and younger ages than other species (e.g., A. lasiocarpa) implies that seed availability may be less limiting for some species under expected declines in abundances of large seed-producing trees.