DOI: 10.5061/dryad.hhmgqnkj9
title: Tradeoffs between leaf cooling and hydraulic safety in a dominant arid land riparian tree species

description: 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 (&gt;40 °C maximum summer
 temperature) and cool-adapted (&lt;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. 
description: <strong>2.1
 Study site</strong> An experimental
 common garden was established in October 2014 with 16 <em>Populus
 fremontii </em>populations (~4,100 propagated cuttings) that
 collectively represent the climatic and elevational range of the species
 (Cooper <em>et al</em>. 2019; Hultine <em>et
 al</em>. 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 (&gt;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
 <em>et al</em>. 2020). Results from Blasini <em>et
 al</em>. (2020) suggest trait expression in <em>P.
 fremontii</em> 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 <em>et
 al</em>. 2011). <strong>2.2 Meteorological data
 </strong> 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. <strong>2.3 Morphological
 traits  </strong>
 <strong>2.3.1 Stomatal anatomy</strong>
 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<sub>stom</sub>) was
 calculated as the number of stomata in eight (0.59
 mm<sup>2</sup>) digital images at 10× magnification per
 genotype. Stomatal size (S<sub>stom</sub>) (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<sup>−2</sup> s<sup>−1</sup>) following Franks and Farquhar (2001).  where d<sub>w</sub> is the diffusivity of water in air (2.43 × 10<sup>−5</sup> m<sup>2/s</sup>), v is the molar volume of air (0.024 m<sup>3</sup>/mol; Jones 2014), D<sub>stom</sub> is the stomatal density, amax is the maximum area of the open stomatal pore, approximated as π(p/2)<sup>2</sup>, where p is stomatal pore length (µm), assumed to be stomatal length divided by two (Franks & Farquhar, 2007) <strong>2.3.2 Leaf traits</strong> Specific leaf area (SLA, cm<sup>2</sup> g<sup>-1</sup>) 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<sub><em>l</em></sub>, cm<sup>2</sup>) and leaf width (w<sub><em>l</em></sub>, cm) were derived from the average leaf size from these measurements (Ackerly <em>et al</em>. 2002). <strong>2.3.3 Whole tree allometry</strong> 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<sub><em>l</em></sub>) and sapwood area (A<sub>s</sub>). 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<sub>c</sub>) 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 <em>et al</em>. 2012). <strong>2.4 Sap-flux-scaled canopy transpiration and stomatal conductance </strong> We installed heat dissipation sensors (Granier 1987) that measured stem sap flux density (J<sub>s</sub>, g H2O m<sup>-2</sup> sapwood s<sup>-1</sup>) on all 56 genotypes from June 2<sup>nd</sup> (day 153) to October 2<sup>nd</sup> (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<sub>s</sub> (g cm<sup>-2</sup> s<sup>-1</sup>), was calculated as   For diffuse porous tree species (e.g., <em>Populus fremontii</em>), 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 <em>et al</em>. 2010): Where ΔT is the difference in temperature between the heated and unheated probes and ΔT<sub>0</sub> 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<sub>0</sub> was used to calculate k over several days. We calculated canopy transpiration (E, g m<sup>- 2</sup> s<sup>- 1</sup>) using the total sap flux density and sapwood area to leaf area ratio (A<sub>s</sub>:A<sub><em>l</em></sub>) according to:   From the sap flux measurements, we also calculated canopy conductance (Gc, mmol m<sup>- 2</sup> s<sup>- 1</sup>) using a simplified version of the Penman–Monteith equation (Campbell and Norman 1998; Hultine <em>et al</em>. 2013; Monteith and Unsworth 2013): where A<sub>s</sub> is conducting sapwood area (m<sup>2</sup>), Al is the total leaf area (m<sup>2</sup>),  is the psychrometric constant (kPa K<sup>-1</sup>), λ is the latent heat of vaporization (J kg<sup>-1</sup>), ρ is the density of moist air (kg m<sup>-3</sup>), and c<sub><em>p</em></sub> is the specific heat of air at constant pressure (J kg<sup>-1</sup> K<sup>-1)</sup>. Internal bark diameter and the depth of hydro-active xylem was estimated to obtain A<sub>s</sub>. Because of the young age of the trees (2.5 years), and because <em>Populus </em>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 <em>et al</em>. 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<sub><em>l</em></sub>) and canopy area (A<sub>c</sub>) were used to calculate intraspecific differences in LAI and therefore canopy boundary layer resistances, where LAI is given by:   Canopy stomatal conductance (G<sub>s</sub>) was extracted from measurements of G<sub>c</sub> by evaluating leaf boundary layer conductance (G<sub>bl</sub>; mmol m<sup>-2</sup>s<sup>-1</sup>), which can be small enough in broadleaf plants to decouple plant canopies from atmospheric conditions. We therefore calculated G<sub>bl</sub> according to Jones, 2014) to compare with calculated values of G<sub>s</sub> (shown below):   Where (d<sub><em>l</em></sub>) is the mean leaf characteristic dimension calculated for each genotype (d<sub><em>l</em></sub> = 0.72 · leaf width) and (m<sub>c</sub>) is the mean canopy wind speed (Campbell & Norman 1998; Jones 2014). Mean mc (m s<sup>-1</sup>) 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<sub>m</sub> is the roughness length (m), and 0.4 is the von Karman constant. We used population specific values of G<sub>c</sub> and G<sub>bl</sub> to estimate canopy stomatal conductance:   We calculated a dimensionless decoupling coefficient (Ω) (Martin 1989; Hultine <em>et al</em>. 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<sub>r</sub> 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.   <strong>2.5 Leaf water potentials </strong> From June to September 2017, leaf water potentials (Ψ) were measured every month at predawn (Ψ<sub>pd</sub>; 0300 – 0500 h local time) and midday (Ψ<sub>md</sub>; 1100 – 1300 h) on each genotype that was instrumented with sap flux probes using the Scholander pressure chamber (PMS Instruments, Corvallis, OR; Scholander <em>et al</em>., 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 Ψ<sub>pd</sub> and Ψ<sub>md</sub> (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.      <strong>2.6 Leaf temperature</strong> 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<sup>th</sup> (day 240) and September 1<sup>st</sup> (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). <strong>2.7 Statistical analysis and leaf energy balance modeling </strong> 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. <strong>2.7.1. Analysis of trait variation between cool- and warm-adapted ecotypes</strong> 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<sub>2</sub> 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 (&lt;1 kPa). This sensitivity of stomatal conductance to changes in vpd can be estimated from (Domec <em>et al</em>. 2009; Hultine <em>et al</em>. 2013; Oren <em>et al</em>. 1999): where G<sub>sref</sub> is the value of G<sub>s</sub> at vpd = 1 kPa in a log-linear relationship and m (the slope of the regression fit) describes the sensitivity of G<sub>s</sub> to changes in vpd (i.e., ln vpd). We also calculated stomatal sensitivity standardized by G<sub>sref</sub>- (S) according to Oren et al (1999) as -m G<sub>sref</sub><sup>-1</sup>. Regression analyses was used to investigate the relationships between G<sub>s</sub> and G<sub>s</sub>:G<sub>Sref</sub> 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<sub>s</sub> and G<sub>s</sub>:G<sub>Sref</sub> between ecotypes in response to log(vpd) and the interaction ecotype*log(vpd) were analyzed using analysis of covariance (ANCOVA). Differences in Ψ<sub>pd</sub>, Ψ<sub>md</sub> 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. <strong>2.7.2 Leaf energy budget model parameterization            </strong> 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<sup>nd</sup> (day 153) to October 2<sup>nd</sup> (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<sup>nd </sup>(day 153) to October 2<sup>nd </sup>(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 <em>et al</em>. 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 <em>P. fremontii</em>.    <strong>2.7.3. Analysis of trait variation among populations</strong> We performed a principal component analysis (PCA) to analyze the relationship between seven morphophysiological traits (S<sub><em>l</em></sub>, A<sub>s</sub>:A<sub><em>l</em></sub>, D<sub>stom</sub>, G<sub>s</sub>, Ψ<sub>md</sub>, SLA and S<sub>stom</sub>) and ΔT at the population level using the ‘factoextra’ and ‘FactoMineR’ packages (Kassambara & Mundt, 2017; Lê <em>et al</em>. 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.