Incorporating photosynthetic acclimation improves stomatal optimisation models

Stomatal opening in plant leaves is regulated through a balance of carbon and water exchange under different environmental conditions. Accurate estimation of stomatal regulation is crucial for understanding how plants respond to changing environmental conditions, particularly under climate change. A new generation of optimality-based modelling schemes determines instantaneous stomatal responses from a balance of trade-offs between carbon gains and hydraulic costs, but most such schemes do not account for biochemical acclimation in response to drought. Here, we compare the performance of seven instantaneous stomatal optimisation models with and without accounting for photosynthetic acclimation. Using experimental data from 38 plant species, we found that accounting for photosynthetic acclimation improves the prediction of carbon assimilation in a majority of the tested models. Non-stomatal mechanisms contributed significantly to the reduction of photosynthesis under drought conditions in all tested models. Drought effects on photosynthesis could not accurately be explained by the hydraulic impairment functions embedded in the stomatal models alone, indicating that photosynthetic acclimation must be considered to improve estimates of carbon assimilation during drought. Summary Statement Accounting for photosynthetic acclimation improves the predictions of carbon assimilation in all the stomatal optimization models evaluated. The influence of drought on photosynthesis cannot be fully explained by the hydraulic impairment function of the stomatal models alone.


Introduction
Accurate estimation of net photosynthetic rates (Anet) in terrestrial ecosystems is crucial for understanding and predicting vegetation dynamics, global carbon budgets, and climate change trends (Arora et al., 2020;Fisher & Koven, 2020;Green et al., 2019;Humphrey et al., 2021;Lombardozzi et al., 2015;Mercado et al., 2018). However, photosynthesis modules within current Ecosystem and Earth System models often lack accuracy and provide biased estimates of carbon assimilation (Prentice et al., 2015;Rogers et al., 2017;Seiler et al., 2022). This is explained, among other factors, by the inability of large-scale models to adequately account for the effects of drought on the regulation of stomatal gas exchange and leaf photosynthetic capacity (Zhou et al., 2013).
Novel hydraulically explicit stomatal optimization models that mechanistically incorporate drought impacts have shown promise in improving the accuracy of instantaneous gas exchange predictions, and therefore Anet, in terrestrial models (Eller et al., 2020;Sabot et al., 2020), while reducing their reliance on poorly defined empirical parameters. Previous studies have thoroughly evaluated the performance of a variety of such models (Sabot et al., 2022b;Wang et al., 2020), and have equivocally advocated the integration of non-stomatal limitations alongside stomatal limitations to photosynthesis, particularly, drought impacts on photosynthetic capacity. Medrano, 2002;Martin-StPaul et al., 2012;Peguero-Pina et al., 2009;Salmon et al., 2020;Zhou et al., 2013). This could be an adaptation to minimise the investment of resources in unprofitable processes and instead diverting them towards growth and development to maximize fitness (Bloom, 1986, Anon, 2005. Downregulation of photosynthetic capacity has also been proposed to be associated with photochemical inhibition by elevated sugar concentrations in mesophyll cells (Hölttä et al., 2017), as well as increased abscisic acid accumulation (ABA) and electrolyte concentrations, inhibiting stromal enzymes (Kaiser, 1987). Although the specific proximate mechanisms governing drought-driven photosynthetic downregulation are not fully understood, it is plausible that its ultimate cause is the cost associated with maintaining photosynthetic capacity. Therefore, exploring its implications for stomatal control provides a basis for improving the accuracy of photosynthetic predictions (Sabot et al., 2022b;Smith et al., 2019;Yang et al., 2019).
Stomatal optimisation models rest on an optimal balance of the trade-offs between the carbon gained by photosynthesis and the risks associated with water lost to transpiration (Cowan & Farquhar, 1977). Inspired by economic thinking, such trade-off optimisation is typically expressed as a carbon profit equation (i.e., carbon gain minus carbon cost), where water loss is substituted for an equivalent carbon-based 'penalty' cost. Thus, plants may take the opportunity for higher carbon gain by allowing more transpiration when water is cheap to use and forego carbon gain when it is costly. A major appeal of optimality-based stomatal models is that, unlike their empirical function (Wang et al., 2020). Following Cowan and Farquhar (1977), stomatal optimisation models originally considered a parametric marginal carbon cost of water use (Katul, Palmroth & Oren, 2009;Konrad, Roth-Nebelsick & Grein, 2008;Mäkelä, 1996;Manzoni et al., 2011;Medlyn et al., 2011). Yet, it was never specified how these costs should vary with the environment, time, and species (Manzoni et al., 2013(Manzoni et al., , 2011Schymanski et al., 2007;Wong, Cowan & Farquhar, 1985). To address this issue, subsequent stomatal optimization models have conceptualised the cost of water use as an impairment to hydraulic function arising from the cavitation of xylem vascular tissues Buckley, Frehner & Bailey, 2023;Eller et al., 2020;Joshi et al., 2022;Lu et al., 2020;Sperry et al., 2017;Wang et al., 2020;Wolf, Anderegg & Pacala, 2016).
Motivated by calls to include non-stomatal limitations to photosynthesis in terrestrial models (Sun et al., 2014;Zhou et al., 2013), some studies have attempted to account for non-stomatal limitations within stomatal optimization criteria (Dewar et al., 2018;Hölttä et al., 2017). These studies attempted to capture gas exchange responses to variation in soil water availability by maximizing photosynthesis under a prescribed linear reduction in instantaneous photosynthetic capacity or mesophyll conductance with decreasing xylem pressure (Dewar et al., 2018) or decreasing turgor pressure (Hölttä et al., 2017). However, this assumption of linear downregulation is little supported by evidence (Wang et al., 2020). Other studies (Drake, 2017;Egea, Verhoef & Vidale, 2011;Keenan et al., 2009;Knauer et al., 2019;Yang et al., 2019;Zhou et al., 2014Zhou et al., , 2013 have relied on calibrated functions for prescribing the reduction in maximum photosynthetic capacity or of the maximum mesophyll conductance. By contrast with the aforementioned approaches, Joshi et al. (2022) suggested to add metabolic impairment to hydraulic impairment within the stomatal optimization procedure, by accounting for a marginal maintenance cost of photosynthetic capacity (α). How α changes across space and through time (e.g., with environmental growing conditions, with changes in other functional traits) is unclear, but metabolic impairments are usually correctly represented in models when they are allowed to vary over timescales of a few days (Caldararu et al., 2020;Haverd 92 et al., 2018;Jiang et al., 2020;Sabot et al., 2022a), whereas stomatal regulation is presumed instantaneous.
In this study, we investigated the effects of including an optimality-based representation of photosynthetic acclimation within six existing models and a new variant of one of them. To that end, we extended their optimality criteria to consider a marginal maintenance cost of photosynthetic capacity α, as proposed in Joshi et al. (2022). We tested the performance of these modified models against observational data stemming from dry-down experiments on >35 C3 species. Our main objective was to assess the ability of the extended models to predict Anet and gs, compared to the original formulations (Objective 1). We also analysed the extent to which each model's estimates of Anet and gs could be explained by hydraulic impairment vs. non-stomatal limitations manifested through the regulation of photosynthetic capacity (Objective 2). Finally, we explored the relationship between the cost parameter α, hydraulic parameters, (sub)species-specific functional traits, and their environmental growth conditions (Objective 3). We expected that accounting for photosynthetic acclimation will improve Anet estimates for all models. Based on previous findings by Wang et al. (2020) and Sabot et al. (2022b), we also expected that hydraulic impairment will affect assimilation to a greater extent than non-stomatal limitations during soil dry-down, albeit to varying degrees among models. Finally, we hypothesized that species characterized by high hydraulic vulnerability (e.g., high P50) or high specific leaf areas would have higher photosynthetic maintenance costs, and that "fast growing but drought sensitive" species would also have higher metabolic sensitivity to drought.

Stomatal optimization models
Hydraulics-enabled stomatal optimization models couple water transport and carbon uptake through gs (see Notes 1 in the Supplementary Material). In these steady-state models, gs links to leaf water potential (ψl) through a mass-balance of water, i.e., root-to-leaf water supply equals transpiration demand (Eq. S5), whilst controlling the rate at which CO2 diffuses within the leaf, and thus, Anet (Eq. S6). A given model's optimality criterion works to adjust the instantaneous gs (consequently, ψl) to maximize profit, i.e., Anet minus a carbon-equivalent hydraulic cost Θ (Eq. 1). Stomatal optimization models can be algebraically rearranged such that they differ only in the mathematical form of the Θ function (see Sabot et al., 2022b andWang et al., 2020 for comprehensive reviews).
The PMAX model was proposed by Sperry et al. (2017) with the aim to express the stomatal cost of water loss based on plant hydraulic theory at a minimal parameter expense. The PMAX model balances a carbon gain function, where Anet is normalised by its maximum potential value (Amax) specific to the instantaneous environmental conditions, with a hydraulic risk function, also normalised, varying with ψl and soil water potential (ψs). A unique pairing of gs and ψl hence maximises the difference between the normalised assimilation and risk functions. This model can also be expressed as in Equation 1 by multiplying both the gain function and hydraulic risk function by Amax, in which case Θ becomes where kψ is the xylem hydraulic conductance at a given water potential (ψ) calculated using Eq. S3 (Notes 1 in Supplementary Material), with the double-subscripts s, l, and crit representing soil, leaf, and critical water potentials, respectively. The hydraulic conductance at the critical water potential (kψcrit) represents the point at which the plant pays the maximum cost of canopy desiccation (Sperry et al. 2017). Here, kψcrit is assumed to be zero for simplicity and generalisation across species.
The PMAX2 model was developed by Wang et al. (2020) in order to meet a specific set of mathematical and biological criteria for a biologically sensible Θ. The hydraulic impairment function of this model is based on penalising the instantaneous transpiration (E(ψl)) by its proximity to Ecrit, i.e., E(ψl) at the critical leaf water potential at which we assume hydraulic conductance is null (ψl → -∞). E(ψl) is weighted by Anet in proportion to Ecrit such that Here, the hydraulic impairment function is equivalent to the ratio of two integrals of the hydraulic vulnerability curve function (P(ψ)). The integral in the numerator ranges from ψs to ψl, and the improper integral in the denominator ranges from ψs to minus infinity (at kψcrit). P(ψ) is calculated using Eq. S4.
The PMAX3 model proposed in this study is a biologically more consistent variant of the PMAX2 model. The model, instead of maximising the differential between Anet and the cost of hydraulic impairment, minimises the opportunity cost of carbon and the hydraulic cost. The model uses a similar hydraulic impairment function as PMAX2. Expressed as in eq. 1, The model assumes a carbon cost (ϖ; μmol m -2 s -1 ) to offset carbon gain, with the former representing an investment into recovering the hydraulic conductance lost to embolized xylem (Sabot et al., 2022b). While the CGAIN model can account for the lagged effects of xylem recovery, as in Sabot et al. (2022b), we assume no lagged costs of xylem recovery. This simplified version of CGAIN is expressed as where Kmax is the maximum whole-plant xylem hydraulic conductance between the roots to the leaves.
The PHYDRO model (Joshi et al. 2022) took a phenomenological approach that relies on an impairment function proportional to the square of the differential between the soil and leaf water potential (Δψ 2 ) where γ (μmol m -2 s -1 MPa -2 ) is a calibrated empirical parameter which determines species-specific hydraulic sensitivity and is potentially related to the isohydricity of a species. This model requires the use of whole-plant vulnerability curve, including the outside-xylem hydraulic pathways in fine roots and leaves.
The SOX stomatal model (Eller et al. 2018) drew on the PMAX model but opted to directly impair Anet by the ratio of the hydraulic vulnerability curve evaluated at ψl to that of the vulnerability curve evaluated at ψs Finally, the SOX2 model proposed by Buckley et al. (2023)  parameter ψ50 (i.e. the water potential at a 50% loss of xylem hydraulic conductance) is tuned by a species-specific stomatal behaviour parameter, δ (dimensionless). The parameter b (dimensionless) is the steepness of the vulnerability curve Acclimation of the photosynthetic capacity To calculate Anet, maximum photosynthetic capacity must be known. Following the photosynthetic model of Farquhar et al. (1980), Anet is determined by the lesser of two biochemical assimilation rates (Aj and Ac) that are functions of the available irradiation (Iabs), leaf intercellular CO2 concentration (Ci), and temperature (Eq. S7 and Eq. S8), minus a day leaf respiration rate. The two rates Aj and Ac depend on the maximum rate of electron transport (Jmax) and maximum rate of carboxylation capacity (Vcmax), respectively, with both Jmax and Vcmax differing in their sensitivities to temperature. The parameter values that set Jmax and Vcmax are typically standardized to 25°C (Jmax,25 and Vcmax,25; Eqs. S14 ad S15) and taken to be species, ecosystem, or functional-type averages, without accounting for variation in their standardized magnitudes (i.e., the ratio of Jmax,25 to Vcmax,25 is a fixed value). However, Jmax and Vcmax have been observed to co-vary in the medium to long term, leading to a coordination among the two photosynthetic rates (Chen et al., 1993).
Assuming coordination, we can establish a single photosynthetic capacity cost, based on only one of the two limiting biochemical photosynthetic rates. Specifically, we let both Jmax and Vcmax acclimate by extending the stomatal optimisation criterion (Eq. 1) to account for the cost of maintaining electron-transport capacity (Joshi et al. 2022).
where α is a marginal maintenance cost per unit of electron-transport capacity under standard conditions (i.e. 25 ºC), and αJmax,25 represents the total carbon cost required to maintain photosynthetic capacity. Here, unlike in Joshi et al. (2022), we optimise Jmax,25 instead of Jmax, to limit the effects of temperature on α.
The optimisation problem presented in Eq. 9 allows us to estimate the values of ψl and Jmax,25 that maximise profit, similar to how Eq. 1 allows us to regulate ψl to maximise Anet. The optimisation of Jmax,25 is performed instantaneously, but using the average environmental conditions over the past seven days to capture the weekly acclimation timescale. We obtain the corresponding coordinated Vcmax,25 via the photosynthetic coordination hypothesis (Eq. S10 -Eq. S13). The optimised acclimated parameter values of Jmax,25 and Vcmax,25 are then adjusted to the instantaneous temperature using an Arrhenius function (Eqs. S14 and S15), and fed back as fixed parameters into the biochemical photosynthetic model for the calculation of instantaneous gs, ψs, and Anet via Eq. 1.
In our approach, photosynthetic capacity is not only acclimated to the prevailing atmospheric conditions of the previous seven-day period, but also to ψs through the averaged hydraulic cost or impairment factor. However, as soil moisture conditions were not available for all seven day periods, we assumed that ψs decreased linearly between the first day of the experiment and the last day. We also assumed that the environmental conditions were constant until the start of the experiment.

Experimental data
We used published data from a series of dry-down experiments to evaluate the performance of the seven stomatal optimization models with or without accounting for the photosynthetic acclimation.

Parameter calibrations and simulation experiments
We conducted a different set of parameter calibrations for each of the three objectives of this study ( Fig. 1). Our first set of calibrations was designed to explore the co-variation of α with hydraulic (or hydraulic cost) parameters in different models, as well as its relation to other functional traits and environmental conditions (Objective 3). To obtain the best estimates of α, we fitted all model parameters that influence the optimality criteria. By contrast, our second set of parameter calibrations was much more conservative, as it aimed to infer the models' performance gain from photosynthetic acclimation in the absence of a species-specific knowledge of α (Objective 1). To that end, we fitted one hydraulic (or hydraulic cost) parameter per model while keeping all other parameters fixed. Our third set of parameter calibrations was used to gain insight into the relative contributions of hydraulic (stomatal) vs. non-stomatal limitations on photosynthesis (Objective 2).
In this case, we fitted two model parameters: one hydraulic (or hydraulic cost) parameter and α. calibrating the Kmax, of these models, their specific hydraulic parameter was calibrated, and the Kmax parameter calibrated using the PMAX model was used. Parameter optimization was achieved by minimizing the error between predicted and measured instantaneous gs, Anet and Ci:Ca ratios using Eq. 10. The parameters box shows the parameters calibrated in each set of calibrations.
In our first set of calibrations ("complete calibration"), we accounted for photosynthetic acclimation as described in the 'Acclimation of the photosynthetic capacity' section (Eq. 9). Using the Nelder-Mead algorithm, we looked to fit the best {Kmax, ψ50, b, α} parameter set for each individual model, species, and dry-down experiment that minimized the normalised root mean square error (NRMSE) between predicted and measured instantaneous rates of gs, Anet and Ci:Ca ratios. Compared to the PMAX, PMAX2, PMAX3, and SOX models which only required calibrating the aforementioned parameters, the SOX2, PHYDRO, and CGAIN models required an additional hydraulic cost parameter to be calibrated (i.e., δ, γ, and ϖ, respectively).
The first set of calibrations was performed on the entire dataset of each species, and it was also repeated without photosynthetic acclimation (Eq. 1). Calibrating four or five parameters on small data subsets (e.g., N=4; Table 1) likely leads to unrealistic predictions of these parameters.
Therefore, the parameter estimates from small data samples were given less weight in subsequent statistical analyses.
where χ is the ratio of Ci:Ca, the prime symbol ' denotes predicted values, the double-subscript i represents the position of the paired predicted and actual values, and N is the total number of actual observations for a given species in an experiment. Our last set of parameter calibrations ("calibrated α acclimated") prescribed the hydraulic parameter values obtained in the second set of calibrations, and α was calibrated on the 80 th wet percentile data of the species and dry-down experiment using the same parameter fitting criterion given by Eq. 10.
Both acclimated (Eq. 9) and non-acclimated (Eq. 1) subsequent model simulations were performed mostly outside their calibration samples, over the entire dry-down experiments data for a given species-experiment. We opted to use the hydraulic (or hydraulic cost) parameters obtained for the non-acclimated models when running their acclimated counterparts, in order to make the simulation outcomes as comparable as possible. In all calibrations and model simulations without acclimation, Jmax,25 and Vcmax,25 were obtained using the average values of the actual gs and Anet measurements for the 80 th wet percentile of each species-study ( Fig. S1-S2), and assuming photosynthetic coordination via Eqs. S6, S11, S12, and S13. Finally, to ensure that using a different parameterisation of hydraulic vulnerability did not give the PHYDRO model an unfair advantage or disadvantage compared to other models, we repeated the second and third set of calibrations (and associated simulations) using outside-xylem vulnerability curve parameters (i.e., ψox50 = ψ50/3 and b=1) for other models and xylem vulnerability curve parameters (i.e., ψ50 and b) for PHYDRO. This incidentally provided us with a test of the sensitivity of the models to the choice of vulnerability curve parameters.

Statistical analyses
The performance of the seven stomatal models in predicting Anet with and without acclimation For each statistical metric, we used linear mixed-effects models ( experiments and species were used as crossed random intercepts to account for potential systematic errors in the experimental design or potential species-specific biases that could be introduced across stomatal models. All performance analyses were repeated for gs. We used model with or without acclimation parametrised according to our third set of calibrations, to evaluate how each stomatal optimisation model divides its predictive capacity into the components of hydraulic impairment vs. limitations due to photosynthetic acclimation (Objective 2). For each stomatal model, we analysed how the explained variability (R 2 ) of Anet was partitioned between the hydraulic penalty (stomatal partition), the cost of maintaining photosynthetic capacity (αJmax,25; non-stomatal partition), and species. We accounted for the explained variability by species because species-level ψ50 and b parameters may be inadequate for a given population within the species. The stomatal variance partition was computed as the marginal R 2 of the actual Anet explained by the estimated Anet without acclimation and adjusting random intercepts for each species in an LMM. We then applied this LMM to the estimated Anet values with acclimation, and we obtained the non-stomatal variance partition as the difference between its marginal R 2 and the previously calculated marginal R 2 for the stomatal variance partition. The species variance partition was obtained as the difference between the conditional and the marginal R 2 of the LMM for the acclimated Anet simulations. To understand whether the importance of these variance components changes as the soil dries, we repeated the analysis on the driest 50 th percentile of data for each species within each dry-down experiment.
Finally, we analysed the relationships between the calibrated α from the first set of calibrations and other species-and experiment-specific traits and growing conditions. The traits considered as possible explanatory variables for α were the species-level Hmax, ψ50 , SLA, Hv, and KL collected or imputed from the literature (Table 1). We selected these traits as they stand for different axes of variation representing different water use strategies and life forms. We also included an estimate of the actual Jmax,25 under well-watered conditions (80 th wet percentile; JmaxWW,25), which was inferred from the observed Anet and gs using the one-point method (Eqs. S6, S11, S12, and S13). The environmental conditions of interest were average Iabs and Ca during the experimental period.
During the experiments, Iabs was strongly positively related to VPD and air temperature (Fig. S3).
All explanatory variables that presented a non-normal distribution (Hmax,SLA,Hv,KL,and JmaxWW,25) were transformed using the natural logarithm to minimize the leverage effect of large values. We performed a multivariate LMM in which the variability of α was explained by the chosen species-level traits, ψ50, JmaxWW, and the growing conditions, with the stomatal models and dry-down experiments used as crossed random intercepts. Because species-level traits were taken as fixed predictors of the LMM, it was important we used the stomatal models as random intercepts to account for the variability in α caused by each stomatal model for each species. Next, we identified the most significant predictors using an Akaike Information Criterion (AIC) stepwise algorithm.
Finally, we calculated the variance inflation factor (VIF) of each predictor to explore possible multicollinearity across explanatory variables. We also repeated the above model and the identification of the most significant predictors using the product of α and JmaxWW,25, which corresponds to a photosynthetic capacity cost under well-watered conditions, as an explanatory variable.
In all LMMs, we used the natural logarithm of each species' sample size as a weighting factor to minimize the effect of a small number of data points (see N, Table 1).

Code
All analyses and calibrations were performed using R 4.2.0 (R Core Team, 2022) and the code is freely available from GitHub at https://github.com/vflo/Acclimated_gs_optimization_models. The following specific R packages were used:  dEoptim (Mullen et al., 2011) to solve the instantaneous optimality criteria presented in Eqs.

Results
All three calibration sets were successful for all the species and stomatal models. Most of the model parameterisations obtained from the first set of calibrations were insensitive to acclimation (Figs. S4a, S4b and S4c), with the exception of the magnitude of -ψ50 and Kmax from SOX and SOX2 which were reduced by acclimation. In the second set of calibrations, where Kmax was the only calibrated parameter for PMAX, PMAX2, PMAX3 and SOX, PMAX and SOX's Kmax were much lower than those from the first set of calibrations in the absence of acclimation (Fig. S4a), and moderately lower for PMAX2 and PMAX3 (p<0.05). For PMAX, Kmax also substantially decreased between the first and second sets of parameter calibrations accounting for acclimation (Fig. S4a).
Model performance with vs. without photosynthetic acclimation All seven stomatal models generated highly and positively correlated in-sample Anet estimates across most species (Fig. 2a, Fig. S5). The mean Pearson's correlation coefficient for the stomatal models without acclimation ranged from 0.761 in PMAX2 to 0.812 in PHYDRO, whilst for the acclimated stomatal models using the average α, it ranged between 0.804 and 0.823 in PMAX2 and PMAX, respectively, and using the calibrated α, it ranged from 0.794 in SOX to 0.817 in PMAX (Fig. 2a). The inclusion of acclimation produced a moderately positive effect on the estimates of PMAX2 and PMAX3 in terms of correlation, whereas the effect was negligible for CGAIN, acclimated CGAIN had significantly higher RMSE than PMAX, SOX and SOX2 (Table S2). The calibrated α acclimated PHYDRO presented the lowest average RMSE for Anet at 3.63 μmolC m -2 s -1 (Fig. 2c, Table S2), whereas the average α acclimated CGAIN presented the highest average RMSE at 5.61 μmolC m -2 s -1 .
All the calibrated α acclimated models presented biases significantly lower than the non-acclimated models and not significantly different to zero, with the exception of the calibrated α acclimated PMAX2 and PMAX3 for which biases were slightly lagger than zero (p<0.05). PHYDRO, CGAIN, PMAX and PMAX2 did not show significant differences in biases between their average α acclimated and their non-acclimated model versions (Fig. 2d). All the non-acclimated stomatal models presented a significant average overestimation (bias higher than zero; p<0.05), with PMAX2 presenting the larger bias (0.486) and PHYDRO the lowest one (0.388; Fig. 2d). Unlike for Anet, the models' performance for gs in terms of Pearson's correlation, m, RMSE, and bias was generally not affected by the inclusion of acclimation (Fig. S7). Yet, acclimation led to a significant improvement of proportionality for CGAIN, PMAX2 and PMAX3 (Fig. S7b). However, m was still significantly higher than 1 for PMAX2 and PMAX3, pointing to these models' failure to produce estimates of gs in good proportionality with the observations (Fig. S7b).
Importantly, all models except the PMAX variety were highly sensitive to different vulnerability curve parameters (Fig. S6), poor predictions result when parameters other than those prescribed in the respective model designs were used. Thus, using the outside-xylem vulnerability curve parameters, both SOX and SOX2 exhibited enhanced sensitivity to changes in ψs, leading to rapid and strict stomatal closure as soil dried ( Fig. S6b and S6d). CGAIN became unstable, and for several species, the model failed to converge. By contrast, both PMAX2 and PMAX3 showed improved performance (Fig. S6) and similar estimates to PMAX. Likewise, using the xylem-level vulnerability curve parameters, PHYDRO was unable to accurately depict the decline of Anet as the soil dried out (Fig. S6b).

Hydraulic impairment vs. limitations due to photosynthetic acclimation
The proportion of variance in Anet explained by hydraulic (stomatal limitation) and metabolic impairment (non-stomatal limitation to photosynthesis) differed among models (Fig. 3). Overall, stomatal limitations were markedly more important than non-stomatal limitations for all models across all conditions (Fig. 3a), except for PMAX2 and PMAX3 in which non-stomatal limitations equally mattered. When analysing the driest half of the data, the variance explained by non-stomatal limitations slightly decreased for all models (Fig. 3b). The acclimated PMAX model explained the greatest amount of variability in Anet, with a marginal R 2 of 0.59 for all data points (stomatal + nonstomatal contributions), whereas the acclimated PMAX, PMAX2 and PMAX3 were the best models for the driest half of the data, achieving a marginal R 2 of 0.55. The proportion of variance explained by the different species was greater for the driest half of the data (Fig. 3a), accounting for nearly half of the variance explained by the stomatal and non-stomatal components taken together, which hints at site-specific adaptations or acclimations and at increased errors introduced by species-level parameterisations when water supply is limiting. Relationship between α, functional traits, and environmental growing conditions Calibrating {Kmax,ψ50,b,α} (and additionally,γ,δ or ϖ,for PHYDRO,SOX2,and CGAIN, respectively) for each model, species, and dry-down experiment (first set of parameter calibrations) enabled a majority of the models to match or exceed a R 2 of 0.6 for both gs and Anet (Table S3)  A multivariate LMM combined with a stepwise AIC variable selection (Table 2)  The results were similar albeit explaining a larger part of the variability (R 2 = 0.59) when the relationships between αJmaxWW,25 and the eight potential explanatory variables were explored (Fig.   S8). In that latter case, six explanatory variables were selected: the SLA and maximum height of the species were associated with lower photosynthetic capacity cost under well-watered conditions, whilst Iabs, species KL , Hv, and ψ50 were positively related to αJmaxWW,25 (Fig. S8).  Table 2. Grey scatter open points show all the model-specific contributions to α, across all species. Black plain lines are partial regressions and dashed lines represent the partial regressions plus or minus their standard error. Size of the points represent the weight applied to each data point, which is equal to the natural logarithm of the sample size (Table 1).

Discussion
This study demonstrates that accounting for a maintenance cost for photosynthetic capacity, as proposed by Joshi et al. (2022), improves predictions of Anet (Fig. 2) across a range of hydraulicsenabled stomatal optimality frameworks, particularly under drought conditions (Fig. 2b). Our approach eliminates the need for species-level knowledge of photosynthetic capacity (e.g., Jmax25 and Vcmax25), for empirical functions describing the sensitivity of photosynthetic capacity to drought (Drake, 2017;Egea et al., 2011;Keenan et al., 2009;Knauer et al., 2019;Yang et al., 2019;Zhou et al., 2014Zhou et al., , 2013, or for semi-empirical functions relating photosynthetic capacity to leaf nitrogen investments (Sabot et al. 2022a). Despite these benefits, our approach requires fitting an additional photosynthetic cost parameter α, whose variable responses to species-specific plant traits and current environmental conditions (let alone future conditions) have not been fully elucidated in this study. Nonetheless, our results suggest that the cost parameter α could be predicted from plant functional traits and parameters, in combination with knowledge of the relevant growth conditions. Despite variability in model-specific α values (Fig. 3), the predictions of the seven stomatal models taken together allowed us to meaningfully explore the relationships between α and other plant traits Whether photosynthetic costs change through time or increase during dry-downs is not resolved by our analysis, but it is plausible that the magnitude of α might depend on hydraulic impairment.
Interestingly, we found a positive relationship between α and Iabs which suggests that species experiencing higher light intensities have higher leaf maintenance costs per unit Jmax25. This seemingly contradicts our expectation, as species growing under higher radiation levels are typically characterized by higher Jmax25 (Peng et al., 2021). However, there is also evidence of more pronounced photo-inhibition for leaves under high light intensities exposed to moderate to severe water deficits, compared to leaves under low radiation (Kaiser, 1987;Epron & Dreyer, 1990). This suggests an enhanced photosynthetic sensitivity to drought (i.e., a higher α) for high Jmax25 leaves, perhaps linked to the production costs of photoprotective compounds whose concentrations acclimate to UV radiation levels (Barnes et al., 2023). Further considering the strong positive correlation between Iabs and VPD (Fig S3), it is likely the increased cost associated with elevated radiation is influenced by the confounding effect of higher atmospheric aridity.
Besides its relations to environmental drivers, α was negatively related with species maximum height but positively related with ψ50 and KL (Fig. S8). This suggests a coordination of the photosynthetic cost with the hydraulic safety-efficiency trade-off, as well as with the life form strategy. Such a pattern would agree with the observational evidence that species with higher KL have smaller hydroscapes and close their stomata at less negative water potentials Santiago et al., 2004;Li et al., 2019), hence, in the context of this study, requiring a higher α to prompt a rapid decrease in photosynthetic capacity as the soil dries out. In the same vein, taller species typically show more 'efficient' water use strategies (Flo et al., 2021) and coordination with ψ50 and KL across species and life forms . This coordination implies that taller plants would typically have stronger stomatal regulation as the soil dries out, and therefore a higher photosynthetic cost may be anticipated. However, note that all woody plants included in this study were seedlings, and many of them were far from their maximum potential heights. Nevertheless, the maximum height of the species still successfully represents the axis of variation of the different life form strategies.
With the caveat that performance evaluation of the models presented here only accounts for model predictions of A net and g s , PMAX exhibited the best performance when photosynthetic acclimation was included ( Fig. 2 and 3). Previous studies have shown this stomatal optimisation model to perform as well or better than empirical stomatal models in the absence of photosynthetic acclimation (Sabot et al., 2022b;Venturas et al., 2018;Wang et al., 2020), and it has also been shown to make robust predictions when accounting for it (Sabot et al. 2022a). The SOX, SOX2, PHYDRO, and CGAIN also exhibited reasonably good performance ( Fig. 2 and 3), especially in terms of improvement to proportionality when acclimation was considered (Fig. 2b), but note that these four models were more sensitive to the choice of vulnerability curve parameters than PMAX (cf. Fig. 2 and S6). In stark contrast, even though the inclusion of photosynthetic acclimation improved the performance of PMAX2 and PMAX3 in terms of correlation, RMSE, and bias, it did not improve the proportionality of their Anet predictions (Fig. 2b) except when outside-xylem vulnerability curve parameters were used (Fig. S6). Overall, our simulations' performance without acclimation did not support the assertion made by Wang et al. (2020) that PMAX2 outperforms other stomatal optimisation models, in line with the findings of Sabot et al. (2022b).
Stomatal optimisation accounting for hydraulic impairment functions alone demonstrated enough flexibility to predict Anet reasonably well when using species-level prescribed parameters of vulnerability curves and calibrating other hydraulic parameters (Table S2, "not acclimated"), particularly for those models that incorporate an additional hydraulic parameter besides just Kmax (i.e., CGAIN, PHYDRO, and SOX2). However, all stomatal optimization models examined here except one (PMAX) were sensitive to the choice of hydraulic vulnerability curve parameterisation.
Therefore, they might not perform so well when directly parameterised based on both hydraulic vulnerability curves and estimates of Kmax for a given plant organ (e.g., leaf, branch, stem, root), especially if those estimates are taken in isolation from parameters related to photosynthetic acclimation, which we have shown to positively affect model outcomes. The difficulty lies in determining how the vulnerability of a plant organ should be scaled to better represent the vulnerability of the entire hydraulic system for each model -this is complex because organ-level traits do not accurately characterize whole-plant hydraulic responses (Wang and Frankenberg 2022) -and integrated with our understanding of other plant functions -this is complex because our understanding of trait coordination at the plant level is incomplete . We recommend strong caution and stress the need for in-depth testing before the application of stomatal optimisation models at large spatial scales, let alone before their integration into global vegetation models and land surface models.
European Research Council (787203 REALM) under the European Union's Horizon 2020380 research programme.

Contributions
VF conceptualized the study with the help of all authors. VF and JJ prepared the database. VF and JJ prepared the code. MS and DS supervised the code. VF conducted the calibrations, simulations, and statistical analyses. VF wrote the first draft. All the authors provided scientific input, and reviewed and edited the final manuscript.