The distribution of woody species in relation to climate and fire in Yosemite National Park, California, USA

Background: The effects of climate on plant species ranges are well appreciated, but the effects of other processes, such as fire, on plant species distribution are less well understood. We used a dataset of 561 plots 0.1 ha in size located throughout Yosemite National Park, in the Sierra Nevada of California, USA, to determine the joint effects of fire and climate on woody plant species. We analyzed the effect of climate (annual actual evapotranspiration [AET], climatic water deficit [Deficit]) and fire characteristics (occurrence [BURN] for all plots, fire return interval departure [FRID] for unburned plots, and severity of the most severe fire [dNBR]) on the distribution of woody plant species. Results: Of 43 species that were present on at least two plots, 38 species occurred on five or more plots. Of those 38 species, models for the distribution of 13 species (34%) were significantly improved by including the variable for fire occurrence (BURN). Models for the distribution of 10 species (26%) were significantly improved by including FRID, and two species (5%) were improved by including dNBR. Species for which distribution models were improved by inclusion of fire variables included some of the most areally extensive woody plants. Species and ecological zones were aligned along an AET-Deficit gradient from cool and moist to hot and dry conditions. Conclusions: In fire-frequent ecosystems, such as those in most of western North America, species distribution models were improved by including variables related to fire. Models for changing species distributions would also be improved by considering potential changes to the fire regime.


AET:
Annual Actual Evapotranspiration AICc: Akaike's Information Criterion with second order correction BURN: the burned status of a plot Deficit: annual climatic water deficit dNBR: differenced Normalized Burn Ratio FRI: Fire Return Interval FRID: Fire Return Interval Departure MTBS: Monitoring Trends in Burn Severity NBR: Normalized Burn Ratio PRISM: Parameter-elevation Regressions on Independent Slopes Model RdNBR: Relative differenced Normalized Burn Ratio

Background
Since the 1800s, scientists have recognized that the distribution of plant species is influenced by climatic factors. Alexander von Humboldt and Aimé Bonpland mapped the vegetation of two volcanoes in Ecuador in relation to elevation and temperature (von Humboldt and Bonpland 1805). Their diagram of the species distribution on the mountains is considered one of the most important contributions to ecology in the 1800s (Egerton 2009). In 1884, Wladimir Köppen developed his first version of a map of climatic zones based on annual and monthly averages of temperature and precipitation (Köppen and Wegener 1924). These maps were used primarily to delineate potential vegetation-the vegetation that could occur over long periods of absence of exogenous factors at broad spatial scales. The potential vegetation concept underlying these early climatic maps consequentially only included the effects of fire in very limited and indirect ways (e.g., only when fire would have eliminated the species entirely). Merriam (1890Merriam ( , 1898 developed life zone maps for Arizona and the United States based on climate and vegetation, and Grinnell and Storer (1924) mapped the same zones in the Yosemite National Park region. Building on these early studies, Botti (2001) described species distributions using ecological zones determined from a vegetation map developed by the National Park Service in the 1930s as part of a state-wide mapping effort (Wieslander 1935). These maps generally depicted existing vegetation, which may or may not have included the effects of fire.
Although ecologists accepted that fire was a factor controlling the establishment, survival, and mortality of forest species in the late 1800s (Pinchot 1899), it was not considered a principal factor governing the distribution of plant species. Clements (1916Clements ( , 1936 viewed fire as a retrogressive process that set back the orderly progression of vegetation succession from its developmental path to its climatic climax. Gleason (1926) felt that the environment had the strongest influence on plant community development and that fire was an unnatural disturbance that limited the duration of the original vegetation but did not necessarily proscribe its range.
The concept of gradient analysis of vegetation patterns in space and time was introduced by Whittaker (1953Whittaker ( , 1967. He considered fire an environmental factor to which some vegetation was adapted. Without fire, those vegetation types might develop into something completely different. Parker (1989) used a gradient analysis to portray patterns of vegetation in Yosemite National Park. He used elevation and a radiation index as axes to plot community distributions. The index was designed to indicate the effects of slope, aspect, and topographic position on received radiation. However, Parker (1989) did not include fire as an environmental variable in his analysis. Stephenson (1998) described the distribution of woody vegetation in the southern Sierra Nevada according to the annual actual evapotranspiration (AET) and annual unmet water demand (climatic water deficit; Deficit). Especially in regions like the Sierra Nevada with winter snowpack and an extended summer period of low precipitation, the simultaneous availability of water and sunlight are key to plant species presence and abundance. Stephenson (1998) showed that AET and Deficit are more meaningful predictors than elevation, synthetic indices, or average annual temperature and annual precipitation alone. Stephenson (1998) acknowledged the role that fire frequency and intensity could have on species composition but did not include fire in his analysis. Lutz et al. (2010) refined the methods of Stephenson (1998) using the vegetation data collected during the 1930s mapping program and explicit calculations of the effect of slope, aspect, and water capacity on AET and Deficit. Although Lutz et al. (2010) did not include fire effects in their calculations, they did discuss the influence that fire or lack of fire might have on species distributions.
Since the initial development of climate-vegetation maps, it has become apparent that fire is responsible for influencing species composition in fire-adapted systems , for maintaining the abundance of large-diameter trees (Lutz et al. 2009a), for structuring the spatial heterogeneity of forests (Roberts et al. 2008;Kolden et al. 2012;Furniss et al. 2019), and for modifying and maintaining vegetation structure over time (Odion et al. 2010). Frequent fire excludes species that have not adapted their life history traits or reproductive cycles to the fire return interval, and infrequent fire permits species with other adaptations (i.e., adaptations to environmental conditions other than fire) to sometimes achieve dominance over species with specific fire adaptions. What is lacking is an analysis that combines fire effects with climate (but see Lutz et al. 2009b). Such an analysis should include the direct effects of fire as well as the effects of excluding fire, bounding species presence and abundance jointly in terms of climate and the frequency and severity of fire.
The direct effect of fire occurrence can be described simply as a binary variable: burned or unburned (BURN). Within a burned area, however, fires can range from having little effect to severe effects that remove all of the vegetation. Fire severity, defined as the proportionate biomass loss after a fire (van Wagtendonk 2018), can be determined at broad spatial scales by the Landsat-derived differenced normalized burn ratio (dNBR; Key and Benson 2005;Key 2006). Fire severity regimes synthesize the patterns of fire severity over long periods of time, multiple fire events, and numerous ecosystem properties (Sugihara et al. 2018). The effect of climate is integrated into fire severity regimes through its influence on weather, fire behavior, and the underlying vegetation productivity. Based on the relative differenced normalized burn ratio (RdNBR; Miller and Thode 2007), the distributions of fire severity values for each ecological zone in Yosemite National Park were defined by Thode et al. (2011) by combining 32 woody species into 19 different fire regime types. Each regime type represents the response of the dominant vegetation to fire as expressed by fire severity. An explicit linking of climate variables and woody species occurrence to the fire severity regime types for each ecological zone would enhance our understanding of the role of fire in Yosemite National Park.
Fire exclusion can be quantified by the fire return interval departure (FRID; Caprio et al. 1997). FRID is generally presented as the number of years a particular area has gone without fire divided by the median fire return interval for that area. It is a metric for how many fires an area has missed compared to what is expected (historic median) for that vegetation type.
Any vegetation classification using extant plots and (recent) historical fire data must necessarily reflect present climatic conditions and those of the recent past. Classifications should therefore be considered in light of the likelihood that they will change in the near future. It has been established that climate in the montane Sierra Nevada will become both warmer and drier in the twenty-first century (Hayhoe et al. 2004). The resultant increasing climatic water deficit in Yosemite National Park (e.g., Lutz et al. 2010) will lead to declines in species at the drought-limited portions of their ranges. Furthermore, fires are projected to become both larger and burn at higher severity (Lutz et al. 2009b). These fires that burn at higher severity are likely to increase tree mortality to higher levels than predicted by existing calibrations , although there is already considerable uncertainty involved in interpreting data from satellite-derived fire severity metrics (e.g. Furniss et al. 2020). However, because fire is likely to become more prevalent, classifications that consider the effects of fire may be more predictive as climate continues to change.
In this study, we use recently acquired vegetation, soil, climate, and fire data to explain the distribution of woody species in Yosemite National Park and the El Portal Administrative Site (Yosemite). Our objectives were: (1) to determine, for all species on all plots, which combinations of evapotranspiration, climatic water deficit, and burned or unburned status best explain their presence; (2) to determine, for species in plots that did not burn, which combinations of actual evapotranspiration, climatic water deficit, and fire return interval departure best explain their presence; (3) to determine, for species in plots that did burn, which combinations of actual evapotranspiration, climatic water deficit, and fire severity best explain their presence; and (4) to understand how ecological zones represent the combined responses of woody species to climate and fire regimes.

Study area
Yosemite National Park, including the El Portal Administrative Site, is a 302 881 ha reserve in the central Sierra Nevada of California, USA (Fig. 1). Elevations range from 590 m along the Merced River in El Portal to 3997 m at Mount Lyell. Yosemite has a mediterranean climate with hot, dry summers and cool, wet winters. Temperatures range from a mean daily minimum of −1°C in January at the high elevations to a mean daily maximum of 32°C in July at the low elevations. Normal annual precipitation also varies with elevation, from 810 mm at the western boundary to 1200 mm at 2600 m elevation. Most precipitation occurs from October through April, primarily as snow at the mid and high elevations. There is little precipitation during the period from mid June to mid September.
The vegetation in Yosemite National Park occurs in ecological zones that generally correspond to elevation and climate (Fig. 1). These zones were determined by reclassifying a new vegetation map based on association, alliance, and mapping unit descriptions as defined by Keeler-Wolf et al. (2012). At the lowest elevations (about 600 to 1200 m), the vegetation is primarily foothill shrublands and woodlands. Lower montane forests occur between 1200 and 1900 m and consist of a mix of conifer and hardwood species. Upper montane forests between 1900 and 2500 m are primarily conifers, as are subalpine forests (2500 to 3300 m). The alpine zone above 3300 m includes meadows and riparian shrubs. Fire frequency and severity, as well as soil and microclimate, can influence vegetation structure within these zones. Patches of recent high-severity fire promote conversion of lower montane forests to lower montane chaparral, and conversion upper montane conifer forest to montane chaparral, each consisting of a suite of broad-leaved shrubby species .

Data collection Vegetation
We used data on woody species occurrence collected from 561 plots within Yosemite. Yosemite National Park research crews collected data from 158 plots between 1991 and 1993 as part of a natural resource inventory in support of the new vegetation map for Yosemite (Keeler-Wolf et al. 2012). Data from an additional 254 plots were collected in 1998 and 1999 by crews from the Nature Conservancy for the vegetation map (Keeler-Wolf et al. 2012). Between 2002 and 2017, National Park Service fire crews added data from 82 fire and fuel inventory plots designed to monitor fire effects (Lutes et al. 2009). Utah State University contributed data from 67 plots collected in 2011 to study the effects of low-severity fire on species composition changes in the Sierra Nevada (Becker and Lutz 2016). All sample site locations were randomly selected from a landscape stratified by elevation, geology, hydrology, and fire history (burned versus unburned), excluding those areas influenced by roads, trails, and other developments. At the time of establishment, data were collected from 0.1 ha plots that varied in shape from circular to rectangular. Each sample represented a complete inventory of woody species (conifers, hardwoods, and shrubs) within plot boundaries.

Fire
We obtained maximum fire severity data of fires that burned prior to a plot's establishment from data compiled by Lutz et al. (2011) and the Monitoring Trends in Burn Severity (MTBS) program (Eidenshink et al. 2007). Because severity data were not available prior to 1984, we could only determine fire severity for fires burning from that date forward. Fire severity is the magnitude of a fire's effect on the environment and can be quantified using Landsat imagery and the differenced normalized burn ratio (dNBR) developed by Key and Benson (2005). The NBR uses the reflectance values for the near infrared (Band4) and the shortwave infrared (Band7) to determine the ratio between burned and unburned areas, and the dNBR contrasts pre-fire NBR (NBR pre ) with 1-year post-fire NBR (NBR post ): If a plot had a dNBR > 0, we considered it to have burned (BURN). In addition to continuous values for dNBR, the MTBS program provides classified values in four categories: unchanged, 0 to 45; low, 46 to 313; moderate, 314 to 599; and high, 600+.
Yosemite has maintained a fire atlas that records the perimeter of every fire occurring in Yosemite National Park since 1930 (Lutz et al. 2011). For the unburned plots, we used these perimeters to determine when each had burned prior to its initial establishment. Van de Water and Safford (2011) summarized fire frequency estimates for California vegetation types and determined mean, median, mean minimum, and mean maximum fire return intervals for each type. We used the median fire return interval (FRI) for the plot's vegetation type, the establishment year (i.e., the year of measurement) of the plot (YR est ), and the year that it burned last (YR last ) to calculate fire return interval departure (FRID). If the plot had not burned since 1930, YR last was assigned as1930. The equation is: In Yosemite, values for FRID have been classified into three categories: low, 0.0 to 1.9; moderate, 2.0 to 3.9; and high, 4.0+. In the analyses, we used continuous values for FRID and dNBR rather than classified values.

Climate
We obtained mean monthly temperature and precipitation data from 30-year climate normals (1981 through 2010) from the PRISM climate mapping project (Daly et al. 2008). The PRISM data consider meteorological phenomena relevant to mountainous terrain using local regressions between available meteorological stations. These 30 arc-second gridded data (~800 m grid cells) were further subdivided by polygons delineating slope, aspect, and soil water-holding capacity according to the methods of Lutz et al. (2010). Soil water capacity in the top 150 cm of the soil profile was taken from Natural Resources Conservation Service (NRCS 2007) maps. In Yosemite, the NRCS mapped contiguous soil types to a resolution of 0.4 ha in developed areas and to 16 ha in remote areas. We calculated AET and Deficit according to a Thornthwaitetype water-balance model (Thornthwaite 1948;Dingman 2002). Thornthwaite-type methods are appropriate when reliable input data are limited to temperature and precipitation (Vörösmarty et al. 1998;Fisher et al. 2010).

Statistical analysis
We created graphical scatterplots to examine distributions of survey plots with respect to gradients of AET and Deficit, previously established by Stephenson (1998) as being predictive of species occurrence. We measured Pearson correlations to test associations between AET and Deficit among the survey plots and, when the association was significant, we used a reduced major axis regression to identify a primary AET-Deficit gradient along which survey plots tended to occur (Legendre 2018).
Because AET and Deficit co-vary, their effects can be confounded. A positive coefficient for AET or Deficit indicates that a species tends to occur at the higher end of the primary gradient shared by AET and Deficit, even when the other coefficient is not necessarily significant, while a negative coefficient indicates greater occurrence at the lower end of the AET-Deficit gradient. A negative coefficient for AET 2 or Deficit 2 represents a concave likelihood of occurrence along the gradient, with the species tending to concentrate midway along the gradient rather than occurring randomly or at extreme ends of the gradient, while a positive coefficient represents a diffuse distribution of species occurrence. We generally expected species to occupy a niche in AET-Deficit space and thus have negative coefficients for AET 2 and Deficit 2 except in cases for which the evidence did not support such a concentration (e.g., in cases of low sample sizes or diffuse distributions).
We conducted three sets of model comparisons corresponding to our first three objectives of assessing BURN, FRID, and dNBR, respectively. For the first comparison, we modeled species presence across the 561 Table 1 Predictor variables used in models of the distribution of woody species in Yosemite National Park, California, USA (1991 to 2017). Each model includes up to five potential predictors of species presence or absence. Predictor variables consist of actual evapotranspiration (AET), climatic water deficit (Deficit), quadratic effects (AET 2 , Deficit 2 ), a dichotomous variable indicating whether the plot had burned (BURN), fire return interval departure (FRID; for plots that had not burned), and differenced normalized burn ratio (dNBR; for plots that burned plots separately for each species that occurred on at least five plots. Species presence on all plots was coded as a dichotomous response variable (1 for presence, 0 for absence) and analyzed with a logistic regression model using four different combinations of the predictors AET, AET 2 , Deficit, Deficit 2 , and BURN. Among the combinations was a null model containing only AET and Deficit. The other three models were constructed by adding either (1) the pair of quadratic effects for actual evapotranspiration and climatic water deficit (AET 2 , Deficit 2 ), (2) fire occurrence (BURN), or (3) all three (AET 2 , Deficit 2 , BURN). For the second comparison, using a subset of species that occurred on at least five unburned plots, we modeled species presence on the subsample of unburned plots using the same four predictor combinations except replacing BURN with FRID. For the third comparison, using a subset of those species that occurred on at least five plots that burned since 1984, we analyzed species presence on the subsample of burned plots by repeating the analysis except we substituted dNBR for BURN, for a cumulative total of eight predictor combinations (Table 1). We re-centered and rescaled AET and Deficit prior to modeling by subtracting their means (AET = 361, Deficit = 207) and dividing by 100; thus, their model coefficients represent the additive change in the log-transformed odds of species occurrence when increasing AET or Deficit by 100 from their average levels.
For each species, we compared the fit of different models to the presence data using Akaike's Information Criterion with second order correction (AICc; Burnham and Anderson 2002), and the best fitting model (i.e., model with lowest AICc) was selected. We evaluated the significance of predictors in the selected model by examining both their P values in the selected model and their variable importance weights (Stephens et al. 2005). For each comparison set, variable importance weights (Burnham and Anderson 2002) were calculated for each predictor addition: (1) the pair of quadratic effects (AET 2 , Deficit 2 ), and (2) the fire-related metric corresponding to that set (BURN, FRID, or dNBR) as the sum total Akaike weight across models with the predictor addition. The predictor additions for each comparison included the pair of quadratic effects (AET 2 , Deficit 2 ) and the fire-related metric corresponding to that set (BURN, FRID, or dNBR). Because each of these additions occurred in half of the models in the set in which they were compared, an importance weight exceeding 0.5 suggested that models with the addition fit better than models without, with 1 as the maximum importance. Because the linear effects of AET and Deficit were in all models, their importance weights were automatically 1 and are not presented. Model coefficients corresponding to significant predictors represent either positive or negative effects on species occurrence, according to the sign of the coefficient. We classified each predictor as having a significantly positive or negative effect (denoted + or -) when its importance weight >0.7 and coefficient P value <0.05, strongly significant effect (++ or --) when its importance weight >0.9 and coefficient P value <0.01, or no significant effect (0) when its importance weight <0.7 or coefficient P value >0.05. We conducted all statistical analyses using R software (R Core Team 2019).

Results
Out of the 561 plots sampled, 43 woody species were present in two or more plots, and 38 were present in at least five plots and analyzed for effects of having burned ( Table 2). The plots were distributed throughout the five ecological zones in Yosemite National Park (Fig. 1). Ninety-one of the 561 plots have burned since 1984, and 470 were unburned. Out of the 38 species that were present on five or more plots, 11 species were unburned, and eight others burned in fewer than five plots, leaving 19 species that burned in at least five plots to be analyzed for effects of fire severity on a sample of 91 burned plots ( Table 2). All but one of the 38 species that were present in at least five plots were also present in at least five unburned plots, leaving 37 species to be analyzed for effects of fire return interval departure on a sample of 470 unburned plots (Table 2).

Conifers
Conifer species were distributed along a gradient from stands of Pinus albicaulis at low AET and Deficit, to mixed species stands at high AET and moderate Deficit, and then to Pinus sabiniana stands at moderate AET and high Deficit (Fig. 2). When all plots were analyzed (Table 3), five species (Abies concolor, Calocedrus decurrens, Pinus lambertiana, Pinus ponderosa, and Sequoiadendron giganteum) had a positive coefficient for AET, and either a positive or non-significant coefficient for Deficit, corresponding to greater occurrence at the upper end of the AET-Deficit gradient. An additional four species (Abies magnifica, Juniperus grandis, Pinus contorta, and Pinus monticola) had a negative coefficient for AET, and either a negative or non-significant coefficient for Deficit, indicating a tendency to occur at the lower end of the gradient. One species (Pinus jeffreyi) had a negative coefficient for AET but a positive coefficient for Deficit, corresponding to plots that were near the center of the gradient but slightly offset below and to the right of the primary gradient. Pinus sabiniana had a positive coefficient for Deficit but a non-significant coefficient for AET, suggesting occurrence at the higher end of the Deficit gradient. Table 2 Scientific names and authority, common names, species codes, ecological zones, and number of unburned, burned, and total plots for woody species in Yosemite National Park, California, USA (1991 to 2017). The numbers of plots of each species that were analyzed are indicated in bold type. Nomenclature is from Jepson Flora Project (2019)  Two of the remaining three species, Pseudotsuga menziesii and Tsuga mertensiana, had non-significant coefficients for both AET and Deficit, while Pinus albicaulis had a non-significant positive coefficient for AET and a significant negative coefficient for Deficit. Pseudotsuga menziesii is at the southern edge of its continuous distribution in Yosemite and generally occurs along watercourses where the distinct soil variables are of lesser importance. We found negative quadratic relationships (i.e., concave functions) with either AET or Deficit for 11 of 14 conifer species, suggesting that most of these species tended to occur in niches within the AET-Deficit gradient.
The BURN variable occurred in the best models of all conifer species except Pinus contorta, Pinus sabiniana, Pseudotsuga menziesii, and Tsuga mertensiana ( Table 3). The coefficient for BURN was significantly positive for Abies concolor, Abies magnifica, Calocedrus decurrens, Pinus lambertiana, Pinus ponderosa, and Sequoiadendron giganteum. The four species that had non-significant BURN terms in their best fit model were Juniperus grandis, Pinus albicaulis, Pinus jeffreyi, and Pinus monticola.
For those species on the unburned plots with FRID in their best supported model (Table 3), the effect was positive for Abies concolor, Calocedrus decurrens, Pinus lambertiana, Pinus ponderosa, and Pseudotsuga menziesii, suggesting higher occurrences at plots with greater FRID. The effect was negative for Abies magnifica and Pinus contorta, and non-significant for Pinus monticola. Out of the eight species on burned plots (Table 4), only four had dNBR in their best supported models. The effect was positive for Calocedrus decurrens and nonsignificant for Abies magnifica, Pinus contorta, and Sequoiadendron giganteum.
Complete tables of model statistics for the best fit models, importance weights for the predictor variables, and predictor variable coefficients for the conifer species are available from the author on request.
As an example, when all plots were analyzed, there was 100.0% importance for both the combined effect of AET 2 and Deficit 2 , and BURN, in predicting Pinus lambertiana occurrence (model weight 100.0%). Model coefficients were positive for AET and Deficit and negative for AET 2 and Deficit 2 , indicating that Pinus lambertiana tends to concentrate in the upper range but perhaps not at the extreme limits of AET and Deficit. When fire did occur on Pinus lambertiana plots, there was a strong positive response. The absence of dNBR in its best model, however, indicates a lack of evidence that fire severity affected presence. A strongly positive model coefficient for FRID suggests that Pinus lambertiana is more likely to occur on plots that have a long interval without fire. Figure 3 depicts the distribution of plots where Pinus lambertiana was present in relation to AET and Deficit. The negative quadratic effects of AET and Deficit can be seen by the clustering of the plots in the upper range of AET and center range of Deficit, with some plots occurring at higher Deficit values. Additional file 1 contains AET-Deficit plots for all conifer species.

Hardwoods
Hardwood species were distributed partially along the upper range of the AET-Deficit gradient, starting with stands of Populus tremuloides at moderate AET and Deficit, to mixed-species stands at high AET and moderate Deficit (Fig. 4). As Deficit increased, additional hardwood species such as Quercus wislizeni departed from this gradient and tended to occupy plots at moderate AET and high Deficit. Hardwoods were largely absent at the lowest values of AET and Deficit.
Eight hardwood species occurred in more than five plots (Table 5). The most common pattern was reflected by seven species (Alnus rhombifolia, Cercis occidentalis, Cornus nuttallii, Quercus chrysolepis, Quercus kelloggii, Quercus wislizeni, and Umbellularia californica), with positive or non-significant coefficients for at least one of AET or Deficit, suggesting a higher likelihood of Table 2 Scientific names and authority, common names, species codes, ecological zones, and number of unburned, burned, and total plots for woody species in Yosemite National Park, California, USA (1991 to 2017). The numbers of plots of each species that were analyzed are indicated in bold type. Nomenclature is from Jepson Flora Project (2019)  occurrence at the higher end of the AET-Deficit gradients (Table 5). The exception to this pattern was Populus tremuloides, which had a negative coefficient for AET and a non-significant coefficient for Deficit and occurred lower on the gradient. We found only negative or non-significant relationships with AET 2 and Deficit 2 .
Out of the eight hardwood species, four had BURN in their best model. The coefficient was significantly positive for Cornus nuttallii and Quercus kelloggii and nonsignificant for Cercis occidentalis and Umbellularia californica (Table 5). The same eight species occurred on unburned plots (Table 5), and only two of them (Quercus chrysolepis and Quercus kelloggii) had FRID in their best model. In both cases, the effect was significantly positive. Three species (Cornus nuttallii, Quercus chrysolepis, and Quercus kelloggii) occurred on at least five previously burned plots and were modeled with dNBR (Table 6). For these species, only Quercus kelloggii had dNBR in its best fit model, and its coefficient was nonsignificant. Complete tables of model statistics for the best fit models, importance weights for the predictor variables, and predictor variable coefficients for the hardwood species are available from the author on request.
As an example, for hardwood species, Quercus kelloggii had support for both the combined effect of AET 2 and Deficit 2 (weight = 0.9801), and BURN (weight = 0.9677), in predicting its occurrence (model weight 94.9%). Model coefficients were positive for AET and Deficit but negative for Deficit 2 , indicating that Quercus kelloggii tends to concentrate in the upper range of AET and Deficit but not at the extremes of Deficit. The positive model coefficient (± SE) for BURN (0.77 ± 0.26) suggests that Quercus kelloggii is more likely to occur on plots that have burned. Figure 5 depicts the distribution of plots where Quercus kelloggii occurred in relation to AET and Deficit. Additional file 2 contains AET-Deficit plots for all hardwood species. and climatic water deficit (Deficit) are shown for the 394 plots where conifers were present and the 167 plots where conifers were absent. The solid line is the reduced major axis regression that identifies the primary AET-Deficit gradient along which survey plots tended to occur. Mean annual AET and mean annual deficit are plotted as triangles for all species. Species symbols for conifers are listed in Table 2 van Wagtendonk et al. Fire Ecology (2020) 16:22

Shrubs
The plots with shrubs present trended from stands of Salix orestera at low AET and Deficit, to mixed species stands at high AET and moderate Deficit, and then to Ceanothus cuneatus stands at moderate AET and high Deficit (Fig. 6). All shrub species met our criterion of occurrence in at least five plots for modeling (Table 7). Six species (Arctostaphylos patula, Arctostaphylos viscida, Ceanothus cordulatus, Ceanothus integerrimus, Chamaebatia foliolosa, and Prunus emarginata) had a positive coefficient for either AET or Deficit, corresponding to their occurrence at the upper range of the AET-Deficit gradient. One species (Quercus vacciniifolia) had a negative coefficient for AET corresponding to its occurrence in Table 3 Best fit models (i.e., models with lowest AICc) and importance weights for predictor variables for conifer presence in all plots and in all unburned plots in Yosemite National Park, California, USA (1991 to 2017). Species listed occurred in five or more plots for all plots and unburned plots, respectively. Models are defined in Table 1, and slope coefficients are displayed for BURN and FRID when in the best-fit model. The coefficients represent the change in log odds of occurrence for every 10-unit increase in FRID. Important effects are indicated by positive (+) and negative (-) symbols according to the sign of their coefficients for AET or Deficit effects, or by asterisks (*) for BURN and FRID effects. Single sign denotes significance (variable importance >0.7 and P value <0.05), double sign denotes strong significance (variable importance >0.9 and P value <0.01), zero (0) or ns denotes non-significance (variable importance <0.7 or P value >0.05), and n/a denotes that the variable did not occur in the best model Pinus sabiniana 1a 0 n/a ++ n/a n/a 1a 0 n/a ++ n/a n/a  Table 4 Best fit models (i.e., models with lowest AICc) and importance weights for predictor variables for conifer presence in burned plots in Yosemite National Park, California, USA (1991 to 2017). Species listed occurred in five or more burned plots. Models are defined in Table 1, and slope coefficients are displayed for dNBR when in the best-fit model. The coefficients represent the change in log odds of occurrence for every 100-unit increase in dNBR. Important effects are indicated by positive (+) and negative (-) symbols according to the sign of their coefficients for AET or Deficit effects, or by asterisks (*) for dNBR effects. Single sign denotes significance (variable importance >0.7 and P value <0.05), double sign denotes strong significance (variable importance >0.9 and P value <0.01), zero (0) or ns denotes non-significance (variable importance <0.7 or P value >0.05), and n/a denotes that the variable did not occur in the best model the lower half of the AET-Deficit gradient (Fig. 6). The remaining nine species with non-significant coefficients for AET and Deficit occurred near the center of the gradient. Six species (Arctostaphylos patula, Ceanothus cordulatus, Chrysolepis sempervirens, Holodiscus discolor, Prunus emarginata, and Quercus vacciniifolia) had negative coefficients for Deficit 2 and non-significant or negative coefficients for AET 2 , indicating a typical concave likelihood along the AET-Deficit gradient. The remaining species had either a positive coefficient for AET 2 , nonsignificant coefficients for AET 2 and Deficit 2 , or did not have the quadratic effects in their best models. Eight shrub species had BURN in their best fit models, but only Arctostaphylos patula, Ceanothus cordulatus, Ceanothus intergerrimus, Chamaebatia foliolosa, and Salix scouleriana had significant coefficients, and they were all positive. Arctostaphylos viscida, Cercocarpus betuloides, and Holodiscus discolor had non-significant coefficients for BURN. Fire return interval departure (FRID) was significant for Arctostaphylos viscida and Chamaebatia foliolosa, not significant for Quercus vacciniifolia, and not included in the best model for the remaining species (Table 7). Fire severity (dNBR) was significantly positive for Ceanothus cordulatus, Ceanothus integerrimus, and Salix scouleriana, and not included in the best model for the remaining species (Table 8).
Complete tables of model statistics for the best fit models, importance weights for the predictor variables, and predictor variable coefficients for the shrub species are available from the author on request.
The shrub plots are exemplified by Ceanothus cordulatus, which had 100.0% support for the combined effect of AET 2 and Deficit 2 , and BURN, in predicting its occurrence (model weight 100.0%). Model coefficients were positive for AET and Deficit and negative for AET 2 and  Table 2 van Wagtendonk et al. Fire Ecology (2020) 16:22 Deficit 2 . This indicates greater occurrence of Ceanothus cordulatus at higher values of the AET-Deficit but not necessarily at highest AET or Deficit (Fig. 7). The positive coefficient for BURN (1.70 ± 0.30) indicates a strong response of Ceanothus cordulatus to fire. AET-Deficit plots for all shrub species can be found in Additional file 3.

Ecological zones
AET and Deficit were positively correlated among plots (r = 0.57; P < 0.0001), and the distribution of plots within the ecological zones in relation to AET and Deficit fell along a primary gradient from the alpine zone at low AET and Deficit through the subalpine, upper montane, and lower montane zones as AET and Deficit increased (Fig. 8). In a departure from the primary gradient, plots in the foothill zone were located at moderate values of AET and the highest values of Deficit (Table 9). Fire return interval departures (FRID) averaged in the moderate category for the foothill and lower montane plots and low in the upper montane, subalpine, and alpine plots ( Table 9). The average fire severity level (dNBR) was in the low category for all zones that contained plots that burned (Table 9). The highest (901.0) and lowest (3.0) dNBR values occurred in the lower montane zone.

Discussion
For some species, recent fire history (i.e., burned or unburned) has provided additional information about their distribution beyond the climate niche space defined by AET and Deficit. The best models for 13 of the 38 woody species analyzed included significant positive terms for fire. Among these were conifers such as Abies concolor, Pinus lambertiana, and Pinus ponderosa, which had strongly significant positive terms. Stevens et al. 2020 considered the two pines to be resistant to frequent fires, and Thode et al. 2011 placed them in the low and low to moderate fire severity regime types. Abies concolor is a shade-tolerant species and has moderate resistance to fire (Stevens et al. (2020). Fires tend to kill the less fire resistant and more shade-tolerant tree species in the understory of pine stands such as smaller diameter Abies concolor, Calocedrus decurrens, and Pseudotsuga menziesii, thereby perpetuating the overstory species. Without recurrent fire, these species will Table 5 Best fit models (i.e., models with lowest AICc) and importance weights for predictor variables for hardwood presence in all plots and in all unburned plots in Yosemite National Park, California, USA (1991to 2017. Species listed occurred in five or more plots for all plots and unburned plots, respectively. Models are defined in Table 1, and slope coefficients are displayed for BURN and FRID when in the best-fit model. The coefficients represent the change in log odds of occurrence for every 10-unit increase in FRID. Important effects are indicated by positive (+) and negative (-) symbols according to the sign of their coefficients for AET or Deficit effects, or by asterisks (*) for BURN and FRID effects. Single sign denotes significance (variable importance >0.7 and P value <0.05), double sign denotes strong significance (variable importance >0.9 and P value <0.01), zero (0) or ns denotes non-significance (variable importance <0.7 or P value >0.05), and n/a denotes that the variable did not occur in the best model Alnus rhombifolia 1a + n/a 0 n/a n/a 1a + n/a 0 n/a n/a Cercis occidentalis 2a 0 n/a + n/a -18.232 ns 1a 0 n/a + n/a n/a Cornus nuttallii 2b ++ -0 0 0.834 * 1a ++ n/a 0 n/a n/a Populus tremuloides 1b Quercus chrysolepis 1a ++ n/a ++ n/a n/a 3a 0 n/a ++ n/a 1.757 **

Quercus kelloggii
Quercus wislizeni 1a 0 n/a ++ n/a n/a 1a 0 n/a ++ n/a n/a Umbellularia californica 2a 0 n/a ++ n/a -1.854 ns 1a + n/a ++ n/a n/a Table 6 Best fit models (i.e., models with lowest AICc) and importance weights for predictor variables for hardwood presence in burned plots in Yosemite National Park, California, USA (1991 to 2017). Species listed occurred in five or more burned plots. Models are defined in Table 1, and slope coefficients are displayed for dNBR when in the best-fit model. The coefficients represent the change in log odds of occurrence for every 100-unit increase in dNBR. Important effects are indicated by positive (+) and negative (-) symbols according to the sign of their coefficients for AET or Deficit effects. Single sign denotes significance (variable importance >0.7 and P value <0.05), double sign denotes strong significance (variable importance >0.9 and P value <0.01), zero (0) or ns denotes nonsignificance (variable importance <0.7 or P value >0.05), and n/a denotes that the variable did not occur in the best model Quercus kelloggii 4a + n/a 0 n/a 0.167 ns eventually replace less shade-tolerant pines (e.g., Becker and Lutz 2016). Quercus kelloggii was the only hardwood to have a strongly significant positive term for fire occurrence in its best model, and Thode et al. (2011) categorized it as having a moderate fire severity distribution. It is a vigorous sprouter and establishes early after moderate-and high-severity fires. The best models for Arctostaphylos patula, Ceanothus cordulatus, Ceanothus intergerrimus, Chamaebatia foliolosa, and Salix scouleriana also had strongly significant positive fire occurrence terms. Stands of Pinus ponderosa with Arctostaphylos patula, Ceanothus intergerrimus, and Chamaebatia foliolosa in the understory are in the a low-to moderate-fire severity regime type (Thode et al. (2011). On poorer soils, Chamaebatia foliolosa forms dense, low-statured stands underneath Pinus ponderosa. Sprouting in Arctostaphylos patula is fire-stimulated and, if fire removes most of the overstory shading, it would continue to occupy the site. Ceanothus integerrimus and Salix scouleriana are sprouters and their seeds remain viable in the soil for decades to centuries, allowing them to reoccupy severely burned stands of overstory trees.
Most of the species with significant positive terms for fire occurred predominately in the lower montane forest zone. These forests experience heterogeneous fire, and patches of unburned, low severity, moderate severity, and high severity are intermixed based on landscape position (Jeronimo et al. 2019), fuel heterogeneity (Cansler et al. 2019), or the stochasticity of fire behavior (Jeronimo et al. 2020). Species can possess either distinct fire adaptations (Stevens et al. 2020) or can persist in fire refugia (sensu Kolden et al. 2012;Meddens et al. 2018). A notable feature of lower montane forests is fire heterogeneity at all scales, thus permitting some species without direct adaptations to fire to persist in small transient unburned areas (Blomdahl et al. 2019). Species that can, by chance, persist through one or two fires in permanent or transient refugia can attain a stature (i.e., sufficient height to live canopy and bark thickness) that allows them to resist future fires. The nine species that had BURN in their best fit model with nonsignificant coefficients included four conifers, two hardwoods, and three shrubs. In each case, the P value exceeded the threshold for significance (>0.05), and seven had importance value less than the threshold (<0.7). Details about the variety of possible reasons that the effect does not satisfy the variable importance or P value criteria can be found in Additional file 4. The best fit models for species on unburned plots included both positive and negative terms for fire return interval departure. Four conifers (Abies concolor, Calocedrus decurrens, Pinus lambertiana, and Pinus ponderosa), two hardwoods (Quercus chrysolepis and Quercus kelloggii), and one shrub (Chamaebatia foliolosa) had strongly significant positive terms. These species all occur in the lower montane zone and experience frequent, low-intensity fires. Long periods without fire, however, allow seedlings and sprouts to become established and increase the presence of these species. In the upper montane and subalpine zones, Abies magnifica and Pinus contorta had strongly significant positive terms for FRID. After long periods without fire in some mesic areas of the upper montane zone, Abies magnifica can invade and replace stands of Pinus contorta. In the subalpine zone, nearly pure all-aged stands of Pinus contorta are perpetuated over long periods of time without fire. Only one conifer and one shrub had non-significant coefficients for BURN in their best model, and both had P values >0.05 and importance values <0.07 (Additional file 4).
Fire severity, as measured by dNBR, had little or no effect on species presence. Out of the eight species that had dNBR in their best model, only Calocedrus decurrens, Ceanothus cordulatus, Ceanothus intergerrimus, and Salix scouleriana had significant terms, and those terms were positive. Except for large mature trees, Calocedrus decurrens has very thin bark and is very susceptible to mortality from low-intensity fires.  (1991 to 2017), in relation to annual evapotranspiration (AET) and climatic water deficit (Deficit) are shown for the 219 plots where shrubs were present and the 342 plots where they were absent. The solid line is the reduced major axis regression that identifies the primary AET-Deficit gradient along which survey plots tended to occur. Shrub species means are shown as triangles for all species. Species symbols for shrubs are listed in Table 2 van Wagtendonk et al. Fire Ecology Table 7 Best fit models (i.e., models with lowest AICc) and importance weights for predictor variables for shrub species presence in all plots and in all unburned plots in Yosemite National Park, California, USA (1991to 2017. Species listed occurred in five or more plots for all plots and unburned plots, respectively. Models are defined in Table 1 and slope coefficients are displayed for BURN and FRID when in the best-fit model. The coefficients represent the change in log odds of occurrence for every 10-unit increase in FRID. Important effects are indicated by positive (+) and negative (-) symbols according to the sign of their coefficients for AET or Deficit effects, or by asterisks (*) for BURN and FRID effects. Single sign denotes significance (variable importance >0.7 and P value <0.05), double sign denotes strong significance (variable importance >0.9 and P value <0.01), zero (0) or ns denotes non-significance (variable importance <0.7 or P value >0.05), and n/a denotes that the variable did not occur in the best model n/a n/a n/a n/a n/a n/a van Wagtendonk et al. Fire Ecology (2020) 16:22 Table 8 Best fit models (i.e., models with lowest AICc) and importance weights for predictor variables for shrub presence in burned plots in Yosemite National Park, California, USA (1991to 2017. Species listed occurred in five or more burned plots. Models are defined in Table 1, and slope coefficients are displayed for dNBR when in the best-fit model. The coefficients represent the change in log odds of occurrence for every 100-unit increase in dNBR. Important effects are indicated by positive (+) and negative (-) symbols according to the sign of their coefficients for AET or Deficit effects, or by asterisks (*) for dNBR effects. Single sign denotes significance (variable importance >0.7 and P value <0.05), double sign denotes strong significance (variable importance >0.9 and P value <0.01), zero (0) or ns denotes non-significance (variable importance <0.7 or P value >0.05), and n/a denotes that the variable did not occur in the best model Abies magnifica, Pinus contorta, Sequoiadendron giganteum, and Quercus kelloggii had dNBR in their best model but had non-significant coefficients with P values greater than 0.05 or importance values less than 0.7 (Additional file 4). The non-significant coefficients are likely due to non-linear relations with more absences in the no change and low dNBR categories than in the Fig. 8 The distribution of plots in each ecological zone in Yosemite National Park, California, USA (1991 to 2017), in relation to annual evapotranspiration (AET) and climatic water deficit (Deficit). The solid line is the reduced major axis regression that identifies the primary AET-Deficit gradient along which survey plots tended to occur. Mean annual AET and mean annual Deficit are plotted as triangles for plots in each zone Table 9 Number and percent of plots, area and percent of the Park, and the mean ± standard error values for annual deficit (Deficit), annual actual evaporation (AET), fire return interval departure (FRID, and fire severity (dNBR) for each ecological zone in Yosemite National Park, California, USA (1991 to 2017) moderate and high categories. This lack of significance can be attributed to both the relatively low severities resulting from Yosemite National Park's managed fire program, which has been in effect since 1972 (van Wagtendonk and Lutz 2007), and the uncertainties inherent in low-to moderate-severity dNBR analysis (Furniss et al. 2020). The distribution of the probability of occurrence of these species in relation to dNBR is depicted in Additional file 5. Based on our statistical modeling, we feel that fire adaptations and life history traits were the reasons why some species appeared after burns. Most shrub species, as well as Cornus nuttallii and Quercus kelloggii, have their aboveground woody structures killed by fire, but then sprout from still-living rootstocks. Thus, we expect that the presence of genets means that occurrence of these species would not be negatively affected by fire, even though metrics such as cover, stem count, or biomass could be reduced by recent fire (Lutz et al. 2017). These species, having generally lower stature and being more shade intolerant, could be negatively affected by long FRIDs that allow other species to outcompete them. Other species' model responses followed from the environmental settings in which they predominate. Slope shape, slope position, extent of exposed rock, fuel bed characteristics, associated species, and proximity to wetlands influence species response to predictor variables. Short-needled conifers at high elevation form shallow fuel beds with higher fuel moistures than lowerelevation types (van Wagtendonk and Moore 2010). As a result, fire frequency and severity are lower and there is reduced selective pressure for species to develop fireadaptive traits. These species could be little affected by long FRIDs but more affected by elevated severity outside historical ranges.
One consequence of the model formulation is that individual woody stems were modeled without respect to stature. This approach works well for smaller species that reach reproductive stature relatively quickly after establishment. However, for some species (e.g., Pinus lambertiana and Sequoiadendron giganteum), it is the persistence of large-diameter individuals that determines the ability of the population to persist in the presence of fire (e.g., Lutz et al. 2018). Larger diameter individuals of these and similar species generally survive all but the most severe fires. Highseverity fires in lower to upper montane zones generally promote the establishment of shrublands, which, in turn, modifies the frequency and severity of subsequent fires (Airey Lauvaux et al. 2016). The historical fire regime in Yosemite (e.g., Scholl and Taylor 2010) generally featured only small patches of high-severity fire. It is unclear whether the twenty-first century restored fire regime will be similar (e.g., Lydersen et al. 2014;Kane et al. 2015).
These models reflect vegetation as it was measured in Yosemite from the late 1980s to the early 2010s, and also the fire regime characteristic during that time (e.g., van Wagtendonk and Lutz 2007). Although we expect changes to both vegetation and fire in the coming decades, these results should be robust into the near future. For example, although the Rim Fire of 2013 burned with patches of very high severity outside Yosemite (Lydersen et al. 2014), within Yosemite the fire burned with a severity similar to recent, smaller fires (Kane et al. 2015). In the longer term, these classifications may need to be reevaluated as climatic water deficit in the Sierra Nevada increases.
The addition of fire information in the description of ecological zones aids in our understanding of how ecological zones represent the combined responses of woody species to climate and fire regimes. Although there were few data to analyze the terms of fire on the distribution of foothill woodland species, the pronounced terms of fire predictors on lower montane tree species underscores the validity of grouping fire regimes into ecological zones (Thode et al. 2011). Similarly, species occurrence in the upper montane and subalpine zones was largely a function of climatic factors, reflecting the geographic nature of where fire predominates. Species and ecological zones were aligned along an AET-Deficit gradient from cool and moist to hot and dry conditions. In a departure from the primary gradient, species in the foothill zone were located at moderate values of AET and the highest values of Deficit. The distributions along the AET-Deficit gradient mirror the findings of Stephenson (1990) and Lutz et al. (2010) and reinforce the validity of their approach.

Conclusion
Fire can shape the presence and absence of woody plant species as much as climate. Changes in the fire regime (the frequency, severity, and spatial patterns of fire) could consequently change the ranges of plants even when the climatic envelope changes little-or not at all. Our work suggests that, in areas of the world where fire is extensive, such as western North America, potential vegetation distribution models that depict vegetation development over long periods of time independent of fire are less conceptually accurate. Because we expect changes in both vegetation distributions and fire due to changing climate, vegetation models that include both elements will likely be more useful. Thus, with fire being restored to Sierra Nevada forests, it is more accurate to seek models of vegetation presence and absence that explicitly consider the frequency and severity of fire.
Additional file 1. Presence and absence of conifer species in Yosemite National Park, California, USA (1991USA ( to 2017, in relation to annual evapotranspiration (AET) and climatic water deficit. The solid line is the reduced major axis regression that identifies the primary AET-Deficit gradient along which survey plots tended to occur. The plots where a species was present and unburned are shown as blue circles, and the plots that burned are shown as circles colored according to their dNBR values. The thresholds for dNBR are: unchanged, 0 to 45, low, 46 to 313, moderate, 314 to 599, and high 600+. Plots where a species was absent are shown as gray dots. The mean for the unburned plots is shown in the blue triangle and for the burned plots in the triangle colored with the mean dNBR level.
Additional file 2. Presence and absence of hardwood species in Yosemite National Park, California, USA (1991 to 2017), in relation to annual evapotranspiration (AET) and climatic water deficit. The solid line is the reduced major axis regression that identifies the primary AET-Deficit gradient along which survey plots tended to occur. The plots where a species was present and unburned are shown as blue circles, and the plots that burned are shown as circles colored according to their dNBR values. The thresholds for dNBR are: unchanged, 0 to 45, low, 46 to 313, moderate, 314 to 599, and high 600+. Plots where a species was absent are shown as gray dots. The mean for the unburned plots is shown in the blue triangle and for the burned plots in the triangle colored with the mean dNBR level.
Additional file 3. Presence and absence of shrub species in Yosemite National Park, California, USA (1991 to 2017), in relation to annual evapotranspiration (AET) and climatic water deficit. The solid line is the reduced major axis regression that identifies the primary AET-Deficit gradient along which survey plots tended to occur. The plots where a species was present and unburned are shown as blue circles, and the plots that burned are shown as circles colored according to their dNBR values. The thresholds for dNBR are: unchanged, 0 to 45, low, 46 to 313, moderate, 314 to 599, and high 600+. Plots where a species was absent are shown as gray dots. The mean for the unburned plots is shown in the blue triangle and for the burned plots in the triangle colored with the mean dNBR level.
Additional file 4. Best models, number of plots, importance values, ΔAICc values, coefficients, standard errors (SE), and P values for species that had significant values of zero (0) and had models that included BURN for all plots, models that included FRID for unburned plots, and models that included dNBR for burned plots in Yosemite National Park, California, USA (1991 to 2017). The ΔAIC of second best model (ΔAICc of 2nd best model) column is the difference in AIC between the best model that had the fire variable and the second best model that, in most cases except where indicated by asterisks (*), did not have the fire variable. The coefficients represent the change in log odds of occurrence for every 10-unit increase in FRID and 100-unit increase in dNBR.
Additional file 5. Logistic regressions (black lines) of species occurrence in relation to dNBR for those species that had dNBR in their best model but did not have a significant coefficient for dNBR in Yosemite National Park, California, USA (1991 to 2017). The black dots represent the probability of occurrence (presence), and have been jittered to distinguish observations. The vertical lines represent lower thresholds for severity categories: no change 0 (dark green); low, 46 (light green); moderate, 314 (yellow); and high, 600+ (red). The 95% confidence interval is shaded in gray.