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

## Cite this dataset

Crockett, Joseph (2023). Data from: Ability of seedlings to survive heat and drought portends future demographic challenges for five southwestern US conifers [Dataset]. Dryad. https://doi.org/10.5061/dryad.280gb5mt5

## Abstract

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.

## README: Ability of seedlings to survive heat and drought portends future range contractions for five southwestern US conifers

This readme file was generated on 2023-01-03 by Joseph L Crockett

### Author/Principal Investigator Information

Name: Joseph L. Crockett

Institution: University of New Mexico

Address: Department of Biology, University of New Mexico, MSC03 2020, Albuquerque, NM 87131-0001, United States

Email: jcrockett@unm.edu

### Associate or Co-investigator Contact Information

Name: Matthew D. Hurteau

Institution:University of New Mexico

Address: Department of Biology, University of New Mexico, MSC03 2020, Albuquerque, NM 87131-0001, United States

Email: mhurteau@unm.edu

### Date of data collection: April 2020 - September 2021

### Geographic location of data collection: New Mexico, CA

Information about funding sources that supported the collection of the data:

This work was supported by the Interagency Carbon Cycle Science program Grant No. 2017-67004-26486/ project accession no. 1012226 from the USDA National Institute of Food and Agriculture, Agriculture and Food Research Initiative program Grant No. 2021-67034-35106/project accession no. 1026366 from the USDA National Institute of Food and Agriculture, and the Joint Fire Science Program under Project JFSP 16-1-05-8 and Project JFSP 20-1-01-9.

### DATA & FILE OVERVIEW

### File list

ABCO_df.csv

ABCO_grid_maps.rda

abco_gridterra_values.rda

abco_gridterra.rda

abco_maca_df.rda

abco_weekly_models_gcb.rda

ABCO.maca_maps.rda

PIED_df.csv

PIED_grid_maps.rda

pied_gridterra_values.rda

pied_gridterra.rda

pied_maca_df.rda

pied_weekly_models_gcb.rda

PIED.maca_maps.rda

PIEN_df.csv

PIEN_grid_maps.rda

pien_gridterra_values.rda

pien_gridterra.rda

pien_maca_df.rda

pien_weekly_models_gcb.rda

PIEN.maca_maps.rda

PIPO_df.csv

PIPO_grid_maps.rda

pipo_gridterra_values.rda

pipo_gridterra.rda

pipo_maca_df.rda

pipo_weekly_models_gcb.rda

PIPO.maca_maps.rda

PSME_df.csv

PSME_grid_maps.rda

psme_gridterra_values.rda

psme_gridterra.rda

psme_maca_df.rda

psme_weekly_models_gcb.rda

PSME.maca_maps.rda

ready_data_frame.csv

d_all.csv

elev_slope_aspect_maca.rda

ibutton_df.csv

aznm.zip

### METHODOLOGICAL INFORMATION

### Description of methods used for collection/generation of data:

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 = 30m) and predicted survival probabilities for these locations. We extracted contemporary climate from GridMET (daily max temperature, precipitation total, and mean VPD, resolution = 4km, Abatzoglou 2013) and soil moisture from Terraclimate (total column soil moisture [mm/m, converted to cm^{3}/cm^{3}], resolution = 4km, 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 21^{st} 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: 4km, 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 cm^{3}/cm3], resolution = 4km, 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 4km 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%.

### Instrument- or software-specific information needed to interpret the data:

R version 4.1.2 (2021-11-01) 'Bird Hippie'

RStudio 2021.09.2 Build 382

### Package List:

dplyr, ggpubr, hexbin, raster, sp, , terra, Survival, brms, car, cowplot, ggplot2, grid, gtable, lubridate, readr, reshape2, sf, sjPlot, stringr, surrosurv, tidybayes, tidyr, wesanderson, broom

### DATA-SPECIFIC INFORMATION

### Specialized formats or other abbreviations used:

.rda files are an R specific data file that can incorporate multiple data structures as well as more rows/cases than .csv formats.

### People involved with sample collection, processing, analysis and/or submission:

Joseph L Crockett: Sample collection, processing, analysis and submission

### DATA-SPECIFIC INFORMATION FOR: ready_data_frame.csv

Number of cases/rows: 34/1814

Variable List:

SPECIES: Sample species

TEMP: Temperature (C)

MOISTURE: Initial soil moisture (%)

REP: Sample id

CHAMBER: Incubator number

GROUP: Sample group

WATER_WEIGHT: weight of initial water (g)

interval: interval period (days)

DATE_START: Start date

DRY_WEIGHT: Soil/container/sample weight (g)

id.x: Sample id

WEIGHT_actual: total weight (g)

DATE: Survey date

STATUS_FIXED: live/dead (0/1)

VPD_25: Vapor pressure deficit (kPa)

TEMP_actual: Chamber temperature (C)

actual_water: derived water weight (weight - dry weight, g)

minimum_weight: minimum measured weight (g)

WEEK.x: Survey week

tstart.x: interval start day (days)

tstop: interval end day (days)

tduration: interval length

event.x: live/dead (1/0)

TEMP_num: Temperature (C)

status_sum: weeks living

id.y: sample id

event.y: live/dead (0/1)

time: interval length (days)

tstart.y: interval start day

WEEK.y: interval start week

tstart: interval start day

event: live/dead (0/1)

WET_WEIGHT_actual: water weight (g)

actual_water_volumetric: water weight (cm3/cm3)

Missing data codes: NA

## DATA-SPECIFIC INFORMATION FOR: d_all.csv

Number of cases/rows: 36/683

Variable List:

SPECIES: sample species

TEMP: chamber maximum daily temperature (C)

MOISTURE: initial soil moisture (%)

REP: sample number

CHAMBER: incubator number

DATE_START: date of incubator start

WEIGHT_WEEK_1: weight (g)

DATE_WEEK_1: sample date, week one

STATUS_FIXED_WEEK_1: live/dead status (0/1)

WEIGHT_WEEK_2: weight (g)

DATE_WEEK_2: sample date, week two

STATUS_FIXED_WEEK_2: live/dead status (0/1)

WEIGHT_WEEK_3: weight (g)

DATE_WEEK_3: sample date, week three

STATUS_FIXED_WEEK_3: live/dead status (0/1)

WEIGHT_WEEK_4: weight (g)

DATE_WEEK_4: sample date, week four

STATUS_FIXED_WEEK_4: live/dead status (0/1)

WEIGHT_WEEK_5: weight (g)

DATE_WEEK_5: sample date, week five

STATUS_FIXED_WEEK_5: live/dead status (0/1)

TARGET_WEIGHT: gravimetric weight for initial soil moisture (g)

DRY_WEIGHT: weight (g)

FIELD_WEIGHT: drained weight (g)

WATER_WEIGHT: weight (g)

WEIGHT_WEEK_6: weight (g)

DATE_WEEK_6: sample date, week six

STATUS_FIXED_WEEK_6: live/dead status (0/1)

WEIGHT_WEEK_7: weight (g)

DATE_WEEK_7: sample date, week six

STATUS_FIXED_WEEK_7: live/dead status (0/1)

LIVE_WEEKS: total weeks living

DEAD_WEEKS: week of death

TEMP_MOISTURE: id by temp/moisture

id: sample id

GROUP: sample group

Missing data codes: NA

## DATA-SPECIFIC INFORMATION FOR: ibutton_df.csv

Number of cases/rows: 6/43897

Variable List:

datetime

GC: Incubator number

temperature: Temperature (C)

humidity: Relative humidity (%)

VPD: Vapor pressure deficit (kPa)

DATE: Date

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: ABCO_df.csv

Number of cases/rows: 12/237

Variable List:

id: Sample id

event: live/dead (0/1)

SPECIES: sample species

TEMP: temperature (C)

MOISTURE: initial soil moisture (%)

CHAMBER: incubator number

GROUP: sample group

time: interval duration

interval: interval start day

tstart: interval start day

tstop: interval end day

tduration: interval duration

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PIED_df.csv

Number of cases/rows: 12/512

Variable List:

id: Sample id

event: live/dead (0/1)

SPECIES: sample species

TEMP: temperature (C)

MOISTURE: initial soil moisture (%)

CHAMBER: incubator number

GROUP: sample group

time: interval duration

interval: interval start day

tstart: interval start day

tstop: interval end day

tduration: interval duration

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PIPO_df.csv

Number of cases/rows: 13/300

Variable List:

id: Sample id

event: live/dead (0/1)

SPECIES: sample species

TEMP: temperature (C)

MOISTURE: initial soil moisture (%)

CHAMBER: incubator number

GROUP: sample group

time: interval duration

interval: interval start day

tstart: interval start day

tstop: interval end day

tduration: interval duration

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PIEN_df.csv

Number of cases/rows: 12/245

Variable List:

id: Sample id

event: live/dead (0/1)

SPECIES: sample species

TEMP: temperature (C)

MOISTURE: initial soil moisture (%)

CHAMBER: incubator number

GROUP: sample group

time: interval duration

interval: interval start day

tstart: interval start day

tstop: interval end day

tduration: interval duration

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PSME_df.csv

Number of cases/rows: 12/248

Variable List:

id: Sample id

event: live/dead (0/1)

SPECIES: sample species

TEMP: temperature (C)

MOISTURE: initial soil moisture (%)

CHAMBER: incubator number

GROUP: sample group

time: interval duration

interval: interval start day

tstart: interval start day

tstop: interval end day

tduration: interval duration

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: abco_maca_df.rda

Number of cases/rows: 6/632765

Variable List:

x: longitude

y: latitude

s_min: minimum survival probability

s_sum: sum of survival probility

year: year

model: CMIP5 model

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: pipo_maca_df.rda

Number of cases/rows: 6/2599569

Variable List:

x: longitude

y: latitude

s_min: minimum survival probability

s_sum: sum of survival probility

year: year

model: CMIP5 model

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: pied_maca_df.rda

Number of cases/rows: 6/3484051

Variable List:

x: longitude

y: latitude

s_min: minimum survival probability

s_sum: sum of survival probility

year: year

model: CMIP5 model

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: pien_maca_df.rda

Number of cases/rows: 6/72304

Variable List:

x: longitude

y: latitude

s_min: minimum survival probability

s_sum: sum of survival probility

year: year

model: CMIP5 model

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: psme_maca_df.rda

Number of cases/rows: 6/1413200

Variable List:

x: longitude

y: latitude

s_min: minimum survival probability

s_sum: sum of survival probility

year: year

model: CMIP5 model

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: elev_slope_aspect_maca.rda

Number of cases/rows: 5/44252

Variable List:

X: longitude

Y : latitude

elevation: elevation (m)

slope: slope (°)

aspect: aspect (°)

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: ABCO.maca_maps.rda

Number of cases/rows: 3/41201

Variable List:

X: longitude

Y : latitude

ABCO: species presence

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PIPO.maca_maps.rda

Number of cases/rows: 3/41201

Variable List:

X: longitude

Y : latitude

PIPO: species presence

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PIED.maca_maps.rda

Number of cases/rows: 3/10426

Variable List:

X: longitude

Y : latitude

PIED: species presence

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PIEN.maca_maps.rda

Number of cases/rows: 3/41201

Variable List:

X: longitude

Y : latitude

PIEN: species presence

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PSME.maca_maps.rda

Number of cases/rows: 3/41201

Variable List:

X: longitude

Y : latitude

PSME: species presence

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: abco_gridterra_values.rda

Number of cases/rows: 4/4389

Variable List:

Year: year

x: longitude

y: latitude

s_min: minimum survival probability

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: pied_gridterra_values.rda

Number of cases/rows: 4/88303

Variable List:

Year: year

x: longitude

y: latitude

s_min: minimum survival probability

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: pien_gridterra_values.rda

Number of cases/rows: 4/115

Variable List:

Year: year

x: longitude

y: latitude

s_min: minimum survival probability

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: pipo_gridterra_values.rda

Number of cases/rows: 4/51001

Variable List:

Year: year

x: longitude

y: latitude

s_min: minimum survival probability

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: psme_gridterra_values.rda

Number of cases/rows: 4/22108

Variable List:

Year: year

x: longitude

y: latitude

s_min: minimum survival probability

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: abco_gridterra.rda

Number of cases/rows: 8/45,627,030

Variable List:

year: Year

x: longitude

y: latitude

tstop: period end day

tduration: period duration

TEMP: temperature (C)

VPD_25: vapor pressure deficit (kPa)

actual_water_volumetric: soil moisture (cm3/cm3)

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: pied_gridterra.rda

Number of cases/rows: 8/124,900,890

Variable List:

year: Year

x: longitude

y: latitude

tstop: period end day

tduration: period duration

TEMP: temperature (C)

VPD_25: vapor pressure deficit (kPa)

actual_water_volumetric: soil moisture (cm3/cm3)

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: pien_gridterra.rda

Number of cases/rows: 8/15,428,160

Variable List:

year: Year

x: longitude

y: latitude

tstop: period end day

tduration: period duration

TEMP: temperature (C)

VPD_25: vapor pressure deficit (kPa)

actual_water_volumetric: soil moisture (cm3/cm3)

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: pipo_gridterra.rda

Number of cases/rows: 8/114,118,710

Variable List:

year: Year

x: longitude

y: latitude

tstop: period end day

tduration: period duration

TEMP: temperature (C)

VPD_25: vapor pressure deficit (kPa)

actual_water_volumetric: soil moisture (cm3/cm3)

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: psme_gridterra.rda

Number of cases/rows: 8/74,598,660

Variable List:

year: Year

x: longitude

y: latitude

tstop: period end day

tduration: period duration

TEMP: temperature (C)

VPD_25: vapor pressure deficit (kPa)

actual_water_volumetric: soil moisture (cm3/cm3)

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: ABCO_grid_maps.rda

Number of cases/rows: 3/3124

Variable List:

X: longitude

Y : latitude

ABCO: species presence

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PIPO_grid_maps.rda

Number of cases/rows: 3/7813

Variable List:

X: longitude

Y : latitude

PIPO: species presence

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PIED_grid_maps.rda

Number of cases/rows: 3/8548

Variable List:

X: longitude

Y : latitude

PIED: species presence

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PIEN_grid_maps.rda

Number of cases/rows: 3/1056

Variable List:

X: longitude

Y : latitude

PIEN: species presence

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: PSME_grid_maps.rda

Number of cases/rows: 3/5110

Variable List:

X: longitude

Y : latitude

PSME: species presence

Missing data codes: NA

### DATA-SPECIFIC INFORMATION FOR: abco_weekly_models_gcb.rda

Includes four objects:

b_r2.list_abco: Coefficient of determiniation from model

cut_events_abco: knot locations for model

fit_BRM_spline_water_weight_VPD_abco: model, requires brms package and cut_events_abco renamed to cut_events (see Final_code_TS)

model.df_abco: dataset used for model

### DATA-SPECIFIC INFORMATION FOR: pied_weekly_models_gcb.rda

Includes four objects:

b_r2.list_pied: Coefficient of determiniation from model

cut_events_pied: knot locations for model

fit_BRM_spline_water_weight_VPD_pied: model, requires brms package and cut_events_pied renamed to cut_events (see Final_code_TS)

model.df_pied: dataset used for model

### DATA-SPECIFIC INFORMATION FOR: pipo_weekly_models_gcb.rda

Includes four objects:

b_r2.list_pipo: Coefficient of determiniation from model

cut_events_pipo: knot locations for model

fit_BRM_spline_water_weight_VPD_pipo: model, requires brms package and cut_events_pipo renamed to cut_events (see Final_code_TS)

model.df_pipo: dataset used for model

### DATA-SPECIFIC INFORMATION FOR: psme_weekly_models_gcb.rda

Includes four objects:

b_r2.list_psme: Coefficient of determiniation from model

cut_events_psme: knot locations for model

fit_BRM_spline_water_weight_VPD_psme: model, requires brms package and cut_events_psme renamed to cut_events (see Final_code_TS)

model.df_psme: dataset used for model

### DATA-SPECIFIC INFORMATION FOR: pien_weekly_models_gcb.rda

Includes four objects:

b_r2.list_pien: Coefficient of determiniation from model

cut_events_pien: knot locations for model

fit_BRM_spline_water_weight_VPD_pien: model, requires brms package and cut_events_pien renamed to cut_events (see Final_code_TS)

model.df_pien: dataset used for model

### DATA-SPECIFIC INFORMATION FOR: aznm

Includes shapefile for combined Arizona and New Mexico state boundaries

## 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 (cm^{3}/cm^{3}) 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 cm^{3}/cm^{3}], 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 21^{st} 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 cm^{3}/cm^{3}], 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%.

## Usage notes

Application: R/Rstudio

Package List: 'dplyr','ggpubr','hexbin','raster','sp','','terra','Survival','brms','car','cowplot','curl','ggplot2','grid','gtable','lubridate','readr','reshape2','sf'

'sjPlot','stringr','surrosurv','tidybayes','tidyr','wesanderson','broom','geomtextpath'

## Funding

National Institute of Food and Agriculture, Award: 2017-67004-26486/project accession no. 1012226

National Institute of Food and Agriculture, Award: 2021-67034-35106/project accession no. 1026366

Joint Fire Science Program, Award: Project JFSP 16-1-05-8

Joint Fire Science Program, Award: Project JFSP 20-1-01-9