C₄ 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 C₄ physiology and biochemistry still has room for improvements. More specifically, the three subtypes of C₄ 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 C₄ families, not just model C₄ 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 C₄ photosynthesis model that can estimate biochemical parameters using gas exchange data while accounting for plant, species, and subtype variation.

To enable a comparative analysis, we compared 6 closely related C₄ 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 CO₂ assimilation rate (A) to intercellular concentration of CO₂ (C_i) and we measured the response of stomatal conductance (g_s) 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 CO₂ concentration at ambient level of 400 ppm, unless otherwise noted. 

We used a hierarchical Bayesian approach to model C₄ photosynthesis to estimate maximum PEPC activity (V_pmax) and maximum Rubisco activity (V_cmax) from the A/C_i data. The hierarchical Bayesian approach followed methods of Patrick et al. (2009) and allowed us to estimate multiple parameters (namely V_pmax and V_cmax) simultaneously, as well as accounting for observational error in the data set while identifying uncertainty in the estimated parameters. We used a C₄ 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 C_i 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 C₄ photosynthesis model (von Caemmerer 2000), and the prior distributions of the variables involved in the deterministic model. These components generated posterior distributions of V_pmax, V_cmax, and the transition C_i 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. V_pmax, V_cmax, 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 CO₂ for PEPC and for O₂ and CO₂ 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.