Increasing vulnerability of an endemic Mediterranean-climate conifer to changing climate and fire regime

Davis, Frank W.; Parkinson, Anne-Marie; Moritz, Max A.; Park, Isaac W.; D’Antonio, Carla M.

doi:10.3389/ffgc.2025.1516623

ORIGINAL RESEARCH article

Front. For. Glob. Change , 30 January 2025

Sec. Fire and Forests

Volume 8 - 2025 | https://doi.org/10.3389/ffgc.2025.1516623

Increasing vulnerability of an endemic Mediterranean-climate conifer to changing climate and fire regime

Frank W. Davis¹^*^†

Anne-Marie Parkinson²

Max A. Moritz^1,3^†

Isaac W. Park⁴^†

Carla M. D’Antonio⁵^†

¹Bren School of Environmental Science and Management, University of California, Santa Barbara, CA, United States
²Santa Barbara County Fire Safe Council, Santa Barbara, CA, United States
³University of California Cooperative Extension, Oakland, CA, United States
⁴Biology Department, Georgia Southern University, Statesboro, GA, United States
⁵Department of Ecology, Evolution and Marine Biology, University of California, Santa Barbara, Santa Barbara, CA, United States

Conifer species of Mediterranean-type ecosystems (MTEs) are especially sensitive to coupled changes in climate and fire regime. Recently California’s MTE conifers have experienced elevated mortality, range shifts, and decreased abundance relative to broadleaf trees. Restricted-range species may be especially vulnerable to environmental change and are a high priority for climate change assessment and adaptation management. We used species distribution modeling and wildfire probability models to assess the potential vulnerability of Pseudotsuga macrocarpa (Bigcone Douglas fir, hereafter BDF), a restricted-range conifer of southern California mountains that is a species of concern because of its ecological importance and recent mortality due to drought and wildfire. We also modeled the distribution of Quercus chrysolepis (CLO), a widespread oak that can be co-dominant with BDF and that can limit fire spread and reduce crown fire risk to BDF compared to the risk from surrounding chaparral vegetation. Ensembles of MAXENT and General Boosted Models were produced for a 2.1 million ha model domain encompassing BDF’s 672,000 ha range for the period 1981–2010, and were projected to 2040–2069 and 2070–2099 using three climate models (warmer-wetter, warmer, and hotter-drier) under a high emission scenario (IPCC 5th Assessment, RCP 8.5). Wildfire projections were analyzed only for mid-century based on warmer-wetter and hotter-drier climate scenarios. Our analysis contributes several important new insights: (1) topoclimatic habitat for BDF could shrink by 19–57% by mid-century, especially for southernmost populations; (2) by mid-century, wildfire probability could increase 2–4 fold in BDF habitat, potentially diminishing the value of mid-to-high elevation sites as topoclimatic refugia; (3) CLO could play a diminishing role in providing refuge from severe wildfires as soon as mid-century but especially by end-of-century; (4) extensive areas of stable mid-century and late-century habitat are confined to national forests in the San Gabriel Mountains and Transverse Ranges. The climate and fire vulnerability of BDF could be reduced by management actions such as mechanical fuel treatments and post-fire restoration, but these are highly constrained by topography, access, and sensitive wildlife species. Our case study illustrates the rapidly increasing vulnerability of endemic conifer species in MTE forests.

Introduction

Mediterranean-type ecosystems (MTEs) are global hotspots of plant diversity and endemism that are especially vulnerable to coupled changes in climate, fire regime and land use (Franklin et al., 2021; Harrison et al., 2024; Peñuelas et al., 2017; Underwood et al., 2009). As water-limited ecosystems, MTEs are highly sensitive to warming which increases evapotranspiration and drought stress (Clark et al., 2016; Peñuelas et al., 2017; Stephenson, 1998). In MTE forests and shrublands, drought stress is often accompanied by increased susceptibility to wildfire, insect pests and pathogens, along with low post-disturbance seedling establishment, resulting in rapid, directional changes in local species and communities (Doblas-Miranda et al., 2015; Franklin et al., 2021b; Peñuelas et al., 2017; Peñuelas and Boada, 2003). Many tree and shrub species in these ecosystems are endemic with restricted ranges, specialized life history requirements, and limited seed dispersal – traits that increase their risk of local and range-wide extinction under rapid climate change (Casazza et al., 2014; Loarie et al., 2008; Rose et al., 2023; Thuiller et al., 2005). Extinction risk is further compounded by land use change and accompanying habitat loss and fragmentation (Franklin et al., 2021b; Rose et al., 2023).

California is an MTE hotspot that supports high endemic plant diversity (Baldwin et al., 2012), including exceptionally high endemic conifer diversity (Rundel, 2019). In recent decades, California’s conifer forests have experienced elevated mortality, disproportionate loss of large trees, species’ elevational shifts, and decreased abundance relative to broadleaf trees (Harrison et al., 2024). These trends have been attributed to the interrelated effects of changing land use and land management, climate warming, extreme drought, severe wildfires, and insect pests and pathogens (e.g., Allen et al., 2010; Cobb, 2022; Harrison et al., 2024; Hill et al., 2023; Kelly and Goulden, 2008; McIntyre et al., 2015; Thorne et al., 2008).

Future climate-related impacts to California’s endemic species are predicted to be greatest on restricted-range species, especially those whose habitats are vulnerable to urbanization and other land use development (Franklin et al., 2021b; Rose et al., 2023). Although some at-risk conifer species remain vulnerable to private land development, many species with high climate exposure risk are mainly or entirely on public lands or private conservation lands not threatened by development, including such noteworthy rare species as giant sequoia (Sequoiadendron giganteum (Lindl.) J. Buchholz) (Nydick et al., 2018), foxtail pine (Pinus balfouriana Grev. and Balf.) (Schoettle et al., 2022), Santa Lucia fir (Abies bracteata (D. Don) A. Poit.) (Ledig et al., 2006), and Brewer spruce (Picea breweriana S. Watson) (Ledig et al., 2012). For these and other endemic conifers in California, conservation management is focused on mitigating effects of climate change, wildfire, and pathogens through actions such as fuels management, prescribed burning, pathogen host management, and climate-adapted reforestation (Maxwell et al., 2022; Millar et al., 2007; North et al., 2019).

Bigcone doulas fir (Pseudotsuga macrocarpa (Vasey) Mayr) (hereafter BDF) is a long-lived conifer endemic to southern California, with a restricted range of less than 672,000 ha (Little and Viereck, 1971)—only 0.6% of the range extent of its widespread, well-studied, parapatric congener, Douglas fir (P. menziesii (Mirb.) Franco). BDF is confined to the Transverse and Peninsular Ranges of southern California, and current patchily distributed populations are almost entirely located on four national forests managed by the U.S. Forest Service (Kauffman et al., 2017; McDonald and Littrell, 1976; Sawyer et al., 2009). Federal land managers consider BDF forests a regionally significant biological resource, in part because they provide important habitat for sensitive wildlife species including the California spotted owl (Strix occidentalis occidentalis) (Hennessy et al., 2024). The species is listed on the IUCN Red List as near threatened (Farjon, 2011).

Several features of BDF make it a compelling case study of climate change exposure for an endemic, range-restricted conifer. Since 1900 southern California has warmed faster than other parts of the state (Cordero et al., 2011). Southern California has also experienced the largest increase in biologically significant drought as measured by climatic water deficit (CWD) (Rapacciuolo et al., 2014) or an ombroxeric drought index (ODI) (González-Pérez et al., 2023). Mean annual precipitation has not changed significantly, but the wet season appears to be shortening, with increasing interannual variability and deeper drought years (Swain, 2021). Recent research has documented BDF’s drought sensitivity and heightened mortality during multi-year droughts (Post-Leon et al., 2022; Wolfgang, 2021), and continued climate warming will likely increase drought exposure of some populations. Moreover, wildfire has increased in southern California Mountains in recent decades (Li and Banerjee, 2021) and both frequency and intensity are projected to increase rapidly by mid- century (Yue et al., 2014). Although unusually fire-adapted for a conifer in the Pinaceae, BDF populations can be locally eliminated in high-intensity wildfires (Howard, 1992; McDonald, 1990; Parkinson et al., 2022), with poor post-fire regeneration (McDonald and Littrell, 1976; Runte et al., 2022).

On the other hand, several mitigating factors could reduce BDF’s vulnerability to climate change. The species mainly occurs in rugged mountainous terrain where there is potential for local topoclimatic buffering and drought refugia (Parkinson et al., 2022; Post-Leon et al., 2022). The risk of high-intensity wildfire in BDF stands varies depending in part on local plant community composition and topography. Populations growing in or adjacent to chaparral suffer higher wildfire mortality than populations growing with canyon live oak (Quercus chrysolepis Liebm.) (CLO), an arborescent evergreen oak that is far less prone to crown fires than chaparral and may limit fire spread and reduce fire intensity in mixed BDF-CLO stands (Parkinson et al., 2022). Very steep slopes with little or no understory vegetation can also serve as topographic fire refugia (Parkinson et al., 2022).

In summary, BDF has some attributes that would appear to make it highly vulnerable to climate change such as small range size, long generation time, and occupancy of sites prone to severe wildfire, and other attributes such as topoclimatic and fire refugia that could serve to reduce its exposure (Franklin et al., 2021b). The ecological importance and restricted range of BDF beg the question of its future distribution under projected climate change. Our objective in this paper is to apply species distribution models (SDM) and wildfire probability models to quantify the bioclimatic niche and wildfire risk of BDF and CLO in BDF’s current range in southern California, and to then apply the SDMs for each species to examine the extent of historic and projected future distributions of BDF and CLO under mid-century and end-of-century climate scenarios. We evaluate changes in wildfire risk only for the mid-century. Our research is focused primarily on BDF and addresses the following questions: (1) What are the most important factors associated with the occurrence of BDF across its range? (2) How will the geographic distribution of topoclimatic habitat be altered by climate change? (3) How is wildfire probability predicted to change in areas projected to be climatically suitable BDF habitat? (4) How will the distribution of CLO habitat change and to what extent will future CLO habitat overlap with BDF habitat? (5) How will climate change impacts on BDF vary among the four national forests and what does that suggest for current management of the species? We hope that our findings will help forest management planners in southern California identify areas that are less exposed to future climate change and wildfire risk in which to focus ongoing fuels management and post-fire restoration efforts (e.g., Hennessy et al., 2024).

Methods

Study area

We modeled the distribution of BDF across the entire range of the species, most of which is confined to the Southern California Mountains and Valley Ecological Section as defined by Goudey and Smith (1994) (Figure 1). We expanded the study area beyond the section boundary in the northern portion of the range to include additional occurrence points on national forest lands in the Santa Ynez Mountains and Sierra Madre Mountains. We excluded two Ecological Subsections that have been extensively developed (Perris Valley and Hills, Fontana Plain and Calimesa Terraces) along with four arid subsections on the edge of the Mojave and Colorado-Sonoran Deserts that are outside of the current range of the species (Figure 1).

Figure 1

Figure 1. Study area map showing the model domain, presence points for BDF, and national forests of southern California. Gray regions are subsections of the Southern California Mountains and Valley Ecological Section that were excluded from the model domain.

The model domain totals 20,856 km² and encompasses most of four national forests that bound the range of BDF. These include the Cleveland National Forest (CNF), Angeles National Forest (ANF), San Bernadino National Forest (SBNF), and Los Padres National Forest (LPNF). Land and Resource Management Plans are maintained for each individual Forest, guided in part by strategic regional planning processes involving all four National Forests.

Regional climate

The study region experiences a Mediterranean climate with hot, dry summers and cool, wet winters. Annual precipitation ranges from 190 mm to 1,390 mm/yr., increasing from north to south but more importantly with elevation, which ranges from 40 m to 3,500 m AMSL (Supplementary Figure S1). Higher elevations in the eastern and southeastern portions of the region can receive significant summer monsoonal precipitation, averaging 20–30 mm/yr. Mean annual climatic water deficit (CWD) – the difference between potential evapotranspiration (PET) and actual evapotranspiration (AET)—varies from <700 mm at high elevations to >1,200 mm at low elevations away from the modifying influence of the coast (Supplementary Figure S1).

Study species

Within its small geographic range, BDF can be locally abundant across a broad elevational range (300–2,700 m, although most populations occur between 1,000 and 1800 m), on many substrates, and within a wide array of vegetation types ranging from chaparral and oak woodlands to mixed conifer forests (Kauffman et al., 2017; McDonald and Littrell, 1976; Sawyer et al., 2009). The species achieves this distribution through a combination of remarkable adaptive traits and the high topoclimatic and edaphic heterogeneity of the southern California mountain ranges. BDF is both long-lived (up to 600+ years), relatively drought tolerant, fire resistant due to its thick bark, and unique among western United States conifers in its ability to resprout profusely from epicormic buds after fires of low to moderate severity (Gause, 1966; Howard, 1992; McDonald, 1990). At elevations below ~1,400 m BDF is patchily distributed on relatively mesic sites on steep north-facing slopes and near the coast where temperature extremes are moderated (Borchert and Hibberd, 1984; Post-Leon et al., 2022; Wolfgang, 2021). Between 1,500–2000 m, density of BDF increases with elevation and the species forms more extensive, continuous groves (Leiberg, 1899). BDF abundance declines above 2000 m.

Site water balance and fire regime both contribute to this distribution pattern at low-to-mid elevations; competition and low temperatures may explain declining abundance at higher elevations. Although a relatively drought tolerant conifer, the species is constrained by site water availability and is vulnerable to drought-induced mortality. Dendroecological studies have documented the sensitivity of the species’ growth to interannual precipitation (Haston et al., 1988; Wolfgang, 2021). Minnich et al. (2016) documented 25% mortality of BDF in the San Bernardino and San Jacinto Mountains during an extreme drought year in 2001–2002. Ecophysiological studies during the historic 2012–2016 drought showed that some northern populations were operating at the edge of critical drought thresholds (turgor loss point, stomatal closure, loss of hydraulic conductivity), particularly at more interior sites on south-facing slopes (Post-Leon et al., 2022). A 2015–2016 survey of BDF stands in the Angeles National Forest documented drought-related tree mortality exceeding 10% (mainly 10—40%) in 19% of sampled stands, mostly at lower elevation sites (Kauffman et al., 2017).

Similarly, despite being fire-adapted, BDF individuals are susceptible to mortality from high severity fire (Howard, 1992; McDonald, 1990; Parkinson et al., 2022) and can experience poor regeneration post-fire (McDonald and Littrell, 1976). Many low-elevation stands are in areas dominated by chaparral where wildfires occur as high-severity crown fires. In those areas BDF stands persist on steep slopes with reduced understory fuel loads or in mesic sites also supporting CLO. CLO thickets can act as a buffer against the spread of chaparral fires due to the higher canopy moisture content and limited understory fuel accumulation in dense stands (McDonald and Littrell, 1976; Minnich, 1977; Parkinson et al., 2022).

CLO is native to Oregon, California, Nevada, Arizona, and Mexico, and is the most widely distributed oak in California as a component of many plant communities (Tollefson, 2008). In the mountains of southern California CLO is often codominant with BDF on mesic, steep, north-facing slopes (McDonald and Littrell, 1976; Minnich, 1999) at elevations between 1,200 m and 2000 m AMSL (Borchert and Hibberd, 1984; Tollefson, 2008). Tree mortality from wildfires decreases with tree size, and top-killed trees resprout vigorously from the root crown or bole (Regelbrugge and Conard, 1993; Tollefson, 2008).

Occurrence data

Occurrence data for BDF were compiled from multiple sources: 2,189 occurrence points were compiled from the Global Biodiversity Information Facility (GBIF, www.gbif.org), Calflora¹, and the CalJep database (Viers et al., 2006). Another 177 locations were obtained from the publicly available USFS Forest Inventory and Analysis (FIA) plot database (Gray et al., 2012). To ensure representation of BDF on very steep and/or remote sites, we compiled another 986 points from vegetation map sources including CALVEG (U.S. Forest Service, 2018; Kauffman et al., 2017; Parkinson et al., 2022). Map polygons labeled as BDF in each dataset were combined and points were randomly selected inside the polygons. These points were then visually inspected using the most recent Google Earth imagery to verify that a BDF tree – recognizable by its distinctive crown—was within 135 m of the point. This distance was selected to match the 270 m resolution of the climate grids used in the SDM analysis. To avoid multiple BDF observations in the same grid cell, the final pooled dataset of 3,352 points was thinned to 1,954 points using the R package spThin and a minimum separation distance of 270 m (Aiello-Lammens et al., 2015) (Figure 1).

Occurrence data for CLO were compiled from FIA plot data, the CalJep database and GBIF. Only 370 points were available for CLO but these points were well distributed across the model domain (Figure 2).

Figure 2

Figure 2. BDF (left) and CLO (right) habitat suitability scores for the period 1981–2010 based on the ensemble SDMs combining GBM and MAXENT models. Habitat scores can range from 0–1,000. BDF presence points are colored to distinguish true positives (yellow) from false negatives (magenta). CLO presence points are colored to distinguish true positives (yellow) from false negatives (blue).

Predictor variables

We used the 2014 California Basin Characterization Model (BCM) grids to represent historic and future topoclimatic conditions http://climate.calcommons.org/bcm. BCM topoclimatic and ecohydrological data gridded at 270 m resolution have been used extensively to model climate change impacts on California topoclimate, hydrology, and species distributions (e.g., Franklin et al., 2013; Rose et al., 2023; Underwood et al., 2018).

We explored multiple species distribution modeling approaches and bioclimatic variables for the period 1981–2010 that were expected to have a strong influence on BDF habitat suitability. Based on the correlation among BCM and topographic variables, their importance in preliminary SDMs, and visual inspection of habitat suitability maps, we narrowed the set to seven variables: mean annual Climatic Water Deficit (CWD_ANN), mean annual precipitation (PPT_ANN), mean summer (June through August) precipitation (PPT_JJA), mean daily minimum temperature for December through February (TMIN_DJF), mean maximum daily temperature for June through August (TMAX_JJA), slope angle (SLOPE_MAX), and total March–May solar radiation (SOLAR_MIN). We used Spring solar radiation as opposed to total annual radiation because tree carbon assimilation is high during this time of year and spatial variation in radiation related to slope orientation and horizon shading is relatively large (Davis and Dozier, 1990; Turner et al., 2020). The variables CWD_ANN, TMIN_DJF and TMAX_JJA were most highly correlated (r² = 0.58). Squared correlations among other variables ranged from 0.01–0.42 (Supplementary Table S1).

Slope angle and seasonal solar radiation were derived from a 90 m digital elevation model from the Shuttle Radar Topography Mission (SRTM3, http://www.webgis.com/srtm3.html). Total March–May solar radiation was modeled using the Solar Radiation tool in ArcGIS Pro v. 2.4. Slope and radiation grids were resampled to 270 m to match the resolution of BCM data. To better capture local steep slope and low radiation environments that have been associated with BDF, we rescaled the 90 m grids to 270 m using the maximum slope and minimum solar radiation for the nine 90 m cells in the 270 m BCM cell.

Climate projections

We used projections from three Global Climate Models (GCMs) from the Fifth IPCC Climate Assessment (CMIP5) under the high emissions scenario (RCP 8.5): CCSM4 (National Center for Atmospheric Research, United States), MIROC-ESM 8.5 (Center for Climate System Research, Japan), and CRNM-CM5 8.5 (Centre National des Recherches Météorologiques, France). CCSM4 is close to the ensemble mean for 18 GCMs used in the California Climate Assessment, CNRM-CM5 is a warmer-wetter projection, and MIROC-ESM is a hotter-drier projection. These models have been identified as well suited to California climate change impact assessments (Flint and Flint, 2014; Underwood et al., 2018). The CMIP5 GCM models were produced at a grid resolution of ~6 km². Grids were statistically downscaled to 270 m as described by Flint and Flint (2012). SDM projections were produced using 30-year means for mid-century (2040–2069) and end-of-century (2070–2099).

SDM modeling

Presence-only data were used to model the topoclimatic habitats of BDF and CLO, with background sample sizes of 15,000 and 10,000, respectively, or roughly 5.2 and 3.4% of the model region (Mapped BDF stands occupy roughly 1% of the model region.) Large background sample sizes to have been shown to improve the performance of presence-only SDMs (Valavi et al., 2022). We explored the effectiveness of four different modeling approaches including Maximum Entropy (MAXENT), General Boosted Models (GBM), Random Forest (RF) and General Additive Models (GAM), as well as model ensembles. Based on two model performance measures—Area under the Receiver Operator Characteristics curve (AUC) and the KAPPA statistic—along with visual examination of spatial patterns in model prediction errors, we ultimately selected a performance-weighted mean Ensemble model combining MAXENT and GBM. These methods have been shown to be especially effective for predictive modeling using presence-only data (Valavi et al., 2022).

Models were fitted using the R package biomod2 v. 4.2–4 (Thuiller et al., 2021, 2009). Each model was fitted using 4 random samples of background points and cross-validated using 4 runs each with 80% of the presence points. Variable importance was estimated based on 3 model permutations. MAXENT models were limited to a maximum of 200 iterations. GBM models were based on 500 trees and a minimum node size of 20 observations. Habitat suitability scores were converted to binary maps (non-habitat, habitat) based on maximizing the sum of model sensitivity and specificity, which is commonly used to threshold SDM suitability scores (Fielding and Bell, 1997; Valavi et al., 2022). We tabulated the area of suitable habitat for BDF and CLO for the modeling domain as a whole as well as for the individual national forests to highlight differences in trends for the different forests and the different portions of the range that they encompass.

Wildfire probability

One km² grids of annual wildfire probability for the periods 1981–2010 and 2040–2069 were obtained from the State of California CalFire Forest and Resource Assessment Program (FRAP).² They estimated the probability of a fire in each grid cell in each year using the approach developed by Park et al. (2021), who fitted a Generalized Additive Model to FRAP wildfire event data for the period 1970–2016 based on BCM grids of modeled AET and CWD, along with geospatial grids of housing density, years since fire, distance to roads, and distance to electrical infrastructure. Grid cells with ≥50% cultivated land, urban cover, water, or barren/rock were excluded from the model. FRAP analysts used these variables to fit a GAM model to a spatially randomized subset of wildfires for the period 1980–2018, and then applied that model to map annual wildfire probability for each year from 2019 to 2099 using historical climate and land use data, climate projections from four GCMs (RCP 8.5) and projected changes in housing density and electrical infrastructure.

We calculated the 30-year mean probability for each cell in each period in our BDF model domain. Wildfire risk models were available for 2 of the GCM models used to produce the SDMs for BDF: CNRM-CM5 and MIROC5, under the 8.5 RCP emissions scenario. We oversampled the 1-km wildfire probability grids to 270 m and extracted those grid cells that were predicted as suitable BDF habitat in 1981–2010 and in 2040–2069 based on the SDM. We elected not to analyze end-of-century wildfire probability projections out of concern that dynamic vegetation change over the longer time period, which was not accounted for in the GAM model projections, could appreciably change the statistical relationships between climate, vegetation, and wildfire probability over the longer term (Keeley and Syphard, 2016).

Results

BDF and CLO habitat models

Maxent, GBM and Ensemble models have good skill (AUC = 0.86–0.88, Table 1). The ensemble model effectively predicts true positives (sensitivity = 0.91) but had only modest skill in the rate of false positives (specificity ~0.68). Most false negatives occur at the lower, drier edges of population occurrences and in the southernmost populations (Figure 2).

Table 1

Table 1. Species distribution model results for BDF for the period 1981–2010, including mean model area under the receiver operating curve (AUC) and variable importance for cross-validated general boosted models (GBM), maximum entropy models (MAXENT), and ensemble models.

Annual CWD (CWD_ANN) and slope angle (SLOPE_MAX) are the most important variables in GBM, MAXENT and Ensemble models (Table 1). Summer precipitation (PPT_JJA) along with summer maximum (TMAX_JJA) and winter minimum (TMIN_DJF) temperatures are moderately important in the MAXENT and Ensemble models. Variable response curves for the ensemble model are provided in Figure 3, and those for the separate GBM and MAXENT models are provided as Supplementary Figures S2, S3. Ensemble habitat scores decline steeply as CWD_ANN values rise above 800 mm and increase with local slope at angles between 20 and 40 degrees. Secondarily, BDF habitat suitability is maximized where annual precipitation is between 500–1,000 mm, mean daily minimum winter temperatures exceed −4°C, and maximum summer temperatures exceed 25°C (Figure 3).

Figure 3

Figure 3. Predictor variable response curves for the BDF species distribution model produced as a performance-weighted mean of GBM and MAXENT models. Predictor variables include (A) CWD_ANN, (B) SLOPE_MAX. (C) PPT_ANN, (D) PPT_JJA, (E) TMAX_JJA, (F) TMIN_DJF. The Y-axis is habitat suitability scores, which can range from 0 to 1. The shaded areas in each plot are 95% confidence intervals.

BDF habitat is predicted to occur over 5,908 km² or 28.5% of the model domain (Supplementary Table S2). Eighty-six percent of predicted BDF habitat occurs on national forest land, and is especially extensive on the Angeles National Forest (53% of administrative area) and Los Padres National Forest (49% of area) (Supplementary Table S2). In these Forests, high suitability scores occur in large contiguous areas at ~900–1700 m elevations on steep north-facing slopes (Figure 2).

Model skill was good but slightly lower for CLO than BDF, irrespective of the method, with AUC values of 0.84 for GBM, MAXENT and ENSEMBLE models (Table 2). Similar to the BDF SDM, annual CWD was the most influential variable in all models. Unlike BDF, the next most important variables were summer maximum daily temperatures, summer precipitation, and spring insolation, followed by winter minimum daily temperatures and slope angle. Visual inspection of habitat suitability scores indicates good agreement between CLO occurrence points and model predictions (Figure 2).

Table 2

Table 2. Species distribution model results for CLO: mean model AUC and variable importance for cross-validated general boosted models (GBM), maximum entropy models (MAXENT), and ensemble models.

Roughly 75% of the area predicted to be CLO habitat is also predicted BDF habitat, consistent with the previously reported strong association of these species (Supplementary Table S2). The overlap in species’ habitats is lower (59–67%) in the lands off the national forests and in the more southern parts of the range such as in the Cleveland National Forest (Supplementary Table S2).

Wildfire probability

For the period 1981–2010, modeled mean annual wildfire probabilities in sites modeled as BDF habitat range from <0.005 (i.e., expected fire return interval (FRI) of 200 years) to as high as 0.04 (FRI = 25 years) (Figure 4). Areas at highest risk occur primarily in the San Gabriel Mountains, as well as the San Rafael Mountains in the Los Padres National Forest and Santa Ana Mountains in the Cleveland National Forest. Areas of lowest wildfire risk occur at lower, drier sites, particularly in the upper Cuyama River Basin of the Los Padres National Forests near BDF’s northern range limit.

Figure 4

Figure 4. Predicted annual wildfire probability for the period 1981–2010 in areas modeled as BDF habitat. For visualization, continuous probabilities are shown in 4 classes. The heavy black line indicates the boundary of the SDM modeling domain.

Mid-century climate projections, wildfire probabilities and SDMs

Projected changes in average annual PPT, annual CWD, summer maximum daily temperatures, and winter minimum daily temperatures in the model domain are tabulated in Supplementary Table S3. By mid-century, annual CWD is projected to increase by 71–123 mm compared to 1981–2010 levels as a result of projected warming and changing precipitation. The dry MIROC ESM 8.5 model projects the largest increase, but even the wetter CCS4 model projects a 74 mm increase in annual CWD due to warmer temperatures. Winter minimum and summer maximum daily temperatures are projected to increase by roughly 1.5–2.3^o C compared to the baseline period 1981–2010.

Based on the hot-dry climate projection (MIROC ESM 8.5), BDF habitat extent is predicted to shrink by 57% from 5,909 km² to 2,584 km² across its entire range due largely to reduced annual precipitation and increased annual CWD (Figure 5; Supplementary Table S3). Modeled habitat contracts dramatically in the southern portion of the range and at lower elevations (Figure 6). For example, habitat in the CNF declines from 29 to 1.5% of the administrative area. Habitat extent in the SBNF drops from 28 to 13% of the area. Substantial habitat loss is also projected for the ANF and LPNF with this hotter and drier scenario. Mean suitability of remaining habitat is also reduced in all of the national forests (Supplementary Table S2).

Figure 5

Figure 5. Predicted extent of BDF habitat for historic, mid-21st century and end-of-century 30-year periods for the model domain and for individual national forests based on (A) a hot-dry projection (MIROC ESM 8.5, RCP 8.5), (B) a warm projection with little rainfall change (CCSM4, RCP 8.5), and (C) a warm-wet projection (CNRM CM5, RCP 8.5).

Figure 6

Figure 6. Mean annual probability of wildfire in predicted BDF habitat for the periods 1981–2010 (HIST) and 2040–2069 based on CNRM 8.5 and MIROC ESM 8.5 climate model projections.

The CCSM4 model mid-century projection (warmer but little change in precipitation) also produces substantial reduction in BDF habitat, from 5,909 km² to 4,608 km² across the model domain (Figures 6B, 7B), although this is not as evenly distributed across the national forests. The Cleveland and San Bernardino National Forests lose the most habitat on a percent basis, whereas increased habitat is predicted in the Los Padres National Forest. The warmer/wetter CNRM-CM5 model results in a slightly lower net loss in habitat from 5,909 km² to 4,831 km², with a slight net habitat gain in the Los Padres National Forest (Figures 6C, 7C).

Figure 7

Figure 7. Distribution of mid-century (left column) and end-of-century (right column) BDF habitat based on (top row) a hot-dry projection (MIROC ESM, RCP 8.5), (middle row) a warm projection with lower rainfall change (CCSM4, RCP 8.5), and (bottom row) a warm-wet projection (CNRM CM5, RCP 8.5), stable areas (predicted historic and future habitat) are shown in blue. Habitats predicted to be lost are shown in orange, and new habitats are shown in green.

Based on the MIROC-ESM 8.5 model, habitat for CLO is predicted to diminish from 7,657 km² to 3,824 km² by mid-century (Supplementary Table S2). The area of joint habitat for BDF and CLO shrinks from 4,433 km² to 1731 km² (from 75 to 67% of BDF habitat) in this hot-dry scenario (Supplementary Table S2). The CCSM4 8.5 and CNRM 8.5 model projections also project net habitat loss, although less so than in the hot-dry scenario. Habitat shrinks to 4,665 km² and 4,679 km², respectively across the model domain. However, the overlap with BDF habitat drops from 75% to 60–65% of overall BDF habitat.

By mid-century the annual wildfire probability in modeled BDF habitat jumps 1–3 fold based on the MIROC5 model and 3–4 fold based on the CNRM model (Figure 6). The results are very similar for all of the national forests (Figure 6) and for joint BDF/CLO habitat (not shown).

End of century climate projections and SDMs

Climate models diverge significantly in end-of century projections for annual precipitation and CWD (Supplementary Table S3). MIROC ESM 8.5 and CCSM4 8.5 models project a slight decrease in annual PPT and increase in CWD compared to the baseline period, whereas the CNRM 8.5 model projects a > 100 mm increase in annual PPT across the model domain. The MIROC ESM 8.5 model projects increases in mean winter minimum daily and summer maximum daily temperatures greater than 4°C, compared to ~2 ^o C for the CCSM4 model and ~ 1 ^o C for the CNRM 8.5 mode.

The ensemble SDM projections for BDF range from continuing loss of habitat during the latter half of the 21st century based on the MIROC ESM 8.5 and CCSM4 8.5 models to little change in total habitat for the CNRM 8.5 model (Figure 5). The MIROC ESM 8.5 projection predicts a reduction in topoclimatic habitat from 5,908 km² in 1981–2010 to 2,290 km², including the near-total loss of habitat in the eastern and southern portions of the range and extensive losses in the northwestern domain and Los Padres National Forest in particular (Figure 7; Supplementary Table S2). More habitat persists in the southern portion of the range under the CCSM4 8.5 projection, partly due to modeled increases in annual PPT and decreases in annual CWD (Figure 7). Projected decreases in annual CWD and increases in annual PPT throughout the southeastern portion of the range in the CNRM 8.5 model are predicted to result in a significant increase in BDF habitat relative to both historic and mid-century model projections (Figure 7). Even under this relatively wet projection, significant portions of historic range are lost in the northern range, notably in the Los Padres, Angeles and San Bernardino National Forests.

Late-century model projections for CLO indicate habit loss from 7,657 km² in 1981–2010 to only 1,413 km² based on the MIROC ESM 8.5 and to 5,137 km² based on the CCSM4 8.5 model. In contrast, CLO habitat is projected to expand relative to historic baseline to 8,468 km² based on the CNRM 8.5 model (Supplementary Table S2). The area that is jointly suitable for BDF and CLO decreases under all 3 model projections, from 4,433 km² to 380–3,122 km² depending on the model (Supplementary Table S2).

Discussion

Numerous studies have used SDMs to investigate climate change impacts on tree species in MTE settings (e.g., Casazza et al., 2014; López-Tirado et al., 2024; Peñuelas and Sardans, 2021; Serra-Diaz et al., 2014). Suitable bioclimatic habitats for northern hemisphere MTE conifers and hardwoods are generally predicted to shift northward and upslope, but the magnitude and the direction of range shifts vary depending on the species, the scale of the analysis, and on whether projections account for disturbances such as land use change, fire, and pests and pathogens (e.g., Franklin et al., 2013; Ruiz-Labourdette et al., 2012; Serra-Diaz et al., 2014). Both empirical data and dynamic ecosystem models tend to show an increase of hardwoods over conifers due to the higher sensitivity of conifers to stand-replacing wildfires, drought, pests and pathogens (Serra-Diaz et al., 2018; Takolander et al., 2019).

Despite its significance as an ecologically important and range-restricted endemic conifer, BDF’s climate vulnerability has received limited attention and only as part of larger multi-species analyses (Franklin et al., 2013; Loarie et al., 2009; Potter et al., 2017; Rose et al., 2023). Based on a moderate mid-century climate scenario (RCP 4.5) and a multi-criterion assessment framework, Potter et al. (2017) classified BDF as having high vulnerability, with moderate adaptation potential based on regeneration capacity, genetic variability and ecological requirements. Rose et al. (2023) estimated that the species is not threatened by land development given its distribution on public lands but has moderately high exposure to climate change by end-of-century under the RCP 8.5 emissions scenario.

We obtained a relatively detailed, spatially explicit view of the species’ climate vulnerability within its current range based on a very large training dataset and using climate and other environmental and biotic factors deliberately chosen based on the known ecology of the species. Our results contribute several important new insights, notably: (1) the projected pace of climate change in southern California is fast enough that the species’ topoclimatic habitat will shrink significantly (19–57%) by mid-century, especially for southern populations in the Cleveland National Forest; (2) by mid-century, wildfire probability will increase significantly in BDF habitat, potentially diminishing the value of mid-to-high elevation sites as topoclimatic refugia; (3) CLO will play a diminishing role in providing refuge from severe, stand-eliminating wildfires as soon as mid-century but especially by end-of-century because of reduced overlap between its future suitable habitat and that of BDF; (4) extensive areas of stable mid-century and late-century habitat are confined to two national forests in the San Gabriel Mountains (ANF) and Transverse Ranges (LPNF).

Topoclimatic controls on BDF distribution

Many studies have found climatic water deficit (CWD) to be a particularly important topoclimatic variable associated with the distribution of California vegetation types (e.g., Ackerly et al., 2015; Stephenson, 1998), plant species establishment and distribution (e.g., Ackerly et al., 2020; Davis et al., 2016), and tree mortality (Das et al., 2013; van Mantgem et al., 2013). The importance of CWD and annual precipitation in our ensemble SDM for BDF is consistent with dendroecological studies showing the sensitivity of the species’ growth to interannual drought (Haston et al., 1988; Wolfgang, 2021). The lesser importance of temperature variables is also consistent with the finding by Haston et al. (1988) of only a weak negative effect of temperature on BDF growth during spring and early summer.

The Basin Characterization Model (BCM) estimates of CWD account for regional temperature and precipitation gradients, local topoclimatic variation in solar radiation and associated evaporative demand, and variation in soil water balance as influenced by the timing of precipitation and snowmelt along with mapped soil properties such as water holding capacity and wilting point (Flint et al., 2013). As such, CWD as modeled in the BCM is a highly integrative topoclimatic variable that is often influential in SDMs for plant species in California and other areas of the western U.S. (e.g., Ackerly et al., 2020; Anderegg et al., 2015; Rose et al., 2023).

Increased wildfire risk and divergence of BDF and CLO habitats

The risk of severe wildfire in existing BDF stands has increased significantly in recent decades and that trend is likely to continue through at least mid-century (Nigro and Molinari, 2019). Thus wildfire refugia will likely continue to play an important role in the distribution and abundance of BDF. A study evaluating habitat features of stands where BDF either survived or recovered from two >81,000 ha wildfires in the LPNF—the” Zaca Fire” in 2007 and the “Thomas Fire” in 2017—concluded that refugial locations include two general types of habitat (1) more arid portions of the species range where the understory is less vegetated or is largely rocky, and (2) sites supporting stands of Q. chrysolepis (Parkinson et al., 2022). Regarding the former, these more arid portions of the LPNF as with CNF and SBNF, will not be suitable for BDF in the future regardless of fire because of increased CWD. Hence these ‘arid’ type fire-refugia may no longer exist for BDF. By contrast, habitat associations with CLO, the second type of fire refugia, will likely continue to be important as both climate and fire refugia into the future. However, our results suggest a diminishing overlap of CLO and BDF habitat by mid-century.

Because BDF seeds are not stimulated by fire and have limited dispersal potential (Minnich, 1977; Vander Wall et al., 2006), the long-term persistence of BDF stands is largely dependent on seed sources from trees within or very near to the burned stand. Unlike other conifers, BDF also does not appear to expand from its existing locations into neighboring chaparral (Russell, 2018). Considering all these factors, the current distribution and quantity of BDF is unlikely to increase, but will only at best remain constant or, more realistically, continuously decline without intervention because of increasing CWD and wildfire. We hasten to add that the wildfire projections used here do not consider fire severity, so it is possible that future fires of lower severity could have lower impact or even contribute to the persistence of BDF stands.

Modeling limitations and caveats

Climate models for southern California show high concordance in projecting 21st century warming, although the magnitude of warming diverges among models over the latter half of the century and depends strongly on the emission scenario (Underwood et al., 2018). Model projections of precipitation means and extremes are highly uncertain for southern California compared to the rest of California or the contiguous U.S., and this uncertainty increases from mid-century to end-of-century. This is true for both the CMIP5 model projections used here as well as the more recent CMIP6 model projections (Petrova et al., 2024; Polade et al., 2017).

Species distribution modeling entails many assumptions and decisions that significantly influence model predictions, as has been extensively discussed in the SDM literature (e.g., Franklin, 2009; Pearson and Dawson, 2003). We will not repeat those concerns here except to address two specific issues: species interactions and non-stationary disturbance regimes.

SDMs often assess climate change impacts on biodiversity by summing independent species models (e.g., Loarie et al., 2008). However, the assumption of species’ independence is problematic when the species influence their respective distributions through competitive or mutualistic interactions (Araújo and Luoto, 2007). We could have included the SDM suitability scores of CLO as a predictor variable in the BDF SDM (e.g., Araújo and Luoto, 2007), or by fitting joint distribution models (Clark et al., 2014), but opted to treat the species independently for several reasons. First, the interaction of BDF and CLO is complicated by the fact that individuals may compete for space and light at the scale of neighbors but CLO may also buffer BDF from severe fires at site-to-landscape scales. Additionally closeness to CLO individuals can increase BDF seedling survival (Runte et al., 2022). In other words, interactions could be positive or negative as a function of scale. Second, it is clear that BDF does not depend on the presence of CLO given its occurrence in a variety of plant communities, arguing for a more individualistic treatment of the two species (Gleason, 1926). Finally, while we might have obtained better predictions for the two species through joint distribution modeling (Clark et al., 2014), model input data would have been limited to plot data where the presence or absence of both species was recorded. By modeling the species separately, we were able to take advantage of a much larger and presumably less bias-prone occurrence data set specifically developed for BDF.

The static and correlative nature of species distribution models limits their power to predict changes in species’ distributions under rapid climate change where dynamic processes like disturbance and dispersal can play critically important roles (Franklin, 2010). Such processes are better explored through dynamic landscape models. Franklin et al. (2005) included BDF in a dynamic landscape simulation of the response of dominant plant species to alternative fire regimes in southern California mountain and foothill environments based on historic fire information and without considering climate trends. BDF – which was parameterized as long lived, shade tolerant, and resprouting—increased under all three fire regimes with fire return intervals ranging from 30–687 years. Additional landscape simulations with current landscape models that incorporate more detailed fire behavior and fire effects information as well as climate change would be worth pursuing in the light of our findings suggesting high vulnerability of the species to climate change as well as increasing wildfire risk and decreasing habitat overlap with CLO.

Restoration and management implications

Because BDF is distributed almost entirely on federal land free from development pressure, is not a commercial timber species, and generally occupies rugged areas that are not readily accessible for recreational use, its vulnerability depends primarily on its climate change exposure and risk of severe wildfire. Existing national forest management plans explicitly consider mechanical fuel treatment in BDF stands as well as reducing fuel in chaparral adjacent to BDF stands to reduce the likelihood of stand replacing wildfires (e.g., USDA Forest Service, 2005), along with reforestation of severely burned stands by outplanting BDF seedlings (Hennessy et al., 2024; Runte et al., 2022).

Forest managers face several challenges in pursuing these activities to increase the climate resilience of BDF forests. Over much of the species’ range, slopes are too steep (> 30%) to allow mechanical fuel treatments or re-planting. Many sites are remote with limited or no vehicular access, and in those areas fuels management will be cost-prohibitive. Although outplanting can be done, watering is not practical in these areas, resulting in high seedling mortality during the prolonged summer drought (Runte et al., 2022). Even areas that are otherwise accessible may be ineligible for active management due to the presence of sensitive species such as the California spotted owl or to areas being designated as ‘wilderness’, which limits some management actions. These concerns place a premium on spatial planning to prioritize fuels management and restoration efforts in areas where there is also strong local management capacity, community support, and a reasonable likelihood of BDF survival under future climate (Hennessy et al., 2024). Our analysis underscores the need to include near-term, climate-driven changes in habitat suitability and wildfire risk in these planning exercises, as has been done in developing the climate-adapted conservation strategy proposed by Hennessey et al. (2024), where forest managers have placed greater priority on protecting BDF populations at higher elevations that are less vulnerable to climate change but at high risk of severe wildfire. Our results provide additional support for that strategy.

We have also shown that each of the four national forests faces different levels of climate exposure and different opportunities to conserve BDF in topoclimatic and fire refugia. A range-wide climate adaptation strategy will need to consider these differences, especially the increasing importance of the two northern forests (ANF and LPNF) in providing extensive contiguous habitat in a hotter, drier world. Prioritizing areas for potentially costly interventions could require close coordination among these administrative units, particularly given that extensive areas of the LPNF are designated wilderness where management options are more constrained.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

FD: Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review & editing, Data curation, Formal analysis, Visualization. A-MP: Conceptualization, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing. MM: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Validation, Writing – original draft, Writing – review & editing. IP: Conceptualization, Methodology, Supervision, Writing – original draft, Writing – review & editing. CD’A: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. Financial support was provided by the National Fish and Wildlife Foundation (NFWF, ID: 64532) grants to CD’A and MM. FD was supported by the La Kretz Research Center at the University of California Sedgwick Reserve.

Acknowledgments

We thank two reviewers for their helpful reviews. We also thank Nicole Molinari for project advice.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Gen AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2025.1516623/full#supplementary-material

Footnotes

1. ^www.calflora.org

2. ^https://www.fire.ca.gov/what-we-do/fire-resource-assessment-program/fire-probability-for-carbon-accounting

References

Ackerly, D. D., Cornwell, W. K., Weiss, S. B., Flint, L. E., and Flint, A. L. (2015). A geographic mosaic of climate change impacts on terrestrial vegetation: which areas are most at risk? PLoS One 10:e0130629. doi: 10.1371/journal.pone.0130629

PubMed Abstract | Crossref Full Text | Google Scholar

Ackerly, D. D., Kling, M. M., Clark, M. L., Papper, P., Oldfather, M. F., Flint, A. L., et al. (2020). Topoclimates, refugia, and biotic responses to climate change. Front. Ecol. Environ. 18, 288–297. doi: 10.1002/fee.2204

Increasing vulnerability of an endemic Mediterranean-climate conifer to changing climate and fire regime

Introduction

Methods

Study area

Regional climate

Study species

Occurrence data

Predictor variables

Climate projections

SDM modeling

Wildfire probability

Results

BDF and CLO habitat models

Wildfire probability

Mid-century climate projections, wildfire probabilities and SDMs

End of century climate projections and SDMs

Discussion

Topoclimatic controls on BDF distribution

Increased wildfire risk and divergence of BDF and CLO habitats

Modeling limitations and caveats

Restoration and management implications

Data availability statement

Author contributions

Funding

Acknowledgments

Conflict of interest

Generative AI statement

Publisher’s note

Supplementary material

Footnotes

References

94% of researchers rate our articles as excellent or good