Climate‐altered fire regimes may increase extirpation risk in an upper subalpine conifer species of management concern

Climate change is transforming forest structure and function by altering the timing, frequency, intensity, and spatial extent of episodic disturbances. Wildland fire regimes in western U.S. coniferous forests are now characterized by longer fire seasons and greater frequency, with further changes expected. Identifying the impacts of altered fire regimes on forest resources may enable land managers to plan mitigation strategies or prepare for novel or altered communities. We created a stochastic, density-dependent, matrix projection model for a whitebark pine (Pinus albicaulis) metapopulation to estimate the impacts of increasing fire frequency on metapopulation persistence. Whitebark pine is a widely distributed foundation species of management concern found in upper subalpine and tree line forests of the Northern Rocky Mountains. We parameterized the model using empirically based demographic data from the Greater Yellowstone Ecosystem (GYE) and validated the model by comparing observed whitebark pine densities to those projected by the model when parameterized with historical demographic rates and fire frequencies. We reparameterized the model with current demographic rates including mortality from insect outbreaks and exotic disease. We compared odds of functional extirpation among six scenarios comprising three altered fire frequencies (fires suppressed, historical fire return interval of 268 yr, and decreasing fire return intervals from current to 97 yr) and two seed dispersal probabilities. Historical parameterization with high dispersal probability projected median whitebark pine densities (40.95 trees/ha, first and third quartiles: 21.89, 67.25), which were similar to empirically estimated densities (40.62 trees/ha, first and third quartiles: 12.04, 114.15). Odds of functional extirpation with increasing fire frequency were 8.26 and 139.91 times higher than historical fire frequency and fire suppression, respectively. In decreasing fire return interval scenarios, odds of functional extirpation were 1.76 times higher in low than high dispersal probability scenarios. These findings suggest that fire suppression may be required to maintain whitebark pine metapopulations in the GYE and that maintaining stand networks connected by high rates of seed dispersal could increase metapopulation resiliency.

Abstract. Climate change is transforming forest structure and function by altering the timing, frequency, intensity, and spatial extent of episodic disturbances. Wildland fire regimes in western U.S. coniferous forests are now characterized by longer fire seasons and greater frequency, with further changes expected. Identifying the impacts of altered fire regimes on forest resources may enable land managers to plan mitigation strategies or prepare for novel or altered communities. We created a stochastic, density-dependent, matrix projection model for a whitebark pine (Pinus albicaulis) metapopulation to estimate the impacts of increasing fire frequency on metapopulation persistence. Whitebark pine is a widely distributed foundation species of management concern found in upper subalpine and tree line forests of the Northern Rocky Mountains. We parameterized the model using empirically based demographic data from the Greater Yellowstone Ecosystem (GYE) and validated the model by comparing observed whitebark pine densities to those projected by the model when parameterized with historical demographic rates and fire frequencies. We reparameterized the model with current demographic rates including mortality from insect outbreaks and exotic disease. We compared odds of functional extirpation among six scenarios comprising three altered fire frequencies (fires suppressed, historical fire return interval of 268 yr, and decreasing fire return intervals from current to 97 yr) and two seed dispersal probabilities. Historical parameterization with high dispersal probability projected median whitebark pine densities (40.95 trees/ha, first and third quartiles: 21.89, 67.25), which were similar to empirically estimated densities (40.62 trees/ha, first and third quartiles: 12.04, 114.15). Odds of functional extirpation with increasing fire frequency were 8.26 and 139.91 times higher than historical fire frequency and fire suppression, respectively. In decreasing fire return interval scenarios, odds of functional extirpation were 1.76 times higher in low than high dispersal probability scenarios. These findings suggest that fire suppression may be required to maintain whitebark pine metapopulations in the GYE and that maintaining stand networks connected by high rates of seed dispersal could increase metapopulation resiliency.

INTRODUCTION
Climate change is expected to transform ecosystem structure and function, not only by altering tree distributions, but also by changing the timing, frequency, intensity, and spatial extent of episodic disturbances such as wildfire and insect outbreaks (Dale et al. 2001, Lenihan et al. 2008, Donato et al. 2009, Westerling et al. 2011, Enright et al. 2015, Buotte et al. 2016, Westerling 2016. In North American coniferous ecosystems, wildland fire is a major form of disturbance that historically renewed seral communities and maintained a mosaic of successional stages on the landscape (Foster et al. 1998, Turner 2010. In recent decades, fire regimes in the U.S. Rocky Mountains have changed in response to warming temperatures and drought; fire seasons are longer, and wildfires have been more frequent and larger (Westerling et al. 2011, Moritz et al. 2012, Dennison et al. 2014, Jolly et al. 2015, Abatzoglou and Williams 2016, Westerling 2016. Since the 1970s, fire seasons have lengthened by nearly three months, annual area burned has increased by almost 3000%, and the number of fires per year has grown by nearly 900% (Westerling 2016). Future climate warming is expected to further increase fire frequency, and thus the time required to burn all areas in a given landscape. In the Greater Yellowstone Ecosystem (GYE), for example, fire rotation time (the time required to burn all area in a given landscape) is expected to decline to fewer than 100 yr and average annual area burned to exceed 1500 ha by the midtwenty-first century (Clark et al. 2017).
Forest ecosystem responses to novel fire regimes are expected to vary; in some cases, changes may disproportionately impact certain species, resulting in altered community composition. In others, conversions to novel forest communities or non-forested ecosystems may occur, impacting ecosystem function and ecosystem services (Enright et al. 2015, Millar andStephenson 2015). Carbon sequestration, energy and water fluxes, nutrient cycling, regulation of runoff, groundwater recharge, disease regulation, timber production, recreation opportunities, and tourism attractions are among the ecosystem services at risk (Costanza et al. 1997, Anderegg et al. 2013 and references therein, Pan et al. 2013, Rocca et al. 2014, IPCC 2018.
Strategic and adaptive forest management will be essential for mitigating the effects of climate change and altered fire and other episodic disturbance regimes (e.g., Canadell andRaupach 2008, IPCC 2018). In fact, some management strategies are shifting from traditional approaches, which aspire to maintain the historical status quo, to more flexible adaptive management strategies that prioritize essential ecosystem services (Millar et al. 2007, Millar and Stephenson 2015, Thom and Seidl 2016). Effective management of forest communities experiencing rapid fire regime changes will require a combination of actions that foster resistance to change, increase system resilience, ensure redundancy, boost representation, and enable adaptation and transition (Millar et al. 2007, Wolf et al. 2015. To identify appropriate management actions, managers require information on the likely trajectories of major forest species under changing fire regimes.
Population-based simulation models are one means to explore population trajectories across a range of scenarios; many have been used for the development of management and recovery plans (Morris et al. 2002, Crone et al. 2011a. For longlived tree species, simulation models may take the place of field studies, which are difficult, if not impossible, to implement because of long generation times. As regeneration failure following fire becomes more common (e.g., Klutsch et al. 2015, Leirfallom et al. 2015, Donato et al. 2016, Stevens-Rumann et al. 2018, it is imperative that we understand consequences for currently vulnerable species.
Whitebark pine (Pinus albicaulis Engelm.) is a long-lived, fire-adapted species of management concern that is widely distributed throughout subalpine and treeline forests of the western United States and Canada (Tomback et al. 2001a, Tomback and. In its Rocky Mountain range, whitebark pine comprises three different but intergrading community types: selfreplacing, treeline, and-on productive upper subalpine sites-successional (Arno 2001). Selfreplacing whitebark pine communities occur on exposed high-elevation sites with poor soils. On these sites, whitebark pine is dominant or codominant, and other occurrence of other conifer species is limited by harsh conditions Hoff 1989, Arno 2001). At treeline, whitebark pine is characterized by dwarfed growth forms (krummholz) and is often the dominant species . In successional stands on productive upper subalpine sites, which are well-represented throughout the GYE and the focus of this paper, whitebark pine is a pioneering species that relies on infrequent high-severity fires characteristic of mixed-severity fire regimes in the GYE (fire return interval 120-400+ yr; Arno 1980, Barrett 1994, Larson et al. 2009) to renew succession and create regeneration opportunities. As the postfire community develops over time, canopy closure reduces regeneration opportunities, and whitebark pine is outcompeted by faster-growing, more shade-tolerant conifers (Arno 1980, Arno and Hoff 1989, Barrett 1994.
Whitebark pine functions ecologically as a keystone and foundation species, and provides a wide range of ecosystem services (Ellison et al. 2005. Its large seeds are a nutritious food source for wildlife, and its high-elevation communities provide important wildlife habitat (Kendall 1983, Mattson et al. 1992, Tomback and Kendall 2001. As a pioneer species following fire on favorable sites, it mitigates harsh postfire conditions and facilitates regeneration of other conifers, thereby fostering community development and biodiversity (Tomback et al. 1993, Tomback 2005, Bansal et al. 2011. Whitebark pine growing at the highest elevations protracts snowmelt and regulates downstream flows, and its roots stabilize soils and prevent erosion (Farnes 1990. A slow-growing, stress-tolerant, and shadeintolerant species (McCune 1988), whitebark pine, rarely reaches reproductive maturity until 30-50 yr of age, in contrast to its forest associates (e.g., Engelmann spruce Picea engelmannii Parry ex Engelm.; subalpine fir; Abies lasiocarpa (Hook) Nutt.; lodgepole pine Pinus contorta Douglas ex Loudon), which can produce seeds at 10-20 yr of age (Krugman and Jenkinson 1974, Lotan and Perry 1983, Alexander et al. 1984. The comparatively late reproductive maturity of whitebark pine leaves it vulnerable to extirpation caused by increased fire frequency (e.g., Enright et al. 2015).
Further, whitebark pine populations in the GYE and elsewhere are faced with disturbance factors in addition to changing fire regimes: the exotic fungal pathogen Cronartium ribicola J.C. Fisher, which causes the often fatal disease white pine blister rust (WPBR), large-scale outbreaks of mountain pine beetle (MPB; Dendroctonus ponderosae Hopkins), and previous fire exclusion practices leading to successional replacement (Arno and Hoff 1989, Arno 2001, Tomback et al. 2001a. Management strategies for whitebark pine populations to prevent and reverse its decline (e.g., Keane et al. 2012) have been developed using a variety of approaches including matrix projection and forest process models, which have considered the demographic influences of WPBR (e.g., Ettl andCottone 2004, Field et al. 2012), MBP (Jules et al. 2016), and climate change (e.g., Ireland et al. 2018). Although some have considered the effects of climate change, few have explicitly considered the effect of altered fire frequency or size on whitebark pine persistence (e.g., Keane et al. 2017). Here, we develop a population-based model to examine whitebark pine ecosystem persistence in the GYE given expected changes to fire regimes and incorporating WPBR and MPB as additional mortality factors. We examine various scenarios to guide potential management mitigation strategies.
Our specific aims were to (1) create a stochastic, density-dependent, matrix projection model simulating a seral whitebark pine metapopulation in the GYE, and use the model to estimate (2) functional extirpation (sapling and mature tree density <0.01 trees/ha; see Population Projections section for rationale) probabilities and median time to functional extirpation under varying fire frequencies, and (3) the effects of varying seed dispersal probabilities (high and low, a proxy for burned areas~20,000 and~100,000 ha in area, respectively) on functional extirpation probability.

METHODS
We created a stage-structured, post-reproduction, stochastic Lefkovitch matrix model (e.g., Caswell 2001) with a one-year time step to estimate the effects of increasing fire frequency and dispersal probability on the functional extirpation probability of a successional whitebark pine metapopulation. The metapopulation is comprised of two distinct populations, each contained within a 1,000-ha area. We then simulated metapopulation size forward in time for 510 yr. We assumed the metapopulation was closed (i.e., no seed input from sources other than the modeled populations) and that populations were separated by unsuitable habitat but connected via bidirectional seed dispersal (e.g., a valley with v www.esajournals.org whitebark pine in the upper subalpine zone on both slopes above the valley floor).
We initially parametrized our model using the historical fire return interval and mature tree survival rate (i.e., not including mortality from WPBR and MPB; see sources in Table 1) to determine whether we could recreate observed whitebark pine densities within the GYE. We then reparametrized the model using recent estimates of mature tree mortality that include the impacts of WPBR and MPB to determine how changing fire regimes will influence whitebark pine's long-term persistence ability, given its current population dynamics (Table 1). To explore metapopulation dynamics, we simulated eight scenarios. We first simulated two historical scenarios differentiated by dispersal probability. Historical simulations were parameterized using the historical fire return interval and low mature tree mortality. High and low dispersal probabilities represented increasing distances between populations. The remaining six scenarios, which we define as recent scenarios, consider mature tree mortality rates including white pine blister rust, three fire return intervals, and two dispersal probabilities. We considered three fire return intervals, (1) historical (μ = 268 yr, σ = 101.8 yr), (2) decreasing, as is expected with climate change (µ =97 yr, σ = 101.8 yr), and (3) suppression (μ = 1000 yr, σ = 750 yr). We considered high and low dispersal probabilities for each fire return interval.

Model description
The model takes the following general form: (1) where n ! f5g ðt þ 1Þ is the stage-structured population vector at time t + 1 comprising population sizes of each of the five life stages included in the model. The matrix Sð θ ! Þ represents the projection matrix comprising survival and transition probabilities, as well as fecundity values (Caswell 2001). We use θ to denote the parameter(s) associated with a particular function. We begin by describing the life stages used in the model, followed by the model structure for one 1,000-ha population, and extensions to incorporate metapopulation dynamics and wildfire.

Whitebark pine life stages
Whitebark pine is a masting species characterized by synchronous and episodic cone production; in the GYE, it usually produces large cone crops once every 3-5 yr (IGBST 2017, Crone et al. 2011a, Sala et al. 2012. Seeds are dispersed up to 32 km by its avian mutualist, Clark's nutcracker (Nucifraga columbiana Wilson). Nutcrackers harvest seeds from whitebark pine cones and store them across the landscape in small caches containing an average of three seeds, and return to retrieve seeds throughout the winter and spring (Tomback 1978, 1982, Hutchins and Lanner 1982. Seed dispersal distances appear to be strongly bimodal, with most caches stored in communal areas and home ranges within a few kilometers of seed sources and a small proportion dispersed over distances greater thañ 25 km, with a median distance of~12 km (Lorenz et al. 2011). Seeds that are unrecovered by nutcrackers and uneaten by granivorous small mammals may germinate (Tomback 1978, 1982, Hutchins and Lanner 1982, Pansing et al. 2017. Cached seeds often germinate during the summer following seed dispersal but some seeds may delay germination, persisting in the soil for two winters or more (McCaughey 1993, Tomback et al. 2001b, Pansing et al. 2017. Seedling mortality is high the year following germination but decreases during the second year (Appendix S1: Table S1). Tree growth is slow, and whitebark pine trees often do not reach reproductive maturity (male and female cones) until 30-50 yr of age, although consistent production of large cone crops is uncommon until 60+ yr (Krugman and Jenkinson 1974). In the GYE, stand-replacing fires that initiate community succession occur, on average, every 268 yr (Arno 1980, Barrett 1994, Larson et al. 2009). For our projection model, we define five whitebark pine life stages: (1) seed caches (CACHE), (2) first-year seedlings (SD1), (3) seedlings (SD), (4) saplings (SAP), and (5) mature adults (MA; Fig. 1). The specifics for each life stage of the model are described below. Sources for survival and transition probabilities are shown in Table 1.
Caches of second-year seeds (CACHE). -We defined an individual in the cache stage to be a seed cache comprised of three seeds (the average number of seeds per cache; Tomback 1978, Hutchins andLanner 1982) that persisted in the v www.esajournals.org soil seed bank for two winters (Fig. 2). We simulated caches rather than seeds contained within caches because multistemmed individuals count as a single tree (e.g., Tomback et al. 2001b). We assumed a CACHE could transition to the SD1 stage (i.e., one or more seeds in the cache germinate) or die.
First-year seedlings (SD1).-The first-year seedling life stage was defined as recently germinated individuals less than one year old. We considered first-year seedlings distinct from older seedlings, because whitebark pine survival rates increase dramatically after the first year (Pansing et al. (2017); D. F. Tomback, unpublished data; Appendix S1: Table S1). We assumed first-year seedlings could survive and transition to the seedling stage (SD) or die.
Seedlings (SD).-We designated seedlings as individuals aged two to 29 yr-the latter the average age at which whitebark pine reach 1.37 m in height (Tomback et al. 1993;D. F. Tomback, unpublished data). Seedlings could (1) survive and remain in the seedling stage, (2) survive and transition to the sapling stage, or (3) die.
Saplings (SAP).-We considered saplings as trees aged 29-50 yr, the average age at which whitebark pine first produces seed cones Mature trees (MA).-We defined mature trees as those more than 40 yr in age and the only individuals that produce seed cones. We incorporated two distinct mortality rates based on the model time period. For historical time periods, we used an average mortality rate consistent with observations of mature pines absent exotic disease or pest infestation (e.g., Franklin 2010, Silver et al. 2013). For recent scenarios, we used estimates from the Greater Yellowstone Monitoring Network that incorporate current mortality from WPBR and MPB  (GRYMN 2018). The fate of mature trees was assumed to be survival or death.

Single population model
The population state vector for a one population model takes the following form: where the elements of n ! f5g depict the number of individuals in each life stage. We use the following projection matrix: where f ð θ ! Þ represents the number of caches created at time t, and is described in the Fecundity section below.
Parameterization of survival and transition probabilities.-Survival probabilities are usually reported as binomial processes where the possible outcomes are survival or death. However, some life stages we considered have three potential outcomes: (1) survive and persist in the same life stage, (2) survive and transition to the next life stage, or (3) die. In these, instances (i.e., SD, SAP stages) survival probabilities were decomposed into the probability of survival and persistence, and the probability of survival and transition.
We used the following equations to decompose the survival and transition probabilities, following Field et al. (2012): where s j,t and t j,t represent the probability of surviving and staying in the same life stage in year t, and surviving and transitioning to the next life stage in year t, respectively, R j represents the life stage-specific residence time (i.e., time an individual persists in any given life stage; Table 1), s * k represents the stochastic survival rate, and k indicates the life stages for which reported survival probabilities needed to be decomposed.
Not all survival and transition probabilities needed to be decomposed, because either the probabilities of every possible outcome for a given life stage have been estimated, or the outcomes were binary (e.g., transition or die). In the former case, for example, mature trees produce caches of seeds that could either (1) germinate, (2) survive (persist as dormant caches), or (3) die. Pansing et al. (2017) estimated the probability of one or more seeds within a cache germinating (t MA-SD1 ), and one or more seeds within a cache surviving as a dormant seed cache (t CACHE ). In the latter case, mature trees could either survive (s MA ) or die. The following probabilities did not require decomposition: Incorporating demographic stochasticity. -We incorporated demographic stochasticity by ran- v www.esajournals.org domly selecting vital rates from probability distributions derived from reported demographic rates (references used are presented in Table 1). We used moment matching to convert means and variances of empirically derived survival and transition probabilities to parameters describing the corresponding beta distribution (α, β), whose domain is defined on the interval (0, 1). This conversion was required because survival probabilities lie on this interval; drawing from other distributions could produce impossible values (i.e., <0 or >1).
At each time step, we selected survival probabilities and transition probabilities not derived from decomposition of survival probabilities (e.g., those for which all possible outcomes been estimated empirically, s * k or t * 1 ) from the beta distribution with parameters defined in Table 1.
The values of s * k were then used as survival and transition probabilities when necessary (i.e., when k = j; Eqs. 4, 5) to create a time step-specific projection matrix.
Density-dependent effects.-Whitebark pine is an early-successional species that is moderately shade-intolerant; germination, seedling survival, and seedling growth rates decline as canopies become more dense Hoff 1989, Perkins 2015). We therefore included density-dependent effects, reducing germination and seedling survival as a function of projected leaf area index (LAI; m 2 /m 2 ). We modified equation 25 used by Field et al. (2012) to include LAI influence from whitebark pine and other conifer species (LAI b ).
We estimated LAI using the following equation: where LA ! T WBP is the row vector comprised of the static leaf area (LA) attributed to an individual whitebark pine in each life stage, area represents the area of a population (1,000 ha), and LAI b represents the LAI contributed by other (heterospecific) conifers.
We defined LAI b as follows: where LA ! T is a row vector specifying the LA of an individual heterospecific conifer in each life stage (defined below), and n ! b is a column vector with elements defining the number of heterospecific conifers in each life stage. We considered seedlings, saplings, and mature trees in calculating LAI b . We did not consider seeds or first-year seedlings because they contribute negligible amounts to LA. We defined the number of heterospecific conifers in a population (n ! b ) as a function of time since fire. We first defined the maximum number of mature trees per 1,000-ha population using data describing late-successional tree densities in the GYE prior to the 1988 Yellowstone fires (x ¼ 931481:5 trees; Tomback et al. 2001b). We assumed progressive mortality would result in more seedlings and saplings than mature trees and defined the maximum number of heterospecific seedlings and saplings as 1.5 and 1.1 times the maximum number of mature trees, respectively. We defined the rate of increase to the maximum number arbitrarily and assumed the number of seedlings increased at a faster rate than the number of saplings, which increased at a faster rate than the number of mature trees. In sum, we defined n b as follows: Relative to recently burned areas, seed germination and seedling survival may be reduced by up to two-thirds in areas with closed canopies (Perkins 2015, Pansing et al. 2017. We therefore reduce germination (r germ ) and seedling survival (r SD ) as a function of LAI (Fig. 3): Lastly, we also incorporate a reduction factor for cone production following the equation used by Field et al. (2012) because cone production can decline with increasing canopy closure (Scott and McCaughey 2006). We start cone availability values at 50% of production because most seed is lost to predispersal seed predation (McKinney and Fiedler 2010).
Incorporation of density-dependent effects modifies the projection matrix as follows: Fecundity.-We simulated fecundity as a combination of seed production and successful dispersal. Because whitebark pine is a masting species, we model cone production using a cosine function with a periodicity of four years Fig. 3. Density-dependent effects are assumed to manifest as reductions in cone production (a), germination (b), and seedling survival (c), as a result of increasing leaf area index (LAI). Proportion of total indicates the amount by which cone production, germination, and seedling survival are reduced given the LAI relative to baseline (i.e., LAI = 0). Cone production begins at 0.5 because we assume half of all seeds are lost to predispersal seed predation.
v www.esajournals.org and normally distributed error (Krugman and Jenkinson 1974, Tomback 1982, Kelly and Sork 2002, McKinney et al. 2009, Crone et al. 2011b): cones t ð Þ ¼ 12:5cos 1:5t ð Þþ14 þ ∈ , ∈ ∼ N 0, 3:5 ð Þ (14) where t is the projection year. The equation was parameterized using the annual average number of cones per tree in transects throughout the Greater Yellowstone Ecosystem from 1980 to 2016 (IGBST 2017). For any model iteration, values lower than zero were set to zero to avoid impossible outcomes. We assumed that all mature trees produce the same number of cones.
The number of cones per tree was used to estimate the number of caches created by nutcrackers in year t. We estimated the number of caches at each time step using Eq. 15, where SpCone = 45 and represents the number of seeds per cone (Tomback 1982), Pcache = 0.7 and is the probability that a seed removed from a cone is cached and not immediately consumed (Keane et al. 1990), Premain = 0.45 (Tomback 1982) and represents the proportion of seeds that are not retrieved from caches by nutcrackers, and SpC = 3 and represents the mean number of seeds per cache (Hutchins andLanner 1982, Tomback 1982), In total, fecundity is defined as

Incorporating metapopulation dynamics
We incorporated metapopulation dynamics for two reasons: (1) assessing the influence of dispersal probability (a proxy for distance to nearest seed source) on functional extirpation probability requires a minimum of two populations, and (2) forest recovery from mixed-severity fires-the fire regimes in upper subalpine environments in the GYE (Morgan et al. 2001)-often begins with seed sources from nearby populations unaffected by fire.
To incorporate metapopulation dynamics, we expanded the population state vector and projection matrix to incorporate the dynamics of the second population and dispersal both between and within populations. Because our metapopulation comprised two distinct populations, the population state vector was comprised of ten elements, five life stages for each population: where stages above the horizontal line above represented population 1, and those below the line represented population 2. We then included parameters allowing dispersal within and between the two populations (Eq. 17; Caswell 2001). We estimated dispersal probabilities using information from Lorenz et al. (2011), which represent the proportion of seeds from population 1 (or 2) that are successfully dispersed to population 2 (or 1). Lorenz et al. (2011) reported median distances traveled by seven nutcrackers dispersing whitebark pine seeds in the Washington Cascade Mountains, the only quantitative data describing distances traveled by nutcrackers to cache seeds. The distribution of distances is strongly bimodal, likely the result of a few nutcrackers dispersing whitebark pine seeds to more distant locations, also observed by Tomback (1978) in the Sierra Nevada Mountain range. We estimated the density of each of the modes of this distribution to define the dispersal v www.esajournals.org probabilities for populations separated by larger distances (>29 km, low probability of dispersal = 0.06) and shorter distances (<10 km, high probability of dispersal = 0.24). Dispersal probabilities remained constant throughout simulations. We assumed all seeds were dispersed to one of the populations (i.e., no seeds are dispersed to unsuitable habitat). Incorporating a second population and dispersal between the population results in Eq. 18, where d represents the probability of dispersal.

Incorporating fire
We considered three different fire return intervals (hereafter, "fire scenarios") to investigate how whitebark pine functional extirpation probabilities may change with different fire frequencies: (1) the historical fire return interval, (2) decreasing fire return interval, and (3) complete fire exclusion.
We determined fire years prior to each iteration of the projection by selecting time between fires from a gamma distribution independently for each population, with parameters defined by the fire scenario. We determined the parameters describing the gamma distribution for each fire scenario by converting the mean and variance using method of moments. Prior to each iteration, we selected values from the gamma distribution and calculated the cumulative sum of the draws to determine the years (t) during which fires would occur: Fire return interval~Gamma μ 2 =σ 2 ,μ=σ 2 À Á where µ represents the mean fire return interval, and σ 2 represents the variance of the fire return interval.
For the historical fire return interval, we estimated the mean and variance of fire return intervals reported for the GYE (Arno 1980, Barrett 1994, Larson et al. 2009). For future changes in fire return internals, we based our distribution on Clark et al. (2017), who suggested that fire rotation times will decrease from current values to fewer than 100 yr. We therefore drew the time between fires from a gamma distribution where the mean fire return interval decreased as a function of time, but variance remained constant as defined in Clark et al. (2017): where µ FRI represents the average fire return interval, t the year, and i the draw number.
We include fire suppression because it is representative of current fire management strategies; 98% of fires in the western United States are contained before reaching 100 ha in size (North et al. 2015). In this scenario, we assume fire suppression efforts are implemented successfully and define a mean fire return interval of 1000 yr, acknowledging that total fire suppression is impossible.
To define whitebark pine mortality during each wildfire event, we scaled mortality rates by the proportion of the population impacted by each fire severity class (high, moderate, low, unburned). This required mortality rate estimates for each fire severity class and the proportional distribution of fire severity classes in whitebark pine forests. We used published mortality rates for each fire severity class (high 0.9-0.95; moderate 0.60-0.89; low 0.2-0.59; unburned 0; Donato et al. 2016, Whittier andGray 2016). We do not consider mortality rates below 0.20 because fire-related mortality can be confounded with background mortality, which we incorporated in the mature tree mortality rate v www.esajournals.org (Whittier and Gray 2016). In each time step during which a fire occurred, we drew a time stepspecific mortality rate for each fire severity class from a uniform distribution with minimum and maximum values defined above (Table 1).
To determine the proportion of the population impacted by each fire severity class, we estimated the proportional fire severity class distribution of all fires in whitebark pine habitat in the GYE between 1984 and 2018 using data from the Monitoring Trends in Burn Severity database (www.mtbs.gov). MTBS maps the extent and severity of all U.S. wildfires larger than 404 ha that burned since 1984. We estimated fire severity using Relative differenced Normalized Burn Ratio (RdNBR), which correlates strongly with fieldbased metrics of fire severity Thode 2007, Parks et al. 2014). RdNBR accounts for changes in prefire vegetation among locations and is therefore better suited for regional comparisons over other commonly used metrics including NBR and dNBR (Miller and Thode 2007, Cansler and Mckenzie 2014, Parks et al. 2014). We downloaded RdNBR layers from MTBS for all fires in the GYE and masked each fire layer to include only areas that burned in whitebark pine habitat (as defined by the most recent whitebark pine range map; WPEF 2019). For each fire, we determined the proportion impacted by high, moderate, and low severity fire, as well as the proportion unburned within the fire perimeter, using thresholds defined by Miller and Thode (2007). A total of 246 wildfires >404 ha occurred in the GYE during the time period 1984-2017, of which 85 burned in whitebark pine habitat. The mean proportion affected by each fire severity class was 0.40, 0.21, 0.20, and 0.19 for high, moderate, low, and unburned classes, respectively (Appendix S1: Fig. S1). During model time steps when a fire occurred, we sampled randomly one fire from our population of 85 and assigned fire severity classes in the same proportion as the selected wildfire. We then calculated the total proportion of trees killed (p mort ) by fire as p mort ¼ ðp high severity Â p mortality in high severity Þ þðp moderate severity Â p mortality in moderate severity Þ þ p low severity Â p mortality in low severity (20) where p indicates a proportion. We assumed all seedlings and saplings were killed.
Although we acknowledge that there are documented associations between fire size and proportion of high-severity fire (e.g., Mckenzie 2014, Harvey et al. 2016b), we do not increase the proportion of high-severity fire with climate change, because the rate of change is still unresolved (e.g., Harvey et al. 2016a, Parks et al. 2018).

Historical and current time period simulations
To model historical conditions, we used an average mature tree mortality rate of 0.01 (survival = 0.99), which is characteristic of mature trees in the absence of disturbance (e.g., Gill 2005, Larson andFranklin 2010). For current time periods, we increased mature tree mortality from 0.01 to 0.13 (GRYMN 2018) to incorporate dynamics driven by WPBR and MPB. We focused on the mature life stage because WPBR infections were absent in samples used to estimate first-year (SD1) and seedling (SD) survival, and MPB rarely attack small diameter individuals (i.e., saplings) and do not attack seedlings (Shanahan et al. 2016, Cole andAmman 1969). We assumed that MPB was present only at endemic levels in the historical time period scenarios to improve interpretation of results, although in reality intermittent outbreaks occur in many regions in the west (Perkins and Swetnam 1996, Gibson et al. 2008). We focus on comparing outcomes among different fire scenarios for recent conditions, because we are interested in assessing the possible impact of decreasing fire return intervals on whitebark pine given current population dynamics.

Population projections
We projected whitebark pine metapopulation size for 510 yr for 5000 iterations for each of eight scenarios (Table 2). We began each model iteration with a starting population state vector ðn ! 10 f g Þ sampled randomly from a population of empirically estimated densities (FIA 2019, Tomback et al. 2001;D. F. Tomback, unpublished data). We estimated seedling, sapling, and mature tree densities using Forest Inventory and Analysis (FIA) data and first-year seedling densities using data from D. F. Tomback et al. (unpublished data).
The FIA program, run by the USDA Forest Service, probabilistically samples all forested land in the United States (FIA 2019). Data collected include the number of trees per species (DBH ≥ 5.0 cm) on each of three 0.0169-ha subplots, and the number of saplings (1.0 cm ≤ DBH < 5.0 cm) and seedlings (height > 15.24 cm and DBH < 1 cm) in each of three 0.00135-ha microplots. We downloaded data for Wyoming, Idaho, and Montana from FIA DataMart and subsequently excluded plots outside of the GYE (FIA 2019). We further restricted plots to those in which whitebark pine was documented. We then estimated seedling, sapling, and mature tree (DBH ≥ 5.0 cm) densities for each plot and year by averaging densities from the nested subplots (mature trees) and microplots (saplings and seedlings) and converted the densities to number of trees per 1000 ha. A total of 1108 plots contained whitebark pine trees.
Because FIA does not count seedlings <15.24 cm in height, we derived densities of first-year seedlings using a 30-yr regeneration dataset following the fate of individual whitebark pine seedlings (see Tomback et al. 2001 for study design). These data were used to estimate the density and, subsequently, count per 1000 ha for each plot and survey year.
At the beginning of each model iteration, we randomly selected two FIA plots and assigned the number of seedlings, saplings, and mature trees from each plot to the elements comprising one population of the state vector. We similarly randomly selected two unique plot/year combinations from the 30-yr dataset to assign the starting number of first-year seedlings in each subpopulation. Lastly, we determined the starting number of caches using Eq. 15, parameterized with the number of mature trees in the starting population state vector. We use Eq. 15 because we are unaware of empirical estimates of seed cache density. Initial time since fire was drawn randomly from a uniform distribution with parameters 0 and 268. To account for potential lack of correlation between density and time since fire that resulted, we eliminated the first 10 time steps of each iteration. For each scenario and iteration, we tracked annually the population size for every life stage in each population and for the metapopulation as a whole, time since the most recent fire, and LAI. We converted the number of individuals in the sapling and mature tree life stages to densities (no. trees/ha) to allow for comparison to densities estimated using FIA plot data from the GYE (FIA 2019).
To estimate functional extirpation probabilities for each scenario, we calculated the proportion of iterations in which density dropped below 0.01 saplings and mature trees/ha for the metapopulation. We chose this threshold based on the probability of nutcracker occurrence in whitebark pine forests, assuming that if nutcrackers do not visit, whitebark pine populations are functionally extirpated, because seeds will not be dispersed. Barringer et al. (2012) found that nutcracker occurrence is associated with the percentage of live trees, and the minimum stand basal area in which nutcrackers were absent was 0.04 m 2 /ha. This basal area corresponds to a mature tree density (as defined by this study) of~3.5 trees/ha. We chose a functional extirpation density threshold more than an order of magnitude lower than 3.5 trees/ha to ensure tree densities are unlikely to result in successful seed dispersal. For each scenario parameterized with the recent mature tree survival rate (i.e., recent time periods), we estimated median time to functional extirpation and the cumulative proportion of model iterations that fell below the functional extirpation threshold across the 500-yr projection period. Subsequently, we estimated odds ratios to assess the effect of changing fire frequency and dispersal probability on functional extirpation (Rita and Komonen 2008). We first estimated the odds ratio comparing functional extirpation in the historical and recent scenarios. Subsequently, for scenarios with recent mature tree mortality, we estimated odds ratios comparing functional extirpation in the decreasing fire return intervals, historical fire return intervals, and suppression, as well as high and low dispersal probabilities within similar fire scenarios. We do not assess outcomes for different fire return intervals for our historical scenarios, because we were interested in assessing the possible impact of decreasing fire return intervals on whitebark pine populations currently on the landscape. The model was built, and all analyses were conducted, using the R programming language (R Core Team 2018; www.github.com/epansing/ WBPdemographicModel).

Historical conditions
Parametrized with historical demographics and fire frequencies (low rates of mature tree mortality and historical fire return interval), we projected a final (t = 500) median whitebark pine sapling and mature tree density of 40.96 trees/ha (first and third quartiles: 21.90, 67.25) and 1.99 trees/ha (first and third quartiles: 0.85, 4.21) for the high and low dispersal probability scenarios, respectively (Fig. 4). Functional extirpation occurred in 2.04% and 2.98% of model iterations in the high and low dispersal scenarios, respectively. Of those iterations that fell below the functional extirpation threshold, the median time to functional extirpation was 371 yr (first and third quartiles: 219, 394) in the high dispersal scenario and 341 yr (first and third quartiles: 218, 424) years in the low dispersal scenario.
Projected densities across the postfire successional trajectory were similar to whitebark pine densities we estimated from FIA plots in the GYE and other published densities , FIA 2019D. F. Tomback, unpublished data): In late-successional forests (longer than the average fire return interval, ≥268 yr postfire), we projected a median sapling and mature tree densities of 62.91 trees/ha (first and third quartiles: 38.71, 95.60) in the high dispersal scenario and 11.50 trees/ha (first and third quartiles: 6.64, 18.48) in the low dispersal scenario, consistent with median density estimated from the GYE FIA monitoring plots used for state vector initialization of 40.62 trees/ha (first and third quartiles: 12.04, 114.15; Fig. 5). During fire recovery (30 yr postfire), we projected a median tree density (first-year seedlings, seedlings, saplings, and mature trees) of 28.41 trees/ha (first and third quartiles: 14.74, 50.00 trees/ha) and 0.78 trees/ha (first and third quartiles: 0.28, 1.89 trees/ha) for high and low dispersal scenarios, respectively. These estimates are consistent with median densities observed in areas impacted by stand-replacing fire during the 1988 Yellowstone fires and adjacent unburned forest following 30 yr of recovery (median = 0 trees/ha; first and third quartiles: 0.0, 500.00 trees/ha; D. F. Tomback, unpublished data).

Recent conditions
The six scenarios incorporating recent population health challenges (i.e., mature tree mortality from WPBR and MPB) all had projected declining median densities across the time horizon (Fig. 6). Final (t = 500) median sapling and mature tree density from high to low were produced by the following scenarios in declining order: (1) The historical fire return interval and high dispersal probability scenario produced the highest median sapling and mature tree density (40.73 trees/ha; first and third quartiles: 22.05, 68.48); the (2) high and then (3) low dispersal probability suppression scenarios (38.36 trees/ha; first and third quartiles: 23.84, 59.95 and 31.95 trees/ha; first and third quartiles: 18.65, 50.32, respectively); (4) the high dispersal decreasing fire return interval (19.39 trees/ha; first and third quartiles: 5.08, 45.77); (5) the low dispersal probability historical fire return interval (1.95 trees/ha; first and third quartiles: 0.81, 4.19), and (6) the low dispersal probability decreasing fire return interval (0.31 trees/ha; first and third quartiles: 0.11, 1.04; Fig. 7).
In the decreasing fire return interval scenarios, 22.6% and 14.3% of model iterations dropped below the functional extirpation threshold in the low dispersal and high dispersal scenarios, respectively, by year 500 (Fig. 8). The median number of years until functional extirpation was 438 yr (first and third quartiles: 317, 453) and 329 yr (first and third quartiles: 188, 445; Table 3) for the low and high dispersal scenarios, respectively. In the historical fire return interval, 2.66% and 2.34% of the model iterations fell below the functional extirpation threshold in the low and high dispersal scenarios by year 500 (Fig. 8). The median number of years before functional extirpation occurred was 306 (first and third quartiles: 232, 389) for the low dispersal scenario and 297 (first and third quartiles: 220, 409) for the high dispersal scenario (Table 3). The low and high  dispersal scenarios for the fire suppression fell below the functional extirpation threshold in 0.16% and 0.14% of iterations, respectively. The median time to functional extirpation was 264 yr (first and third quartiles: 241, 379) for the low dispersal scenario and 245 yr (first and third quartiles: 193, 390) for the high dispersal scenario.

Comparing odds of functional extirpation among scenarios
Higher mature tree mortality rates (i.e., recent time periods) resulted in similar odds of functional extirpation to historical time periods (OR, 1.00; 95% CI: 0.837, 1.20) with historical fire return intervals. When considering recent scenarios only, the odds of functional extirpation were 10.7 (95% CI: 8.87, 13.0) and 6.94 (95% CI: 5.68, 8.55) times higher in the decreasing fire return interval scenario than the historical fire return interval in the low and high dispersal probability scenarios, respectively. The odds of functional extirpation in the decreasing fire return interval scenarios were 182.12 (95% CI: 91.87, 419.08) and 118.58 (95% CI: 57.13, 295.06) times higher than fire suppression scenarios in the low and high dispersal scenarios, respectively. Lastly, odds of functional extirpation were 17.05 (95% CI: 8.39, 40.37) and 17.09 (95% CI: 8.03, 43.48) times higher in the historical fire return interval relative to fire suppression in the low and high dispersal probability scenarios, respectively.
Comparing odds between dispersal probabilities within the decreasing fire return interval scenarios for the recent scenarios, odds of functional extirpation were 1.76 (95% CI: 1.58, 1.95) times higher in low relative to high dispersal scenario. In contrast, the odds of functional extirpation between dispersal probabilities in the historical fire return interval and suppression scenarios were similar; the odds of functional extirpation in low dispersal probability scenarios were 1.14 (95% CI: 0.88, 1.48) and 1.14 (95% CI: 0.36, 3.71) times higher than high dispersal probability in the historical and suppression scenarios, respectively.

DISCUSSION
We created a matrix projection model to estimate the effects of climate-induced increases in fire frequency and dispersal probability on the functional extirpation probability of a successional whitebark pine metapopulation in the GYE. To ensure the model realistically captured whitebark pine population dynamics similar to those in the GYE, we initially parametrized the model with historical fire return intervals of upper subalpine forests in the GYE and mature tree survival rates consistent with healthy pine stands not impacted by WPBR or MPB outbreaks. In this historical parameterization, median-projected densities were comparable to empirical density estimates derived from FIA plots in the GYE, but projected third quartiles were lower (FIA 2019). However, FIA data were collected from a probabilistic sample of U.S. forest ecosystems (Burrill et al. 2018), and sampled stands will necessarily include multiple whitebark pine community types. The bimodal distribution of whitebark pine density derived from FIA data shown in Fig. 5 likely reflects two common whitebark pine community types: successional and self-replacing Hoff 1989, Arno 2001). Whitebark pine densities in successional communities are usually lower than in self-replacing communities, suggesting that the lower mode represents successional communities, and the higher mode, self-replacing communities (e.g., GYCCWPS 2011, data from Barringer et al. 2012: Appendix S2). Because we parameterized the model to reflect a successional whitebark pine stand, we would not expect the model to project the higher densities found in self-replacing stands.
In scenarios parameterized with mature tree mortality rates incorporating WPBR and MPB (i.e., recent conditions), odds of functional extirpation using the historical fire return interval were 13.26 and 17.09 times higher than with fire suppression in low and high dispersal scenarios, consistent with expectations that fire suppression will not lead to complete replacement by shadetolerant competitors (Campbell andAntos 2003, Amberson et al. 2018). However, final median sapling and mature tree densities were not discernably different among historical fire return interval and suppression scenarios, which is inconsistent with expectations that whitebark pine densities will decline over time in seral communities where wildfires have been suppressed (Murray et al. 2000, Amberson et al. 2018. Our findings suggest that mechanisms other than reduced regeneration success may explain declining densities with prolonged fire suppression, but are otherwise consistent with other modeling efforts of whitebark pine response to climate change, which have found higher cover and basal area in suppression scenarios . Declining densities may result from Clark's nutcracker caching preferences; nutcrackers may prefer caching in open and recently burned areas relative to closedcanopy forest (Tomback et al. 2001b). The model was not parameterized to reflect caching preference (i.e., we assumed equal caching probabilities over all densities); it is possible that incorporating caching preference would result in decreased densities in late-successional forests.

Increasing wildfire frequency
Although odds of functional extirpation were higher with more frequent than historical fire, sapling and mature tree densities at year 500 in surviving iterations did not differ between the decreasing fire return interval high dispersal probability scenario and historical fire frequency scenarios (based on overlapping IQRs). In contrast, our simulations projected drastically lower final sapling and mature tree densities in the decreasing fire return interval low dispersal probability scenario relative to the historical fire frequency scenarios. Keane et al. (2017) created a landscape simulation model for the East Fork of the Bitterroot River, Montana, and the Crown of the Continent Ecosystem, Montana, to assess whether restoration actions are likely to be successful in future climates. Two of their scenarios comprised future and current time period simulations of whitebark pine forests without fire management or restoration actions. Basal area increased in the future relative to current time periods in the Bitterroot and declined in the Crown. Keane et al. (2017) attributed basal area declines in the Crown to increased WPBR incidence and mortality that decreased seed production. In contrast, basal  area in the East Fork of the Bitterroot increased despite increased area burned annually. Our findings that, relative to historical scenarios, tree density decreased in the decreasing fire return interval low dispersal probability scenario and remained similar in the decreasing fire return interval high dispersal scenario (as determined by overlapping IQRs) highlight the importance of seed availability and successful dispersal, and also suggest that given adequate seed dispersal, whitebark pine may be minimally affected by increasing fire frequency. However, declining seed production resulting from WPBR mortality and increasing distance to nearby seed sources caused by increasing fire size may work together to increase the functional extirpation risk and highlight the importance of maintaining a network of healthy seed sources.
Empirical findings and simulation models show that distance to seed source governs forest resiliency and the likelihood and rate of recovery through its influence on postfire tree regeneration, which our findings support (Donato et al. 2009, Leirfallom et al. 2015, Kemp et al. 2016, Harvey et al. 2016b, Urza and Sibold 2017, Stevens-Rumann et al. 2018. Further, temporal trends in landscape patterns of stand-replacing fire suggest that a higher proportion of burned areas may become isolated from nearby seed sources because of increased fire size and altered burn shapes (Harvey et al. 2016b). Consequently, many studies conclude that altered fire regimes will lead to declining forested area (Kemp et al. 2016, Harvey et al. ,2016b, c, Keane et al. 2017, Stevens-Rumann et al. 2018). However, the impacts of these changes have been investigated primarily in wind-dispersed species, where 95% of seeds are transported fewer than 180 m (Alexander et al. 1984, McCaughey andSchmidt 1987). In simulations with shortened fire return intervals, Hansen et al. (2018) found nearly complete regeneration failure (<50 stems/ha) for common wind-dispersed species in the GYE, except serotinous lodgepole pine, when distance to the nearest seed source was more than 1 km. In contrast, we found >75% chance of persistence in the high dispersal/ decreasing fire return interval scenario, representing a distance to seed source of~10 km, order of magnitude farther than the distance for which Hansen et al. (2018) found regeneration failure was a near certainty. This suggests that long-distance seed dispersal through avian seed transport and caching may increase confer resiliency following large-scale and severe disturbances. As a result, predictions of forest loss caused by changing fire regimes may be overestimated in communities with avian-dispersed conifer species. Colonization of disturbed areas may readily occur with seeds from distant unburned forest communities through long-distance seed dispersal (e.g., Tomback 2005, Pesendorfer et al. 2016; this factor may shift community types toward those dominated by species with long-distance seed dispersal and/or mitigate future forest loss in the face of unprecedented ecosystem perturbation. It further indicates that multiple species with different life-history traits should be considered before making conclusions about future forest community structure and composition, because species with unique lifehistory traits may be more resilient to climate-altered disturbance regimes. Defining an appropriate threshold below which populations are considered extinct or extirpated is a long-debated topic in the ecological community. Extinction has been defined variably, such as the threshold below which population declines result in changes to community function and stability, the population size leading to loss of other species in the community, or when reproductive failure occurs (Bull et al. 2009, Säterberg et al. 2013, Baumsteiger and Moyle 2017, all of which are further complicated in plant species by lag times between extinction events and loss of the last individual (Vellend et al. 2006, Cronk 2016. Because functional loss can occur before loss of the last individual, we chose to focus on functional extirpation. We selected our functional extirpation threshold (sapling and mature tree density <0.01 trees/ha), because low probability of nutcracker occurrence and high probability of seed predation mean seed dispersal are unlikely. In reality, loss of keystone and foundational functions provided by whitebark pine will occur at higher densities than our threshold because functional roles provided by whitebark pine including food and habitat provisioning for wildlife and facilitation of postfire recovery may cease at higher densities.
Other species-specific studies investigating the effect of climate-altered fire regimes have had to make similar decisions about thresholds for change, whether it be based on regeneration density, tree cover, or another metric . However, most of these decisions have lacked linkage to provisioning of ecosystem services or ecosystem function and more studies linking tree density to functional roles are necessary. Defining these thresholds using functional definitions would improve inference, promote maintenance of ecosystem function and provisioning of ecosystem services, aid management strategy development, and provide a metric for management success.
Because this model was nonspatial in nature, we could not incorporate fuel-fire feedbacks. There remains a great amount of uncertainty relating to how feedbacks may be impacted by climate change. Although fire reduces fuels and therefore limits the severity and spread of future fires and reburns (Harris andTaylor 2017, Parks et al. 2018), these limiting effects may be shortlived or overcome by hot, dry conditions (Harvey et al. 2016a, Parks et al. 2018. Our findings are predicated on fire severity distributions remaining static. If fires become more severe in the future, odds ratios comparing functional extirpation in decreasing relative to historical fire return interval scenarios may be underestimated because of higher mortality rates associated with more severe fires. A similar bias would be present if fire spread increases. In contrast, decreasing fire severity and/or spread would result in reduced mortality, and creation of novel fuel breaks could produce local refugia to which the spread of fire is unlikely (Meddens et al. 2018). We also do not directly consider fire size. We assume that each fire perimeter encompasses the entire whitebark pine population and that the fire severity distribution of a fire is an appropriate fire severity distribution for the whitebark pine population. Although wildfires of mixed fire regimes characteristically produce heterogenous patches of fire severity, the size of highseverity patches tends to increase with fire size (Harvey et al. 2016c). As such, if the fire perimeter becomes much larger than the areal extent of the population, the fire severity distribution and therefore mortality rates may not match what we would expect to find on the ground.
Further, our results are contingent upon demographic rates remaining similar. Postfire regeneration failure is expected to become more common over time (Harvey et al. 2016b, Tepley et al. 2017, Stevens-Rumann et al. 2018), but growth rates of established trees may increase as a result of conditions favorable to growth (Powell and Hansen 2007, Hansen et al. 2016. Decreasing regeneration success could increase extirpation probability, even if those individuals that do establish have higher growth rates.

Management implications
Current fire management actions in whitebark pine ecosystems focus on prescribed fire, wildland fire use, and mechanical thinning to mediate future fire severity and slow successional advance (Keane 2018). However, management actions that incorporate burning may be counterproductive as fire regimes continue to change, because they focus on increasing wildfire in an era when the fire problem in whitebark pine ecosystems is shifting from too little to potentially too much. Flexible management approaches will be required as fire frequency increases and whitebark pine mortality from WPBR accelerates.
Our projections suggest decreasing fire return intervals will increase the odds of functional extirpation of whitebark pine in the GYE by 8.26 times relative to historical fire return intervals, assuming fire regimes change as predicted by Clark et al. (2017) and demographic rates in the GYE remain similar to those we used. Fire suppression, and in particular spot suppression techniques (i.e., suppression tactics that are used on a small area, often to protect valuable assets), may become increasingly important to whitebark pine population persistence in the GYE to mitigate the increased functional extirpation risk caused by increasing fire frequency and protect valuable whitebark stands that contain trees putatively resistant to WPBR that are integral to restoration activities (Keane 2018).
Efforts are being made to identify the spatial distribution of management actions for many species at risk from climate change (e.g., Hennon et al. 2012, Ireland et al. 2018. Ireland et al. (2018) created climate-informed management strategies and identified areas within the GYE that may increase whitebark pine persistence but did not consider population connectivity. Our findings suggest that priority locations could be v www.esajournals.org refined by developing a spatial network of populations connected by seed dispersal, potentially with a focus on locations likely to be persistent fire refugia (e.g., stands with geophysical barriers that prevent fire spread; Meddens et al. 2018). Identifying locations within the GYE to employ targeted suppression efforts would improve resiliency to altered fire regimes by reducing functional extirpation risk, provision ecosystem services by maintaining landscape spatial heterogeneity (Turner et al. 2013), and facilitate adaptation by promoting seed dispersal to areas that may become newly suitable as a result of climate change (Parmesan and Yohe 2003, Millar et al. 2004, Chang et al. 2014.