Increased whitebark pine (Pinus albicaulis) growth and defense under a warmer and regionally drier climate

1. Introduction

Widespread tree mortality associated with changing climate and disturbance regimes predisposes many western United States forests to future, unpredictable changes in forest stand structure, function, and species composition (van Mantgem et al., 2009; Allen et al., 2010, 2015; Smith et al., 2015; Goeking and Windmuller-Campione, 2021). Long-term drying trends coupled with increased fire activity, pathogen activity, and insect outbreaks, all highlight a need to understand physiological responses of trees to biotic and abiotic stressors (Stephens et al., 2013; Jenkins et al., 2014; Bentz et al., 2016; Dudney et al., 2021). Currently, the majority of research linking drought- and disturbance-related mortality centers on lower-elevation, moisture-limited forest systems (van Mantgem and Stephenson, 2007; McDowell et al., 2008; Das et al., 2013; Zhang et al., 2017; Stephenson et al., 2019) with little attention to higher-elevation, energy-limited systems (i.e., systems limited by growing season mean temperature and length). There is a critical need to evaluate how climate trends and severe drought that are causing widespread tree mortality in low-elevation, moisture-limited systems will affect cooler, short-growing season subalpine forests in the western U.S. More specifically, there is a need to understand how tree growth and defense strategies in energy-limited systems are consistent with, or differ, from moisture-limited systems. This is particularly important for informing ecosystem models that are often parameterized based on data from lower-elevation forest taxa.

Whitebark pine (Pinus albicaulis) is an ecologically and culturally important high-elevation conifer with a split distribution. The western populations range from the southern Sierra Nevada in California (36° N) north through the Coast and Cascade Ranges of Oregon and Washington into west-central British Columbia (55° N). The eastern distribution tracks the northern Rocky Mountains, extending from Colorado poleward through Idaho and Montana into western Alberta (Little, 1971; Arno and Hoff, 1989; Fryer, 2002). Whitebark pine has been recognized as a foundation species in high-elevation montane ecosystems due to its importance in promoting snowpack accumulation and retention, facilitating vegetative establishment, and providing a critical food source for a variety of wildlife (Tomback et al., 2001). However, there is great concern about the future viability of whitebark pine forests due to extensive tree mortality associated with a highly aggressive exotic pathogen, white pine blister rust (Cronartium ribicola), and extensive outbreaks of the native mountain pine beetle (Dendroctonus ponderosae) linked to a warming climate. These factors, along with future projections of continued population declines due to sustained warming (Chang et al., 2014), led the U.S. Fish and Wildlife Service to formally list whitebark pine as a “threatened” species under the Endangered Species Act in December 2022. Managers tasked with supporting and maintaining the viability of whitebark pine into the future need information on how whitebark pine forests respond to climate and disturbance synergies across biophysical gradients throughout its geographic range.

In whitebark pine, and other conifers, defense characteristics play a crucial role in moderating bark-beetle interactions, with trees persisting through bark beetle outbreaks often exhibiting distinct resin duct characteristics (Kane and Kolb, 2010; Ferrenberg et al., 2014; Hood et al., 2015; Hood and Sala, 2015; Kichas et al., 2020). However, research demonstrating linkages to warming temperatures and increased bark beetle activity in the western U.S. (Bentz et al., 2005, 2010, 2016) underscores a need to understand how tree defenses may respond to changing biotic and abiotic pressures. While a subset of studies have looked at relationships between resin duct morphology and climate, these are almost exclusively from lower-elevation moisture-limited systems (Rosner and Hannrup, 2004; Gaylord et al., 2013, 2015; Hood et al., 2015; Rodríguez-García et al., 2015; Slack et al., 2016, 2017; Reed and Hood, 2020; Vázquez-González et al., 2021), with a gap in knowledge for trees in high-elevation, energy-limited systems. There is also a lack of uniformity in methods employed to disentangle anatomical changes related to geometric size-age relationships vs. climate and other factors across studies (Hood and Sala, 2015; Vázquez-González et al., 2020). Because mature trees typically experience increased mortality relative to smaller, younger trees during a bark beetle outbreak (Safranyik and Carroll, 2006; Boone et al., 2011; Audley et al., 2020), it is important to remove the influence of geometric size-age relationships before investigating climate-resin duct relationships and comparing such associations across trees of differing ages in natural, forested ecosystems. Longer time periods outside of the juvenile growth curve are optimal for estimating and separating geometric growth curves from climate-driven multi-decadal patterns and long-term trends in resin duct anatomy, particularly in long-lived tree species such as whitebark pine.

Based largely on the observed bioclimatic envelope of modern whitebark pine populations, species distribution models project that increased warming will amplify whitebark pine mortality through direct and indirect effects of drought, decreased tree resilience, and increased susceptibility to insects and pathogens (Millar et al., 2007, 2012; Chang et al., 2014; Wong and Daniels, 2017; Goeking et al., 2018; Millar and Delany, 2019; Goeking and Windmuller-Campione, 2021). However, while numerous studies have investigated climate-growth relationships in conifers, there are exceedingly limited data on how climate influences resin duct morphology, which is important for understanding growth-defense interactions and within-tree resource allocation patterns. Here we present an evaluation of high-elevation whitebark pine growth and defense in response to climate variability by aggregating mean growth and anatomical defense structure measurements across two U.S. northern Rocky Mountain sites spanning 1924 to 2015. Our objectives were: (1) to evaluate relationships between age- and growth-related trends in resin duct morphology (duct production, size, area, density, relative area) with climatic variability (i.e., seasonal and annual temperature and moisture) and (2) to determine if regional whitebark pine growth and resin duct morphology is constrained by changing climatic conditions, including warming temperatures and regional drying. From past work assessing climate controls on tree growth (Fritts, 1971; Wettstein et al., 2011), we expect growth will be positively related to warmer growing season temperatures and regional drought if they are energy-limited trees (Kipfmueller and Salzer, 2010), in contrast to moisture-limited systems where tree growth correlates negatively with warmer and drier conditions. In contrast to previous research from mostly moisture-limited systems (Gaylord et al., 2013, 2015; Vázquez-González et al., 2021), we also speculate that with increased growth, whitebark pine trees will produce more and larger resin ducts, due to more favorable growing season conditions facilitating increased cambial cell division, resource acquisition, and carbohydrate accumulation.

2. Materials and methods

2.1. Site description

Data for this study were collected across two high-elevation (1,905–2,120 m a.s.l.) whitebark pine sites on the Flathead Indian Reservation as part of a larger fire history reconstruction project for the Confederated Salish and Kootenai Tribes (CSKT; Figure 1A). We selected sites where whitebark pine was an important component of the canopy (>15% canopy cover). The Boulder (BLD) site (1,905 m a.s.l.) is located in the western Mission Range on the eastern edge of Flathead Lake (Figure 1B). Mean annual temperatures are generally cool, ranging from 1 to 6°C with a 20- to 70-day frost-free period. This area receives more precipitation than the Three Lakes (3LK) site, with an average of 1,397 mm annually occurring primarily as snowfall from November through March. The Three Lakes (3LK) site (2,120 m a.s.l.) is located approximately 74 kilometers away in the Reservation Divide Mountains between the Flathead Indian Reservation and Missoula, Montana (Figure 1C). Mean annual temperatures at 3LK are similar to BLD, ranging from 2 to 5°C with a 30- to 60-day frost-free period. Precipitation averages 1,143 mm per year, which occurs primarily as snowfall from November through March.

Figure 1. Flathead Indian reservation in northwestern Montana with the two sites (A) and a detailed map of the Boulder (BLD) study site (B) and three Lakes (3LK) study site (C). Permanent plots (blue circles) were located on a 200 m grid as part of a larger fire history reconstruction project for the confederated Salish and Kootenai Tribes. Plots cover approximately 100 ha at 3LK and 88 ha at BLD.

Forest stands at both sites are relatively open-canopy, mixed-conifer forests consisting of lodgepole pine (P. contorta var latifolia) and whitebark pine, with smaller components of subalpine fir (Abies lasiocarpa) and Engelmann spruce (Picea engelmannii). Evidence of disturbance legacies, including wildfires, insect, and pathogen outbreaks are common at both sites. Whitebark pine mortality, due to cumulative impacts from mountain pine beetle and white pine blister rust, was prevalent, with the majority of trees (∼82%) being dead or dying at the time of sampling. The broader region has been affected by numerous large-scale bark beetle outbreaks, occurring in the 1930s, 1960s–1980s, and most recently 2002–2009 (Arno and Hoff, 1989; Kipfmueller and Swetnam, 2002; Buotte et al., 2017; van de Gevel et al., 2017), while white pine blister rust was introduced to this region of the northern Rocky Mountains circa 1950 (Mielke, 1943; Samman et al., 2003; Geils et al., 2010). Despite high observed whitebark pine mortality, whitebark pine comprises a relatively low percentage of overall forest composition across our study sites (∼15–30%) and stand/canopy structure has remained largely unchanged post-mortality due to rugged topography and the relatively open-grown conditions of these high-elevation forests.

2.2. Data collection and preparation

As part of the CSKT fire history project, a 200 m² grid was overlain across each site using a combination of remote sensing and field surveying for optimal placement (Figures 1B, C). Within the 200 m² grid, individual plots (hereafter “macroplots”) were located 200 m apart along cardinal directions, with a total of 38 plots for the 3LK site (covering approximately 100 ha) and 34 plots for the BLD site (covering approximately 88 ha) (Kichas et al., 2020). Each macroplot was spatially referenced with Garmin Rhino 650 GPS units and plot centers were permanently marked and tagged. From the macroplot center, four 10 m wide, 100 m long belt transects were established in each cardinal direction (north, south, east, and west) (Supplementary Figure 1). The closest ten trees within each 100 m transect were sampled (to a maximum distance of 100 m) with a diameter threshold of 15 cm diameter at breast height (DBH) (Heyerdahl et al., 2014). Data collected for each tree at each macroplot included species, condition (live or dead), size (DBH), canopy base height, crown class (dominant, codominant, intermediate, or suppressed), tree height, and evidence of fire (e.g., fire scarring and/or presence of charcoal on tree stem). Using Haglöf 4.3 mm diameter increment borers, we collected increment cores from 40 sample trees per macroplot near the base of each tree (at or below 15 cm from the root collar), as we required cores that were as close to the pith as possible to accurately determine tree establishment dates for the fire history project. For this study, we utilized tree-ring data from 58 live whitebark pine trees at the 3LK site and 54 live whitebark pine at the BLD site for a total of 112 trees. The majority of trees across the two sites established in the early 20th century, with the mean size (diameter at breast height) of whitebark pine trees sampled being 21.2 cm (± 0.57 cm SE), with a mean age of 102 years (± 1.5 years SE).

Increment cores were air-dried and processed using standard dendrochronological techniques (Stokes and Smiley, 1996). We measured ring-widths to the nearest 0.001 mm using a Velmex measuring stage interfaced with the recording software Measure J2X (Voortech Consulting 2005). We crossdated tree-ring series visually and graphically using ring-width patterns and frost rings (LaMarche and Hirschboeck, 1984). Crossdating was facilitated using site-specific chronologies from ring-width measurements, which were compiled as part of previous studies (Kichas et al., 2020, 2021) and statistically checked using the program COFECHA (Holmes et al., 1986).

Individual tree cores were scanned at high resolution (1200 d.p.i.) on an Epson V550 flatbed scanner and annual resin duct features were measured across all years of growth using the program ImageJ (version 1.46r, National Institutes of Health, Bethesda, MD, USA). All resin ducts were measured to the nearest 1 mm × 10^–7mm using the ellipse tool and the calendar year in which they formed was documented. For defense metrics we followed the methods outlined by Hood and Sala (2015) and collected measurements across five categories, including non-standardized defense measures of resin duct size (mean size of all resin ducts per annual ring; mm²), duct production (number of ducts year^–1), and total duct area (sum of resin duct size: mm² year^–1), as well as resin duct measures that have been standardized to annual tree growth, including duct density (number of resin ducts mm^–2 year^–1) and relative duct area (% annual ring).

2.3. Data analysis

All analyses were conducted using the statistical package R (v. 3.5.0) (R Development Core Team, 2008). We used multivariate, correlation, and regression analyses to examine relationships between radial growth, resin duct morphology and climate. We first investigated the potential influence of local fire activity on whitebark pine tree growth and resin duct morphology, using spatially explicit 20th century wildfire records, to determine dates of fire events adjacent to our study sites. We then conducted superposed epoch analysis (SEA) to evaluate whether there were differences in growth or defense in response to documented regional fire activity. Results from the SEA suggested no significant differences in growth and defense to wildfire activity, allowing us to evaluate the relationships of tree growth and resin duct anatomy to climatic conditions without considering fires (Supplementary Figures 2, 3). The sample plots were also not located directly in or below any avalanche slide paths, nor were site slope angles conducive for generating damaging or destructive snow avalanche events. Accordingly, we did not detect any synchronicity in growth reduction, cambial scarring, reaction wood or resin duct development (e.g., traumatic resin ducts) in our sampled trees that would indicate a widespread avalanche event (or active insect outbreaks) that could obscure the relationships between tree growth, resin duct development, and climate (DeRose et al., 2017; Peitzsch et al., 2021). We measured growth and defense metrics over the full record (1908–2015) but limited our climate analysis to a period of 1924–2015 to both increase our sample depth (n = 80 trees at 1924) and reduce inclusion of growth and defense data over the earliest part of the record when juvenile growth vigor may obscure climate relationships. The expressed population signal for the individual and combined BLD and 3LK tree-ring and resin duct chronologies was above 0.85, limiting additional error attributable to decreasing sample number through time to <15% (Wigley et al., 1984).

For climate analyses, tree growth and resin duct data were conservatively and interactively detrended using the dplR package (version 1.7.2) (Bunn, 2008) with the goal of removing geometric age-related growth trends while retaining as much interannual to decadal+ scale variability, as well as preserving any long-term increasing growth trends (if present) as possible. The interactive detrending allowed us to deal with the infrequent and unique growth curves or responses that arise due to topography, soils, or individual tree disturbance events. For sites with energy-limited tree growth, the use of detrending methods that allow for positive slopes (e.g., flexible splines or smoothers) on most or all the series would result in the removal of an increasing trend in growth. Thus, each tree-ring and resin duct series was individually fit with either a negative exponential curve (79.7%), Hugershoff curve (1.6%), Friedman super smoother (5.1%), or a spline with a 0.5 frequency response at 75% the series length to remove the geometric growth trend (13.6%) (Fritts, 1976; Supplementary Table 1). The resulting series were then standardized by dividing by the fitted trend line, transformed into unitless ratio (index), and then combined into a composite chronology by calculating the robust weighted mean across all samples (Cook and Peters, 1997).

Since no single detrending method is perfect, and different methods excel at preventing different problems or potential biases, we performed a sensitivity test of our analyses by detrending using several different methods. Though positive series-length trends in growth will be removed, we also detrend our data using (1) individual splines configured to 2/3 of the individual series length (0.5 frequency cutoff) as well as (2) a 75-year (∼75% of our 100-year series) spline (0.5 frequency cutoff) fit to all series. Both methods are fairly conservative in preserving interannual to decadal+ scale variability, and perform well in reducing variability in ring-width series arising from dense forest dynamics and biased end effects from overly inflexible detrending methods (Klesse, 2021). Other than the expected removal of a positive trend in growth (i.e., ring width index, RWI), the major interannual to decadal+ features of each growth and defense metric chronology were highly similar (r > 0.8) (Supplementary Figure 4). As such, we feature our iteratively detrended growth and defense chronologies for use in subsequent analyses, with results from the sensitivity tests provided in Supplementary material.

Estimates of local hydroclimate within the study areas were obtained at a spatial resolution of 4 km² from the Wolock and McCabe (2018) monthly water-balance model (MWBM) output (Data available).¹ The MWBM uses precipitation and temperature input data from the Parameter-elevation Regressions on Independent Slopes Model (PRISM, from the PRISM Climate Group, Oregon State University)² (Daly et al., 2008) to estimate monthly evapotranspiration (potential and actual), runoff, snowfall, and soil moisture storage. We opted to use this version of the PRISM-based water balance dataset because (1) instrumental weather data were not available at these study sites (Daly et al., 1994; Kipfmueller and Salzer, 2010; Abatzoglou et al., 2017), (2) its accuracy has already been evaluated at watershed scales across the contiguous U.S. (McCabe and Wolock, 2011), and (3) both PRISM and the MWBM are commonly utilized dataset for many dendroecological and climatic studies, particularly in areas with complex topography (Millar et al., 2007; Kipfmueller and Salzer, 2010; Pederson et al., 2013; Guiterman et al., 2015; Hood et al., 2015; Carnwath and Nelson, 2016; Ferrenberg et al., 2017; Cansler et al., 2018; Goeking and Izlar, 2018; Martin et al., 2019, 2020).

To obtain the best representation of climatic conditions affecting trees in our study we averaged climate data from the nearest four grid points for each site. We grouped monthly climate data into seasonal and annual windows to evaluate potential responses of trees to seasonal variability. We also included climate data from the year prior to growth to account for potential lagged influences of climate on tree growth and duct morphology. Monthly climate variables obtained for this analysis included temperature, precipitation, vapor pressure deficit (VPD), soil moisture storage, as well as the Palmer Drought-Severity Index (PDSI) from the western Montana climate division in the WestWide Drought Tracker product (Abatzoglou et al., 2017) (Data available),³ and April 1 snow-water equivalence (SWE) (Pederson et al., 2011) (Data available).⁴ We included the PDSI data because we wanted to investigate possible relationships of whitebark pine growth and defense to regional (Montana climate division 1) warming and drying patterns, in addition to the local, site-specific climate data obtained through the MWBM. In addition, PDSI is a measure of long-term drought, integrating the current and lagged effects of precipitation, temperature, soil moisture storage and evapotranspiration over multiple months. This aspect of PDSI generally corresponds well with the physiological response of tree growth to changing environmental conditions.

To determine climate variables most strongly associated with tree growth and resin duct anatomy, we performed regression “best” subset selection using the LEAPS package (version 3.1) (Lumley, 2020). The best subsets regression procedure selects for an optimized subset of predictor variables from a larger pool of candidate variables by comparing all possible model combinations and using objective criterion, such as reduced mean square error and increased proportion of variance explained, to identify top-performing models. For our model building procedure, we utilized Mallow’s Cp, the coefficient of determination (R²), and Bayesian information criterion (BIC) to inform the most parsimonious selection of predictor climate variables, and in addition, top models were run through a forward-backward stepwise procedure to ensure the elimination of overfit models. Variables were allowed to enter the model with an alpha-to-enter significance level α_E = 0.15 and removed from the model with an alpha-to-remove significance level α_R = 0.1.

We also partitioned climate variables into upper and lower quartiles (e.g., Supplementary Figure 5) and used the non-parametric Wilcoxon Rank-Sum Test to evaluate potential differences in whitebark pine growth and resin duct morphology across the extremes of the respective distributions of climate variables related to warmer and drier, or cooler and wetter conditions. For this analysis, we included warm season (May-September) average temperature, soil moisture, and maximum vapor pressure deficit, cool season (October-April) precipitation, April 1st snow-water equivalence, and annual average regional Palmer Drought Severity Index.

3. Results

3.1. Growth and defense increase with regional warming and drying

Assessments of relationships between tree size, growth, and resin duct morphology indicated the need to remove the biological influence of early vigorous tree growth from tree anatomical defense metrics using dendrochronological detrending methods prior to performing climate analyses (Figure 2 and Supplementary Figure 6). Relationships between the detrended whitebark pine growth and defense metrics and the upper and lower quartiles of climate variables show that increased temperature, reduced precipitation, drier soils, and increased vapor pressure deficit collectively drive increased tree growth and production of anatomical resin defenses (Figure 3 and Supplementary Figure 7). Significant climate-related differences in tree growth, resin duct production, duct size, and duct area indicate that more favorable growing conditions and/or extended growing season length allowed for both increased growth and defense in these high-elevation whitebark pine trees (Figure 3 and Supplementary Figure 7).

Figure 2. Time series plots of raw tree-ring data (top row), detrended data (middle), and duct metrics that have been standardized to tree growth (bottom row). Tree growth (ring-widths: mm), resin duct production (no. ducts year^{– 1}), duct size (mm²), duct area (mm² year^{– 1}), and resin duct metrics that have been standardized to tree growth: duct density (no. ducts mm^{– 2} year^{– 1}) and relative duct area (% annual ring). Smoothed lines in the time series plots represent 5-year moving averages.

Figure 3. Comparison of detrended whitebark pine tree growth and resin duct morphology across the upper and lower quartiles of reported climate variables from 1924 to 2015: annual average palmer drought severity index (PDSI), warm season (May-September) average temperature (°C), maximum vapor pressure deficit (VPD), soil moisture storage, cool season (October–April) precipitation (mm), and April 1st snow-water equivalence (SWE). The blue- and red-colored plots indicates significance at the p < 0.05 level. Reported p-values are from the non-parametric Wilcoxon’s Rank-Sum Test.

The relationship between raw (non-detrended) data shows positive or neutral correlations between whitebark pine tree size (diameter at breast height), radial growth, and resin duct morphology, suggesting that trees at our sites are both growing well and are well-defended (Table 1). In contrast, the negative correlations with resin duct metrics that have been standardized to tree growth (resin duct density and relative duct area), reflect a dilution effect, as relative duct area and duct density geometrically decline with increased tree size even if there is no actual decrease in duct production and/or total resin duct area.

Table 1. Spearman correlation coefficients for raw, non-detrended whitebark pine tree size (diameter at breast height; DBH), tree growth (ring-width: mm), resin duct production (no. ducts year^–1), duct size (mm²), duct area (mm² year^–1), and resin duct metrics that have been standardized to tree growth: duct density (no. ducts mm^–2 year^–1) and relative duct area (% annual ring).

3.2. Whitebark pine growth and resin duct morphology have unique climate responses

Utilizing the detrended datasets in a “best” subsets regression framework to explore the major climate drivers of growth and defense, we found that climate explains a significant proportion of variation in whitebark pine tree growth and resin duct morphology, with R² values of models ranging from 20 to 31% (Table 2). Cool-season (October–April) total precipitation is significantly negatively correlated with all detrended growth and defense metrics (Table 2 and Figure 4), which corresponds to increased growth occurring under regionally drier (more negative) PDSI values (Figure 3). In addition to cool-season precipitation, average spring temperature and spring PDSI are positively related to whitebark pine tree growth, while annual average maximum vapor pressure deficit is positively associated with resin duct production and duct size (Table 2 and Figure 4). Annual average soil moisture is also positively associated with resin duct size while annual average temperature is positively related to resin duct area. For the resin duct metrics that were standardized to tree growth (duct density and relative duct area), we found negative relationships with fall PDSI and positive associations with annual total precipitation and the previous year’s snow-water equivalence (Table 2 and Figure 4). Our sensitivity tests using the spline detrended chronologies yielded similar results, though the influence of temperature is generally reduced as expected (Supplementary Figure 7 and Supplementary Table 2). In general, regardless of detrending method, we observed increased whitebark pine growth and defense in response to warming (increased temperature and vapor pressure deficit) and regional drying (decreased precipitation, soil moisture storage, and snow-water equivalence corresponding with more negative PDSI values) across our study sites (Figure 3 and Supplementary Figure 7).

Table 2. Top climate models predicting tree growth (ring width index, RWI), detrended duct production/size/area, and standardized duct metrics (duct density and relative duct area) as determined by best subsets regression.

Figure 4. Scatterplots of whitebark pine tree growth, detrended resin duct production, duct size, and duct area, as well as duct metrics standardized to tree growth (resin duct density and relative duct area) and the primary climate variables derived through the “best” subsets regression framework (Table 2). PDSI, palmer drought severity index; VPD, vapor pressure deficit; SM, soil moisture; PPT, precipitation; SWE, soil-water equivalence; RWI, ring width index.

4. Discussion

4.1. Growth and defense increase with regional warming and drying

Unlike many lower-elevation, moisture-limited forests, high-elevation whitebark pine trees in the U.S. Northern Rocky Mountains have experienced increased growth and produced more and larger resin ducts with greater total duct area under warmer and regionally drier conditions (Figures 3, 4 and Table 2). As documented with previous research, drier, lower-elevation forests experience recent warming and drying as a stressor, resulting in decreased growth and defense (Allen et al., 2010, 2015; Anderegg et al., 2015; DeSoto et al., 2020). In contrast, at our high-elevation, energy-limited systems, warming and drying appear to release the physical constraints of snowpack and the generally cool-wet conditions that can limit photosynthetic and meristem activity restrictions on tree growth (Korner, 2015). Increased warming in these high-elevation, subalpine systems results in earlier snowmelt and ground thaw, which increases biologically available water and provides temperatures favorable for xylogenesis and nutrient exchange (Danby and Hik, 2007; Chhin et al., 2008). We speculate that increased temperature allows increased secondary xylem cell division, which, coupled with adequate soil moisture, leads to increased radial growth (Korner, 2015; Cabon et al., 2020). Warmer temperatures may also increase photosynthesis that can positively influence growth, although this relationship is weak and assimilation is less constrained by temperature than photosynthesis (Cabon et al., 2022). Temperature is also positively correlated with resin duct formation, although the effects of precipitation are more nuanced, with slight water deficits increasing duct production while severe droughts lower production (see review by Vázquez-González et al., 2020). It is likely that a variety of mechanisms could be working to collectively drive the observed increase in whitebark pine growth and defense (Cabon et al., 2022; Dow et al., 2022). For instance, increased biologically available water early in the growing season, facilitated by earlier snowpack melt and runoff (Pederson et al., 2009, 2011, 2013; Barichivich et al., 2014), could promote cambial cell division and expansion (Cabon et al., 2020, 2022). Decreased precipitation (and presumably increased insolation due to less cloud cover) could also facilitate increased spring and autumn photosynthetic activity, promoting greater carbon sequestration annually, resulting in the observed increases in tree growth and resin duct defenses (Körner, 1998; Goldblum and Rigg, 2005). Our results align with Runyon et al. (2022), one of the only other climate studies of pine defenses in high-elevation systems, which found the highest concentration of both constitutive and induced terpene defenses in five-needle pines occurred at the coldest and driest sites. In these high-elevation systems, whitebark pine trees don’t appear to be experiencing substantive moisture stress in response to recent regional warming and drying, although this dynamic may change with continued warming and drying. More research on the environmental controls on resin duct formation is needed.

Our results also provide evidence showing age as a primary control on resin defenses with greater duct production and area over the juvenile growth years (Figure 2 top; DeAngelis et al., 1986). This supports field observations, whereby young or juvenile trees typically experience lower mortality during insect outbreaks compared to larger, more mature trees (Safranyik and Carroll, 2006; Boone et al., 2011; Audley et al., 2020), although the specific mechanisms modulating these interactions remain unclear (e.g., beetle larvae may have poorer or simply less feeding substrate in younger, less mature trees and/or younger trees may have increased flexibility in the allocation of resources to resin-based defenses relative to older individuals). Accordingly, to assess the potential secondary climate drivers of growth and defense, detrending is required (Figure 2 and Supplementary Figure 6). Individual series, biological growth curve detrending methods (e.g., negative exponential curves, negative to zero slope regression lines, and rigid splines) designed to retain decadal- to centennial-scale variability and long-term trends in tree-ring records were successfully utilized in removing the early juvenile growth curve present in our data and producing detrended, standardized indices of growth and defense suitable for use in climate analyses (Figure 2 and Supplementary Figures 4, 6).

The degree to which plants allocate resources toward defenses relative to growth remains an area of active debate with a number of competing models put forth to conceptualize these complex physiological processes (Stamp, 2003). While our study is not designed to test any of these hypotheses directly, the observed positive relationships between tree growth and resin duct defenses we found in whitebark pine, as well as enhanced resin chemistry reported for other high-elevation conifers (Runyon et al., 2022), suggest no apparent tradeoffs but rather increases in both growth and defense with warming conditions in these energy-limited systems. While we did detect negative relationships between tree growth and the resin duct metrics that have been standardized to tree growth (duct density and relative duct area), this is in agreement with other studies (Kane and Kolb, 2010; Ferrenberg et al., 2014; Hood and Sala, 2015; Moreira et al., 2015; Mason et al., 2019; Zhao and Erbilgin, 2019) and reflects a dilution effect, as the metrics standardized to tree growth geometrically decline with increased tree size, even if there is no actual decrease in resin duct production and/or total resin duct area. Given the observed positive increase to both tree growth and defense, we propose that warmer temperatures, leading to an extended growing season, reduces constraints to meristem activity and allows these high-elevation trees to simultaneously allocate resources toward both tree growth and resin duct development, with no discernible growth-defense tradeoffs (Table 1, Figure 3 and Supplementary Figure 7). While similar patterns of increased tree growth under warmer and regionally drier conditions have been reported for other conifers in high-elevation, energy-limited systems throughout western North America (Graumlich and Brubaker, 1986; Peterson et al., 2002; Salzer et al., 2009, 2014; Miyamoto et al., 2010; Millar et al., 2020; Kemp-Jennings et al., 2021), these results are the first to also show enhanced anatomical resin duct defenses.

4.2. Whitebark pine growth and resin duct morphology have unique climate responses

The combination of resin duct characteristics collectively influences a trees ability to produce, store, and mobilize oleoresin. We found that key climate variables (cool-season precipitation, annual maximum vapor pressure deficit, spring PDSI and average temperature) collectively explain a considerable proportion of variability in resin duct anatomy (duct production, duct size, and duct area) (Table 2, Figure 4 and Supplementary Figure 2) after the primary control (size-age relationship) has been removed through detrending (Figure 2 and Supplementary Figure 6). Together, our findings suggest that these high-elevation trees are responding to climatic controls in nuanced and complex ways by actively adjusting components of tree growth and resin duct morphology in response to changing environmental conditions.

Prior research investigating growth and defense relationships in whitebark pine trees growing at our study sites found that individuals that persisted through late 20th century mountain pine beetle outbreaks had significantly greater resin duct defenses than trees that died during these outbreaks (Kichas et al., 2020). This research also found that increased winter precipitation negatively correlated with tree growth and resin duct area in live whitebark pine trees. Together, and supported by our findings here, this suggests that warming and regional drying has created environmental conditions favorable for enhanced photosynthesis, meristem activity and simultaneous resource allocation to both growth and defense. However, the degree to which these enhanced defenses may mitigate against future projected increases in the magnitude/frequency of beetle outbreaks, due to enhanced thermal suitability in high-elevation whitebark pine forests, remains unclear. Also, while duct anatomy is important for tree survival during bark beetle outbreaks, resin chemistry and resin viscosity are critical components modulating tree-insect-pathogen interactions (Kichas et al., 2021). Sampling temporal variation in resin chemistry, viscosity, and flow is not possible in a retrospective analysis in this study; however, relatively small increases in resin duct area would lead to a greater capacity to deliver oleoresin to insect attack sites, as resin flow is proportional to the fourth power of resin duct diameter in accordance with Hagen-Poiseuille’s law (Schopmeyer et al., 1954; Vázquez-González et al., 2020; Yi et al., 2021). Increasing temperatures also correlate with reduced resin viscosity and increased resin flow in some species, which may further enhance resin mobilization in whitebark pine under warmer conditions (Zas et al., 2020; Wermelinger et al., 2021), as well as increase terpene concentrations (Runyon et al., 2022). In addition, warming temperatures stimulate the production of ethylene, a phytohormone that stimulates the differentiation of resin duct epithelial cells from tracheid development in annual growth rings of trees (Yamamoto and Kozlowski, 1987; Franceschi et al., 2002; Hudgins and Franceschi, 2004; Vázquez-González et al., 2020). This could further prime resin duct defenses of whitebark pine and other conifers as conditions become warmer and regionally drier in these high-elevation environments.

We found that warmer temperatures and regionally drier conditions are associated with increased tree growth, duct production, duct size, and duct area (Table 2 and Figures 3, 4). Ultimately, these high-elevation, energy-limited systems experience warming and regional drought as a release from cold and wet conditions that limit photosynthetic and meristematic activity. As these high-elevation ecosystems continue to experience increased warming, energy-limited whitebark pine trees appear poised to be better defended against insects and pathogens via increased vigor and enhanced resin duct anatomy (Tables 1, 2 and Figures 3, 4). These metrics of resin-based defensive structures are important characteristics influencing tree mortality to insect-pathogen complexes (Kane and Kolb, 2010; Hood et al., 2015; Kichas et al., 2020). However, the benefit of enhanced growth and defense in response to warming and drying conditions is likely a non-linear relationship with critical thresholds at which further warming and drying would eventually result in these systems becoming moisture-limited with subsequent physiological implications related to moisture stress (Andrus et al., 2021; Pedlar et al., 2021). Also, since climate-controlled production of defense structures is a secondary response relative to the primary size-age control, there are likely thresholds of tree size beyond which the additional production of defensive structures during drought offers limited additional benefit.

Whitebark pine mortality across our study sites was high prior to sampling (∼82% mortality) due to extensive mountain pine beetle activity in the broader region throughout the 20th century. Trees sampled for this analysis are part of a surviving cohort and, as such, our results may not accurately capture climate-resin duct relationships for phenotypes beyond the region of our analysis. While the whitebark pine trees sampled for this study have responded to these biotic and abiotic pressures with increased resin duct defenses, it is unknown how future disturbance synergies may push these systems past a threshold at which regional warming and drying may inhibit rather than enhance growth and defense responses. In addition, there are numerous other interacting factors that could influence tree growth and mortality, which may obscure these relationships. For instance, projected enhanced thermal suitability for mountain pine beetle might amplify outbreaks, potentially overwhelming even well-defended trees (Boone et al., 2011; Bentz et al., 2016; Reed and Hood, 2020). White pine blister rust can also cause extensive mortality in healthy, vigorous trees; however, the extent to which this pathogen affects resin duct formation and resin flow is unknown, but limited research suggests that resin-duct defenses are unrelated to genetic resistance of white pine blister rust (Holtz and Schoettle, 2018). Other physiological processes, such as resource allocation to reproduction (e.g., masting and cone production), can further complicate growth and defense relationships in mature trees. For instance, Redmond et al. (2019) found that resin duct area decreased during years of increased cone production in piñon pine (Pinus edulis), while no discernable tradeoffs with growth were detected. Growth form could also be an important factor for long-term whitebark pine persistence as research has demonstrated lower mortality rates (Maher et al., 2020) and differing growth responses to climate (Millar et al., 2020) in krummholz whitebark pine trees relative to full-stem trees growing at treeline. While our study was not explicitly designed to test for such interactions, future studies could investigate relationships between exposure to white pine blister rust and insect-pathogen complexes on resin duct development and within-tree resource allocation patterns, and to quantify resin duct characteristics of krummholz vs. full-stem whitebark pine trees.

5. Conclusion

We found that high-elevation whitebark pine trees in the northern Rocky Mountains have experienced increased growth and defense under warming temperatures, reduced precipitation and snowpack, drier soils and increased vapor pressure deficit. In response to these changing conditions, these subalpine forest ecosystems are experiencing a release from conditions that historically limited radial growth (e.g., short growing season, presence of snowpack into growing season). Contrary to previous research conducted in moisture-limited systems, our results show that regional warming and drying trends may enhance both growth and defense in the energy-limited subalpine forests where many whitebark pine stands occurs. This relationship highlights an important distinction between high-elevation systems and lower-elevation systems - where drought stress, as reflected by traditional climate variables and indices, beneficially translates to an extended annual growing season. We posit that a lengthened growing season in these high-elevation systems promotes accumulation of carbon and increased meristem activity, ultimately, increasing tree growth and defense structures. Of critical importance, we expect this positive relationship to cross a threshold at which even energy-limited high-elevation systems experience drought-stress and limited growth and defense as warming continues. The dynamic nature of these responses further emphasize how relationships with tree growth, resin ducts, and climate will vary in complex ways, within and across species, as well as over topo-climatic gradients (Mäkinen et al., 2002; Miyamoto et al., 2010; Ackerman and Goldblum, 2021). Future studies should investigate these patterns and consider climate-resin duct relationships across elevational transects to assess if and how these relationships change across important biophysical gradients and investigate at what point the benefits of regional warming and drying in these systems will become liabilities for subalpine forests.

Data availability statement

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

Author contributions

NK: conceptualization, investigation, data production, formal analysis, methodology, visualization, writing–original draft, and funding acquisition. GP: conceptualization, investigation, formal analysis, methodology, resources, and writing–review and editing. SH: investigation, methodology, and writing–review and editing. RE: project administration and resources. DM: methodology, supervision, funding acquisition, and writing–review and editing. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by grants from the Joint Fire Science Program Graduate Research Innovation (GRIN) Program (JFSP 17-2-0114), National Science Foundation (BCS-1539820/BCS-2029775/BCS-1832486), and USDA National Institute of Food and Agriculture (NIFA; 2013-3842420992).

Acknowledgments

We are grateful to the Confederated Salish and Kootenai Tribes, specifically Rich Jannsen of the Department of Natural Resources, Jim Durglo of CSKT Tribal Forestry, and Tony Incashola Jr. of CSKT Tribal Forestry, who allowed us to conduct this research on tribal lands. The Salish Cultural Committee, and Salish Kootenai College provided additional support. We also thank Robert E. Keane (U.S. Forest Service Rocky Mountain Research Station), Jeremy Littell (U.S. Geological Survey, Alaska Climate Adaptation Science Center), and two reviewers who provided valuable comments that greatly improved the manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

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.

The handling editor declared a past co-authorship with the author GP.

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.2023.1089138/full#supplementary-material

Footnotes

1. ^ https://doi.org/10.5066/F71V5CWN
2. ^ http://prism.oregonstate.edu
3. ^ https://wrcc.dri.edu/wwdt/
4. ^ https://www.ncei.noaa.gov/access/paleo-search/study/11205

References