Ability of seedlings to survive heat and drought portends future demographic challenges for five southwestern US conifers

Climate change and disturbance are altering forests and the rates and locations of tree regeneration. We examined seedling survival of five southwestern United States (US) conifer species found in warmer and drier woodlands (Pinus edulis, P. ponderosa) and cooler and wetter subalpine forests (Pseudotsuga menziesii, Abies concolor, and Picea engelmanii) under hot and dry conditions in incubators. We constructed models that explained 53% to 76% of the species-specific survival variability, then applied these to recent climate (1980-2019) and projected climate (1980-2099) for the southwestern US. We found that lower elevations within species’ range would have low survival under projected climate and that range contraction would be greatest for species that currently occupy warm-dry conditions. These results demonstrate that empirically derived physiological limitations can be used to identify where species composition or vegetation type change are likely to occur in the southwest US.

Methods
We obtained 135 one growth year old seedlings of each species (675 seedlings total) from the New Mexico State University John T. Harrington Forestry Research Center in Mora, NM between April 2020 and September 2021. Seedlings were grown from locally sourced seeds from mature trees in northern New Mexico to ensure that they represented local adaptation to conditions. Seedlings were grown in a greenhouse in 10 cm containers at staggered intervals to ensure that they were of similar age and size when they were placed in the incubators.

We transplanted seedlings into 22 cm deep pots (volume: 590 cm3) with well-drained soil (2 parts sphagnum moss, 1.5 parts vermiculate, 1.5 parts sand) brought to field capacity following transplanting. We allowed soil moisture to draw down to the treatment level, measuring soil moisture gravimetrically. Once soil moisture matched the treatment condition, we randomly assigned seedlings to one of two Percival Model E-36L1 incubators. We intended to use 15 seedlings per species per temperature/moisture combination, but several seedlings died during moisture drawdown resulting in several treatments using fewer than 15. Incubator rack positions were adjusted to ensure that seedlings in each incubator received equivalent photosynthetically active radiation (~260 mol). We programmed temperature treatments to follow a diurnal cycle with lower temperatures at night (15C) and progressive steps to treatment temperatures during the day.

We set photoperiods at 15/9 hours to reflect growing season conditions. Incubators controlled temperature and light. We placed iButtons (Model number DS1923; Temperature accuracy +/- 0.5C; humidity resolution 0.6%; https://www.maximintegrated.com/en/products/ibutton-one-wire/data-loggers/DS1923.html) within each incubator to record the actual temperature and humidity at hourly intervals. Temperatures matched the programmed values and humidity was highest at the start of each stage of the experiment and decreased as moisture was lost to evapotranspiration. We calculated Vapor Pressure Deficit (VPD, kPa) at each time step as the difference between saturated and effective water pressure of the air. We assessed seedling health weekly with visual assessments of needle coloration and by measuring leaf fluorescence with a MultispeQ v2.0 fluorometer (Guadagno et al. 2017). The efficiency of light adapted photosynthetic reaction centers (measured as a ratio of Fv* to Fm*) corresponds well to destructive measures of cell conductance yet provides a non-destructive, rapid assessment of plant death with greater accuracy than visual assessment of foliage color. We determined plant death as either 95% brown/grey needle coloration or below 0.1 Fv*/Fm*. At plant death, we recorded time between treatment start date and seedling death to express results in days until death.

Methods for processing the data:

To test the physiological tolerances of seedlings from a variety of climates, we subjected seedlings to temperature and soil moisture combinations ranging from those commonly found in burned landscapes to those projected with ongoing climate change. We used five species whose southwestern distributions range from warmer and drier woodlands to cooler and wetter subalpine forests (Supp. Fig. 1). *Pinus edulis* Engelm. is a widespread conifer in the southwestern US, considered drought-hardy and commonly found between 1370 and 2440 m (Burns & Honkala 1990). *Pinus ponderosa* Douglas ex C. Lawson has an extensive range in the western US, is fire-tolerant as an adult, and in the southwest and southern Rockies is found up to 3050 m (Burns & Honkala 1990). Due to a legacy of fire suppression and resultant forest densification in the southwest US, *Pseudotsuga menziesii* (Mirb.) Franco has colonized forests previously dominated by *P. ponderosa*, though is less fire tolerant, climate tolerant, and is generally found at a higher elevation range in the southwestern US (2440m to 3290m) (Burns & Honkala 1990). *Abies concolor* (Gord. & Glend.) Lindl. Ex Hildebr is found up to 3400 m in the central Rockies and is sensitive to heat and drought but generally tolerant of a range of soil conditions (Burns & Honkala 1990). *Picea engelmannii* Parry ex Engelm. is the least widely distributed species in the southwestern US of the five species we examined, occupies the coolest and wettest areas, and is found between 2400 m and 3700 m elevation (Burns & Honkala 1990).

Because heat and drought effects vary by species, we used a 3x3 full factorial design, with three levels of temperature (34°C, 39°C, 44°C (based on growing season air temperature measurements in a high-severity burn area of the 2011 Las Conchas fire in northern New Mexico and the maximum temperature limits of the incubators) and three levels of soil moisture (5%, 10%, 15% of soil moisture at field capacity, measured gravimetrically). We calculated VPD from chamber relative humidity to use as a predictor variable because it is an integrated measure of temperature and moisture, but because chambers were unable to control rH levels, VPD varied over time.

Data analysis:

To analyze the species-specific relationships between temperature, soil moisture, VPD, and time-to-death, we first used two-way ANOVAs to compare temperature and moisture treatments in R (R core team 2021) using a Type II sum of squares implemented in the car package (Smith & Cribbie 2014, Fox & Weisberg 2019). We converted soil moisture weights (g) to volumetric by calculating the ratio of moisture to soil volume (cm3/cm3) so that we could use models to examine projected climate with volumetric soil moisture. We then used a Bayesian framework to construct species-specific discrete-time proportional-hazard models in R with the brms package, which fits models using ‘Stan’ (Tutz & Schmid 2016, Bürkner 2017). These models allow for an event to be modeled if it occurs between regular observation intervals as well as incorporate time-varying covariates as predictors. These models present the hazard of an event (here, death) occurring. Models took the form of

Y_i  ~ bernoulli(μ_i)
logit (μ_i )  ~ a_ij+ β_1 x_1i+ β_2 x_2i  + β_3 x_3i+s(x_4i)
s(x_4i )  ~ β_4 x_4i+ z_k,for 1,…,k knots

a_ij  ~ Normal(0,4)
β_1  ~ Normal(1,1) 
β_2~Normal(-1,1)
β_3  ~Normal(1,1)
β_4  ~ Uniform(-inf,inf)
z_k  ~ Normal(0,σ_τ)
σ_τ  ~ Students-t(3,0,2.5)

with descriptions of coefficients and priors in table 1.

Where logit (μ_i )  is the logit of death occurring, a_ij  is the intercept of seedling j; β_1, β_2, and β_3   are the coefficients of temperature, initial soil moisture, and vapor pressure deficit, respectively; s(x_4i )  is the spline function for time since start, with coefficient β_4  and intercept z_k, for each 1:k knots. Errors have a Bernoulli distribution. Based on a literature search for likely effects of variables, we generated weakly informative, skeptical priors for each covariate (Table 1 and Supp. Table 1) and visually examined prior predictive distributions to ensure they generated realistic-looking data in the absence of observations. Models were fit with a Bernoulli family with a logit link and with a random intercept of plant ID to account for the repeated measures of each plant during the experiment. During model development, we determined that scaling and centering temperature, VPD, and initial soil moisture reduced divergent transition, and following scaling/centering these variables, we extracted the scale and center factors to apply to projected climates. We ran six chains with a 2000 iteration burn-in followed by 4000 iterations, and a thinning rate of 1, totaling 12000 post-warmup draws. We adjusted sampling algorithm settings (i.e., changing the adapt_delta value) where needed to achieve convergence of chains. To validate model performance, we conducted Gelman-Rubin diagnostic tests and checked that MCMC chain trace plots achieved stationarity and demonstrated mixing without autocorrelation between iterations. (Table 1, Supp. Fig 4). We then compared the posterior predictive distributions to the expected observations using the bayesplot package and Bayesian R2, which included both total variance explained and the marginal variance attributed to fixed effects, as well as calculated the root mean squared error (RMSE) from 10-fold cross validations. Following model assessment, we extracted 1000 posterior draws from the linear predictors and calculated the probability of surviving to time t given the hazard of an event:

Eq. 5    S(t)= exp(-∑0^t (μ(t〗))  S(t)=exp(-∑0^t μ_i ))  

In which the survival to time t is the exponentiated negative sum from 0 to time t of the hazard Ui. 

Present-day and future climate scenarios

To determine how present-day species ranges compare to the modelled survival probability, we extracted climate data from 1980 to 2019 for modeled species ranges from the National Individual Tree Species Atlas (Ellenwood et al. 2015, resolution = 30 m) and predicted survival probabilities for these locations. We extracted contemporary climate from GridMET (daily max temperature, precipitation total, and mean VPD, resolution = 4 km, Abatzoglou 2013) and soil moisture from Terraclimate (total column soil moisture [mm/m, converted to cm3/cm3], resolution = 4 km, Abatzoglou et al. 2018). We then calculated species presence using the modeled species ranges from the National Individual Tree Species Atlas as pixels with > 0 basal area and upscaled these data to match the resolution of GridMET.

From GridMET, we first calculated the pixel-wise precipitation-free period. For each day in the period, we calculated the mean of the daily temperature and VPD maximums for all days up to that day. We used this approach rather than calculating the mean of the entire period because a single mean daily temperature/VPD maximum for an entire precipitation-free period could obscure shorter heat waves or droughts that occur during that period. We scaled temperature, VPD, and soil moisture with the scale factors used to process data for our models and took 100 draws from the linear predictor to calculate the mean survival for each day of each year using eq 5. We chose 100 draws for projections to avoid computation limitations stemming from size of the area/days/years we analyzed. For each pixel, year, and species, we calculated the minimum survival value. We then calculated the annual percent area for each species’ range that exceeded our experimental thresholds (i.e., >34°C). We calculated the number of days that pixels in each bin experienced conditions likely to result in less than a 10% probability of survival. To determine whether area at risk or survival changed during the 1980-2019 period, we compared area at risk and mean survival between the 1980-1999 and 2000-2019 periods with T tests using a 0.05 significance level.

To examine how the modelled survival probability may change within present-day species ranges during the 21st century, we used Multivariate Adaptive Constructed Analogs (MACA) downscaled CMIP5 projections forced with the RCP8.5 emissions scenario to calculate the pixel-wise precipitation-free periods and mean daily maximum temperature for each period (MACAv2-METDATA, resolution: 4 km, daily, Abatzoglou & Brown, 2012). In lieu of projected soil moisture, which at the time of writing was not available at a similar scale as MACA, we incorporated monthly climatologies of Terraclimate that were calculated using a 4°C temperature increase (monthly normal, total column soil moisture [mm/m, converted to cm3/cm3], resolution = 4 km, Abatzoglou et al. 2018). We calculated daily temperature maximums and periods with less than 1mm of precipitation from five downscaled CMIP5 models (CCSM4, bcc-csm1-1-m, ACCESS1-3, GFDL-ESM2G, and CESM1-CAM5; Supp. Table 2). For each day in the period, we calculated the mean of the daily temperature and VPD maximums for all days up to that day. As with the present-day thresholds, we calculated the number of aggregate days that pixels in each bin experienced conditions likely to result in less than 10% probability of survival. Using elevation and slope values extracted from a 4 km DEM provided with the gridMET data, we calculated the mean elevation and median aspect per year per species of pixels in which survival probability is less than 10%.