Leaf structure and ultrastructure changes induced by heat stress and drought during seed fi lling in fi eld-grown soybean and their relationship with grain yield

Studies focusing on terminal drought combined with heat impacts on plants of agronomic value remain scarce, and even less under fi eld conditions. The objective of this study was to investigate leaf structural and ultrastructural changes induced by heat stress (HS) and drought stress (DS) during seed fi lling and their relationship with physiological variables and yield determination. Two soybean cultivars were grown in fi eld conditions. During seed fi lling four treatments were applied, including a control (without manipulation, at ambient temperature and fi eld capacity), HS (episodes exceeding 32°C for 6 h d) during 21-d, DS (20% of fi eld capacity soil water content) during 35-d, and HS×DS. Drought principally reduced leaf area, whereas heat decreased leaf thickness, possible as acclimation strategies, but also irreversible reducing CO2 assimilation sites. Both stresses damaged the outer and inner membranes of chloroplasts, causing swollen chloroplasts and accumulation of plastoglobules, loss of chlorophyll content, and negatively affecting chlorophyll fl uorescence. Thus, the performance and integrity of the photosynthetic machinery were reduced. Through a morpho-functional perspective and a holistic multiscale approach, our results provide evidence of photosynthesis impairment and yield drops under stressful conditions which were associated with structural and ultrastructural (particularly at the level of chloroplasts) modifi cations of


INTRODUCTION
Climate change is one of the main risk factors involved in agricultural production worldwide (Rolla et al. 2018), and therefore threatens global food security. Soybean [(Glycine max (L.) Merr.] is one of the most important commodity crops globally cultivated, and about 90% of the world's production occurs under rainfed conditions, which are often characterized by high temperatures and low or erratic rainfall (Thuzar et al. 2010). These problems are expected to further accentuate in the future (Stocker et al. 2013). Although it is well known that drought stress (DS) along with the occurrence of short episodes of high temperatures (HT) will be important characteristics of future climate scenarios, there are limited studies elucidating their impact on plants of agronomic importance (Sehgal et al. 2019, Mahalingam & Bregitzer 2019, Correia et al. 2018, and even less on soybean under field conditions (Ergo et al. 2018). This is of particular interest since growth, yield, and seed quality are strongly controlled by the environment.
Heat and drought stress during late reproductive development in soybean trigger large negative effects on the crop productivity, decreasing yield (Brevedan & Egli 2003, Dornbos & Mullen 1991. This reduction is in part a consequence of several physiological and biochemical modifications, including a detrimental impact on leaf primary metabolism and redox state (Egli & Bruening 2007, Siebers et al. 2015. Furthermore, in a recent study, we found that heat and drought stress during seed filling, concomitantly induced low levels of primary metabolism variables (i.e. quantum efficiency and yield of photosystem II photochemistry, total leaf soluble sugars, and starch) and elevated levels of oxidative stress (i.e. low values of the ferric reducing ability of plasma, FRAP) and high canopy temperature (Ergo et al. 2018). In addition, leaf chlorophyll dropped causing loss of photosynthetic capacity, depicting stressinduced leaf senescence (Ergo et al. 2018). These changes resulted in yield reduction and highlighted the obvious dependence of grain yield on leaf primary metabolism (Ergo et al. 2018) towards the end of the growth cycle in this species which is strongly limited by assimilating source capacity during grain filling (Borrás et al. 2004).
Biological functions are dependent on the structure and ultrastructure of cells, organelles, and molecules; in consequence, the structural leaf traits changes are major determinants of plant responses to stress conditions. These morpho-functional characters are often overlooked in field studies. Leaf thickness and stomatal density, structural variables widely analyzed under heat stress or drought stress principally in controlled environments have shown dissimilar patterns. For instance, thicker leaves were observed in heated soybean (Djanaguiraman et al. 2011) and droughtstressed olive (Olea europaea L.) (Guerfel et al. 2009, Ennajeh et al. 2010, and common bean (Phaseolus vulgaris L.) leaves (Silva et al. 1999) exhibiting no variation or higher stomatal density (Ciha & Brun 1975, Buttery et al. 1993. The opposite was observed by Makbul et al. (2011) andJumrani et al. (2017) in heated soybean, and by Chartzoulakis et al. (2002) in drought-stressed avocado (Persea americana Mill.). Taken as a whole, none of the previously cited studies have altogether examined the links between the structure and ultrastructure modifications underpinning the mechanisms of the physiological-biochemical traits responses to heat or drought stress conditions and grain yield determination.
The morpho-functional responses to abiotic stress conditions are strongly dependent on the severity, which in turn results from the combination of the duration and intensity of the stress, as well as on the crop developmental stage. Besides, results to date were obtained from experiments carried out in pots under controlled environments, thus arising the old question about how accurate are the observed physiological and structural leaf responses to abiotic stress imposed in controlled environments and how extrapolative are them to field-grown plants, especially in soybean, where it has been shown that leaf structure was more variable and sensitive in plants grown in controlled environments than those grown at the field (Van Volkenburgh & Davies 1977). Ultrastructural changes at the subcellular level under heat and/or drought stress are generally similar and comprise major modifications in chloroplasts (Wahid et al. 2007). Such changes include loss of grana stacking and damage of thylakoid membranes, as indicated by an increased Fo/Fm ratio, leading to significant photosynthesis reductions in wheat (Triticum aestivum L.) (Tambussi et al. 2000) and soybean (Djanaguiraman et al. 2011). However, to our knowledge, a holistic study integrating relationships between leaf photosynthetic performance, structure and ultrastructure, and yield generation in field stressed soybean has not been documented. This knowledge is essential to complement our understanding of the functional relationships under stressful productive scenarios and to identify new traits for breeding programs.
The objective of this work was to study in field-grown soybean, leaf structure and ultrastructure modifications under heat and drought stress during seed fill and their association with physiological variables and grain yield determination through a morphofunctional perspective and a multiscale approach. This objective was supported by the following working hypothesis: the interactive effect of HS and DS cause differential alterations of leaf morpho-functional characters, which are in close relation with the physiologicalbiochemical responses and grain yield reduction.

Plant husbandry
In the growing season 2012-2013, an experiment was conducted at the research station of the National Institute of Agricultural Technology (INTA), located in Manfredi (31°49ʹS, 63°46ʹW), Córdoba province, Argentina. The soybean cultivars used for the experiment were SPS4×4 RR and SPS4×99 RR from Syngenta Company, being both of maturity group IV and indeterminate growth habit. These two cultivars were selected because they are commonly sown across the Argentinean soybean region, exhibit similar cycle length (i.e. crop duration from sowing to full maturity) but differ in yield stability across environments according to unpublished data from Syngenta Development Area. The planting date was 16 Nov. and the plant density was 25 plants m -2 . Seeds were inoculated with commercial strains of Bradyrhizobium japonicum. Plots were four rows four m long with 0.52 m distance between rows. The soil is silty loam Typic Haplustoll (USDA Soil Taxonomy). More details about crop husbandry can be found in Ergo et al. (2018).

Experimental design and treatments
Four treatments were applied at R5.5 phenological stage (Fig. 1), defined as seeds of six mm length in a pod at one of the four uppermost nodes, according to the scale of Fehr & Caviness (1977). Treatments consisted of manipulations of temperature and water availability during seed filling in soybeans. The treatments involved a control (without manipulation, at ambient temperature and field capacity, by drip irrigation), heat stress (episodes that exceed 32°C for 6 h d -1 around noon), drought stress (20% of available water content), and heat stress × drought stress.
Treatments were distributed in a split-split plot design with water regimes, cultivars, and thermal regimes in the principal plot, sub plot, and sub-sub plot, respectively. Two replicates were always used. Water regimes resulted in two types of plots: non-drought-stressed (NDS) and drought-stressed (DS). In NDS, soil water was kept near field capacity during the entire growing using drip irrigation; whereas in DS, the target available soil water content from R5.5 to R8 phenological stage (when 95% of the pods have reached their full mature color, i.e. full maturity) was 20%. To avoid rainfall water incomes treatment areas were emplaced under independent mobile rainout shelters (8 m long, 5 m wide, and 6 m height) covered with transparent polyethylene film (100 μm thickness). The mobile rainout shelters were only in place during precipitation events. The evolution of soil moisture content (percent of fi eld capacity) was surveyed up to 200 cm depth on each plot at 20 cm depth intervals using the gravimetric method. Gravimetric water content was measured four times in control plots and DS plots: i) before the onset of the stresses (at R3 phenological stage, defi ned as pods having 5 mm long at one of the four uppermost nodes), ii) at the onset of the stresses (at R5.5 phenological stage), iii) at the end of the heat stress treatment (i.e. 21 d from R5.5 phenological stage), and iv) at full maturity (at R8 phenological stage). In those plots exposed to heat stress treatments (irrigated HS and HS combined with DS) soil water content was measured at two moments: at the end of the heat stress treatment, and full maturity. Details of this determination and soil available water content calculation were given in Ergo et al. (2018).
Thermal regimes also resulted in two types of plots: non-heat-stressed (NHS) and heatstressed (HS). In the NHS, plot temperature was kept at ambient levels, whereas in the HS plots, plants were exposed to brief episodes (6 h) of high temperatures (above 32°C) during 21 d from R5.5 (Fig. 1). Heat treatment was imposed through transparent polyethylene fi lm made covers (100 μm thickness) fi xed to rigid metal shelters (4 m long, 2 m wide, and 2 m high) opened at the Figure 1. Scheme indicating soybean phenology and the duration of heat stress (black line) and drought stress (red line) treatments. End of seed fi ll in drought-stressed plots and heat-stressed combined with drought-stressed plots = 106 days after emergence (DAE); end of seed fi ll in irrigated heat-stressed plots and control plots =116 DAE. Samples for leaf structure, ultrastructure, leaf area, and specifi c leaf weight determinations were collected 95 DAE (coincidently with the end of heat stress imposition). Indicated developmental stages (Fehr & Caviness 1977) are V1, V2, V3, VN: 1st, 2nd, 3rd, Nth node; R3: beginning pod (pods in one of the uppermost nodes is 0.5 cm). ; R5.5: mid grain fi lling (seeds of 6mm length in a pod at one of the four uppermost nodes); R7: beginning of physiological maturity (fi rst pod on the main stem is matured to a brown pod color); R8: full maturity (when 95% of the pods have reached their full mature color). Illustration: Oliveira Junior et al. (2016). bottom (0.60 m above ground level) to enable gas exchange. These covers remained closed only 6 hours a day (from 1000 h to 1600 h) during the days of treatment, after which the polyethylene film was rolled completely on the middle of the shelter allowing plants of HS plots to be exposed to ambient conditions.
The polyethylene enclosure generated a greenhouse effect promoting temperature increases. Button-type digital loggers (DS1923L-F5, resolution: 0.5°C, I-buttons data loggers, Digi-Key Co. Ltd., USA) located in the center of each shelter at the top of the canopy (at 1.2 m from the ground level) recorded air temperature and relative humidity (RH) inside each shelter every 30 min during the treatment period. The degree of heat stress was computed as the difference in mean temperatures between control and heated plots. The framework proposed by Neiff et al. (2016) was used to calculate heat-stressful temperatures (HST, °C h -1 ), with a modification of the optimum temperature value, which was set at 32°C. This adjustment takes into account that above this temperature both growth rate and the duration of seed filling are reduced in soybean (Egli & Wardlaw 1980, Gibson & Mullen 1996. During the daily hours of heating treatment, the average vapour pressure deficit (VPD, kPa) was calculated based on the framework developed by Allen et al. (1998). More details about the heat stress characteristics in terms of temperature and VPD, as well as the average daily maximum, minimum and mean air local temperature (obtained from INTA Manfredi Meteorological Station) during the seed filling, can be found in Ergo et al. (2018).
Measurement of light intensity (n=3) was accessed above the canopy inside and outside the shelters with a one m long line-quantum sensor (Cavadevices, Buenos Aires, Argentina) between 1100 and 1400 h on clear days. Additionally, we had also measured the red/far-red ratio (R:FR) using the Skye SKR 110 Red/Far Red sensor (Skye Instruments Ltd, Llandrindod Wells, UK) at noon during sunny days. Measurements of R:FR (n=3) were performed 10 cm above the canopy with the sensor parallel to the soil surface. During the period of heating treatments CO 2 was measured (n=3) from a buffer box placed at 0.4, 0.8, and 1.2 m height above ground using a portable, openflow gas-exchange system LI-6400 (LI-COR).

Data collection
Leaf structure, ultrastructure, leaf area, and specific leaf weight determinations were made on the central fully expanded leaflet of the third trifoliate leaf from the main stem apex of plants from the central rows of the control, HS, DS, and HS × DS plots. Samples were collected 21 d after the plot reached R5.5, coincidently with the end of heat stress imposition ( Fig. 1), to avoid recovery processes (especially in irrigated heated plots) and also because drought stress intensity was sufficiently advanced (6.48±11.26% in DS and 4.30±2.68% in DS×HS plots, of soil available water content), as shown by the soil water gravimetric data (see Fig. 2 in Ergo et al. 2018). It is important to highlight that the stress treatments were imposed during the late reproductive period of soybean, a crop characterized by accelerated senescence processes and loss of leaves.
For the structure and ultrastructure determinations (except the ones processed for the transmission electron microscopy), leaflet discs of tagged trifoliate leaves were collected and immediately fixed in FAA (formalin, acetic acid, and ethanol) until further procedures. The fixation technique for observations with transmission electron microscopy is described below.
Leaf clearings: Stomata count was conducted in semi-permanent slides based on the technique of diaphanization proposed by Payne (1969). Ten random samples of adaxial and abaxial epidermis were obtained. Leafl et discs were cleared and stained with safranin. Stomatal index was calculated based on the formula proposed by Salisbury (Wilkinson 1980) (Eqn. (1)): Where I s is the stomatal index, D s is the stomatal density, and D ec is the epidermal cell density is equal to the number of stomata and epidermal cell per unit leaf area (1 mm -2 ), respectively; each stoma consists of the stomatal pore and two fl anking guard cells.
Light microscopy: For the preparation of histologic sections, the material was processed by dehydration through an ethanol series with a pre-impregnant rinsing of tertiary butyl alcohol (Gonzalez & Cristóbal 1997) and infi ltration in Histoplast® paraffin (Biopack, Buenos Aires, Argentina), according to Johansen (1940). Leaves were sectioned transversely (10-12 μm thick) with a Microm HM 350 rotary microtome (Microm International, Walldorf, Germany). Sections were stained with astra blue-safranin (Luque et al. 1996) and were permanently mounted with synthetic Canada balsam (Biopur, Buenos Aires, Argentina). The cross-sections were made at the level of the main vein in the middle part of the leafl ets. A total of 15 measurements of Figure 2. Biplot displaying relations between soybean leaf structure (black circles) and physiological traits (white circles). Data correspond to the average of two soybean cultivars grown under four conditions during seed fi lling (grey triangles): control (non-heat-stressed and non-drought-stressed), heat-stressed (HS), drought-stressed (DS), and HS × DS plots. Measurements were performed on the central leafl et of the third trifoliate leaf from the main stem apex, 21 d after the plot reached R5.5 (coincidently with the end of HT imposition), except the canopy temperature which was measured also during the HS treatment. ᶲPSII= quantum effi ciency of photosystem II photochemistry; Fv/Fm= maximum quantum yield of photosystem II; Fo/Fm= thylakoid membrane damage, SPAD= leaf chlorophyll content, AD= adaxial, AB= abaxial.
the thickness (µm) in the upper and lower epidermis and cuticle, parenchymal palisade and mesophyll, were performed in each experimental unit. The measurement of the cuticle thickness was made using the ultrafine sections (1 µm), derived from the transmission electron microscopy technique (explained below) (N=20). Total leaf thickness (TLT) was calculated as the sum of the thickness of the cuticle, the epidermis, and the mesophyll. Anatomical features were observed and photographed using a Leica MZ6 stereomicroscope and a Leica DM LB2 compound microscope (Leica, Wetzlar, Germany), both equipped with a digital camera (Leica digital, Wetzlar, Germany).
Scanning electron microscopy: The fixed leaf discs were dehydrated through a series of increasing acetone solutions. The material was then critical point-dried with solventsubstituted liquid carbon dioxide and sputtercoated with gold-palladium. Photomicrographs were obtained using a JEOL 5800 L Vat 20 kV (JEOL USA, Peabody, MA, USA).
Transmission electron microscopy: Small leaflet discs of 2 mm diameter were pre-fixed in 1% glutaraldehyde, 4% formaldehyde in phosphate buffer (pH 7.2) for 2 h. The samples were then post-fixed in 1.5% OsO 4 at 2°C in the same buffer for 3 h. Then, they were dehydrated in ascending graded series of acetone and embedded in Spurr's resin. Sections of 1 mm thick were made on a Reichert-Jung ultramicrotome and stained with toluidine blue. Ultrathin sections of 70 nm were made using an ultramicrotome (Reichert-Jung, Vienna, Austria) and stained with uranyl acetate and lead citrate (Zarlavsky 2014). The sections were examined under a JEOL-JEM 1200 Ex II TEM at 85 kV (Jeol Ltd., Tokyo, Japan). The photomicrographs obtained were used to assess the ultrastructure of the leaf stomatal apparatus, as well as chloroplasts of palisade mesophyll and bundle sheath cells. Chloroplasts dimensions in mesophyll cells as well as in the bundle sheath cells (i.e. width and length) were measured, and the size of the chloroplast expressed in μm 2 was estimated. Detail of chloroplast number analyzed in mesophyll cells were 14 in control, 9 in HS, and 12 in DS, as well as in HS × DS plots, whereas in bundle sheath cells were 20 in control, 9 in HS, 10 in DS, 7 in HS × DS. In all cases, the sample size was reduced to the optical field.
Leaf area (LA, cm 2 ) determination was made through digital photography processed with Image-Pro plus (Ipwin32, Media Cybernetics, Inc., Rockville, USA) (N=12). Specific leaf weight (SLW, mg cm -2 ) was determined on 12 leaf discs of 2 cm 2 taken from each of six collected leaves. Leaf discs were placed in an oven at 70°C for 72 h, until constant weight. Dry weights were computed and SLW was calculated using the following formula: SLW = LDW/LA, where LDW is the leaf dry weight and LA is the leaf area of the dried leaves discs.
Relative chlorophyll content, chlorophyll a fluorescence, and canopy temperature measurements were made during the SF on the central leaflet of the third trifoliate leaf from the main stem apex of six plants per plot from the central rows of the control, HS, DS, and HS × DS plots. Relative chlorophyll content was measured with a portable self-calibrating chlorophyll meter (Minolta SPAD-502, Osaka, Japan) five times during the SF. Measurements of chlorophyll a fluorescence were performed on cloudless and windless (< 5 km h -1 ) days and between 1100 h and 1400 h using a pulse amplitude modulated fluorometer (FMS2, Hansatech Instruments, Pentney King's Lynn, UK) three times during the SF. Maximum quantum yield (Fv/Fm = (Fm -Fo)/ Fm) and quantum efficiency of photosystem II photochemistry (ᶲPSII = (Fm'-Ft)/ Fm') were measured at the upper surface on six plants per plot, being the leaves dark-adapted (30 min) for Fv/Fm measurements. Where Fm' represents the maximal fluorescence yield of an illuminated leaf, Ft is the steady-state fluorescence yield, and Fo and Fm are the minimal and maximal fluorescence yields of dark-adapted leaves, respectively. These two last variables were used to calculate the ratio Fo/Fm, indicative of the thylakoid membrane stability. Canopy temperature was measured using a hand-held infrared thermometer (Testo 845, Barcelona, España) two times throughout SF (CT SF ). Readings were made at an angle less than 45 degrees to the horizontal plane integrating many leaves. Five readings were taken from each experimental unit between 1200 h and 1500 h, only on clear, sunny days with minimal wind.
Seed yield: Yield was measured after physiological maturity by hand-harvesting a bordered area of 1.13 to 2.08 m 2 in each replicate. Pods were threshed, and seeds were weighed and counted. Seed number (SN) per square meter (seeds m -2 ) was calculated as the quotient between seed weight per square meter (g m -2 ) and individual seed weight (SW mg seed -1 ), both expressed on a dry moisture basis.

Data analyses
Statistical analysis were performed using InfoStat (Di Rienzo et al. 2015). To access a global and integrated view of all variables determined, a principal component analysis (PCA) was performed (Johnson & Wichern 2002). This multivariate statistic technique for dimensionality reduction allowed us to explore associations between yield, leaf physiological, and structural traits taking into account several factors simultaneously (temperature and water regimes), better representing the real situation in crop systems. In our previous study, the effect of cultivars or the interaction between cultivars and heat or drought stress were not statistically significant for yield, SN, and SW (see Table II in Ergo et al. 2018). Thus, for further analyses the effect of cultivars was not considered as a factor, instead, it doubled the number of repetitions from two to four. A two-way ANOVA was conducted to test the individual and interactive effects of heat and water stress at the 5% significance level (P≤0.05) for all parameters evaluated. Whenever the ANOVA test indicated a significant difference, a pairwise comparison of means by Fisher's least significant difference (LSD) test (Sokal & Rohlf 1995) was performed. InfoStat statistical software was used for the analysis (Di Rienzo et al. 2015).

Seed filling environmental conditions and grain yield components
Average daily mean air temperatures during the 21 d HS treatment were 27.5°C ± 1.0, 27.4°C ± 0.9, and 22.5°C ± 0.7 for HS, HS × DS, and control plots, respectively. Average daily maximum air temperatures in the same period were 35.8°C ± 1.5, 38.0°C ± 1.4 and 28.1°C ± 1.0 for HS, HS × DS, and control plots, respectively. Daily heating during six-hours increased the average maximum air temperature from 31.6°C ± 1.3 in control plots to 35.5°C ± 0.8 and 37.8°C ± 0.7 in HS and HS × DS plots, respectively. During the sixhour heating the average mean air temperature in HS and HS × DS plots was also higher, 3.2 and 5.0°C, respectively than the control (28.1°C). This resulted in a slightly higher intensity of HS (quantified as HST) for HS × DS plots compared to HS plots (13.1°C h −1 ± 0.5 vs 10.7°C h −1 ± 0.5).
Artificial heating increased the VPD in HS (2.3 kPa ± 0.4) and HS × DS plots (3.5 kPa ± 0.5) during the treatment period compared to control (1.6 kPa ± 0.2). Noticeably, differences in this variable between heated and control plots were entirely attributed to higher temperatures during the heating period since the relative humidity inside the polyethylene shelters was not affected (see Table I in Ergo et al. 2018). The light intensity was slightly reduced by the HS semi-closed polyethylene shelters in a magnitude equal to 12.9% in HS and 13.2% in HS × DS compared to ambient values (1349.6 μm m -2 s -1 ). Whereas, during the heating treatment the R:FR above the plant canopy was 1.11 ± 0.02 for control, 1.08 ± 0.02 for HS, and 1.08 ± 0.03 for HS × DS, being the differences in the R:FR between control and heated shelters plots negligible (lesser than 3%). The CO 2 concentration (μmol mol −1 ) at the top of the plant canopy (1.2 m above ground) was 379.6 ± 4.5 for control, 373.7 ± 6.6 for HS, and 375.7 ± 14.3 for HS × DS. At the middle of the plant canopy (0.8 m above ground), the CO 2 concentration (μmol mol −1 ) was 378.9 ± 2.6 for control, 369.6 ± 10.7 for HS, and 373.8 ± 6.5 for HS × DS. At the bottom of the plant canopy (0.4 m above ground) the CO 2 concentration (μmol mol −1 ) was 382.5 ± 3.5 for control, 373.6 ± 11.4 for HS, and 382.9 ± 5.6 for HS × DS. As it emerges from the registered values, the differences in the CO 2 concentrations between control and heated (i.e. inside the semi-closed polyethylene shelters) plots during the heating period were less than 4%. The volumetric soil moisture contents measured up to 2 m depth in irrigated plots at the onset of treatments (R5.5 phenological stage) for both control and HS plots were similar (75 ± 10%) and this condition was maintained until physiological maturity ending with 84.0 ± 6.01% and 82.6 ± 3.02%, respectively. In plots under water deficit, soil available moisture content started to decrease and was about 16-20% at R5.5, thereafter it declined throughout the stress period ending with 2.58 ± 7.82% in DS plots and 3.02 ± 8.79% in HS × DS plots, causing a severe terminal drought stress (see Fig. 2 in Ergo et al. 2018).
In our previous work, we showed that HS combined with DS reduced yield and its components (SN and SW). Yield and SN were also reduced in irrigated HS treatments; however, SW was similar to that exhibited by the control. This highlighted the capacity of soybean crops to benefit from improved environmental conditions late during the reproductive stage through SW modifications (Ergo et al. 2018).
Heat and drought stress effects on leaf tissue structure and ultrastructure: correlations with physiological and yield traits The biplot obtained from the PCA indicated that the first two principal component vectors (PC1 and PC2) explained a high proportion (92.3%) of the total variation of the data (Fig. 2). Thus, reflecting the fundamental patterns of the relationships among the measured traits. The correlation coefficient between any two traits is approximated by the cosine of the angle between their vectors. Acute and obtuse angles indicate a positive and negative correlation, respectively. One of the most prominent relations revealed by this biplot (Fig. 2) was the strong positive associations among yield, ᶲPSII, and Fv/Fm indicated by the acute angles among variables, being the three vectors oriented towards irrigated control plots. Likewise, SN, specific leaf weight, and the dimensions of almost all cell layers contributing to total leaf thickness were positively associated and oriented also towards irrigated control plots. In contrast, seed weight, leaf area, and SPAD trait vectors were in the direction of the irrigated HS plots. In addition, the biplot also indicates a negative correlation of yield, SN, and SW with epidermal cell density, stomatal density and stomatal index (adaxial and abaxial), canopy temperature during seed filling (both during and after HS treatment), thylakoid membrane damage (Fo/ Fm), mesophyll chloroplasts and bundle sheath size which were located in the opposite side of the plane towards DS plots alone, or combined with HS (Fig. 2). Tables I and II show the average and standard error for each of the 17 analyzed leaf structure traits. The information contained in both Tables supports the correlations between leaf morpho-anatomical features and physiological traits under different treatments observed in Fig. 2. Leaf area (LA) was significantly decreased by drought (22% in DS and 25% in HS × DS) compared to well-watered treatments (Table I). In contrast, the specific leaf weight (expressed as unit weight per unit leaf area) exhibited greater and significant decreases under heated treatments independently or in combination with DS compared to the individual DS. In these plots, SLW was reduced by 27 and 30%, respectively compared to controls (Table  I). Overall, HS, DS, and HS × DS treatments significantly decreased the thickness of almost all cell layers contributing to total leaf thickness ( Table I). The effects were more pronounced in the irrigated plots exposed to brief periods of HS, where the thickness of the adaxial and abaxial epidermis, palisade, mesophyll, and cuticle adaxial and abaxial were decreased by 17, 15, 29, 30, 31, and 35%, respectively (Table I).
In drought-stressed plots alone, or combined with HS significant decreases of the mentioned layers were also observed but of less magnitude compared to irrigated heated plots. As a result, these changes lead to a noticeable reduction of the TLT being 29% thinner in HS and 13% thinner in DS plots alone or combined with HS compared to control plots (Table I). Interestingly, across treatments, the different cell layers preserved the proportion representing the epidermis (AD and AB), the mesophyll, and the cuticle (AD and AB) between 10-12%, 85-88%, and < 3%, respectively of the TLT (Table I).
The reduction observed in LA was reflected in the significant increase of the density of abaxial and adaxial epidermal cells. The stomatal density and index showed a similar trend, being significantly increased by drought compared to controls in both surfaces of the leaf, with the only exception being the upper stomatal density and index which were similar to controls (Table II). Although soybean leaf is amphistomatous (stomata are present on both sides), and as expected from a mesophytic species, the proportion of stomata per unit area was greater in the abaxial than the adaxial side, and this trend was maintained across treatments (202, 235, 249 and 170% more stomata mm -2 abaxial than adaxial in control, HS, DS, and HS × DS, respectively). Greater increases of all the mentioned variables (adaxial and abaxial) were observed, under the combination of both stresses, coincidently with the higher leaf area reduction in this treatment than in the individual DS plots, although these differences were not statistically significant (Table II). The light micrographs of leaf anatomy in cross-sections showed a diminution of the epidermis adaxial and abaxial, palisade parenchyma, mesophyll, cuticle adaxial, and abaxial, resulting in a notable reduction of the TLT in HS (Fig. 3c-d), DS (Fig. 3e-f), and the combined stresses (Fig. 3g-h) compared to control (Fig. 3a-b) ( Table I). The midrib of control leaves showed turgent parenchyma cells (Fig. 3a) and the dorsiventral mesophyll exhibited three to four cellular layers of elongated columnar palisade parenchyma and three to four of spongy parenchyma cells, which were well-rounded or well-lobed shaped, thereby allowing gases to circulate through abundant air spaces (Fig. 3b). Under stressful conditions, the most noticeable changes involved cell size and their distribution patterns. For instance, leaves from irrigated plots exposed to brief periods of HS displayed more obliterate and smaller parenchymal cells at the level of the midrib (Fig. 3c), as well as larger intercellular spaces in the mesophyll, in both palisade and spongy parenchyma (Fig.  3d) compared to the control ( Fig. 3a-b). It is noteworthy, that under DS alone or combined with HS, cells of palisade decreased in diameter and were more compactly arranged, and cells of spongy parenchyma were smaller (Fig. 3e-h) compared to those from control.
Overall, across treatments, the scanning electron micrographs of leaf epidermis showed that the stomatal apparatus were predominantly paracytic or rubiaceous being both subsidiary cells of different sizes (Fig. 4a-h). Micrographs from control leaves (Fig. 4a-b) disclosed a larger size of epidermal cells relative to those observed in leaves under independent or combined stresses (Fig. 4c-h), leading to a greater density of both epidermal and stomatal cells under these treatments (Table II). A detailed level of the impact of the stresses over the stomatal apparatus was obtained by ultrastructure analysis of the guard cells. Relative to control (Fig. 5a) these cells displayed clear membrane damage in chloroplasts, loss of grana and stroma thylakoid and plasma membrane, as well as an increase of the cytoplasm vacuolization when plants were subjected to heat, drought, and the combined stress (Fig. 5b-d).
Ultrastructural images of control mesophyll cells showed chloroplasts with intact outer membrane, well-formed stacked thylakoids and stroma thylakoids, several grains of starch, and a few plastoglobules (Fig. 6a). Remarkable changes were observed in chloroplasts ultrastructure when analyzing micrographs obtained from HS, DS, and HS × DS (Fig. 6b-d). Under single HS and DS, chloroplasts showed notable unstacking of membranes and accumulation of starch and plastoglobules (Fig. 6b-c). Chloroplasts subjected to the combined HS × DS exhibited reduced and less stacked thylakoids, as well as decreases in starch grains and abundant plastoglobules (Fig.  6d). Moreover, the number of plastoglobules per chloroplast increased 43% in HS, 471% in DS, and 441% in HS × DS compared to chloroplasts of control micrographs (0.88 ± 0.64 plastoglobules per chloroplast). The chloroplasts of the bundle sheath cells in the control micrographs (Fig. 7a) were smaller than those of cells affected by the three stress treatments, which also revealed severe membrane damages similar to those observed in the chloroplasts of mesophyll cells (Fig. 7b-d). Heat stress, drought, and their combination also caused severe damage to the mitochondria outer and inner membranes resulting in a diffuse outline, in mitochondria of the mesophyll cells (Fig. 6b-d) as well as of the bundle sheath cells (Fig. 7b-d).   Heat and drought stress effects on mesophyll chloroplasts and bundle sheath cells size, thylakoid membrane damage, and chlorophyll content Significant (p = 0.0009) increases of mesophyll chloroplast size were observed under drought stress both alone (37%) or in combination with episodes of HS (48%) (Fig. 8a). The independent or combined effect of HS and DS was even more pronounced in the case of the chloroplasts of the bundle sheath cells, which exhibited a large size increase of 80, 92, and 70% in HS, DS, and HS × DS, respectively (Fig. 8b).
During the first 10 days of treatments (from day 75 to day 84 after emergence), there were non-detectable changes in the thylakoid membrane integrity as indicated by Fo/Fm ratio, which remained unaffected (Fig. 9a). However, after this 10 d period, the damage of thylakoid membranes significantly increased under DS and HS × DS. Furthermore, measurements of Fo/Fm ratio 29 d after the onset of the stresses (day 103 after emergence) indicated significantly (p = 0.0001) increases of thylakoid membrane damage by 21% and 39% in DS and HS × DS, respectively compared to control plots (0.229) ( Fig. 9a). Across treatments, chlorophyll content exhibited an increasing trend from day 68 (6 d before the onset of HS treatment) to 83 days after emergence (10 d after the onset of HS treatment), i.e. after seed fill has begun, being this increase 28% for control and HS, 21% for DS and 25% for HS × DS. Afterward, chlorophyll content decreased in different magnitude depending on the treatment. For instance, 10 days after finishing the HS treatment chlorophyll content decreased 28% and 23% in irrigated control and HS plots, respectively, being greater the decrease caused by DS separately (51%) or in combination with HS (54%) (Fig. 9b). Thus, plants from DS and HS × DS plots ended up with less chlorophyll content, 37% and 43% less compared to plants grown under control conditions (34.9 SPAD value).

DISCUSSION
In this research, we studied how leaf structure and ultrastructure of field-grown soybean respond to heat stress, drought, and their interactive effect during seed filling and their association with physiological variables and grain yield determination. Thermal and water regimes were manipulated through heat and rainout shelters both made of polyethylene film. During the heating hours, the semi-closed polyethylene shelters slightly reduced CO 2 concentration and did not increase relative humidity compared to ambient levels. Actually, the RH was 54.6% and 45.7% in HS and HS × DS, respectively, whereas in ambient plots it was 67.2% (Ergo et al. 2018). The reductions in light intensity and the R:FR ratio due to the semi-closed polyethylene shelters were negligible and did not produce confounding effects. Although the light intensity and R:FR ratio have been shown to affect photosynthetic metabolism and chlorophyll content (Yang et al. 2018, Wherley et al. 2005) no significant effects on leaf photosynthetic performance (Bowes et al. 1972) as well as yield and its components (Wahua & Miller 1978) were detected in soybean exposed to 20% of light intensity reductions. Besides, 35% in light intensity reductions were reported to not affect the percentage of open stomata and stomatal aperture length, maintaining these plants a higher chlorophyll a, b, and total chlorophyll concentration compared to control plants (Zhang et al. 2016).
Most of the available information about the physiological, structure, and ultrastructure modifications in response to HS or/and DS have come mainly from studies conducting experiments under controlled environments, difficulting the extrapolation to field-grown plants. For instance, growth chamber and phytotron experiments may not provide a realistic approach, since i) the same temperature regime is often imposed throughout the whole period of interest, resulting in more severe temperature treatments than similar temperature regimes in the field where temperatures vary diurnally (Egli et al. 2005), and ii) plants are grown in pots which limit root depth, and thus induce a rapid rather than gradual rate of stress development. To our knowledge, there are currently no studies of the interactive effect of HS and DS on the relationship between plant structure and function in soybean plants under field-cropped conditions. In recent research, we showed an integrating approach combining physiologicalbiochemical parameters (describing the primary metabolism and redox state) to explain yield variations in field-grown soybean under HS and DS conditions (Ergo et al. 2018).
The plasticity in leaf area is an important feature by which plants control the lightharvesting capacity as well as the gas exchange surface, and consequently the photosynthetic potentials given by the relationship between the carbon and water economy (Ting 1982). Cell elongation growth is essentially motorized by water entering into the cell and it is the most sensitive process to water deficit induced by drought or high evaporative demand given by heat stress (Boyer 1970, Blum 1996. Thus, the observed increased densities of both epidermal and stomatal cells could be a consequence of cell expansion reduction (as indicated by the decrease in LA) in response to DS alone or combined with HS, indicating that acclimation to drought occurred. In line with these results, the biplot highlighted a strong negative correlation between LA and the mentioned traits, as displayed by the obtuse angles formed between their vectors (Fig. 2). Furthermore, the microphotographs of leaves from DS treatments alone or combined with HS ( Fig. 4e-h) clearly showed a decrease in the size of the epidermal cells, which in turn promoted closer packing of stomata. Earlier studies in field-grown soybean found that drought-stressed plants had significantly greater stomatal density and smaller leaf area than non-drought-stressed ones (Ciha & Brun 1975). Matching responses were observed later not only in soybean (Buttery et al. 1993) but also in other species as well, such as olive (Guerfel et al. 2009, Bosabalidis & Kofidis 2002, wheat , and in Lotus creticus L. (Bañon et al. 2004) when subjected to drought compared to well-watered plants.
Unlike stomatal density, which is a function of both stomatal initiation and epidermal cell expansion, the stomatal index normalizes for the effects of the latter (i.e. density of epidermal cells), being affected only by stomatal initiation as described by Royer (2001). Thus, even when several studies support the positive correlation between DS and stomatal density; this response seems not to be universal for the stomatal index as shown by Estiarte et al. (1994), who reported no significant changes in SI in wheat plants grown in dry treatments compared to those from the wet ones. Likewise, Clifford et al. (1995) reported drought-induced large increases in SD in ambient air (35% on the adaxial and 28% on abaxial surfaces) with a concomitant no significant effect on SI. In our work, SI increased showing a trend response to drought similar to that of SD (regardless of the heat stress), suggesting an influence on stomatal initiation. Stomatal index is primarily a function of CO 2 concentration, with different lines of evidence supporting the inverse relationship with CO 2 atmospheric levels in most modern vascular C3 plants (Woodward 1987, Woodward & Bazzaz 1988, Royer et al. 2001). On the other hand, in water-limiting environments, it has been shown that meristemoids, a key type of epidermal cell identified to be stomatal precursors (Geisler et al. 2000), are irreversible pushed into stomata differentiation when drought persists (Skirycz et al. 2011). The activity of meristemoids, often underestimated, might account for an important portion of cell division especially at the later stages of leaf development and even under permanent stress (Skirycz et al. 2011). Therefore, we postulate that at least in part, drought may exert an indirect effect on SI, through a decrease in the internal CO 2 concentration (Ci) in the leaf during stomata closure as water availability becomes scarce and that this decrease in Ci may trigger meristemoid differentiation into stomata, being such morphological response a part of a more complex acclimation strategy. Undoubtedly, this hypothesis needs further investigation to be corroborated.
Leaf thickness plays an important role in plant functioning and is related to strategies of resource acquisition and use (Jumrani et al. 2017). We observed that TLT was markedly decreased by individual stresses (HS and DS) and their interaction through decreases of size in all cell and extracellular layers (i.e. adaxial and abaxial epidermis, palisade, mesophyll, and cuticle adaxial and abaxial) ( Table I). As expected from the prominent positive association between the specific leaf weight and the TLT revealed by the biplot (Fig. 2), the responses of SLW were similar to TLT, with greater and significant decreases for the heated treatments alone or combined with DS than for controls. In greenhouses research, significant reductions of soybean leaf thickness and SLW were also observed in heated plants (Jumrani et al. 2017) et al. 2002). In well-watered crops, leaves must dissipate vast quantities of heat to avoid damage, being the sensible heat loss (through conduction and convection) to the cool air one of the most important processes in the regulation of leaf temperature (Taiz & Zeiger 2002). Canopy temperature measurements during HS treatment imposition showed that leaves from HS, DS, and HS × DS plants were 4.4, 5.3, and 6.2°C above that from control plants (23.6°C). We hypothesize that the thinning of the leaves under stressful conditions might be regarded as an acclimation mechanism. On the other hand, we also observed a strong negative association between canopy temperature during HS treatment and the TLT highlighted by the biplot (Fig. 2). Nevertheless, the response direction of TLT and SLW to heat or drought seems not always to be the same (reduction), as indicated by Djanaguiraman et al. (2011) who found that day/night temperatures of 38/28°C significantly increased the soybean TLT (21%) above that of control growing at 28/18°C. A considerably thicker mesophyll, and consequently an increased specific leaf weight were also reported for common bean when less irrigation was applied (Silva et al. 1999). Even when the responses observed by these authors (Djanaguiraman et al. 2011, Silva et al. 1999 are not obvious, none of them provided any explanation nor did they propose possible physiological mechanisms involved in leaf structural changes to heat or drought stress. Even when the structural changes involving an increase in stomatal density (Table II) and abundant intercellular spaces in the mesophyll (Fig. 3c-h) may imply an improvement for the CO 2 diffusion, the decrease in LA and TLT inexorable decreased the number of CO 2 fixation sites, detrimentally affecting photosynthesis. Interestingly, a significant and strong (p=0.0238, R 2 = 0.60) positive linear association was detected between ᶲPSII and specific leaf weight across treatments (data not shown), the latter positively correlated with TLT (Fig. 2). . Thylakoid membrane damage (Fo/Fm ratio) (a), and chlorophyll content (SPAD value) (b) of the third trifoliate leaves from the main stem apex of soybean plants during seed filling exposed to four water and temperature treatments: control (non-heatstressed and nondrought-stressed), heat-stressed, drought-stressed and heat-stressed × drought-stressed plots. Vertical bars denote ± SE of means. Black arrow = onset of heat treatment, white arrow = end of heat treatment.
At the subcellular level, ultrastructure micrographs of chloroplast of the guard (Fig.  5b-d), mesophyll (Fig. 6b-d), and bundle sheath (Fig. 7b-d) cells showed that both the single and simultaneous application of HS and DS caused adverse impact on chloroplasts membranes integrity. Inducing high disorganization of thylakoid membrane system composed of stacking grana and stroma thylakoids. Additionally, membranes of mitochondria and plasmalemma were disrupted exhibiting diffuse outlines, indicative of severe damage. The thylakoid membrane damage observed through imagines was supported by measurements at the field of the Fo/Fm ratio, which exhibited significant increases under DS alone or combined with HS (Fig. 9a). Similar ultrastructural chloroplast damage was observed in soybean and turfgrass under HS (Xu et al. 2006, Djanaguiraman et al. 2011, in Cynodon dactylon (L.) Pers. under DS (Utrillas & Alegre 1997), and Zea mays L. and Helianthus annuus L. under the combination of HS and DS (Dekov et al. 2000). Besides, abnormal chloroplasts with a swollen shape and expanded volume were also observed in the micrographs of the mesophyll (Fig. 6b-d) and bundle sheath (Fig. 7b-d) cells, also confirmed by measurements of chloroplasts dimensions (i.e. width and length) as shown in Fig. 8a-b. These changes at least in part may be related to the disruption of the chloroplast envelope as well as the distortion of the thylakoid membrane system abovementioned, and the accumulation of plastoglobules (thylakoid-associated lipid droplets). Interestingly, plastoglobules have been demonstrated to be involved in thylakoids lipid remodeling repairment during stress and in the conversion from chloroplasts to gerontoplasts, senescing plastids that appear during leaf senescence (Rottet et al. 2015). Under our experimental conditions, the involvement of plastoglobules in both processes cannot be discarded, since the exposure of soybean plants late in the reproductive period to prolonged water deficit combined with HS promoted leaf senescence as shown in our recent study (Ergo et al. 2018).
The different stress conditions induce excitation energy excess (EEE) processes, given principally by the uncoupling between the light and carbon reaction of photosynthesis. Dependent on the severity of the stress condition, the EEE could induce acclimation or deleterious responses, where the reactive oxygen species (ROS) are major modulators acting as a signal or toxic molecules depending on a delicate equilibrium between their production and scavenging. It is well established that organelles such as chloroplast, mitochondria, and peroxisomes are a major source of ROS in plant cells, being the rate of ROS production increased under abiotic stress conditions (Gill & Tuteja 2010). Moreover, it has been found that linoleic and linolenic acids, the major polyunsaturated fatty acids in the plant membrane galactolipids (thylakoid membrane) and phospholipids (all other membranes) are particularly susceptible to ROS attack (Møller et al. 2007). Therefore, the remarkable changes observed in the membranes of chloroplasts, mitochondria, and plasmalemma may partly be related to an increased accumulation of intracellular levels of ROS, as observed in our previous research where FRAP values were lowest in HS and/or DS plots, evidencing elevated levels of oxidative stress compared to control plots (Ergo et al. 2018). Since more than 80% of the chlorophyll is located in the grana thylakoids, also detected to be highly enriched with photosystem II (Albertsson 2001), it may be inevitable the interplay between the ultrastructural chloroplast injury of stressed leaves, the chlorophyll loss and the damage to the photosynthetic machinery. These observations explained the strong association between leaf morphology and function. Indeed, chlorophyll fluorescence variables, which are mainly related to PSII function, showed a significant reduction of ᶲPSII and Fv/Fm in response to HS and DS, indicative of a decline in the photosynthetic performance and integrity, respectively (Ergo et al. 2018). In addition, the chlorophyll content decreased as shown by drops in SPAD value (Fig.  9b). This significant loss of chlorophyll pointed out to senescence processes accelerated by DS alone or combined with HS, in a crop characterized by its high seed protein content, where redistribution of N is an important source for protein biosynthesis (Masclaux et al. 2001).
Overall, our results indicate that leaf area is mainly affected by DS, leading to a reduction of water loss to the atmosphere (i.e. water economy); meanwhile, leaf thickness is mainly affected by HS. However, these acclimation strategies conduct to an irreversible reduction in CO 2 assimilation sites. On the other hand, at the subcellular level, both stresses caused the accumulation of plastoglobules and damaged membrane integrity of plasmalemma, mitochondria, chloroplasts, as well as thylakoids. These effects triggered the loss of chlorophyll content and negatively affected chlorophyll fluorescence. The decrease of the CO 2 fixation sites together with the photosynthetic machinery impairment resulted in less assimilate production detrimentally impacting crop productivity. Besides, the differential effects of each type of stress could be because in our field experimental conditions drought developed gradually and heat rapidly. Although both stresses were applied several days they presented contrasting development patterns, which are common in field-grown crops since air temperature changes are likely to occur more rapidly than soil moisture variations. Interestingly and in contrast with what has been expected, not only the gradualism in the development of water deficit but also the rapidity of the heat stress progress were translated into changes in the structure and ultrastructure of leaf components in favor of a new homeostasis state, allowing plants to acclimate.

CONCLUSIONS
In summary, DS and HS induced structural and ultrastructural leaf changes in close relation with physiological-biochemical variables associated with yield and its components in field-grown soybean. In addition to these strong associations between physiological-biochemical variables and yield, the present multiscale approach disclosed other variables, namely, alterations in morpho-functional characters, often overlooked in field studies, contributing to understanding the effects of combined stresses on soybean yield reduction.