We selected species and climate conditions local to the Jemez Mountains, of northern New Mexico, United States (36°N, 106°W), a region that has been the site of numerous studies on drought-and heat-related tree mortality (Allen and Breshears, 1998; Breshears et al., 2005; Rich et al., 2008; Breshears et al., 2009; McDowell et al., 2019). The landscape has also experienced extensive change from uncharacteristic high-severity wildfires over the past several decades (Allen, 2007; Stevens et al., 2021; Guiterman et al., 2022). The Jemez Mountains have a semi-arid continental climate: across the Jemez Mountains, mean annual precipitation (1991–2020 norm) is 564 mm with a largest proportion (214 mm) falling in the summer monsoon season from July to September (Hegewisch and Abatzoglou, 2022). Highest average summer temperatures (7.8 to 24.3°C) and lowest rainfall (27.4 mm) are in June, the most likely time for prolonged drought and heatwave conditions (Hegewisch and Abatzoglou, 2022). Forest soils are composed of volcanic rock and gravel with a loam texture, dating mostly to Pleistocene eruptions that created the dominant physiographic feature of the range, a prominent crater known as the Valles Caldera (Soil Survey Staff, 2002; Muldavin and Tonne, 2003). The elevational range of the Jemez Mountains extends from 1,800 m to 3,500 m, supporting ecological communities across elevations from piñon-juniper woodlands to dense spruce-fir forests (Johnson, 2001). Ecological communities across the mountain range have been exposed to an extensive history of low-intensity fires, documented by tree ring records across the Jemez (Dewar et al., 2021), before being interrupted in the late 1800s by a century of overgrazing and fire suppression (Margolis and Malevich, 2016). Fire deficit in this region and subsequent accumulation of woody fuels, in combination with drier and hotter climate, increased the risk of severe fires (Parks et al., 2015; Singleton et al., 2019). High-severity fires have characterized the landscape in recent decades, raising uncertainty about patterns of forest regrowth and recovery (Stevens et al., 2021; Guiterman et al., 2022).

To compare juvenile tree vulnerability as an indicator of forest recovery potential, we selected five study species that are dominant or co-dominant in mid-to upper-elevation forest communities across the Jemez Mountains. Species were selected both for their relative dominance within elevational zones and for their availability in tree nurseries in the western United States, which are a primary source of plant material for active reforestation management (Landis et al., 2010). As a result, the trees were not sourced from our specific study area, and we could not evaluate the climate conditions of the seed source as a factor in our study. Our selected species were Pinus ponderosa Douglas ex Lawson & C. Lawson (ponderosa pine), Pinus edulis Engelmann (piñon pine), Pinus flexilis E. James (limber pine), Pseudotsuga menziesii Mirbel Franco (Douglas-fir), and Picea engelmanni Engelmann (Engelmann spruce).¹ Within the Jemez Mountains, Picea engelmanni occurs most frequently in cooler, wetter sites at upper elevations of 3,000 meters (10,000 ft). Pinus flexilis and Pseudotsuga menziesii occur in slightly drier sites above 2,900 meters (9,500 ft). At the lowest forest elevational range, Pinus ponderosa and Pinus edulis occur in open forests and forest-woodland ecotones, at mean elevations of 2,590 meters (8,400 ft) (Muldavin and Tonne, 2003). Pinus edulis and Pseudotsuga menziesii were purchased from the University of Idaho’s Franklin H. Pitkin Forest Nursery in Moscow Idaho, United States. Picea engelmanni, Pinus ponderosa, and Pinus flexilis were purchased from the Colorado State Forest Service Nursery in Fort Collins Colorado, United States. Trees had well-established root systems and are consequently referred to in this study as juveniles, which represent individuals that have progressed beyond the seedling stage (first one to 2 years of initial growth following seed germination and initial root and shoot development), but prior to the sapling stage (established trees >1 m tall) (Fenner, 2000; National Park Service, 2003; Burrill, 2021). Juvenile trees used in this experiment averaged 4 ± 0.9 years old at the beginning of the experiment (Table 1).

We developed temperature profiles for each species to establish ambient temperature conditions for experimental growth chambers. To calibrate these, we used detailed 1987 vegetation surveys across the Jemez landscape (Allen, 1989) to determine a precise latitudes, longitudes, and elevations for each species. These latitudes and longitudes were then entered into Prism Explorer to retrieve our estimates of the daily mean, minimum, and maximum temperatures of those locations (Climate PRISM Group, Oregon State University); using all the location occurrences in that survey, for each of the species, gave us acceptable estimates of daily temperature profiles that were mimicked in the growth chambers. The accepted 30-year historical baseline that PRISM weather data are modeled from werebetween the years 1970 to 2000, with a spatial resolution of 4 km, for the month of June, historically one of the warmest months for this landscape (Bowen, 1990). We defined our lower-elevation ambient temperature regime where Pinus ponderosa and Pinus edulis co-occur as warm, dry forest types in the Jemez. Where cooler, wetter forests support a mix of Picea engelmanni, Pseudotsuga menziesii, and Pinus flexilis, we defined our higher-elevation ambient temperature regime (Bowen, 1990). Experimental temperatures were set to fluctuate between minimum and maximum temperatures, following a diurnal pattern, and heatwaves were calculated as a 10°C increase over these ambient temperature patterns (Supplemental Figure S1).

Three temperature-controlled plant growth chambers (Conviron, Model E7/2, Pembina North Dakota, United States) were set to temperatures that matched the PRISM data for each species or species group (Table 2). The growth chambers maintained 12 h of visible light (400 to 700 nm wavelengths) at a total output of 7,180 lumens from Agrobrite compact fluorescent light bulbs in each 24-h diurnal cycle. The fourth chamber was used to implement a one-week heatwave for a subset of plants in each chamber (staggered in time as detailed below in Experimental Design). The high elevation group was divided into two chambers, one containing Pinus flexilis and Pseudotsuga menziesii, and the other containing Picea engelmanni. We calibrated chamber temperatures using two Kestrel weather instruments (Kestrel 5,500, Nielsen-Kellerman, United States) placed on either end of each chamber to measure temperature accurately and ensure consistency across chambers. Our Kestrel measurements show the average lower-elevation temperatures ranged from 11.3°C (min) to 27.6°C (max), while the average higher-elevation temperatures ranged from 4.9°C (min) to 19.4°C (max). During the experimental heatwave of 10°C increase, the average lower-elevation temperatures ranged from 18.7°C (min) to 36.3°C (max), and the average higher-elevation temperatures ranged from 13.6°C (min) to 29.1°C (max). Both chambers achieved near identical temperature fluctuations and heatwave increases, with average temperatures deviating 0.3–0.5°C from the target temperatures. We maintained the heatwave chamber to +10°C (actual range 8.9–10°C) above each species’ ambient temperature.

We hypothesized that juveniles of different species would die at varying rates under climate conditions of prolonged drought, and that mortality rates would increase when exposed to an acute heatwave. Timing, duration, and magnitude of heatwave are all variables of potential interest but must be constrained within a single study. Our heatwave parameters (+10°C above ambient temperature for 7 days) are meant to test the extremes of modeled heatwave probabilities in the Southwest (Guirguis et al., 2018). The duration and intensity of our heatwave is consistent with experimental propositions to test heatwaves effects, severe enough to elicit responses in plant mortality but within the realm of realistic probability for expected heatwave events (Perkins-Kirkpatrick and Lewis, 2020; Breshears et al., 2021). We selected a magnitude of 10°C above ambient temperatures with a duration of 1 week, which has been proposed as a target for initial studies of heatwaves to enable cross-study comparisons (Breshears et al., 2021). Currently, heatwaves of this magnitude are occurring across elevations in the Southwest (Zachariah et al., 2023). We implement the heatwave 7 weeks into the study, which approximates a typical interval between peak soil moisture in late spring and when hotter summer temperatures begin to occur in the study area (e.g., for soil moisture dynamics in Pinus edulis woodlands; see Figure A1 in Breshears et al., 2009). This heatwave study expands upon previous studies that tested heatwaves of 5–7°C increase over ambient for 5–6 days (Notarnicola et al., 2021; Marchin et al., 2022). We expand our scope beyond these prior studies to explore heatwave effects on higher-elevation species with lower ambient temperatures (Table 2).

We acquired 50 individuals per species as a target sample size, plus additional plants to account for mortality during transportation and transplanting. All plants arrived in April 2021 and were placed into a temporary electric growth chamber (Conviron, Model CMP3000, Pembina North Dakota, United States) set to a constant 15°C with 12 h of fluorescent light per 24-h diurnal cycle. Plants were watered regularly to maintain non-limiting soil moisture. All individuals were moved from the growth chamber to the greenhouse in May 2021 with a daytime high temperature of 29°C and a nighttime low temperature of 10°C. After 1 month of acclimation to the greenhouse, the plants were repotted into 10 × 10 × 35 cm pots in June 2021. Soil was mixed to roughly replicate an intermediate textured soil type, in terms of sand, silt and clay composition, of the Jemez Mountains, using sandy loam soil (Acme Sand and Gravel, Tucson, AZ, United States) sieved to a gravel size of 2 mm and amended with 30% pure sand and 25% sphagnum peat moss based potting soil (Sunshine Mix #4, Sungro Horticulture, Agawam, MA, United States). To begin the experiment, plants were placed into the experimental growth chambers (Conviron, Model E7/2, Pembina North Dakota, United States) in August and September 2021, using a staggered-in-time design as described below.

To implement the experiment, 50 individual plants per species were assigned randomly to two treatment groups: a temperature treatment, and a water treatment. Within the temperature treatment group, 25 individuals per species were assigned randomly to absence or occurrence of a heatwave. Plants remained under ambient temperature conditions within their growth chamber, except that individuals within the heatwave treatment were exposed to a 10°C increase above ambient temperature for 7 days during week seven of the experiment before being returned to ambient temperature conditions. In each temperature treatment, the majority (n = 20) received no additional water for the duration of the experiment, while a subset of individuals (n = 5) was watered to saturation every 2 weeks to enable us to confirm that drought was driving juvenile mortality.

We used a staggered-in-time design to expose each species to a heatwave during week seven, given our restriction that only a single chamber was available for implementing heatwaves. The start of the experiment was staggered by 2 weeks between each of the other three chambers, meaning that “Week 1” of the experiment corresponds to different calendar dates for each chamber. This allowed sufficient time and space to prepare the heatwave chamber for each round of heatwaves. Plants were placed into the growth chambers in August and September 2021. All plants were watered and fertilized at the beginning of the experiment (Week 0), and measurements began on Week 1. Fertilizer was applied as instructed: about 29.6 grams of Liquinox All-Purpose Fish Fertilizer (NPK 5–1-1) mixed with 3.78 liters of water. Each plant was watered to saturation and then 0.3 liters of fertilizer mixture was sprayed onto the soil with a chemical sprayer. During Week 3, containers with 297 grams of calcium chloride desiccant (Damp Rid, W.M. Barr, Memphis, TN, United States) were used to offset the humidity increase associated with watering treatments. Two desiccant containers per species were placed in each chamber, monitored weekly and refilled as needed. At Week 5, the foliage of all plants was saturated with insecticidal soap (Bondie Products LLC, Oriskany, NY, United States), widely used in greenhouse and growth chamber experiments, to treat a mild occurrence of scale insect (Adelges sp.) [Hemiptera: Adelgidae] mostly confined to the Pinus ponderosa juveniles.

Tree size was characterized by taking single measurements of height (mm), basal diameter (mm), and dry weight biomass (g); these were recorded over the course of the experiment (Table 1). Age was estimated using branch whorl counts (Husch et al., 2003). Height was measured from the soil surface to the apical bud, and basal stem diameter was measured with a metric caliper at the base of each tree. Biomass was measured at the end of the experiment once all trees had died. Individuals were clipped at the soil surface, placed into labeled paper bags, set in a drying oven (Yamato, Model DX 402, Santa Clara, California, United States) set at 75°C for 24 h, and weighed.

We took weekly measurements capturing three key potential indicators of plant stress: (1) water availability, (2) stomatal conductance, and (3) leaf browning assessed using ocular estimates and digital photography.

First, to measure water availability and soil water loss, we weighed plants gravimetrically (grams) each week (Eamus, 2006). Weight data was used to calculate the amount of water available to the plant each week, graphically represented as a curve showing water weight in grams over time. Water stress was calculated as an inflection point along this curve, reflecting an estimate of when water availability transitions from relatively non-limiting to relatively limiting (Van Genuchten, 1980). Inflection points were calculated using R to compute the second derivative of the curve to find changes in concavity. Curves were inspected visually to verify the date of the inflection point.

Second, we measured stomatal conductance (mmol m⁻² s⁻¹) weekly on the uppermost needles of a subset of juveniles (n = 3 watered, 10 droughted) using a leaf porometer (Model SC-1, Meter group, Pullman, Washington, USA) as an indicator of plant transpiration. Low conductance measurements indicate reduced plant transpiration (Marchin et al., 2022). Because conductance measurements of non-conducting tissue (dry filter paper) recorded a minimum value of 90 mmol m⁻² s⁻¹, we assumed that plants had reached permanent stomatal closure when conductance was consistently at or below this minimum value.

Third, to establish a mortality date, we visually inspected leaf color twice per week to record leaf browning. Species tolerance to mortality drivers can be detected observationally by leaf browning trends, reflecting chlorophyll degradation and loss of photosynthetic capacity (Myneni et al., 1995; Jong et al., 2012). As a metric of mortality, foliar color has been shown to relate to hydraulic failure (Blackman et al., 2019). Ocular estimates of 90% brown thresholds have been used as a benchmark to determine plant death in other mortality studies (Adams et al., 2009, 2017). We categorized percent brown as 10%, 25%, 50%, 75%, or 90% brown. Studies have also found that leaf browning among Pinus species can lag behind hydraulic failure-induced mortality, showing that leaf color change is a delayed indicator of tree death (Blackman et al., 2019; Hammond et al., 2019). We similarly observed that many plants retained their green color even while exhibiting other signs of severe tissue degradation. To best establish time of death, we supplemented our percent brown estimates with criteria related to needle loss and brittleness (needles that broke under slight pressure). Plants were declared dead when they reached a 90% threshold for percent brown, needle loss, or brittle needles. We continued to measure individuals for 2 weeks after their death date to confirm each mortality event; the experiment concluded when all droughted plants died, given that watered plants were not expected to die. As a visual record of plant death, photographs were taken weekly against a grey background in the same position and with similar lighting (white balance of 3,200–3,400 K) using a digital camera (Sony Alpha 6,000 mirrorless, 24.3-megapixel 16–50 mm zoom lens). Further image analysis might enable more consistent, quantitative, repeatable, and precise metrics to assess color continuum change from green to brown (e.g., Hammond et al., 2019).

We analyzed significant differences in time to mortality within each species using two-sample t-tests in two contexts. First, we verified our control group to ensure that juvenile death was attributed to drought rather than chamber conditions, testing for differences between watered and droughted juveniles within each species. Second, we tested for differences in heatwave effects among droughted juveniles within each species. We visually inspected for normality and used Shapiro–Wilk tests to compare raw and transformed data; we determined that untransformed data met the assumptions of normality (Shapiro and Wilk, 1965). We performed statistical tests with and without outlier data points; because outliers did not change statistical conclusions, we included outliers in the final analysis. Using a two-sample t-test, we evaluated the null hypothesis that juveniles within a species would respond equally to drought, with or without heatwave exposure. The alternative hypothesis is that juveniles with exposure to a heatwave would have more rapid mortality rates. Because juveniles with heatwave exposure were expected to die more quickly on average, a one-sided value of p was used to evaluate significance (alpha = 0.05).

We used a one-way ANOVA to test significance of time to mortality among species. We compared individuals within drought-only treatment groups, and separately compared individuals within drought-heatwave treatment groups. We compared sample sizes and variance of data to determine that our data fit ANOVA assumptions. Data were analyzed using three measures of time to mortality, as indicators of drought stress: time to mortality (in weeks) since (1) water restriction, at the beginning of the experiment; (2) water stress, calculated from weight data at the water loss inflection point; and (3) permanent stomatal closure, when conductance data consistently reached a minimum. Across these three measures, we evaluated the null hypothesis that all species would respond equally to drought and have the same time to mortality. The alternative hypothesis is that at least one species mortality rate would differ. Given that mortality rates differed among species, further post-hoc analysis of species differences were tested using Tukey’s honest significant difference (HSD). We graphed predicted survival probability against observed values using the R package survival to produce a Kaplan–Meier curve, a nonparametric estimation of survival (Kaplan and Meier, 1958).

All statistical analyzes and visualizations were completed using R Statistical Software, version 4.2.1 (R Core Team, 2022; Therneau, 2022) and figures were produced using the package ggplot2 (Wickham, 2016). All code is available on publicly available GitHub pages.²

Our experimental focus was on drought effects with and without a heatwave, so we used our set of watered juveniles as a control group to confirm drought effects as the driver of mortality. We followed the progression of juvenile death until 100% of droughted individuals reached mortality. In the time that all droughted individuals had died, watered juveniles of four species (Picea engelmannii, Pinus edulis, Pinus flexilis, and Pseudotsuga menziesii) had 0% mortality. One of five watered Pinus ponderosa individuals died in both the ambient and heatwave groups. The number of juveniles that died during the experiment was statistically different between watered and droughted Pinus ponderosa juveniles (two-sample t-test, t₄₈ = −12.39, one-sided p < 0.0001). Collectively, these results validate that conditions within the growth chambers were suitable for juvenile growth if water was not limiting.

We focused our primary analysis on species responses to prolonged drought. Under ambient temperatures and drought conditions, average time to mortality differed among species (analysis of variance F-test, F_4,94 = 256.6, p < 0.0001). Pinus ponderosa tree juveniles died most quickly (ӯ = 10.7 weeks, SD = 1.71, n = 20), followed by Pinus edulis (ӯ = 14.6 weeks, SD = 2.17, n = 19), Picea engelmannii (ӯ = 19.2 weeks, SD = 1.51, n = 20), Pseudotsuga menziesii (ӯ = 19.4 weeks, SD = 2.31, n = 20), and Pinus flexilis (ӯ = 30.3 weeks, SD = 2.39, n = 20). All species were significantly different from each other, except for Picea engelmannii and Pseudotsuga menziesii (post-hoc Tukey HSD test; Figure 1). Using a Kaplan Meier survival function, we plotted a predicted survival probability curve against observed values for each species (Figure 1). These contrasts did not change when analyzing species under the drought-heatwave treatment (post-hoc Tukey HSD test; Supplemental Figure S2; Supplemental Table 1A), indicating that drought was the primary driver of mortality rather than heatwave exposure.

Our secondary analysis was focused on species responses to a weeklong heatwave event overlapping the experimentally imposed drought. Among all species, the mean time to mortality was earlier among juveniles with heatwave exposure (0.3 to 2.7 weeks earlier than drought-only juveniles), although this difference was not statistically significant for either of the lower-elevation species (Pinus ponderosa, Pinus edulis; Figure 2).

Two of three higher-elevation species (Pinus flexilis, Pseudotsuga menziesii) exposed to a heatwave had a statistically significant shorter time to mortality (Table 3; Figure 2). Pinus flexilis individuals with a heatwave treatment (ӯ = 27.6, SD = 5.05, n = 20) died an average of 2.7 weeks before those with no heatwave treatment (ӯ = 30.3, SD = 2.39, n = 20; two-sample t-test, t₃₈ = −2.16, one-sided p = 0.018), and Pseudotsuga menziesii juveniles with a heatwave treatment (ӯ = 17.4, SD = 2.39, n = 20) died an average of 2.0 weeks before those with no heatwave treatment (ӯ = 19.4, SD = 2.31, n = 20) (two-sample t-test, t₃₈ = −2.69, one-sided p = 0.005).

Across all species, we found that water stress was the earliest stress indicator, followed by stomatal closure 2 to 4 weeks later (Figures 3, 4). Under ambient temperatures and drought conditions, time to mortality from water stress differed significantly among species (analysis of variance F-test, F_4,94 = 103.1, p < 0.0001, post-hoc Tukey HSD test). Following the water stress inflection point, Pinus ponderosa juveniles died most quickly (ӯ = 6.5 weeks, SD = 1.77, n = 20), followed by Pinus edulis (ӯ = 10.5 weeks, SD = 2.89, n = 19), Picea engelmannii (ӯ = 11.7 weeks, SD = 1.41, n = 20), Pseudotsuga menziesii (ӯ = 11.9 weeks, SD = 1.45, n = 20), and Pinus flexilis (ӯ = 19.0 weeks, SD = 2.15, n = 20). Using a Kaplan Meier survival function, we plotted a predicted survival probability curve against observed values for each species (Supplemental Figure S2). The order species died did not change when adjusted for water stress, although Pinus edulis was not significantly different when compared to Picea engelmannii and Pseudotsuga menziesii (post-hoc Tukey HSD test; Supplemental Table 1B). By normalizing water stress, we expected species differences to diminish between lower-elevation species and higher-elevation species, anticipating that lower-elevation species would show greater drought tolerance. This effect is shown for Pinus edulis but not Pinus ponderosa, which maintained significantly earlier death dates compared to all other species.

Under ambient temperatures and drought conditions, the time to mortality from stomatal closure also differed significantly among species (analysis of variance F-test, F_4,45 = 27.1, p < 0.0001). Pinus ponderosa tree juveniles died most quickly (ӯ = 5.0 weeks, SD = 1.83, n = 10), followed by Pinus edulis (ӯ = 6.8 weeks, SD = 2.44, n = 10), Picea engelmannii (ӯ = 7.0 weeks, SD = 1.94, n = 10), Pseudotsuga menziesii (ӯ = 8.1 weeks, SD = 1.29, n = 10), and finally Pinus flexilis (ӯ = 14.2 weeks, SD = 2.86, n = 10). Using a Kaplan Meier survival function, we plotted a predicted survival probability curve against observed values for each species (Supplemental Figure S3). The order in which species died did not change when adjusted for stomatal closure stress, although Pinus ponderosa became non-significantly different compared to Pinus edulis and Picea engelmannii (post-hoc Tukey HSD test; Supplemental Table 1C). By normalizing for stomatal closure stress, we expected species differences to diminish even further. This effect is shown for Pinus edulis and Pinus ponderosa as they became more indistinguishable from Picea engelmannii and Pseudotsuga menziesii. However, these four species remained significantly different from Pinus flexilis after we adjusted for metrics of stress onset.

Notably, species contrasts did not change when analyzing species under the drought-heatwave treatment (post-hoc Tukey HSD test; Supplemental Tables 1B,C). Except for stomatal closure stress comparisons between Pinus ponderosa and Pinus edulis, all significant species differences under drought conditions were also significant under heatwave conditions. This further indicates that drought was the primary driver of mortality rather than heatwave exposure.

Across all levels of analysis–with and without adjusting to time of onset of stress–Pinus ponderosa died significantly earlier, and Pinus flexilis died significantly later, than most other species (post-hoc Tukey HSD test; Supplemental Tables 1A,C). Adjusting the time to mortality by stress points or stomatal closure did not change the order in which species died. The species most vulnerable to drought was Pinus ponderosa, followed by Pinus edulis, Picea engelmannii, Pseudotsuga menziesii, and Pinus flexilis.

The ability of juvenile trees to survive prolonged drought is a key determinant of forest regeneration capacity, informing future forest composition under extreme climate conditions and following severe disturbances (Tepley et al., 2017; Stevens-Rumann and Morgan, 2019; Stevens-Rumann et al., 2022). We implemented an experimental drought which lasted until all individuals had died, exploring the bounds of how long juvenile trees can persist under conditions of extended drought compounded by a heatwave. Our study indicates substantial differences in species vulnerabilities to drought at the juvenile life stage. The length of a lethal drought for lower-elevation species under warmer ambient conditions was short: 11 weeks for juvenile Pinus ponderosa and 14 weeks for Pinus edulis. Higher-elevation species under cooler ambient conditions did not reach mortality until much later: 20 weeks for Pseudotsuga menziesii and Picea engelmannii, and 30 weeks for Pinus flexilis.

Droughts of this duration are within the range of climatic episodes that have occurred in recent decades (e.g., Late 2001–2002; Breshears et al., 2009), and are projected for an altered-climate future. A 10-week drought strong enough to induce widespread mortality of Pinus ponderosa in the Southwest occurred approximately 30 times in the last century (Figure 2B in Adams et al., 2017). Yet shorter droughts are more frequent than longer droughts, and a 30-week drought, such as needed to induce widespread mortality of our Pinus flexilis species, might occur once a century (Figure 2B in Lauenroth and Bradford, 2009; Adams et al., 2017). Thus, species such as Pinus ponderosa may be more vulnerable to drought-induced mortality due to the greater likelihood of shorter droughts, especially at the juvenile life stage. Notably, the drought conditions leading up to the die-off event documented in our study area for adult trees in 2002 (Breshears et al., 2005) had anomalously low precipitation for almost all 52 weeks of that year (Figure A1 in Breshears et al., 2009), highlighting the potential for a protracted drought as in our experiment. Droughts within a single year may additionally be compounded by interannual drought such as in the Southwest U.S., which is currently in a 23-year megadrought matching both the length and duration of the worst drought in the past 1,200 years (Williams et al., 2022).

Heatwave impacts on juvenile trees in our study were relatively minor, yet consistent. Across all species, the average time to mortality was faster among individuals with heatwave exposure. For two high elevation species–Pinus flexilis and Pseudotsuga menziesii–this difference was statistically significant, accelerating average time to mortality by 2 to 3 weeks. Despite observations of heatwaves having initiated large tree mortality events (Ruthrof et al., 2018), crown dieback of trees (Matusick et al., 2018), and reduced productivity (Ciais et al., 2005; Young et al., 2022), this study found that the impact of a heatwave was small relative to the duration of the drought, perhaps because some individuals were already experiencing physiological stress from drought at that point. We initiated our heatwave at week seven during the dry-down process; at this time point, only Pinus ponderosa had progressed to stomatal closure, and the three higher-elevation species had not reached the soil moisture inflection point. Thus, our results suggest that the timing and duration of a heatwave could determine its effect on time to mortality among species (Breshears et al., 2021). A heatwave alone did not trigger mortality in watered juveniles.

Hotter temperatures during drought pose a triple threat to plants: amplification of atmospheric drought (or increased vapor pressure deficit), intensified soil drought, and direct effects of heat stress (Hammond et al., 2022). Related research on survival of juveniles in a post-fire landscape of the Jemez Mountains found that temperature-related metrics such as Heat Load Index were particularly influential in predicting seedling survival among species including Pinus ponderosa, Pinus edulis, and Pseudotsuga menziesii (Marsh et al., 2022). Other studies report that heatwaves do not directly damage leaf tissue, but instead intensify drought conditions (De Boeck et al., 2010; Salomón et al., 2022) which is consistent with our findings.

Differences in time to mortality among species at the juvenile stage are a complex emergent response resulting from a multitude of factors. Three intertwining factors–those which vary concurrently across elevation–deserve attention as to how they played a role in this study: (1) species-specific vulnerability to embolism during drought, (2) the associated ecohydrological response (for example, how that climate and the juvenile canopy structure affect soil evaporation), and (3) the climate conditions to which juvenile trees are exposed. Our experiment was not designed to isolate these three factors; rather, they inherently covary as part of our design and merit discussion.

Embolism in response to protracted drought is a key driver of plant mortality (Choat et al., 2012; Adams et al., 2017; McDowell et al., 2019, 2022). Drought-stressed plants may undergo cavitation in their xylem cells, disrupting their water transport system (Choat et al., 2018; McDowell et al., 2022). A common metric of drought tolerance is Ψ₅₀, the water potential of a plant when 50% of its water transport system has been disrupted by cavitation (Choat et al., 2012; Klein, 2014). Species with a more negative Ψ₅₀ value can tolerate more extreme drought conditions before losing hydraulic function. Research on Ψ₅₀ values suggests that, among our study species, Pinus ponderosa would have the greatest vulnerability to drought, followed by Pseudotsuga menziesii, Pinus flexilis, Pinus edulis, and finally Picea engelmannii (Domec and Gartner, 2001; Choat et al., 2012). In our study, Pinus ponderosa had the shortest time to mortality, consistent with this prediction. However, time to mortality among the other four species did not match this expected embolism order, even when we adjusted our time to mortality based on stress metrics. A potentially co-occurring driver of mortality among plants is carbon starvation via stomatal closure (McDowell et al., 2008, 2022; Peltier et al., 2021). In this study, stomatal closure followed water stress points by 2 to 4 weeks across all species. The time from stomatal closure to mortality, on the other hand, varied widely across species. While most species reached mortality 5 to 8 weeks following stomatal closure, Pinus flexilis survived for 13 weeks. These differences in time to mortality could possibly be attributed to higher tolerance to carbon limitation or even greater carbon reserves to withstand longer periods of stomatal closure, although our study was not designed to elucidate this mechanism. The early mortality in Pinus ponderosa is consistent with mechanisms of both embolism and carbon starvation, but leaves open questions about the mechanisms of drought and temperature vulnerability in the other species at the juvenile stage.

Ecohydrological factors–how plant structure interacts with environmental conditions to influence water availability–may help to elucidate the role of evapotranspiration in this study (Rodriguez-Iturbe, 2000; D’Odorico et al., 2010). A field study in the Jemez of planted seedling survival under shrub canopy found that planting efforts were more successful in cooler microclimates where denser shrub canopy provided more shading for seedlings and reduced vapor pressure deficit (Crockett and Hurteau, 2022). These authors found that shading provided by denser shrub canopy compared to more open shrub canopy delayed the time to mortality among Pinus ponderosa, Pinus edulis, and Pseudotsuga menziesii by 2 weeks. Other studies have found that increased tree basal area is correlated with increased mortality, suggesting that thinning forests may reduce resource competition among adult trees (Bradford and Bell, 2017; Bradford et al., 2022; McCauley et al., 2022), although it is unclear how juvenile trees may interact with this competition given differences in plant-water availability at various soil depths (Breshears et al., 2009). More open tree canopy architecture can allow greater evapotranspiration and water loss, while closed canopy structures can shade soil and retain moisture (Royer et al., 2012). This effect may have been operating in our study at the individual tree scale, so that leaf arrangement interacted with evaporative loss in our potted plants. In our study, Pinus flexilis was observed to have more compact needles with branches concentrated near the soil, which may have led to greater water retention. Even after soil moisture became more limiting, lower evaporation rates may have extended the time Pinus flexilis was able to survive.

Current species elevational distribution, represented in this study with corresponding temperature ranges, substantially influenced differing mortality rates among species. We expected a priori that lower-elevation species growing in warmer and drier environments would be more adapted to withstand drought stress. Instead, we found that warm ambient temperatures led to high evapotranspiration rates and water loss, followed by more rapid species death. Cooler ambient temperatures of higher-elevation species, conversely, acted as a buffer against drought conditions presumably by slowing the loss of water through evaporation; longer droughts were necessary before higher-elevation species would reach water stress, perhaps overriding expected differences in drought tolerance due to tree physiology.

All three interacting factors of species physiology, ecohydrological influences, and temperature effects could influence our experimental results. Nonetheless, the outcomes of this study indicate the potential consequences of prolonged drought and heatwave exposure at the juvenile state across species, adding to our overall understanding of these species’ drought tolerances at critical demographic life stages.

We acknowledge several limitations in this preliminary study. First, we used growth chambers to implement a controlled heatwave event and normalize drought stress across species to compare stress tolerance, which would be difficult to implement in a replicated field setting. Second, as is typical with growth chamber experiments, individuals within each chamber may be spatially correlated in their responses, and results can be difficult to translate to the landscape scale. Third, conductance measurements are sensitive to fluctuations in humidity and temperature; to maintain measurement accuracy, we calibrated the porometer prior to each measurement period. Fourth, ecohydrological influences with landscape topography, microclimate, and the broader plant community can attenuate climate drivers of heat and drought, influencing water stress among species. Fifth, as is often the case in reforestation efforts, plant material did not originate from the specific study region, and regional species adaptations were therefore not reflected in this study. While there is mixed evidence about the effect of seed sources from different local climates on survival rates of planted seedlings (Pinto et al., 2023), local variation could be important when making predictions for how species might respond across elevational gradients in specific landscapes. Recent studies have suggested that epigenetic responses (heritable traits in juveniles in response to stresses experienced by the parent organism) may help mitigate localized stress to drought heat stress (Carneros et al., 2017; Amaral et al., 2020; Backs and Ashley, 2021). Overall, we employed temperature conditions within the typical elevation range for ecologically co-occurring species the Jemez Mountains, and this study design was appropriate for imposing consistent environmental conditions across species.

Potential limitations in our experimental design suggest promising future lines of research. More mechanistic methods to complement our insights into drought tolerance would require a larger sample size to allow for weekly destructive sampling (e.g., of biomass, plant water potential and non-structural carbohydrates) or rewatering plants at different time intervals to estimate when plants cannot recover (e.g., Hammond et al., 2019). Estimates of plant death and percent brown are consistent with other plant mortality studies as relative measures, yet these are still whole-plant subjective measurements with coarse definitions. We mitigated observer variation by having one technician consistently evaluating tree mortality metrics. Multispectral or hyperspectral image data include a wider and more finely resolved array of wavelengths beyond what would be detected by ocular estimates, such as near infrared wavelengths that are longer than the visible wavelengths. Although results have not been reported in this study, initial progress has been made to develop and test an objective, quantifiable, and repeatable method of image analysis using weekly photographs tracking plant color change.

To uncover the impacts of heatwave at various stress points, future growth chamber studies could adjust heatwave timing, duration (days), and intensity (degrees of increase) after the onset of drought. A differently parameterized heatwave experiment (e.g., Breshears et al., 2021) would help to explore the dimensions of heatwave impacts. Our drought experimental framework could be expanded to include re-watering of trees following varying intervals of drought, to examine the recovery potential of individuals. These results could be complemented with further studies of the effects of warming baseline temperatures on species vulnerability (e.g., Adams et al., 2009, 2017) for a more complete picture of future forest change. Additionally, exposing trees of multiple species to a range of ambient temperatures could help to differentiate proximate responses from the tolerance to environmental conditions. In this experiment, the co-occurrence of species in a single climate represents a narrow view of the area a species may occupy. For instance, the co-occurrence of Pinus ponderosa and Pinus edulis represents the lower elevational boundary of Pinus ponderosa and the upper elevational boundary of Pinus edulis. A more complete factorial experiment in which species are exposed to a gradient of conditions would further elucidate species tolerances across their elevational range.

Drought and temperature regimes, and widespread disturbances, will almost certainly have a prominent role in shaping future species distributions (Thuiller, 2004; Loehman et al., 2018; Parks et al., 2019; Davis et al., 2023). Landscape-scale droughts and periodic heatwaves are likely to become persistent climate factors in the Southwest (Guirguis et al., 2018; Williams et al., 2020). The backdrop of drought, even within the context of species-specific climatic adaptations, will be a determining factor in the survival of tree species to a rapidly changing climate. In this dynamic environment, young trees represent a critical demographic stage forcing species range shifts, especially toward the leading edge of the species distribution (Talluto et al., 2017; Parks et al., 2019; Shriver et al., 2022). Young trees play an outsized role in post-fire restoration, where naturally occurring seed dispersal is either limited or absent (Shive et al., 2018; North et al., 2019; Stevens-Rumann and Morgan, 2019).

This study expands the scope of forest recovery beyond initial seedling establishment following disturbance, to include persistence against subsequent climate pressures during early tree development. Our results uncover substantial species-specific responses to drought and heat at the juvenile growth stage. When normalizing species to different levels of stress, Pinus flexilis survived three times longer than Pinus ponderosa at the juvenile growth stage, while Pinus edulis, Picea engelmannii, and Pseudotsuga menziesii had comparable tolerance with intermediate survival rates. We also saw a pronounced temperature effect, where warmer baseline temperatures allow drought conditions to be realized more quickly and cooler baseline temperatures delayed the onset of drought conditions. Thus, higher-elevation species at cooler temperatures may be able to withstand longer droughts. This could be a mechanism to draw species up in elevation over time in response to climate change (Monleon and Lintz, 2015).Importantly, our results complement observations of recruitment failure in post-fire landscapes that are projected to become more pronounced under future climate (Davis et al., 2023) by experimentally isolating the effects of drought with and without a heatwave.

The relatively high vulnerability of Pinus ponderosa to drought seems surprising given its dominance in dry landscapes across the Southwest (Oliver and Ryker, 1990). Pinus ponderosa is highly adapted to dry, fire-frequent landscapes, with traits such as thick bark and elevated open crowns to avoid fire damage, although these traits are expressed primarily in the adult stage (Pausas et al., 2017; Keeley and Pausas, 2022). As a mature tree, Pinus ponderosa can survive low intensity fire disturbance, but juveniles are considerably less fire resistant (Partelli-Feltrin et al., 2020). Juvenile regeneration in the species is episodic and highly dependent on climate (White, 1985; Mast et al., 1999); several years of favorable climate are needed to support juvenile survival and growth (Kemp et al., 2019). Although mature Pinus ponderosa may be favored at lower elevations due to its fire adaptations, juveniles facing climate extremes may be vulnerable in the warm-dry tails of their climatic range (Kolb et al., 2007; Remke et al., 2020). In contrast, Pinus flexilis has a relatively small abundance in the Jemez Mountains despite its drought tolerance (Muldavin and Tonne, 2003). Pinus flexilis is concentrated on rocky dry sites in cool environments where they hold a competitive advantage over other conifers (Menon et al., 2021; Hankin et al., 2023). In the southern Rocky Mountains, they are concentrated in the mixed conifer and spruce-fir zones, and are successionally replaced by spruce and fir species on most sites (Rebertus et al., 1991). In an experimental study of seedling emergence under low soil moisture conditions, Pinus flexilis had the highest emergence rates of the three high-elevation five-needle pines in the study (Hankin et al., 2023). As this study suggests, Pinus flexilis may gain a competitive advantage in sites where forest species must contend with lengthening droughts.

Natural and managed forest regeneration is dependent on the ability of tree juveniles to survive expected climate anomalies. Many forest and woodland ecosystems in the Southwest are forecast to be highly vulnerable to vegetation conversion, based on species distributions and shifting climate envelope models (Triepke et al., 2019). Strategic tree planting, especially in the context of prolonged drought, is an active management technique that can mitigate regeneration failure (Field et al., 2020; Pinto et al., 2023). Planting in refugia, where climate conditions are most suitable (McDowell et al., 2019; Krawchuk et al., 2020) can facilitate regeneration among disturbed landscapes (Stevens et al., 2021). The juvenile age class, as examined in this study, is especially relevant to managed reforestation strategies (Landis et al., 2010).

Our study employed juvenile trees acquired from nurseries, and our results can thus inform managers concerned with potential success of outplanting under emerging climate conditions. This demographic stage is key to implementing mechanistic models of forest ecological resilience (Falk, 2017; Falk et al., 2022). In turn, ecological resilience is the foundation for emerging management and adaptation paradigms, including the Resist-Resilience-Transition (Peterson St-Laurent et al., 2021) or Resist-Accept-Direct (Thompson et al., 2021) management frameworks. Resisting change through recovery would be focused on planting a species according to its historic distribution, which may require ongoing active management to resist persistent drought conditions. On the other hand, planting species at higher elevations with cooler temperatures, beyond their typical climatic zones, would represent a transition approach, facilitating upward migration of species (Brusca et al., 2013; Parks et al., 2019). Other transition approaches might include planting more drought-tolerant species which are not currently dominant on the landscape but may become more widespread in the future. Adaptive management will be especially important in the face of uncertain climate futures (Jackson, 2021).

In conclusion, our incorporation of drought and temperature tolerance at the juvenile demographic stage fills a key information gap regarding potential species persistence, recovery, and range shifts over time. Future management actions and planting strategies can consider these experimental outcomes to anticipate landscape-level climate impacts on juvenile tree survival as the key to species adaption in a rapidly changing world.