Can the Responses of Photosynthesis and Stomatal Conductance to Water and Nitrogen Stress Combinations Be Modeled Using a Single Set of Parameters?

Accurately predicting photosynthesis in response to water and nitrogen stress is the first step toward predicting crop growth, yield and many quality traits under fluctuating environmental conditions. While mechanistic models are capable of predicting photosynthesis under fluctuating environmental conditions, simplifying the parameterization procedure is important toward a wide range of model applications. In this study, the biochemical photosynthesis model of Farquhar, von Caemmerer and Berry (the FvCB model) and the stomatal conductance model of Ball, Woodrow and Berry which was revised by Leuning and Yin (the BWB-Leuning-Yin model) were parameterized for Lilium (L. auratum × speciosum “Sorbonne”) grown under different water and nitrogen conditions. Linear relationships were found between biochemical parameters of the FvCB model and leaf nitrogen content per unit leaf area (Na), and between mesophyll conductance and Na under different water and nitrogen conditions. By incorporating these Na-dependent linear relationships, the FvCB model was able to predict the net photosynthetic rate (An) in response to all water and nitrogen conditions. In contrast, stomatal conductance (gs) can be accurately predicted if parameters in the BWB-Leuning-Yin model were adjusted specifically to water conditions; otherwise gs was underestimated by 9% under well-watered conditions and was overestimated by 13% under water-deficit conditions. However, the 13% overestimation of gs under water-deficit conditions led to only 9% overestimation of An by the coupled FvCB and BWB-Leuning-Yin model whereas the 9% underestimation of gs under well-watered conditions affected little the prediction of An. Our results indicate that to accurately predict An and gs under different water and nitrogen conditions, only a few parameters in the BWB-Leuning-Yin model need to be adjusted according to water conditions whereas all other parameters are either conservative or can be adjusted according to their linear relationships with Na. Our study exemplifies a simplified procedure of parameterizing the coupled FvCB and gs model that is widely used for various modeling purposes.


INTRODUCTION
In the past decades, many crop models have been developed for predicting yield in response to changing environments. Some studies evaluated the performance of different crop models under different growth conditions such as different temperature, water supply and soil fertility (Jamieson et al., 1998;Adam et al., 2011;Palosuo et al., 2011). Surprisingly, when testing these models under a large land scale or long time span, the yield predictions in most models turned out to be an artifact of the balance between incorrect predictions of assimilation and leaf area index (Jamieson et al., 1998) or between biomass production and harvest index (Palosuo et al., 2011). The radiation-use efficiency approach that was taken in many crop models may over-simplify underlying processes and a more detailed approach, based on quantitative functional relationships for underlying processes, is needed in order to capture the effects of high temperature and high radiation intensities on crop growth under changing environments (Challinor et al., 2009;Adam et al., 2011). While detailed models usually require more effort in terms of model parameterization, some parameters and functional relationships are found to change very little (i.e., are conservative) among crop types (von Caemmerer et al., 2009) and environmental conditions (Yin, 2013). Therefore, it is important to test the conservative level of commonly used functional relationships, so as to balance between the level of detail in these models and the efforts needed for model parameterization.
Photosynthesis is the primary physiological process that drives crop growth and productivity and influences many plant quality traits, and is strongly affected by environmental factors. Accurately predicting photosynthesis is the first step toward predicting crop growth, yield and quality in response to environmental changes. Water and nitrogen variations frequently occur in crop fields. The effects of water and nitrogen on photosynthesis have been extensively and separately studied (Grassi et al., 2002;Xu and Baldocchi, 2003;Gu et al., 2012). The combined effect of water and nitrogen on photosynthesis, however, has received less attention.
Previous modeling studies have shown that the use of empirical factors to capture the effect of stresses, does not model photosynthesis reliably in many cases (Jamieson et al., 1998). The effects of environmental factors on leaf photosynthesis can be best investigated by use of the biochemical model of Farquhar, von Caemmerer and Berry (the FvCB model hereafter) (Farquhar et al., 1980) combined with diffusion models. The FvCB model has been widely used to describe photosynthesis in response to multiple environmental changes (Harley et al., 1992;Grassi et al., 2002;Xu and Baldocchi, 2003;Monti, 2006;Qian et al., 2012). The model describes photosynthesis as the minimum of the Rubisco-limited rate and the electron transportlimited rate. Major parameters in this model are the maximum Rubisco carboxylation rate (V cmax , definitions of all model variables hereafter are listed in Table 1), the maximum electron transport rate (J max ) and the mitochondrial day respiration (R d ). These biochemical parameters have been found to be linearly correlated with leaf nitrogen content per unit leaf area (N a ) under environmental changes such as various nitrogen supply (Grassi et al., 2002;Yin et al., 2009) and elevated CO 2 (Harley et al., 1992;Yin, 2013), as well as seasonal changes (Zhu et al., 2011). However, whether or not the linear relationships between these biochemical parameters and N a exist under drought is debatable, mainly due to inconsistent effect of drought on N a (Díaz-Espejo et al., 2006;Damour et al., 2008Damour et al., , 2009).
The FvCB model itself requires the CO 2 concentration in the chloroplast (C c ) as an input variable. To this end, estimating stomatal conductance (g s ) and mesophyll conductance (g m ) is necessary to enable the FvCB model to predict photosynthesis using the atmospheric CO 2 level (C a ) as input. The stomatal conductance model of Ball et al. (1987) (the BWB-type model hereafter), as one of the most commonly used models of g s , is often coupled with the FvCB model (Harley et al., 1992;Kosugi et al., 2003). In the BWB-type model, g s responds to net photosynthetic rate, relative humidity and CO 2 concentration at the leaf surface. Although it is phenomenological, the BWBtype model is widely used to model g s at leaf level (e.g., Leuning, 1995) and is the most feasible yet biologically robust tool for extrapolating g s at the field or forest stand level (Misson et al., 2002;Alton et al., 2007). The original BWBtype model does not capture stomatal responses to soil water status, thus some efforts were made toward modifying the BWB-type model to predict g s under drought. Either the slope used in the BWB-type model (describing the response of g s to photosynthetic rate, relative humidity or vapor pressure deficit (VPD) and CO 2 concentration) (Tuzet et al., 2003;Maseyk et al., 2008;Héroult et al., 2013) or the residual stomatal conductance (the value of g s when irradiance approaches to zero) (Misson et al., 2004) was reported to decrease under drought, and was related to soil moisture or leaf water potential (Baldocchi, 1997;Wang and Leuning, 1998;Misson et al., 2004;Keenan et al., 2010;Egea et al., 2011;Li et al., 2012;Zhou et al., 2013;Müller et al., 2014). In another study, however, neither of these two parameters was affected by drought (Xu and Baldocchi, 2003). So far, there is no consensus as to how to adjust the BWB-type model parameters to properly model g s under drought. Moreover, there are very few studies that investigated the responses of these parameters to nitrogen supply and to the combination of water and nitrogen supply.
g m has been considered as infinite in most early studies, in which intercellular CO 2 concentration (C i ) was used to substitute C c in the FvCB model (Harley et al., 1992;Kosugi et al., 2003). However, this assumption has later been proved not true since C c is lower than C i (Warren, 2004). Ignoring g m leads to the underestimation of V cmax , especially under stress conditions such as drought (Monti, 2006). g m has been found to decrease under water-deficit conditions and low nitrogen availability in many previous studies (reviewed in Flexas et al., 2008). There have been only a few attempts to incorporate the effect of drought on g m in the photosynthesis model by using a dependence of g m on g s (Cai et al., 2008) based on the observation of a close correlation between g s and g m in response to water-deficit conditions (Flexas et al., 2002;Warren, 2008;Perez-Martin et al., 2009) or by including an empirical soil moisture dependent function for g m (Keenan et al., 2010). Given that so far no consensus exists, more Factor used to calculate electron transport rate from chlorophyll fluorescence - Vapor pressure deficit kPa χ J Slope of the linear relationship between J max25 and N a µmol e − (g N) −1 s −1 χ V Slope of the linear relationship between V cmax25 and N a µmol CO 2 (g N) −1 s −1 φ 2 Apparent operating efficiency of PSII photochemistry mol e − (mol photon) −1 investigations are needed to incorporate the responses of g m to water and nitrogen variations into the photosynthesis model. When applying the combined FvCB, g s and g m model for predicting photosynthetic responses to fluctuating environmental variables, inevitably many parameters need to be quantified. Information about which parameters are conservative and which are variable depending on the treatment is extremely useful for predicting photosynthesis under diverse environmental conditions. Given the previous experience that the FvCB model parameters, once expressed as a function of N a , are not altered by environmental variables such as elevated [CO 2 ] (Yin, 2013), we are particularly interested in examining whether the responses of FvCB, g s and g m model parameters to water and nitrogen stress can be modeled using a single set of parameters when they are related to leaf nitrogen content. The objectives of this study are (i) to test whether or not water and nitrogen stress combinations change the linear relationships between photosynthetic biochemical parameters and leaf nitrogen content, and (ii) to investigate the responses of stomatal conductance model parameters and mesophyll conductance to different water and nitrogen conditions and to quantify these responses for the purpose of model simplicity. To this end, we used Lilium (L. auratum × speciosum "Sorbonne") as the test plant, as this plant is commonly grown under low-investment greenhouses where plants are frequently subject to different water and nitrogen regimes.

Plant Materials and Experimental Design
Four experiments with the same type of water and nitrogen treatments were conducted in different growth seasons in a plastic greenhouse located at Nanjing, China (32 • N, 118 • E) during 2009 to 2011 ( Table 2). The greenhouse, covered by antidrop polyvinyl chloride film, was composed of two spans and east-west oriented with a length of 28 m, span width of 8 m, gutter height of 3 m and arch height of 5 m. Heating pipes were installed during winter season. During summer season, the greenhouse was cooled through natural ventilation and an inner shading screen installed at the position with a distance of 1.0-1.4 m to the top. Temperature, VPD and photosynthetically active radiation are shown in the Supplementary Data (Figures S1-S3). No CO 2 enrichment was applied, and standard cultivation practices for disease and pest control were used as is common for commercial Lilium production in China. Lilium bulbs, with a circumference of 14-16 cm, were planted in plastic pots filled with substrates of sand, turf and soil (3:1:1). The physicochemical properties of the substrate are shown in Table 2. The pots, with a depth of 14 cm, upper diameter of 18 cm and bottom diameter of 12 cm, were put on seedling beds (l × w × h = 25.0 m × 1.7 m × 1.0 m) and arranged at a density of 36 plants m −2 .
Two water levels were used: well-watered conditions, with a soil water potential (SWP) of −4 to −15 kPa according to Li et al. (2012), and water-deficit conditions, with a SWP of −20 to −40 kPa. The SWP at 0.1 m below the soil surface was monitored using tensiometers (SWP-100, Institute of Soil Science, Chinese Academy of Sciences) with three replicates per water level. When the SWP reached its designed lower limit value, plants were irrigated until it reached the designed upper limit value. The SWP at 0.1 m below the soil surface and the corresponding gravimetric soil water content were measured to establish calibration curves. These curves were then used to determine the amount of water required for irrigation. The dates of starting water treatment in the four experiments are shown in Table 2.
At each water level, there were four levels of nitrogen supply: 25, 45, 65, and 85 mg available nitrogen per kg substrate (hereafter N25, N45, N65, and N85, respectively). Nitrogen was added in the substrate as urea taking into account that urea can be converted into nitrate within 1 or 2 days (Harper, 1984). The amount of urea needed was calculated based on the targeted treatment level and the amount of available nitrogen in the substrate ( Table 2), and urea was directly spread in the substrate, with the dates shown in Table 2. According to Sun (2013), 65 mg available nitrogen per kg substrate is the optimal level of nitrogen supply in commercial Lilium production for the cultivar used in this study. Treatments, with a plot area of 2.0 ×1.5 m 2 and three replicates per treatment, were arranged in a split-plot design with water level assigned to the main plots and nitrogen level to the sub-plots.

Gas Exchange and Chlorophyll Fluorescence Measurements
Gas exchange was measured on newly fully expanded leaf (the 4th leaf counting from the top downward) at flower bud visible stage using the LI-6400 Portable Photosynthesis System (Li-Cor BioScience, Lincoln, NE, USA) under 21% O 2 . In Experiment 1, both light response curves and C i response curves were measured in order to identify any differences in photosynthesis parameter estimation by using these two types of curves. For light response curves, incident irradiance (I inc ) in the leaf cuvette was decreased in the series of 1, 500,1,200,1,000,600,400,200,100,50,20, and 0 µmol m −2 s −1 , while keeping C a at 370 µmol mol −1 . For C i response curves, C a was increased stepwise: 50,100,150,200,250,380,650, 1,000, and 1,500, while keeping I inc at 800 µmol m −2 s −1 . The microclimate conditions in the leaf chamber were automatically controlled. The CO 2 concentration and water vapor between leaf and the reference chamber were automatically matched before data were recorded. We found that photosynthesis parameters estimated from A n -I inc curves and A n -C i curves were similar (see Results). Therefore, in Experiments 2-4, only A n -I inc curves were measured, as measurement of A n -C i curves inevitably involves the problem of CO 2 leakage into and out of the leaf cuvette, which would require additional measurements to correct for.
Chlorophyll fluorescence was simultaneously measured using FMS2 (Hansatech Instruments Ltd, UK) at a similar position on the leaf where gas exchange was measured. The steady-state fluorescence (F s ) was measured under natural radiation level (ranged from 0 to 1,200 µmol m −2 s −1 ) and saturating I inc (at 1,500 µmol m −2 s −1 ) after 3-5 min light adaptation, followed by applying a light pulse > 7,000 µmol m −2 s −1 for <1 s to measure maximum fluorescence F  photosystem II photochemistry (Φ 2 ) was calculated as Genty et al., 1989). Due to inadequate environmental control in the lowinvestment greenhouse, air temperature and VPD hardly stayed constant although they were kept within the range suitable for Lilium growth (Figures S1, S2). Therefore, all gas exchange and chlorophyll fluorescence measurements in the four experiments were subjected to variations of temperature and VPD.
In order to convert chlorophyll fluorescence data on Φ 2 into electron transport rate, combined measurement of gas exchange and chlorophyll fluorescence was conducted using the LI-6400XT Portable Photosynthesis System (Li-Cor BioScience, Lincoln, NE, USA) at low oxygen using a gas blend of 2% O 2 and 98% N 2 in the leaf chamber at flower bud visible stage (Experiment 1). A n −I inc curves were measured while keeping C a at 1,000 µmol mol −1 , to create non-photorespiratory conditions. A n at high C a levels (i.e., 650, 1,000, and 1,500 µmol mol −1 ) at 2% O 2 was also measured while keeping I inc at 800 µmol m −2 s −1 . Φ 2 was assessed using the same procedure as described above. In order to establish the correlation of estimating R d using different methods, combined measurement of gas exchange and chlorophyll fluorescence for A n −I inc curves (under 21% O 2 , keeping C a at 370 µmol mol −1 ) was also conducted in Experiment 1. All gas exchange data wherever the set-point C a differed from the ambient CO 2 level were corrected for CO 2 leakage from measurements using thermally killed leaves.

Leaf Characteristics
After gas exchange and chlorophyll fluorescence measurements, the leaves were cut, and leaf area was measured before being put in the oven at 105 • C for 30 min and subsequently at 80 • C until constant weight. Leaf nitrogen concentration (for organic nitrogen) was measured by using the Kjeldahl digestion method (Sun, 2013). Briefly, leaf dry samples were ground, and a 0.5 g of ground sample was digested with 30% hydrogen peroxide and 5 mL of concentrated sulphuric acid at 340 • C. 10 mL of 10 mol L −1 sodium hydroxide was then added for distilling the digested solution. The distillate was titrated using 0.02 mol L −1 sulfuric acid, and bromocresol green-methyl red was used as the indicator. Leaf nitrogen content per unit leaf area (N a , g m −2) was calculated based on leaf nitrogen concentration, leaf dry weight and leaf area.

Estimation of Photosynthetic Model Parameters
The FvCB model (Farquhar et al., 1980) predicts net photosynthetic rate (A n ) as the minimum of the Rubisco carboxylation-limited rate (A c ) and the electron transport-limited rate (A j ): where C c and O are the chloroplast partial pressures of CO 2 and O 2 , respectively; K mC and K mO are the Michaelis-Menten coefficients of Rubisco for CO 2 and O 2 , respectively; R d is day respiration; Γ * is the CO 2 compensation point in the absence of R d and was calculated as 0.5O K mC K mO [exp(−3.3801 + 5,220 298R(T+273) )] (Yin et al., 2004), derived from the parameter values of Bernacchi et al. (2001); J is the photosystem II electron transport rate that is used for CO 2 fixation and photorespiration.
R d was firstly estimated as the y-axis intercepts of the linear regression plots of A n against I inc (the Kok method hereafter) (Sharp et al., 1984). The Kok method tends to underestimate R d (Sharp et al., 1984;Yin et al., 2011). Therefore, R d was also estimated from the linear regression of A n against (I inc 2 /4) (the Yin method hereafter) (Yin et al., 2009(Yin et al., , 2011 using data available from the combined measurement of gas exchange and chlorophyll fluorescence, in order to establish the calibration relationship between values of R d estimated by the two methods. As combined gas exchange and chlorophyll fluorescence data were used only in part of our measurements, all the R d estimated based on the Kok method was then corrected according to the established calibration relationship to obtain R d estimates for all treatments. The calculation of A c or A j in the FvCB model requires C c , which is unknown beforehand. Therefore, A j relevant parameters were estimated based on Yin et al. (2009) using chlorophyll fluorescence data. To convert fluorescence-based data on Φ 2 into electron transport rate J, a calibration needs to be made for each water and nitrogen treatment. This was done by linear regression plot of A j against (I inc Φ 2 /4), using data obtained under nonphotorespiratory conditions from low light levels of the A n −I inc curve and three high CO 2 levels. The slope s of this linear regression was used as a calibration factor to calculate values of electron transport rate under all conditions: J = sI inc 2 (Yin et al., 2009). The obtained J was then fitted to the following equation to obtain electron transport parameters of the FvCB model: where κ 2LL is the conversion efficiency of incident light into J at strictly limiting light; J max is the asymptotic maximum value of J when I inc approaches to saturating level; θ is a convexity factor for response of J to I inc , and was assumed to have a constant value of 0.8 (Yin and Struik, 2015). Since chlorophyll fluorescence measurement was conducted under fluctuating temperature, the value of J max at 25 • C (J max25 ) and κ 2LL were calculated by combining Equation (4) with Equation (7) (see later) that describes the temperature response of J max . With J max25 and κ 2LL calculated as described above, J max for each A n −I inc curve from gas exchange measurement was derived according to the temperature level during each measurement using Equation (7) (see later). J at each light level in the A n −I inc curve was then derived using Equation (4) based on J max and κ 2LL calculated before.
With J and R d calculated, g m was then estimated assuming that g m was constant across the entire light response curve. Whether or not g m is constant across light or CO 2 levels remains debatable, but this assumption allows the identification of any differences among water and nitrogen treatments in the actual average g m . For that purpose, a relatively less measurement error-sensitive method, the NRH-A method (Yin and Struik, 2009a), was used to estimate the value of g m as constant, by fitting the following non-rectangular hyperbolic (NRH) equation for the A j part of the C i -based FvCB model: where x 1 = J/4 and x 2 = 2Γ * ; C i is the intercellular CO 2 level. According to our experimental data, A j -limitation in a light response curve of Lilium usually occurred at or below 1,000 µmol m −2 s −1 , as a good linear relationship between A n and J was observed within this range ( Figure S4). The advantages of the NRH-A method over other existing methods including the most widely used variable-J method in deriving the average g m was fully illustrated by Yin and Struik (2009a). Equation (5) can also be applied to calculate A c by setting: x 1 = V cmax and x 2 = K mC (1 + O/K mO ). V cmax was then estimated by fitting the combined Eqs. (1), (4) and (5) to the entire light response curve or C i response curve using the already estimated values of J max , κ 2LL , R d and g m as input.

Temperature Responses of Photosynthesis Parameters
To account for the effect of the varying temperature during measurement, temperature response functions were introduced so that the estimation of key parameters could be adjusted to the same reference temperature for the comparison among treatments. The temperature responses of R d and Rubisco kinetic properties (V cmax , κ mC and κ mO ) were described by an Arrhenius function Equation (6), and the temperature responses of J max and g m were described by a peaked Arrhenius function Equation (7), normalized with respect to their values at 25 C: where X stands for each parameter; X 25 is the value of each parameter at 25 • C (R d25 , V cmax25 , κ mC25 , κ mO25 , J max25 , and g m25 ); E x is the activation energy of each parameter (E Rd , E Vcmax , E KmC , E KmO , E Jmax , and E gm ); S x and D x are the entropy term and the deactivation energy, respectively (applying to J max and g m ); T is the leaf temperature; R is the universal gas constant. Since Rubisco kinetic properties are generally assumed conserved among C 3 species (von Caemmerer et al., 2009), values of κ mC25 , κ mO25 , E KmC , and E KmO were fixed at 272.4 µbar, 165.8 mbar, 80,990 J mol −1 , and 23,720 J mol −1 , respectively, according to Bernacchi et al. (2002). To avoid over-parameterization, E Rd was fixed at 46,390 J mol −1 (Bernacchi et al., 2001); S Jmax and D Jmax were fixed at 650 J K −1 mol −1 (Harley et al., 1992) and 200,000 J mol −1 (Medlyn et al., 2002), respectively; E gm , S gm , and D gm were fixed at 49,600 J mol −1 , 1,400 J K −1 mol −1 , and 437,400 J mol −1 , respectively (Bernacchi et al., 2002).

The Relationships between Biochemical Parameters and Leaf Nitrogen Content
The photosynthetic capacity parameters V cmax25 and J max25 are linearly related to N a (Harley et al., 1992;Braune et al., 2009): where N b is the base leaf nitrogen content at or below which A n is zero, and a value of 0.35 g N (m 2 leaf) −1 was used in this study (Archontoulis et al., 2012); χ V is the slope of V cmax25 against N a , and χ J is the slope of J max25 against N a .

Parameterization of the Stomatal Conductance Model
A phenomenological model for stomatal conductance for CO 2 transfer was first described by Ball et al. (1987), revised by Leuning (1995), and further revised by Yin and Struik (2009b). Li et al. (2012) called this model the BWB-Leuning-Yin model. In the model, stomatal conductance was described by: where g 0 is the residual stomatal conductance when the irradiance approaches to zero; C i * is the C i -based CO 2 compensation point in the absence of R d and was calculated as (Γ * − R d /g m ) using Γ * , R d and g m calculated before as input; f vpd is a function describing the effect of VPD, which is not yet understood sufficiently and may be described empirically as Yin and Struik (2009b): where a 1 represents the ratio of C i to C a for vapor saturated air, and b 1 represents the decreasing slope of this ratio with increasing VPD, if g 0 approaches to zero. Because of this obvious meaning of a 1 and b 1 , we chose Equation (11), instead of the equation of Leuning (1995), for our analysis of the effect of VPD on g s . Combining Equations (10) and (11), g 0 , a 1 and b 1 can be estimated by using the data of A n , C i and VPD obtained from gas exchange measurement. For that, measured stomatal conductance for water vapor transfer was divided by a factor 1.6 to convert it to g s for CO 2 transfer that is required for Equation (10).

Statistical and Model Analyses
Using a non-linear regression with the GAUSS method in PROC NLIN of SAS (SAS Institute Inc., Cary, NC, USA), FvCB model parameters (V cmax25 , J max25 , κ 2LL , R d25 , g m25 , E Vcmax , and E Jmax ) and BWB-Leuning-Yin model parameters (g 0 , a 1 , and b 1 ) were estimated. Whether or not the treatment effect on each estimated parameter was significant was tested using an F-test. Following that, conserved parameter values across treatment classes were also estimated. With these estimated parameters available, we aimed to test to what extent conserved parameter values could be used to predict A n and g s under water and nitrogen stress combinations, for the purpose of simplifying model parameterization. For such, a step-wise procedure was followed. First, we analyzed whether or not water and nitrogen stress combinations change the linear relationships between biochemical parameters and N a , and tested to what extent conserved parameter values in the C i -based FvCB model (Equation 5) could be used to predict A n under different water and nitrogen conditions. Second, we tested to what extent conserved parameter values could be used in the BWB-Leuning-Yin model to predict g s under different water and nitrogen conditions. Third, we explored the coupled FvCB and BWB-Leuning-Yin model (for the analytical solution for this coupled model, see Yin and Struik, 2009b), which allows using C a as input to predict A n . We used this coupled model to assess to what extent conserved parameter values in both the FvCB model and the BWB-Leuning-Yin model could be used to predict A n (using C a as input) across various water and nitrogen treatment regimes.

Model Parameterization
Data of A n −I inc curves showed that both water-deficit conditions and low nitrogen supply decreased A n (Figure 1). The initial linear part of these curves was explored to estimate R d . Values of R d estimated by the Kok method were generally lower than those estimated by the Yin method (Figure 2). The linear correlation between values of R d estimated by the two methods (Figure 2) was used to correct all R d estimated by the Kok method.
The plot of A j against (I inc Φ 2 /4) using data obtained under low O 2 condition from low light levels of the A n −I inc curves and three high CO 2 levels was essentially linear (Figure 3). Both water and nitrogen conditions affected the value of the linear slope s, the calibration factor used to convert Φ 2 into J. The factor decreased by low nitrogen supply and by water-deficit conditions.
V cmax estimated from A n -I inc curves and from available A n -C i curves under the same measurement conditions were very similar when using the same input values of J max , κ 2LL , R d , and g m (Figure 4). This suggested the reliability of using A n -I inc curves to estimate V cmax . The estimated parameter values of the FvCB model for each treatment are listed in Table 3, and those of the BWB-Leuning-Yin model and g m are listed in Table 4. All parameters were reliably estimated, as the standard error values of the estimates were relatively small (Tables 3, 4).

The Response of Estimated Parameter Values to Water and Nitrogen Treatments
Water-deficit conditions significantly decreased V cmax25 , J max25 , k 2LL , and R d25 at all nitrogen levels ( Table 3). V cmax25 , J max25 , and R d25 decreased with decreasing of nitrogen availability whereas κ 2LL showed such a response to a much less clear extent under both water-deficit conditions and well-watered conditions ( Table 3). V cmax25 , J max25 , and κ 2LL were significantly lower in the combined water deficit and low nitrogen availability treatments than in other treatments ( Table 3). Neither E Jmax nor E Vcmax was significantly affected by water and nitrogen treatments ( Table S1).
Water-deficit conditions significantly decreased g 0 , a 1 , b 1 , and g m25 at all nitrogen levels ( Table 4). g 0 and g m25 decreased with decreasing nitrogen availability whereas a 1 and b 1 responded little to nitrogen treatments under both waterdeficit conditions and well-watered conditions ( Table 4). g m25 was significantly lower under the combined water deficit and the lowest nitrogen availability treatment than in other treatments ( Table 4).

The Relationships between Estimated Parameter Values and Leaf Nitrogen Content
Under both water-deficit conditions and well-watered conditions, V cmax25 , J max25 , κ 2LL , R d25 , g m25 , and g 0 linearly increased with increasing N a (Figure 5). X V and X J were determined as 62 µmol (g N) −1 s −1 and 93 µmol (g N) −1 s −1 , respectively (Figures 5A,C). The N adependent relationship was relatively less clear for other parameters (Figures 5B,D-F), but an F-test revealed that water and nitrogen treatments did not significantly alter the linear relationships in all the six parameters. Linear relationship existed between V cmax25 and J max25 with a slope of 1.49 under different water and nitrogen treatments (Figure 6).

Comparison between Model Predictions and Measured Values for A n and g s
Since the linear relationships between biochemical parameters and N a were found to exist under different treatment combinations (Figure 5), we further tested to what extent conserved parameter values could be used in the FvCB model to predict A n under different water and nitrogen conditions. Two sets of comparisons between the measured A n and the predicted A n were conducted, (i) using treatment-specific parameter values (i.e., using specific parameter values obtained under each treatment) (Figures 7A,B), and (ii) using shared parameter values (i.e., incorporating the N a -dependent linear relationships and using overall E Jmax and E Vcmax values) (Figures 7C,D). For this second set of comparison, the overall values of E Jmax and E Vcmax for all treatments were estimated (Table S1) by incorporating the linear relationships between parameters (V cmax25 , J max25 , κ 2LL , R d25 , and g m25 ) and N a . The coefficient of determination (r 2 ) between estimated and measured A n in both comparisons ranged from 0.85 to 0.94 (Figure 7).
We also tested to what extent conserved parameter values could be used in the BWB-Leuning-Yin model (Equations 10 and 11) to predict g s under different water and nitrogen conditions. Since nitrogen had been found to have little effect on a 1 and b 1 ( Table 4) and g 0 could be linearly correlated with N a under both well-watered conditions and water-deficit conditions (Figure 5F), we tested to what extent conserved values of a 1 , b 1 , and g 0 can be used. For this purpose, we incorporated the linear relationships between model parameters (g 0 and g m25 ) and N a , and estimated overall values of a 1 and b 1 for all treatments (Table 4). Three sets of comparisons between the measured g s and the predicted g s were conducted, (i) using treatment-specific parameter values (Figures 8A,B), (ii) using shared parameter values for each water treatment (i.e., incorporating the N a -dependent linear relationships and using overall values of a 1 and b 1 for each water treatment group given in Table 4) (Figures 8C,D), and (iii) using shared parameter values for all treatments (i.e., incorporating the N adependent linear relationships and using overall values of a 1 and b 1 for all treatments given Table 4) (Figures 8E,F). Using treatment-specific parameter values in the BWB-Leuning-Yin model, the r 2 between estimated and measured g s was 0.61 under well-watered conditions and 0.57 under a water deficit (Figures 8A,B); using shared parameter values for each water treatment, the r 2 was 0.55 for well-watered plants and 0.43 under a water deficit (Figures 8C,D). When shared parameters were used for all treatments, g s was appreciably underestimated under FIGURE 4 | Comparison of V cmax estimated from A n -I inc curves and A n -C i curves using the same input values of J max , κ 2LL , R d and g m (Each point represents value of V cmax estimated from A n -I inc curve or A n -C i curve measured on the same leaf).
well-watered conditions (Figure 8E), but overestimated under a water deficit (Figure 8F).This third set of predictions of g s , when compared with the first set of predictions, underestimated g s by 9% under well-watered conditions and overestimated g s by 13% under water-deficit conditions.
As g s was either underestimated or overestimated by the BWB-Leuning-Yin model using shared parameter values for all treatments (Figures 8E,F), we further assessed the impact of this inaccurate estimation of g s on the prediction of A n . Two sets of comparisons between the measured A n and the predicted A n were conducted. In the first comparison, shared values of the FvCB model parameters for all treatments and shared values of the BWB-Leuning-Yin model parameters for each water treatment were used in the coupled model; the r 2 between estimated and measured A n was 0.89 under wellwatered conditions and 0.80 under water-deficit conditions (Figures 9A,B). In the second comparison, shared values of both the FvCB model parameters and the BWB-Leuning-Yin model parameters for all treatments were used; the r 2 was 0.89 under well-watered conditions (Figure 9C), but A n was overestimated by 9% under water-deficit conditions ( Figure 9D).

Methodology to Estimate Photosynthetic Parameters
In our study, all model parameters were estimated based on the A n -I inc curves, instead of A n -C i curves, for estimating the FvCB parameters. We tested that the estimated V cmax values by using these two types of curves were quite similar (Figure 4), as also shown in a previous study (Archontoulis et al., 2012). The approach of using A n -I inc curves provides an alternative to the prevailing approach of using A n -C i curves and has its own advantages. First, the FvCB model is commonly used to predict leaf photosynthesis in canopies under field conditions, where it is the light level, not the CO 2 level, that fluctuates most significantly in space and in time. This suggests that the FvCB parameters estimated from A n -I inc curves should more closely represent field situations, relative to those based on A n -C i curves. Second, using A n -C i curve is known to have problems of CO 2 leakage and down-regulation of Rubisco at the low level of CO 2 during the measurement. The A n -I inc curve-based approach avoids these problems since the whole response curve is measured under ambient CO 2 level. However, using A n -I inc curves also tends to have problems. First, V cmax cannot always be estimated from A n -I inc curves since the entire A n -I inc curve can be A j limited sometimes (Archontoulis et al., 2012), especially for field crops that have high light saturating point (e.g., cotton, Wise et al., 2004). Second, the rate of TPU (triose phosphate utilization), if exerting a limitation on photosynthesis, cannot be estimated using A n -I inc curves since like Rubisco limitation, any TPU limitation on A n -I inc curves also happens at high irradiance levels (Archontoulis et al., 2012). Nevertheless, our limited data (Figure 4) show the evidence in support of using A n -I inc curves as an alternative approach to estimate V cmax . More comparisons between the two approaches using A n -I inc and using A n -C i curves are needed for different crop types and environments. We adopted some parameter values from the literature as input to avoid over-parameterization of the FvCB model. First, θ (the convexity factor for response of electron transport rate to incident light) was set to a constant value of 0.8 according to Yin and Struik (2015). It is worthy to notice that the actual value of θ could vary across species and environments. In our experiment, θ may be affected by different water and nitrogen treatments, as well as different light environment caused by different growth season. Initial analyses showed that letting θ be fitted as well resulted in enormous unrealistic variation of the estimated J max and κ 2LL . Since the biological meaning of θ is less obvious than that of J max and κ 2LL , we decided to set θ as a constant value to avoid biased estimations of J max and κ 2LL . Equation (4) with θ of 0.8 generates a very similar light response shape as given by the other widely used quadratic hyperbolic equation initially used by Harley et al. (1992). Second, in line with some previous studies (Xu and Baldocchi, 2003;Li et al., 2012), we adopted the activation energy of R d (E Rd ) and g m (E gm ), the deactivation energy of J max (D Jmax ) and g m (D gm ), and the entropy term of J max (S Jmax ) and g m (S gm ) from literature (Bernacchi et al., 2001(Bernacchi et al., , 2002. Whether or not these temperature response parameters change with water and nitrogen conditions is still not clear and further studies are needed. Third, Rubisco kinetic properties (κ mC25 , κ mO25 , E KmC , and E KmO ) were adopted from Bernacchi et al. (2002). Despite the generally assumption that Rubisco kinetic properties are conserved among C 3 species (von Caemmerer et al., 2009), values of these constants reported in the literature are different (Bernacchi et al., 2001(Bernacchi et al., , 2002Dreyer et al., 2001). The choice of Rubisco parameters also affected our FvCB parameter estimation. Since all parameters in the FvCB model are interrelated with each other, potential errors in our parameter estimation exist  Treatment if parameter values we adopted from the literature were not applicable in our study.

Photosynthetic Biochemical Parameters in Response to Water and Nitrogen Conditions
Our study showed that a long-term mild water deficit and water and nitrogen stress combinations did not have significant effects on the linear relationships between biochemical parameters of the FvCB model (i.e., J max25 , κ 2LL , V cmax25 ) and leaf nitrogen content per unit area (N a ) (Figure 5). Previous studies showed that a short-term water deficit did not change the linear relationships between biochemical parameters and N a (Díaz-Espejo et al., 2006;Gu et al., 2012), whereas under longterm drought, either the slopes of the relationships between biochemical parameters and N a were changed (Wilson et al., 2000;Díaz-Espejo et al., 2006) or considering the effect of leaf mass per area (LMA) in the linear regressions was needed (Xu and Baldocchi, 2003). A few other studies (Damour et al., 2008(Damour et al., , 2009 found that drought totally modified the fundamental relationships between J max and N a since N a was either increasing (Damour et al., 2008) or not affected (Damour et al., 2009) under drought whereas J max decreased. The discrepancy of the response of N a to drought found in different studies may be caused by different species. Damour et al. (2008) worked with lychee tree and Damour et al. (2009) worked with mango tree, whereas we focused on herbaceous species Lilium and Gu et al. (2012) worked with rice. Besides, different approaches used to estimate FvCB parameters could also affect the results in different studies. First, as stated earlier, we used A n -I inc curves to parameterize the FvCB model. Whether or not the approach of using A n -I inc curves and the approach of using A n -C i curves yield similar results under drought still requires more comparisons. Second, early studies tend to ignore s (the calibration factor for converting fluorescence-based efficiency of photosystem II photochemistry Φ 2 into electron transport rate J) and g m (mesophyll conductance) during the estimation of biochemical parameters. This could lead to inaccurate estimation of biochemical parameters since both s (Figure 3) and g m ( Table 4; also reviewed in Flexas et al., 2008) decreased under drought.
The calibration factor s used to convert Φ 2 into J is actually a lumped physiological parameter (s = ρ 2 β[1-f pseudo /(1−f cyc )]) that includes the absorptance of light by leaf photosynthetic pigments (β), the proportion of absorbed light partitioned to photosystem (PS) II (ρ 2 ), and the fraction of electrons at PSI following the cyclic transport around PSI (f cyc ) and following the pseudocyclic transport (f pseudo ) (Yin et al., 2009;Yin and Struik, 2009a). s was found to decrease by low nitrogen supply in previous study (Yin et al., 2009), which is also found in our study (Figure 3). This decrease may be explained by the decreasing of β as a result of the decreased photosynthetic pigments in lownitrogen leaves (Evans and Terashima, 1987). Interestingly, we found that s was smaller under water-deficit conditions compared to that under well-watered conditions despite the similar N a (e.g., s in N65 under well-watered conditions compared with s in N85 under water-deficit conditions). It has been reported that drought did not change the partitioning of electrons between PSI and PSII (Genty et al., 1987). However, stomatal closure caused by drought results in the decreasing of CO 2 concentration in the leaf, and consequently the amount of electrons used for CO 2 fixation decreases (Cornic and Briantais, 1991). Excessive electrons need to be consumed by other sinks apart from CO 2 fixation by following pseudo-cyclic electron transport (Cornic and Briantais, 1991;Biehler and Fock, 1996), or electrons need to follow cyclic flow around PSI (Kohzuma et al., 2009). Our results for the decreased s under water-deficit conditions independent on N a suggest that drought induced an increase of f pseudo or f cyc or both in our experimental conditions. Associated with estimating the factor s, mitochondrial day respiration (R d ) was estimated. Water-deficit conditions did not affect R d in all N treatments, and there were non-significant effects of nitrogen on R d under both well-watered conditions and water-deficit conditions ( Table 3). Nevertheless, water-deficit conditions significantly decreased R d in N85 and N45 treatments and generally there was a trend showing that drought and decreasing of nitrogen level decreased R d (Table 3), as also revealed in some previous studies (González-Meler et al., 1997;Huang and Fu, 2000). Therefore, we established an N a -dependent relationship of R d ( Figure 5F) and applied this relationship to capture the changes of R d under different water and nitrogen conditions. The linear relationship between respiration rate and leaf nitrogen content was also found under different light conditions (Ryan, 1995) and growth locations (Reich et al., 1998).
A relatively stable J max25 /V cmax25 ratio among different water and nitrogen treatments was found in our study (Figure 6), in line with some previous studies (Makino et al., 1992;Walcroft et al., 1997;Díaz-Espejo et al., 2006). Some studies simplified the parameterization of the FvCB model by using a fixed value for either the J max /V cmax ratio (Kosugi et al., 2003) or the J max25 /V cmax25 ratio (Müller et al., 2005). However, care needs to be taken in setting a constant J max /V cmax ratio. First, when temperature varies, this ratio cannot be constant because J max and V cmax have different temperature response curves. In fact, the J max /V cmax ratio was found to decrease with temperature increase (Walcroft et al., 1997;Medlyn et al., 2002;Díaz-Espejo et al., 2006). When scaled to a common temperature, a better correlation between J max and V cmax was found (Leuning, 1997). Second, g m has a strong influence on this J max /V cmax ratio. In early studies (Grassi et al., 2002) when g m was not considered, a J max /V cmax ratio of ca 2.0 was obtained (Leuning, 1997), which is higher than our estimate where g m was considered (ca 1.5, Figure 6). Finally, some studies found that water and nitrogen conditions also affected the J max /V cmax ratio (Grassi et al., 2002;Gu et al., 2012). Therefore, the approach using a fixed value for the J max /V cmax ratio to parameterize the FvCB model should receive critical reservation (Xu and Baldocchi, 2003;Archontoulis et al., 2012).
In short, our study suggested that it is feasible to incorporate linear relationships between biochemical parameters and N a in the FvCB model to predict photosynthesis under different water and nitrogen conditions, since the FvCB model using shared parameter values for all treatments gave satisfactory predictions of A n under different water and nitrogen conditions (Figures 7C,D).

Stomatal Conductance Parameters and Mesophyll Conductance in Response to Water and Nitrogen Conditions
Accurately modeling stomatal conductance (g s ) and mesophyll conductance (g m ) are necessary steps toward predicting A n under changing environments. The BWB-type model of g s takes into account the effects of both environments and plant physiological status on g s , and has been widely tested able to satisfactorily predict g s for well-watered plants (Leuning, 1995;Li et al., 2012). Some efforts have been devoted to predict g s under drought conditions using the BWB-type model by introducing proper approaches to adjust parameter values used in the model. In general, most studies kept g 0 (residual stomatal conductance when the irradiance approaches to zero) as a fixed value and adjusted the value for the slope (roughly represents a 1 and b 1 in the BWB-Leuning-Yin model used in our study) by introducing a modifying factor of soil moisture (Egea et al., 2011;Li et al., 2012), or precipitation and evaporation (Baldocchi, 1997), or predawn xylem water potential (Sala and Tenhunen, 1996), or leaf nitrogen content and leaf water potential (Müller et al., 2014). Leuning (1995) suggested that the BWB-type model should be able to predict g s under water-deficit conditions by only adjusting the value for a 1 . We found that both a 1 and b 1 decreased with the decreasing of SWP (Table 4), and without considering these decreases, g s was overestimated under waterdeficit conditions ( Figure 8F). Further estimation of a 1 under water-deficit conditions by using the value for b 1 obtained under well-watered conditions resulted in a value of 0.586 for a 1 , which is much larger than the original value of 0.262 obtained under water-deficit conditions ( Table 4). Therefore, values for both a 1 and b 1 need to be adjusted to properly predict g s under waterdeficit conditions. However, a 1 and b 1 were little affected by nitrogen availability (Table 4) and no correlation between a 1 and N a , nor between b 1 and N a , under different water and nitrogen conditions was found in our study. The approach introducing a modifying factor of leaf nitrogen content on the slope (Müller et al., 2014) is able to predict g s in response to drought, and this could merely be due to similar responses of leaf nitrogen content and the slope to soil water condition rather than because a functional relationship exists between the slope and leaf nitrogen content. Our study did not present a quantitative relationship of a 1 and b 1 with water supply conditions since there were only two water-level treatments. Further studies including more water levels would be needed to quantify changes of a 1 and b 1 under different water and nitrogen conditions. g 0 was affected by both water conditions and nitrogen availability (Table 4), and a linear relationship between g 0 and N a (Figure 5F) was used in our study to take into account the changes of g 0 under different water and nitrogen conditions. Although this linear relationship is less clear compared to linear relationships between N a and biochemical parameters (e.g., J max25 and V cmax25 ) (Figure 5), an F test showed that there is no significant difference between using a conserved linear relationship and using separate relationships to describe the N a dependence of g 0 in response to water-deficit conditions. Under drought, plants tend to reserve water by reducing water loss, which makes it unlikely that g 0 is unaffected by waterdeficit conditions. However, few modeling studies considered the change of g 0 under drought condition (Misson et al., 2004;Keenan et al., 2010). The reason for using a fixed value for g 0 in previous studies could be that changing the value of g 0 should not affect the prediction of g s very much for plants with relatively high g s since the value of g 0 itself is normally very small and approaches to zero. However, this may not hold true for plants with low g s , as is the case in our study, since the value of g 0 may have relatively larger impact on predicting g s . g m has received growing attentions in modeling photosynthesis , since g m has been found to be finite and vary greatly among environments (Flexas et al., 2008;Yin et al., 2009). Previous studies found that g m decreased under drought and low nitrogen availability (reviewed in Flexas et al., 2008). We found that g m was enhanced by high nitrogen level and strongly decreased by the combination of water deficit and low nitrogen availability (Table 4). A relatively strong linear correlation between g m and N a was found in our study (Figure 5E), as also found in previous studies (von Caemmerer and Evans, 1991;Warren, 2004). Such a correlation may be explained by the surface area of the chloroplasts facing the cell walls, an anatomical determinant of g m (von Caemmerer and Evans, 1991;Evans et al., 1994), which depends on N a .
Our results showed that the relation between g m and N a was hardly changed by water-deficit conditions ( Figure 5E). In contrast, Gu et al. (2012) found that the change of g m by water-deficit conditions was not explained by the change of N a but was negatively correlated with LMA. Nevertheless, LMA is generally considered as setting a limitation for the maximum g m (Flexas et al., 2008;Perez-Martin et al., 2009) rather than is used to model g m in response to environments, mainly because the change of LMA results from the long-term environmental adaptation of the plants (Poorter et al., 2009) whereas g m can vary quickly in response to environmental changes (Flexas et al., 2006). This is supported by our result of using the N a -dependent linear relationship to take into account the effects of water and nitrogen on g m . Together with the incorporation of other N adependent relationships of biochemical parameters, the model yielded similar results of A n prediction compared to those using treatment-specific parameter values (Figure 7).
Some studies incorporated a dependence of g m on g s in the photosynthesis model (Cai et al., 2008) as a close correlation between g s and g m in response to soil water deficit was commonly observed (Flexas et al., 2002;Warren, 2008;Perez-Martin et al., 2009). An approach incorporating the dependence of g m on g s was shown to give better prediction of A n of different genotypes than the one incorporating the dependence of g m on leaf nitrogen (Ohsumi et al., 2007). However, the approach has been criticized as having no physiological justification  since g m and g s respond differently to other environmental factors such as VPD (Warren, 2008;Perez-Martin et al., 2009). As there is not yet sufficient physiological knowledge to reliably quantify the variability of g m , some studies merely used a modifying factor of soil water conditions to take into account the effect of water deficit on g m (Keenan et al., 2010;Egea et al., 2011). Whether or not the linear relationship between g m and N a could be a promising step toward modeling the variation of g m needs to be further tested.

The Effect of g s Estimation on the Prediction of A n
The coupled FvCB and BWB model has been increasingly used to model photosynthesis in response to environmental changes such as elevated CO 2 (Harley et al., 1992) and drought stress (Keenan et al., 2010;Müller et al., 2014) and seasonal changes (Kosugi et al., 2003). Normally in those previous studies, values of the biochemical parameters were related to the leaf nitrogen content and values of the stomatal conductance model parameters were changed according to the CO 2 level (Harley et al., 1992), leaf water potential (Müller et al., 2014), or growth season (Kosugi et al., 2003).
Our study showed that considering the decreases of the stomatal conductance model parameters (a 1 and b 1 ) by drought was needed, otherwise, the coupled FvCB and BWB-Leuning-Yin model overestimated A n under drought ( Figure 9D) due to an overestimation of g s (Figure 8F). The strong decrease of a 1 by drought (Table 4) indicates the decreasing of C i /C a ratio for vapor saturated air. The decrease of b 1 by drought (Table 4) suggests a negligible control of VPD on g s under drought condition. These results are in line with previous studies that under drought condition, g s at vapor nearly saturated air tended to be lower and g s was less sensitive to VPD (Forseth and Ehleringer, 1983;Perez-Martin et al., 2009). However, an exceptional case, which g s showed much stronger sensitivity to VPD under drought, was also found in the previous study without an explanation provided (Perez-Martin et al., 2009).
The BWB-Leuning-Yin model without considering the effect of water level on a 1 and b 1 also underestimated g s under wellwatered conditions ( Figure 8E). But the subsequent prediction of A n was not affected much ( Figure 9C). This is probably explained by that under well-watered conditions, C i is generally high and changing C i at its high level only slightly affects A n according to the diminishing-return relationship of A n vs. C i . Therefore, as shown in Figure 9, the estimation of g s had more effect on the prediction of A n under water-deficit conditions than under well-watered conditions.

CONCLUDING REMARKS
A previous analysis (Yin, 2013) showed that the relationship of many crop model parameters (including those FvCB biochemical parameters) as a function of plant nitrogen status was little altered by elevated CO 2 concentration. Our present study examined whether this assertion could be extended for the water and nitrogen stress combinations. We showed that the N a dependence of biochemical parameters of the FvCB model, g 0 of the BWB-Leuning-Yin model and the g m value were little altered by water and nitrogen stress combinations (Figure 5). By incorporating these N a -dependent relationships with the FvCB model and BWB-Leuning-Yin model, parameterization of these models could be simplified while maintaining satisfactory predictions. The obvious exception is parameters a 1 and b 1 of the BWB-Leuning-Yin model, which depended little on nitrogen treatments but greatly on water treatments ( Table 4). This is probably because the BWB-Leuning-Yin model is largely phenomenological, and its related conclusions are only valid for the specific species and conditions examined in this study. While the variation of parameters a 1 and b 1 had a great impact on the prediction of stomatal conductance, it had a considerably lower impact on the prediction of leaf photosynthesis. Nevertheless, a further study is needed to quantify how these two parameters vary with water-deficit conditions, as they have a stronger bearing on modeling leaf transpiration.

AUTHOR CONTRIBUTIONS
NZ analyzed the data and drafted the manuscript. GL, SY, DA, and QS carried out the measurements. WL and XY made substantial contributions to conception and experimental design, and critically revised the manuscript.

ACKNOWLEDGMENTS
This research was funded by China Natural Science Foundation (31171453), and was conducted in the framework of collaboration between the College of Agriculture, Nanjing Agricultural University and the Centre for Crop Systems Analysis, Wageningen University & Research. NZ thanks the China Scholarship Council for awarding her a fellowship to conduct this work. We thank Dr. Chuang Cai for the discussion about data analysis and the two reviewers for their critical comments.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2017. 00328/full#supplementary-material Table S1 | Activation energy of J max and V cmax (standard error of estimate in brackets) estimated for each water and nitrogen treatments and their shared values for all treatments. Different letters following the data in the same column indicate significant difference (P < 0.05).    . All data points were chosen from light levels at or below 1,000 µmol m −2 s −1 and leaf temperature at 20 ± 2 • C (Vertical error bar indicates standard error of measured A n ; horizontal error bar indicates standard error of calculated J).