Moisture and vegetation cover limit ponderosa pine regeneration in high-severity burn patches in the southwestern US

Fire regimes are shifting in ponderosa pine (Pinus ponderosa Lawson & C. Lawson)-dominated forests, raising concern regarding future vegetation patterns and forest resilience, particularly within high-severity burn patches. The southwestern US has recently experienced a marked increase in large fires that produce large, high-severity patch interiors, with few surviving trees. These areas could be more susceptible for forest loss and conversions to alternative vegetation types than areas closer to the forest edge with more available seed sources. To better understand forest recovery, we surveyed ponderosa pine regeneration within edge and core areas (>200 m from edge) of high-severity patches in ten fires that burned between 1996 to 2008 across Arizona and New Mexico, USA. Specifically, we compared regeneration density, height, and canopy cover in patch edge and core areas and used generalized linear models to investigate the abiotic and biotic factors that contribute to ponderosa pine seedling establishment and density. High-severity burn-patch edge and core plots were not significantly different in seedling density, height, or canopy cover across fires. Seedling establishment was more likely at higher-elevation mesic sites and less likely when Gambel oak (Quercus gambelii Nutt.) was more abundant. Seedling density was negatively impacted by shrub, grass, and Gambel oak cover. Regeneration density varied among fires but analysis of regeneration in aggregated edge and core plots showed that abundance of seed availability was not the sole factor that limited ponderosa pine regeneration, probably because of surviving tree refugia within high-severity burn patches. Furthermore, our findings emphasize that ponderosa pine regeneration in our study area was significantly impacted by xeric topographic environments and vegetation competition. Continued warm and dry conditions and increased wildfire activity may delay the natural recovery of ponderosa pine forests, underscoring the importance of restoration efforts in large, high-severity burn patches.


Introduction
In the last century, fire regimes have shifted drastically in dry forests across the southwestern US, particularly those dominated by ponderosa pine (Pinus ponderosa Lawson & C. Lawson). Prior to Euro-American settlement, frequent, low-severity fire limited fuel accumulations and tree densities and promoted open, unevenaged stands (Swetnam and Baisan 1996;Covington et al. 1997). However, human-induced disruptions to fire regimes (e.g., livestock grazing, fire exclusion) have resulted in changes in forest structure and composition such that they now support anomalous crown fire (Covington and Moore 1994;Fulé et al. 1997;Moore et al. 1999;Agee and Skinner 2005;Savage and Mast 2005). Current forest conditions, in combination with a warmer and drier climate, have given rise to an increase in high-severity fire in Southwestern forests over the past several decades (Singleton et al. 2019;Mueller et al. 2020). Given that ponderosa pines are poorly adapted to such rapid shifts in fire regime and climate, there is heightened concern that forests will not return to prior composition, structure, and function. Furthermore, expected increases in fire activity, fire severity, and temperatures may further limit recovery and threaten the sustainability of ponderosa pine-dominated forest.
One key component of forest resilience (i.e., the ability of a system to absorb disturbance and return to a similar state following disturbance) is the recovery of natural regeneration (Johnstone et al. 2016). Ponderosa pine regeneration is governed by life history traits that facilitate persistence under a regime of frequent, low-severity fire (Pearson 1923). Various adaptations allow mature pines to survive surface fires including thick, insulating bark, self-pruning limbs, long needles that protect buds, and thick bud scales that protect meristem tissue. (Schubert 1974;He et al. 2012;Keeley 2012). Ponderosa pine is an obligate-seeding species, with heavy, small-winged seed that does not facilitate long-distance wind dispersal or rapid colonization of a site after stand-replacing fire. Following severe wildfire, seedling recruitment is heavily dependent upon seed availability from adjacent mature conifers where areas near the edges of high-severity patches (<200 m from edge) tend to have greater seedling densities (Haire and McGarigal 2010;Chambers et al. 2016;Harvey et al. 2016;Kemp et al. 2016;Haffey et al. 2018). However, these adaptive traits are at odds with rapid changes in fire severity and climate, making continued survival of ponderosa pine in post-highseverity fire landscapes uncertain.
Contemporary increases in fire size in the Southwest are associated with large high-severity burn-patch sizes outside the historical range of variation (HRV) for dry forests (Singleton 2020). Over the past three decades, large high-severity patches >100 ha account for 61% of the landscape burned at high severity across all forest types (M. Singleton, Northern Arizona University, Flagstaff, Arizona, USA, unpublished data). In contrast, historical high-severity burn-patch sizes in dry forests were typically less than 1 ha (Pearson 1923) but up to 60 ha (Iniguez et al. 2009). Increasing fire size also corresponds with large high-severity burn-patch interiors >200 m from forest edges (i.e., total core area, Singleton 2020) leaving expansive areas without surviving trees. Since ponderosa pine regeneration is limited to shorter distances from seed trees, it is unclear if pine regeneration will occur within these large, high-severity burn-patch interiors (hereafter, core areas). Under a warming climate, altered fire regimes, and a lack of abundant seed sources, high-severity core areas are predicted to be increasingly susceptible to forest loss and conversion to non-forest states (Walker et al. 2018;Parks et al. 2019;Coop et al. 2020). Therefore, understanding regeneration patterns in large, severely burned areas is critical in managing for forest resilience.
In addition to seed dispersal, successful ponderosa pine seedling establishment also depends on site conditions including biotic (i.e., fungal communities, grass and shrub competition) and abiotic (i.e, topography, slope, moisture environments) factors (Marcolin et al. 2019;Owen et al. 2019). Site conditions in large, highseverity burn patches can potentially vary between edge and core areas, thereby influencing conifer regeneration. For example, recent greenhouse research points to slower growth rates of ponderosa pine in high-severity burn-patch core areas (J. Lippert, Northern Arizona University, Flagstaff, Arizona, USA, unpublished data), suggesting that seedbed conditions may be less favorable or may have less suitable site conditions (i.e., lack of shade, less moisture) for successful seedling regeneration than areas closer to patch edges (but see Owen et al. 2020). Furthermore, Owen et al. (2019) found that highseverity burn-patch core areas have lower density and richness of ectomycorrhizal fungi, which play a role in nutrient uptake and plant growth, compared to edge plots in two Southwestern fires. Since burn-patch edges contain more available seed trees, shade, and potentially different moisture environments and fungal communities, we would expect to see increased seedling densities, more vigorous growth, and decreased abundance of other plant species compared to the potentially harsher environmental conditions within patch core areas. However, more research is needed in assessing differences between burn-patch edge and core areas across many fires in the Southwest, particularly if ponderosa pine regeneration density and growth characteristics (i.e., height, canopy cover) and vegetation abundance differ. Understanding such differences will enable managers to make better informed decisions to understand and manage post-fire landscapes.
Given expected projections of climate and fire activity, it is also important to understand which biotic and abiotic factors contribute to regeneration success and densities in the Southwest. Previous research assessing drivers of post-fire conifer regeneration in high-severity burn patches demonstrated highly variable results between studies (Korb et al. 2019). Thus, it remains unclear if this variation is natural (i.e., geographic) or related to differences in methods. For example, most of these studies were based on a small number of fires (Haire and McGarigal 2010;Owen et al. 2017), on small plot sizes (Haffey et al. 2018), or were conducted outside the Southwest (Bonnet et al. 2005;Dodson and Root 2013;Chambers et al. 2016;Kemp et al. 2016;Rother and Veblen 2016). Since there is potential site variation within burn patches and large geographic variation in factors that drive conifer regeneration between regions in the US (Korb et al. 2019), the goal of this study was to assess ponderosa pine regeneration in high-severity burn-patch edge and core areas and provide a regionally focused assessment of the biotic and abiotic factors that affect ponderosa pine forest recovery across the Southwest. Here we quantified natural regeneration in ten wildfires in northern Arizona and New Mexico, USA, in high-severity areas next to the forest edge and in core areas (>200 m from forest edges) and asked: (1) Do site conditions such as regeneration characteristics (density, height, canopy cover) and vegetation competition differ between edge and core plots? (2) What are the abiotic (i.e., topography and climate) and biotic (i.e., plant competition) drivers of ponderosa pine establishment? (3) Given successful regeneration, what are the drivers of seedling density?

Study area and site selection
We collected field data during the summers of 2018 and 2019 in high-severity burn patches from ten recent fires (1996 to 2008) across northern Arizona and New Mexico (Fig. 1). These fires burned in five national forests and ranged from approximately 1650 to 181 000 ha in area burned, of which 430 to 60 000 ha burned severely. We selected fires and burn patches based on the following criteria: (1) fire must have occurred prior to 2010 to allow for enough time for natural regeneration; (2) fire must have burned in ponderosa pinedominated vegetation type; (3) burn patches must not be reburned or had occurrence of management activities such as planting or thinning; (4) burn patches must contain core area, defined as area exceeding 200 m distance from forest edge; and (5) burn patches must be accessible by foot. The number of sites per fire were based on these criteria. Additional tree species found in the study area were most commonly Gambel oak (Quercus gambelii Nutt. The climate within the study area is warm and semiarid with monsoon rains in late summer (July to September) and winter snow providing the primary sources of moisture for the region (Sheppard et al. 2002). Within fire perimeters, mean winter (January) and summer (July) temperatures ranged from −4.6 to 3.5°C and 13.9 to 23.2°C, respectively. Mean annual precipitation ranged from 362 to 878 mm (Prism Climate Group, https://prism.oregonstate.edu/). Elevation ranged from 1539 to 3173 m as recorded by the lowest and highest pixels within fire perimeters in digital elevation models, and fires spanned approximately 2.5 degrees in latitude (~320 km)

Field data collection
We used satellite-derived relative differenced normalized burn ratio (RdNBR) from the Monitoring Trends in Burn Severity Project and classified burned areas as high severity using a field-calibrated threshold (≥643; Eidenshink et al. 2007;Miller and Thode 2007;Singleton et al. 2019). We used the Southwest Ecological Restoration Unit (ERU) vegetation layer to identify burn patches that were ponderosa pine-dominated forests. The ERU framework identifies pre-fire vegetation by grouping homogenous vegetation classes based on similar site potential and historical disturbance regimes (Wahlberg et al. 2014). Using ArcGIS 10.7 (ESRI, Redlands, California, USA), we generated core areas by buffering classified high-severity 30 m pixels 200 meters inward (Haire and McGarigal 2010). We define a highseverity patch as contiguous area of classified highseverity pixels within the ponderosa pine vegetation class that is large enough to contain core area. Within mapped high-severity burn patches, we established two transects in a random direction (azimuth 0°to 359°) from the centroid of the burn patch. Along each transect, we established two circular 30 m radius plots (0.28 ha), one near the forest edge (≤60 m) and the other within the core area (≥200 m). One transect was established when a burn patch was smaller and its area did not permit two transects. (Table 1). We did not exclude establishing core plots near small groups of surviving trees that were undetected as "unburned islands" during burn severity classification (generally <10 trees), as this represents the heterogeneity of burned areas and small islands are important for regenerating high-severity patch interiors Coop et al. 2020). Plot sizes were sufficiently large to minimize the likelihood of seedling absences, to capture variability in all other biotic variables, and to better match the coarse-scale satellite-derived variables (Table 2).
Within each plot, we tallied the number of tree seedlings >15 cm tall, categorized seedlings into six height size classes (15 to 30 cm, 30 to 50 cm, 50 to 75 cm, 75 to 100 cm, 100 to 137 cm, >137 cm), and recorded tree species. If seedlings were greater than 1.37 m tall, we recorded diameter at breast height (dbh) or diameter at root collar for juniper species. We also established a nested 15 m radius subplot in which we measured the canopy cover for each seedling by taking perpendicular measurements of the long and short side of the canopy.
To quantify vegetation competition, we estimated shrub and grass cover within each plot by establishing four 30 m transects in each cardinal direction from the plot center. Along each transect, we sampled shrubs using the line-intercept method and recorded species. We also used this method to estimate Quercus gambelii Nutt. (hereafter, Gambel oak) cover since stems tended be less than 1.37 m dbh and were most often in dense thickets. We measured Gambel oak separately from shrub cover because of its potential facilitative effect on ponderosa pine regeneration in Southwestern fires (Owen et al. 2017). We measured grass cover by making ocular estimates within two 1 m 2 quadrats along each 30 m transect spaced at 10 m and 20 m (n = 8 per plot). From the plot center, we recorded elevation using a global positioning system and recorded the distance to the nearest living Table 1 Summary information of the ten fires surveyed for ponderosa pine regeneration across Arizona and New Mexico, USA, that burned severely in ponderosa pine forests. Field data were collected in 2018 and 2019. Within each fire perimeter, table includes total area burned in forests and woodlands, percent of area burned severely, the number of plots and transects established in highseverity burn patches, and ponderosa pine regeneration density (mean and median). High-severity fire was mapped using a RdNBR threshold of ≥643 ( mature ponderosa pine using a laser rangefinder. If the nearest tree exceeded the limits of the rangefinder (>500 m) or was obstructed, we used 60 cm National Agriculture Imagery Program (NAIP; www.usgs.gov) aerial imagery to locate and measure the distance to the nearest seed tree. We also recorded slope and aspect for each plot.
To quantify climate in each plot, we used the TerraClimate (Abatzoglou et al. 2018) gridded Table 2 Description of abiotic and biotic predictor variables used in ponderosa pine seedling establishment and seedling density models, using field data that were collected in high-severity burn patches in ten fires across Arizona and New Mexico, USA, in 2018 and 2019. Topographic variables were derived from 30 m resolution digital elevation models (DEM). Climate variables were calculated from 4 km resolution TerraClimate data (Abatzoglou et al. 2018  dataset (4 km resolution) and calculated 30 yr average (1981 to 2010) and 3 yr post-fire average of precipitation, actual evapotranspiration (AET), and climate water deficit (CWD) (Abatzoglou et al. 2018;Rodman et al. 2020a, b). We also calculated a 1 yr post-fire average of precipitation (Welch et al. 2016). Actual evapotranspiration is defined as the evaporative water loss from land surfaces given the supply of water (precipitation and snowmelt ;Stephenson 1998;Abatzoglou et al. 2018). Since AET represents the interaction of available energy and water, it is an index of site productivity (Dobrowski et al. 2015). Climate water deficit is an index of drought stress and represents the evaporative demand of a given site that is not met by the available water supply (Stephenson 1998).
We expected that higher densities of regeneration will correspond with high site productivity (i.e., high values of AET represent high availability of water and energy) and low drought stress (lower values of CWD).
In addition to biotic variables data collected in the field, we also calculated abiotic topographical variables including the slope cosine and sine index (Stage 1976), abs-aspect (Welch et al. 2016), northness, eastness, and curvature (Table 2). Following methods in Rodman et al. (2020b), we used an image-classification approach to quantify post-fire canopy cover surrounding each plot. Specifically, we used the 2017 60 cm 4-band color infrared aerial photography from NAIP and calculated percent canopy cover of mature conifers within 60 to 300 m circular windows, at 30 m increments from each plot center (sensu Rodman et al. 2020b). We performed an accuracy assessment of these methods by generating 950 random points stratified equally in non-forest and forest classified areas. We received an overall accuracy of 91% and a kappa value of 0.83 (Additional file 1)

Statistical analysis Differences between edge and core plots
We compared three regeneration characteristics and three vegetation competition variables between edge and core plots using the one-way Kruskal-Wallis rank sum test (α = 0.05; Kruskal and Wallis 1952) since data were not normally distributed. Specifically, we compared ponderosa pine seedling count, height, and canopy cover, as well as percent grass, Gambel oak, and shrub cover. We calculated ponderosa pine seedling height by taking the mean of each seedling's height class, multiplying it by the number of seedlings within the height class, and summing by plot. To compare canopy cover data, we calculated the area of canopy (i.e., area of an ellipse) for each seedling and summed by plot. For cover estimates, we averaged percent grass cover (n = 8 1 m 2 quadrats) and averaged shrub and Gambel oak cover (n = 4 30 m transects) for each plot.
Logistic regression for seedling presence or absence (hereafter, presence/absence) To determine which factors influenced ponderosa pine regeneration, we fit a binomial generalized mixed model and included fire as a random effect to account for within-fire variation and potential spatial autocorrelation. Generalized linear mixed effect models are an extension of linear mixed models that allow for response variables to have non-normal distributions and allow for random effects (Bolker et al. 2009). We completed all statistical analyses in R version 3.6.1 (R Development Core Team 2019) and developed the logistic mixed model in the lme4 package version 1.1-23 using the glmer function with the family specification set to "binomial" and the link specification set to "logit" (Bates et al. 2015).
To improve model stability and increase the likelihood of model convergence, we scaled all predictor variables prior to model fitting. We used a stepwise model selection process based on minimizing the Akaike Information Criterion (AIC) statistic and selected models with the lowest AIC score (Burnham and Anderson 2004). In cases where competing models had ΔAIC < 2, the most parsimonious model was selected. In addition to minimizing AIC, we also minimized predictor collinearity by calculating the variance inflation factor (VIF) for each candidate model and only retained predictors if the VIF < 3 among terms.
We assessed logistic model fit using the area under curve (AUC) of the Receiver Operating Curve (ROC) in the AUC package version 0.3.0 (Ballings and Van den Poel 2013), which gives the probability that the model correctly classifies observed values. AUC values closer to 1 (AUC values >0.8) indicate excellent model prediction while values near 0.5 indicate that the model is no better than random (Hosmer and Lemeshow 2000). We also assessed model fit by calculating the marginal and conditional R 2 using the performance package version 0.5.1 (Lüdecke et al. 2020). Marginal R 2 describes the proportion of variance explained by fixed effects while conditional R 2 describes the proportion of variance explained by both fixed and random effects (Nakagawa and Schielzeth 2013).

Count models for regeneration abundance
To assess drivers of seedling density, we fit ponderosa pine stems per hectare using a generalized mixed model (GLMM) using the glmmTMB package version 1.0.1 (Brooks et al. 2017). Since conifer density data were right skewed (28% of plots did not contain regeneration) and were overdispersed, we modeled the data with a negative binomial distribution (family = nbinom2). We scaled all predictors and used the buildmer package version 1.6 (Voeten 2020) using the AIC selection criterion to select the best fit model from all possible predictor variables. To mitigate multicollinearity, we excluded correlated variables in the model builder (>|0.70|), but added them independently to candidate models, and made selections that minimized AIC scores (Burnham and Anderson 2004). We did not include fire as a random effect since there was no improvement in AIC score and since the random effect variance was close to zero. Thus, to assess model fit, we calculated only the marginal R 2 . Residual diagnostics, as well as over-dispersion and zero-inflation tests were completed using the DHARMa package version 0.3.2.0 (Hartig 2020).

Differences between edge and core plots
From the 54 plots that we sampled, we counted a total of 4136 conifer seedlings. Of these post-fire coniferous seedlings, 3615 seedlings were ponderosa pine (87%). Although highly variable, we observed ponderosa pine regeneration across all ten fires with a mean regeneration density of 237 (±76 SE) and median density of 19.5 seedlings ha −1 . Mean (median) seedling density was higher in Arizona fires 325 (57) stems ha −1 compared to New Mexico fires 27 (9) stems ha −1 . Mean percent grass cover was the dominant cover type across all fires (18.1%), followed by Gambel oak (15.2%) and then shrub cover (11.9%).
There were no significant differences between edge and core plots in all six regeneration characteristics and vegetation competition variables assessed (Kruskal-Wallis 0.064 ≤ χ 2 ≤ 1.177, df = 1, P ≥ 0.278). Median densities in edge plots were 28.3 stems ha −1 and 10.6 stems ha −1 in core plots. Edge and core plots had similar median cover estimates: edge plots had a median percent shrub, Gambel oak, and grass cover of 6.2%, 14.5%, and 15.1%, respectively, while core plots had 7.4%, 13.2%, and 16.8%, respectively. Median canopy cover in edge plots was 0.29 m 2 compared to 0.05 m 2 in core plots.

Seedling presence/absence
We found that a combination of topographic and biotic variables predicted the presence of ponderosa pine seedlings within high-severity burn patches in the study area. The most important variable predicting seedling presence was the slope cosine index (β = 2.160, P = 0.006; Table 3; Fig. 2A), indicating that the likelihood of pine establishment increased as the interaction of aspect and slope increased (i.e., on steeper northern slopes). Curvature was another important topographic variable in the model and was positively related to pine regeneration, when seedling presence was more likely on convex terrain (β = 1.547, P = 0.021; Table 3; Fig. 2C). Elevation was also a significant topographic predictor in the model and had a positive relationship with seedling establishment, for which increasing elevation had a positive effect on predicting the likelihood of ponderosa pine presence (β = 18.477, P = 0.029; Table 3; Fig. 2D). Percent Gambel oak cover was the only significant biotic variable in the model and had a negative effect on the likelihood of seedling presence (β = −2.315, P = 0.019; Table 3), with a lower probability of seedling presence with increasing Gambel oak cover (Fig. 2B).
The area under the ROC curve was 0.93 which indicates that our model had high predictive performance in discriminating sites that are likely for ponderosa pine regeneration to establish. Marginal R 2 was 0.738 while conditional R 2 was 0.777.

Seedling density
We found that biotic variables largely predicted the abundance of ponderosa pine regeneration in the study area. Percent Gambel oak cover was the only variable that was significant in both the seedling presence/absence model and the seedling abundance model and with the same directionality (i.e., negative coefficient). Percent Gambel oak cover had a negative relationship with regeneration density and was the most important variable in predicting seedling density (P = 0.0001; Table 3), indicating that, as percent Gambel oak cover increased, the predicted density of regeneration decreased (Fig. 3D). Percent grass cover and shrub cover, and distance to nearest living mature conifer were also negatively related to ponderosa pine regeneration densities (P ≤ 0.048; Table 3, Fig. 3B, C). The only important topographic variable in the model was the northness index, which had a positive effect on density, with more northern aspects positively effecting tree densities (P = 0.0005; Fig. 3A). Our model overpredicted seedling density and had a median of 76 seedlings ha −1 compared to the median of 19 seedlings ha −1 for the observed data. However, 72% of the observed values fell within two standard errors of the predicted values and marginal R 2 was 0.862.

Discussion
Following severe wildfire, forest resilience depends on the recovery of tree species. One goal of this study was to investigate ponderosa pine forest recovery within large high-severity core areas >200 m from residual live forest edges, where seedling densities rapidly decline (Haire and McGarigal 2010). To understand differences in recovery patterns within high-severity burn patches, we compared seedling and site characteristics between core and edge areas. We also identified factors that contribute to ponderosa pine establishment and abundance in the Southwest. Our results highlighted three important findings with respect to high-severity burn patches: (1) ponderosa pine seedling densities, height, and canopy cover do not significantly differ between patch core and edge areas across fires; (2) the likelihood of ponderosa pine establishment is largely controlled by topographic variables; and (3) vegetation competition limits seedling abundance.

Are patch interiors more vulnerable?
Our results suggest that ponderosa pine regeneration characteristics as well as competitive abiotic environments were similar across high-severity burn patches. That is, we found that ponderosa pine seedling density, height, and canopy cover did not significantly differ between areas within the interior of high-severity burn patches and areas adjacent to forest edges. Similarly, we found that vegetative competition within high-severity core areas did not significantly differ in percent grass, shrub, and Gambel oak cover between patch edge and core areas. Within high-severity burn patches, many studies have documented that seedling densities significantly decline past 150 to 200 m from mature conifers (Haire and McGarigal 2010;Haffey et al. 2018). Although we documented fewer seedlings in core plots, median densities were similar to edge plots across our study area, thus suggesting that core areas may not be experiencing complete regeneration failures. Our results were likely due to fire "refugia," or unburned islands of surviving trees remaining after wildfire (Krawchuk et al. 2016). Fire refugia are important in providing seed sources and are critical for natural regeneration, especially within large contiguous patches of tree mortality . Within our dataset, 31% of core plots had at least one living mature conifer that was < 150 m from plot center, which highlights the heterogeneity of the burn mosaic within severe fires and patch interiors. Interestingly, within four core plots, we observed seedling densities ranging from 170 622 stems ha −1 at distances 272 to 513 m from the nearest living pine within the Bridgerknoll and Rodeo-Chediski fires. This finding suggests that, despite increasing distances from seed sources, other long-distance dispersal mechanisms, such as animals, may play an important role in regenerating some burn patches (Haire and McGarigal 2010;Owen et al. 2017). In contrast to our findings, Owen et al. (2017) and Owen et al. (2020) found that edge plots within the Pumpkin and Rodeo-Chediski fires in Arizona had significantly more seedling densities than Fig. 3 Prediction probabilities of generalized linear mixed model used to predict ponderosa pine seedling densities in ten fires surveyed in Arizona and New Mexico, USA, in 2018 and 2019. Solid blue lines are the marginal effects of each predictor while other covariates are held at their mean. Blue shaded ribbon represents the 95% confidence interval of predicted values. Ponderosa pine seedling densities are greater on (A) northern aspects, but are limited in areas with greater (B) grass cover, (C) shrub cover, and (D) Gambel oak cover, and (E) at greater distances from nearest seed tree. See Additional file 3 for graphs with observed data core plots, and found no differences in seedling height or other growth characteristics between plot types. These contrasts in findings were likely due to the increase in fires sampled in our study potentially capturing more between-fire variability across the Southwest.
Additionally, we found that seedling canopy cover and height in core areas were not significantly different from those in edge areas. These findings suggest that, since seedling vigor (i.e., canopy cover) and growth rate (i.e., height) are similar, site conditions and resources to assist regeneration success, such as moisture availability or soil nutrients, are also similar. Furthermore, vegetation cover was similar between plot types, suggesting that a similar level of competition was also present. It is important to note, however, that seedling densities and all other site characteristics were highly variable across fires. Although our goal was to assess the effects of high-severity fire on ponderosa pine regeneration across the Southwest, it will be beneficial to assess post-fire recovery within fires for more targeted vulnerability assessments.

Seedling presence/absence
Overall, our model suggested that ponderosa pine seedling presence in high-severity burn patches was most strongly influenced by topographic variables. Similar to other studies, we found that increasing elevation increased the likelihood of seedling presence (Haire and McGarigal 2010;Dodson and Root 2013;Chambers et al. 2016;Rother and Veblen 2016;Haffey et al. 2018;Lopez Ortiz et al. 2019). Environmental gradients, such as increasing elevation, also reflect gradients in resource availability, particularly soil moisture and decreased heat load that are important for seedling establishment (Dodson and Root 2013;Chambers et al. 2016;Kemp et al. 2016). Additionally, we found that increasing slope cosine index and curvature had a positive effect on seedling presence in our logistic regression model. Although areas with greater slopes tend to have dry, shallow soils that may limit regeneration density, a study in inland Northwest USA found that shade-intolerant coniferous species are more likely to occur on rising slopes due to the increased light availability (Shen and Nelson 2018). In our study, the likelihood of ponderosa pine occurrence increased on more northern aspects on greater slopes, highlighting the importance of moisture (north aspects) and micro-climates for seedling establishment (Ziegler et al. 2017). Furthermore, according to our model, the shape of the slope was also important. Topographic curvature is commonly used to identify convex, concave, or flat surfaces to assess the movement of water and subsequent erosion processes. Upwardly convex (positive curvature) surfaces such as hills or ridges indicate areas that decelerate water movement and have divergent water flow, while concave curvature represent depressions and valleys that have convergent flow and accelerate water movement. Seedling establishment may be more probable on convex surfaces because post-fire soil erosion is lessened and, thus, may be more favorable for seed germination. However, these results may also be attributed to an interactive effect of curvature and Gambel oak cover. Gambel oak generally dominated drainage bottoms (i.e., concave terrain) throughout our study area, so it is likely that competition is driving seedling establishment on convex slopes. Given that unburned patches are associated with moister, cooler land surfaces, increased Gambel oak in drainages over time may have important ecological implications for future fire refugia and, in turn, further decrease ponderosa pine regeneration (Dwire and Kauffman 2003).

Seedling density
We found that post-fire ponderosa pine regeneration density in high-severity burn patches across the Southwest was most influenced by vegetation competition. Although the composition of understory varied across the study area, we found that increasing percent grass, shrub, and Gambel oak cover negatively affected ponderosa pine seedling density. Competition between ponderosa pine seedling and coexisting vegetation can vary greatly depending on species and other factors such as soil type, moisture availability, or plant root structure (Abella 2008;Abella and Fulé 2008;Puhlick et al. 2012). For example, grasses quickly develop superficial roots following fire and can deplete soil moisture to outcompete ponderosa pine regeneration, particularly during the pre-monsoon period (Pearson 1942). Gambel oak also has an immediate post-fire competitive advantage over ponderosa pine because it vigorously resprouts following severe crown fire and grows rhizomatous roots that absorb water deep within the soil profile (Abella and Fulé 2008;Kaufmann et al. 2016).
Previous studies have demonstrated similar competitive (i.e., negative) relationships between post-fire conifer regeneration and grasses, Gambel oak, and other resprouting shrubs (Bonnet et al. 2005;Savage and Mast 2005;Haire and McGarigal 2010;Roccaforte et al. 2012;Welch et al. 2016;Rodman et al. 2020b). However, vegetation competition was not a ubiquitous factor that affected conifer regeneration across all studies in the western US, as results from these studies varied both regionally and by fire (see Owen et al. 2017). The Southwest is a region where resprouting shrubs and grasses are rapidly recolonizing sites after high-severity fire (Barton 2002;Savage and Mast 2005), and our results provide further evidence that herbaceous and woody plants are limiting ponderosa pine tree densities.
Unexpectedly, we found that seedling density was not affected by percent canopy cover of surrounding mature conifers, suggesting that higher (or lower) seed availability does not yield higher (or lower) densities, at least within the study region. This may be in part due to the metric we used, whereas metrics that incorporate both proximity and abundance of multiple seed sources (i.e., distance-weighted density) show stronger correlations with tree densities (Haire and McGarigal 2010;Coop et al. 2019;Downing et al. 2019). However, our finding suggests that seed abundance is not the only factor that limits post-fire tree recovery, and that not all mature trees produce abundant seed. For example, Puhlick et al. (2012) found a strong positive correlation of seed tree superiority (i.e., dominant seed tree) and ponderosa pine seedling densities in unburned areas across Arizona, suggesting that the size and productivity of mature trees plays a larger role in regeneration than maturity alone. Future studies could use Puhlick et al. (2012) methods for identifying superior seed trees in high-severity burn patches to better understand their influence compared to seed tree abundance on post-fire recovery.
Although ponderosa pine seed germination is highly dependent on climate, our models did not identify climate as a significant factor in seedling establishment or density. Annual climate conditions have become warmer and drier in the Southwest in recent decades (Mueller et al. 2020). As a result, post-fire environments are increasingly moisture stressed and have less suitable climatic windows for regeneration (Savage et al. 2013). Studies in the Southwest implicate precipitation as an important factor in influencing increased tree densities (Pearson 1923;Cooper 1960;Puhlick et al. 2012) and drought conditions that favor the establishment of dense shrubfields or grasses (Savage et al. 2013). Although there is growing evidence that climatic conditions, particularly drought stress, reduces conifer establishment and densities in montane coniferous forests (Rother and Veblen 2016;Stevens-Rumann et al. 2018;Rodman et al. 2020a, b), climate variables were not an influential driver in our study. This may be in part due to seasonal or annual climate conditions playing a larger role in episodic tree recruitment (Iniguez et al. 2016), which were not thoroughly assessed in this study. Furthermore, studies that document climate factors as strong predictors of tree densities have utilized fine-resolution (250 m) gridded datasets that may better capture smaller topo-climatic variations compared to coarser resolutions (4 km or 800 m; Rodman et al. 2020a, b). Thus, further research is needed in applying downscaled climatic datasets to better assess impacts of climate on seedling regeneration in the Southwest. It is well known, however, that warmer, drier conditions can stress trees and can trigger higher rates of tree mortality or inhibit conifer germination (Breshears et al. 2009;Savage et al. 2013). Under adverse climatic conditions, in conjunction with contemporary shifts in fire regime, Southwestern ponderosa pine recruitment rates may be increasingly constrained.

Management implications and conclusions
Assessing regeneration in high-severity burn patches within ten fires across the US Southwest allowed us to describe general patterns. First, we documented that aggregated core areas of high-severity burn patches are not regenerating differently or have different abundance of vegetation competition than aggregated forest edge areas. Although patch interiors are farther away from a greater abundance of seed sources, seed abundance is not the sole cause of limited regeneration. Importantly, our results suggest that ponderosa pine regeneration in edge and core areas are equally constrained by the abiotic and biotic factors assessed in this study. Second, we found that factors driving both edge-and core-area seedling establishment and densities are largely related to topographic moisture gradients and competition, respectively. Areas that are less xeric (i.e., higher in elevation and on north aspects) were more favorable to successful regeneration, while shrub, grass, and Gambel oak cover were increasing their competitive pressure and limiting tree densities. Although we found strong relationships of predictor variables with regeneration in this study, seedling regeneration may also be due to other factors that were not assessed, such as soil type and conditions, ectomycorrhizal associations, nurse logs, plant associations, and pre-fire basal area. Further research is needed in examining these and other potential factors.
Overall, we observed regeneration in all fires, but abundance was highly variable among edge and core plots across different patches in fires that burned between 1996 and 2008. Overall, we observed that seven out of ten fires had median conifer densities that do not meet the minimum 125 stems ha −1 target that will produce mature trees that are consistent with historical range of variation (HRV) and ecological restoration targets (Reynolds et al. 2013;Ouzts et al. 2015). Although we documented low seedling densities across fires, the seedlings currently found at these sites may be the first of several regeneration pulses over time and, thus, ponderosa pine could potentially recover to its previous density (Haire and McGarigal 2010). However, it remains unclear exactly how much time and the minimum number of seedlings needed to reforest an area to its HRV, especially under contemporary conditions. It is likely that subsequent reburning and increasing competitive pressure from other plant species, along with warmer and drier conditions, will continue to affect mechanisms of ponderosa pine forest recovery and resilience in Southwestern post-high-severity fire landscapes .
Wildfires are becoming larger and more severe in the Southwest, resulting in large patches of high-severity burn areas (Singleton 2020) that can potentially threaten the sustainability of dry, forested systems (Coop et al. 2020). Given this reality, managers need to develop strategies to create more resilient conditions. To this end, we suggest increased use of managed fires during shoulder fire seasons when weather and fuel conditions are milder  and therefore less likely to create more large high-severity burn core areas (Singleton 2020). Managed fires are an effective tool for reducing fuel loads at landscape scales and, thus, reducing the size and severity of future fires (North et al. 2012;Hunter et al. 2014). In addition to avoidance of high-severity burn patches, we also recommend that managers identify and map areas more likely to develop regeneration based on the important topographic variables described in this study. Given that post-fire seedlings are often vulnerable to other disturbances, these maps could then be used to protect reforested areas from future disturbances including future fires that can reset the system and contribute to type-conversion of forested systems (Savage et al. 2013;Coop et al. 2020).
Additional file 1. Results of accuracy assessment of 950 randomly generated points in classified 2017 National Agriculture Imagery Program (NAIP) imagery. Image classification was performed to map forest and non-forest classes to quantify mature conifer canopy cover in ten fires surveyed for post-fire ponderosa pine regeneration across Arizona and New Mexico, USA (sensu Rodman et al. 2020b). Columns indicate class types determined from reference data (i.e., NAIP imagery) and rows indicate class types determined from classified maps. Table includes producer's accuracy (errors of omission) and user's accuracy (errors of commission). The overall map accuracy was 91.4% and the kappa value was 0.83. The kappa statistic reflects the difference between observed accuracy and agreement expected by chance (Congalton and Green 2008). Kappa values close to 1 indicate higher agreement between reference and classified data.
Additional file 2. Prediction probabilities of generalized linear mixed model plotted with observed values. Model was used to predict ponderosa pine seedling establishment (presence/absence) in ten fires surveyed in Arizona and New Mexico, USA, in 2018 and 2019. Solid blue lines are the marginal effects of each predictor while other covariates are held at their mean. Black dots represent observed values, and the blue shaded ribbon represents the 95% confidence interval of predicted values. Ponderosa pine seedling presence was more likely on (A) steep northern aspects, (B) areas with less oak cover, (C) convex terrain, and (D) higher elevations.
Additional file 3. Prediction probabilities of generalized linear mixed model plotted with observed values. Model was used to predict ponderosa pine seedling densities in ten fires surveyed in Arizona and New Mexico, USA, in 2018 and 2019. Solid blue lines are the marginal effects of each predictor while other covariates are held at their mean. Black dots represent observed values, and blue shaded ribbon represents the 95% confidence interval of predicted values. Ponderosa pine seedling densities were greater on (A) northern aspects, but are limited in areas with greater (B) grass cover, (C) shrub cover, and (D) Gambel oak cover, and (E) at greater distances from nearest seed tree. Density values >1000 seedlings ha −1 are not shown to improve graph interpretability.