Radiation and Drought Impact Residual Leaf Conductance in Two Oak Species With Implications for Water Use Models

Stomatal closure is one of the earliest responses to water stress but residual water losses may continue through the cuticle and incomplete stomatal closure. Residual conductance (gres) plays a large role in determining time to mortality but we currently do not understand how do drought and shade interact to alter gres because the underlying drivers are largely unknown. Furthermore, gres may play an important role in models of water use, but the exact form in which gres should be incorporated into modeling schemes is currently being discussed. Here we report the results of a study where two different oak species were experimentally subjected to highly contrasting levels of drought (resulting in 0, 50 and 80% losses of hydraulic conductivity) and radiation (photosynthetic photon flux density at 1,500 μmol m–2 s–1 or 35–45 μmol m–2 s–1). We observed that the effects of radiation and drought were interactive and species-specific and gres correlated positively with concentrations of leaf non-structural carbohydrates and negatively with leaf nitrogen. We observed that different forms of measuring gres, based on either nocturnal conductance under high atmospheric water demand or on the water mass loss of detached leaves, exerted only a small influence on a model of stomatal conductance and also on a coupled leaf gas exchange model. Our results indicate that, while understanding the drivers of gres and the effects of different stressors may be important to better understand mortality, small differences in gres across treatments and measurements exert only a minor impact on stomatal models in two closely related species.

Stomatal closure is one of the earliest responses to water stress but residual water losses may continue through the cuticle and incomplete stomatal closure. Residual conductance (g res ) plays a large role in determining time to mortality but we currently do not understand how do drought and shade interact to alter g res because the underlying drivers are largely unknown. Furthermore, g res may play an important role in models of water use, but the exact form in which g res should be incorporated into modeling schemes is currently being discussed. Here we report the results of a study where two different oak species were experimentally subjected to highly contrasting levels of drought (resulting in 0, 50 and 80% losses of hydraulic conductivity) and radiation (photosynthetic photon flux density at 1,500 µmol m −2 s −1 or 35-45 µmol m −2 s −1 ). We observed that the effects of radiation and drought were interactive and species-specific and g res correlated positively with concentrations of leaf non-structural carbohydrates and negatively with leaf nitrogen. We observed that different forms of measuring g res , based on either nocturnal conductance under high atmospheric water demand or on the water mass loss of detached leaves, exerted only a small influence on a model of stomatal conductance and also on a coupled leaf gas exchange model. Our results indicate that, while understanding the drivers of g res and the effects of different stressors may be important to better understand mortality, small differences in g res across treatments and measurements exert only a minor impact on stomatal models in two closely related species.

INTRODUCTION
Plant transpiration through stomatal pores and leaf cuticles dominates global evapotranspiration (Hetherington and Woodward, 2003). As water stress intensifies under global warming, there is an increasing interest toward understanding ecological variation in residual leaf conductance (g res ).
After stomatal closure, water loss continues until mortality due to a mixture of cuticular water loss and incomplete stomatal closure (residual conductance; Blackman et al., 2016;Martin-Stpaul et al., 2017).
Studies addressing ecological and physiological variation in the drivers of residual conductance are currently rare (Heredia-Guerrero et al., 2018). According to a recent review on this topic (Duursma et al., 2019), only 10 studies have addressed the effect of drought on g res and, from those, only 4 had been performed on trees. Consequently, multifactorial studies addressing ecological variation in residual conductance are much needed to understand its variation. For instance, while shade and drought are both known to decrease residual conductance (Boyer et al., 1997;Shepherd and Wynne Griffiths, 2006), it is currently unknown whether the effect of both stressors would be additive or interactive. However, the effects of residual conductance on mortality have been documented to be dramatic: time to mortality nearly doubles if g res declines from 4 to 2 mmol m −2 s −1 (Duursma et al., 2019).
Understanding the physiological and ecological drivers of g res has been the topic of some discussion (Riederer and Müller, 2007;Fernández et al., 2017). Some studies report that variations in the degree of sclerophylly (as indicated by leaf mass area) would increase g res because leaves that are more scleromorphic will show thicker cuticles, but other work has demonstrated that changes in wax composition may compensate for such effect (Bueno et al., 2020). Another alternative, explored to a lesser degree, is that further reductions in g res may be inhibited by changing carbohydrate allocation priorities (Zhang et al., 2020). In other words, as non-structural carbohydrate reserve pools deplete, cuticle production to prevent cuticular water losses may be limited by NSC availability.
Understanding variation in residual conductance is also necessary for models of water use (Leuning, 1995;Barnard and Bauerle, 2013;De Kauwe et al., 2015), where residual conductance acts as the intercept of commonly used stomatal models (g int ). The most common stomatal models being used in land surface models are Ball-Berry model types, which have the general form: Where g s is stomatal conductance, A, C a , and D represent photosynthesis, ambient CO 2 concentration and vapor pressure deficit, respectively, and m is the slope parameter. When g int is estimated through regression fitting, it may either be equal to 0, which creates problems because then the ratio of intercellular to ambient CO 2 (C i /C a ) does not vary with light (Collatz et al., 1991;Leuning, 1995;Duursma et al., 2019), or it may be negative, which is nonsensical.
There are at least two possible definitions of g int : (1) g 0 , which represents the lowest conductance reached as photosynthesis tends to 0 because light declines (Leuning, 1995;Barnard and Bauerle, 2013); (2) g min , which refers to the residual conductance after (complete or not) stomatal closure under strong water stress (Duursma et al., 2019). We note that some studies use g min and g res interchangeably but, for clarity, we will differentiate them here as previously defined.
The problem then becomes how to measure g 0 and g min . g 0 could simply be measured as daytime conductance (g d ) under low light in non-droughted plants and, similarly, g min could similarly be measured from g d in droughted plants (for as long as photosynthesis tends to zero, in both cases; Barnard and Bauerle, 2013;Duarte et al., 2016). Additionally, residual conductance has most often been measured by monitoring the water mass loss in detached leaves (g MLD ; Kerstiens, 1996;Schuster et al., 2017). g 0 and g min could thus be measured with this method by comparing g MLD in plants that have grown under strong light limitation or under strong water limitation, respectively. The problem with this approach, however, is that some acclimation responses (particularly in response to low radiation) could alter leaf morphology and it is unclear whether g 0 measured through g MLD after low light acclimation would be representative of that in plants without acclimation to low radiation.
An alternative would be to use nocturnal conductance (g n ; Lombardozzi et al., 2017) in non-droughted and droughted plants. An advantage would be that photosynthesis would always be zero in this case. Duursma et al. (2019), however, proposed that g n should not be used given the evidence of active regulation of stomatal conductance overnight (Resco De Dios et al., 2019), and that the drivers of nocturnal conductance could differ from those driving daytime conductance (Ogle et al., 2012). Amongst other processes, g n varies through time due to circadian regulation (Resco De Dios et al., 2015). However, g n often retains some sensitivity to D such that maximum stomatal closure and, potentially, residual conductance, may be achieved at lower D than during the daytime (Barbour and Buckley, 2007). One could thus hypothesize that measurements of g n under high D may be indicative of g res .
Regardless of how g 0 and g min are estimated, Duursma et al. (2019) proposed to replace Eq. 1 by: (2) That is, according to Eq. 2, residual conductance would not be added to the right-hand term of Eq. 1. Instead, one would use measured residual conductance (the maximum between g 0 and g min ) as an actual minimum (De Kauwe et al., 2015). However, this formulation has not yet been tested against data and, therefore, we do not yet know whether it enhances the predictive power of stomatal models.
Here we evaluate the effects of shade and water stress on g res across two different oak species, the deciduous Quercus faginea and the sclerophyll Q. ilex. These two species are common in the calcareous soils from Spain and the Western Mediterranean Basin and we expected conductance to be significantly lower in Q. ilex, a species with a more conservative water use. More specifically, we sought to test: (1) how do drought and shade interact to affect g res ? and (2) what are some of the possible mechanisms underlying variation in g res across drought and shade treatments? Because g MLD is probably the most accepted method to measure residual conductance, here we focused on g MLD . In particular, we addressed whether g MLD would be driven by water stress (as indicated by water potential), NSC, LMA, or nitrogen concentration (N mass , an indicator of photosynthetic capacity), and whether g res could limit respiration. We also sought to understand: (3) whether we obtain different values of g res depending upon whether it is measured from g MLD , g n , and g d ; and (4) how do we incorporate residual conductance into Ball-Berry type stomatal models and what are the consequences of variation in g res across treatments and types of measurements for coupled leaf gas exchange models?

Experimental Design and Growing Conditions
The experiment was performed at the experimental fields from the University of Lleida (Spain; 41.62 N, 0.59 E). We built a rain-out shelter covered by clear polyethylene plastic, which is commonly used in greenhouse building. Half of the structure received solar radiation (sun treatment), with a maximum photosynthetically active radiation (PAR) of 1,500 µmol m −2 s −1 . The other half was covered by a dense shading cloth (shade treatment) with a maximum PAR of 35-45 µmol m −2 s −1 , which was near the light compensation point in this species (data not shown). The structure had openings on both sides to increase ventilation. Temperature inside the rain out shelter was 3 • C higher than outside, but differences between the sun and shade treatment were negligible. Full details on the infrastructure have been provided elsewhere (Resco De Dios et al., 2020).
For this study we sourced 2 year-old seedlings from local nurseries (n = 120). The ecotypes for both species were original from the mountain range of the Iberian System. Plants were grown in 11 L cubic pots (20 cm × 20 cm × 27.5 cm). The substrate used was Humin Substrat Neuhaus N6 [Klasman-Deilmann GmbH, Geeste, Germany], a commercial potting mix. Pots were regularly fertilized with a slow release NPK MgO fertilizer (17-09-11-2, Osmocote Universal, KB, Ecully, France) and daily watered to field capacity until treatment implementation. The position of the pots was randomly shifted every other week.
The plants grew for 4 months into the rainout shelter before experiment inception in July 2017. That is, they developed new leaves under the assigned experimental light conditions. Although we cannot discard legacy effects from the previous growing season in the nursery (Aranda et al., 2001), all plants were treated equally.
We performed a full factorial experiment with the plants experiencing two light treatments crossed with three water stress treatments. Half of the plants grew under the sun treatment and the other half under the shade treatment, as previously described. We implemented three different water stress treatments using three different levels of percent loss of hydraulic conductivity (PLC): (i) P 0 , where plants were irrigated at field capacity; (ii) P 50 , where plants experienced 50% losses in hydraulic conductivity and which represents an important stress; and (iii) P 80 , where the plant experienced 80% losses in hydraulic conductivity, which represents a major stress and potentially mortality (Resco et al., 2009).
We kept plants at field capacity until treatment implementation. We then stopped watering and allowed plants to dehydrate and we measured midday stem water potential ( md ) every other day in a subset of plants (n = 5). The levels of PLC were controlled from the relation between midday shoot water potential ( md ) and PLC values reported previously in vulnerability curves from Quercus faginea Lam. (Esteso-Martínez et al., 2006) and Quercus ilex L. (Peguero-Pina et al., 2014). Shoot md was regularly measured during treatment implementation with a pressure bomb (PMS 1000, PMS Instruments, Albany, Oregon) after clipping the sample and allowing for equilibration in the dark for ∼30 min. Once plants reached the target PLC, we kept soil moisture constant at that level for 2 weeks. This was achieved by weighing a subset of pots (n = 5 per each treatment) and adding back the water that had evaporated every day. We also measured native embolism to test the actual levels of PLC that we achieved in every treatment, as previously published (Resco De Dios et al., 2020). It is important to note that we did not always reach the target PLC levels (see Supplementary  Table 1), but treatment implementation was successful in that we created a gradient in water availability with our treatments. Full details have been provided by Resco De Dios et al. (2020).

Gas Exchange Measurements
Leaf gas exchange was measured with a portable photosynthesis system (LI-6400XT, Li-Cor Inc., Lincoln, NE, United States). We measured 3-5 plants in each treatment at two different periods during the night: between 23:00 h and 01:00 h and between 03:00 h and 05:00 h, and also during the day (10:00-13:00 h). We did not observe significant differences between the stomatal conductance measured over night at the different times (p = 0.79) and measurements were pooled together in subsequent analyses. To understand if measurement errors arising from low flux rates affected our measurements, we also conducted measurements with an empty chamber for 4-5 h, following previously published protocols (Resco De Dios et al., 2013). Results were always one or more orders of magnitude lower or negative. Given these results, we concluded that leaf observations were reliable and that a general correction was not required.
Block temperature was set at 25 • C during the night and at 30 • C during the day, CO 2 at 400 ppm and relative humidity at ∼30%. This meant that D during nighttime measurements was at ∼2.2 kPa, which was substantially higher than that naturally occurring during the night (Resco De Dios et al., 2020). We chose this design to induce nocturnal stomatal closure and test whether g n indicates g res .
During the daytime, we performed measurements at two different levels of PAR: 1,500 µmol m −2 s −1 and 40 µmol m −2 s −1 . We first measured under growth PAR (1,500 µmol m −2 s −1 for plants in the sun treatment and 40 µmol m −2 s −1 for plants in the shadow treatment) and then at the other PAR level (40 µmol m −2 s −1 for plants in the sun treatment and 1,500 µmol m −2 s −1 for plants in the shadow treatment). The leaves were exposed for 10-20 min under the different light intensities until acclimation to the new light level. We only used data measured under growth PAR for analyses, and the rest was reserved for model validation. We note that a sudden exposure to 1,500 µmol m −2 s −1 for a plant growing in the shade would represent a sunfleck, and this could affect the performance of steady-state stomatal models (Way and Pearcy, 2012). When the leaf did not cover the chamber completely, it was scanned and we corrected measurements for leaf area.
In order to parameterize the photosynthesis component of the coupled leaf gas exchange model, we also measured the response of photosynthesis (A) to different internal CO 2 concentrations (C i ) following the protocols from Long and Bernacchi (2003). Briefly, we started measurements with an ambient CO 2 concentration (C a ) of 400 ppm and, after 5 min of acclimation, we sequentially changed C a to 300,250,200,150,100,50,0,400,500,650,800,1000,1250, and 1500 ppm. These measurements were performed at saturating light (1,500 µmol m −2 s −1 ), setting block temperature at 30 • C and with RH as high as we could achieve, which was ∼50%.
Measurements of g MLD g MLD was measured as the mass loss of detached leaves following Phillips et al. (2010) in five leaves per treatment, weighting the leaves every 5 min during 2 h after collection. We wrapped the petiole with paraffilm so that only water lost through the leaf was measured. We performed the measurements in the laboratory, briefly after collection, where we monitored the temperature and relative humidity. Residual conductance was then calculated as: where E MLD is mass loss per projected leaf area (mol m −2 s −1 ), P is atmospheric pressure (kPa) and D is the vapor pressure deficit (kPa). g 0 was defined as g MLD when the leaf originated from the shade treatment at P 0 (P 0 _shade) and g min was defined as g MLD at P 80 in the sun treatment (P 80 _sun).

Analyses of Non-structural Carbohydrates and Elemental Composition
To better understand the physiological mechanisms explaining variations in g res with treatments, we analyzed the concentrations of non-structural carbohydrates, changes in leaf mass per area (LMA) and nitrogen concentrations (N mass ). We collected all the leaves in five plants for each treatment. Immediately after collection, we scanned the leaves to measure the total area and they were then microwaved for 30 s and 700W to stop further metabolic processes. We then oven dried the samples (48 h in 105 • C) and recorded the dry mass. Leaf area and dry weight was used to estimate LMA. We followed previously developed protocols for extracting the percentage for sugar and starch (Palacio et al., 2007). This method consists of grinding the dried leaves with a mill (IKA A10, IKA-Werke, Staufen, Denmark) and making two extractions: one for extracting soluble sugars (sugars from now on) and a second extraction for starch. The first step of the sugar extraction consisted of adding 10 ml of ethanol (80% v/v) to 50 mg of sample, which we then left for 30 min at 60 • C in water bath, and then we centrifuged (NEYA 8, REMI ELEKTROTECHNIK LTD., Vasai, India) the sample for 10 min at 3200 rpm. In the second step we added 50 µl of the supernatant, 450 µl of ethanol (80%), 500 µl of phenol (28%), and 2500 µl of sulfuric acid (96%), we shook the mix and let it stand for 30 min. In the third step we read the absorbance at 490 nm with spectrophotometer (Spectrophotometer UV-1600PC, VWR, Radnor, PA, United States) after removing the supernatant and drying the sample at 70 • C during 16 h.
In the starch extraction, we added 4 ml of sodium acetate (pH 4.5) to the dry sample and left it for 60 min in a water bath (60 • C). Once the sample cooled down, we added 1 ml of Amyloglucosidase (0.5% w/v) and we incubated the mix in the stove for 16 h at 50 • C. We then added sample 50 µl of supernatant, 450 µl of sodium acetate (pH 4.5), 500 µl of phenol (28%), and 2,500 µl of sulfuric acid (96%). We then mixed it and let sit for 30 min, and then we measured the absorbance at 490 nm with the spectrophotometer.
We analyzed nitrogen concentration in an elemental analyzer (Carlo Erba 1110 Elemental Analyzer) at the University of Wyoming following previously published procedures (Hoffman et al., 2019).

Statistical Analyses
We examined statistical differences across treatments in g MLD , g n , and g d using an ANOVA (followed by Tukey's HSD test) with species, light and water treatments as explanatory variables. Measurements of g MLD , were conducted on different individuals within a treatment. Consequently, we examined whether values were comparable within a given treatment by examining variation in the mean ± 95% CI in g MLD , g n , and g d .
To examine potential drivers of variation in g res , we additionally performed correlation analyses between conductance and NSC, LMA, gas exchange parameters and md . All data was analyzed with R 3.6.3 (R Core Team, 2020) using base packages and, additionally, "corrplot" for plotting the correlation table (Wei and Simko, 2017).

Modeling
In order to examine the effects of the different forms of measuring residual conductance over stomatal predictions and coupled photosynthetic responses, we performed two exercises. First, we examined the effects on stomatal predictions on different implementations of Eqs 1, 2. Second, we examined the effects of the different measured values of g res on a photosynthesis-stomatal conductance coupled model.
For the first exercise, we compared the performance of different versions of the Ball-Berry (BB) model (Ball et al., 1987). First, we examined the version proposed by Duursma et al. (2019, BBD): and we used three different forms of the left hand term [max (g 0 , g min )]. That is, we compared model performance when the left hand term used g 0 and g min estimated from g MLD (BBD MLD ), g n (BBD n ), and g d (BBD d ). In all cases, g 0 was defined as conductance (g MLD , g n , or g d , depending on the case) in the shade treatment without water stress (P 0 _shade) and g min as conductance in the sun under strong water stress (P 80 _sun). We compared these results with the original version of the Ball-Berry model (BB): where g int and m were both estimated through least squares fitting. Finally, we used an intermediate option where we used Eq. 5 but where g int was replaced by actual g MLD measurements (BB_meas_g MLD ), instead of being estimated through least squares. We also tried with g n , in addition to g MLD , but differences were negligible, as will be discussed later in more detail.
Model calibration was performed with data collected under growth PAR (1,500 µmol m −2 s −1 for sun treatment and 40 µmol m −2 s −1 for the shade treatment). Model validation was performed with data collected under different PAR levels. That is, with PAR at 40 µmol m −2 s −1 for the sun treatment and at 1,500 µmol m −2 s −1 for the shade treatment. Model comparison was performed calculating the Akaike Information Criterion (AIC) and a model was considered more plausible when the AIC was smaller by a difference of 2 or more units (Burnham and Anderson, 2002). We also examined the variation in the slope, intercept and R 2 of the observed vs predicted relationship.
For the second exercise, we simulated the effects of the different values of g MLD , g n and g d on predictions of C i with varying PAR and on the effect temperature on leaf evaporation. We used the A/C i curves to parameterize a coupled photosynthesis model (Duursma, 2015) and we conducted the simulation following previously published protocols (Duursma et al., 2019). We note that differences in mesophyll conductance across species and treatments could affect estimates of photosynthetic parameters (Flexas et al., 2012).

RESULTS
Effects of Shade and Drought on g MLD and g n We observed that g MLD varied significantly with species and light and also with light and water (Table 1 and Figure 1). The interactions between species and light resulted in g MLD significantly declining from 6.9 in the sun to 3.4 mmol m −2 s −1 shadow in Q. faginea. However, g MLD in Q. ilex did not differ across light levels (5.6 in the sun and 4.4 mmol m −2 s −1 in the shade). The interaction between light and water was such that g MLD declined with drought in the sun treatment (from 7.4 at P 0 to 5.5 mmol m −2 s −1 at P 80 ), but g MLD increased with drought in the shade from 3.1 at P 50 to 5.0 mmol m −2 s −1 at P 80 ).
Variation in g n followed a pattern of variation similar to that of g MLD in that it also varied significantly with species and light treatments (Table 1 and Figure 2). g n was not different between species at the shade treatment (4.5 and 5.6 mmol m −2 s −1 in Q. faginea and Q. ilex, respectively), but there was a significant increase in g n in Q. faginea (7.8 mmol m −2 s −1 ) in the sun  Table on the effects of species, light treatment, water treatment on residual conductance measured from the mass loss of detached leaves (g MLD ), from nocturnal conductance (g n ), and also from daytime conductance (g d ).

Factor
Df F P-value treatment. Instead, g n in Q. ilex at the sun treatment was similar to that in the shade (4.0 mmol m −2 s −1 ; Figure 1B). Differences across water treatments were not significant.

Effects of Shade and Drought on g d
Of particular relevance for this study is to examine g d when A net approaches zero (Figure 3B), so that one can test the potential use of g d as an indicator of residual conductance. There are different definitions in the literature as to what is meant by photosynthesis approaching zero (Leuning, 1995;Barnard and Bauerle, 2013).
Here we used g d when A net was at, or below, 1 µmol m −2 s −1 . In Q. faginea, this occurred under the shade treatments at all water stress levels, where g d varied between 14.6 and 29.5 mmol m −2 s −1 (Figure 3). In Q. ilex, A net was always below 1 µmol m −2 s −1 in the shadow treatments at all water stress levels. However, there was significant variation in g d as it varied from 58 mmol m −2 s −1 in P 0 to 14 and 4 mmol m −2 s −1 in P 50 and P 80 , respectively. Within the sun treatments, A net was always below 1 under water stress (at P 50 and P 80 ) where g d varied between 4 and 1 mmol m −2 s −1 , respectively. g d under water stress (P 50 and P 80 ) was not different between shadow and sun treatments (Figures 3A,B).
Differences Between g MLD , g n , and g d Within a given treatment, g n was indistinguishable from g MLD : 95% CI error bars always overlapped ( Figure 3C). In Q. faginea, values of g MLD were usually below those of g n , but the absolute difference was less than 4 mmol m −2 s −1 . In Q. ilex, the difference between g n and g MLD was less than 1 mmol m −2 s −1 .
In contrast, g d was consistently and significantly above both, g MLD and g n in Q. faginea. It should be noted that, for this comparison, we only used g d when A net was below 1 µmol m −2 s −1 . That is, we did not seek to compare values of g d with g MLD and g n if A net was above 1 µmol m −2 s −1 because, in that case, photosynthesis does not tend to zero. The average difference of g d with g MLD was 17.7 mmol m −2 s −1 and the average difference of g d with g n was 10 mmol m −2 s −1 . The only case in which g d was not different from g MLD and g n was in the sun treatments in Q. ilex.

Correlates Explaining Variation in g MLD
Overall, the relationships between the different indicators of g MLD and other physiological parameters were species-specific ( Figure 4A). The only exceptions were N mass and NSC concentrations which had a negative and a positive correlation, respectively, with g MLD in both species (Figure 4). In turn, N mass correlated negatively with NSC concentrations and LMA in both species. NSC also correlated with LMA in both species, albeit positively. In Q. faginea, g MLD and g n also correlated positively with LMA and g MLD also correlated positively. In Q. ilex, g n showed a negative correlation with respiration (R) and a positive correlation with md and with A net. NSC concentrations were negatively affected by the shade treatment (Supplementary Figure 1A).

Modeling g d : Comparing Different Formulations of the BB Model
We first compared the performance of the model proposed by Duursma et al. (2019) when g 0 and g min had been defined on the basis of g MLD (BBD MLD ), of g n (BBD n ), and of g d (BBD d ). In all cases, the original g 0 and g min were defined as the level conductance (g MLD , g n , or g d , depending on case) in the P 0 _shade treatment (low light) and in the P 80 _sun treatment, respectively, (high water stress).
Model performance was superior when the model was based on g MLD (BBD MLD ), but differences with the model based on g n were minor ( AIC = 0.3 for Q. faginea and 2.3 for Q. ilex). However, the model based on g d showed consistently a larger AIC, indicating smaller plausibility ( Table 2).
We compared the performance of these three models against the original Ball-Berry (BB) and we observed that BBD MLD and BBD n performed better only in Q. faginea, where the difference in AIC was bigger than 4. For Q. ilex, however, the AIC was similar across models although the intercept of the observed vs predicted relationship was significantly different from 0 only in the BB model.
Finally, we compared the performance of the Ball-Berry model but where, instead of fitting g int through least squares, we use actual g MLD measurements (BB_meas_g MLD ), which we defined originally as g MLD under water stress (P 80 _sun). We observed that this was the best model in Q. ilex as it had the smallest AIC although the difference was not significant with BBD MLD ( AIC = 1.64). In Q. faginea, BB_meas_g MLD peformed worst than BBD MLD ( AIC = 2.6).
Differences between BB_meas_g MLD were significant with the BB model (AIC = 2 for both species) and it was also more plausible than the BBD d model in Q. faginea (AIC > 2). Differences between BB_meas_g int and the other models were not significant. We tried fitting BB_meas_g int with different values of g int (e.g., using values under shade, or from g n ), but differences were not significant (data not shown).

Modeling g d : Coupled Leaf Gas Exchange Model
Depending on how g res was measured, we found significant differences of simulated gas exchange. In particular, when g d was used we always observed higher values of C i at any PAR level and also higher leaf transpiration rates (E l ) as temperatures increased because g d was often larger than g MLD and g n (Figure 5).
Generally speaking, there was little difference in simulated C i and E l regardless of whether g MLD or g n were used, and whether they were defined from g 0 or from g min . The only exception was that, in Q. ilex, there were some differences in predicted C i (particularly at low PAR levels) and in predicted leaf transpiration (particularly at peak E l ) when g res was defined from g n : using g n from the P 0 _shade (g 0 ) treatment led to higher predicted C i and E l than using g n from the P 80 _sun treatment (g min ). It should be noticed that g n from the P 0 _shade treatment was one order of magnitude larger than g n from the P 80 _sun treatment (9.1 vs 0.6 mmol m −2 s −1 , respectively).

DISCUSSION
We observed that residual conductance varied significantly across light and water treatments in an interactive (nonadditive) fashion and the responses differed across species. There were no significant differences as to whether residual conductance was measured from g MLD or from g n , but the values were significantly higher when using g d . g MLD was positively correlated with NSC concentrations, suggesting that further FIGURE 3 | Variation in daytime leaf conductance (A, g d ) and net assimilation (B, A net ) across different light treatments (shadow vs sun) and water treatments (P 0 , P 50 , and P 80 ). Bars indicate mean values per treatment and error bars indicate SE. Hatched bars in (A) indicate that A net in that treatment is significantly higher than 1. (C) Differences in residual conductance as measured by from the mass loss of detached leaves (g MLD ), from nocturnal conductance under high D (g n ), or from g d when A net was smaller than 1 [note that some values of g d are missing in (C) if A net for that treatment was higher than 1]. Error bars indicate 95% CI. reductions in g MLD under drought may be limited by low NSC availability. From a modeling perspective, the small measured differences between g MLD and g n generally did not impact model performance. Although residual conductance differed significantly under experimental treatment, such differences in residual conductance showed only a moderate impact on model performance. That is, model performance did not critically depend upon whether residual conductance was measured under strong shade or under strong water stress. There was also little difference in model fit when either g MLD or g n were used as an absolute minimum in Eq. 4 (BBD MLD or BBD n ), or when they were used as the intercept of the BB equation (BB_meas_g MLD ).

Shade and Drought Interact as Drivers of g MLD Although Responses Are Species-Specific
We observed that g MLD declined under increasing drought in the sun treatment. In the shade treatment, however, g MLD remained low and constant, regardless of the water treatment. This result indicates that drought only affects g MLD under high light because, under shade, light limitations lower g MLD to a minimum that is not affected by water stress. It is worth noting that, at least for some species, full acclimation after changes in the light growth environment may require more than one growing season (Aranda et al., 2001). In other words, the strong limitation imposed by the low light over g MLD may increase even more in subsequent years.
Previous studies had identified how g MLD often decreases under exposure to water stress and light, as a result of changes in wax composition, when each effect is examined in isolation (Shepherd and Wynne Griffiths, 2006). However, our experiment may be the first to examine g MLD responses in a multifactorial experiment. Interestingly, light and water effects were not additive. That is, we did not observe a lower g MLD under low light and high water stress, as would be expected from an additive effect of both factors.
The response to shade was, however, species-specific. g MLD increased in the sun treatment only in the deciduous  2 | Model comparison. We compared for Q. faginea (QF) and Q. ilex (QI) different models based on the Akaike Information Criterion (AIC), the change in AIC relative to the lowest ( AIC) and the R 2 , slope and intercept of the observed vs predicted relationship. For the slope and intercept we show the mean value (and SE). The models compared include the Ball-Berry model (Eq. 5, BB), the BB model where the intercept is measured, rather than estimated (BB_meas_g MLD ) and the modification of the Ball-Berry model proposed by Duursma (Eq. 4,BBD) using g MLD (BBD MLD ), g n (BBD n ), or g d (BBD d ) as g 0 and g min . We also provide the values of the intercept, g int (the maximum between g 0 and g min in the first three models), that were used in each case. Values in bold (or italics) in slope and intercept indicate values significantly (or marginally) different from 1 and 0, respectively. 1 Note that g int in this case is estimated through least squares, rather than measured as in the other models.
FIGURE 5 | Effects of different values of g res on a coupled photosynthesis model as in Duursma et al. (2019). Models are for Q. faginea (A,B) and Q. ilex (C,D). g 0 and g min represent g res under low light (P 0 _shadow treatment) and water stress (P 80 _sun treatment), respectively. Subscripts MLD , n , and d indicate that g 0 or g min (depending on case) were estimated from mass loss of detached leaves, from nocturnal conductance, or from daytime conductance.
Q. faginea, whereas the increase in Q. ilex g MLD under sun was not significant.

g MLD Correlated With Low NSC Concentrations
We can speculate that the reason why g MLD was not lower under the high water stress and shade treatment (relative to other shade treatments under less water stress) is related to carbohydrate limitations. We observed a significant and positive correlation between NSC and g MLD across species. A synthesis of variation in NSC across species reports that a minimum NSC of 46% is always conserved (Martínez-Vilalta et al., 2016). In our results we also observed a minimum NSC that was close to the 46% of the maximum NSC that we measured (Figures 4C,D). A possible explanation on why g MLD did not decrease further under the joint drought and shade stress is related to a lack of NSC to feed the building of additional wax and/or epidermal layers. That is, once plants have reached the minimum NSC threshold of 46% relative to maximum, they will seek to preserve their NSC for other functions, such as osmoregulation, at the expense of building thicker cuticles or additional wax layers. We note that osmoregulation under shade may be impaired in oaks (Aranda et al., 2005;Rodriguez-Calcerrada et al., 2010).
At any rate, this is the first study, to our knowledge to raise this possibility. This result should thus be interpreted with caution. We acknowledge that the correlation between g MLD and NSC may have been affected by jointly considering plants under different light and water regimes. Subsequent work would thus be needed to confirm this hypothesis.
Residual Conductance in Relation to Respiration, LMA and N mass Despite stomatal closure, g res did not limit CO 2 diffusion out of the leaf. In fact, there was a negative correlation between nocturnal conductance and respiration in one of our (Q. ilex) species, indicating higher CO 2 efflux at lower g n and, consequently, that reduced g n levels were far from limiting respiration. This contradicts earlier studies that cytotoxic CO 2 build-up could occur under nocturnal stomatal closure (Fricke, 2019) but it aligns along with the results of modeling, indicating that only under conductances that are orders of magnitude lower to those reported here could a cytotoxic CO 2 build-up occur (Resco De Dios et al., 2019).
g MLD increased in the sun in Q. faginea. LMA also increased with light (data not shown) and it was significantly correlated with g n and with g MLD in Q. faginea. LMA is an indicator of the degree of sclerophylly, which could serve to decrease residual conductance by increasing cuticle thickness. However, LMA also increased with light in the sclerophyll Q. ilex, where LMA did not correlate with g n or g MLD . This result matches with previous studies in Quercus indicating that any effects of LMA in g n and g MLD may be modified by changes in the cuticle composition (Bueno et al., 2020). We note that this argument is speculative and based only on circumstantial evidence. N mass showed a negative correlation with g MLD in both species, and with g n in Q. faginea. This result indicates that species with a higher photosynthetic investment will decrease the investment in residual conductance. This points toward a potential mechanism underlying the trade-off between investment for C uptake (higher N mass ) and preventing catastrophic water losses (reduced g MLD ). Further studies will be necessary to test the generality of this hypothesis.
Residual Conductance May Be Measured via g MLD or via g n Under High D Measurements of nocturnal conductance under the relatively high D from this experiment were statistically indistinguishable from independent measurements of residual conductance indicating that the latter was driving the former. It had been previously argued that measurements of g n would not be valid indicators of residual conductance, because g n is actively regulated (Duursma et al., 2019). Our results suggest that this argument from Duursma et al. (2019) may only be valid when g n is measured under low D.
Previous studies document that stomata often reach complete closure (or as complete as it can be) under lower D in the night, than in the day (Barbour and Buckley, 2007). This phenomenon would explain why g n was much lower than g d although D was comparable, and it is likely explained by the capacity of stomata to sense and open in response to light.
We also show how modeling results were not affected by either using g MLD or g n . This result, however, needs to be interpreted with caution. We only focused on BB-type stomatal models and other results may be obtained in different model types. For instance, as Duursma et al. (2019) noted, changing minimum conductance from 2 to 4 mmol m −2 s −1 halved the time to reach mortality in a hydraulics model because it doubled the water losses (Duursma et al., 2019). Although differences were statistically not significant between g MLD and g n , we still observed differences in mean residual conductance of 4 mmol m −2 s −1 across measurements, indicating that measurement errors and other sources of uncertainty may play a large for other model types, such as mortality models.

Modeling g s
We acknowledge our dataset was limited to thoroughly test the best form of the BB model: we sampled under highly contrasting light and water conditions, but only once in time. We would thus need data over more time periods for a more thorough evaluation. However, our dataset allows for the development of some hypotheses, which may be expanded in subsequent studies.
We observed that there were only little differences between Eq. 5, where the original BB function was used but including measured g MLD (BB_meas_g MLD ), instead of the version proposed by Duursma et al. (2019;Eq. 4). Duursma et al. (2019) note that residual conductance acts as an actual minimum in the function they propose. However, if the goal is to use residual conductance as an actual stomatal minimum, one could consider the following equation instead: g s = max min g 0 , g min , m A RH/Ca (6) where the minimum between g 0 and g min is chosen [not the maximum, as proposed by Duursma et al. (2019)]. At any rate, we did not observe major differences in model performance between the BBD model or Eq. 5. This result indicates that it is unlikely that losses in model performance will derive from the adoption of the alternative model formulations, as proposed by previous studies (De Kauwe et al., 2015;Duursma et al., 2019).
Our results also indicate that g MLD and g n can both be interchangeably, and that the choice between g 0 and g min exerts negligible consequences for model fitting. Earlier studies indicate a major effect of g res (De Kauwe et al., 2015;Duursma et al., 2019). This is because those studies used a wide range of g res values (10-40 mmol m −2 s −1 , depending on the case), much higher that the variability we reported here when using g MLD and g n across treatments ( Table 2). Synthesis studies similarly indicate limited variation in g res within a family (Duursma et al., 2019). After discarding g d as a reliable indicator of g res , our results indicate a minor effect of different methods and approaches used for measuring g res and for modeling water use, at least in our two closely related species.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
VR designed the experiment. CA performed research with the help of FC. VR analyzed the data and wrote the manuscript with important feedback from all co-authors. All authors contributed to the article and approved the submitted version.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2020. 603581/full#supplementary-material Supplementary Figure 1A | Variation in leaf non-structural carbohydrate concentrations measured across different light treatments (shadow vs sun) and water treatments (P 0 , P 50 , and P 80 ) in Quercus faginea (QF) and Q. ilex (QI). Bars indicate mean values per treatment and error bars indicate SE.
Supplementary Table 1A | Target midday water potential ( md ) to reach the desired PLC according to Esteso-Martínez et al. (2006) for Q. faginea and to Peguero-Pina et al. (2014) for Q. ilex, and actual values. Mean (and SE) actual values are presented. The letters in "Actual PLC" indicate the results of post hoc analyses (Tukey HSD). This is a reproduction (with permission from the publisher) of Table 1 originally published in Resco De Dios et al. (2020).