Stomatal response to VPD in C4 plants with different biochemical sub-pathways

C4 plants are integral to many ecosystems around the world and are abundantly found in grasslands and savannas in North America, Australia, and Africa, and deserts in Central Asia among other ecosystems in and around Europe (Edwards and Still 2008; Pyankov et al. 2010; Rudov et al. 2020; Wan et al. 2001). Despite making up 20% of the global plant biomass (Ehleringer et al. 1997), our understanding of C4 physiology and biochemistry still has room for improvements. More specifically, the three subtypes of C4 plants (NADP-malic enzyme, NAD-malic enzyme, and PEP carboxykinase) should be better characterized in terms of their biochemistry. To do that, we require a comprehensive characterization of species from across C4 families, not just model C4 species. However, it is often labor-intensive to collect the full gamut of physiological and biochemical data for large numbers of species. Here, I describe a hierarchical Bayesian approach in parametrizing the C4 photosynthesis model that can estimate biochemical parameters using gas exchange data while accounting for plant, species, and subtype variation.

Methods
To enable a comparative analysis, we compared 6 closely related C4 lineages that each contain an independent origin of NAD-ME and NADP-ME subtypes, for a total of 12 species. 

We analyzed the response of net CO2 assimilation rate (A) to intercellular concentration of CO2 (Ci) and we measured the response of stomatal conductance (gs) to VPD. All gas exchange measurements were conducted with two Li-6400 Portable Photosynthesis Systems (Li-Cor Biosciences, Lincoln, NE, USA) over a course of five months. We conducted measurements at a constant leaf temperature of 30 °C, PPFD of 1800 µmol m-2s-1, and CO2 concentration at ambient level of 400 ppm, unless otherwise noted. 

We used a hierarchical Bayesian approach to model C4 photosynthesis to estimate maximum PEPC activity (Vpmax) and maximum Rubisco activity (Vcmax) from the A/Ci data. The hierarchical Bayesian approach followed methods of Patrick et al. (2009) and allowed us to estimate multiple parameters (namely Vpmax and Vcmax) simultaneously, as well as accounting for observational error in the data set while identifying uncertainty in the estimated parameters. We used a C4 photosynthesis model (von Caemmerer 2000) that was simplified to only estimate PEPC-limited and Rubisco-limited values of A. From the model, we estimated the Ci where the transition from PEPC-limited to Rubisco-limited A occurred. Our hierarchical model had three main components: the likelihood equation of observing the data, the deterministic model which is the C4 photosynthesis model (von Caemmerer 2000), and the prior distributions of the variables involved in the deterministic model. These components generated posterior distributions of Vpmax, Vcmax, and the transition Ci between the limitations. Because this was a hierarchical model, we were able to set parameters to vary at either the species level or the plant level. Vpmax, Vcmax, and the transition point were all set to vary at the plant level such that each plant can have its own distinct value. The Michaelis-Menten constants for CO2 for PEPC and for O2 and CO2 for Rubisco were allowed to vary at the species level such that all the replicates within a species shared the same value. All other parameters in the model were included as constants. Because there were small variations in temperature during gas exchange measurements (28.3 – 32.1 °C) and the oxygenation and carboxylation rates of Rubisco and carboxylation rates of PEPC are temperature dependent, we used an Arrhenius function to standardize all values to 30 °C (von Caemmerer, 2000). The hierarchical Bayesian photosynthesis model was implemented using JAGS (Plummer 2003) interfaced with RStudio using the package “rJAGS” (Plummer et al. 2022). We ran three parallel Markov chain Monte Carlo (MCMC) chains for 30 000 iterations each and the chains were evaluated for convergence using the Gelman-Rubin diagnostic (Gelman and Rubin 1992). The burn-in samples (first 5000) were discarded. The model goodness-of-fit was evaluated by generating predicted A values and comparing these with observed A values. If the model was a good fit, then the predicted A and observed A should fall along the 1:1 line.