Introduction

Quantifying the CO2 diffusion inside leaves of C3 plants is important in both physiological and ecological contexts. Physiologists assess leaf photosynthetic efficiency and capacity, and both of them depend on how CO2 from the atmosphere travel to the chloroplast stroma and how much CO2 released by respiration and photorespiration [“(photo)respired CO2” hereafter] can be refixed by Rubisco (Busch et al. 2013; von Caemmerer 2013). Ecologists often project the impact of global land CO2 fertilization (Sun et al. 2014). The model of Farquhar, von Caemmerer and Berry (1980; “the FvCB model” hereafter), which is widely used as a component for this projection, requires the CO2 level at carboxylation sites of Rubisco (Cc) as its input. The drawdown of Cc, relative to the CO2 level in the ambient air (Ca), depends not only on stomatal conductance for CO2 transfer (gsc) but also on mesophyll conductance (gm), such that (von Caemmerer and Evans 1991):

$${C}_{\mathrm{c}}={C}_{\mathrm{i}}-A/{g}_{\mathrm{m}}$$
(1)

where Ci is the intercellular air space (IAS) CO2 level and A is the net photosynthesis rate.

The FvCB model calculates A as the minimum of the Rubisco activity limited rate (Ac) and electron transport-limited rate (Aj) of photosynthesis, and Sharkey (1985) added a third limitation, accounting for the rate set by triose phosphate utilization (Ap) (see Supplementary Text S1). Equation (1) has been combined with the FvCB model to estimate gm from combined data of gas exchange and chlorophyll fluorescence measurements on photosystem II (PSII) electron transport efficiency Φ2 (Harley et al. 1992; Yin and Struik 2009). The most commonly used method to estimate gm is the ‘variable J method’ (Harley et al. 1992), derived from the Aj part of the FvCB model, using measurements that have to include photorespiratory conditions (Laisk et al. 2006; see Supplementary Text S2). Equation (1) has also been used to estimate gm from online carbon isotope discrimination measurements (e.g. Evans et al. 1994; Tazoe et al. 2011; Barbour et al. 2016a) or oxygen isotope techniques (Barbour et al. 2016b).

Equation (1), as the classical gm model, treats the (photo)respired CO2 in the same way as it treats the CO2 flux that comes from the IAS. Mesophyll resistance (the inverse of mesophyll conductance) consists of components imposed by IAS, cell wall, plasmalemma, cytosol, chloroplast envelope and stroma (Evans et al. 2009; Terashima et al. 2011). Unlike the CO2 from the IAS, the (photo)respired CO2, mainly coming from the mitochondria, does not need to cross the cell wall and plasmalemma, and thus experiences a different resistance. For this reason, Tholen et al. (2012) developed an Equation for the drawdown of Cc, relative to Ci:

$${C}_{\mathrm{c}}={ C}_{\mathrm{i}}-A\left({r}_{\mathrm{w}\mathrm{p}}+{r}_{\mathrm{c}\mathrm{h}}\right)-(F+{R}_{\mathrm{d}}){r}_{\mathrm{c}\mathrm{h}}$$
(2)

where F and Rd are CO2 fluxes from photorespiration and respiration, respectively, rwp is the combined cell wall and plasma membrane resistance, and rch is the chloroplast envelope and stroma resistance (rch). Combining Eqs. (1) and (2) results in gm = 1/[rwp + rch + rch(F + Rd)/A]. Tholen et al. (2012) concluded that mesophyll conductance, as defined by Eq. (1), is influenced by the ratio of (photo)respired CO2 release to net CO2 uptake, (F + Rd)/A, thereby resulting in an apparent sensitivity of mesophyll conductance to [CO2] and [O2]. As this sensitivity does not imply a change in intrinsic diffusion properties, gm as defined by Eq. (1) is an apparent parameter. We shall call it the apparent mesophyll conductance (gm,app). In developing their model, Tholen et al. (2012) assumed a negligible IAS and cytosol resistance, but Eq. (2) still holds if the IAS resistance is lumped into rwp, and part of cytosol resistance is lumped into rwp, and the remaining part is lumped into rch (Berghuijs et al. 2015). If rwp and rch both represent physical resistances, the total mesophyll diffusion resistance (rm,dif) is rwp + rch, and the model of Tholen et al. can be rewritten as

$${g}_{\mathrm{m},\mathrm{a}\mathrm{p}\mathrm{p}}= \frac{1}{{r}_{\mathrm{m},\mathrm{d}\mathrm{i}\mathrm{f}}[1+\omega (F+{R}_{\mathrm{d}})/A]}$$
(2a)

where ω is the fraction of rch in rm,dif.

However, the relative position of mitochondria and chloroplasts is underrepresented in the model of Tholen et al. (2012). Considering six scenarios of the arrangement of these organelles, Yin and Struik (2017) derived the model: gm,app = 1/{rm,dif[1 + ω(1 − λk)(F + Rd)/A]}, where λ is the fraction of mitochondria located in the inner cytosol (i.e. the cytosol area between chloroplasts and vacuole), and k is a factor allowing an increase (k > 1), no change (k = 1), and a decrease (0 ≤ k < 1) in the fraction of inner (photo)respired CO2, caused by gaps when chloroplasts are not continuously aligned. The gaps largely depend on the anatomical parameter Sc/Sm, the ratio of chloroplast area to the mesophyll area exposed to IAS (Sage and Sage 2009). As (1 − λk) is between 0 and 1, the model predicts that the sensitivity of gm,app to (F + Rd)/A is lower than Tholen et al. (2012) initially stated (Yin and Struik 2017). The model of Tholen et al. applies to an extreme case, either where mitochondria are located exclusively in the outer cytosol between plasmalemma and chloroplasts (λ = 0) or where (photo)respired CO2 are completely mixed in cytosol if cytosol resistance is negligible and there are chloroplast gaps (k → 0). In another extreme case where mitochondria are located exclusively in the inner cytosol (λ = 1) and chloroplasts cover completely the cell periphery (k = 1), the model predicts no sensitivity of gm,app to (F + Rd)/A, and Eq. (1) would work well as gm,app becomes gm,dif (= 1/rm,dif). Equation (1) also works when rch is negligible compared to rwp (ω = 0) as if (photo)respired CO2 is released in the same organelle where RuBP carboxylation occurs. Either situation (λk = 1 or ω = 0) can be approximately represented by leaves where mitochondria lies only in the inner cytosol, intimately behind chloroplasts that form a continuum.

Most likely scenarios are somewhere between the two extremes defined by Eqs. (1) and (2), that is, 0 < λk < 1 and 0 < ω < 1. All these scenarios result in different fractions of re-assimilation of (photo)respired CO2 (Yin and Struik 2017), both within mesophyll cells and via IAS (see Supplementary Text S3). It would be useful if ω, λ and k can be measured. One way to derive ω is to use individual resistances that can be calculated from microscopic measurements on leaf anatomy (Evans et al. 1994; Peguero-Pino et al. 2012; Tosen et al. 2012a, b; Tomas et al. 2013; Berghuijs et al. 2015), despite uncertainties in the value of gas diffusion coefficients. Another possible method to estimate ω is to first estimate rwp from oxygen isotope techniques assuming that the outer limit of carbonic anhydrase activity represents the cytosol immediately adjacent to the cell wall (Barbour 2017). Parameter λ can be assessed using electron microscope images for mitochondria distribution (Hatakeyama and Ueno 2016). Most difficult is to measure k, which depends on Sc/Sm. However, whether a high Sc/Sm would make k > 1 or < 1 would depend on the λ value as well as on cytosol resistance, and such a complex relationship is hard to quantify with a simple resistance model. However, because ω, λ and k lump together co-defining the sensitivity of gm,app to (F + Rd)/A, the model of Yin and Struik (2017) can be rewritten to

$${g}_{\mathrm{m},\mathrm{a}\mathrm{p}\mathrm{p}}= \frac{1}{{r}_{\mathrm{m},\mathrm{d}\mathrm{i}\mathrm{f}}[1+m(F+{R}_{\mathrm{d}})/A]}$$
(3)

where m = ω (1 − λk). Although Eq. (3) looks the same as Eq. (2a), their underlying intracellular fluxes for CO2 gradient and re-assimilation differ (see Supplementary Text S3). Equation (3) may be used for estimating m from noninvasive gas exchange measurements where (F + Rd)/A varies.

Many reports (e.g. Flexas et al. 2007a; Vrábl et al. 2009; Yin et al. 2009; Tazoe et al. 2011) showed that gm,app responds to changes in [CO2] or irradiance levels. gm,app was shown in tobacco to increase when [O2] was decreased from 21 to 1% (Tholen et al. 2012). All these responses can be described using a phenomenological equation (Yin et al. 2009). Tholen et al. (2012) explained the O2 response and the commonly observed decline of gm,app with decreasing CO2 below the ambient level (e.g. Flexas et al. 2007a; Vrábl et al. 2009; Yin et al. 2009), based on the earlier introduced sensitivity of gm,app to (F + Rd)/A, because both increasing O2 and decreasing Ci increase (F + Rd)/A. However, the sensitivity of gm,app to (F + Rd)/A cannot explain the observed response of gm,app to irradiances. Moreover, it is unknown whether gm,dif would be conserved across irradiance, CO2 and O2 levels.

In this study, we described a method that explores varying (F + Rd)/A ratios to analyse mesophyll resistance from combined gas exchange and chlorophyll fluorescence measurements. The varying (F + Rd)/A ratios were mainly created using five levels of O2, on two contrasting species tomato and rice. Using these data, we assessed (i) the value of the m factor and whether it differs between species, (ii) whether gm,dif responds to [CO2], irradiance and [O2], and (iii) how the re-assimilation of (photo)respired CO2 is affected by the m factor.

Materials and methods

Experiments and growth conditions

Seeds of tomato and rice were sown, and uniform seedlings were transplanted into pots 2 weeks after sowing, in glasshouse compartments. Pots were filled with soil, and after assessing initial soil nutrient contents, extra nutrients were applied (Table 1). Tomato plants were watered regularly, while rice plants were maintained submerged.

Table 1 Growth and measurement conditions during the experiments with tomato and rice

About 60% of the radiation incident on the glasshouse was transmitted to the plant level. During daytime supplemental light from 600 W HPS Hortilux Schréder lamps (Monster, NL) was automatically switched on when the incident solar flux dropped below a threshold and off when it exceeded a threshold outside glasshouse. These threshold levels were set different for tomato and rice (Table 1), to mimic growth environments of the two species.

Simultaneous gas exchange and chlorophyll fluorescence measurements

We used the Li-Cor-6400XT open gas exchange system with an integrated fluorescence head enclosing a 2-cm2 area (Li-Cor Inc, Lincoln-NE, USA). Young but fully expanded leaves of four replicated plants from staggered sowings were measured for incident irradiance (Iinc) and Ca response curves in each species (Table 1).

Curves were measured at five O2 concentrations (Table 1). Additional light response curves were obtained at 1000 µmol mol−1Ca and 2% O2 to establish nearly nonphotorespiratory conditions for calibration (see later). Gas from a cylinder containing a mixture of O2 and N2 was humidified and supplied via an overflow tube to the air inlet of the Li-Cor where CO2 was blended with the gas, and the IRGA was adjusted for O2 composition of the gas mixture according to the manufacturer’s instructions. Based on pre-test measurements, we used 7–8 min for each step of an A − Iinc curve, and 3–4 min for each step of an A − Ci curve, to reach a steady state. All CO2 exchange data were corrected for CO2 leakage into and out of the leaf cuvette, using measurements on boiled leaves (Flexas et al. 2007b), and then Ci was re-calculated.

When A reached steady state at each light or CO2 step, steady-state fluorescence (\({F}_{\mathrm{s}}\)) was recorded. Maximum fluorescence (\({F}_{\mathrm{m}}^{^{\prime}}\)) was measured using a 0.8 s light pulse of > 8000 µmol m−2 s−1, or the multiphase flash with each phase of 300 ms and ramp depth of 40% (Loriaux et al. 2013). The PSII operating efficiency (\(\Delta F/{F}_{\mathrm{m}}^{^{\prime}}\)) was set as \(({F}_{\mathrm{m}}^{^{\prime}}-{F}_{\mathrm{s}})/{F}_{\mathrm{m}}^{^{\prime}}\) (Genty et al. 1989).

Calibration and pre-determination of Rd and Rubisco parameters

Setting that \({\Phi }_{2}=\Delta F/{F}_{\mathrm{m}}^{^{\prime}}\), Rd was estimated as the negative intercept of a linear regression of A against (IincΦ2/4) using data of A − Iinc curves within the electron transport-limited range for the nonphotorespiratory condition (Yin et al. 2009, 2011). The slope of the regression yields a calibration factor (s), which lumps (1) absorptance by leaf photosynthetic pigments, (2) the factor for excitation partitioning to PSII, (3) basal forms of alternative electron transport, (4) any difference between real efficiency of PSII electron transport (Φ2) and \(\Delta F/{F}_{\mathrm{m}}^{^{\prime}}\), and (5) possibly difference in chloroplast populations sampled by gas exchange and by chlorophyll fluorescence (van der Putten et al. 2018). The electron transport rate J can then be obtained as \(J=s{I}_{\mathrm{i}\mathrm{n}\mathrm{c}}(\Delta F/{F}_{\mathrm{m}}^{^{\prime}})\) (Yin et al. 2009). Like other calibration methods, this procedure assumes that the calibration factor is the same for photorespiratory and nonphotorespiratory conditions, for which photosynthetic rates differ by a factor of (Cc − Γ*)/(Cc + 2 Γ*) (see Eqs. S1.1 and S1.3 in Supplementary Text S1; but with cautions from recent literature, Busch et al. 2018; Tcherkez and Limami 2019).

The parameter Γ* was calculated as 0.5O2/Sc/o, where Sc/o is the relative CO2/O2 specificity of Rubisco (von Caemmerer et al. 1994). Values from in vitro measurements of Cousins et al. (2010) on Sc/o (= 3.022 mbar μbar−1) and Michaelis–Menten coefficients of Rubisco for CO2 (KmC = 291 μbar) and for O2 (KmO = 194 mbar) were taken, assuming that Rubisco kinetic constants are conserved among C3 species. This assumption was checked by in vivo estimates of Sc/o from the lower parts of A − Ci curves of five O2 levels (see “Results”).

Model method

After the above parameters were quantified, we first checked whether gm,dif was variable based on the combined data of gas exchange and chlorophyll fluorescence. Using measured A, Ci and a tentative value for m across its range (0 ≤ m ≤ 1), gm,dif was calculated as

$${g}_{\mathrm{m},\mathrm{d}\mathrm{i}\mathrm{f}}=\frac{A + m\left(F +{ R}_{\mathrm{d}}\right)}{{C}_{\mathrm{i}} -{ C}_{\mathrm{c}}}$$
(4)

where F and Cc can be solved from the Aj equation of the FvCB model, see Eq. (S1.6) in Supplementary Text S1 and Eq. (S2.1) in Supplementary Text S2, respectively. Equation (4) was derived by Yin and Struik (2017, see their Eq. 19), in analogy to the variable J method of Harley et al. (1992; also see Eq. S2.2 in Supplementary Text S2).

The obtained gm,dif responded to a change in both Ci and irradiance (see “Results”). Explaining these responses would need a separate study; to estimate m, here we adopted the generic phenomenological equation of Yin et al. (2009) to describe this response:

$${g}_{\mathrm{m},\mathrm{d}\mathrm{i}\mathrm{f}}={g}_{\mathrm{m}\mathrm{o},\mathrm{d}\mathrm{i}\mathrm{f}}+\delta (A+{R}_{\mathrm{d}})/({C}_{\mathrm{c}}-{\Gamma }_{*})$$
(5)

where gmo,dif and δ are parameters. If δ = 0, Eq. (5) becomes a constant gm,dif mode (= gmo,dif). Any nonzero δ would predict a variable gm,dif in response to CO2, O2 and irradiance levels, and if gmo,dif = 0, parameter δ, as discussed later, represents the carboxylation: mesophyll resistance ratio. Equation (5) was combined with the FvCB and other equations to solve for A (Supplementary Text S1, where reasons for using Eq. 5 are also explained). This results in an equation expressing A as a function of Ci and other variables:

$$A=(-b\pm \sqrt{{b}^{2}-4ac})/(2a)$$
(6)

where

$$a={x}_{2}+{\Gamma }_{*}\left(1-m\right)+\delta ({C}_{\mathrm{i}}+{x}_{2})$$
$$b=m\left({R}_{\mathrm{d}}{x}_{2}+{\Gamma }_{*}{x}_{1}\right)-\left[{x}_{2}+{\Gamma }_{*}\left(1-m\right)\right]\left({x}_{1}-{R}_{\mathrm{d}}\right)-\left({C}_{\mathrm{i}}+{x}_{2}\right)\left[{g}_{\mathrm{m}\mathrm{o},\mathrm{d}\mathrm{i}\mathrm{f}}\left({x}_{2}+{\Gamma }_{*}\right)+\delta \left({x}_{1}-{R}_{\mathrm{d}}\right)\right]-\delta [{x}_{1}\left({C}_{\mathrm{i}}-{\Gamma }_{*}\right)-{R}_{\mathrm{d}}\left({C}_{\mathrm{i}}+{x}_{2}\right)]$$
$$c=-m\left({R}_{\mathrm{d}}{x}_{2}+{\Gamma }_{*}{x}_{1}\right)\left({x}_{1}-{R}_{\mathrm{d}}\right)+\left[{g}_{\mathrm{m}\mathrm{o},\mathrm{d}\mathrm{i}\mathrm{f}}\left({x}_{2}+{\Gamma }_{*}\right)+\delta \left({x}_{1}-{R}_{\mathrm{d}}\right)\right][{x}_{1}\left({C}_{\mathrm{i}}-{\Gamma }_{*}\right)-{R}_{\mathrm{d}}\left({C}_{\mathrm{i}}+{x}_{2}\right)]$$

where x1 = Vcmax (maximum carboxylation activity of Rubisco) and x2 = KmC(1 + O2/KmO) for the Ac-limited conditions; x1 = J/4 and x2 = 2Γ* for the Aj-limited conditions, and for the Ap-limited conditions: x1 = 3Tp (where Tp is the rate of triose phosphate export from the chloroplast) and x2 = − (1 + 3α) Γ* (where α is the fraction of glycolate carbon not returned to the chloroplast).

We found that the \(\sqrt{{b}^{2}-4ac}\) term of Eq. (6) should always take the – sign for either Ac- or Aj-limited rate, but the solution for Ap is mathematically complicated if α > 0 (see Supplementary Text S4). Our data showed that A often declined with increasing Ci within high Ci ranges (see “Results”), suggesting the limitation by triose phosphate utilization with α > 0 (Harley and Sharkey 1991). We conducted sensitivity analyses to choose a value of α although metabolic flux data (Abadie et al. 2018) suggest that its value might be small. We then used Eq. (6) to estimate four parameters: m (0 ≤ m ≤ 1), δ, Vcmax and Tp, by a nonlinear fitting to all data of A − Ci and A − Iinc curves of the five O2 levels (gmo,dif was set to zero, see “Results”). For that, J, as defined earlier as\(s{I}_{\mathrm{i}\mathrm{n}\mathrm{c}}(\Delta F/{F}_{\mathrm{m}}^{^{\prime}})\), were used as input. Our method assumed that Rd does not vary with [O2], and was based on the expectation that neither Vcmax nor Tp varies with [O2], as confirmed experimentally for Vcmax (von Caemmerer et al. 1994). The fitting minimizes the sum of squared differences between estimated and measured A values, using the GAUSS method in PROC NLIN (SAS Institute, NC, USA). SAS scripts can be obtained upon request.

Once A was calculated from Eq. (6), Cc could be solved from Eq. (S2.1) in Supplementary Text S2. Then, gm,dif was re-calculated from Eq. (4) using the estimated m and measured A and Ci, where x1 and x2 terms were chosen according to whether the modelled A was Ac-, Aj- or Ap-limited. This showed gm,dif in response to CO2, irradiance, and O2 levels.

With gm,dif and other parameters, we calculated the fraction of (photo)respired CO2 being refixed (frefix), the fraction of (photo)respired CO2 being refixed within the mesophyll cells (frefix,cell), and the fraction of (photo)respired CO2 being refixed via IAS (frefix,ias), using Eqs. (S3.4), (S3.5) and (S3.6), respectively, in Supplementary Text S3. In these equations, rsc is the stomatal resistance to CO2 diffusion (being 1.6 times measured stomatal resistance to water vapour), and rcx is the carboxylation resistance (which is \(({C}_{\mathrm{c}}+{x}_{2})/{x}_{1}\), von Caemmerer 2000). As discussed in Supplementary Text S3, these calculations need ω and λk as inputs. The estimate for m was 0.3 for tomato and 0.0 for rice (see “Results”). For tomato, we measured ω (0.65) for leaves of the same age in the same cultivar “Growdena” (see Berghuijs et al. 2015) for calculating λk, from m = ω (1 − λk). For rice, λk was set to 1.0 to agree with the estimate that m = 0. In such a case, ω is not needed as Eqs. (S3.4) and (S3.5) become simplified as Eqs. (S3.3) and (S3.9) in Supplementary Text S3, respectively.

Results

Use of the five O2 levels generated diverse shapes of photosynthetic responses to irradiance and CO2 levels (Fig. 1). Our model approach, combined with data for A (Fig. 1) and for \(\Delta F/{F}_{\mathrm{m}}^{^{\prime}}\) (Fig. S1), yielded an estimation of a set of parameters as described below.

Fig. 1
figure 1

Measured (points) and modelled (curves) net CO2 assimilation rate A of tomato (filled circle, solid curves) and rice (open circle, dashed curves) as a function of incident irradiance Iinc (left panels) and of intercellular CO2 concentration Ci (right panels) at different O2 percentages as shown in individual panels. Each point represents the mean of four replicated plants. The A − Iinc curve under nonphotorespiratory (NPR) condition was obtained at 2% O2 combined with ambient CO2 level of 1000 μmol mol−1. Curves were drawn from connecting two nearby values calculated by the model

Estimated Rd and s

Data of A − Iinc curves within the range of Iinc ≤ 200 µmol m−2 s−1 showed that the relationship between A and (IincΦ2/4) was linear for the conditions with a gas mixture of 2% O2 with 1000 μmol mol−1Ca (Fig. 2), where Φ2 was set to be \(\Delta F/{F}_{\mathrm{m}}^{^{\prime}}\). The value of Rd estimated from this linear relationship was 1.2 (standard error or s.e. 0.1) µmol m−2 s−1 for tomato and 1.1 (s.e. 0.1) µmol m−2 s−1 for rice. The slope of the A − (IincΦ2/4) linearity (i.e. calibration factor s) was 0.4570 (s.e. 0.0076) for tomato and 0.5488 (s.e. 0.0076) for rice. Values of s were also re-estimated, together with other parameters, in fitting Eq. (6) to all data; but the re-estimated s remained the same, suggesting that we reached a nonphotorespiratory condition using the gas mixture.

Fig. 2
figure 2

Linear relationship between net CO2 assimilation rate A and IincΦ2/4, where Φ2 is set to be \(\Delta F/{F}_{\mathrm{m}}^{^{\prime}}\) and Iinc is ≤ 200 μmol m−2 s−1 (each point represents the mean of measurements on leaves from four replicated plants), for nonphotorespiratory condition (2% O2 combined with Ca = 1000 μmol mol−1). The intercept of regression lines gives an estimate of −Rd (see Yin et al., 2011), and the slope gives an estimate of the calibration factor s for converting \(\Delta F/{F}_{\mathrm{m}}^{^{\prime}}\) into the linear electron transport rates (see the text)

The first few data points of the A − Ci curves were linear, and gross leaf photosynthesis values A + Rd were plotted versus Ci within this linear range. The intercept of this line with the Ci-axis gives the estimate of the Ci-based CO2 compensation point, commonly noted as Ci*. The value of Ci* increased linearly with increasing O2 levels (Fig. 3). Half of the reciprocal of this linear slope gives an in vivo estimate of Sc/o, which was 2.71 mbar μbar−1 for tomato and 3.13 mbar μbar−1 for rice. Using the method of Yin et al. (2009) gave similar in vivo estimates of Sc/o (results not shown). These values are close to 3.02 mbar μbar−1 measured in vitro for wheat by Cousins et al. (2010), confirming that Sc/o is conserved among C3 species. We will use 3.02 mbar μbar−1 for further analysis (but see sensitivity analysis later).

Fig. 3
figure 3

Values of CO2 compensation point Ci* [identified as the intercept at the Ci-axis of the initial strictly linear part of leaf gross CO2 assimilation rate (A + Rd) versus Ci] plotted as a function of the O2 levels, for tomato and rice leaves

Dependence of gm,dif on CO2 and irradiance level

Equation (4) assuming an electron transport limitation, was applied to check the pattern of gm,dif across a range of Iinc and Ci levels, by setting m either to 0 (equivalent to the variable J method of Harley et al. (1992) for gm,app) or to a value between 0 and 1. A similar response was obtained for various O2 levels, except for 2% O2. At that oxygen concentration, Eq. (4), like the variable J method, cannot be reliably applied due to insufficient photorespiration (see Supplementary Text S2). Although the obtained gm,dif sometimes had unrealistic values largely due to unrealistic values of Cc (as often occurs when using the variable J method, see Yin and Struik 2009), an overall trend of gm,dif in response to Iinc and to Ci was obtained. An example of the response is shown in Fig. 4 for the case of 10% O2 level for tomato. gm,dif increased monotonically with increasing Iinc (Fig. 4a), and decreased gradually with an increase in Ci (Fig. 4b). Changing m did not change the response pattern, but only the absolute value of gm,dif, and a nonzero m resulted in higher gm,dif than the value obtained from setting m = 0 (Fig. 4).

Fig. 4
figure 4

Calculated gm,app using the variable J method of Harley et al. (1992) (open square) or gm,dif using Eq. (4) where parameter m is set to 0.29 (filled circle), as a function of a incident irradiance Iinc or b intercellular CO2 level Ci, under the condition of 10% O2 for tomato leaves. Points were obtained, based on the Aj part of the FvCB model, using measured A and J that was derived from chlorophyll fluorescence with the calibration as described in the text. The monotonically descending curve in panel (b) is drawn from values of the modelled gm,dif using the full FvCB model of three limited rates

Estimates of parameters δ, m, Vcmax and Tp

Equation (6) for describing A was applied to estimate gmo,dif, δ and m, using data of both A − Iinc and A − Ci curves. The obtained gmo,dif did not differ significantly from zero (p > 0.05), which is supported by the result that the calculated gm,dif by Eq. (4) at low Iinc was close to zero (Fig. 4a). Also, model fit became worse if δ was fixed to zero than if gmo,dif was fixed to zero, supporting the variable gm,dif mode. Sensitivity analysis with respect to α suggested that a change within its relevant range had no impact on the estimates of parameters other than Tp (see below). We set gmo,dif to zero, and α to 0.3 (Busch and Sage 2017), in the subsequent analysis.

Equation (6) describes well both A − Iinc and A − Ci curves (Fig. 1), with an overall R2 being > 0.99 for either species (Table 2). Most of the data points (> 80%) were Aj-limited, indicating that chlorophyll fluorescence signals generally echoed gas exchange data since we calculated J from chlorophyll fluorescence measurements as \(s{I}_{\mathrm{i}\mathrm{n}\mathrm{c}}(\Delta F/{F}_{\mathrm{m}}^{^{\prime}})\). Only a few points at low Ci of A − Ci curves or at high Iinc of A − Iinc curves were Ac-limited, and a few points at high Ci of A − Ci curves under low O2 conditions were Ap-limited. The estimated m was ca 0.3 for tomato but was 0.0 for rice (Table 2). The estimated δ was also higher for tomato (1.4) than for rice (1.0) (Table 2). Other parameter values were similar for the two species: 113.7 and 111.0 µmol m−2 s−1 for Vcmax, and 8.3 and 7.8 µmol m−2 s−1 for Tp, for tomato and rice, respectively.

Table 2 Estimates (standard errors in brackets) of two major parameters (δ and m), and Vcmax and Tp, from fitting Eq. (6) to irradiance- and CO2 response curves of five O2 levels for leaves of tomato and rice

Sensitivity analysis

Given that any uncertainty in estimated s and Rd and in other parameters (Sc/o, KmC, KmO and α) may have an impact on the major estimated parameters (m and δ in this study), we carried out sensitivity analyses. The estimation of δ and m was very sensitive to s and Sc/o, and less sensitive to Rd (Fig. S2), but virtually insensitive to KmC, KmO and α (results not shown). Both δ and m decreased monotonically with increasing s (Fig. S2a). The estimate of δ decreased with increasing Sc/o, whereas that of m changed in an opposite direction (Fig. S2b). The obtained response of δ (the parameter in Eq. (5) on mesophyll conductance) to both Sc/o and s is expected in the same way as gm,app responds to these parameters (Harley et al. 1992). The opposite response of m to Sc/o and s is probably because photorespiration, i.e. the F term in Eq. (3), which is relevant to determining m, has an opposite response to Sc/o and s. As Rd has the same effect as the F term has (see Eq. 3), the estimated m decreased with increasing Rd, whereas δ changed in an opposite direction (Fig. S2c). As expected, any sensitivity to KmC and KmO occurred with the estimated Vcmax, whereas a sensitivity to α occurred with Tp (results not shown).

Calculated fractions for re-assimilation of (photo)respired CO2

The calculated fractions of (photo)respired CO2 being refixed, using Eqs. (S3.4–S3.6) in Supplementary Text S3, are shown in Fig. 5, using the result at 21% O2 as the example. The trends were similar for O2 levels above 2%. Except for very low Iinc or Ci levels, the refixed fractions were quite consistent over a wide range of conditions. frefix,cell was lower in tomato (0.25) than in rice (0.49) (Fig. 5), largely due to the fact that the estimated m was 0.3 for tomato but 0.0 for rice (Table 2). In contrast, frefix,ias was higher in tomato than in rice. As a result, the total re-fixation fraction frefix was comparable for the two species, i.e. up to ca 0.6.

Fig. 5
figure 5

Calculated fractions of total re-assimilation (filled circle, frefix), of re-assimilation within mesophyll cells (open square, frefix,cell), and of re-assimilation via the intercellular air spaces (open triangle, frefix,ias) at different incident irradiance (a, c) or intercellular CO2 (b, d) levels, in leaves of tomato (a, b) and rice (c, d), when the O2 level was 21%. The horizontal dashed line represents the calculated frefix,cell using the model predicted A values. In the calculation for tomato, we used the value of ω (the proportion of rch in total rm,dif) of 0.65 that we measured, as reported by Berghuijs et al. (2015, see the text)

Responses of stomatal and mesophyll conductance to O2

Except for a few cases, gsc generally decreased with increasing [O2], and was lower in tomato than in rice (Fig. 6). The calculated value of gm,dif also decreased with increasing [O2], except for very high CO2 conditions which lowered gm,dif to the extent that the O2 response of gm,dif was no longer significant (Fig. 6d,j). gm,dif was higher in tomato than in rice.

Fig. 6
figure 6

Stomatal conductance for CO2 diffusion gsc (open symbols) and mesophyll conductance gm,dif (closed symbols) of tomato (af) and rice (gl) leaves in response to O2 level, at high (left panels), medium (middle panels) and low (right panels) Iinc levels (ac, gi) or Ca levels (df, jl). Values of Iinc or Ca are shown at each corresponding panels, where units of Iinc and Ca are μmol m−2 s−1 and μmol mol−1, respectively

Discussion

Analysing mesophyll resistance

Compared with the CO2 flux coming from IAS, (photo)respired CO2 experiences different resistances. This suggests the need to dissect rm,dif into sub-components. Anatomical measurements can partition mesophyll resistance into individual sub-components (Peguero-Pina et al. 2012; Tosens et al. 2012a, b; Tomas et al. 2013; Carriquí et al. 2019). The calculation of these sub-components relies on many assumed diffusion or permeability coefficients that are uncertain (Berghuijs et al. 2015). Furthermore, this approach does not quantify the effect of the arrangement of mitochondria and chloroplasts on the intracellular CO2 diffusion.

In line with anatomical measurements, Eq. (2) dissects rm,dif into two sub-components rwp and rch (Tholen et al. 2012). Oxygen isotope techniques may estimate rwp based on certain assumptions (Barbour 2017), but so far have been explored to separate rwp and rch within the framework of the classical gm model, Eq. (1) (Barbour et al. 2016b). Equations (1) and (2) both underrepresent the intracellular arrangements of organelles. In contrast, the model of Yin and Struik (2017), Eq. (3), has a factor lumping (i) the rch:rm,dif ratio (ω), (ii) the fraction of (photo)respired CO2 that are released in the inner cytosol (λ), and (iii) k, the factor for the change in λ as a result of the chloroplast gaps. The factor k is particularly hard to assess. Since ω, λ and k lump as such that m = ω (1 − λk), Eq. (3) provides an approach by exploring nondestructive gas exchange and chlorophyll fluorescence measurements under different levels of O2 that created large variations in photorespiration. Instead of estimating individual resistances, our nonlinear fitting approach estimates rm,dif as a whole, as well as the m factor. Common nonlinear procedures, typically by fitting A − Ci curves, estimate four or even more parameters of the FvCB model, such as Vcmax, J, Tp, and gm (e.g. Sharkey et al. 2007). In our method, J was measured from chlorophyll fluorescence. Despite a wide range of O2 levels exploited, we restricted the number of estimated parameters to four from Eq. (6)-based nonlinear fitting. All four were reliably estimated for both species with small standard errors (Table 2).

Our approach is still a simplified representation of complex diffusion pathways. Some respiratory flux may originate in the chloroplasts, in cytosol, and in the heterotrophic tissues such as epidermis, vasculature, and bundle sheath (Tcherkez et al. 2017). These components of Rd could be incorporated as additional terms into the CiCc gradient equation, Eq. (S1.7) in Supplementary Text S1. However, they are ignored here as fractions of these components in Rd are generally unknown. There may also be some activity of phosphoenolpyruvate carboxylase (Vpepc) in cytosol (Douthe et al. 2012; Abadie and Tcherkez 2019), which would counteract the effect of (F + Rd) on gm,app. But our procedure of estimating Rd may have accounted for this, i.e. the estimated Rd represents the net rate of true Rd minus Vpepc. Tholen et al. (2012) showed that small amounts of Vpepc have little impact on gm,app.

Variation of gm,dif with CO2, irradiance and O2 levels

Reports using chlorophyll fluorescence data consistently showed that gm,app initially increases and then decreases with increasing Ci and increases monotonically with increasing Iinc (e.g. Flexas et al. 2007a; Yin et al. 2009). Similar results for gm,app in response to Ci (Vrábl et al. 2009; Tazoe et al. 2011) and to Iinc (Douthe et al. 2012) were sometimes reported, using the carbon isotope discrimination method. No change in anatomical arrangements was observed that could explain the variable gm,app (Carriquí et al. 2019). Gu and Sun (2014) showed that the reported response of gm,app to a change in CO2 or in Iinc may be due to the artefact of errors in experimental measurements. Although resolving experimental uncertainties is urgently needed, consistent variations of gm,app cannot be ascribed only to experimental errors because responses due to random errors would be irregular and inconsistent among various reports. Théroux-Rancourt and Gilbert (2017) demonstrated that changing patterns of light penetration within the leaf 3D-structure leads to different contributions of each cell layer to bulk-leaf mesophyll conductance, resulting in an apparent response of the bulk-leaf gm,app to light intensity. However, their theory cannot explain the response of gm,app to Ci.

Most results using the variable J method of Harley et al. (1992) showed that within a low Ci range, gm,app typically decreases with decreasing Ci (e.g. Flexas et al. 2007a; Vrábl et al. 2009; Yin et al. 2009; Fig. 4b). Tholen et al. (2012) also noted a decrease of gm,app with increasing O2 level. Tholen et al. (2012) explained these responses to Ci and O2 as a consequence of the fact that gm,app as an apparent parameter decreases with an increase in the (F + Rd)/A ratio. If this is the only explanation of variable gm,app, one would expect that gm,dif would be independent of Ci because Eq. (4) for gm,dif already accounts for the (F + Rd)/A ratio. However, gm,dif still declined with decreasing Ci within its low range, albeit to a lesser extent (Fig. 4b). In fact, within the low Ci range, A is limited by Rubisco activity; so, as noted by Yin et al. (2009), using the variable J method or Eq. (4) assuming an electron transport limitation results in an artefactual decline of estimated gm because the occurrence of additional alternative electron transport is wrongly attributed to the mesophyll diffusional limitation. This assertion was supported by the result that the decrease of gm,dif with decreasing Ci within the low Ci range was no longer obtained once gm,dif was calculated from the full FvCB model (Fig. 4b). This suggests that the decrease of gm,app with decreasing Ci is explained more by the occurrence of alternative electron transport than by the theory of Tholen et al. Anyway, the theory does not explain the decreases of gm,app with increasing Ci within its high range, or with decreasing Iinc, or with lowering temperature as reported previously (Bernacchi et al. 2002; Warren and Dreyer 2006; Yamori et al. 2006; Scafaro et al. 2011; Evans and von Caemmerer 2013).

We found that gm,dif increased with Iinc (Fig. 4a). This increase continued within the high Iinc range, where some additional alternative electron transport is also expected, suggesting that the increase of gm,dif with Iinc overrode any artefactual decline caused by alternative electron fluxes. Also, gm,dif decreased with increasing O2 (Fig. 6), in the same direction as the O2 response of gm,app reported by Tholen et al. (2012). Literature on O2 responses of diffusional conductance is scarce (Farquhar and Wong 1984; Buckley et al. 2003). Our data showed that both gsc and gm,dif generally declined with increasing O2. So, gm,dif is variable, in response to Ci, Iinc, and O2, in a similar pattern as gsc responds to these variables (e.g. Morison and Gifford 1983; Farquhar and Wong 1984; Buckley et al. 2003).

The identified variable gm,dif was based on the assumption that m (= ω (1-λk)) is constant, independent of short-term changes (within 3–8 min) in irradiance or [CO2]. This is supported by Carriquí et al. (2019), who reported that anatomical parameters determining ω and k hardly vary with short-term changes in irradiance or [CO2]. Chloroplasts and mitochondria in some plants may move under varying light, but they always colocalize (Islam et al. 2009), suggesting that λ also hardly varies. We are unable to find evidences supporting quantitative changes that m or its components must have to obtain invariable gm,dif with irradiance, [CO2] and [O2].

gm,dif defined here is still a bulk-leaf trait. Like bulk-leaf gm,app, it may not represent intrinsic transport properties. Also, our result on the variable gm,dif is subject to experimental confirmation by other methods. If proven true, future studies are needed to examine if the variable gm,dif can emerge from fluxes and concentrations across the real 3D-structure of leaves, as well as in relation to membrane permeability and other properties. Here we only describe the response from bulk-leaf equations themselves. gm,dif can be formulated from Eq. (4) and A = V − F − Rd (where V is the carboxylation rate) as

$${g}_{\mathrm{m},\mathrm{d}\mathrm{i}\mathrm{f}}=\frac{V-(1-m)(F+{R}_{\mathrm{d}})}{{C}_{\mathrm{i}}(1-{C}_{\mathrm{c}}:{C}_{\mathrm{i}})}$$
(7)

When Iinc increases, only the numerator increases significantly; so Eq. (7) predicts that gm,dif increases with increasing Iinc. If the CO2 gradient from Ci and Cc is regulated such that the Cc:Ci ratio is roughly constant for a given O2 level (results not shown), Eq. (7) also predicts that gm,dif will decrease monotonically with Ci because according to the FvCB model, the V increment per Ci increment decreases with increasing Ci. Finally, the F term increases when O2 increases; as a result, Eq. (7) predicts that gm,dif decreased with increasing O2 (Fig. 6).

Interpretation of the model and estimated parameter values

Our method is based on Eq. (3), the equation summarized by Yin and Struik (2017) from considering six possible scenarios for the intracellular organelle arrangement. Recently, Ubierna et al. (2019) came up with the same model but formulated gm,app in a Michaelis–Menten-like equation, i.e. gm,app = A·gm,dif/[A + m(F + Rd)] (see their Eq. 15; note that gm,app was written as gm in their notations). The maximum value of gm,app is gm,dif, while the Michaelis–Menten constant “Km” is m(F + Rd). For the case of tomato where m = 0.3 and Rd = 1.2, the “Km” occurs at A ≈ 2.0 μmol m−2 s−1 for the ambient O2 condition. This suggests that gm,app and gm,dif only differ significantly when A is low, which our results (Fig. 4) confirmed.

In view of the variation of gm,dif shown in Fig. 4, we adopted Eq. (5), which accommodates either constant or variable gm,dif in relation to Ci, Iinc and O2 levels. Although the equation is phenomenological and has an a priori assumption that gm,dif grows with relative carboxylation and the estimates of its parameters are expectedly sensitive to the pre-input values of s, Sc/o and Rd (Fig. S2), the model generated useful insights.

Our results supported no constant gm,dif, but a variable gm,dif with parameter δ being 1.0 for rice and 1.4 for tomato (Table 2). Equation (5) with gmo,dif = 0 for our variable gm,dif mode can be rewritten to \({r}_{\mathrm{m},\mathrm{d}\mathrm{i}\mathrm{f}}=({C}_{\mathrm{c}}-{\Gamma }_{*})/[\delta \left(A+{R}_{\mathrm{d}}\right)]\). As \(\left(A+{R}_{\mathrm{d}}\right)\) can be calculated from the FvCB model as \(({C}_{\mathrm{c}}-{\Gamma }_{*}){x}_{1}/\left({C}_{\mathrm{c}}+{x}_{2}\right)\), the above equation becomes \({r}_{\mathrm{m},\mathrm{d}\mathrm{i}\mathrm{f}}=({C}_{\mathrm{c}}+{x}_{2})/\left(\delta {x}_{1}\right)\). As \(({C}_{\mathrm{c}}+{x}_{2})/{x}_{1}\) is defined as carboxylation resistance rcx (von Caemmerer 2000), it follows that

$$\delta ={{r}_{\mathrm{c}\mathrm{x}}/r}_{\mathrm{m},\mathrm{d}\mathrm{i}\mathrm{f}.}$$
(8)

Thus, parameter δ of Eq. (5) has a meaning, representing the carboxylation: mesophyll resistance ratio. Our estimates for δ (Table 2) suggest that \({{r}_{\mathrm{c}\mathrm{x}}\, \mathrm{a}\mathrm{n}\mathrm{d}\, r}_{\mathrm{m},\mathrm{d}\mathrm{i}\mathrm{f}}\) had similar values in rice leaves, whereas rcx was ca 40% higher than rm,dif in tomato leaves.

Our estimate of the factor m was ca 0.3 for tomato and 0.0 for rice (Table 2). Thus, using Eq. (1), which is the special case of the generalized model when m = 0, actually suits for rice leaves but does not work for tomato leaves when (F + Rd)/A is high. As stated in Introduction, the classical model works well if mitochondria are located exclusively in the inner cytosol (λ = 1) and chloroplasts cover fully the mesophyll periphery that k = 1. Sage and Sage (2009) and Busch et al. (2013) showed that compared with other species, in rice leaves, there are stromules that effectively extend chloroplast coverage of the cell periphery and mitochondria locate in the cell interior and are intimately associated with chloroplasts/stromules. These features engender such a structure as if (photo)respired CO2 is released in the same compartment where RuBP carboxylation occurs. This is the case when Eq. (1) works well. Therefore, our results with curve-fitting to gas exchange data actually agree with anatomical differences between species.

Such differences are also shown in the fractions of re-fixation of (photo)respired CO2 calculated from resistance components (Fig. 5). With the distinct anatomical feature of rice leaves, (photo)respired CO2, if to exit mesophyll cells, will have to travel via the stroma, thereby maximizing the re-fixation of (photo)respired CO2 within the cell. Therefore, rice had higher values of frefix,cell than tomato (Fig. 5). For a given set of resistance values, the organelle arrangements as in rice leaves that make the highest frefix,cell can result in low frefix,ias (see Supplementary Text S3). Moreover, in line with the observation of Ouyang et al. (2017) on rice ‘IR64′, the cultivar we used, rice had high stomatal conductance, compared with tomato (Fig. 6). A low gsc would make (photo)respired CO2 more difficult to exit into the atmosphere via IAS. This also contributed to higher values of frefix,ias in tomato than in rice (Fig. 5). As a result, the two species had similar values (up to 60%) of the total re-fixation, frefix. The calculated frefix,ias and frefix varied with Iinc or Ci levels (Fig. 5), because resistance components rsc and rcx varied with these variables. The calculated frefix,cell was more constant (Fig. 5). Substituting Eq. (8) into Eq. (S3.5) in Supplementary Text S3 gives

$${f}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{i}\mathrm{x},\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}=\frac{(1-\omega )(\delta +\omega \lambda k)}{\left(1+\delta \right)[\delta -\omega \lambda k(\delta +\omega -1)]}$$
(9)

As all terms are constant, Eq. (9) describes why frefix,cell stayed invariant. Using isotope mass spectrometry and gas exchange measurements, Busch et al. (2013) determined frefix,cell, frefix,ias and frefix, being 0.29, 0.22, and 0.51 for rice under ambient CO2 and high light conditions. Our estimates for rice somewhat differed from their values for comparable conditions (Fig. 5d).

The difference in the value of factor m between the species also has implications on values of Ci* and the relationship between Ci* and Γ*. Ci* was lower in rice than in tomato at a given O2 level (Fig. 3), and Ci* at the lowest O2 in rice was even negative (Fig. 3b). A negative Ci* could be due to measurement noises, uncertainties in assuming constant Rd, and the influence of varying amounts of Vpepc (see earlier discussions). However, for the case where m = 0, it can be seen from Eq. (1) that Ci* = Γ*Rd/gm (von Caemmerer et al. 1994); so, Ci* is always lower than Γ*. This agrees with our linear relation for rice in Fig. 3b (where the term 0.16O2 can be considered as Γ*, given that Γ* = 0.5O2/Sc/o). Setting the intercept of this relation equal to –Rd/gm and knowing that Rd = 1.064 µmol m−2 s−1 (Fig. 2b), gm for rice can be solved as ca 0.1 mol m−2 s−1 bar−1, comparable with its value calculated in the other way for low CO2 conditions (Fig. 6l). So, a negative Ci* for low O2 conditions could actually represent biological realities, i.e. high intracellular re-fixation of both respired CO2 and photorespired CO2 sufficed to (over)compensate for photorespiratory losses. In contrast, for cases where m ≥ 0, the relation between Ci* and Γ* can be formulated from Eq. (S1.7) in Supplementary Text S1 as

$${C}_{\mathrm{i}\mathrm{*}}={ \Gamma }_{\mathrm{*}}-[\left(1-m\right){R}_{\mathrm{d}}-mF]{/g}_{\mathrm{m},\mathrm{d}\mathrm{i}\mathrm{f}}$$
(10)

Equation (10) means that Ci* is no longer necessarily lower than Γ* (see also Tholen et al. 2012), depending on relative values of (1–m)Rd versus mF. Our result in Fig. 3a suggests that Ci* is 2.134 μbar higher than Γ* for tomato. Thus, different m values estimated by curve-fitting for the two species are supported by the Ci* vs Γ* relationships in Fig. 3, suggesting that our approach is internally consistent.