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 anti-drop 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₂ 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⁻².

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² 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 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₂. 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⁻² s⁻¹, while keeping C_a at 370 μmol mol⁻¹. 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⁻² s⁻¹. The microclimate conditions in the leaf chamber were automatically controlled. The CO₂ 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₂ 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⁻² s⁻¹) and saturating I_inc (at 1,500 μmol m⁻² s⁻¹) after 3–5 min light adaptation, followed by applying a light pulse > 7,000 μmol m⁻² s⁻¹ for <1 s to measure maximum fluorescence Fm′. The apparent operating efficiency of photosystem II photochemistry (Φ₂) was calculated as Φ2=1-Fs/Fm′ (Genty et al., 1989).

Due to inadequate environmental control in the low-investment 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 Φ₂ 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₂ and 98% N₂ 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⁻¹, to create non-photorespiratory conditions. A_n at high C_a levels (i.e., 650, 1,000, and 1,500 μmol mol⁻¹) at 2% O₂ was also measured while keeping I_inc at 800 μmol m⁻² s⁻¹. Φ₂ 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₂, keeping C_a at 370 μmol mol⁻¹) was also conducted in Experiment 1. All gas exchange data wherever the set-point C_a differed from the ambient CO₂ level were corrected for CO₂ leakage from measurements using thermally killed leaves.

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⁻¹ sodium hydroxide was then added for distilling the digested solution. The distillate was titrated using 0.02 mol L⁻¹ sulfuric acid, and bromocresol green-methyl red was used as the indicator. Leaf nitrogen content per unit leaf area (N_a, g m⁻²⁾ was calculated based on leaf nitrogen concentration, leaf dry weight and leaf area.

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₂ and O₂, respectively; K_mC and K_mO are the Michaelis-Menten coefficients of Rubisco for CO₂ and O₂, respectively; R_d is day respiration; Γ_* is the CO₂ compensation point in the absence of R_d and was calculated as 0.5OKmCKmO[exp(-3.3801+5,220298R(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₂ 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Φ₂/4) (the Yin method hereafter) (Yin et al., 2009, 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 Φ₂ 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Φ₂/4), using data obtained under non-photorespiratory conditions from low light levels of the A_n−I_inc curve and three high CO₂ 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Φ₂ (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₂ 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₁ = J/4 and x₂ = 2Γ_*; C_i is the intercellular CO₂ 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⁻² s⁻¹, 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₁ = V_cmax and x₂ = 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.

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₂₅ 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₃ 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⁻¹, and 23,720 J mol⁻¹, respectively, according to Bernacchi et al. (2002). To avoid over-parameterization, E_Rd was fixed at 46,390 J mol⁻¹ (Bernacchi et al., 2001); S_Jmax and D_Jmax were fixed at 650 J K⁻¹ mol⁻¹ (Harley et al., 1992) and 200,000 J mol⁻¹ (Medlyn et al., 2002), respectively; E_gm, S_gm, and D_gm were fixed at 49,600 J mol⁻¹, 1,400 J K⁻¹ mol⁻¹, and 437,400 J mol⁻¹, respectively (Bernacchi et al., 2002).

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² leaf)⁻¹ 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.

A phenomenological model for stomatal conductance for CO₂ 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₀ is the residual stomatal conductance when the irradiance approaches to zero; C_i* is the C_i-based CO₂ 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₁ represents the ratio of C_i to C_a for vapor saturated air, and b₁ represents the decreasing slope of this ratio with increasing VPD, if g₀ approaches to zero. Because of this obvious meaning of a₁ and b₁, 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₀, a₁ and b₁ 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₂ transfer that is required for Equation (10).

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₀, a₁, and b₁) 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.

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 Φ₂/4) using data obtained under low O₂ condition from low light levels of the A_n−I_inc curves and three high CO₂ 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 Φ₂ 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).

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₀, a₁, b₁, and g_m25 at all nitrogen levels (Table 4). g₀ and g_m25 decreased with decreasing nitrogen availability whereas a₁ and b₁ responded little to nitrogen treatments under both water-deficit 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).

Under both water-deficit conditions and well-watered conditions, V_cmax25, J_max25, κ_2LL, R_d25, g_m25, and g₀ linearly increased with increasing N_a (Figure 5). X_V and X_J were determined as 62 μmol (g N)⁻¹ s⁻¹ and 93 μmol (g N)⁻¹ s⁻¹, respectively (Figures 5A,C). The N_a-dependent 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).

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²) 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₁ and b₁ (Table 4) and g₀ 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₁, b₁, and g₀ can be used. For this purpose, we incorporated the linear relationships between model parameters (g₀ and g_m25) and N_a, and estimated overall values of a₁ and b₁ 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₁ and b₁ 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_a-dependent linear relationships and using overall values of a₁ and b₁ for all treatments given Table 4) (Figures 8E,F). Using treatment-specific parameter values in the BWB-Leuning-Yin model, the r² 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² 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 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² between estimated and measured A_n was 0.89 under well-watered 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² was 0.89 under well-watered conditions (Figure 9C), but A_n was overestimated by 9% under water-deficit conditions (Figure 9D).

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₂ 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₂ leakage and down-regulation of Rubisco at the low level of CO₂ during the measurement. The A_n-I_inc curve-based approach avoids these problems since the whole response curve is measured under ambient CO₂ 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, 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₃ species (von Caemmerer et al., 2009), values of these constants reported in the literature are different (Bernacchi et al., 2001, 2002; Dreyer 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 if parameter values we adopted from the literature were not applicable in our study.

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 long-term 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, 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 Φ₂ 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 Φ₂ into J is actually a lumped physiological parameter (s = ρ₂β[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 (ρ₂), 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 low-nitrogen 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₂ concentration in the leaf, and consequently the amount of electrons used for CO₂ fixation decreases (Cornic and Briantais, 1991). Excessive electrons need to be consumed by other sinks apart from CO₂ 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).

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₀ (residual stomatal conductance when the irradiance approaches to zero) as a fixed value and adjusted the value for the slope (roughly represents a₁ and b₁ 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₁. We found that both a₁ and b₁ decreased with the decreasing of SWP (Table 4), and without considering these decreases, g_s was overestimated under water-deficit conditions (Figure 8F). Further estimation of a₁ under water-deficit conditions by using the value for b₁ obtained under well-watered conditions resulted in a value of 0.586 for a₁, which is much larger than the original value of 0.262 obtained under water-deficit conditions (Table 4). Therefore, values for both a₁ and b₁ need to be adjusted to properly predict g_s under water-deficit conditions. However, a₁ and b₁ were little affected by nitrogen availability (Table 4) and no correlation between a₁ and N_a, nor between b₁ 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₁ and b₁ 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₁ and b₁ under different water and nitrogen conditions.

g₀ was affected by both water conditions and nitrogen availability (Table 4), and a linear relationship between g₀ and N_a (Figure 5F) was used in our study to take into account the changes of g₀ 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₀ in response to water-deficit conditions. Under drought, plants tend to reserve water by reducing water loss, which makes it unlikely that g₀ is unaffected by water-deficit conditions. However, few modeling studies considered the change of g₀ under drought condition (Misson et al., 2004; Keenan et al., 2010). The reason for using a fixed value for g₀ in previous studies could be that changing the value of g₀ should not affect the prediction of g_s very much for plants with relatively high g_s since the value of g₀ 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₀ may have relatively larger impact on predicting g_s.

g_m has received growing attentions in modeling photosynthesis (Niinemets et al., 2009), 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_a-dependent 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 (Niinemets et al., 2009) 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 coupled FvCB and BWB model has been increasingly used to model photosynthesis in response to environmental changes such as elevated CO₂ (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₂ 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₁ and b₁) 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₁ by drought (Table 4) indicates the decreasing of C_i/C_a ratio for vapor saturated air. The decrease of b₁ 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₁ and b₁ also underestimated g_s under well-watered 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.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer DT and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.