Novel multimodel ensemble approach to evaluate the sole effect of elevated CO2 on winter wheat productivity

Elevated carbon-dioxide concentration [eCO2] is a key climate change factor affecting plant growth and yield. Conventionally, crop modeling work has evaluated the effect of climatic parameters on crop growth, without considering CO2. It is conjectured that a novel multimodal ensemble approach may improve the accuracy of modelled responses to eCO2. To demonstrate the applicability of a multimodel ensemble of crop models to simulation of eCO2, APSIM, CropSyst, DSSAT, EPIC and STICS were calibrated to observed data for crop phenology, biomass and yield. Significant variability in simulated biomass production was shown among the models particularly at dryland sites (44%) compared to the irrigated site (22%). Increased yield was observed for all models with the highest average yield at dryland site by EPIC (49%) and lowest under irrigated conditions (17%) by APSIM and CropSyst. For the ensemble, maximum yield was 45% for the dryland site and a minimum 22% at the irrigated site. We concluded from our study that process-based crop models have variability in the simulation of crop response to [eCO2] with greater difference under water-stressed conditions. We recommend the use of ensembles to improve accuracy in modeled responses to [eCO2].

better growth and yield 14 . Elevated CO 2 (CO 2 fertilization effect) will be beneficial for Ribulose 1,5-bisphosphate carboxylase/oxygenase (Rubisco) and may inhibit photorespiration and increase photosynthesis. Ainsworth and Rogers 15 concluded that [eCO 2 ] stimulated light-saturated photosynthesis by 31% and reduced stomatal conductance by 22% in free-air CO 2 enrichment (FACE) experiments. Meanwhile, Kruijt et al. 16 concluded that at [eCO 2 ], stomatal activity is reduced. This change in stomatal activity resulted in a 50% increase in water use efficiency. Varga et al. 10 reported increased water use efficiency in winter wheat under stress conditions due to [eCO 2 ] (700-1000 ppm). Photosynthesis increases with [eCO 2 ] following a Michaelis-Menten curve. The Michaelis-Menten constant (K m ) could be used to quantify the [eCO 2 ] effects under different temperatures 17 . Almost 23% of the carbon fixed by photosynthesis is lost due to photorespiration and if it is stopped completely, the carboxylation reaction could increase to 53%. It is suggested that with the future rise in [CO 2 ] Rubisco will have higher Km (6.3 to 15 µM) resulting in higher photosynthetic rate and efficiency 15 . Tausz et al. 14 concluded that photosynthetic acclimation or photosynthetic downregulation could be inhibited by [eCO 2 ]. Similarly, [eCO 2 ] resulted in a decrease in evapotranspiration (ET) 18 . Transpiration efficiency could be improved by increased net photosynthesis and reduced stomatal conductance 19 . Radiation use efficiency (RUE) can be defined as biomass produced per unit of light energy used by crops. Different process-based crop models use RUE-based functions to moderate the effect of [eCO 2 ] on biomass accumulation. Yin and Struik 20 proposed a new framework to quantify the conversion efficiency of incident solar radiation into phytoenergy by annual crops. They indicated that for C 3 crops the overall efficiency of converting incident solar radiation into phytoenergy was 2.2% (RUE = 1.22 g MJ −1 ) under 400 μmol mol −1 [CO 2 ] which could be increased to 3.6% (RUE = 1.75 g MJ −1 ).
Different approaches have been used to study the effect of [eCO 2 ] on plants growth, development and yield. These include FACE experiments 9,12,21 , open top chamber (OTC) [22][23][24] , temperature gradient tunnel (TGT) 25 and crop modeling. The FACE approach is considered more appropriate compared to other experimental approaches as it can provide data that better resemble field conditions. In general, C 3 cereal crop response to [eCO 2 ] under water stress is comparatively higher (22%) than under irrigated conditions (16%) 26 . Similarly, increased photosynthesis (10-45%) in C 3 crops with increased canopy temperature, yield, biomass and water use efficiency and decreased stomatal conductance and evapotranspiration have been reported under FACE experiments 9 . Wheat, the main C 3 cereal crop, showed reduced stomatal conductance and evapotranspiration with increased photosynthesis and canopy temperature under [eCO 2 ]. This resulted in higher biomass and yield in wheat even under water stress conditions. Hocking and Meyer 27 reported doubled drymatter production in wheat under [eCO 2 ] treatments compared to control. Meanwhile, higher water use efficiency (19-23%) under the high N treatment was reported in wheat under FACE 28 . FACE experiments from Australia and China reported a 21-23% increase in biomass and 24.8% increase in wheat grain yield under [eCO 2 ] 29 .
Crop simulation models often used to study crop behavior under changing climate. Many researchers have used crop modeling under different climatic scenarios 16,26,[30][31][32][33][34][35][36][37][38][39][40][41][42] . The Agricultural Model Intercomparison and Improvement Project (AgMIP) studied the impact of climate change on agricultural production and food security using process-based crop models 43 . Crop model comparison under the AgMIP framework revealed that uncertainties in wheat yield simulation increased with increased temperature-by-CO 2 interactions 44 . Similarly, Asseng et al. 39 tested 30 different wheat crop models in response to elevated temperature and predicted that most of the models simulated yield well under baseline conditions, but with increasing spread at the higher future temperatures. Furthermore, Rosenzweig et al. 45 found strong negative effects of climate change particularly at higher temperature. However, they recommended further research to minimize uncertainties related to the representation of carbon dioxide, nitrogen, and high temperature. Most of the earlier crop modeling work was focused more on studying and quantifying the impact of temperature on crop growth, development and yield [46][47][48][49][50] . They, in general, found a reduction in grain yield with some level of uncertainty under higher temperature. Similarly, earlier modeling studies focused on the combined effect of climatic parameters i.e. [eCO 2 ], temperature, nitrogen and drou ght 9,24,34,51,52 . Some of the earlier work studied the interaction of increased temperature and CO 2 [53][54][55] .
The interaction between increasing temperature and [eCO 2 ] is difficult to isolate, and the interpretation of projections tend to focus on crop model performance under warming, with limited attempts to understand model performance solely in response to [eCO 2 ]. For example 56 , evaluated the integrated effect of temperature and CO 2 on wheat phenology and yield using CERES and N-Wheat. Similarly, interactive effect of CO 2 and temperature on soybean [Glycine max (L.) Merr.].
water use efficiency (WUE), foliage temperature, canopy resistance and evapotranspiration were studied earlier 57 . Root Zone Water Quality Model (RZWQM2) was used to model current and future climate change effects on winter wheat production but again they studied CO 2 fertilization and warming effects in combination 58 . However, a systematic comparison of [eCO 2 ] responses, independent of temperature, of crop models often used for climate change projections has not been attempted to our knowledge. The present study was designed to evaluate the performance of five process-based crop models (APSIM-Wheat, CropSyst, DSSAT-CERES-Wheat, EPIC and STICS) under different levels of [eCO 2 ]. The objectives of the present study were to (i) quantify impact of [eCO 2 ] on winter wheat biomass and yield and (ii) bring/suggest accuracy in the models' response to [eCO 2 ] and determine whether a multi-model ensemble approach would minimize uncertainty in climate change simulation.

Results
Biomass response to eCo 2 . The simulated biomass results at the high rainfall site (Pullman) depicted bias among models in response to [eCO 2 ]. All models provided different standard errors (Table 1). Significant difference among models for simulated biomass at ambient carbon dioxide concentration was observed with highest value simulated by CropSyst and lowest by DSSAT. However, with carbon dioxide concentration at 700 µmol mol −1 the simualted biomass started increasing with the largest response in STICS, followed by CropSyst and APSIM. A smaller response for biomass was simulated in EPIC followed by DSSAT at 700 µmol mol −1 . Nevertheless, the overall response to [eCO 2 ] was similar among models (Table 1). [eCO 2 ] increased the simulated biomass on average by 34% among all models. The percentage increase of simulated biomass from ambient CO 2 concentration (400 µmol mol −1 ) to 1000 ppm was 39, 37, 34, 33 and 25% for STICS, DSSAT, CropSyst, APSIM and EPIC, respectively.
Biomass accumulation at the low rainfall site (Lind) was more responsive to [eCO 2 ] compared to the high rainfall site among all models (Table 1). Overall, biomass accumulation among models was similar but a higher biomass was simulated in the EPIC crop model. Table 1 showed that at [aCO 2 ], STICS responded lowest for biomass accumulation followed by DSSAT while on increasing [CO 2 ] to 1000 ppm, biomass accumulation increased significantly. The standard error for biomass accumulation in response to [eCO 2 ] remained highest in DSSAT followed by EPIC. The range of standard error for DSSAT was 541 to 884 while for EPIC it was 519 to 727. On average, all models simulated 44% increase in biomass. The highest biomass percentage change from [aCO 2 ] (400 µmol mol −1 ) to to [eCO 2 ] (1000 µmol mol −1 ) was shown by EPIC (48%) followed by STICS (46%). However, overall order of increase was EPIC > STICS > CropSyst > APSIM > DSSAT.
The effect of [eCO 2 ] on biomass accumulation at the irrigated site (Moses Lake) revealed that it did not increase significantly with increased [CO 2 ] compared to the other two sites. The highest response to [eCO 2 ] was observed for CropSyst which ranged from 23198 to 27860 kg ha −1 . APSIM simulation for biomass accumulation in response to [eCO 2 ] ranged from 19301-22947 kg ha −1 with standard error of 312-350. Similarly, DSSAT simulated wheat crop biomass was 18216 kg ha −1 at [aCO 2 ] (400 µmol mol −1 ) which increased to 26615 kg ha −1 at 1000 ppm, [CO 2 ]. The standard error for simulated biomass was highest for DSSAT (Table 1). Genearlly, all five crop models agreed with the assumptions that elevated atmsospheric CO 2 concentrations increases crop biomass but the effect was more prominent at the low rainfall site (Lind) compared to the well watered site (Moses Lake). The average percentage change in response to [eCO 2 ] by all models at the irrigated site was 22%. The highest percentage change in biomass from ambient to [eCO 2 ] concentration was observed in DSSAT (32%) followed by STICS (27%). However, the lowest change was simulated in APSIM. The overall order of increase in biomass was DSSA T > STICS > EPIC > CropSyst > APSIM.
The outcome of biomass ratio by APSIM against increased CO 2 concentration showed an increasing trend with the highest response at the Lind compared to the irrigated one (Moses Lake). The biomass ratio remained similar between low and high rainfall sites at 500 µmol mol −1 but at higher concentrations it became significantly different (Fig. 1a). However, the trend was higher at the low rainfall site and lower at the high rainfall site after 900 µmol mol −1 . The biomass ratio increase at the irrigated site remained close to 1.2 while at the low and high rainfall sites it reaches to 1.5 and 1.8, respectively (Fig. 1a).
The response of CropSyst simulation to elevated CO 2 for biomass ratio revealed that it remained highest at the dryland site compared to the irrigated and high rainfall sites. The overall trend among all locations was increasing but at the irrigated site the increase was not too high (Fig. 1b). Overall, CropSyst and APSIM depicted similar trends for biomass ratio at all three locations.
With the increased CO 2 concentrations from baseline 400 µmol mol −1 to 1000 µmol mol −1 , the DSSAT crop model depicted a linear increase in biomass ratio among all locations (Fig. 1c). However, like APSIM and CropSyst, the biomass ratio remained highest at dryland site followed by the high rainfall and irrigated sites. The difference among sites in response to [eCO 2 ] was not too much as seen earlier in APSIM and CropSyst (Figs 1a aCO 2 (µmol mol −1 ) eCO 2 (µmol mol −1 ) www.nature.com/scientificreports www.nature.com/scientificreports/ and 2c). Similarly, in contrast to APSIM and CropSyst, the DSSAT response to [eCO 2 ] increased linearly. To some extent, APSIM and CropSyst depicted sigmoid trends for simulated biomass ratio in response to [eCO 2 ].
EPIC depicted an increasing trend for biomass ratio in response to [eCO 2 ]. The trend was more prominent at the low rainfall site. The biomass ratio was similar at high rainfall and irrigated sites from 400 to 500 µmol mol −1 while at higher concentrations it increased at the high rainfall site compared to the irrigated site (Fig. 1d). for APSIM, CropSyst, DSSAT, EPIC, and STICS at three Pacific Northwest sites. Three coloured lines represents sites at Pacific Northwest with ± standard errors and each box represents Model outputs for Biomass and Yield. Means are averaged over three replicates.
However, the EPIC biomass ratio response for high rainfall and irrigated sites remained significantly different from the other models.
The STICS response to [eCO 2 ] showed a similar increasing trend like other models for the dryland site followed by the high rainfall site (Fig. 1e). However, the response remained similar among high and low rainfall site up to 500 µmol mol −1 while at higher concentrations it increased sharply at the low rainfall site. Biomass ratio at the irrigated site increased linearly but remained lower than the other two sites.
Yield response to eCo 2 . The impact of [eCO 2 ] on grain yield showed increasing trends by all five crop models at the high rainfall site (Pullman). The highest yield at [aCO 2 ] 400 µmol mol −1 was simulated by STICS (7101 kg ha −1 ) followed by APSIM (6368 kg ha −1 ). The lowest yield at [aCO 2 ] was predicted by DSSAT (6214 kg ha −1 ). The range of standard errors at 400 µmol mol −1 was 209-239 with a highest standard error with EPIC. However, on moving from ambient CO 2 concentration (400 µmol mol −1 ) to 700 µmol mol −1 , the highest yield was simulated by STICS (10013 kg ha −1 ) followed by APSIM and CropSyst. The standard error range for grain yield at 700 µmol mol −1 was 255-323 with the highest standard error shown by STICS. The yield response at 1000 µmol mol −1 CO 2 concentration showed that STICS simulated the highest yield (11353 kg ha −1 ) while the lowest was depicted by DSSAT (8835 kg ha −1 ). The standard error at 1000 µmol mol −1 CO 2 ranged from 274-343 (Table 2). Overall, the average percentage increase in grain yield by all models from 400 to 1000 was 34%. However, the highest increase was observed in CropSyst and STICS (38%) followed by APSIM where it remained 37%. The lowest percentage change was observed for EPIC (26%) followed by DSSAT.
The simulated grain yield increased significantly as a function of CO 2 concentration at the low rainfall site ( Table 2). The average yield increase among all models from [eCO 2 ] to 1000 µmol mol −1 was 45%, comparatively higher than high rainfall site where it was 34%. The highest percentage increase was shown by EPIC (49%) and STICS (46%) followed by CropSyst (44%), APSIM (43%) and DSSAT (42%). The grain yield at [aCO 2 ] was in the range of 2339-3159 kg ha −1 with the highest yield simulated by EPIC. Similarly, standard error ranged from 134-249 with the highest error depicted by EPIC. On increasing CO 2 concentration, yield increased linearly with the highest response shown by EPIC. The yield response trend at 700 µmol mol −1 CO 2 revealed that EPIC (4937 kg ha −1 ) predicted highest wheat yield followed by CropSyst (4189 kg ha −1 ), STICS (4017 kg ha −1 ), DSSAT (3675 kg ha −1 ) and APSIM (3503 kg ha −1 ). The standard error at 700 µmol mol −1 CO 2 concentration was in the range of 158-304 with highest standard error depicted by DSSAT. Similarly, at 1000 µmol mol −1 CO 2 the yield remained highest (6162 kg ha −1 ) for EPIC followed by CropSyst (5058 kg ha −1 ). The standard error ranged from 195-376 with the highest error depicted by DSSAT.
The models' simulated results of grain yield for the irrigated site (Moses Lake) in response to elevated [eCO 2 ] showed a linear relationship ( Table 2). The highest response was with CropSyst at ambient as well as eCO 2 compared to other models. The highest yield at [aCO 2 ] was 10559 kg ha −1 with a standard error of 184. However, the highest standard error (350) was observed for DSSAT with the grain yield of 8322 kg ha −1 . Overall, the average grain yield at [aCO 2 ] combind over all models was 8608 kg ha −1 . The range of standard error at 400 µmol mol −1 CO 2 was 142-350. On increasing CO 2 from 400 to 700 µmol mol −1 , the yield increased among all models. The maximum yield was simulated by CropSyst (12213 kg ha −1 ) followed by DSSAT at 700 µmol mol −1 CO 2 . The range of standard error at 700 µmol mol −1 CO 2 was 165-367 with highest standard error for DSSAT. A similar trend was Asymptotic DSSAT response to elevated CO 2 . Wheat uses an asymptotic look-up multiplier on RUE for the relative response to elevated CO 2 to produce biomass. The asymptotic look-up multiplier for modeled effects of elevated CO 2 on RUE is given in the WHCER045.spe file. Regression equation was fitted to show the response of CO 2 factor, photosynthesis with CO 2 reference and validated by coefficient of determination (R 2 ). www.nature.com/scientificreports www.nature.com/scientificreports/ observed at 1000 µmol mol −1 CO 2 with standard error in the range of 169-395. The average percentage increase in grain yield among all models from 400 to 1000 µmol mol −1 CO 2 was 22% which was almost 50% less than the Lind. The highest percentage increase was with DSSAT (30%) followed by STICS (26%  2 ] at all sites and it was 12-19% for Pullman, 11-20% for Lind and 14-18% for Moses Lake respectively ( Table 3).
The yield ratio for APSIM at three sites showed that the value of this ratio was a function of CO 2 (Fig. 1f). However, this response was significantly different among the several climatic locations. The highest yield ratio was obtained for water stress conditions where it increased from 1 to 1.8 followed by the high rainfall site (Pullman) where the change was from 1 to 1.6. However, the lowest yield ratio was obtained at the irrigated site where it remained in the range of 1-1.2.
The yield ratio response to [eCO 2 ] for CropSyst depicted higher increasing trend under water stress conditions compared to irrigated (Fig. 1g). The yield ratio response between high and low rainfall was similar at 500 µmol mol −1 CO 2 but at higher concentrations, the increase was more prominent at the low rainfall site. The yield increase ratio at the water stress site was from 1 to 1.8 while under high rainfall conditions it was from 1 to 1.5. The response to [CO 2 ] at the irrigated site remained in the range of 1 to 1.2. The yield ratio response of APSIM was similar to CropSyst (Fig. 2d-f).
The simulated yield ratio response of DSSAT revealed that it was more linear and directly related to [CO 2 ] under water stress conditions than APSIM and CropSyst (Fig. 1h). However, the increased ratio (1-1.7) was a little bit lower than DSSAT. The other difference in DSSAT simulated yield ratio response compared to APSIM and CropSyst was that it showed similar trends for high rainfall and irrigated conditions. The range was from 1-1.4.
The yield ratio trend simulated by EPIC was significantly different among locations (Fig. 1i). The highest yield ratio (1-2.0) was observed under water stress conditions. EPIC response to [eCO 2 ] was similar to other crop models (APSIM, CropSyst and DSSAT) particularly under water stress conditions. However, under high rainfall and irrigated conditions, the yield ratio remained similar untill 500 µmol mol −1 above which it increased at the high rainfall site. The range of increase at the high rainfall site was 1-1.4. In contrast to other crop models, EPIC depicted a decreasing trend under irrigated conditions. aCO 2 (µmol mol −1 ) eCO 2 (µmol mol −1 )  www.nature.com/scientificreports www.nature.com/scientificreports/ With the increase in CO 2 from baseline 400 µmol mol −1 to 1000 µmol mol −1 , the STICS model depicted yield increase as a function of CO 2 concentration at all locations (Fig. 1j). The dominant trend was observed at the dryland site similar to previous model results followed by yield ratio outcomes at the high rainfall site. The lowest ratio was obtained under non-stressed conditions (Fig. 1j). The range of yield ratio increase for the dryland site was 1-1.8 while for high rainfall site it was 1-1.6. The increased yield ratio at the irrigated site was in the range of 1-1.3.

Discussion
Crop models can be used to see the impact of various climatic variables on crop biomass and yield under changing climate. These models have been continuously refined to give a realistic picture of different climatic variables impacts on crop production. However, models have limitations in the prediction of crop responses to global change factors 59 . Challinor et al. 60 and Soussana et al. 33 emphasized the need to improve crop models for assessments of climate change. Since most of the models need improvements in climate impact assessment, we evaluated five different process-based models under eCO 2 . The models we tested simulated increased biomass on average by 34% but this trend was greatest under dryland conditions as compared to irrigated. This difference could be because of lower stomatal conductance and reduced transpiration which might result in higher water use efficiency. Earlier work concluded that [eCO 2 ] alone will increase biomass and yield in C3 crops as photosynthesis of C3 plants are not CO 2 -saturated and photosynthesis rates increases under [eCO 2 ]. The two major plant responses e[CO 2 ] are (i) increased net photosynthesis with a consequent increase in crop growth and yield, and (ii) decreased stomatal conductance and increase crop water use efficiency 12,61 . Similar to our work some studies suggest that these responses become very important under water-limited conditions and they reported greater e[CO 2 ] response by plants under drier conditions because of greater water use efficiency 61 . Furthermore, buffering action of e[CO 2 ] has been reported against heat waves which resulted in the increased crop production under semi-arid environments. Specific genotypic adaptation strategies have been suggested to capture the positive effects of elevated [CO 2 ] under drier conditions 62 . Thus, under dryland conditions, small amounts of water could contribute to enhanced photosynthate production and its translocation to the grain 63 . Manderscheid et al. 64 reported increased biomass under [eCO 2 ] by 17% with significant CO 2 × N interaction and this can be increased further up to 30% by +200 ppm of CO 2 . Hence, similar to our modeling outcomes it has been suggested that under eCO 2 crop water use of well-fertilized wheat will improve due to reduction in seasonal evapotranspiration. Fitzgerald et al. 65 reported a large response to e[CO 2 ] and suggested field level research to provide a detailed mechanistic understanding for adapting crops to climate change. O'Leary et al. 29 reported that elevated CO 2 (700 μmol mol −1 ) increased leaf area index (21%), photosynthetic area index (25%) and biomass (23%). However, our biomass response to eCO 2 was slightly higher than reported by Ainsworth and Long 66 and was consistent with Kimball 9 . e[CO 2 ] experiments in Nanjing China who reported a 13.6% increase in biomass, but their maximum CO 2 was 550 μmol mol −1 . The relatively small increase in their experiment was because of a decrease in biomass between heading and maturity while positive effects were found during pre-anthesis growth phases which resulted in higher grain yield 67   www.nature.com/scientificreports www.nature.com/scientificreports/ which results in higher photosynthesis and biomass production. According to Long et al. 68 a 200 ppm rise in CO 2 could theoretically result in a 40% increase in photosynthesis per unit of leaf area which is considerably higher than experimental results 69 . Fitzgerald et al. 70 tested the hypothesis that biomass and yield response to e[CO 2 ] is greater in semi-arid agroecosystems or not and they measured biomass and yield increase up to 79% well above the measured highest response (34%) by Liu et al. 71 . Shimono et al. 72 suggested different prescreening techniques which could be used for identification of elevated CO 2 -responsive genotypes. However, among these techniques, we suggest process-based crop models as a good option to test CO 2 -responsive genotypes among different plant species. Kumagai et al. 73 proposed that the Finlay-Wilkinson relationship could be used as a pre-screening criterion for e[CO 2 ] responsiveness, which is regression based, while crop models provide a detailed mechanistic approach.
Variable yield responses were simulated by models under e[CO 2 ], however, the largest response was observed at Lind compared to Moses Lake, where it was lowest ( Table 2). This magnitude of yield response was higher than the findings of O'Leary et al. 29 and Yang et al. 74 who reported an average increase of 26% and 24.8% respectively. Deryng et al. 75 reported increased crop yield under rising CO 2 which might be because of enhanced photosynthesis and reduced water use. Simulation modeling is a complex process as it involves interactions of multiple factors thus simulation by individual crop model might contain uncertainties. Quantifying the impact of climate change (e.g. elevated CO 2 , rising temperature and variability in rainfall) on crop yield by single crop model is problematic as suggested by different researchers 44,76,77 . Therefore, ensemble multi-models for simulation of climate change impacts could be a good option as proposed in the present study and in earlier work. Transpiration efficiency (TE) needs to be considered as an important mechanism in models for explaining response of increased crop yield under dryland conditions and elevated CO 2 78 . Similar to our findings yield stimulation due to e[CO 2 ] was reported by Fitzgerald et al. 70 but it was in the higher range of 24% to 70% depending upon environment (Fig. 1). Similarly, they reported that heat wave effects could be ameliorated under e[CO 2 ]. Long et al. 68 reported that [CO 2 ] of 200 ppm above ambient could result in a 10-20% increase in crop yield. However, Ziska et al. 79 studied the impact of [CO 2 ] on quantitative and qualitative traits of wheat varieties and concluded that elevated CO 2 resulted in increased seed yield while grain and flour protein declined. Similarly, lower nutritional quality in grains of non-legume crops was reported by Jin et al. 80 due to e[CO 2 ]. C-N-P stoichiometry of terrestrial plants to the rising CO 2 concentration showed that concentrations of N, P and N:P will decrease by −9.73%, −3.23% and −7.23% while C and C:N will increase by +2.19 and +13.29% respectively 81 .
Water use efficiency is the amount of grain produced per unit of water used by the crop. Elevated CO 2 concentration improved the water use efficiency with more positive effect under stress condition compared to the well-watered plants. This could be due to the fact that under water-limiting and [eCO 2 ] conditions photosynthetic CO 2 uptake response increases resulting in the higher CO 2 fixation. In present studies, the highest WUE was simulated by the crop model EPIC at low rainfall site while for all other models increasing trend was observed. Different studies 10,19,29,82 discussed the positive effect of [eCO 2 ] on WUE which could be due to reduced stomatal conductance resulting to lower canopy transpiration and crop water use 12,83,84 . However, negative WUE were reported by earlier researchers and they suggested upgradation in the code of models to have a sufficiently strong effect of CO 2 on stomatal conductance and on transpiration 85 . Similarly, stomatal resistance also regulates photosynthesis and transpiration and it affected by CO 2 and vapor pressure deficit (VPD) 86 . Thus, leaf level CO 2 exchange rate and stomatal closure have association with VPD which increases with temperature and have a strong relationship with radiation use efficiency (RUE) [87][88][89][90] . Our study suggests integrating the effects of all these crucial factors in the models so that they can simulate results in a biologically realistic manner. Since most of the models are unable to accommodate these factors particularly under [eCO 2 ] thus we recommend using model ensemble or adapt physiological mechanisms in the model. Therefore, it is necessary to consider these findings in process-based models to have better response under eCO 2 . The response to e[CO 2 ] in our studies suggested further evaluation/improvement of models, particularly under stressed conditions. Models could be improved by local calibration with consideration of radiation use efficiency (RUE) and transpiration efficiency(TE) methods of biomass accumulation. Similarly, temperature and light intensity interactions with eCO 2 should be considered, which will render models more effective for future climate change studies.

Conclusion
Model response to [eCO 2 ] showed significant increases in biomass and yield of wheat. Overall models were able to capture [eCO 2 ] response but with differences in response to environmental conditions. The response was higher under dryland conditions compared to irrigated which could be because of lower stomatal conductance and transpiration resulting in higher water use efficiency. However, to have more accurate simulation results from models it is important to calibrate the model under local dryland conditions and consider the interactive effect of light intensity with [eCO 2 ]. In the future models could be used to pre-screen large numbers of germplasms for [eCO 2 ] responsiveness at relatively low cost. Process-based crop models have variability in the simulation of crop response to elevated CO 2 with a greater difference under water-stressed conditions. An ensemble approach will increase the accuracy of model response to elevated CO 2 .

Materials and Methods
Five process-based crop models were evaluated in the present study: APSIM-Wheat 91,92 , CropSyst 93,94 , DSSAT-CERES-Wheat 95 , EPIC 96 and STICS 97,98 . These models were chosen because of their wide use in climate change studies and their ready availability. The focus of this evaluation was on biomass and yield responses to [eCO 2 ], with attention to changes in crop transpiration, biomass and grain yield. The relevant details of the approaches used in these models to simulate the response to [eCO 2 ] are described below.
www.nature.com/scientificreports www.nature.com/scientificreports/ APSIM. APSIM (Agricultural Production system SIMulator) is a dynamic daily time step model that can simulate crop growth, development and yield using different management and biophysical modules. It is capable of simulating soil C, N, P and water dynamics in interactions with different management/crops systems driven by daily meteorological data (Solar radiation, maximum and minimum temperatures, rainfall). Radiation-use efficiency (RUE) approach is used to calculate daily potential production of crops which is then limited to actual above ground biomass production on a daily time step basis by N, P and soil water availability 92 . Soil water balance is simulated by cascading bucket approach of CERES 99 through SOILWAT and SWIM3 modules 100 2 ] from 350 to 700 ppm at a mean temperature of 20 °C will bring a 21% increase in RUE 101 . RUE is scaled by the ratio of light limited photosynthetic response at the elevated CO 2 compared with CO 2 at 350 ppm. CropSyst. CropSyst is multi-year, multi-crop, daily time-step cropping system model. It can simulate the effect of management, soil and climate on crop growth, development and yield. Detail of CropSyst model is available on the website (http://modeling.bsyse.wsu.edu/CS_Suite_4/CropSyst/index.html). The [eCO 2 ] effect on crop biomass can be simulated by CropSyst which relies on biomass accumulation under experimental elevated CO 2 (C a ). C a in CropSyst is expressed as the ratio of growth under elevated C a (C eo ) to growth under control baseline C a (C bo ). The subscript, "o", represents experimental conditions. A Michaelis-Menten-type equation models relative biomass growth (G r ) in response to C a (Eq. 2).
r a x a G r (relative biomass growth similar to leaf photosynthesis response to intercellular CO 2 concentration) will be less than 1.0 if C a < C bo and vice versa. G x (maximum growth increases relative to baseline conditions) and s can be obtained after considering G r = 1.0 (C a = C bo ) and equal to G ro when C a = C eo (Eqs 3-4).
The value of e for a given elevated C a (e e ) is calculated by considering e e = eGr. Stomatal resistance (u) is also considered in CropSyst under elevated C a as it reduced transpiration ratio F (elevated to baseline Ca crop transpiration per unit leaf area). Thus, u under elevated Ca (ue) is given by ue = uGr/F. Both e and u are parameters specified for baseline Ca conditions. Stockle et al. 86 , presented calculation of stomatal resistance as a function of C a . Meanwhile, F could be calculated as the ratio of elevated to baseline Ca considering aerodynamic and canopy resistance 94 .

CERES-Wheat. CERES-Wheat embedded in Decision Support Systems for Agrotechnology Transfer
(DSSAT, v 4.7) uses daily time step from planting to maturity to simulate the growth and development of crops 105 . Potential growth (G p ) is a function of photosynthetically active solar radiation (PASr) and its interception (Si) by crops but limited by suboptimal temperature, soil water, N and P deficits. Cardinal temperature approach has been used in CERES-Wheat to simulate temperature effects on crop growth and grain filling with an optimum temperature of 34 °C 106,107 . CERES-Wheat uses an asymptotic look-up multiplier on RUE for the relative response to elevated CO 2 to produce biomass. The asymptotic look-up multiplier for modeled effects of elevated CO 2 on RUE is given in the WHCER045.spe file (Fig. 2). The CERES-Wheat model simulates the effect of elevated CO 2 on actual transpiration by increasing stomatal resistance as a function of CO 2 concentration.
An approach for reducing transpiration as a function of rising CO 2 was developed for the DSSAT models in the early 1990s by J.W. Jones and L.H. Allen (personal communication, see TRANS routine of DSSAT code). The computations include equations for leaf stomatal resistance (R s ) response to 330 ppm or current CO 2 , whole canopy stomatal resistance (R c ) to reference or current CO 2 (dividing R s by total LAI), and canopy boundary layer resistance (R a ) as a function of LAI. Finally, a ratio effect (T ratio ) of CO 2 (current CO 2 versus 330 ppm reference CO 2 ) to reduce daily transpiration is computed in the following equation, considering the psychrometric constant (δ), gamma (γ), canopy resistances (R c ), and boundary resistance (R a ) (Eq. 5). Soil data. The soil of the Pullman site was silty clay loam with bulk density of 1.35 g cm −3 . Soil texture at Lind was coarse silt loam with bulk density of 1.31 g cm −3 . Field capacity in the specific root zone at Pullman was 0.30 mm mm −1 while wilting point water in specific root zone was 0.12 mm mm −1 . The soil series of Pullman was Palouse which is a deep well drained soil. The textue of Lind was coarse silt loam having sand, silt and clay percentages of 21.7, 70.8 and 7.5, respectively. The bulk density at Lind was 1.31 g cm −3 with field capcity of 0.33 mm mm −1 and wilting point of 0.007 mm mm −1 .The soil series at Lind was Lind which is a deep poorly drained soil. The irrigated site, Moses Lake, was in the Ephrata soil series, and had sandy loam soil texture with drain upper limit of 0.40 mm mm −1 while the lower limit was 0.27 mm mm −1 . Bulk density of soil at Moses Lake was 1.41 g cm −3 . www.nature.com/scientificreports www.nature.com/scientificreports/ biomass and yield at 400 ppm were considered as baseline. The ratio in biomass and yield change were calculated for all other concentrations of CO 2 using 400 ppm as baseline. The ratio was plotted against CO 2 concentration to see models' response to [eCO 2 ] and uncertainty among the models. Average with standard error for WUE was calculated in response to [eCO 2 ].   Table 4. Input data of study sites for model calibration.