Data from: Ability of seedlings to survive heat and drought portends future demographic challenges for five southwestern US conifers
Data files
Nov 29, 2023 version files 1.87 GB
-
ABCO_df.csv
7.27 KB
-
ABCO_grid_maps.rda
26.42 KB
-
abco_gridterra_values.rda
24.01 KB
-
abco_gridterra.rda
63.51 MB
-
abco_maca_df.rda
9.17 MB
-
abco_weekly_models_gcb.rda
30.57 MB
-
ABCO.maca_maps.rda
118.29 KB
-
aznm.zip
202.05 KB
-
d_all.csv
126.38 KB
-
elev_slope_aspect_maca.rda
818.96 KB
-
ibutton_df.csv
2.94 MB
-
PIED_df.csv
35.82 KB
-
PIED_grid_maps.rda
70.32 KB
-
pied_gridterra_values.rda
446 KB
-
pied_gridterra.rda
971.48 MB
-
pied_maca_df.rda
51.48 MB
-
pied_weekly_models_gcb.rda
30.57 MB
-
PIED.maca_maps.rda
84.48 KB
-
PIEN_df.csv
7.52 KB
-
PIEN_grid_maps.rda
9.12 KB
-
pien_gridterra_values.rda
1.10 KB
-
pien_gridterra.rda
24.74 MB
-
pien_maca_df.rda
1.08 MB
-
pien_weekly_models_gcb.rda
30.29 MB
-
PIEN.maca_maps.rda
106.74 KB
-
PIPO_df.csv
10.36 KB
-
PIPO_grid_maps.rda
59.37 KB
-
pipo_gridterra_values.rda
263.85 KB
-
pipo_gridterra.rda
20 B
-
pipo_maca_df.rda
37.63 MB
-
pipo_weekly_models_gcb.rda
29.72 MB
-
PIPO.maca_maps.rda
143.62 KB
-
PSME_df.csv
7.61 KB
-
PSME_grid_maps.rda
42.69 KB
-
psme_gridterra_values.rda
114.01 KB
-
psme_gridterra.rda
537.43 MB
-
psme_maca_df.rda
19.51 MB
-
psme_weekly_models_gcb.rda
30.86 MB
-
PSME.maca_maps.rda
127.13 KB
-
README.md
28.32 KB
-
ready_data_frame.csv
351.41 KB
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 cm3/cm3], 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 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: 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 cm3/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 (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%.
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'