Using δ13C and δ18O to analyze loblolly pine (Pinus taeda L.) response to experimental drought and fertilization

1Department of Forestry and Environmental Resources, North Carolina State University, 2820 Faucette Drive, Raleigh, NC 27606, USA; 2College of Life Sciences and Oceanography, Shenzhen University, 3688 Nanhai Boulevard, Nanshan District, Shenzhen, Guangdong 518060, China; 3Bordeaux Sciences Agro, UMR 1391 INRA-ISPA, 33195 Gradignan Cedex, France; US Geological Survey, Wetland and Aquatic Research Center, 700 Cajundome Boulevard, Lafayette, LA 70501, USA; Department of Forest Ecology and Management, Swedish University of Agricultural Sciences, Skogens ekologi och skötsel, 901 83 Umeå, Sweden; 6Department of Forest Resources and Environmental Conservation, Virginia Polytechnic Institute and State University, 310 West Campus Drive, Blacksburg, VA 24061, USA; 7Rayonier Inc., 851582 Highway 17N, Yulee, FL 32097, USA; Department of Ecosystem Science and Management, Texas A&M University, 495 Horticulture Street, College Station, TX 77843, USA; Eastern Forest Environmental Threat Assessment Center, United States Department of Agriculture Forest Service, 3041 East Cornwallis Road, Research Triangle Park, NC 27709, USA; 10Corresponding author (wlin17@foxmail.com) http://orcid.org/0000-0003-3240-1389


Introduction
Projections of future water availability suggest that the southeastern USA will likely experience an overall decrease in precipitation, aswellasincreasingtemporalvariabilitycharacterizedbyfewerbut more intense rain events (Ryan 2011, Intergovernmental Panel on Climate Change (IPCC), 2014, Melillo et al. 2014). The increasing likelihood of droughts between these events will stress forests in a temporally and spatially complex manner (Anderegg et al. 2015, Bansal et al. 2015, Wolf et al. 2016. Loblolly pine (Pinus taeda L.) is the predominant tree species on over 13.4 million ha of southern forest lands and represents one-half of the standing pine volume (Schultz 1997, Wear andGreis 2012).
Due to the large area and high productivity, this species plays a crucial role in the regional carbon budget and economic vitality of the southeastern USA (Schultz1999).To understandthe potential physiological and economic impacts of these projected changes on loblolly pine ecosystem function and management implications, the Pine Integrated Network: Education, Mitigation and Adaptation Project (PINEMAP; Will et al. 2015) was implemented.
Commercial loblolly pine has been bred and selected for maximum economic productivity, aided by intensive management practices. External nutrient inputs, competition and pest control, and site preparation techniques all represent substantial inputs of resources and labor responsible for maximizing final yields and decreasing rotation length (Fox et al. 2007). In particular, mid-rotation fertilization (as opposed to fertilization at stand establishment) has been found to increase growth and productivity in P. taeda (Albaugh et al. 1998, Samuelson et al. 2001 and is a staple of commercial forestry. Given the fundamental controls of plant allometry (Chen et al. 2013), nutrient supplementation simultaneously increases leaf area and decreases root area, which are expected to narrow the trees' safety margin under severe water deficit, and increases the risk of hydraulic failure due to catastrophic cavitation-induced embolism (Bakker et al. 2009, McNulty et al. 2014. Thus, quantitative understanding of the interaction between nutrient and water availability on plant physiological processes is critical for understanding carbon sequestration responses and projecting plant resilience in a drier environment.
Both drought and fertilization may influence the leaf-level gas-exchange processes (Farquhar et al. 1989, Brooks andCoulombe 2009). An early indicator of plant water stress is a decrease in stomatal conductance (g s , see Table S1 available as Supplementary Data at Tree Physiology Online for symbols and abbreviations), which responds to the water potential gradient along the plant hydraulic pathway (Chaves et al. 2002, Flexas andMedrano 2002), whereas its limitation on photosynthesis (A) is indirect and occurs under greater water stress (Flexas et al. 2004). However, physiological (A) and structural (canopy volume and leaf area) changes in loblolly pine after fertilization have been shown to vary by site, time since fertilization and genotype (King et al. 2008, Maier et al. 2008, Pell 2015. Similarly, the reported change of canopy conductance (G s ; canopyaveraged g s ) to fertilization has been variable (Ward et al. 2013, 2015, Bartkowiak et al. 2015, Wightman et al. 2016. Such physiological responses are documented in the isotopic composition of tree ring cellulose. Intrinsic water-use efficiency (iWUE), the leaf-level ratio of A to g s , can be inferred from stable carbon isotope ratio (δ 13 C), assuming that mesophyll conductance is relatively constant (Farquhar et al. 1989, Warren and Adams 2006, Flexas et al. 2008). In addition, oxygen isotope composition provides information about g s (Barbour 2007). Thus, pairing δ 13 C and δ 18 O data would allow differentiation between photosynthesis-and stomatal conductance-driven shifts in iWUE (dual isotope approach, Scheidegger et al. 2000). However, the applicability of the dual isotope method remains a matter of debate as the assumptions may not always be met (Scheidegger et al. 2000, Roden andSiegwolf 2012). Alternatively, mechanistic models have allowed direct simulations of tree ring cellulose oxygen isotope discrimination ( 18 O c ) with an array of climate and physiological variables (Roden et al. 2000, Barbour and Farquhar 2000. Such models provide yet another way of evaluating the environmental impacts on 18 O c through sensitivity analysis (e.g., Brooks and Mitchell 2011).
Sensitivity analysis has been utilized widely and has shed light on the ecophysiological processes involved and the environmental changes that lead to the change in 18 O c (Roden et al. 2000, Brooks and Coulombe 2009, Ogée et al. 2009, Brooks and Mitchell 2011. It is performed to evaluate the effects of individual parameters while keeping the others unchanged (one-dimensional sensitivity analysis thereafter). However, due to the complexity of the equations and interactions between the parameters, the relative importance of parameters defined can depend on the values of other parameters. Therefore, it is necessary to understand model behavior across all possible parameter space, and multi-dimensional sensitivity analysis, a way to explore model behavior by varying all parameters together (Williams and Yanai 1996), can then be utilized to achieve this goal.
In the current study, we evaluated the physiological responses of loblolly pine to drought and fertilization at the Virginia replicate of the PINEMAP Tier III experimental sites , at the northern edge of this species' distribution range. With 13 C compositions of tree ring cellulose, in conjunction with trunk sap flow and allometry data, we partitioned the productivity responses to structural (leaf area) versus physiological (A) components. Multi-dimensional sensitivity analysis for the tree ring cellulose oxygen isotope discrimination models, on the other hand, explored the model behavior across growing scenarios and provided information about source water change.

Site description and experimental design
The study site was located in Buckingham County, Virginia (37 • 27 37 N, 78 • 39 50 W), in the Piedmont physiographic region, with 0-15% slopes. This site represented the northernmost of the four drought and fertilization experiments established under the PINEMAP study . The study was established in 2012 in a production forest, and the understory vegetation was mechanically cleared. A thick and uniform litter layer (about 5 cm) of dead needles remains across the site. The experimental set-up was a randomized complete block factorial design with two levels of water availability and two levels of nutrient availability with four blocks. There are 30 to Tree Physiology Volume 00, 2019 Downloaded from https://academic.oup.com/treephys/advance-article-abstract/doi/10.1093/treephys/tpz096/5636870 by DigiTop USDA's Digital Desktop Library user on 07 January 2020 53 trees for measurements within each of the 16 plots. All trees survived by 2016 in seven plots, while no more than two trees died in each of the rest of the plots. The treatments were control (C, ambient throughfall without fertilization), drought (D, 30% throughfall reduction without fertilization), fertilization (F, ambient throughfall with fertilization of a complete suite of nutrients) and fertilization plus drought (FD, throughfall reduction with fertilization). The drought treatment consisted of removing ∼30% of precipitation throughfall by troughs positioned below the live crown, alongside tree rows. The troughs were installed on both sides of a row of trees, each 0.6 m wide and about 0.3 m from each other. They ran the length of the entire plot, including border rows, and the effluent was directed downhill from all plots using drainage pipes. Given that the troughs were ∼1 m above the ground, plots were small (0.11 ha square plot, with 32.9 m at each side), and the drought treatment plots were adjacent to and intermixed with ambient rainfall plots, vapor pressure deficit (VPD) and wind speed were similar among plots. The radiation level that reached the ground was not significantly affected because the canopy was closed when the experiment started. Therefore, we assume that the effect of troughs on soil evaporation was low, even though it was not directly quantified. Fertilizer was applied at a rate of 224 kg N ha −1 , 27kg P ha −1 , 52kg K ha −1 and 1.12 kg ha −1 of micronutrient mix (6% sulfur, 5% boron, 2% copper, 6% manganese and 5% zinc; Cameron Chemical, Virginia Beach, VA, Southeast Blend) at the start of the study in April 2012. Trees represented the local open-pollinated seed source. Additional information on the experimental set-up was reported by Will et al. (2015).
The soil is a well drained, fine, mixed, subactive, mesic Typic Hapludult of the Littlejoe soil series, with a silt loam surface soil and a clay loam subsoil. The depth to ground water is > 200 cm. Mean annual precipitation is 1120 mm and mean annual temperature is 13.6 • C. The loblolly pine trees were planted as 2-year-old seedlings in 2003, and the treatments started in April 2012. Precipitation, photosynthetically active radiation, relative humidity (RH) and air temperature (T a ) were measured using micrometeorological sensors mounted above the canopy at the center of the study site. The VPD was calculated using RH and T a after Campbell and Norman (1998). Additional information on site micrometeorology was reported by Ward et al. (2015). Daily averages of wind speed and atmospheric pressure recorded at the Lynchburg International Airport weather station, VA, ∼50 km away from the study site, were obtained from National Centers for Environmental Information (https://www.ncdc.noaa.gov/data-access).

Estimation of aboveground biomass increment and mean daily canopy conductance
Tree height and diameter at breast height (DBH) of all trees in each plot were measured annually during dormant seasons beginning in 2012. The aboveground biomass increment (AGBI) of each measurement tree was calculated using tree height, DBH, age and stand density after Gonzalez-Benecke et al. (2014).
Leaf area index (LAI) was measured using a LAI-2200C Plant Canopy Analyzer (Li-Cor, Lincoln, NE, USA) on 1 April and 25 August 2015 in each plot. Monthly LAI and DBH between April and September were estimated by linear interpolation, approximating the seasonality reported for earlier years at the same site . Monthly tree transpiration on a daily basis (kg_water day −1 ) in 2015 was estimated from sap flow measurements using thermal dissipation probes after Granier (1985). The details on the approach were reported by Ward et al. (2015). As leaf area within a given stand is roughly proportional to sapwood area (Domec et al. 2009, Tor-ngern et al. 2017, the monthly leaf area of the measurement tree was estimated using plot size, monthly LAI and monthly ratio of the sapwood area to the sum of sapwood area from the whole plot. Sapwood area was estimated based on DBH using the allometric relationship reported by Gonzalez-Benecke and Martin (2010). Finally, the mean daily canopy conductance (mean daily G s , mol m −2 day −1 , leaf area basis) of each measurement tree from May through September was estimated using tree transpiration, monthly VPD and monthly leaf area .

Core sampling, cellulose extraction, precipitation collection and stable isotope analysis
Increment cores of 5 mm and 12 mm diameter were collected from three out of the five representative trees per plot that were equipped with sap flow sensors. The cores were collected from the northern side of the trees at breast height (1.3 m) in January 2016, four growing seasons after the start of the drought and fertilization treatments. The cores were oven dried at 60 • C, mounted and sanded progressively down to 1200 grit. The ring width was measured to the nearest 0.001 mm using a tree ring measuring system (Velmex Inc., Bloomfield, NY, USA) and the cores were cross-dated using the program COFECHA (Holmes 1983).
Latewood produced in 2015 was separated from the 12 mm cores, and α-cellulose was extracted (Wieloch et al. 2011, Lin et al. 2017. The 13 C and 18 O stable isotope ratios of αcellulose were determined at the Cornell University Stable Isotope Laboratory (http://www.cobsil.com) using a Thermo Delta V isotope ratio mass spectrometer (Thermo Fisher Scientific, Bremen, Germany) interfaced to a NC2500 elemental analyzer (Carlo Erba, Italy) and to a Temperature Conversion Elemental Analyzer (Thermo Fisher Scientific, Bremen, Germany). The within-run isotopic precision of the methodology using quality control standards was 0.2 for carbon and 0.4 for oxygen (K. Sparks 2015, personal communication).
Rainwater was collected at the study site on a monthly basis during the experimental period using a carboy fitted with a funnel. A small amount of mineral oil was applied to Tree Physiology Online at http://www.treephys.oxfordjournals.org reduce evaporation in the field. Oxygen isotope ratios (δ 18 O) were determined at the Stable Isotopes for Biosphere Science laboratory at Texas A&M University (http://sibs.tamu.edu) using a high temperature conversion elemental analyzer coupled to an isotope ratio mass spectrometer (TC/EA-IRMS; Delta V, Thermo Scientific, Waltham, MA, USA). Given the deep ground water table (>200 cm), precipitation was assumed to be the major source water for this site. Since latewood cellulose was used for isotopic analysis, we considered precipitation from May through September a good proxy for source water in our site. However, the water sample collected from 29 May to, and 16 June was missing. Thus we calculated precipitation δ 18 O (δ 18 O p ) from four periods in 2015: 17 June-16 July, 17 July-21 August, 22 August-16 September and 16 September-6 October. The average δ 18 O p was weighed by the amount of rain water collected during each period.
Carbon isotope discrimination ( 13 C c , ) was calculated as: where δ 13 C c is the 13 C isotopic ratio of the extracted α-cellulose and δ 13 C a is the isotopic ratio of the air. Given that VPD at the site was mostly within low to medium levels (<2 kPa), boundary layer conductance was relatively high (6.4 mol m −2 s −1 ), and that evidence for variable mesophyll conductance across genotypes and water availabilities is lacking for this species (Wilson 2014), iWUE was calculated as: where a is the discrimination against 13 CO 2 during diffusion through the stomata (4.4 ), b is the net discrimination due to carboxylation (27 ) and c a is ambient CO 2 concentration. The annual averages of c a and δ 13 C a were calculated from monthly atmospheric CO 2 concentration and δ 13 C a values obtained from the sampling station of Scripps CO 2 Program at La Jolla Pier, CA, USA (http://scrippsco2.ucsd.edu/data/atmospheric_co2/ljo). Oxygen isotope discrimination of cellulose ( 18 O c , ) was calculated as: where δ 18 O c is the 18 O isotopic ratio of the extracted α-cellulose and δ 18 O p is that of the collected precipitation at the study site. Because the plots were spatially collocated and trees were of the same species and age, we began by assuming that the source water δ 18 O across treatments was similar, and that the differences in 18 O c reflected the difference in the evaporative process at the leaves affected by different treatments. We then use sensitivity analyses to test this assumption.

Sensitivity analyses of 18 O c models
Similar to Brooks and Mitchell (2011), we evaluated the sensitivity of tree ring cellulose oxygen isotope discrimination ( 18 O c ) contributed by three component processes: evaporative enrichment at evaporative sites (Farquhar et al. 2007), the Péclet effect describing 18 O in the bulk leaf water ) and the fractionations in cellulose biosynthesis (Barbour and Farquhar 2000).
The evaporative enrichment at the evaporative sites ( 18 O e ) was calculated with the modified Craig-Gordon model (Farquhar et al. 2007): (4) where + is equilibrium fractionation coefficient and a function of leaf temperature (T l ), k is the kinetic fractionation coefficient and a function of the ratio of g s and boundary layer conductance (g b ), w a is the vapor pressure of the atmosphere and a function of T a and RH, w i is the saturated vapor pressure inside leaves and a function of T l and 18 O v is the enrichment of atmospheric water vapor above source water and is equal to − + under wellmixed outdoor conditions (Barbour 2007, Song et al. 2014a. As the variation of 18 O v across the site was considered negligibly small due to the small plot size and mixing of air, it was calculated using T a instead of T l for each of the simulation scenarios in the sensitivity analysis. The average isotopic enrichment of bulk leaf water at steady state ( 18 O L ) was modeled by considering the Péclet effect (℘): where ℘ is the Péclet number, a function of effective path length (L), transpiration (E), molar density of water (C ) and diffusivity of H 18 2 O (D ; Eq. 6). E is a function of g s , g b , atmospheric pressure (p a ) and leaf-to-air VPD (Song et al. 2014a, Eq. 7). D is dependent on T l , while L is usually an empirically fitted value (Barbour 2007).
Finally, 18 O c was estimated following Barbour and Farquhar (2000) as: tion of source water at the site of cellulose formation (thus its value equals 1 for trunk wood) and ε wc is the isotope fractionation of 27 , associated with the oxygen atom exchange between water and carbonyl groups.
Equations 4-8 require T a , T l , RH, g s , g b , L and p ex as input parameters. The parameters L and p ex were assumed constant across the treatments. L was set to 33 mm, corresponding to mean values reported for conifers by Brooks and Mitchell (2011), while p ex was set as 0.4, similar to Cernusak et al. (2005) and Keel and Weigt (2016). Using an exponential L-E relationship as reported by Song et al. (2013) did not improve the fit between simulated values with observations, except when paired with the p ex value used for Pinus rigida (Song et al. 2014a). Considering the lack of measurements in loblolly pine and the similarity of the trends exhibited in sensitivity analyses by these two approaches, we decided to use the values of L as reported by Mitchell (2011) andD. Yakir (2017, personal communication). Although p ex has been reported to vary within and across species (Gessler et al. 2009, Song et al. 2014b, Cheesman and Cernusak 2017, we used a fixed value considering that: (i) the majority of studies report p ex ≈0.4 (Gessler et al. 2014), (ii) current evidence for the effect of longterm water supply conditions on p ex of the same species is weak (Song et al. 2014b, Cheesman andCernusak 2017) and (iii) to the extent that it could vary, the within-treatment variability of p ex probably exceeded that between treatments (Song et al. 2014b;X. Song, 2018, personal communication). g b was calculated with the growing season mean wind speed (2.0 m s −1 ) and needle diameter (d) for loblolly pine (1.5 mm; Campbell and Norman 1998). As a result, the parameters T a, RH, g s and leaf-to-air temperature difference (T l -T a , to decouple T l from T a ) were selected for one-and multi-dimensional sensitivity analyses.
One-dimensional sensitivity analysis was performed first to evaluate the effects of individual parameters. The reference conditions for this analysis were chosen as the growing season (May-September 2015) average daytime (6:00 h to 21:00 h) T a (24.3 • C), RH (69%) and G s measured in the control plots in 2013 (67 mmol m −2 s −1 , unpublished data). Although T a and RH are expected to be similar across the plots, all four parameters were allowed to vary for better comparison to that of the multi-dimensional sensitivity analysis.
Multi-dimensional sensitivity analysis was conducted to explore model behavior by varying all four parameters (T a, RH, g s and T l -T a ) together. We limited the parameter space, in this case a 4-dimensional hypervolume, to realistically observed or modeled ranges. Firstly, only the observed joint parameter space of T a and RH, as measured at the site (Figure 1a), was used. For instance, at air temperature of 30 • C, RH ranged from 30% to 80% in the sensitivity analysis. A second effort to limit the parameter space was to set an upper limit for g s using the relationship between G s and VPD (Eq. 9, the solid line in  Conifer needles are usually well coupled with the atmosphere (Jarvis and McNaughton 1986), yet T l -T a as high as 5 • C have been noted (Martin et al. 1999, Aubrecht et al. 2016). Thus we simulated the ranges of T l -T a using Eq. 10, following Campbell and Norman (1998).
where γ * is apparent psychrometric constant and a function of g s , g b , wind speed (u) and d; s is the slope of saturation mole fraction function and a function of T a and p a , R ni is radiation load and a function of absorbed radiation, emissivity of surfaces and the Stefan-Boltzmann constant, g Ha is the boundary layer conductance for heat and a function of u and d, g r is radiative conductance and a function of T a and c p is specific heat of air at constant pressure. Our analysis of Eq. 10 indicated that T l -T a was sensitive to radiation load and low wind speed. That is, the range of T l -T a was calculated using minimal and maximal R ni as well as minimal u. The details of the simulation are included in Note S1 available as Supplementary Data at Tree Physiology Online.
Tree Physiology Online at http://www.treephys.oxfordjournals.org Given the small plot size and proximity of different treatment plots, T a and RH are considered the same among plots and treatments. Thus, any change in 18 O c in the treatment plots was attributed to changes in g s and T l -T a due to our treatments.

Statistical analysis
Statistical analyses were performed in R software (version 3.3.2, R Development Core Team 2016). Two-way ANOVA was used to evaluate the effects of drought and fertilization on AGBI, iWUE and 18 O c . Due to a single outlier in iWUE, which was more than three times the interquartile range above the third quantile, Type III sum of squares was considered in ANOVA of iWUE, without considering the block effect. As the number of trees with available sap flux measurements varied among plots, we used a mixed-effect model to evaluate the effects of drought and fertilization on mean daily canopy conductance using 'lme4' package (Bates et al. 2015). Drought and fertilization were treated as fixed factors, while block and month were considered random factors. The P-values were obtained using Satterthwaite's approximation in 'lmerTest' package (Kuznetsova et al. 2017). Given that LAI is a stand-level variable and was measured twice in 2015, the same procedure was applied to LAI measurements, except that measurement date was treated as the only random factor.

Aboveground biomass and leaf area
Four growing seasons after implementing the experimental treatments, fertilization significantly increased growth of loblolly pine. Although tree heights were similar across treatments (P = 0.38 for drought effect and P = 0.24 for fertilization effect, n = 12), LAI ( Figure 2a) and AGBI (n = 12, Figure 2b) were consistently greater in fertilized plots at both ambient and reduced water availability (P = 0.05 for both LAI and AGBI). Their response to drought were not statistically significant (P = 0.60 and P = 0.34 for LAI and AGBI, respectively), and there were no significant interactions between treatments. On average, peak LAI under fertilization was 3.45 m 2 m −2 , or 12% higher than in unfertilized plots, while the average AGBI in 2015 was 21.9 kg per tree under fertilization, 29% higher compared with the reference treatment (Figure 2b).

Canopy conductance, iWUE and observed 18 O c
Mean daily canopy conductance (mean daily G s ), calculated from tree transpiration, ranged from 0.5 to 11.3 kmol H 2 O m −2 day −1 from May through September. Due to the large variance among individual trees (e.g., the coefficient of variation was 0.52 for D treatment and 0.66 for the control), the treatments did not significantly differ. However, Cohen's d was 0.45  (iWUE) and (e) the observed tree ring cellulose oxygen isotope discrimination ( 18 O c ). * Represents significance at 5% level, * * * represents significance at 1 level, while ns represents not significant.
Intrinsic water-use efficiency was 85.6 ± 1.2 μmol mol −1 (mean ± SE, n = 24, see Table S2 available as Supplementary Data at Tree Physiology Online) in C and F treatments and 91.1 ± 1.0 μmol mol −1 in drought treatments (D and FD, n = 24, Figure 2d). Fertilization had no effect on iWUE (P = 0.99, n = 12), while the effect of drought was significant (P = 0.001, n=12). There was no interaction between the fertilization and drought treatments, with drought affecting iWUE similarly at both ambient and elevated nutrient availability. The observed oxygen isotope discrimination ranged from 33.4 to 36.3 , with a significant drought-induced decrease at both ambient and operational nutrient supply (n = 12, Figure 2e, see Table S2 available as Supplementary Data at Tree Physiology Online).

O c simulation and sensitivity analyses
The simulated tree ring cellulose oxygen isotope discrimination was calculated for 2015 growing season mean T a (24.3 • C) and RH (69%), and G s from the control plots in 2013 (67 mmol m −2 s −1 ), yielding 33.4 . With summer rains (from June through September of 2015) as source water, which had a weighted average δ 18 O of −5.0 , the observed 18 O c in the control treatment was 35.4 or 2.0 higher than the modeled estimate.
The one-dimensional sensitivity analysis (see Figure S1 available as Supplementary Data at Tree Physiology Online) indicated that except for T l -T a , any increases of T a , RH and g s would decrease the simulated 18 O c . The effects of T l -T a and g s on 18 O c were similar in magnitude but opposite in sign.
The responses of 18 O c to changes in g s and T l under a range of scenarios of possible T a and RH combinations are shown in the multi-dimensional sensitivity analysis (Figure 3). The parameter space was then restricted to observed environmental conditions and corresponding ranges of g s and T l -T a (Eq. 9 and 10, respectively); the unlikely parameter combinations are not discussed. Similar to the one-dimensional sensitivity analysis, the relationship between 18 O c and g s decreased monotonically, although 18 O c was more sensitive to changes in g s when RH was low and T l -T a was high (Figure 3d and i). In contrast, 18 O c increased monotonically with T l -T a in most scenarios. Notably, this relationship could become reversed at low RH and high g s (i.e., the areas next to those shaded in dark gray in Figure 3d, e, i and j). However, such conditions occurred less than 4% of the time during the growing season of 2015.

Discussion
No interactions between treatments were detected for any of the response variables. Therefore, the direct effects of drought and fertilization are discussed separately.

Tree responses to drought
The 30% reduction of throughfall led to an average of 8.4% (24.6 ± 0.32% vs 16.2 ± 0.02%) decline in soil water content in the top 15 cm (Bracho et al. 2018) and increased the iWUE of loblolly pine by 7.6% (P < 0.01) at both ambient and elevated nutrient availability (Figure 2d). Although ANOVA did not deem the drought response in g s statistically significant due to high among-tree variability, Cohen's d = 0.45 showed the treatment response (21% below reference) to be of intermediate magnitude. More importantly, the increasing of iWUE indicates a decline in g s in response to drought. Given the muted responses of LAI and AGBI (Figure 2a and b) to drought, photosynthesis (A) must have remained relatively unaffected. As iWUE represents the leaf-level ratio of A to g s , the increased iWUE under drought suggests a decline in g s . Therefore, we consider the decline in g s observed in the current study as physiologically significant.
The multi-dimensional sensitivity analysis of 18 O c models ( Figure 3) helps us understand the physiological responses of loblolly pine to drought. The relationship between 18 O c and T l -T a increases monotonically, except at low RH and high g s (Figure 3d). However, given that such combination of T a and RH occurred only rarely (<4% of the time in the 2015 growing season), it is unlikely thatsuch a reversedrelationshipwouldaffectthesimulated 18 O c substantially. Therefore, when we focused on a representative responsesurfaceof 18 O c tog s andT l -T a (Figure 4),thedecrease of g s was expected to increase 18 O c (vertical yellow arrow in Figure 4). As T l -T a probably increases due to the decrease of g s , 18 O c is also expected to increase (horizontal yellow arrow in Figure 4). Yet, the observed 18 O c under drought was actually significantly lower than in non-drought plots (Figure 2e), directly opposite to the model prediction (Eqs 4-8, Figure 4). Given that the effective path length has a limited influence on 18 O c (Ogée et al. 2009) and that drought had a limited effect on tree growth in our current study, the most plausible explanation for the observed change in 18 O c is a change in the δ 18 O of source water. As to where exactly the 18 O-depleted water in the drought treatment came from, we can only speculate, because soil water was not sampled during this study. The winter precipitation (366 mm from November 2014 through March 2015), released mostly as snow thaw, amounted to 180% of the water-holding capacity of the soil (as per gSSURGO (Gridded Soil Survey Geographic Database, https:// gdg.sc.egov.usda.gov/)), and could thus easily have reset the soil water δ 18 O values by the start of the 2015 growing season, including the drought treatment plots. Due to the input of more enriched summer precipitation (see Figure S2 available as Supplementary Data at Tree Physiology Online) and evaporation of surface soil water, it is expected that surface soil water was more 18 O-enriched than subsoil water at our study site during the 2015 growing season. In particular, surface soil water Tree Physiology Online at http://www.treephys.oxfordjournals.org Figure 3. The multi-dimensional sensitivity analysis for simulated oxygen isotope discrimination of tree ring cellulose (simulated 18 O c ), varying with air temperature (T a ), relative humidity (RH), leaf-to-air temperature difference (T l -T a ) and stomatal conductance (g s ). The area shaded in dark gray indicates unlikely simulation results due to the canopy conductance-VPD relationship described by Eq. 9, while the area shaded in light gray indicates simulation beyond the limits of T l -T a . Figure 4. The simulated oxygen isotope discrimination of tree ring cellulose (simulated 18 O c ) with varying leaf-to-air temperature difference (T l -T a ) and stomatal conductance (g s ) at site mean air temperature (T a = 24.3 • C) and relative humidity (RH = 69%) of 2015 growing season. The arrows in the plot indicate the directions of changes in g s and T l -T a under shifting water and nutrient availability. of the drought plots was likely more enriched in 18 O due to evaporation, by a simple calculation using Rayleigh equation. With comparable leaf area and 30% less new water input, it is easy to see that the drought plots would deplete the newly available surface water sooner than non-drought treatments, and depend more on water stored from previous rain events, as reported previously (Tang and Feng 2001). It is also likely that water from the earlier precipitation events, with lower δ 18 O, is closer to the surface in drought than in non-drought plots, further contributing to a possible shift in source water. In a similar study with precipitation reduction from an arid site in New Mexico, USA, Grossiord et al. (2017) reported seasonal water source change for both piñon pine (Pinus edulis) and oneseed juniper (Juniperus monosperma) by sampling xylem and soil water directly. Even though the current study was conducted in a mesic ecosystem, in a much wetter climate (see Site description and experimental design), the δ 18 O of latewood αcellulose suggested a similar change of drought-exposed trees to a different water source.

Tree responses to fertilization
The reported effects of nutrient availability on iWUE in trees have been mixed (Elhani et al. 2005, Balster et al. 2009, Brooks and Coulombe 2009, Brooks and Mitchell 2011, Cornejo-Oviedo et al. 2017. Of the limited measurements of δ 13 C in loblolly pine foliage, the most relevant in the current context is a companion study at a similar experimental site in Georgia, USA that showed a short-lived increase in foliar δ 13 C (Samuelson et al. 2014). This effect was detected only in the year of treatment, which is consistent with our current observation of no response to fertilization in iWUE 4 years into the experiment. In the current study, the mean daily G s did not respond to fertilization (Figure 2c), with a negligible effect size (Cohen's d = 0.007). When combined with the lack of an iWUE response, this suggests that leaf-level A was also not affected by fertilization. Given that both AGBI and LAI increased under fertilization (at both ambient and reduced water availability; Figure 2a and  we conclude that the increase in AGBI was driven by LAI rather than by leaf-level photosynthetic capacity, consistent with earlier reports (Samuelson et al. 2001, Gough et al. 2004. As discussed earlier, g s did not respond to fertilization. However, a smaller mean T l -T a is expected due to greater self-shading from increased LAI under fertilization (Figure 2a). Such change may lead to a decrease in 18 O c (green arrow in Figure 4). Therefore, we expect the simulated 18 O c to remain unchanged or decrease slightly under fertilization. As the observed change in 18 O c was not statistically significant (Figure 2e), the model prediction was consistent with observations.

Implications for tree ring cellulose 18 O studies
While uncertainties remain in capturing fractionations in the post-photosynthesis processes for tree ring cellulose 18 O (Gessler et al. 2014), good overall agreement between simulation and observation of latewood cellulose of evergreen species from temperate zones has been reported (Ogée et al. 2009). In the current study, given the uncertainties in estimating L and p ex , and the assumptions of T l = T a , 18 O v = −ε + and source water δ 18 O equals to precipitation δ 18 O, the difference between simulated and observed 18 O c in the control plots was around 2 , indicating a reasonable simulation. However, the lack of actual source water measurements did not allow us to test our estimate that droughted trees accessed more 18 O-depleted water, and we were not able to evaluate the response dynamics because only one year of tree rings were sampled. Therefore, we advise future investigations to measure all potential sources of water such as soil and xylem water, and to take advantage of the full time series of tree ring isotopic records, including that of pretreatment, when embarking on research involving tree ring stable isotopes like ours.
Although a limited number of parameters are needed for 18 O c models, their impacts on simulation are not immediately obvious. Our one-dimensional sensitivity analyses indicated similar responses to those reported by others (Roden et al. 2000, Brooks and Coulombe 2009, Brooks and Mitchell 2011. In contrast, the multi-dimensional sensitivity analysis explored model behavior across the possible parameter space and shows relationships that are not self-evident (Figure 3). For example, the effect of g s on 18 O c varies with RH, being stronger at low RH (despite being partly counteracted by the negative effect of VPD on g s ). This relationship offers a potential mechanistic explanation for the observations of decoupling between canopy gas exchange and tree ring isotope composition (e.g., Offermann et al. 2011). Second, our current analysis highlights the importance of explicit consideration of leaf temperature in the 18 O c simulation. In earlier studies, T l was assumed equal to T a (Brooks and Mitchell 2011), incorporated into the study using observations (Song et al. 2013), simulated empirically (Cernusak et al. 2009, Munksgaard et al. 2016 or given hypothetical values (West et al. 2008). But the effect of T l on 18 O c has not been discussed. Similar to g s , the sensitivity of 18 O c to T l varies with RH, especially at low g s . Thus, ignoring fluctuation of T l in high RH environments may bias estimates of 18 O c . In this light, the new generation of thermal infrared cameras may simplify measurements of leaf temperature (Aubrecht et al. 2016), providing a more empirical means of assessing the leaf temperature effects that were modeled here. Finally, as source water may be different between drought treatment and the control, dual-isotope approach could not be applied and the multi-dimensional sensitivity analysis provides an alternative tool. Such quantitative analysis of 18 O between different environmental drivers as well as their interactions is novel and directly relevant for ecosystem models.
As the parameter space in this study is specific to our site conditions, we expect that multi-dimensional sensitivity analysis would expand the inference space of the 18 O models for plant ecophysiological studies when applied in other scenarios.

Supplementary data
Supplementary data for this article are available at Tree Physiology Online.