Integration of high-throughput phenotyping with anatomical traits of leaves to help understanding lettuce acclimation to a changing environment

Main conclusion The combination of image-based phenotyping with in-depth anatomical analysis allows for a thorough investigation of plant physiological plasticity in acclimation, which is driven by environmental conditions and mediated by anatomical traits. Abstract Understanding the ability of plants to respond to fluctuations in environmental conditions is critical to addressing climate change and unlocking the agricultural potential of crops both indoor and in the field. Recent studies have revealed that the degree of eco-physiological acclimation depends on leaf anatomical traits, which show stress-induced alterations during organogenesis. Indeed, it is still a matter of debate whether plant anatomy is the bottleneck for optimal plant physiology or vice versa. Here, we cultivated ‘Salanova’ lettuces in a phenotyping chamber under two different vapor pressure deficits (VPDs; low, high) and watering levels (well-watered, low-watered); then, plants underwent short-term changes in VPD. We aimed to combine high-throughput phenotyping with leaf anatomical analysis to evaluate their capability in detecting the early stress signals in lettuces and to highlight the different degrees of plants’ eco-physiological acclimation to the change in VPD, as influenced by anatomical traits. The results demonstrate that well-watered plants under low VPD developed a morpho-anatomical structure in terms of mesophyll organization, stomatal and vein density, which more efficiently guided the acclimation to sudden changes in environmental conditions and which was not detected by image-based phenotyping alone. Therefore, we emphasized the need to complement high-throughput phenotyping with anatomical trait analysis to unveil crop acclimation mechanisms and predict possible physiological behaviors after sudden environmental fluctuations due to climate changes. Supplementary Information The online version contains supplementary material available at 10.1007/s00425-022-03984-2.


Introduction
It is estimated that nowadays about 50% of the global crop yield loss is due to climate change which has increased temperature, severity, and frequency of drought, also causing long periods of severe air vapor pressure deficit (VPD), threatening the sustainability of future agriculture (IPCC 2017).
Furthermore, climate change imposes sudden fluctuations in environmental conditions (Cramer et al. 2011;Will et al. 2013) which threaten agricultural production, affecting crop yield, biomass, biochemical composition, and visual quality of plants, thus impacting the marketability of products (Rouphael et al. 2018).
Plants react to environmental changes by implementing a series of morpho-physiological strategies, the extent of which depends on the plasticity of the species in acclimation (Amitrano et al. 2021a;Choat et al. 2018). When subjected to any environmental stress, plant cells activate a signaling Communicated by Dorothea Bartels. cascade, including hormones, secondary messengers, and signal transducers, which converge to regulate the expression of genes that encode for proteins and enzymes involved in the stress response. This leads to a biochemical and physiological reprogramming, e.g., the accumulation of reactive oxygen species (ROS) or the closure of stomata (Zandalinas et al. 2018).
For example, since most crops perceive high air VPD as an "atmospheric water stress" (Vadez et al. 2013;Will et al. 2013), they try to react by modifying stomatal regulation, with consequences on photosynthesis, to maintain the water balance (Amitrano et al. 2021b). However, while adjusting their physiological processes, plants cannot exceed the limits imposed by their anatomical structure (stomatal size and density, organization of the mesophyll, vein density, etc.), which has, therefore, been hypothesized to be a bottleneck for optimal functioning (Korner 2015). Indeed, due to its ability to influence plant tissue development, VPD has been recognized as the "hidden driver" behind morpho-physiological traits of plants growing in controlled environment (Amitrano et al. 2019). For instance, small leaves are more common in dry environments (under high VPD), as plants can react to adverse conditions by reducing transpiration costs (Sack and Scoffoni 2013), and where turgor-driven cell enlargement is reduced. Indeed, the thickness of the boundary layer increases with the size of the leaf, so it would be difficult for larger leaves to reduce the heat loss in a dry environment (Wang et al. 2020). Therefore, the evaporative demand (expressed as VPD) and the water availability to which the plants are subjected are closely linked, influencing not only the entire plant-water relationships, but also the size of the leaves and other morpho-anatomical traits (Scoffoni et al. 2015). For these reasons, the plasticity of water-related leaf traits, such as the ability to limit transpiration in conditions of high VPD, is promising to alleviate the stress, while allowing water saving, and increasing yield under limiting environmental conditions (Wang et al. 2020). Furthermore, narrower leaves allow to save "construction costs" attributed to more expanded vascular and cell-wall fractions, with a competitive advantage in harsh environments, where the availability of carbon and water is already insufficient (Carins Murphy et al. 2012;Niinemets et al. 2007). The morpho-anatomical traits of plants, measurable at the level of single organism (Kattge et al. 2020), depend on the genetic properties of the species and reflect their evolutionary lineage, but are strongly influenced by the surrounding environment Violle et al. 2007). Recently, plant traits have been used in ecology studies and vegetation modelling as they can act as a proxy to predict how plants behave when subjected to different environmental stimuli, how they will affect other trophic levels, and how they will influence the whole ecosystem (Garnier and Navas 2012;Wright et al. 2017). In a climate change scenario, crop research focused on plant traits is becoming essential to achieve the optimization of resource use and maximization of production, ultimately ensuring food security and the livelihood of a growing population, also thanks to technological innovation in agriculture (Ray et al. 2013;Rockström et al. 2009). In this motivating context, highthroughput phenotyping technologies have been developed to characterize plant responses using image analysis in a non-destructive and non-invasive way (Furbank and Tester 2011;Tardieu et al. 2017), and to highlight the coordination of morpho-physiological and metabolic traits (Poorter et al. 2009). The response to abiotic stresses, including water shortage, heat, and salt stress, has been extensively studied in crops using phenotyping platforms that quantitatively analyze geometric-, color-, fluorescence-and NIR-related traits (Dobrowski et al. 2005;Woo et al. 2008;Mehta et al. 2010;Klukas et al. 2014). As an example, Kim et al. (2020) used image-based phenotyping to determine the water use efficiency, water loss rate and transpiration rate of rice cultivars under drought stress and to detect the difference between tolerant and susceptible plants. However, there is still a gap between the whole-plant phenotyping level and the molecular and anatomical level (Costa et al. 2019;Langstroff et al. 2021), therefore, as many traits as possible should be collected to integrate phenotypic data points and facilitate the integration of knowledge also through modelling (Van Eeuwijk et al. 2019).
To the best of our knowledge, no studies have attempted to comprehensively correlate the effects of VPD and water stress on morpho-anatomical and physiological development of crops, considering the combination of high-throughput phenotyping with observation at microscopic level. Therefore, this study proposes the integration of automated image-based phenotyping with in-depth anatomical analysis of plant tissue to detect the early stress symptoms of plants subjected to different environmental conditions and to investigate the anatomy-driven mechanisms of plant plasticity in acclimation, highlighting the structure-to-function relationships. We aimed to evaluate whether the analysis of anatomical traits can represent an added value to highthroughput phenotyping allowing to predict the degree of eco-physiological acclimation of plants to sudden changes in VPD, as influenced by anatomical traits.
Indeed, establishing the appropriate combination of air and water supply would depend on the environmental conditions, such as radiation, temperature, and relative humidity, and it is essential for maintaining the desired photosynthetic performance. To this end, we cultivated green and red 'Salanova' lettuce (Lactuca sativa L. var. capitata) plants in a phenotyping growth chamber in two cycles with low and high VPD and two irrigation levels (well-watered and low-watered) and then subjected them to short-term changes in the VPD condition. The results obtained provided useful information to understand how to improve the performance of lettuce which is the most common vegetable cultivated in controlled environment (Inoue et al. 2021) and could be a starting point for extending the results to field crops to model their responses when subjected to sudden weather fluctuations.
After 10 days, the more uniform 64 green and 64 red plants at the stage of 4 true leaves were moved to the Phenotyping Growth Chamber (PGC) and placed in the Lem-naTec carriers (Fig. 1a-c). To prevent evaporation from the soil and to provide a uniform background for the phenotyping cameras, the soil surface of all pots was covered with a blue rubber mat (Fig. 1a-c). A few days before the experiments, the soil water content corresponding to 100% of field capacity (FC) was determined by weighing soil-filled pots after full watering and after drying for 3 days at 80 °C as reported in Junker et al. (2015).
Once in the PGC, half of the plants were irrigated at 100% FC and were referred to as WW (well-watered) and the other half at 30% FC and were referred to as LW (low-watered). The irrigation took place via an automatic weighing/watering station with a pump that irrigates the lower container, to avoid any interference with the growth of the lettuce.
Two experimental trials were subsequently conducted in the PGC at 2 different environmental VPDs. VPD values were obtained by keeping temperature fixed and accordingly modifying the relative humidity (RH %) as reported in Amitrano et al. (2021a, b): (i) the first trial was carried out at a VPD of 0.7 kPa (Low VPD; L); (ii) the second one at 1.4 kPa (High VPD; H). Average values and amplitude of VPD levels were decided according to the values typically monitored in plant factories and greenhouses (Garcia et al. 2011;Lu et al. 2015;Zhang et al. 2015). Twelve days after the transferring (DAT) of plants in the PGC, the environmental conditions in the chamber were changed and plants were kept for 5 days at the opposite VPD (from Low to High, LH; and from High to Low, LH) ( Fig. 2) to simulate a short term exposure to the opposite VPD and test their short-term acclimation ability, following the approach reported in Amitrano et al. (2021a, b). These experiments focused on lettuce plants at the beginning of the cultivation cycle, because we were interested to evidence possible plant sensitivity to changing environmental factors at an early stage of plant development.

Image acquisition and analysis
All images were automatically captured using the Lem-naTec system. All plants were imaged daily from the top view using three imaging procedures VIS, NIR, FLUO ( Fig. 1d-i). VIS images were acquired in the visible light spectrum (~ 390-750 nm) using a Basler (Basler AG, Ahrensburg, Germany) Pilot piA2400-17gc (RGB) camera with a resolution of 2454 × 2056 pixels. A fluorescence imaging system (FLUO, excitation: 400-500 nm, emission: 520-750 nm) using a Basler Scout scA1400-17gc and (RGB) camera with a resolution of 1624 × 1234 pixels, allowed the quantification of static fluorescence signals of plants.
Near-infrared (NIR) imaging was performed in the wavelength range 1450-1550 nm using a Nir 300 PGE sensor (Allied Vision Technologies Gmb Hformer VDS Vosskühler GmbH, Stadtroda, Germany, monochrome camera) with a resolution of 320 × 254 pixels.
After transferring all plants to PGC, on days 12 and 22 (before and after the switch in environmental conditions), imaging of chlorophyll "a" fluorescence was performed using the FluorCam imaging fluorimeters (Photon Systems Instruments, Brno, Czech Republic).
Measurements of Φ PSII were made after light adaptation of the plants in the adaptation tunnel and after a period of illumination in the FluorCam chamber as indicated in the "Results" section. The duration of the saturating light pulse to induce Fm was 800 ms with intensity of 4100 μmol m −2 s −1 (white light).
After 1 h of lights off, the dark-adapted plants were subjected to eight pulses of saturated light (800 ms; 4100 μmol m −2 s −1 white light) over the course of 145.8 s. After the first pulse of saturating light, the light was turned off (5 s) to measure F 0 , followed by the second pulse of saturating light and actinic light to measure maximum fluorescence (F m ) followed by 10 s of dark relaxation to measure the Kautsky kinetics.
For the F m quenching analysis, white light and actinic light were turned on for 120 s and supplemented by six saturating light pulses every 20 s followed by further 10 s dark relaxation measurements.
Image analysis was performed using the IAP (Integrated Analysis Platform) software 9 in different steps: image acquisition, pre-processing, feature extraction and post-processing. From the analysis of the images, the most important traits were extracted and classified as geometric (TA in the RGB channel), color-related (Y2G, Lab_a, Lab_b, R2G, hsv_h and hsv_v in the VIS channel and hsv_s in the FLUO channel), physiological-related (Int in the NIR channel) and (Fv/Fm, ΦPSII, NPQ in the chlorophyll fluorescence channel). All the above-mentioned measured traits and their biological meanings are summarized in Table S1 and reported in Tables 3 and 4 and in Figs. 1, 5 and 7.

Biomass traits
Plants were harvested at 23 DAT and determination of fresh (FW) and dry (DW) weights of above-ground biomass was performed at the single plant level. After weighting the biomass for fresh weight, the dry biomass was determined after drying all the plants in oven at 80 °C for 72 h (until stable weight).

Light microscopy and quantitative leaf anatomical traits
On harvest days (16 days after transferring into the PGC, still in the early stage of the cultivation cycle), a single fully expanded leaf from 6 plants per treatment was collected and fixed in F.A.A. solution (40% formaldehyde, glacial acetic acid, 50% ethanol, 5:5:90 by volume) and transferred to the PWA (plant and wood anatomy) laboratory of the University of Naples Federico II for light microscopy analysis. Each leaf was dissected to obtain sub-samples of about 5 × 5 mm in the median region thus avoiding any possible bias due to comparing portion of leaves with different position, age/degree of development. Samples were dehydrated in ethanol series up to 95% and then embedded in JB4 acrylic resin (Polysciences, Germany). Samples were cut into thin cross sections of about . Short-term exposure under the opposite VPD is also represented 5 µm by means of a rotating microtome and stained with 0.5% Toluidine blue in water, as reported in Feder and O'brien, (1968). Sections were observed under the BX51 light microscope (Olympus, Germany) and digital images were collected and analyzed using the Olympus EP50 digital camera and the Olympus CellSens 3.2 software to characterize the leaf lamina by measuring the thickness of upper (UET) and lower epidermis (LET) as well as the thickness of the palisade (PT) and spongy (ST) parenchyma, and the quantity of intercellular spaces, used as a proxy for the compactness of the parenchyma and expressed as percentage of tissue occupied by the intercellular spaces on a given surface (IS %). All measurements were made in the same three regions per cross section, taking care to avoid main veins, and then averaged to obtain a complete view of the traits under investigation.
Stomatal traits were determined on abaxial peels of the lamina, taken centrally in each leaflet, avoiding the main vein and the margin. Peels were analyzed under the abovementioned microscope and five measurements from 3 different peels were performed.
Stomatal density (SD) was calculated as the number of stomata per mm 2 , measured with 20 × magnification, while the stomatal area (SA), expressed in µm 2 , was measured with 50 × magnification in 10 stomata per leaf, considering both the guard cell major (pole to pole) and minor axes to calculate the area of an imaginary ellipse, as reported in Amitrano et al. (2021b). Vein traits on the same leaves were also detected following the protocol by Sack and Scoffoni (2013).
Briefly, the leaves were chemically cleared in in series of dilutions with ethanol (up to 100%) and then double-stained using safranin (1% safranin in 100% ethanol) and fast green (1% fast green in 100% ethanol). Under light microscope, each leaf was imaged in three different areas and used to calculate the vein density (VD) as the ratio between the sum of the lengths of the veins of all orders and the area of the image (Amitrano et al. 2021a).

Statistical analyses
To test the influence of the different independent factors: (i) VPD, (ii) cultivar (C), (iii) water (W) on dependent variables, a three-way analysis of variance (ANOVA) was performed with the IBM SPSS Statistics software (SPSS, Chicago, IL, USA). In case of significant interactions, the data were then subjected to one-way analysis of variance (ANOVA) and the mean values were separated according to Tukey's test with p < 0.05. Line plots, correlation plots (corrplot package with the Spearman's method) and the principal component analysis (PCA; prcomp() function) were all performed using the R software environment for statistical computing and graphics (version 4.4.1).

Growth and biometry
At the end of the growth cycle, significant differences were found between 'Salanova' plants grown under low (L) and high (H) VPD at the two watering levels (well-watered, WW; low-watered, LW). As reported in Table 1, VPD (V), cultivar (C) and water level (W) had a significant effect alone and in combination (VPD × C × W) on plant biomass. Both fresh and dry weight followed the same trend with weight gains at low VPD in both WW and LW (L WW and L LW).
In addition, both weights were enhanced in WW plants with increases of 66-77% compared with LW. No significant differences were found between green (G) and red (R) plants under WW and LW. In Fig. 3a, the growth in term of plant area expansion showed significant differences between the conditions. Plant area increased at low VPD in both WW and LW watering levels. At 12 DAT the switch in environmental conditions only affected plants grown at LW, which underwent a slight bending, to promptly resume their exponential growth (Fig. 3b).
At the end of the growing period, only red plants under high VPD (H WW R) and then subjected to short-term low VPD exposure were able to reach values of plant area comparable with plants grown under low VPD (L WW G and R) but the same did not happen for plant subjected to low-watering levels. 13.57 ± 0.44 b 0.83 ± 0.08 c Significance VPD *** *** C ** ** W *** *** VPD × C × W *** ***

Quantitative anatomy
Microscopy observations showed that the growth in the two VPD conditions and the two different watering levels did not affect the overall organization of the leaf tissues in qualitative terms. Indeed, as shown in Fig. 4, green and red plants under all the condition of VPD and watering levels maintained the typical dorsiventral structure, without any alterations in tissue organization due to the imposed environmental stress. However, the quantification of anatomical traits evidenced the occurrence of significant differences among the treatments. VPD (V), cultivar (C) and water level (W) had a significant effect alone and in combination (VPD × C × W) on the morpho-anatomical traits (Table 2). In particular, the interaction between the factors was significant for all the measured traits, except for UET and ST.
Otherwise, for cultivar alone was significant only for PT (p < 0.01), LET (p < 0.005), ST and LT (p < 0.05). As for the organization of leaf lamina, no significant differences were found with respect to the UET. By contrast, the thickness of lower epidermis was always enhanced under L compared with H VPD plants, except for red plants under low-watered condition (RLW), where no significant differences were observed between L and H plants.
Many differences were found in the thickness of palisade (PT) and spongy (ST) parenchyma. PT followed a different trend under L and H VPD. More specifically, in L VPD plants the thickness increased by 25-29% under LW, compared with WW, whereas in the H VPD plants, PT was reduced by 29-31% under LW. ST followed a similar trend compared with PT with thickness increases in L VPD of 32-36% under LW compared with WW and thickness Concerning stomatal traits, both SD and SA varied greatly between conditions. SD was improved under L compared with H VPD plants, except for the red plants in well-watered condition, where no significant differences were found between L and H plants. Furthermore, SD was always greater under WW than LW with increases in L WW of 15-16% compared with L LW and 32% in H WW compared with H LW. SA was always enhanced in LW compared with WW plants with increases in L of 19-22% and in H of 31-32%. Moreover, SA was higher under H than in L VPD among all treatments.
As for the organization of the veins, VD was always higher in WW plants compared with LW with increases of 25-44%. The highest values were found in both green and red H WW plants and the lowest in LW plants despite the VPD condition.

Principal component analysis and correlation analysis
PCA showed that for phenotypic traits the first two components cumulatively accounted for 78.3% of the total variance, with PC1 accounting for 50.5% and PC2 for 27.8%. PC1 was highly positively correlated with all variables, except for the two variables hsv_s and lab_b, while PC2 was positively correlated with hsv_s, lab_b and hsv_h (Fig. 5).
The PCA scatterplot in Fig. 5a clearly separated the cultivars (i.e., red and green 'Salanova') into two groups, revealing a strong clustering of green and red plants. This clustering is confirmed by the scatterplots in Fig. 5b, c, where each group is divided into sub-groups depending on the watering level (Fig. 5b) and on the watering level in combination with the cultivar (Fig. 5c), whereas no clustering was observed when the data were divided by VPD (Fig. 5d).
The PCA in Fig. 6 showed that for anatomical traits the first two components explained a cumulative 59.9% of the total variance with PC1 accounting for 39.3% and PC2 for 20.6%. Unlike the phenotypic traits, in Fig. 6 the anatomical traits were not separated based on the cultivar (Fig. 6a) but there was a clear distinction into clusters, when data were separated by watering levels (Fig. 6b, c) and VPD conditions (Fig. 6d).
Many positive and negative correlations at different significant levels (p < 0.001 ***, p < 0.005**, p < 0.05*) were found between anatomical, biomass and phenotypic traits in plants grown at low (Fig. 7a-d) and high (Fig. 7e-h) VPD. In all treatments a strong positive correlation (p < 0.001, intensity level = 1) between weights (FW and DW) was always observed. Under L VPD a strong positive relationship (p < 0.001, intensity level = 1; Fig. 7a-d), was found between the palisade thickness (PT) and the lower epidermis thickness (LET), while the same was not observed under high VPD (Fig. 7e-h). A negative correlation was scored between NPQ and QY in all LW treatments (Fig. 7a, b, e, f), showing highest strength (p < 0.001) in H VPD treatments (Fig. 7e, f).
Furthermore, in G plants under both watering and VPD levels (Fig. 7a, c, e, g) it was possible to find a negative correlation between the phenotypic traits yellow to green (Y2G) and total area (TA), while in R plants (Fig. 7b, d, f, h) this correlation was positive. Table 3 reports the phenotyping results of imaging fluorimeters showing PSII performance before and after the change in environmental conditions. As main factors, VPD (V), cultivar (C) and water level (W) always had a significant effect, also showing a significant effect in combination (VPD × C × W) with a p value that was always < 0.01. F v /F m values have always been higher under L VPD despite cultivar (G and R) and water regimens (WW and LW). Overall highest values were found in L WW plants compared with L LW (reduction of 1.2-3.6% compared with WW). LH plants showed reduced Fv/Fm values compared with L (reduction of 2.4%), but in any case, higher or at least comparable with H and HL. H VPD plants showed lowest F v /F m values, especially under LW (reduction of 2.5-3.7% compared with WW). HL plants showed higher values than H plants (approximately 2% increase). As for the differences in cultivars, compared with L WW, G WW plants showed higher values under LH and lower H and HL, whereas no differences were found between green and red plants grown under low VPD and low-watered condition (LGWW and    PSII performance before and after changing environmental conditions. As main factors, VPD (V), cultivar (C) and water level (W) always had a significant effect, except for hsv_s, hsv_v and hsv_h were VPD alone was not significant. Moreover, their interaction (VPD × C × W) also showed a significant effect on all parameters, except for hsv_s, hsv_v and hsv_h and R2G, with a p value < 0.001 in Y2G and Int and < 0.05 in lab_a and lab_b. More specifically, Y2G, lab_b and Int exhibited higher values in green plants grown at high VPD in low-watered conditions. Hsv_s and hsv_v exhibited higher values in green plants at well-watered conditions in both L and H VPDs before and after the change in conditions. As for the R2G trait, higher values were found in red plants under low watered after changing the environment from high to low VPD.

The combination of VPD and water levels influences the development and the morpho-anatomical traits of the leaves.
Changes in VPD and water levels during plant growth strongly influenced biomass allocation and the morphoanatomical organization of the leaves in green and red 'Salanova' plants. In previous high-throughput phenotyping studies (Yang et al. 2009;Dodig et al. 2019), RGB imaging, including plant area and convex hull area, were used to predict plant biomass in many crops (especially corn, barley, and wheat). Here, biomass data in terms of fresh and dry weight measured manually in all plants at the end of the experiment (Table 1) were in good agreement with the plant area estimated by RGB cameras (Fig. 3), showing the highest area and weight (both fresh and dry) in wellwatered plants under low VPD (L WW) and the lowest in low-watered plants under high VPD (H LW). Indeed, it is known that the reduction of VPD in the air by humidification significantly increases the biomass allocation in the leaves and in the fruits (Zhang et al. 2017). Furthermore, a larger leaf area positively affects the transpiration rate during growth, also favoring a balanced absorption of nutrients and a correct water supply to the plant (Pons et al. 2001).
Interestingly, in the present study after the switch in environmental conditions, only the red cultivar under high VPD (HL R WW) managed to reach values of plant area similar to the red cultivar at low VPD, in well-watered conditions (Fig. 3b, c). The short-time exposure to low VPD probably boosted their development in terms of area expansion, but it was not yet sufficient to improve the fresh and dry weight, which in high VPD plants was still lower than the low VPD condition. Increases in plant area and in biomass in the red 'Salanova' cultivar compared with the green one has been reported in earlier studies and explained as the red 'Salanova' reached maturity earlier than the green one (Amitrano et al. 2021b;El-Nakhel et al. 2019).
In the present study, the red cultivar under well-watered conditions also showed a higher stomatal density than the green one (Table 2). Indeed, a positive correlation (p < 0.001) was observed between the stomatal density (SD) and the fresh weight (FW) (Fig. 7d). This probably explains the improvements in leaf area when red plants developed at high VPD, where subjected to short-term low-VPD (more favorable environmental conditions). Higher stomatal density is known to improve carbon gain and maintain carbon gain/water loss homeostasis in plants (Lawson and Blatt 2014). However, in WWR plants between SD and FW there was a tendency toward positive correlation that, however, was not statistically significant (Fig. 7h).
Otherwise, the lack of fresh and dry weight gain in wellwatered red plants at high VPD (H R WW) (together with the area) could probably be ascribed to the greater development of spongy tissue in these plants, compared with other treatments. As also confirmed by the negative correlation (p < 0.05) between SD and DW (Fig. 7h). The spongy parenchyma, in fact, is richer in intercellular spaces with cells that develop slightly thinner cell walls, compared with the palisade parenchyma (Fan et al. 2013). In these plants it was also possible to score a negative correlation between the spongy thickness (ST) and the DW (Fig. 7h). Table 3 PSII performance of green (G) and red (R) 'Salanova' plants grown at the two watering levels (well-watered, WW; low-watered, LW) under Low (L) and High (H) VPD before and after the change in environmental conditions (from low to high, LH; from high to low, HL) Maximum photochemical efficiency (F v /F m ), the quantum yield of PSII electron transport (Φ PSII ), non-photochemical quenching (NPQ). Mean values and standard errors are reported. NS, *; **, and ***not significant or significant at p < 0.05, 0.01, and 0.001, respectively. Different letters correspond to significantly different values (p < 0.05) *** *** *** C *** ** ** W *** *** *** VPD × C × W *** ** ** Overall, in the present study the combination of elevated VPD and water stress induced the development of an anatomical leaf structure in both green and red plants that probably increased the resistance of the mesophyll to water flow and was responsible for a reduced ability in acclimation of these plants when subjected to short-term changes in environmental conditions. For example, a higher percentage of intercellular spaces is a clear sign of greater resistance to water flows in leaves, due to the reduced connection of the cells and the greater presence of airspaces, which increased the length of the path, also reducing the speed of the flows (Amitrano et al. 2019).
Moreover, the lower stomatal density and larger stomatal area could be one of the main reasons of the lower yield of these plants and of the reduced acclimation capacity when subjected to changing conditions, since they did not allow the plants to efficiently regulate the transport and transpiration of water (Dickison 2000). Indeed, the proportion of anatomical traits and in particular of the palisade and spongy parenchyma and of the stomatal and vein density can vary in relation to plant species and habitat, being strictly influenced by the surrounding environment and can affect plant photosynthesis and the entire eco-physiological behavior (Esau 1965).
Here, together with a larger size of the fewer stomata and a higher percentage of intercellular spaces, a lower vein density was also observed in all plants subjected to water stress, irrespective of the VPD condition. The veining of the leaves is essential for water supply; an extended vascular system manages to maintain water balance even when transpiration is high (Brodribb and McAdam 2017;Liu et al. 2016). Different authors proposed that in species with high leaf plasticity, which promptly adapt to changes in environmental conditions, stomatal and vein densities should be coordinated to allow the maintenance of a high physiological efficiency (photosynthesis, conductance, water use efficiency) of the species under sub-optimal environmental conditions (Brodribb et al. 2007;Carins Murphy et al. 2012). However, contrasting results were found for leaf size, vein, and stomatal density under high and low VPDs. Here, at high VPD there was a different response between well-and low-watered plants: while in well-watered plants, the stomatal and vein densities increased, probably to try to compensate for the air drought, in the low-watered plants the reduced stomatal density together with the reduced vein density might be the cause behind the lower biomass production and the reduced capacity in acclimation of these plants when subjected to the changing VPD.
In the present study it was possible to observe statistically significant differences in the PSII performance (F v /F m , Φ PSII , NPQ) between well-and low-watered plants grown at low and high VPD (Table 3). Differently, in previous studies using measurement of chlorophyll "a" fluorescence, F v /F m was relatively insensitive to drought stress. Just to mention a few, Munns et al. (2010) found no differences in the F v /F m time course of both well-irrigated and drought-prone barley plants.
In rice, the F v /F m parameter did not significantly change in the early response to drought stress (Kim et al. 2020). Similarly, in our study the F v /F m values were always above 0.8 (even in LW plants), which is considered a threshold for stressed plants (Ogaya et al. 2011). This may be because plants in the first phase of development showed resistance to water stress, implementing morphological adjustments that led to the reduction of their yield and biomass but allowed them to survive without permanent damage (Yang et al. 2021).
After the switch in environmental conditions, both F v /F m and Φ PSII showed the same trend, with reduction moving from low to high VPD (L to HL) and enhancement in values from high to low VPD (H to LH); however, this trend was more evident in F v /F m, where the differences were always statistically significant (Table 3). Conversely, the NPQ, an indicator of the light energy dissipated in heat to preserve the integrity of photosystem II (Ruban 2016), followed an opposite trend, being higher under LW and high VPD (H and LH).
We attributed these differences in the photosynthetic apparatus to the different organization of the leaf tissue. Indeed, the reduction of stomatal density with increases in water deficit are common in several species (Orsini et al. 2011;Waqas et al. 2017) and together with a reduced opening of the pores are the most common causes of reduction of the photosynthesis and conductance under drought stresses (Domec et al. 2017;Kelly et al. 2016). By contrast, leaf vein traits are even less studied and their development under drought stress has not been investigated much, so far. In addition, in our study, low-VPD plants, especially under limiting water, developed a thicker palisade mesophyll. This trait together with a high stomatal density (compared with H plants) was probably the reason for the maintaining a greater photochemical efficiency even after the switch in environmental conditions (LH). In fact, looking at the related physiological traits, LH plants maintained F v /F m and Φ PSII values higher or at least comparable with H and LH plants.
However, to the best of our knowledge, this different VPD-driven acclimation of lettuces under water stress has never been reported before and is likely to have a major influence on photosynthesis and plant acclimation to changing environmental conditions that cannot be easily detected with RGB cameras.

The limits of phenotyping in detecting the anatomical differences of plants
Non-destructive assessment of plants can provide valuable information on their growth and acclimation to environmental conditions. In this study, phenotypic traits revealed early stress in LW plants, especially under high VPD. For example, the Y2G (yellow to green) increased in low-watered plants under H and LH VPD compared with well-watered plants (Table 4), remarking a stressful condition.
Yellow to green (Y2G) is one of the color-related traits most sensitive to drought stress, which can also reveal symptoms of wilting and senescence (Neumann et al. 2015). In agreement with Y2G, the same plants showed improved Lab_b. High value of this trait indicates the yellow color and Lab_b is, therefore, a proxy for the stress level. Furthermore, several studies (Dodig et al. 2019;Ibraheem et al. 2012) have found a relationship between increase in Lab_b and decrease in dry weight in different cultivars.
The same relationship is reported in this study, particularly for H plants. The main response of plants to drought and high evaporative demand was undoubtedly a reduction in plant area, which likely translates into a decrease in water content of plants. Furthermore, a previous study on rice found a strong correlation between plant area and NIR intensity (indicator of plant water content) (Kim et al. 2020). In agreement with this study, here the Int trait (NIR intensity) was always positively correlated with plant total area (TA) (Fig. 7) and it resulted lower in HVPD and LW plants, indicating lower water content in 'Salanova' lettuces, compared with the other treatments (Table 4).
Also, as expected, the Lab_a trait was improved in R plants as a high value of this trait indicates the red color and small values indicate the green color. Lab_a was also boosted under LW, indicating increased redness of waterstressed plants.
However, as highlighted by the multivariate analysis in Fig. 5 on phenotypic (RGB, NIR, FLUO) traits, the samples grouped by color (Fig. 5a) and by water availability (Fig. 5b), while the different VPD seems to have no effect on the grouping of data in the PCA (Fig. 5d).
Otherwise, by performing a second multivariate analysis on the anatomical traits, it was evident how the samples were clearly grouped by VPD (Fig. 6d) and by water availability (Fig. 6b) and that the influence of color was minimal (Fig. 6a). Thus, bearing in mind that in the present study fluorescence imaging of the chlorophyll "a" revealed many differences between plants grown under different VPDs both before and after the change of conditions, we can, therefore, conclude that these physiological differences in PSII performance were due to anatomical differences between treatments, which cameras (with the exception of chlorophyll "a" fluorescence) were not actually able to detect. Indeed, these cameras could easily detect the stress already in progress in plants but were not able to detect differences in their structural organization, which were responsible for the greater or fewer capacity in acclimation in case of sudden environmental fluctuation.
Therefore, these combined analyses (phenotyping and tissue anatomy) could serve as a proxy for early diagnosis of environmental stress even before symptoms occur (Chen et al. 2012;Humplík et al. 2015;Janka et al. 2018;Tschiersch et al. 2017) and fill the still existing knowledge gap in functional phenomics concerning the relationships between plant phenotype and functional anatomical traits (Kim et al. 2020) which could, therefore, be used as key indicators of photosynthetic efficiency (Baker and Rosenqvist 2004).

Conclusions
In the current context of climate change, very rapid shifts in environmental conditions are expected; therefore, the understanding of the morpho-physiological acclimation mechanisms of crops is necessary to acquire knowledge to improve the cultivation techniques both in the open-field and in controlled-environment agriculture, where sudden climate fluctuation are due to environmental control systems.
In this context, many high-throughput phenotyping experiments have investigated plant responses to various environmental stresses, such as drought, salinity, extreme heat, and frost events (Campbell et al. 2015;Hairmansis et al. 2014). However, to the best of our knowledge, none have used high-throughput phenotyping in combination with microscopy analysis of plant tissue anatomy to study the combined effect of VPD (low and high) and of two different water levels (well-watered and low-watered) especially with the aim of detecting early stress signals in lettuce plants during the early stages of development, focusing on the morphophysiological acclimation of the plant to the change in VPD.
The main findings of this study showed that different VPDs and watering levels determine a different anatomical leaf structure in lettuce, which influenced the extent of crop acclimation. Phenotyping imaging techniques, although effective and reliable, should be complemented with structural or molecular measurements to get the complete picture of crops' plasticity, which is driven by crop anatomical development and influenced by environmental conditions. Indeed, while a combination of RGB, fluorescence and thermal imaging can provide a powerful tool for the early diagnosis of stress symptoms (including air drought and water stress), it is, therefore, essential to support high-throughput phenotyping with the quantification of leaf anatomical traits at specific stages of plant development to forecast the anatomical-mediated ranges for plant physiological acclimation, finally modelling the transient adaptation of the plant to a changing microclimate. Funding Open access funding provided by Università degli Studi di Napoli Federico II within the CRUI-CARE Agreement. Financial support by the Access to Research Infrastructures activity in the Horizon 2020 Programme of the EU (EPPN2020 Grant Agreement 731013) is gratefully acknowledged. Chiara Amitrano's research work was conducted as part of the PhD programme promoted by the Italian Ministry of Education (Programma Operativo Nazionale FSE-FESR "Ricerca e Innovazione 2014-2020", Asse I "CapitaleUmano", Azione I.1 "Dottorati Innovativi con caratterizzazione industriale".

Data availability
The data supporting the findings of this study are available from the corresponding author, upon reasonable request.

Conflict of interest The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.