Leaf carbon gain optimization in hot environments requires balancing leaf thermoregulation with avoiding excessive water loss via transpiration and hydraulic failure. The tradeoffs between leaf thermoregulation and transpirational water loss can determine the ecological consequences of heat waves that are increasing in frequency and intensity. We evaluated leaf thermoregulation strategies in warm (>40 °C maximum summer temperature) and cool-adapted (<40 °C maximum summer temperature) genotypes of the foundation tree species, Populus fremontii using a common garden near the mid-elevational point of its distribution. We measured leaf temperatures and assessed three modes of leaf thermoregulation: leaf morphology, midday canopy stomatal conductance, and stomatal sensitivity to vapor pressure deficit. Data were used to parameterize a leaf energy balance model to estimate contrasts in midday leaf temperature in warm- and cool-adapted genotypes. Warm-adapted genotypes had 39% smaller leaves and 38% higher midday stomatal conductance, reflecting a 3.8 °C cooler mean leaf temperature than cool adapted genotypes. Leaf temperatures modeled over the warmest months were on average 1.1 °C cooler in warm- relative to cool-adapted genotypes. Results show that plants adapted to warm environments are predisposed to tightly regulate leaf temperatures during heat waves, potentially at an increased risk of hydraulic failure. 

2.1 Study site

An experimental common garden was established in October 2014 with 16 Populus fremontii populations (~4,100 propagated cuttings) that collectively represent the climatic and elevational range of the species (Cooper et al. 2019; Hultine et al. 2020a). The garden was located within the Agua Fria National Monument in central Arizona (34° 15’ 34.42” N; 112° 03’ 29.39” W; elevation 988 m) (Figure 1) and was established on a 1.2 Ha portion of former cropland next to the intermittently flowing Agua Fria River. During the winter of 2013 – 2014, cuttings were collected from a total of 12 genotypes per population. Genotypes were collected at least 20 meters apart to avoid using clones within each population. The individual cuttings were treated with root hormone and potted in the Northern Arizona University greenhouse for 4 months. In the garden, 0.3 m tall saplings were planted 2 m apart from each other in a randomized block design with a total of four replicated blocks, each with 16 populations comprising 64 genotypes each. A drip irrigation system was used to water each tree with approximately 20 L, two to three times per week during the growing season.

From the original 16 populations established in the garden, we studied 7 populations with a total of 56 genotypes (n = 8 genotypes per population) representing the broadest possible range in mean annual temperature of the source populations, from 10.7 to 22.6 °C, and an elevation gradient from 72 to 1,940 m (Figure 1). In addition to the local Agua Fria National Monument population, three populations, respectively from sites with higher and lower mean maximum summer temperatures than the common garden location were selected. The three populations from the lower Sonoran Desert were defined as a warm-adapted ecotype because the extreme mean maximum summer temperatures (>40 °C) they experience at their source sites is above that of the common garden location (Figure 2). The three populations from higher elevation provenances in the Sonoran Desert and Colorado Plateau were sourced from locations with similar or lower mean maximum summer temperatures than the common garden location and therefore were categorized as a cool adapted ecotype (Figure 2).

The experimental common garden and some of the genotypes used in this research were part of a previous investigation that studied morpho-physiological trait variability at multiple trait spectra in relation to local temperatures at the population source sites (Blasini et al. 2020). Results from Blasini et al. (2020) suggest trait expression in P. fremontii is highly coordinated and reflect local adaptation to either exposure to freeze-thaw conditions in high elevation populations or exposure to extreme thermal stress in low-elevation populations. In the present investigation, we evaluated intraspecific differences in leaf thermoregulation in relation with mean annual temperature (MAT) and mean maximum summer temperatures (MMST) transfer distances (Table 2), defined as the MAT and MMST of the source population location subtracted from the MAT and MMST of the common garden location (Grady et al. 2011).

2.2 Meteorological data

A micrometeorological station installed at the garden measured relative humidity, air temperature, photosynthetically active radiation (Q, mmol m-2 s-1), and wind speed continuously every 30 seconds and stored as 30 minutes means from May 30 (day 151) to Oct 23 (day 296) of 2017 with a Campbell CR10X-2M data logger (Campbell Scientific, Logan, UT, USA). Air temperature and relative humidity were measured with a shielded Vaisala HMP 60 AC temperature/humidity probe (Vaisala, Woburn, MA, USA) placed 3 m above the ground surface. Photosynthetic active radiation was measured with an Apogee SQ-110-SS sun calibration quantum Sensor (Apogee Instruments, Logan, UT, USA). We used air temperature and relative humidity to calculate air vapor pressure deficit (vpd, kPa) using both half-hourly and daily averages.

2.3 Morphological traits 

2.3.1 Stomatal anatomy

In 2016, we randomly collected fully expanded leaves from mid-height and south-facing canopy of each genotype (n = 56) to determine stomatal density, length, width and area. We followed the nail polish impression method (Hilu & Randall 1984) to obtain four impressions per leaf, two impressions in both the abaxial and adaxial sides of the leaves (n = 560 impressions). An Olympus CX41 light microscope (Olympus Corp., Center Valley, PA, USA) was used to observe and obtain two images from each impression with a Moticam Pro 282A camera (Motic Instruments, Richmond, BC, Canada). Stomatal density (D_stom) was calculated as the number of stomata in eight (0.59 mm²) digital images at 10× magnification per genotype. Stomatal size (S_stom) (length × width) was observed on 700 stomata from digital images at 40× magnification (n = 100 per population) using an open-source imaging program, ImageJ (https://imagej.nih.gov/ij/). We calculated maximum theoretical stomatal conductance to water vapor (Gsmax, mmol m⁻² s⁻¹) following Franks and Farquhar (2001).

 where d_w is the diffusivity of water in air (2.43 × 10⁻⁵ m^2/s), v is the molar volume of air (0.024 m³/mol; Jones 2014), D_stom is the stomatal density, amax is the maximum area of the open stomatal pore, approximated as π(p/2)², where p is stomatal pore length (µm), assumed to be stomatal length divided by two (Franks & Farquhar, 2007)

2.3.2 Leaf traits

Specific leaf area (SLA, cm² g^-1) was calculated as the one-sided area of a fresh leaf, divided by its oven-dry mass (Wright & Westoby 2002). SLA was measured in June, July, and September 2017. A subset of 12 to 20 fully expanded leaves from the mid and south facing section of the canopy were collected per genotype and scanned with a high-resolution computer scanner, and one-sided leaf area was measured with Image J. The scanned leaves were then oven-dried for 72 hours at 75 °C and weighed to calculate SLA. Leaf size (S_l, cm²) and leaf width (w_l, cm) were derived from the average leaf size from these measurements (Ackerly et al. 2002).

2.3.3 Whole tree allometry

In July 2017, we used allometric relationships between whole-tree stem diameter and leaf area through a branch summation approach to estimate whole-tree canopy leaf area (A_l) and sapwood area (A_s). The diameter of all leaf-bearing branches from the main stem in each of the 56 genotypes were measured with a digital caliper. To calculate whole canopy leaf area, a subset of the collected leaves per genotype was scanned with a high-resolution computer scanner, and one-sided leaf area was measured with Image J. Then, we generated a regression of branch diameter to leaf area from a subset of branches per genotype. Scanned leaves were oven-dried for 72 hours at 75 °C within each subset of branches, and then their weight was multiplied by SLA to determine total leaf area of the branch. (see section 2.3.2). Whole-tree height (H), canopy diameters (4–8 measurements per genotype) and their respective canopy areas were measured five times during the 2017 growing season with a telescoping measuring pole. Canopy area (A_c) was determined using the ellipse equation, πab, where a is the mean radius of longest canopy axis and b is the radius of two perpendicular canopy axes (Ansley et al. 2012).

2.4 Sap-flux-scaled canopy transpiration and stomatal conductance

We installed heat dissipation sensors (Granier 1987) that measured stem sap flux density (J_s, g H2O m^-2 sapwood s^-1) on all 56 genotypes from June 2^nd (day 153) to October 2^nd (day 275) 2017. Each sensor consisted of a pair of 20 mm long, 2 mm diameter stainless steel probes inserted approximately 15 cm apart along the axis of the hydro-active xylem. The sap flux density was calculated from the differences in temperature between the heated and unheated reference probes. Sap flux density, J_s (g cm^-2 s^-1), was calculated as

 

For diffuse porous tree species (e.g., Populus fremontii), q is the prefactor coefficient (0.0119), p is the scaling exponent (1.23), and k is related to the temperature difference between the two probes (Granier 1987; Bush et al. 2010):

Where ΔT is the difference in temperature between the heated and unheated probes and ΔT₀ is the temperature difference during hydrostatic conditions (data provided in repository). We assumed that hydrostatic conditions only occurred during evening periods when vpd was at or near zero. Thus, in some cases a single value for ΔT₀ was used to calculate k over several days.

We calculated canopy transpiration (E, g m^{- 2} s^{- 1}) using the total sap flux density and sapwood area to leaf area ratio (A_s:A_l) according to:

 

From the sap flux measurements, we also calculated canopy conductance (Gc, mmol m^{- 2} s^{- 1}) using a simplified version of the Penman–Monteith equation (Campbell and Norman 1998; Hultine et al. 2013; Monteith and Unsworth 2013):

where A_s is conducting sapwood area (m²), Al is the total leaf area (m²),  is the psychrometric constant (kPa K^-1), λ is the latent heat of vaporization (J kg^-1), ρ is the density of moist air (kg m^-3), and c_p is the specific heat of air at constant pressure (J kg^-1 K^-1).

Internal bark diameter and the depth of hydro-active xylem was estimated to obtain A_s. Because of the young age of the trees (2.5 years), and because Populus trees tend to have large active sapwood depths with uniform sap velocities (Lambs & Muller 2002), we assumed the active sapwood included the entire cross-sectional area beneath the bark. Measurements of leaf area index (LAI), the projected leaf area per unit of ground area (Watson 1947; Bréda 2003; Chapin et al. 2011), provided a way to estimate the physical boundaries between the whole-tree canopy and the surrounding atmosphere (Bréda 2003). Therefore, whole tree leaf area (A_l) and canopy area (A_c) were used to calculate intraspecific differences in LAI and therefore canopy boundary layer resistances, where LAI is given by:

 

Canopy stomatal conductance (G_s) was extracted from measurements of G_c by evaluating leaf boundary layer conductance (G_bl; mmol m^-2s^-1), which can be small enough in broadleaf plants to decouple plant canopies from atmospheric conditions. We therefore calculated G_bl according to Jones, 2014) to compare with calculated values of G_s (shown below):

 

Where (d_l) is the mean leaf characteristic dimension calculated for each genotype (d_l = 0.72 · leaf width) and (m_c) is the mean canopy wind speed (Campbell & Norman 1998; Jones 2014). Mean mc (m s^-1) was estimated from wind speed (m) measured at 3 m above the ground level and by multiplying canopy frictional velocity (mv) to m. mv was calculated following Campbell & Norman (1998) as:

 

Where z is the genotype canopy height (m), d is the zero-plane displacement (m), z_m is the roughness length (m), and 0.4 is the von Karman constant. We used population specific values of G_c and G_bl to estimate canopy stomatal conductance:

 

We calculated a dimensionless decoupling coefficient (Ω) (Martin 1989; Hultine et al. 2013) to evaluate the sensitivity of transpiration to changes in boundary layer conductance.

 

Where ε is the change of latent heat to the change in sensible heat of saturated air and G_r is the long-wave radiative transfer conductance. It  is expected to reach its upper limit (1.0) as the influence of stomatal resistance over transpiration decline.

 

2.5 Leaf water potentials

From June to September 2017, leaf water potentials (Ψ) were measured every month at predawn (Ψ_pd; 0300 – 0500 h local time) and midday (Ψ_md; 1100 – 1300 h) on each genotype that was instrumented with sap flux probes using the Scholander pressure chamber (PMS Instruments, Corvallis, OR; Scholander et al., 1965; Turner, 1988). To take these measurements, a single shoot tip from each of the 56 genotypes was cut with a sharp razor blade at mid-height and south-facing canopy. Differences between Ψ_pd and Ψ_md (DY) were calculated for each genotype, population, and ecotype over each measurement period, to provide an index of the transpiration-induced changes in water potential gradients from the roots to the leaves.     

2.6 Leaf temperature

We measured leaf temperature on 17 warm- and 24 cool-adapted genotypes (total n = 41) instrumented with the sap flux probes in the common garden between 13:00 and 15:00 h of August 28^th (day 240) and September 1^st (day 244) of 2017: two of the warmest days during the study. We evaluated leaf temperature on three to four separate leaves in each individual genotype using a thermal imaging camera ThermaCam (Flir One, Flir Systems, Wilsonville, OR USA). This handheld device integrates a thermal and visual sensor of 80x60 and 1440x1080 pixels, respectively, to a smartphone, with a typical accuracy range of ± 3 °C or ± 5% (https://www.flir.com/products/flir-one-gen-3/). Leaf temperatures were taken 30 cm away from a full expanded leaf at three to four different locations in the canopy. Only leaves located in the middle (1.5 to 2.0 m) and sun-lit areas of each tree canopies were used to measure leaf temperature. This resulted in a total of 81 and 92 leaf temperature measurements for day 240 and day 244 respectively (total n = 173 measurements, 99 measurements for cool adapted genotype and 74 for warm adapted genotypes on both days). We used air temperature data collected from the on-site micrometeorological station to calculate the difference between air and leaf temperatures (ΔT).

2.7 Statistical analysis and leaf energy balance modeling

All statistical analyses were conducted in R version 4.0.2 (R Development Core Team 2011). Prior to analyzing the data, we examined whether each variable met the assumptions of normality and homogeneity of variance, using a Shapiro and Barlett test.

2.7.1. Analysis of trait variation between cool- and warm-adapted ecotypes

Morphological trait comparisons between cool- and warm-adapted ecotypes, including all measurements of leaf morphology, leaf area to sapwood area ratios, and mean daily sap-flux-scaled canopy fluxes were conducted using a standard student’s t-test.

Because riparian tree species are found exclusively in places with abundant water available, stomatal conductance and whole-tree water use in this species are intrinsically influenced by atmospheric characteristics like irradiance, atmospheric CO₂ concentrations and atmospheric vapor pressure deficit (Landsberg, Waring & Ryan 2017). Specifically, increases in vpd have been found to correlate with decreases in stomatal conductance while the stomatal sensitivity to changes in vpd has been described to be proportional to the stomatal conductance at low vpd levels (<1 kPa). This sensitivity of stomatal conductance to changes in vpd can be estimated from (Domec et al. 2009; Hultine et al. 2013; Oren et al. 1999):

where G_sref is the value of G_s at vpd = 1 kPa in a log-linear relationship and m (the slope of the regression fit) describes the sensitivity of G_s to changes in vpd (i.e., ln vpd). We also calculated stomatal sensitivity standardized by G_sref- (S) according to Oren et al (1999) as -m G_sref^-1. Regression analyses was used to investigate the relationships between G_s and G_s:G_Sref of the cool- and warm-adapted ecotypes to log(vpd) during the hottest time of the day (11:00 to 19:00). Comparisons in mean G_s and G_s:G_Sref between ecotypes in response to log(vpd) and the interaction ecotype*log(vpd) were analyzed using analysis of covariance (ANCOVA).

Differences in Ψ_pd, Ψ_md and ΔΨ between cool- and warm-adapted genotypes were analyzed by individual mixed-effects repeated measures ANOVA (type III with Satterthwaite’s method) using the ‘lmer’ R package (Bates, Mächler, Bolker, & Walker 2015; Kuznetsova, Brockhoff & Christensen 2017). In this test, individual traits were represented as response variables while the group (cool and warm-adapted ecotypes) and month were treated as categorical fixed effects with two and three levels, respectively. Individual genotype nested within ecotypes was incorporated as a random effect.

2.7.2 Leaf energy budget model parameterization           

A leaf energy balance model was executed in R (R Core Team, 2018) through the ‘tealeaves’ package (Muir 2019). The model calculates leaf temperature from a suite of leaf traits, environmental parameters, and physical constants. Leaf traits included leaf size, stomatal ratio (stomata density adaxial: stomatal density abaxial), and mean canopy stomatal conductance during the hottest time of the day (11:00 to 19:00) from June 2^nd (day 153) to October 2^nd (day 275) 2017 were included in the model. The environmental parameters used in the model included air temperature, relative humidity, and wind speed collected at the common garden during the hottest time of the day from June 2^nd (day 153) to October 2^nd (day 275) 2017. Other environmental parameters included in the model were atmospheric pressure at 998 meter above sea level (90.0 kPa), reflectance for shortwave irradiance (albedo) (0.2, unitless), and incident short-wave (solar) radiation flux density (1000 W/m2) (Muir 2019, Okajima et al. 2012).

A sensitivity analysis was performed on three and two environmental and morphophysiological variables, respectively, to determine their overall effect on leaf temperature resulted from the leaf energy balance model. These variables were air temperature, relative humidity, wind speed (environmental) and stomatal conductance and leaf size (morphophysiological). We used the ‘konfound’ package in R to run the sensitivity analysis (Frank, Maroulis, Duong & Kelsey 2013) to determine the influence that environmental variables (relative humidity, wind speed and air temperature) and morphophysiological traits (leaf size and stomatal conductance) have on the modeled leaf temperatures in P. fremontii.   

2.7.3. Analysis of trait variation among populations

We performed a principal component analysis (PCA) to analyze the relationship between seven morphophysiological traits (S_l, A_s:A_l, D_stom, G_s, Ψ_md, SLA and S_stom) and ΔT at the population level using the ‘factoextra’ and ‘FactoMineR’ packages (Kassambara & Mundt, 2017; Lê et al. 2008). We determined the number of meaningful PCA axes using the Kaiser criterion and the Broken Stick Model in the ‘vegan’ and ‘biodiversity’ R package. Trait representation in the principal component biplot was based on the magnitude of the correlation between each trait and the principal component. Thus, traits in this biplot were represented as vectors with a length and direction indicating the strength and trend of a given trait’s relationship among other traits. Specific location of the vector in the biplot indicates the positive or negative impact that a trait has on each of the two components x-axis, first component (PC1) and y-axis, second component (PC2). To analyze the relationship between the seven populations and the traits distribution in the PCA biplot, we constructed seven 95% confidence ellipses based on the PCA scores of each population. Subsequently, we performed ANOVA Tukey’s HSD tests to assess significant differences in PC axes scores at the population level.