Short-term thermal photosynthetic responses of C4 grasses are independent of the biochemical subtype

In vitro and in vivo thermal photosynthetic responses vary among C4-grasses independently of biochemical subtype, and leakiness does not correlate with PEPC/Rubisco or initial slope/CO2-saturated photosynthesis rate.


Introduction
Understanding how photosynthesis responds to temperature is critical to our ability to predict the responses of natural and cropping ecosystems to climate change. Modelling the photosynthetic responses of C 3 plants to temperature is a routine task where we have a growing appreciation of the natural variation in the underlying mechanisms, such as the temperature dependence of mesophyll conductance (g m ) and the kinetics of the CO 2 -fixing enzyme ribulose 1,5-bisphosphate carboxylase/oxygenase (Rubisco) (Farquhar et al., 2001;Bernacchi et al., 2013;Walker et al., 2013;von Caemmerer and Evans, 2015;Sharwood et al., 2016a). Unlike its C 3 counterpart, the semi-mechanistic C 4 photosynthesis model (Farquhar, 1983;von Caemmerer, 2000) has not been comprehensively tested across a wide range of species and temperatures. Addressing this gap is critical given that C 4 plants include some of the world's most important food, feed, and biofuel crops (e.g. maize, sorghum, and sugarcane), dominate the understory of warm-climate grasslands and savannas (Ehleringer and Monson, 1993;Brown, 1999;Sage et al., 2013), and account for 20-25% of terrestrial productivity (Lloyd and Farquhar, 1994).
C 4 photosynthesis has evolved independently many times, resulting in multiple biochemical pathways named after the primary C 4 -acid decarboxylase enzyme found in the bundle sheath cells (BSCs), which are NADP-malic enzyme (-ME), NAD-ME, and phosphoenolpyruvate carboxykinase (PEP-CK) (Hatch 1987). Although PEP-CK operates as a secondary decarboxylase in many C 4 species (Leegood and Walker, 2003;Furbank, 2011;Sharwood et al., 2014;Wang et al., 2014), the primary decarboxylase is generally associated with a suite of anatomical, biochemical, and physiological features (Gutierrez et al., 1974;Hattersley, 1992;Kanai and Edwards, 1999;Ghannoum et al., 2011). The grass family includes species from all biochemical subtypes, but little is known about how these different pathways respond to temperature.
C 4 photosynthesis is characterized by the operation of a CO 2 -concentrating mechanism (CCM) comprising a C 4 cycle that fixes atmospheric CO 2 into a C 4 acid in mesophyll cells (MCs) that is then transported to, and decarboxylated in, BSCs where Rubisco and the C 3 cycle are localized (Hatch, 1987). The C 4 cycle is faster than the C 3 cycle, resulting in elevated CO 2 in the BSCs and thus minimizing photorespiration and maximizing CO 2 assimilation rates at low intercellular CO 2 (C i ) (Hatch, 1987;Kanai and Edwards, 1999). Some CO 2 leaks out of the BSCs into the surrounding MCs, costing additional ATP that is required for the phosphorylation of phospenolpyruvate (PEP) in the MCs (Hatch, 1987).
Leakiness (ϕ) describes the efficiency of C 4 photosynthesis and is defined as the rate of CO 2 leakage out of the BSCs as a fraction of the rate of PEP carboxylation or C 4 cycle rate (V p ) (Farquhar, 1983). According to the C 4 model, CO 2 leakage (L) depends on the BSC wall conductance to CO 2 (g bs ) and the gradient between the CO 2 concentration in the BSCs (C bs ) and MCs (C m ), such that L=g bs (C bs -C m ). In turn, the BSC-MC CO 2 gradient depends on the balance between the activity of the C 4 (e.g. PEP carboxylase, PEPC) and C 3 (e.g. Rubisco) cycles (Henderson et al., 1992;von Caemmerer, 2000). The C 4 photosynthesis model stipulates that the initial slope of the A-C i curve (IS) depends on the maximal PEPC activity (V pmax ) while the CO 2 -saturated rate (CSR) depends on maximal Rubisco activity (V cmax ), in addition to other factors (von Caemmerer, 2000). Consequently, the thermal responses of V pmax and V cmax are expected to influence the thermal response of both IS and CSR as well as the BSC-MC CO 2 gradient, and hence leakiness. Therefore, in this study we asked whether the thermal responses of these key parameters (V pmax , V cmax , IS, CSR, and ϕ) depend on the biochemical subtype of the C 4 species.
In a recent study, Sharwood et al. (2016a) showed that the declining response of Rubisco's in vitro CO 2 /O 2 specificity (S c/o ) with increasing temperature was more pronounced for NAD-ME relative to PEP-CK and NADP-ME species. Given that the CSR of C 4 leaves depends strongly on Rubisco and its catalytic properties, we predict that, relative to the other two subtypes, V cmax of NAD-ME species may increase less with temperature, with possible consequences on V pmax /V cmax and ϕ. These predictions remain untested, especially in that thermal responses of V pmax have been less widely studied than V cmax (Tieszen and Sigurdson, 1973;Pittermann and Sage, 2000;Sage, 2002;Boyd et al., 2015;Sharwood et al., 2016a).
In an earlier study, high leaf temperature was found to increase O 2 uptake in NAD-ME but not NADP-ME species, pointing to increased Mehler reaction and/or Rubisco oxygenation at higher temperature in NAD-ME species (Siebke et al., 2003). Higher Rubisco oxygenation may be symptomatic of a less favorable CO 2 /O 2 concentration ratio in the BSCs, which is corroborated by the less relaxed Rubisco affinity for CO 2 (higher K m(CO2) ) in NAD-ME relative to NADP-ME subtypes (Sharwood et al., 2016a). The suberin lamella lining the BSC walls of NADP-ME and PEP-CK grasses may reduce g bs , which could impact C bs (Hatch, 1987;von Caemmerer and Furbank, 2003). Consequently, we predict that C bs and hence leakiness may be affected differently by temperature depending on the C 4 subtypes. Previous measurements revealed a slight decrease in leakiness in response to increasing temperature between 21 and 35 °C for the C 4 monocot Sorghum bicolor (Henderson et al., 1992). To our knowledge, the thermal response of leakiness has not yet been compared in C 4 grasses with different biochemical subtypes.
The temperature dependence of key parameters of the C 4 photosynthesis model has primarily been determined for a few representative species, such as maize and Setaria (Henderson et al., 1992;Chen et al., 1994;Massad et al., 2007;Boyd et al., 2015). In addition, parameterization can be done in vivo (derived from gas exchange) or in vitro (derived from enzyme assays), but comparisons between these two parameterization methods remains unexplored for different C 4 species and temperatures. Relationships between ϕ and the ratio of in vivo (IS/CSR) or in vitro (V pmax /V cmax ) activities of C 4 and C 3 carboxylases have been demonstrated for single species under different growth environments (Ranjith et al., 1995;Saliendra et al., 1996;Meinzer and Zhu, 1998;Gong et al., 2017). However, such relationships have rarely been explored for diverse C 4 grasses at different temperatures. Addressing these knowledge gaps is critical to our ability to interpret gas exchange data and relate them to underlying biochemistry, as has been done for C 3 plants (von Caemmerer and Farquhar, 1981;Hudson et al., 1992;Bernacchi et al., 2001).
Consequently, in this study we measured the photosynthetic thermal responses of eight diverse C 4 grasses belonging to the three biochemical subtypes. We aimed at determining whether the C 4 biochemical subtype influences the thermal responses of (i) in vitro PEPC (V pmax ) and Rubisco (V cmax ) maximal activities, (ii) initial slope (IS) and CO 2 -saturated rate (CSR) derived from the A-C i curves, and (iii) CO 2 leakage out of the bundle sheath as estimated from carbon isotope discrimination. Our results showed that the thermal photosynthetic responses varied among the C 4 species independently of the biochemical subtype. We also observed that the ratios of IS/CS and V pmax /V cmax were uncorrelated with each other or with leakiness when determined for a range of C 4 species. Finally, we derived constants for thermal dependency that can be incorporated in the C 4 photosynthesis model.

Plant culture
The experiment was conducted in a naturally lit glasshouse chamber (5 m 3 ) with a mean noon photosynthetic photon flux density (PPFD) of 1230 ± 22 µmol quanta m −2 s −2 (±SE) and a 10-h photoperiod (see Supplementary Fig. S1 at JXB online). The air temperature inside the glasshouse compartment was regulated and day/night temperatures averaged 28/22 °C. Relative humidity was monitored and ranged between 60-80% during the day. Grass seeds (Table 1) were obtained from the Australian Plant Genetic Resources Information System, Australia. Seeds were sown in germination trays containing a common germination mixture. Seedlings at 3-4 weeks old were transplanted into the experimental pots (2 l) containing Osmocote ® Professional -Seed Raising & Cutting Mix (Scotts; www.scottsaustralia.com.au). Nutrients were supplied through the addition of Osmocote ® Plus Trace Element: Total All Purpose (N:P:K=19.4:1.6:5) (Scotts) and periodic watering with soluble Aquasol (N:P:K=23.3:3.95:14) (Yates; www.yates.com.au/). For each species there were eight pots, each of which contained a single plant. Pots were wellwatered and regularly rotated within the glasshouse chamber.

Leaf gas exchange measurements
Leaf gas exchange measurements were carried out using a portable open photosynthesis system (LI-6400XT, LI-COR, Lincoln, USA). Measurements were conducted between 10.00 h and 14.00 h, 7-8 weeks after transplanting, on the last fully expanded leaf (LFEL) attached on the main stem. Prior to gas exchange measurements, plants were moved to another chamber where air temperature could be changed separately. Measurements were made at leaf temperatures of 18, 25, 34, and 40 °C by using the internal heating system of the photosynthesis unit in conjunction with the glasshouse chamber heating system, whilst maintaining a relatively constant humidity inside the leaf chamber. Prior to measurements, each leaf was allowed to reach a steady state of CO 2 uptake at ambient CO 2 (Reference=400 μl l −1 ), PPFD of 1800 μmol m −2 s −1 , and relative humidity of 50-70%. A steady-state measurement was taken at each of the four leaf temperatures, and this was followed by measuring the responses of CO 2 assimilation rate (A) to step increases of intercellular CO 2 (C i ) by raising the LI-6400XT leaf chamber [CO 2 ] in 10 steps (i.e. 50, 100, 150, 200, 250, 325, 400, 650, 1200, and 1500 μl l −1 ) with 2 and 3 min as the minimum and maximum waiting times during each step change, respectively. For dark respiration, light in the LI-6400XT leaf chamber was switched off for 20 min before measurements were made. There were three or four biological replicates, i.e. plants per species. The initial slope (IS) of each A-C i curve was estimated by fitting a linear model to the initial 3-4 linear data points such that maximum C i was about 55 μl l −1 at all leaf temperatures. The maximum CO 2 -saturated rate of each A-C i curve at each leaf temperature was considered as the CO 2 -saturated rate (CSR).

Calculation of photosynthetic carbon isotope discrimination, leakiness, and C 4 cycle rate
Bundle-sheath leakiness was determined by measuring realtime 13 CO 2 / 12 CO 2 carbon isotope discrimination using a LI-6400XT attached to a tunable diode laser (TDL, model TGA100, Campbell Scientific, Inc., Logan, Utah, USA) under similar conditions to the steady-state measurements. For estimation of TDL precision, we used the δ 13 C value of the LI-6400XT reference gas for S. bicolor. The mean SD of repeated measurements for δ 13 C were 0.24, 0.14, 0.20, and 0.28‰ at 18, 25, 32, and 40 o C, respectively, with an overall SD of 0.21‰. Photosynthetic carbon isotope discrimination (∆) was calculated according to Evans et al. (1986): where δ e , δ o , C e , and C o are the δ 13 C (δ) and CO 2 mol fraction (C) measured with the TDL of the air entering (e) and leaving (o) the leaf chamber. Leakiness (ϕ) was calculated using the model of Farquhar (1983) as modified by Pengelly et al (2010) and Pengelly et al (2012). The slightly modified equation used here is: where Δ is the photosynthetic carbon isotope discrimination measured by the TDL. a i is the fractionation factor associated with the dissolution of CO 2 and its diffusion through water.
Here, we assume that s=a i . The term t, which represents ternary effects of transpiration rate on the carbon isotope discrimination during CO 2 assimilation, is defined according to Farquhar and Cernusak (2012) as: where E is the transpiration rate and g t ac is the total conductance to CO 2 diffusion including boundary layer and stomatal conductance (von Caemmerer and Farquhar, 1981).
The combined fractionation factor through the leaf boundary layer and stomata is denoted by a′: where: C a , C i , and C ls are the ambient, intercellular, and leaf surface CO 2 partial pressures, respectively; a b (2.9‰) is the fractionation occurring through diffusion in the boundary layer; s (1.8‰) is the fractionation during leakage of CO 2 out of the bundle sheath; and a (4.4‰) is the fractionation due to diffusion in air (Evans et al., 1986).
where: b 3 is the fractionation by Rubisco (30‰); b 4 is the combined fractionation of the conversion of CO 2 to HCO 3 − and PEP carboxylation (-5.74‰ at 25 °C); f is the fraction associated with photorespiration; Г* is the CO 2 partial pressure where rate of photorespiratory CO 2 release balances the rate of carboxylation; and C s is the CO 2 partial pressure in the BSC. The fractionation factor e associated with respiration was calculated from the difference between δ 13 C in the CO 2 cylinder used during experiments (-5.6‰) and that in the atmosphere under growth conditions (-8‰) (Tazoe et al., 2008). A and R d denote the CO 2 assimilation rate and daytime respiration, respectively; R d was assumed equal to the measured dark respiration (Atkin et al., 1997). We considered mesophyll conductance (g m ) to be 1.78 mol m −2 s −1 at 30 o C (Barbour et al., 2016;Ubierna et al., 2017). The temperature dependency of g m was accounted for by using the Arrhenius function: where T k is the leaf temperature in K; g m25 is the mesophyll conductance at 25 °C; E a is the activation energy, taken as 40.6 kJ mol −1 (the value for for Zea mays) (Ubierna et al., 2017); and R is the universal gas constant (0.008314 kJ K −1 ). In this study, leaf gas exchange was measured at high light and it was assumed that f C Γ * s =0 (Pengelly et al., 2010(Pengelly et al., , 2012Ubierna et al., 2013;von Caemmerer et al., 2014).
We used ϕ, A and R d to calculate the C 4 cycle rate (V p ), assuming that Rubisco oxygenation approximated to zero, i.e. V o ≈0 (Pengelly et al., 2012).

Determination of Rubisco and soluble protein contents
Following gas exchange measurements, replicate leaf discs were cut and rapidly frozen in liquid nitrogen and then stored at -80 °C until analysis. Each leaf disc was extracted in 0.8 ml of ice-cold extraction buffer [50 mM EPPS-NaOH (pH 7.8), 5 mM DTT, 5 mM MgCl 2 , 1 mM EDTA, 10 μl protease inhibitor cocktail (Sigma), 1% (w/v) polyvinyl polypyrrolidone] using a 2-ml Tenbroeck glass homogenizer kept on ice. The extract was centrifuged at 15 000 rpm for 1 min and the supernatant was used for enzyme activity (see below), Rubisco content, and soluble protein assays. For Rubisco content, subsamples were activated in buffer [50 mM EPPS (pH 8.0), 10 mM MgCl 2 , 2 mM EDTA, 20 mM NaHCO 3 ], and content was estimated by the irreversible binding of [ 14 C]-CABP to the fully carbamylated enzyme (Sharwood et al., 2008). Extractable soluble proteins were measured using the Pierce Coomassie Plus (Bradford) protein assay kit (Thermo scientific, Rockford, USA).

Temperature dependency
The temperature response of A, CSR, IS, V pmax , and V cmax were fitted in the R software (R Development Core Team, 2015) using two different equations, as follows. A modified form of the Arrhenius equation was used to fit the temperature dependence, which yields a peak function (Harley et al., 1992;Crous et al., 2013), and is given by the following equation: where, T k is the measurement temperature (either leaf or assay buffer) in K; E a is the activation energy (kJ mol −1 ) and represents the expansion for the initial part of the temperature response curve; k 25 is the parameter at 25 °C; H d is the deactivation energy (kJ mol −1 ); R is the universal gas constant; and ∆S (kJ mol −1 K −1 ) is the entropy term, which describes the peak part of the curve. H d and ∆S together describe the rate of decrease in the function above the optimum. To avoid over-parameterization, the H d of all parameters was set as constant (200 kJ mol −1 ) for model fitting, as has been done previously (Medlyn et al., 2002;Crous et al., 2013). The outputs from the modified Arrhenius and the June et al. (2004) equations are compared in Supplementary Fig.  S2. The modified form of the Arrhenius equation over-estimates the optimum parameter rate, P opt , at the temperature optimum, T opt . Consequently, T opt and P opt were derived by following the June et al. (2004) where T is the measurement temperature (either leaf or assay buffer) of parameter P in °C; P opt is the optimum parameter rate at the optimum temperature, T opt ; and Ω is the difference in temperature from T opt at which the parameter falls to e −1 (0.37) of its value at T opt . A smaller value of Ω means a narrower peak. This equation effectively assumes that the reversible processes are symmetrical around the optimum temperature.

Statistical analyses
The statistical design included species, subtypes, and leaf temperature as factors. Due to the limited number of species used, the statistical analyses could not incorporate phylogenetic or evolutionary effects. Gas exchange measurements and enzyme activity measurements were performed on three or four replicates at each temperature. The coefficients derived by fitting equations (10) and (11) for each parameter were used to test for differences between the thermal responses of species and subtypes. The effect of species was compared using a linear model with type-II ANOVA. The effect of subtype was compared using a linear mixed-effect model using the lme4 package (https://rdrr.io/cran/lme4/) in R (R Development CoreTeam, 2015). Significance tests were performed using type-II ANOVA. Coefficient means were ranked using post hoc Tukey tests.
Recent studies have highlighted the flexibility of the decarboxylases; in particular, the expression of PEP-CK activity in NADP-ME subtypes (Wingler et al., 1999;Furbank, 2011;Bellasio and Griffiths, 2014). Except for a relatively high activity in Z. mays (one-third of total decarboxylase activity), PEP-CK activity for most NADP-ME subtypes including those used in the current study does not exceed 5-25% of the total decarboxylase activity (see Supplementary Fig. S3). In addition, among the three species in the NADP-ME subtype, two are highly domesticated (S. bicolor and Z. mays), and Z. mays and Cenchrus ciliaris had appreciable secondary PEP-CK activity. However, the three NADP-ME species did not separate according to their known PEP-CK activity or domestication for the parameters collected ( Supplementary  Fig. S3). Consequently, for the purposes of the current study, it was valid to analyse the NADP-ME species as a single group rather than as one made up of multiple subgroups.

Leaf gas exchange at 25 °C
Leaf gas exchange parameters were measured at ambient CO 2 and near-saturating light intensity for all the C 4 grasses concurrently with stable carbon isotope discrimination. At 25 °C, net CO 2 assimilation rate (A), stomatal conductance (g s ), the ratio of internal to atmospheric CO 2 partial pressure (C i /C a ), dark respiration (R d ), and photosynthetic carbon isotope discrimination (∆) all varied significantly among the species (P<0.05; Fig. 1, Supplementary Table S1). Only A varied significantly (P<0.05) according to the C 4 subtype, with NAD-ME species having lower A relative to NADP-ME and PEP-CK species. Photosynthetic carbon isotope discrimination was lowest in Megathyrsus maximus and Z. mays, and highest in C. ciliaris, Chloris gayana, and Leptochloa fusca.

Temperature response of leaf gas exchange
The initial slope (IS) and CO 2 -saturated rate (CSR) of photosynthetic response to intercellular CO 2 partial pressure (from A-C i curves) were estimated at four temperatures (see Supplementary Fig. S4). At 25 °C, IS, CSR, and IS/CSR varied significantly with species but not with subtypes. Zea mays and L. fusca had the highest IS and IS/CSR, while M. maximus had the lowest CSR relative to the other C 4 species (Fig. 2, Supplementary Table S2).

Rubisco and PEPC measurements at 25 °C
Maximal activities of PEPC (V pmax ) and Rubisco (V cmax ) were measured on the same leaves used for gas exchange. At 25 °C, V cmax and V pmax were highest in the two NADP-ME subtypes S. bicolor and Z. mays. V pmax and V pmax /V cmax were lowest in the two NAD-ME subtypes along with M. maximus (PEP-CK) relative to the other species (Fig. 3, Supplementary  Table S2). In line with earlier studies, we obtained an average extraction yield of 80% for V cmax (maximal CO 2 assimilation/ Rubisco activity) for wild C 4 grasses (Meinzer and Zhu, 1998;Kingston-Smith et al., 1999;Crafts-Brandner and Salvucci, 2002;Dwyer et al., 2007).
Leaf Rubisco content varied among the species (3.8-7.8 µmol sites m -2 ), constituting about 6% of the soluble protein fraction (Supplementary Table S3). Among all the species, Rubisco and soluble protein contents were highest in Z. mays and Panicum coloratum and lowest in C. ciliaris, Eriochloa meyeriana, and L. fusca. Rubisco activation measured at 25 °C tended to be higher in NADP-ME (69%) compared to PEP-CK (53%) and NAD-ME (54%) subtypes.

Thermal responses of photosynthetic parameters
The short-term thermal responses of A, IS, CSR, V pmax , and V cmax were well characterized by the modified Arrhenius and June et al. (2004) equations (equations 10 and 11, Supplementary Figs S2 and S5-S8). No subtype effect was observed for the activation energy (E a ), entropy factor (∆S), parameter at 25 °C (k 25 ), the optimum parameter (P opt ) at the optimum temperature (T opt ), or the width of curvature around T opt (Ω) derived for IS, CSR, V pmax , and V cmax (Table 2).
For CO 2 assimilation rate, the range of variation observed for T opt (33-43 °C), E a (34-54 kJ mol −1 ), and Ω (18-24 °C) was not significantly different among all the species (Fig. 1A-C, Table 2). E a was significantly lower (P<0.05) in NADP-ME subtypes (37 kJ mol −1 ) relative to NAD-ME (52 kJ mol −1 ). Conversely, k 25 tended to higher (P=0.1) in NADP-ME subtypes relative to NAD-ME. For CSR, C. ciliaris had the highest P opt , T opt , and Ω, while P. coloratum and M. maximus had the lowest P opt , and Z. mays and P. coloratum had the Fig. 1. Thermal responses of leaf gas exchange and photosynthetic carbon isotope discrimination in eight C 4 grasses. (A-C) CO 2 assimilation rate, A, (D-F) stomatal conductance, g s , (G-I) ratio of intercellular to ambient CO 2 , C i /C a , and (J-L) photosynthetic carbon isotope discrimination, Δ, as a function of leaf temperature for the C 4 subtypes NADP-ME, PEP-CK, and NAD-ME (as indicated). The grasses were grown in a common glasshouse. Data in (A-F) were fitted according to June et al. (2004) and the derived constants are shown in Table 2. Leaves were measured at 1800 μmol m -2 s -1 PPFD and 400 μl l -1 CO 2. Values are means of 3-4 replicates ±SE.
lowest T opt . For IS, C. ciliaris had the highest T opt , P opt , and Ω, while E a was highest in S. bicolor and lowest in Z. mays and M. maximus relative to the other species (Fig. 2, Table 2). Overall, the ratio IS/CSR was not affected by leaf temperature, except in S. bicolor where IS/CSR increased with temperature (Fig. 2, Supplementary Table S2).
For the temperature response of in vitro V cmax , no parameter varied significantly according to the biochemical subtype, while E a , k 25 , and P opt varied among all the species (Table 2). Rubisco E a was highest in L. fusca (58 kJ mol −1 ) and lowest in S. bicolor (36 kJ mol −1 ), while Rubisco k 25 and P opt were highest in S. bicolor and Z. mays and lowest in C. ciliaris ( Fig. 3A-C, Table 2).
All parameters describing the thermal response of in vitro V pmax varied significantly among all the species. In addition, P opt and k 25 tended to be lowest (P<0.07) in NAD-ME, intermediate in PEP-CK, and highest in NADP-ME subtypes ( Fig. 3D-F, Table 2). Overall, the ratio IS/CSR was unaffected by leaf temperature (Fig. 3, Supplementary Table S2).

Photosynthetic carbon isotope discrimination and bundle sheath leakiness
Photosynthetic carbon isotope discrimination (∆) was unchanged between 25 and 40 °C and increased significantly at 18 °C for most of the species, except in S. bicolor and E. meyeriana where it decreased at 18 °C ( Fig. 1J-L Supplementary  Table S1). Overall, the NAD-ME subtypes showed higher Δ compared to NADP-ME and PEP-CK. Leakiness (ϕ) was higher at the two lowest temperatures (18 and 25 °C) relative to the two highest temperatures (34 and 40 °C), with ϕ at 40 °C being similar among all the species (Fig. 4A-  Thermal responses of the CO 2 -saturated rate (CSR) and the initial slope of the A-C i curve (IS) in eight C 4 grasses. (A-C) CO 2 -saturated rate, CSR, (D-F) initial slope of the CO 2 response curve, IS, and (G-I) the IS/CSR ratio as a function of leaf temperature for the C 4 subtypes NADP-ME, PEP-CK, and NAD-ME (as indicated). Data in (A-F) were fitted according to June et al. (2004) and the derived parameters are shown in Table 2. Leaves were measured at 1800 μmol m -2 s -1 PPFD. Values are means of 3-4 replicates ±SE both 18 and 34 °C. Generally, most leakiness values ranged between 35% at 18 °C and 10% at 40 °C ( Supplementary Fig.  S11). The estimated C 4 cycle rate (V p ) increased between 18 and 35 °C and was similar between 35 and 40 °C for all the species (Fig. 4D-F, Supplementary Table S1).

Relationships among photosynthetic parameters
In vivo estimates of CO 2 assimilation rate (A), stomatal conductance (g s ), CO 2 -saturated rate (CSR), initial slope (IS), and leakiness (ϕ), and in vitro measurements of V cmax and V pmax were used to assess how in vivo and in vitro photosynthetic parameters correlated with short-term changes in temperature.
Weak relationships were observed between IS and V pmax (r 2 =0.21, Fig. 5A) and between CSR and V cmax (r 2 =0.48, Fig. 5B), and a strong relationship was observed between A and g s (r 2 =0.78, Supplementary Fig. S10). Leakiness (ϕ) was neither correlated to V pmax /V cmax nor to IS/CSR (Fig. 6A, B). The ratios of IS/CSR and V pmax /V cmax were also not correlated (Fig. 6C).
To shed further light on what controls IS in diverse C 4 grasses, we used the simplified expression linking V pmax with IS derived by Pfeffer and Peisker (1998): where K p is the Michaelis-Menten constant for CO 2 . A weak fit was observed between measured V pmax values and model predictions when published K p values from Z. mays were used ( Supplementary Fig. S12A, dotted and dashed lines). When equation (12) was solved to simultaneously predict g m and K p using measured V pmax and IS values for Z. mays at all four temperatures together with published temperature dependencies of g m and K p (Bauwe, 1986;Boyd et al., 2015;Fig. 3. Thermal responses of photosynthetic enzyme activities in eight C 4 grasses. (A-C) Rubisco activity, (D-F) PEPC activity, and (G-I) the PEPC/ Rubisco activity ratio as a function of temperature for the C 4 subtypes NADP-ME, PEP-CK, and NAD-ME (as indicated). For each extract, the temperature responses of both enzymes were measured at 18, 25, 34, and 40 °C. Data in (A-F) are fitted according to June et al. (2004) and the derived parameters are shown in Table 2. Values are means of 3-4 replicates ±SE.

Table 2. Summary of thermal responses of photosynthetic parameters for eight C 4 grasses.
Coefficients are derived by fitting equation (10) (modified Arrhenius) and equation (11) (June et al., 2004). E a is the activation energy (kJ mol −1 ), ∆S is an entropy term (kJ mol −1 ), k 25 Ubierna et al., 2017), the observed fit between the measured and modelled data was much improved (Supplementary Fig.  S12A, continuous line, and S12B). Fitted values for Z. mays at 25 °C for g m and K p were then found to be 1.1 mol m −2 s −1 and 115 µbar, respectively (Fig. S12A).

Discussion
This study analysed the thermal photosynthetic responses of eight diverse C 4 grasses, deriving constants for thermal dependency that can be incorporated in the C 4 photosynthesis model (von Caemmerer, 2000). The main aim of the study was to determine whether the C 4 biochemical subtype influenced the thermal responses of in vitro (V pmax and V cmax ) and in vivo (IS, CSR, and ϕ) photosynthetic parameters that influence the efficiency of the C 4 CO 2 -concentrating mechanism (CCM). The study also explored whether the ratios of IS/CSR and V pmax /V cmax were correlated with each other or with leakiness across a range of C 4 species. The answers to these questions were largely negative, as discussed below.

Thermal photosynthetic responses varied among the C 4 grasses independently of the biochemical subtype
This study revealed large interspecific variations for most parameters derived from the photosynthetic temperature responses; however, these variations were largely independent of the C 4 subtype, disagreeing with our prediction (Table 2).
There was one exception to this generalization, with the activation energy (E a ) of CO 2 assimilation being lowest in the NADP-ME subtype. In addition, V pmax at 25 °C (k 25 ) and at T opt (P opt ) were marginally (P=0.07) lowest in NAD-ME and highest in NADP-ME subtypes. The latter observation is consistent with our previous reports of NAD-ME subtypes generally having lower V pmax and V pmax /V cmax (Pinto et al., 2016). The E a for in vitro V cmax (36-58 kJ mol −1 ) and V pmax (39-65 kJ mol −1 ) and their corresponding in vivo CSR (36-58 kJ mol −1 ) and IS (27-70 kJ mol −1 ) varied by 2-fold among the species examined here. Although these values are in line with published estimates for C 3 (V cmax ) and C 4 (V cmax and V pmax ) species, they tended to lie at the lower end of the spectrum reported for most of the C 4 grasses used in our study (Ishii et al., 1977;Bernacchi et al., 2001;Medlyn et al., 2002;Sage, 2002;Galmés et al., 2005Galmés et al., , 2015Massad et al., 2007;Walker et al., 2013;Perdomo et al., 2015;Sharwood et al., 2016a). For example, higher E a was reported for V cmax (78 kJ mol −1 ) and V pmax (95 kJ mol −1 ) activities in the NADP-ME grass Setaria viridis (Boyd et al., 2015), while relatively higher values of E a were reported for the in vivo V pmax of the C 4 species Z. mays and Andropogon gerardii (60-77.9 kJ mol −1 ) (Wu and Wedding, 1987;Chen et al., 1994). Interestingly, E a was higher for V cmax than V pmax in C. ciliaris, Ch. gayana, E. meyeriana, and L. fusca, while the opposite was observed in S. bicolor, Z. mays, M. maximus, and P. coloratum. Similarly, lower E a for V cmax than V pmax was reported in S. viridis (in vitro) (Boyd et al., 2015), and the opposite trend was reported in A. gerardii (in vivo) (Chen et al., 1994).
We predicted that thermal photosynthetic responses including leakiness may depend on the biochemical subtypes, based on known differences in Rubisco catalytic properties, PSII activity in the BSCs, suberization of the BSC walls, and possibly other anatomical and biochemical traits. In this study, similar thermal responses for leakiness were observed between the C 4 subtypes, and ϕ at 40 °C was similar for all the species. Comparable results were obtained in an earlier study using a smaller set of C 4 grasses measured at a common temperature (Cousins et al., 2008). These results suggest that CCM efficiency is similar among the C 4 subtypes despite their biochemical and anatomical differences. This was evident despite the potential for increased oxygenation in NAD-ME and PEP-CK species with increases in temperature (Siebke et al., 2003) as a result of significant PSII activity in the BSCs of these two subtypes (Edwards et al., 1976;Ghannoum et al., 2005). It has also been suggested that the absence of suberin in the BSC walls of NAD-ME subtypes is counterbalanced by the greater cytosolic barrier that exists for CO 2 diffusion through the centripetal arrangement of BSC chloroplasts surrounded by mitochondria towards the vascular bundle side (von . However, little is known about the thermal dependence of the CO 2 diffusion path from the BSCs to MCs in the different C 4 subtypes. In the current study, leakiness tended to be higher at the two lowest (18 and 25 °C) relative to the two highest (34 and 40 °C) temperatures for four out of the eight C 4 grasses (see Fig. 6. Relationships between leakiness and the activity ratio of in vitro or in vivo C 4 and C 3 cycle carboxylases. (A) Leakiness (ϕ) versus the ratio of the initial slope of the CO 2 response curve (IS) and the CO 2 -saturated rate (CSR), (B) ϕ versus the ratio of PEPC and Rubisco activities, and (C) the IS/CSR ratio versus the ratio of PEPC and Rubisco activities,measured at 18,25,34,and 40 °C (as indicated) in the C 4 subtypes NADP-ME (circles), PEP-CK (squares), and NAD-ME (triangles). The solid lines represent linear regressions of all data points. Supplementary Fig. S11). Similar results were previously reported for a C 4 monocot and a C 4 dicot species (Henderson et al., 1992). To gain further insights about the factors leading to higher leakiness at lower temperatures, we calculated the C 4 cycle rate (V p ). For all the species, the calculated V p increased with temperature, which indicates that increased leakiness at lower temperatures in some of the species is not related to increased pumping of CO 2 into the BSCs, but is probably due to increased Rubisco limitation. In agreement with this, Sage (2002) showed that C 4 plants are limited by V cmax at low temperature. Taken together, these findings indicate a greater Rubisco limitation at low temperature, which may lead to a greater proportion of CO 2 leakage out of the BSCs in some C 4 species (Kubien et al., 2003;Kubien and Sage, 2004), in a manner not explained by the biochemical subtype.
Maximal in vitro activities of PEPC and Rubisco, the initial slope and maximal rate of the A-C i curves, and leakiness are not correlated across a range of C 4 grasses and leaf temperatures Major progress has been made in our understanding of C 3 photosynthesis following the development of a fully mechanistic model (Farquhar et al., 1980) that was subsequently validated by combining leaf gas exchange measurements with estimates of Rubisco activity and mesophyll conductance (g m ) in wild-type (von Caemmerer and Farquhar, 1981;von Caemmerer et al., 1994;Bernacchi et al., 2002;von Caemmerer and Evans, 2015) and genetically altered plants (Hudson et al., 1992). For C 4 photosynthesis, model validation has been undertaken mostly under standard conditions using a genetically altered C 4 dicot, Flaveria bidentis (von Caemmerer et al., 1997;Kubien et al., 2003;Pengelly et al., 2012), or under a range of conditions using a limited number of model C 4 grass species, such as S. bicolor (Henderson et al., 1992), Z. mays (Yin et al., 2016) and S. viridis (Boyd et al., 2015). This study presented an ideal opportunity to undertake model validation using in vitro enzyme assays and leaf gas exchange measurements for a range of C 4 species.
A weak correlation was observed between CSR and V cmax , while none was found between IS and V pmax (Fig. 5), or between E a of V cmax (or V pmax ) and E a of CSR (or IS) across the various C 4 grasses. This discrepancy is not surprising and reflects the multitude of factors controlling the A-C i curves besides the maximal in vitro activity of Rubisco and PEPC. Other than Rubisco activity and level of activation, CSR is also dependent on RuBP and PEP regeneration, and both limitations are generally indistinguishable at high light (von Caemmerer, 2000). The IS is determined by the maximal PEPC activity (V pmax ), its Michaelis-Menten constant for CO 2 (K p ), and the mesophyll CO 2 partial pressure (C m ), which depends on mesophyll conductance (g m ) (von Caemmerer, 2000). The parameters g m and K p were not measured in this study, but using published values at 25 °C together with their temperature dependencies for Z. mays (Bauwe, 1986;Boyd et al., 2015;Ubierna et al., 2017) to model the relationship between IS and V pmax , a weak fit was observed between measured and predicted V pmax values ( Supplementary Fig. S12A, dotted and dashed lines). In contrast, when model predictions were made using simultaneously solved g m and K p values using measured V pmax and IS for Z. mays at all four temperatures, the observed fit between measured and modelled V pmax was much improved (Supplementary Fig. S12A, continuous line, and S12B). This exercise indicated that g m and K p as well as their temperature responses may vary with growth environment and species (or genotype), as has recently been shown for Rubisco kinetics among these C 4 grasses (Sharwood et al., 2016a). The variation in g m and its temperature response in C 3 species is well established (von Caemmerer and Evans, 2015), and discrepancies for g m in C 4 Z. mays have been reported (Barbour et al., 2016;Ubierna et al., 2017).
In C 3 species, IS (Rubisco-limited photosynthesis) is largely insensitive to temperature because both the K m and V cmax of Rubisco have a Q 10 of 2 (Berry and Raison, 1981;Sage and Sharkey, 1987). In contrast, we observed a strong temperature response of IS (PEPC-limited) for the C 4 species examined, indicating that the E a for V pmax is greater than the E a of K p (Table 2, Fig. 2D-F). A similar trend has been shown in a recent study (Boyd et al., 2015). In addition, variation among C 4 species for the E a of IS suggests that the ratio between E a for V pmax and K p also differs. This has been partly supported by our data in the case of the E a for V pmax (Table 2, Fig. 3D-F). These findings indicate that there is a significant diversity for PEPC kinetics among C 4 species, and warrants further investigation.
The improved relationship between measured and modelled V pmax relative to that observed between IS and measured V pmax ( Supplementary Fig. S12A versus S12B) demonstrates that the C 4 photosynthesis model, when well parameterized, can accurately predict CO 2 assimilation rates at limiting C i in C 4 leaves. Given the rapidly improving methodologies for measuring the exchange of 13 C and 18 O stable isotopes that allow us to estimate both g m and K p (von Caemmerer et al., 2004;Boyd et al., 2015;Barbour et al., 2016;Ubierna et al., 2017), the C 4 photosynthesis model will form a valuable tool for interpreting leaf gas exchange data for any C 4 species under a wide range of environmental conditions, similarly to what has been widely achieved for C 3 plants. Further parameterization and validation of the C 4 photosynthesis model for diverse C 4 species remains a focus of our current research.
Bundle sheath leakiness (ϕ) depends on g bs and the (C bs -C m ) gradient, which in turn depends on the balance between the activities of PEPC and Rubisco (Henderson et al., 1992;von Caemmerer, 2000). In accordance with model predictions, earlier studies reported a relationship between leakiness and in vitro V pmax /V cmax (Ranjith et al., 1995;Saliendra et al., 1996;Meinzer and Zhu, 1998) and in vivo IS/CSR (Gong et al., 2017). These studies used a single species (sugarcane or Cleistogenes squarrosa) exposed to soil or atmospheric water deficit, and estimated ϕ from the C-isotope composition of leaf dry matter rather than photosynthetic C-isotope discrimination as done here. In the current study, leakiness was neither correlated to in vitro V pmax /V cmax nor to IS/CSR when all species and leaf temperatures were considered (Fig. 6). There were three exceptions at the species level: leakiness was positively correlated with V pmax /V cmax (but not with IS/CSR) in Ch. gayana, L. fusca, and Z. mays (see Supplementary Figs S13 and S14). Hence, we argue that, in response to short-term changes in leaf temperature, CCM efficiency was balanced by biochemical factors, such as enzyme activity, as well as by physical factors, such as g bs (von Caemmerer, 2000). A contribution by g m is unlikely because very little influence of g m on the carbon isotope discrimination (∆) has been predicted (von Caemmerer et al., 2014).

Conclusions
Modelling thermal photosynthetic responses at the leaf level is critical for predicting canopy-scale gas exchange in response to diurnal and seasonal changes in leaf temperature (Harley and Baldocchi, 1995). Thermal sensitivities of parameters used by the C 4 photosynthesis model are needed to accurately predict CO 2 exchange in response to temperature. The findings from the current study demonstrated that, like C 3 photosynthesis (Bernacchi et al., 2001;Walker et al., 2013), in vivo and in vitro thermal responses of key photosynthetic parameters (e.g. PEPC and Rubisco activities) differ across C 4 species. Hence, relying on the thermal responses of selected species to model C 4 photosynthesis cannot accurately describe ecosystem responses.
The current study also demonstrated that variations in the thermal responses of leakiness (ϕ) among the C 4 grasses is not aligned with the C 4 subtypes. This indicates that various biochemical and anatomical trade-offs operate to maintain similar CCM efficiencies in the various C 4 biochemical pathways. In addition, no correlations were observed among leakiness, V pmax /V cmax , and IS/ CSR across a range of C 4 grasses and leaf temperatures. Hence, more work is needed to characterize the thermal responses of g m and g bs in diverse C 4 species by combining stable isotope and chlorophyll fluorescence studies (Yin et al., 2016).

Supplementary data
Supplementary data are available at JXB online. Table S1. Summary of leaf gas exchange parameters for eight C 4 grasses. Table S2. Summary of A-C i -derived parameters and enzyme activities for eight C 4 grasses. Table S3. Rubisco and protein content of leaves. Fig. S1. Light environment in the glasshouse during plant growth. Fig. S2. Comparison of the June et al. (2004) and modified Arrhenius models for temperature response. Fig. S3. Principle component analysis plots for species of the NADP-ME subtype, and all three C 4 subtypes. Fig. S4. Photosynthetic CO 2 response curves (A-C i ) measured at four-leaf temperatures in eight C 4 grasses. Fig. S5. Thermal responses of the CO 2 assimilation rate (A), CO 2 -saturated rate (CSR), and initial slope of the A-C i curve (IS) in eight C 4 grasses fitted using the June et al. (2004) model. Fig. S6. Thermal responses of photosynthetic enzyme activities in eight C 4 grasses fitted using the June et al. (2004) model. Fig. S7. Thermal responses of the CO 2 assimilation rate (A), CO 2 -saturated rate (CSR), and initial slope of the A-C i curve (IS) in eight C 4 grasses fitted using the modified Arrhenius model. Fig. S8. Thermal responses of photosynthetic enzyme activities in eight C 4 grasses fitted using the modified Arrhenius model. Fig. S9. Photosynthetic carbon isotope discrimination, Δ, as a function of C i /C a measured during the gas exchange for eight C 4 grasses. Fig. S10. Relationship between CO 2 assimilation rates and stomatal conductance in eight C 4 grasses. Fig. S11. Thermal response of leakiness in eight C 4 grasses. Fig. S12. Relationship between measured and modelled PEPC activity at various leaf temperatures in Z. mays. Fig.S13. Relationship between leakiness and ratio of PEPC to Rubisco activity in eight C 4 grasses. Fig. S14. Relationship between leakiness and ratio of initial slope of A-C i (IS) to CO 2 -saturated rates in eight C 4 grasses.