In deriving the following kinetic equations for this mechanism (Figure 1) we assumed only that both k₁₀ and k₁₆ are very much smaller than any of the remaining rate constants (effectively, k₁₀ = k₁₆ = 0). We emphasize that no such approximations (k_i = 0) were made anywhere else in the derivation. The Michaelis constants (K_M) for carboxylation and oxygenation are then given, respectively, by equations of the form (Equations A23, A24; see Appendix in Supplementary Materials for details of derivations)

The general equation for the specificity of carboxylation relative to that of oxygenation (relative specificity) is then (Equation A25).

In Equation (3), the relative specificity (S_C/O) is formally a function of 10 rate constants (k₅..k₉, k₁₁..k₁₅), five for each of the carboxylation and oxygenation reactions. KR=k3(k3+k4) is a function only of rate constants for the enolization step (Equation A22), i.e., independent of carboxylation or oxygenation, and 0 < γ < 1. Both kcatC and γ_C are formally functions of k₃, k₇, k₈ and k₉ (Equation A26). It is evident (Equation A26) that if k₇ is the slow step that determines the maximum catalytic rate (kcatC=k7), then γ_C = 1. Similarly (Equation A27), if kcatO=k13, then γ_O = 1. However, we need not make these types of assumptions here, and simply regard γ_Ck₆ and γ_Ok₁₂ as effective dissociation rate constants.

The results of linear regression analysis performed on a number of data sets are summarized in Table 2. The green algae, bacteria and cyanobacteria data in Table 1 and other plant species (Galmés et al., 2014) could not be considered individually for analysis due to the small numbers of observations (N < 3). The log-scale plots (Figures 2A,B) of K_C over the full range of kcatC values in Table 1 suggest a linear correlation and hence regression analysis of ln(K_C) on kcatC (“All data” sets in Table 2). P < 0.05 for both coefficients were obtained only for carboxylation using the “All data” sets (Figures 2A,B), carboxylation using a subset of the C₃ plants (Galmés et al., 2014), oxygenation using Triticeae data (Prins et al., 2016) and oxygenation using only the higher plant data (Figure 2C). The residuals were found to be near-normally distributed (Figure 3). Reliable expected values for effective CO₂ and O₂ dissociation rate constants can be derived from the coefficients in regressions (Table 2) that yield P < 0.05 for both coefficients (i.e., both the gradient and intercept). The results are given in Table 3. For the regression of ln(K_C) on kcatC, equating the first terms (a1+a1b1kcatC) in the expansion of the exponential form (a1eb1kcatC) with Equation (1) we find that the value of K_C at kcatC=0 is given by 〈γCk6KRk5〉=a1. From the regression analysis carried out using the full data set in Table 1 (Figure 2A) and the subset utilized by Savir et al. (2010) (Figure 2B), we obtain values of a₁ = 9.7 μM and a₁ = 4.5 μM, respectively. From the expansion of the exponential we also find that 〈1KRk5〉≈a1b1 at kcatC=0, where the two estimates are a₁b₁ = 2.2 μM.s and a₁b₁ = 1.5 μM.s, respectively. In Figures 2A,B, Equation (1), which will obviously deviate from the trend line as kcatC increases, has been graphed using these values. Combining these results obtained for 〈γCk6KRk5〉 and 〈1KRk5〉 we estimate (at kcatC=0) expected effective rate constants for CO₂ dissociation (〈γ_Ck₆〉) of 4.3s⁻¹ and 3.0s⁻¹, respectively. Assuming the scheme (Figure 1) correctly describes the kinetic mechanism, the deviation from linear behavior suggests there exists at least one type of correlation between rate constants. From Equation (1), the expected value of K_Rk₅ conditional on kcatC in terms of regression parameters a₁ and b₁ is then given by (Figure 4A).

Therefore, we may also use Equation (4) to define the expected effective dissociation constant conditional on kcatC as (Figure 4B).

In Figure 4 it is assumed (Tcherkez et al., 2006) that the exponential increase in 〈K_C〉 conditional on kcatC arises from 〈K_Rk₅〉 (one correlation effect, i.e., due to CO₂ binding) while 〈γ_Ck₆〉 is a constant in Equation (4). Alternatively, in Figure 5 we have assumed that variation arises from 〈γ_Ck₆〉 (another correlation effect i.e., due to CO₂ dissociation) while 〈K_Rk₅〉 is now the constant. Here the respective constants are the values of 〈γ_Ck₆〉 and 〈K_Rk₅〉 at kcatC=0 as determined from the regression (Figure 2B). There is, of course, also the possibility that variability in both 〈K_Rk₅〉 and 〈γ_Ck₆〉 contribute to the non-linear behavior of 〈K_C〉, i.e., both 〈K_Rk₅〉 and 〈γ_Ck₆〉 are conditional on kcatC. In general, therefore, we could ascribe any functional dependence for either 〈K_Rk₅〉 or 〈γ_Ck₆〉 to this non-linear behavior.

For the regression of K_O on kcatO (Figure 2C), we have included only the data for all higher plants (Table 1). Unlike the above regressions of K_C on kcatC there are no indications of any deviations from non-linear behavior. The graph of K_O on kcatO for the higher plants in particular clearly conforms to a linear function, and the residuals of regressed K_O data are near normally distributed (Figure 3). From the intercept we find the expected value of the dissociation constant

and from the gradient we obtain the constant

From Equations (6, 7) we estimate the expected value of the effective O₂ dissociation rate constant, 〈γOk12〉≈0.3 s-1. Finally, from the above determinations of 〈1KRk5〉 (from Figure 2B) and 〈1KRk11〉 we can estimate the expected CO₂ to O₂ ratio of the rate constants for binding at kcatC=0 as 〈k5k11〉≈190.

The graph of reciprocal relative specificity, SO/C=1SC/O, against kcatC (Figure 2D) suggests a linear dependence. The residuals of regressed S_O/C data are near normally distributed (Figure 3). We first consider the expected value of S_C/O conditional on kcatC as the reciprocal of the equation for the straight line that describes 〈S_O/C〉, i.e.,

where a2=7.4×10-3 mol/mol and b2=1.2×10-3s are the regression parameters (Figure 2D). Although Equation (8) generally provides a good fit to the data (Figure 6), it clearly does not display the correct limiting behavior as kcatC approaches zero Equation (3). However, defining the expected value as the ratio 〈SC/O〉=〈SC〉〈SO〉 and substituting 〈SC〉=kcatCa1eb1kcatC (Figure 2A), the expected value of S_C/O conditional on kcatC can be written as

As there are no correlations between kcatC and kcatO (Figure 7A) or K_O (Figure 7B), S_O is also not correlated (Figure 7C), and so the best possible approximation for Equation (9) takes the form S_C/O ∝ S_C. The constant 〈SO〉-1 in Equation (9) can therefore be estimated by a linear regression of S_C/O (excluding the outlier, R. rubrum, Figure 2D) on kcatCa1eb1kcatC subject to the constraint S_C/O = 0 at kcatC=0 to obtain the correct general equation for the expected value of S_C/O conditional on kcatC as (Figure 6).

Assuming correlation (Figure 2B) arises from CO₂ binding, the factor implicit in Equation (10) corresponding to 〈k5k11〉 (Figure 6) that is also conditional on kcatC is estimated by (Equations 4, 7, Figure 4A).

We use Equation (3) to rationalize the in vitro kinetic data for the Leu to Val mutation at position 335 (L335V) in tobacco (Whitney et al., 1999). The decrease in kcatC from 3.43 s⁻¹ in the wild type to 0.81 s⁻¹ in the mutant is accompanied by a large decrease also in S_C/O from 81 to 20 mol/mol. In Figure 8, S_C/O is plotted against kcatC assuming that in Equation (3) the term k5(kcatO+γOk12)k11kcatO is constant on the curve, i.e.,

We determine the constant factor such that S_C/O = 81mol/mol for the wild-type tobacco at the two limits (kcatC≫γCk6 and kcatC≪γCk6) for specific values of γ_Ck₆ = 1, 2, 3 and 4s⁻¹. Note that in the limit kcatC≫γCk6 we obtain SC/O=k5k11=81 mol/mol, while the lower limit for kcatC≪γCk6 gives S_C/O = 0. Noting that 〈k5k11〉=〈KRk5〉〈KRk11〉, the remaining kinetic parameters [K_C = 10.7 μM, kcatO=1.17s-1, K_O = 295 μM for wild type, and K_C = 5.1 μM, kcatO=0.39s-1, K_O = 48.9 μM for the mutant] (Whitney et al., 1999) can be used to simply determine the expected value of the ratio k5k11 as.

where Δ is the difference between wild type and mutant.

The trend lines (Figure 2) clearly intercept the vertical axes well above zero, indicating significant expected values for the dissociation constants γ_Ck₆ and γ_Ok₁₂. However, the rate constant for CO₂ dissociation has been previously estimated as not more than about 5% of kcatC (Pierce et al., 1986; McNevin et al., 2007), so that it has generally been assumed that kcatC(kcatC+γCk6)≈1. Our estimates (Figures 2A,B, Table 3) of the expected value (at least for low kcatC) are much higher, and find support in the kinetics modeling study of RuBisCO from spinach. We find that the expected values of dissociation rate constants (γ_Ck₆) for the binding of the substrate CO₂ are 4.3 s⁻¹ (Figure 2A), 3.0 s⁻¹ (Figure 2B), and 3.1 s⁻¹ for a subset of C₃ plants (Galmés et al., 2014; Prins et al., 2016), noting that the differences are not statistically significant (Table 3). These values can be compared with 1.6 ±1.1 s⁻¹ estimated for the CO₂ dissociation rate constant in spinach (McNevin et al., 2006, 2007), and the 5−10 μM range of 〈KDC〉 for lower values of kcatC (Figures 4B, 5B) is also consistent with a KDC=k6k5 of 3 μM for spinach RuBisCO (McNevin et al., 2006). The effective CO₂ dissociation rate constant, γ_Ck₆, impacts the kcatC dependence of S_C/O (Figure 8). As kcatC approaches γ_Ck₆ Equation (12) describes the rapid decline in S_C/O due to increasing probability that the CO₂ will dissociate from RuBP before catalysis takes place. The observed values of kcatC and S_C/O for the L335V mutant (Whitney et al., 1999) are entirely consistent with a γ_Ck₆ greater than kcatC. The expected value of k5k11 as given by Equation (13) is also consistent with the value obtained when averaged over a larger number of RuBisCOs with lower kcatC (Figure 6). Thus, changes in the gas-substrate binding in the mutant RuBisCO appear to be minimal, the bulk of the effect being described by Equation (12). The dissociation rate constant of O₂ is generally considered effectively zero (Tcherkez, 2013, 2016). However, although the expected value of 0.4 ± 0.4 s⁻¹ for γ_Ok₁₂ in higher plants obtained here (Table 3) is significantly lower than the mean kcatO of 1.3 ± 0.2 s⁻¹ (from data in Table 1) it is still sufficient to have an impact on K_O (Equation 2). Additionally, the expected value of 2.3 ± 1.9 s⁻¹ for γ_Ok₁₂ in Triticeae (Table 3) and the corresponding mean kcatO of 0.83 ± 0.16 s⁻¹ (Prins et al., 2016) are not significantly different. Statistical analysis of the available data therefore suggests the expected (or average) value of the dissociation rate is not significantly lower than that of the catalytic rate. Moreover, a knowledge of rate differences in any particular RuBisCO requires more kinetic data than is currently available. Consequently, there is no justification for generally neglecting either of the dissociation rate constants, γ_Ck₆ or γ_Ok₁₂, i.e., assuming they are an order of magnitude or more lower than the corresponding catalytic rates, as has been done previously (Tcherkez, 2013, 2016).

Assuming 〈k5k11〉 decreases with increasing kcatC (Equation 11, Figure 6), it could be regarded as a proxy for S_C/O (Tcherkez, 2013). Also as the specificity of oxygenation, S_O, is not correlated with kcatC (Figure 7C), the variation in 〈k5k11〉 would be largely constrained to the dependence of 〈k₅〉 on kcatC (Figure 4A). It has been hypothesized (Tcherkez et al., 2006; Tcherkez, 2013) that such a constraint is to be expected from the predicted energetics of the reaction as tighter binding of CO₂ to ribulose bisphosphate (increasing k₅) would necessarily raise the activation free energy (decreasing kcatC) required for the subsequent steps leading to turnover of product. However, the generality of this tight-binding (TB) hypothesis has come under question (Hanson, 2016) for its inability to explain the variations in S_C/O that have been observed in some RuBisCOs (Young et al., 2016). It would seem that the TB hypothesis suffers from a more fundamental problem in that it is based on an incomplete and unrepresentative data distribution. In the present analysis, Equation (8) provides the better fit R² = 0.63 to the selected data (Figure 2D), although it is not the more general equation for 〈S_C/O〉 (Equation 10, Figure 6). Similar types of relationships that provide an even tighter fit to the data have been reported elsewhere: KC∝(kcatC)2 [R² = 0.90] and SC/O∝(kcatC)-0.51 [R² = 0.79] (Savir et al., 2010). The TB hypothesis is posited on k5k11 determining the dependence of S_C/O on kcatC. Significantly, all of these analyses are in fact conditional on kcatC≫γCk6, i.e, neglect of the CO₂ dissociation rate constant, k₆. However, the high level of variance in K_C and S_C/O (Figures 2A, 6, respectively) argues for a more cautious data interpretation in the regression analysis. Statistically, the quadratic (Savir et al., 2010) and exponential (Figure 2B) forms both describe the dependence of K_C on kcatC equally well, but only the latter, more general case (Equations 3, 10), allows nonzero values for γ_Ck₆ (Figure 5A).

The deviation of any given data point (Figure 6) from the expected value (Equation 10) can be attributed to variations in the parameters of Equation (3). We expect that S_O will generally produce random variations in S_C/O (Figure 7C), although, possibly lower k₁₁ (higher K_O, Figure 2C) for the cyanobacteria may in part account for a systematic reduction in S_C/O. The CO₂ dissociation term, γ_Ck₆, will certainly become apparent at low enough kcatC values (Figures 4, 5). In particular, variations in γ_Ck₆ may contribute significantly to the large variance seen in the non-green algae (Figures 2A, 6). If the catalytic rate correlates with k₅, regression analysis defines only the first moment, 〈k₅〉, of the distribution (Figure 4A and Equation 11, Figure 6), and provides no information on the variance. In the absence of any coupling, mutations produce random changes in the underlying rate constants, k_i. Irrespective of whether rate constants are correlated, the expected value of k_i is given by ∑snkisn where kis is the value of a rate constant for a given sequence (s). In reality, the composition of the sequence space, Ω (i.e., any number of known sequences), will be determined in varying degrees by genetic drift and natural selection, as these determine the probability that a mutation becomes fixed. If the variations in kis themselves are entirely random (zero correlation), we might expect both S_C/O and kcatC at the high end of their observed values, as there is nothing to constrain them and the combined effect should have become fixed in some species by positive selection. The TB hypothesis attempts to explain this absence of both high S_C/O and high kcatC by positive selection processes occurring within particular constraints (Figures 4A, 6) imposed on the chemical reaction steps (Tcherkez et al., 2006), but it may also be explained by competing selection pressures. The essential difference is that the origin of the evolutionary constraints is shifted from kis to Ω.

From a biophysical perspective, thermodynamic stability is recognized as the most important constraint on the evolution of proteins and their ability to acquire new function (Tokuriki and Tawfik, 2009; Sikosek and Chan, 2014). The necessity of a protein to maintain the integrity of its folded structure despite the destabilizing effects of accumulated mutations results in only a small percentage being fixed by positive selection. Consequently, in the evolution of C₃ to C₄ plants, destabilizing mutations that are selected on the basis of improved activity are followed by mutations that restore stability with little impact on activity (Studer et al., 2014). This leads to an apparent tradeoff between activity and stability that may well limit the ability of RuBisCO to fix the number of mutations required to increase both S_C/O and kcatC. Depending on the sub-cellular CO₂/O₂ ratio, the fixed mutations increase specificity (for low ratio) or catalytic rate (for high ratio), or a varying combination of both, whichever best optimizes photosynthesis.

The origin of the constraint(s) has significant implications for the optimization of RuBisCO activity. If the constraint is on Ω (i.e., from competing selection pressures) rather than kis, greater variability may be exhibited. To what extent the functional limits of RuBisCO are reflected in the minimum and maximum values of kinetic parameters is not yet clear for RuBisCOs with higher kcatC because of the absence of empirical data. Much effort has been directed toward research on higher plants with particular emphasis on the evolution of C₃ to C₄ plants with their associated CCMs, although the recent work on diatoms may now help stimulate investigations into a more diverse range of photosynthetic organisms (Hanson, 2016; Young et al., 2016). Diatoms and C₃ plants share very similar kcatC, although the variance, var(S_C/O), for diatoms is relatively large (with corresponding variations in CCM expression), whereas for C₃ plants var(S_C/O) is barely significant (Figure 6, Table 1). This could raise the possibility of improving specificity, if not kcatC, in higher plants. It is perhaps not surprising that the non-green (red) algae, from which diatoms have evolved with somewhat lower kcatC values, also exhibit high var(S_C/O) (Figure 6). The data distributions are incomplete (Figures 2A, 6, Table 1); there is a scarcity of data for green algae, photosynthetic bacteria and cyanobacteria, with kcatC values between 6 s⁻¹ and 14 s⁻¹. Discoveries of significant variance among these also may provide important clues on how to achieve increases in both kcatC and S_C/O in higher plants.

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.