Deriving C4 photosynthesis parameters by fitting intensive A/Ci curves

Measurements of photosynthetic assimilation rate as a function of intercellular CO2 (A/Ci curves) are widely used to estimate photosynthetic parameters for C3 species, yet few parameters have been reported for C4 plants, because of a lack of estimation methods. Here, we extend the framework of widely-used estimation methods for C3 plants to build estimation tools by exclusively fitting intensive A/Ci curves (6-8 more sampling points) for C4 using three versions of photosynthesis models with different assumptions about carbonic anhydrase processes and ATP distribution. We use simulation-analysis, out-of-sample tests, existing in vitro measurements and chlorophyll-fluorescence-measurements to validate the new estimation methods. Of the five/six photosynthetic parameters obtained, sensitivity analyses show that maximal-Rubisco-carboxylation-rate, electron-transport-rate, maximal-PEP-carboxylation-rate and carbonic-anhydrase were robust to variation in the input parameters, while day-respiration and mesophyll-conductance varied. Our method provides a way to estimate carbonic anhydrase activity, a new parameter, from A/Ci curves, yet also shows that models that do not explicitly consider carbonic anhydrase yield approximate results. The two photosynthesis models, differing in whether ATP could freely transport between RuBP and PEP regeneration processes yielded consistent results under high light, but they may diverge under low light intensities. Modeling results show selection for Rubisco of low specificity and high catalytic rate, low leakage of bundle sheath and high PEPC affinity, which may further increase C4 efficiency.

existing in vitro measurements and chlorophyll-fluorescence-measurements to validate the 23 new estimation methods. Of the five/six photosynthetic parameters obtained, sensitivity 24 analyses show that maximal-Rubisco-carboxylation-rate, electron-transport-rate, maximal-25 PEP-carboxylation-rate and carbonic-anhydrase were robust to variation in the input 26 parameters, while day-respiration and mesophyll-conductance varied. Our method provides 27 a way to estimate carbonic anhydrase activity, a new parameter, from A/C i curves, yet also 28 shows that models that do not explicitly consider carbonic anhydrase yield approximate 29 results. The two photosynthesis models, differing in whether ATP could freely transport 30 between RuBP and PEP regeneration processes yielded consistent results under high light, 31 but they may diverge under low light intensities. Modeling results show selection for 32 Δ S, entropy for temperature 50 dependence for parameters; φ PSII , quantum yield; γ*(25), the specificity of Rubisco at 25°C; 51 g bs , bundle sheath conductance for CO 2 ; g bso , bundle sheath conductance for O 2 ; g m , 52 mesophyll conductance for CO 2 ; I, light intensity; J max (25), maximum rate of electron 53 transport at 25°C; K c (25), Michaelis-Menten constant of Rubisco activity for CO 2 at 25°C; 54 K o (25), Michaelis-Menten constants of Rubisco activity for O 2 ; K p (25), Michaelis-Menten 55 constants of PEP carboxylation for CO 2 ; O bs , O 2 concentration in the bundle sheath cells; 56 Q 10 for K p , temperature sensitivity parameter for K p ; R, the molar gas constant; R d , daytime 57 respiration; R dbs , daytime respiration in bundle sheath cells; R dm , daytime respiration in 58 mesophyll cells; Rubisco, ribulose-1,5-bisphosphate carboxylase/oxygenase; RuBP, 59 ribulose-1,5-bisphosphate; T k , leaf absolute temperature; V c , velocity of Rubisco 60 8 determining transition points instead of assigning in advance, and checking for inadmissible 135 fits. Second, since both RuBP regeneration and PEP regeneration need ATP (Hatch, 1987), 136 we also examine two different assumptions about ATP distribution between RuBP 137 regeneration and PEP regeneration in C 4 photosynthesis models. Third, we validate the 138 estimation methods in four independent ways, using: (i) simulation tests using A/C i curves 139 generated using our model with known parameters and adding random errors, (ii) out of 140 sample test, (iii) existing in vitro measurements and (iv) Chlorophyll fluorescence 141 measurement. Finally, we used the C 4 photosynthesis model to perform sensitivity analyses 142 and simulation analyses for important physiological input parameters. These analyses allow 143 us to illustrate the underlying physiological significance of these parameters to the ecology 144 and evolution of the C 4 photosynthesis pathway. 145 146

C 4 Mechanism 148
The CO 2 concentrating mechanism of C 4 pathway increases CO 2 in the bundle sheath cells 149 to eliminate photorespiration. Like the C 3 pathway, the diffusion of CO 2 starts from the 150 ambient atmosphere through stomata into intercellular spaces, and then into the mesophyll 151 cells. In the mesophyll cells, the first step is the hydration of CO 2 into HCO 3 by carbonic 152 anhydrase. PEPC, then, catalyze HCO 3 and PEP into C 4 acids and the C 4 acids are 153 transported to the bundle sheath cells. In the bundle sheath cell, C 4 acids are decarboxylated 154 to create a high CO 2 environment for the C 3 photosynthetic cycle, and PEP is regenerated. 155 All the modeling equations and mechanistic processes used for our estimation method are 156 from von Caemmerer (2000), Hatch and Burnell (1990), Boyd et al. (2015)  RuBP regeneration and PEP carboxylation limited assimilation (AET) and RuBP 165 regeneration and PEP regeneration limited assimilation (ATT). Since the C 4 cycle operates 166 before the C 3 cycle and provides substrates for the C 3 cycle, the determination process of A n 167 is as follows: 168 (1) 169 (2) 170 which we used for our estimation method. 171 172

Plant Material 173
We performed intensive A/C i curves on nine different C 4 species to develop and examine the 174 efficacy of our estimation tools: Zea mays L., Eragrostis trichodes (Nutt.) Alph. Wood, 175 Andropogon virginicus L., Schizachyrium scoparium (Michx.) Nash, Panicum virgatum L., 176 Panicum amarum Elliott, Setaria faberi Herrm., Sorghastrum nutans (L.) Nash 177 andTripsacum dactyloides (L.) L. The intensive A/C i curves contain more sample points 178 under more CO 2 concentrations than the default curve used for C 3 species. Here we set the 179 CO 2 concentrations as 400, 200, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300, 325, 1 0 350, 400, 500, 600, 700, 800, 1000, 1200, 1400 ppm under light intensity of 1500 μ molm -2 s -181 1 (light intensity encountered by the plants in greenhouse). At each point, data was recorded 182 when the intercellular CO 2 concentration equilibrated within 2-5 minutes. The datasets were 183 obtained using a standard 2 x 3 cm 2 leaf chamber with a red/blue LED light source of LI-184 6400 (LI-COR Inc., Lincoln, NE, USA). If the stomatal conductance of a species does not 185 decrease quickly at high CO 2 , then the sample points at the high CO 2 level can be increased. 186 Fluorescence was measured along with A/C i curves for seven C 4 species (CO 2 concentration 187 is similar with above). After each change of CO 2 concentration and A reached steady state, 188 the quantum yield was measured by multiphase flash using a 2 cm 2 fluorescence chamber 189 head (Bellasio et al., 2014). All the measurements are conducted at 25 o C and VPD is 190 controlled at 1-1.7kPa. The cuvette was covered by Fun-Tak to avoid and correct for the 191 leakiness (Chi et al., 2013). 192 193

Estimation Protocol 194
We implemented the estimation methods using the non-linear curve-fitting routine in MS 195 Excel (Supplementary Material I, II, III) and independently in R ("C4Estimation") to get 196 solutions that minimize the squared difference between observed and predicted assimilation 197 rates (A). Five (or six when considering carbonic anhydrase) parameters will be estimated by 198 fitting the A/C i curve: V cmax , J, R d , V pmax , g m , and k CA. Other input parameters for C 4 are in 199 Table S1. 200 201 Input data sets and preliminary calculations. The input data sets are the leaf temperature 202 during measurements, atmosphere pressure, two CO 2 bounds (CaL and CaH discussed in the 1 1 following section), and the assimilation rates (A) and the C i s (in ppm) in the A/C i curve. 204 Also, reasonable initial values of output parameters need to be given in the output section to 205 initiate the non-linear curve fitting (Supplementary Material IV). C i will be adjusted from 206 the unit of ppm to the unit of Pa inside the program as suggested by Sharkey et al. (2007). 207 208 Estimating limitation states. We set upper and lower limits to the value of C i between 209 which the assimilation rates are freely determined by limitation states. Also, we can avoid 210 over-parameterization by pre-assigning limitation states at the lower and upper ends of the 211 C i range. We assumed that under very low C i (CaL), CO 2 is the limiting substrate; thus, V p is 212 limited by V pc and A is given by A c (AEE); under very high C i (CaH) electron transport is 213 limiting, thus, V p is limited by V pr and A is given by A j (ATT) (Fig. 1). The points between 214 CaL to CaH are freely determined by AEE, ATE, AET or ATT from eq. (16) and (17) to 215 minimize the cost function. We suggest setting CaL as 10 Pa initially, then adjusting based 216 on the preliminary results. The points of constant A at high C i end can initially be set as 217 being limited by ATT primarily (based on the three points, we can CaH) or use 65 Pa as the 218 first trial. The range of freely determined points can be adjusted by users by setting 219 appropriate CaL and CaH. In the column of "Estimate Limitation", whether the data points 220 are limited by AEE (represented by "1"), ATT (represented by "4") or freely vary 221 (represented by "0"), all the assignments of "1", "4" and "0" are determined automatically 222 by the given values of CaL and CaH. One can input "-1" to disregard a data pair. Users can 223 adjust limitation states according to how many points and the range of C i they have in their 224 We assume different processes in the C 4 photosynthesis are coordinated with each other and 227 co-limit the assimilation rate (Sharkey et al., 2007;Yin et al., 2011b;Ubierna et al., 2013;228 Bellasio et al., 2015). Thus, the estimation parameters allow the limitation states to be 229 compactly clustered with each other (Fig. 1). However, if there were only a few points under 230 CaL, the estimation results will depend heavily on the given initial values and unbalanced 231 results would be more likely. Fig. S1 shows an example of unbalanced estimation results by 232 deleting some points under 10 Pa or setting a very low CaL: in the estimation results, A n is 233 limited by AEE at very low C i and is mostly limited by A j (shown by AET and ATT) in the 234 C 3 cycle. In this case, A c (shown by AEE and ATE) has a clear redundancy at higher C i . 235 to minimize the following joint cost function (eq. 3 and 4) by varying the above five or six 252 output parameters (V cmax , J, R d , V pmax , g m , and k CA ): 253 n is the total number of observations, A ci is determined by AEE and ATE and A ji is 256 determined by AET and ATT from eq. (1), A mi is the observed net assimilation rate. 257 In this calculation, we take Michaelis-Menten constant of Rubisco activity for CO 2 (K c ), 258 Michaelis-Menten constant of Rubisco activity for O 2 (K o ), the specificity of Rubisco (γ*), 259 Michaelis-Menten constants of PEP carboxylation for CO 2 or HCO3 -(K p ), the fraction of O 2 260 evolution occurring in the bundle sheath (α) and bundle sheath conductance (g bs ) as given 261 (input parameters), similar to Sharkey et al. (2007). We conduct further sensitivity analyses 262 in the following section to determine the effects of variability of these inputs parameters on 263 the estimation results. anhydrase indeed shows limitation to V pc , which confirms its potential role as a limiting step 293 in the C 4 cycle. However, V pc calculated from CO 2 are only a little higher than V pc calculated 294 from HCO 3 -, which resulted in the similar estimation results. In addition, the estimation 295 errors and true errors from Yin's equations are quite small (average<1), and also similar 296 between models with and without carbonic anhydrase. 297 298 Estimation methods based on the two equations of different assumptions about electron 299 transport between RuBP regeneration and PEP regeneration yield consistent parameter 300 estimates and assimilation-CO 2 response curves (Fig. 2), but there were minor differences. 301 The second assumption that ATP, resulting from electron transport, is freely allocated 302 between PEP carboxylation-regeneration and RuBP regeneration leads to a bump at low 303 CO 2 when estimating ATE. The two assumptions produce different ATE under low CO 2 ; but 304 this is largely inconsequential because, under low CO 2 , assimilation is usually limited by 305 AEE.

Sensitivity analysis 318
The parameters K c , K o , γ*, K p , α, and g bs can vary among species in nature (Cousins et al., 319 2010) and it is therefore important to know how sensitive our results are to variation in these 320 parameters. We conducted a sensitivity analysis for variation in these parameters on the 1 7 estimated V cmax , J, R d , V pmax , g m and k CA (Fig. 3). This analysis shows all the estimated 322 parameters are robust under the variation of α (Fig. 3A) and showed little variation 323 responding to the change of γ* (Fig. 3E) and K o (Fig. 3C); however, the estimated 324 parameters are dependent on the other input parameters to different extents (Fig. 3B, D, F). 325 We calculate the average percentage change of estimated parameters along with the 50 % 326 decrease and 100 % increase of the input parameters. V cmax showed some medium extent of 327 sensitivity for K c , K p , and g bs with the average percentage change of 23.11, 7.54 and 17.69 % 328 respectively. J is robust in the variations of K c , and g bs (the average change is less than 2%) 329 and with a medium 6.96 % change for K p . k CA is robust in the variations of K c , K p , and g bs 330 (average change less than 5%). V pmax is sensitive for K p with the average change of 27.34%, 331 moderately sensitive to the change of g bs with 4.01 % and 13.38% change and is robust for 332 K c . R d is sensitive to K c , K p , and g bs with the change of 6.73, 43.88 and 13.38%. g m is 333 strongly sensitive to K c , K p , and g bs with the average percentage changes of 22.95, 107.04 334 and 23.19 %. This results suggest that V cmax , J, V pmax , and k CA estimated using our method 335 are relatively robust.

Physiological significance for assimilation rate of the input parameters 346
In addition to the sensitivity analysis, we performed a simulation analysis to illustrate the 347 physiological importance of input parameters further, and to indicate further the importance 348 of physiological properties in maintaining the efficiency of C 4 photosynthesis pathway. We 349 chose the estimation parameter set of T. dactyloides as an example, held photosynthetic 350 parameters constant V cmax (28 μmol m -2 s -1 ), J (134 μmol m -2 s -1 ), R d (0.78 μmol m -2 s -1 ), g m 351 (30.00 μmol m -2 s -1 Pa -1 ) and V pmax (41.91 μmol m -2 s -1 ), while changing the values of α, γ*, 352 g bs , and K p (as half or twice of the original parameters) to see their effects on the 353 assimilation rate, C bs and the O 2 concentration in bundle sheath (O bs ) (Fig. 4, Table 1). 354 Using photosynthetic parameter sets of other species to perform the simulation analysis 355 yielded similar results (data not shown). The change of α did not lead to changes in 356 assimilation rate (Fig. 4A) and led to small changes in O bs (Table 1). The decrease of γ* to 357 half of the current value led to a small change of C bs and assimilation rate (less than 0.5 358 μmol m -2 s -1 ) while doubling γ* led to a larger, but still not significant, change (less than 1 359 μmol m -2 s -1 ) (Fig. 4B, Table 1). Importantly, the changes of assimilation rates were less 360 than 0.3 μmol m -2 s -1 when C i was less than 20 Pa, which is the regular range of C i under 361 current ambient CO 2 . However, the change of g bs significantly changed the assimilation rate 362 2 0 and C bs (Fig. 4C, Table 1). The change of K p significantly affected the assimilation rate and 363 C bs to a large degree under low C i (Fig. 4D, Table 1). 364 365 Table 1 The average change of percentage of CO 2 concentration (C bs ) and O 2 concentration at 366 bundle sheath (O bs ) compared to the reference value of α0, γ*0, g bs 0 and K p . Simulation

Validating the estimation methods 380
In order to test our estimation methods, we first conducted a simulation test with 381 manipulated error terms. We use the estimated results of the nine species as known 382 parameters (the known values in Fig. 5) to generate new datasets using the C 4 383 2 2 photosynthesis equations based the first assumption of electron transport and adding error 384 terms to the assimilation rates. The error terms were randomly drawn from a normal 385 distribution of mean zero and standard deviation of 0.1 or 0.2 in an effort to simulate the 386 inevitable random errors in the real measurements. Estimating simulated data sets gave us an 387 idea about how likely we can capture the real parameters of the species given unavoidable 388 errors in measurements. The results show that robust estimation results for V cmax , J, V pmax , 389 and R d can be obtained (Fig. 5A, B, C, D). However, some estimation results of g m and k CA 390 show some deviation from the real values (Fig. 5E, F). 391

392
To test whether our estimation method could give accurate predictions across typical 393 prediction scenarios, (CO 2 ranging from 20 Pa to 60 Pa), we performed out of sample tests 394 for our nine target species. To perform these tests, we removed five points of CO 2 395 concentrations between 20 and 60 Pa range out of the A/C i curves and used the rest of the 396 A/C i curves to estimate parameters. And then we used these parameters to predict the 397 assimilation rate under the CO 2 concentrations we took out before and calculated the 398 estimation errors. In general, the estimation errors for all our species were small (Table 2). 399 400 Table 2 Out of sample test results. Five measured points from 20 Pa-60 Pa were taken out when we 401 conducted the estimation process. Then the calculated assimilation rates under these five CO 2 402 concentrations were compared with the measured ones. The data shows estimated error between the 403 calculated and measured assimilation rates (data show mean (standard error)).

Species
A  We tried to compare our estimation methods with in vitro measurements or other estimation 417 methods using isotopic analysis, especially for Zea. Our estimation results for Zea obtained 418 similar V cmax with the in vitro estimated Rubisco activity of Pinto et al. (2014); however, the 419 estimated value for V pmax is a little lower than the in vitro PEPC activity measurement with a 420 difference of around 10 μmol m -2 s -1 . For species of the Panicum family with NAD-ME 421 subtype, P. virgatum and P. amarum in the current study and P. coloratum in Pinto et al. 422 (2014), the estimated V cmax and V pmax are quite consistent with the in vitro measurements. 423 Ubierna et al. (2017) reported the g m for Zea ranged from 1.69 ± 0.17 to 8.19 ± 0.80 μmol 424 m -2 s -1 Pa -1 using 18 O and in vitro V pmax . Our estimation method fitted a g m for Zea of 7.34 425 μmol m -2 s -1 Pa -1 , which falls into the range of their measurements. Barbour et al. (2016) (Genty, Briantais & Baker 1989). Fluorescence analysis (Baker et al. 2007) is a 433 powerful tool for identifying the limitation states of C 3 species (Sharkey et al. 2007). If 434 Chlorophyll fluorescence is increasing with increasing CO 2 , A n is limited by Rubisco 435 carboxylation limited; when Chlorophyll fluorescence stays constant with increasing CO 2 , 436 A n is limited by RuBP regeneration. For C 4 species, however, the situation is more 437 complicated. Since V p could be limited by V pr and V pc (eq. (9)). Part of the RuBP 438 carboxylation limited condition and RuBP regeneration limited condition for the C 3 cycle 439 will mix together, leading to a linear increase of fluorescence with increasing of CO 2 , but of 440 a small slope (Fig. S2) increase the carbon gain. However, there is a trade-off between the specificity of Rubisco 521 for CO 2 and its catalytic rate (Savir et al., 2010;Studer et al., 2014). Based on this trade-off, 522 we can hypothesize that since C 4 elevates CO 2 around Rubisco relative to the O 2 523 concentration, maintaining low specificity might be optimal, in order to get high catalytic 524 rate of the enzyme to reach higher assimilation rate as shown by the empirical measurements 525 of Sage (2002) and Savir et al. (2010). Our simulation analysis showed the increase of 526 specificity for CO 2 (decrease of γ*) did not increase the assimilation rate much, which 527 indicates the selection upon Rubisco specificity in C 4 plants should be relaxed. (3) g bs 528 represents CO 2 leakage from bundle sheath to the mesophyll cell, and changes in g bs 529 significantly change the assimilation rate and C bs . Therefore, avoiding CO 2 leakage was of 530 great importance for the evolution and efficiency of C 4 photosynthesis pathway (Brown and 531 Byrd, 1993;Ubierna et al., 2013;Kromdijk et al., 2014). 532

533
Although we have shown that parameter estimation can be achieved solely with A/C i curves, 534 it is easy to combine our methods with ancillary measurements to yield more accurate 535 estimation results by defining the parameters as estimated or known or add additional 536 constraints (Supplementary Material IV). Yin et al. (2011b) proposed a method to obtain R d 537 from the fluorescence-light curve, since the method used for C 3 species, the Laisk method, is 538 inappropriate (Yin et al., 2011a). Additional measurement of dark respiration could be an 539 approximation for R d or could help to provide a constraint for estimating R d in our 540 estimating method. Ubierna et al. (2017) discussed the estimation method of g m using 541 instantaneous carbon isotope discrimination. With external measurement results, one can 542 change estimated parameters (such as R d , g m and J) as input parameters, instead of output 543