The experiment was conducted in two environmentally controlled greenhouses with the same characteristics (15 m in length, 10 m in width, and 3.5 m in height, north-south oriented) under spring-summer climatic conditions from May to August 2018. Two widely grown tomato cultivars (JinPeng NO.1, CV1, JinPeng & Co., Ltd., China; FenGuan, CV2, ZhongYa & Co., Ltd., China) with relatively distinct VPD responses were examined (Du et al., 2020). Seeds were sown in plugs for germination and transplanted at the four-leaf stage to 4.5 L plastic pots containing the same amount of organic substrate and perlite mixture in a 3:1 proportion (v/v). Soil moisture was maintained at ~90% container capacity according to a previous method (Zhang et al., 2017). Plants were periodically trimmed to maintain rapid vegetative growth throughout experiments. Plants were grown in two environmentally controlled greenhouses and maintained under the same growth conditions but contrasting VPD. A high VPD was achieved in a natural greenhouse environment, with a VPD of ~3–5 kPa around midday, while low VPD was maintained in the range of.5–1.5 kPa by humidification. A high-pressure micro-fog system was activated when the VPD exceeded the target values, and the characteristics of the system were described in detail in a previous study by Zhang et al. (2018). The average daily meteorological data inside the greenhouse during the growth period were ~maintained at a temperature of 20–32°C, relative humidity of 50–75%, and photosynthetically active radiation of 45–65 Wm⁻².

The effects of VPD perturbations on leaf photosynthetic performance and plant water status were investigated ~50 days after treatments. Afterward, 15 uniform plants from each treatment were selected as samples and transferred to growth cabinets in the evening prior to photosynthetic measurements. The light and temperature of the growth cabinets were controlled steadily at normal levels throughout the experiment.

Leaf gas exchange and chlorophyll fluorescence were measured simultaneously on healthy and expanded leaflets at the same nodes by portable gas exchange systems equipped with a leaf chamber fluorometer (LI-6400, Li-Cor, Inc., Lincoln, NE, USA). All portable gas exchange systems were enclosed in growth cabinets. The VPD inside cabinets and the leaf chamber was simultaneously controlled across a series gradient of 0.5, 1.5, 2.5, 3.5, and 4.5 kPa. The temperature, light, and CO₂ concentrations were controlled at the following constant and steady conditions throughout the experiment: temperature of 28 ± 1°C; saturating photosynthetic photon flux density (PPFD) of 1,100 μmol m⁻² s⁻¹; CO₂ concentration of 400 μmol mol⁻¹. The VPD was increased stepwise across the gradients for at least 60 min until photosynthesis and the plant water status achieved a new steady state.

The curve of the photosynthetic rate (P_n) vs. intercellular CO₂ concentration (C_i) was determined using a previous procedure (Li et al., 2019b), across the VPD range of 0.5–4.5 kPa. Briefly, a P_n-C_i curve was generated by controlling the ambient CO₂ concentration (C_a) from 400 to 300, 200, 150, 100, and 50 μmol mol⁻¹ and then increased to 400 μmol mol⁻¹. After re-achieving a steady-state at 400 μmol mol⁻¹, C_a was increased gradually from 400 μmol mol⁻¹ to 1,200 μmol mol⁻¹. The carboxylation efficiency (CE) was estimated according to linear regression of the P_n-C_i curve in the range of C_a ≤ 200 μmol mol⁻¹ (Sun et al., 2016). The maximum rate of Rubisco carboxylation capacity (V_cmax) and maximal rate of electron transport (J_max) were determined according to the FvCB model (Farquhar et al., 1980).

Carbon dioxide diffuses via stomatal and mesophyll barriers in a series circuit, which was driven by the CO₂ partial pressure gradient (Li et al., 2019b). Stomatal conductance for CO₂ diffusion (g_sc) was determined according to the water diffusion conductance (g_sw) and the ratio between molecular diffusivities of water and CO₂ in gas (Giuliani et al., 2013). The mesophyll conductance (g_m) was estimated by the variable J method (Harley et al., 1992):

where P_n is the net photosynthetic rate and C_i is the intercellular CO₂ concentration. P_n and C_i were measured by steady-state gas exchange. R_d is the mitochondrial respiration rate in the light, and Γ^* is the CO₂ compensation point inside the chloroplast. R_d and Γ^* were calculated according to a previous study by Laisk and Oja (1998). Briefly, P_n-C_i curves were measured at two light intensities (75 and 500 μmol m⁻² s⁻¹) at CO₂ concentrations of 30–120 μmol mol⁻¹. Γ^* (x-axis) and R_d (y-axis) were derived according to the intersection point of the P_n-C_i curves. J is the electron transport rate, which was calculated as described by a previous study (Tomas et al., 2013).

According to the series circuit, the total CO₂ diffusion resistance (1/g_tot) can be determined as 1/g_tot = 1/g_s + 1/g_m (Niinemets et al., 2009). Therefore, g_tot can be determined as:

The photosynthetic limitation was divided into the components of stomatal limitation (L_s), mesophyll limitation (L_m), and biochemical limitation (L_m). The proportions of individual components imposed on photosynthesis were determined as follows (Muir et al., 2014; Li et al., 2019b):

By definition, L_s + L_m + L_b =1; ∂A/∂C_c was determined as the slope of the P_n-C_c curves at CO₂ concentrations of 40–110 μmol mol⁻¹.

Once photosynthetic measurements were completed, the adjacent leaflets were harvested for the determination of the water status of the plant. The leaf water potential (Ψ_leaf) was measured by a pressure chamber (PMS-1000, PMS Instruments Inc., Corvallis, OR, USA). An extra test was performed where Ψ_leaf of adjacent leaflets was compared, and no differences in Ψ_leaf were detected between two adjacent leaflets. Some plants were kept in dark conditions for ~8–10 h for the determination of the soil water potential (Ψ_soil) (Tsuda and Tyree, 2000). Since water movement was ~zero under dark conditions, Ψ_soil remained relatively constant and can be assumed to equal the xylem pressure potential of leaves under dark conditions.

After determination of the plant water status, the leaflet area was measured by a leaf area metre. The leaflet samples were dried at 80°C in an oven to a constant dry mass and weighed. The leaf mass area (LMA) was determined as the ratio of leaf dry mass to leaf area.

After reaching the steady-state of photosynthesis, leaflets were harvested for ABA detection and transcriptome sequencing. Phytohormone contents were determined by a liquid chromatography electrospray ionisation tandem mass spectrometry (LC-ESI-MS/MS) system (HPLC, Shim-pack UFLC SHIMADZU CBM30A system, www.shimadzu.com.cn/; MS, Applied Biosystems 6500 Triple Quadrupole, www.appliedbiosystems.com.cn/). Briefly, the leaflets for photosynthetic measurements were harvested and frozen in liquid nitrogen. The samples were extracted with methanol/water/formic acid and filtered before LC-MS/MS analysis. The detailed protocol was described on MetWare (http://www.metware.cn/) based on the AB Sciex QTRAP 6500 LC-MS/MS platform. Samples were detected with three biological replicates.

Total RNA was extracted from leaflet samples for transcriptome sequencing, according to a previous study (Zhang et al., 2019b). Sequencing libraries were constructed using the Ultra^TM RNA Library Prep Kit for Illumina (NEB, USA) according to the instructions of the manufacturer. The detailed protocol was described in a previous study, which was briefly described in a simplified diagram.

The clean data were obtained by processing the raw data through in-house Perl scripts. The low-quality data and sequencing adapters were trimmed. The fragments per kilobase of transcript per million fragments mapped (FPKM) was calculated based on gene length and read counts. Differentially expressed genes (DEGs) were assigned according to the adjusted P < 0.05. Gene ontology (GO) enrichment was determined by submitting DEGs to the GO database to classify the genes. Kyoto Encyclopaedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg) was used to perform pathway enrichment analysis. Terms with corrected P < 0.05 were identified as significantly enriched by DEGs. To confirm the reliability of transcriptome sequencing, 10 candidates expressed genes in RNA-seq were simultaneously evaluated by qRT–PC analysis. The qRT–PCR values were linearly correlated with the RNA-seq FPKM values (P < 0.001).

All statistical analyses were performed using SPSS 19. One-way ANOVA was used to determine the significant difference of average values according to Tukey's test (P < 0.05). Regression analysis was performed by Microsoft Excel.

Vapour pressure deficit significantly affected the distribution of water potential along the soil-plant-atmospheric pathway (Figure 1). Atmospheric evaporative demand increased with VPD elevation, which triggered plant water stress and a linear decline in the leaf water potential (Figure 1A). With VPD elevation, the drawdown of Ψ_leaf in high-VPD-grown plants was less than that in low-VPD-treated plants according to the slope of linear regression (Figure 1A). The driving force for passive water flow between the soil and leaf (ΔΨ _soil−leaf) increased with the VPD, and the magnitude of the increase was greater in low-VPD-grown plants than high-VPD-grown plants (Figure 1B). Since Ψ_leaf was negligible compared with the large negative air potential, the water driving force at the leaf-air boundary (ΔΨ_leaf−air) increased dramatically with an increase in the VPD (Figure 1C). The magnitude of the increase in ΔΨ_leaf−air was considerably greater than that of ΔΨ_soil−leaf, and the difference between ΔΨ_leaf−air and ΔΨ _soil−leaf was enlarged with the VPD: the ratio of ΔΨ_leaf−air to ΔΨ_soil−leaf increased logarithmically from ~50 at 0.5 kPa to 150 at 1.5 kPa and then maintained at a steady level (Figure 1D). The statistical analyses of the plant water status are shown in Supplementary Table 1.

The photosynthetic rate responded to CO₂ elevation in similar patterns regardless of cultivar and VPD growth conditions: the photosynthetic rate rose rapidly across low CO₂ concentrations and then reached a steady state (Supplementary Figure 1). The maximum steady-state photosynthetic rate declined as the VPD increased from 0.5 to 4.5 kPa (Supplementary Figure 1). The maximum carboxylation rate (V_cmax), maximum electron transport rate (J_max), and CE declined linearly with VPD elevation (Supplementary Figure 2). The drawdown of V_cmax, J_max, and CE with VPD elevation was moderated in high-VDP-grown plants compared with low-VPD-grown plants according to the slope of linear regression (Supplementary Figure 2). The statistical analyses of photosynthetic parameters across VPD ranges are shown in Supplementary Table 2.

The stomatal, mesophyll, and total conductance for CO₂ diffusion decreased linearly with VPD elevation, regardless of the cultivar and VPD growth conditions (Figure 2). The magnitudes of drawdown in the stomatal, mesophyll, and total conductance were lower in high-VPD-grown plants than in low-VPD-grown plants for two cultivars (Figure 2).

The CO₂ concentration along the “source-path-sink” was reduced to different extents with VPD elevation (Figures 3A,B). The drawdowns of C_i and C_c caused by VPD elevation were relatively lower in high-VPD-grown plants than in low-VPD-grown plants (Figures 3A,B). Consequently, the CO₂ transport efficiency of C_i/C_a, C_c/C_a and C_c/C_i decreased linearly with VPD elevation. The declining slopes of C_i/C_a, C_c/C_a, and C_c/C_i vs. VPD were lower in high-VPD-grown plants than in low-VPD-grown plants (Figures 3C–E). The statistical analyses of CO₂ concentrations along the “source-path-sink” across VPD ranges are shown in Supplementary Table 3.

The fractions of stomatal, mesophyll, and biochemical limitations imposed on photosynthesis varied with VPD elevation (Figure 4). Under low VPD conditions, the stomatal and mesophyll conductance for CO₂ diffusion were high and imposed relatively minor limitations on photosynthesis. The stomatal and mesophyll conductance accounted for a low proportion of photosynthetic limitation, while biochemical carboxylation for carbon fixation was the most significant limitation for photosynthetic processes under low VPD conditions (Figure 4). The fraction of stomatal limitation increased linearly with the VPD, from ~15% at 0.5 kPa to 35% at 4.5 kPa (Figure 4A). A similar pattern was observed in the mesophyll limitation: the fraction of mesophyll limitation also increased linearly with VPD elevation, from ~23% at 0.5 kPa to 33% at 4.5 kPa (Figure 4B). The increments in the fractions of stomatal and mesophyll limitations tended to be less marked in high-VPD-grown plants. In contrast, the fraction of total limitations attributed to the biochemical limitation of carbon fixation gradually decreased linearly with VPD elevation, from ~65% at.5 kPa to 35% at 4.5 kPa (Figure 4C). The statistical analyses of stomatal, mesophyll, and biochemical limitation fractions across VPD ranges are shown in Supplementary Table 4.

Biochemical limitation accounted for the greatest limitation on photosynthesis under low VPD conditions, regardless of the cultivar and growth conditions (Figure 5). The limitations that stomatal and mesophyll conductance imposed on photosynthesis gradually increased and predominated under high VPD stress (Figure 5). Diffusion limitations, i.e., the sum of the stomatal and mesophyll resistance, were the rate-limiting step for the photosynthetic process under high VPD conditions, which imposed the greatest limitation on photosynthesis in tomato plants (Figure 5).

The mesophyll conductance was significantly and positively correlated with the stomatal conductance (Supplementary Figure 3A). Meanwhile, the stomatal and mesophyll conductance for CO₂ diffusion were closely linked to the leaf water status, wherein significant and positive correlations were found in the leaf water potential vs. the stomatal and mesophyll conductance (Supplementary Figures 3B,C). Acclimation to VPD modified leaf structural traits, wherein LMA tended to be slightly greater in high-VPD-grown plants than in low-VPD-grown plants (Supplementary Figure 4A). A significant and negative correlation between g_m and LMA was observed (Supplementary Figure 4B).

With VPD elevation, the foliar ABA content increased exponentially (Figure 6A). The leaf ABA content was linearly and negatively correlated with the CO₂ diffusion conductance of g_s and g_m (Figure 6B). The slope of linear regression in g_s was more negative than that in g_m, indicating that g_s was more sensitive to ABA in response to VPD stress (Figure 6B).

Kyoto Encyclopaedia of Genes and Genomes analysis showed that physiological processes of “metabolic pathway” and “plant hormone signal transduction” were involved in the response to VPD and were potentially associated with ABA biosynthesis and signal transduction (Figure 7). “Metabolic pathway” was the dominant pathway in response to VPD elevation. “Plant hormone signal transduction” also exhibited a significant pathway in response to VPD stress in the range from 1.5 to 3.5 kPa (Figure 7).

The enriched genes in the comparison between different VPD treatments were annotated in three main GO categories: biological process, cellular component, and molecular function. A great difference in the top 50 GO enrichment terms was observed between different VPD conditions (Figure 8). In the comparison of 0 with 1.5 kPa and 0.5 kPa with 3.5 kPa, the “ABA-activated signalling pathway” was significantly enriched and was involved in the VPD response (Figures 9A,C). In the comparison of 0.5 kPa with. 2.5 kPa, ABA was potentially involved according to the enriched terms of “cellular hormone metabolic process,” “hormone biosynthetic process,” and “hormone metabolic process”. Gene expression with rising VPD can be classified into 10 patterns according to K-means analysis (Figure 9). Gene expression patterns were mostly classified into the pattern of “Subclass 6” with 768 genes, wherein gene expression remained relatively stable under mild VPD stress and increased dramatically under high VPD water stress (Figure 9). The genes associated with ABA biosynthesis and signal transduction followed different patterns.

Passive water movement was driven by the gradient of free energy along the soil-plant-atmospheric continuum, which could be quantified as the gradient in water potential in the liquid phase. Water movement at the leaf-air boundary in the gas phase was driven by the difference in the VPD. Based on physical principles, excessive air desiccation triggered a high VPD and great negative air-water potential. Δψ_leaf−air was substantially >Δψ_soil−leaf, which drove transpiration. The substantial difference between Δψ_leaf−air and Δψ_soil−leaf was logarithmically enlarged with an increase in the VPD (Figure 1). Quantitatively, the atmospheric driving force at the leaf-air boundary could be >100-fold larger than the soil-leaf component under high VPD conditions (Figure 1). The great asymmetry between the atmospheric evaporative demand and soil water supply triggered disruption in the water balance despite plants being well irrigated. Root water uptake and supply were inadequate to keep pace with the great atmospheric driving force under high VPD conditions, which consequently triggered leaf dehydration and decline in water potential. Therefore, the VPD is a crucial external stimulus moving water through a soil-plant-environment continuum. VPD fluctuates dramatically over the diurnal course in crop production, especially for greenhouse cultivation. Soil moisture is relatively stable over the short term compared with the VPD (Caldeira et al., 2014). Plant-water relations are regulated to a greater extent by the VPD and to a lesser extent by soil moisture. Similar to soil drought, VPD-induced plant water stress is also an important factor triggering photosynthetic depression.

The stoma is the “gatekeeper” for the exchange of water vapour and CO₂. Guard cells surrounding the stomatal pore respond to perturbations of the soil-plant-atmospheric hydraulic continuum, which is putatively transduced into stomatal movements by feedback and feedforward mechanisms (Buckley, 2005, 2017, 2019). Stomatal control of transpired water loss is critical for sustaining physiological processes, such as leaf water status and photosynthetic CO₂ uptake. It has been recognised that plants respond to drought by closing guard cells to prevent the development of water dehydration in plant tissues (Novick et al., 2016). In the present study, the atmospheric driving force was an order of magnitude greater than the water supply, which led to a great dissymmetry between the water supply and evaporative demand. The dissymmetry between the water supply and evaporative demand triggered declines in the leaf water potential and stomatal closure. However, the mechanism of VPD-triggered stomatal closure is still uncertain and is a “black box” (Buckley, 2016). Some hypotheses hold that stomatal closure in angiosperms under high VPD conditions is an active process that is regulated by hormonal and hydraulic signals (Merilo et al., 2018; Pantin and Blatt, 2018). The plant stress hormone ABA is continuously produced and delivered with a transpiration stream to guard cells (Qiu et al., 2017; Merilo et al., 2018). In the present study, leaf ABA rapidly accumulated with the rise in the VPD. Transcriptome analysis suggested that ABA biosynthesis and signal transduction were potentially involved in the response to the VPD. Based on the theory of ion channel-mediated guard cell signal transduction (Julian et al., 2001), hypothetical mechanisms of the ABA-mediated stomatal closure response to high VPD-induced water stress are proposed in Supplementary Figure 6. Ions and water flowed into guard cells under low VPD conditions and sustained turgor for stomatal openness. Under high VPD-induced water stress, ABA rapidly accumulated and promoted stomatal closure by altering ion channel activities.

Although stomatal closure prevented excess water loss to maintain physiological processes by passive or active mechanisms, closed “gatekeepers” simultaneously increased stomatal resistance for photosynthetic CO₂ uptake from air to intercellular. The intercellular CO₂ concentration was gradually reduced with VPD elevation (Figure 3). Consequently, the stomatal limitation imposed on photosynthesis gradually became pronounced with VPD elevation (Figure 4). The declines in the leaf water potential and stomatal conductance with VPD elevation were less marked in high-VPD-grown plants in the present research. The distinct response to the VPD is potentially modulated by physiological acclimation to growth conditions (Fanourakis et al., 2019). In contrary to the previous study, no significant differences were observed in the photosynthetic parameters across VPD ranges between the two examined cultivars. The distinct responses of two cultivars between the previous and present study were potentially caused by the examined VPD conditions. The dynamic VPD response of the present and previous studies were performed under cabinets and greenhouse conditions, respectively.

In addition to the first barrier of stomata, CO₂ movement from intercellular to carboxylation sites is constrained by mesophyll resistance. The present study demonstrated that mesophyll resistance was a significant component of diffusion resistance from air to Rubisco in tomato plants. A strong positive correlation between the mesophyll and stomatal conductance was observed among treatments (Figure 6). Similar to the stomatal conductance, the mesophyll conductance of tomato plants was also linearly reduced with VPD elevation (Figure 2). Under low VPD conditions, stomatal conductance coupled with mesophyll conductance was high for efficient CO₂ transport to carboxylation sites within chloroplasts. High diffusion conductance facilitated high chloroplast CO₂ concentrations for carbon fixation (Figure 3). With VPD elevation, the CO₂ concentration inside chloroplasts was substantially reduced under high VPD conditions. The limitation of mesophyll conductance imposed on photosynthesis gradually dominated with VPD elevation (Figure 5). Leaf anatomical traits from the substomatal cavity to the carbon fixation site determine the maximum potential of mesophyll conductance (Muir et al., 2014; Xiong et al., 2017; Earles et al., 2018; Han et al., 2018; Carriqui et al., 2019). The LMA is a composite of underlying traits such as the lamina thickness, mesophyll thickness, cell wall thickness, cell shape, and bulk leaf density, which anatomically regulate the mesophyll conductance (Muir et al., 2014). The LMA determines the upper limit on mesophyll conductance. Meanwhile, the LMA is closely linked to abiotic stress tolerance (Xiong and Flexas, 2018; Xiong et al., 2018). Generally, a higher LMA is a good indicator of greater stress tolerance. In the present study, the LMA of high-VPD-grown plants was lower but higher than that of low-VPD-grown plants (Supplementary Figure 4). Long-term acclimation to a high VPD facilitates enhanced drought tolerance to prevent dehydration by regulating the leaf thickness, cuticular permeability, stomatal morphology, and other anatomical features (Fanourakis et al., 2016, 2020). Long-term exposure to a VPD also affects stomatal sensitivity and morphological features such as the stomatal size, density, index, and spacing, which consequently modulate transpired water loss (Fanourakis et al., 2016, 2020). As mentioned above, root water uptake and supply are inadequate to keep pacing with the great atmospheric driving force under high VPD conditions. A higher LMA indicated dense structural traits, which buffered cellular transpired water loss and prevented leaf tissue dehydration under high VPD conditions. However, CO₂ and water transport share pathways through the mesophyll cell walls and perhaps plasma membranes within leaves (Barbour, 2017; Groszmann et al., 2017; Zhao et al., 2017; Drake et al., 2019). Although dense structural traits improved drought tolerance, the resistance of CO₂ diffusion was simultaneously increased. The LMA was negatively correlated with mesophyll conductance in the present study, which is consistent with previous studies (Hassiotou et al., 2009).

In addition to anatomical determinations, biochemical regulations such as ABA, carbonic anhydrase, and aquaporin facilitate rapid mesophyll conductance responses to short-term changing environmental factors (Momayyezi et al., 2020). The mesophyll conductance is negatively correlated with the leaf ABA content in tomato plants, which is in accordance with a previous study (Sorrentino et al., 2016; Qiu et al., 2017). The foliar ABA content rapidly increased upon long-term and short-term exposure to a high VPD, which is in accordance with a previous study (McAdam and Brodribb, 2016). However, the ABA-mediated regulatory mechanism has rarely been reported. CO₂ entering from intercellular to carboxylation sites inside chloroplasts must pass through plasma membranes. The resistance of transport across the membrane accounts for a great proportion of the mesophyll resistance. It is now well established that aquaporins function as water pores for water transport across membranes and play significant roles in maintaining water homeostasis in response to drought and salinity (Qian et al., 2015; Zhang et al., 2019a). There is increasing evidence that some specific aquaporins (which localise to the plasma membrane and chloroplast inner envelope membrane) are permeable to CO₂ and contribute to the mesophyll conductance (Uehlein et al., 2012; Groszmann et al., 2017; Zhao et al., 2017). Similar to guard cells, some PIPs pores mediate CO₂ uptake and water transport across the plasma membrane (Zhang et al., 2021). Therefore, the specific PIPs potentially reconciled the trade-off between carbon gain and water loss in response to VPD-induced water stress. PIPs were sensitive to drought signals and responded rapidly to enclose gating and inhibit activity.

Gating is a general mechanism of membrane-mediated channels for controlling the permeability of water and CO₂. Although the inhibition of PIPs channels prevents water loss under high VPD stress, CO₂ uptake across the membrane is also restricted by the gating enclosure. ABA has been reported as a signal-inducing variation in the aquaporin content and activity (Fang et al., 2019). Therefore, VPD potentially modulated PIPs gating for CO₂ and water permeability via ABA signalling, which contributed to mesophyll conductance and the photosynthetic rate.