Selecting Drought-Tolerance Markers: An Exploratory Analysis in Contrasting Soybeans

Laila Toum Instituto de Tecnología Agroindustrial del Noroeste Argentino, Estación Experimental Agroindustrial Obispo ColombresConsejo Nacional de Investigaciones Cientí cas y Técnicas (CONICET) Lucía Sandra Pérez-Borroto Instituto de Tecnología Agroindustrial del Noroeste Argentino, Estación Experimental Agroindustrial Obispo ColombresConsejo Nacional de Investigaciones Cientí cas y Técnicas (CONICET) Andrea Natalia Peña-Malavera Instituto de Tecnología Agroindustrial del Noroeste Argentino, Estación Experimental Agroindustrial Obispo ColombresConsejo Nacional de Investigaciones Cientí cas y Técnicas (CONICET) Catalina Luque Cátedra de Anatomía Vegetal. Facultad de Ciencias Naturales e IML, Universidad Nacional de Tucumán, Miguel Lillo 205, Tucumán, Argentina. Bjorn Welin Instituto de Tecnología Agroindustrial del Noroeste Argentino, Estación Experimental Agroindustrial Obispo ColombresConsejo Nacional de Investigaciones Cientí cas y Técnicas (CONICET) Ariel Berenstein Instituto de Química Biológica (IQUIBICEN) de la Facultad de Ciencias Exactas y Naturales (FCEyN) de la Universidad de Buenos Aires, Intendente Guiraldes 2160, Buenos Aires, Argentina. Darío Fernandez Do Porto Instituto de Química Biológica (IQUIBICEN) de la Facultad de Ciencias Exactas y Naturales (FCEyN) de la Universidad de Buenos Aires, Intendente Guiraldes 2160, Buenos Aires, Argentina. Adrian Vojnov Instituto de Ciencia y Tecnología “Dr. César Milstein”, Fundación Pablo Cassará-Consejo Nacional de Investigaciones Cientí cas y Técnicas (CONICET) Atilio Pedro Castagnaro Instituto de Tecnología Agroindustrial del Noroeste Argentino, Estación Experimental Agroindustrial Obispo ColombresConsejo Nacional de Investigaciones Cientí cas y Técnicas (CONICET) Esteban Mariano Pardo (  mpardokarate@gmail.com ) Instituto de Tecnología Agroindustrial del Noroeste Argentino, Estación Experimental Agroindustrial Obispo ColombresConsejo Nacional de Investigaciones Cientí cas y Técnicas (CONICET)


Introduction
Soybean [Glycine max (L.) Merrill] represents one of the most important sources of vegetable oil and protein in the world 1 . Calculation models based on the growing global population and current agricultural production suggest that crop yields, including soybean, must be doubled to provide enough food in 2050 2 . Yield, the principal breeding target for most crop plants, is massively affected by suboptimal growth conditions primarily due to climate factors such as drought and extreme temperatures. In addition, the progressive climate change will reduce water availability for many rainfed crops like soybean, affecting their growth and productivity 3 . Hence, breeding for yield stabilization and drought tolerance in soybean and other crops is essential for sustainable agriculture and food production 4 .
The development of cultivars with improved yield under water de cit has had relatively limited success for several reasons. First, the direct selection for yield improvement under drought is expensive, time-consuming, laborious and complex due to intrinsic genotype by environmental interactions 5 . Moreover, inherent stress susceptibility is often masked by spillover effects of high yield potential when determining plant performance under drought conditions. Consequently, a high-yield variety will often give signi cant yields during drought, even though the relative yield reduction can be substantial 6 . Instead, analytical approaches that emphasise breeding for yield stabilization through an indirect selection strategy, using morphophysiological or biochemical markers, have gained increasing attention 7,8 . However, the challenge for effectively using targeted breeding approaches lies in developing reliable and reproducible markers. These markers should be (i) strongly related with yield and stress-tolerance traits, where possible, (ii) be non-destructive, (iii) be easily measurable in early phenological stages and (iv) have a high narrow-sense heritability to facilitate the selection in breeding populations 9 . Therefore, identifying and validating drought-tolerance traits are essential steps to obtaining valuable markers for breeding programs and selecting superior genotypes. Soybean drought tolerance has been evaluated through markers such as water use e ciency (WUE), root morphology and penetrability of hardpan, leaf wilting, excised leaf water loss and relative water content (RWC) with varying degrees of success 10 .
Advances in next-generation sequencing (NGS) and subsequent evolution in multi-omics technologies have contributed to understanding some underlying mechanisms of response to water de cit in soybean and other crops 11 . Omics studies, however, must be complemented with morphophysiological and biochemical approaches to ensure an integrative perspective of plant adaptation to water scarcity and accurately assess the role of individual traits regarding stress tolerance and yield. Usually, the main issue in such studies is the lack of a well-de ned and reliable phenotyping methodology that validate the traits accuracy 12 . It is safe to say that, currently, phenotyping systems are the major operational bottleneck in plants breeding, limiting the translation of genetic and genomic analysis into stress-tolerant phenotypes.
In a previous evaluation of drought tolerance under greenhouse and eld conditions, two contrasting soybean cultivars were identi ed from our germplasm: MUNASQA (drought-tolerant) and TJ2049 (drought-susceptible) 13 . In this study, in-depth phenotyping was conducted for both genotypes to identify and statistically validate markers associated with drought-tolerance and yield-stabilization traits, aiming for their future use in soybean breeding programs, including ours, to select genotypes better adapted to low water availability.

Results And Discussion
Here, we show differences in the molecular, morphophysiological and biochemical responses of two contrasting soybean genotypes, MUNASQA (drought tolerant) and TJ2049 (drought susceptible), subjected to mild water de cit treatments in V3 and R5 phenological stages.

Molecular insights
Transcriptional changes were assessed in MUNASQA and TJ2049 after 72 h of exposure to mild water de cit. From 38.658 transcripts analysed, droughtstressed MUNASQA and TJ2049 plants exhibited 2952 and 1126 transcripts with signi cant changes (P < 0.05) in their expression levels (Suppl. Data 1), respectively. After an FDR=0.1, 399 and 15 transcripts were assigned as MUNASQA and TJ2049 DEGs (Suppl. Data 2). The transcript loss detected in TJ2049 might be explained by the larger variation among replicates observed in the water de cit samples and measured as SD/mean ratio for each gene (Suppl. Table  1, Fig. S2). However, in a previous exploratory 454 sequencing experiment, TJ2049 showed signi cantly fewer DEGs than MUNASQA (data not shown).
Large-scale transcriptional reprogramming has long been recognised as the rst response to drought, initiating stress mitigation pathways and metabolic changes 14 . Moreover, the quickness to sense and respond to stresses could be an essential factor to differentiate tolerant and susceptible genotypes. Here, after a mild drought, far from normal eld levels, the genotypes exhibited dramatic changes on their transcriptional pro les and no observable phenotypic alteration. In the heat map (Fig.1A), DEGs expression pro les were classi ed in 8 clusters (Suppl. Data 3), and almost 50% of MUNASQA DEGs showed repression under stress conditions, a difference reinforced by the Venn diagram (Fig.1C).
Generally, in response to drought, plants initially trigger transcriptional control and hormone signalling, leading to metabolic adjustment for coping with low water availability. Cellular mechanisms such as water/ion uptake and transport, redox homeostasis, ROS scavenging, osmoregulation and membrane protection are accompanied by physiological responses such as stomata regulation, root development and protection of the photosynthesis machinery 15 .
The expression of ten randomly selected DEGs by qRT-PCR (Fig. 1D) agreed with the MACE pro le data for MUNASQA and TJ2049, showing a high similarity between both methodologies. Moreover, genes encoding detoxifying proteins like SOD, CAT and APX were differentially regulated in both genotypes (Suppl. Table 3). Enzymes like these are stress-response indicators, frequently regulated at both transcriptional and post-transcriptional levels under drought 16 . Here, we detected two up-regulated SOD genes under stress conditions in MUNASQA and, contrary to TJ2049, three down-regulated genes for APX and CAT. In accordance, SOD activity measurements in R5 corroborated these differences in transcriptional activity between the two genotypes (Suppl. Table 4).

Morphophysiological and biochemical phenotyping
MACE assay, by itself, does not fully explain the extent of MUNASQA and TJ2049 responses to the applied stress. Thus, extensive morphophysiological and biochemical measurements were performed to minimise the gap between the transcriptional regulation and phenotypical alterations observed. Here, 22 markers associated with stress-response, growth and water-use were evaluated in both genotypes exposed to mild water de cit in V3 and R5 stages.

Stress-response markers
An e cient antioxidant activity is crucial for withstanding low cellular water content, and its importance in plant drought tolerance has been extensively reported 17 . Here, substantial differences were found in MUNASQA and TJ2049 SOD, APX, POX, and CAT enzymes regulation over time, under water de cit and independently of the phenological stage (Table 2). MUNASQA reached maximum activity for all enzymes after 4 d of water de cit, while the highest activity for TJ2049 occurred to 8 d after the stress onset, except for CAT. Noticeably, under non-stressed conditions, all enzymatic activities, excluding CAT, showed higher activity in the tolerant genotype than in the susceptible one, strengthening the hypothesis that TJ2049, compared to MUNASQA, presents delayed stress perception and response mechanisms. In fact, and according to Carvalho et al. (2019) 16 , a successful response to water de cit may depend not only on which enzymes are activated but also on the activation timing. Moreover, the general gene expression pro le regarding these enzymes regulation was highly consistent with the biochemical results.
Regarding PRO, one of the most common osmoprotectants in plants 18 , MUNASQA and TJ2049 accumulated the osmolyte in response to water de cit over time and in both V3 and R5 stages ( Table 2). Several studies have demonstrated a direct correlation between high osmoprotectant accumulation and drought tolerance in many crops 19 . Here, the tolerant genotype MUNASQA exhibited a higher and more rapid accumulation of PRO after 4 d of water de cit at both phenological stages. In agreement with our results, a recent study in soybean reported higher PRO accumulation in the drought-tolerant genotype A5009 RG, compared with the drought-susceptible ADM50048 20 .
When analysing MDA production, an indicator of lipid peroxidation and stress severity 21 , a higher and more rapid accumulation was detected in TJ2049 plants in response to water de cit, compared to MUNASQA (Table 2).
Drought also affects leaf pigments content 21 . Changes in photosynthetic pigments can alter various light-harvesting processes, while the accumulation of photoprotective compounds plays an essential role in preventing photo-oxidative damage 22 . Here, MUNASQA and TJ2049 showed alterations in pigment content under water de cit ( Table 2). The CHL was signi cantly reduced over time due to stress in both genotypes, but this reduction was signi cantly lower in the tolerant one. Similar results were found in drought-tolerant maise that showed lower CHL reductions under stress than susceptible genotypes 23 . Regarding CAR levels, TJ2049 and MUNASQA showed an increased synthesis in response to drought, although the last one exhibited greater and faster accumulation.
According to our results, MUNASQA is more e cient at removing reactive oxygen species (ROS) under drought stress conditions, an essential mechanism for stress tolerance 24 .

Growth and yield markers
The ability to produce high seed yield or biomass under limited water access is considered the optimal indicator of drought tolerance in crops [25][26][27] . We evaluated the effect of mild water de cit in MUNASQA and TJ2049 growth and yield by monitoring changes in various markers related to biomass and seed production, including LAI, LAR and NAR, RGR, CGR, relative yield and yield-DSI (Table 3).
Total leaf surface area and LAI are strongly associated with canopy interception, evapotranspiration and photosynthesis 28 . Here, independently of the phenological stage or water availability, MUNASQA plants exhibited a higher LAI compared to TJ2049, indicating a larger assimilatory capacity and, as a consequence, photosynthetic potential. In vegetative stages, a higher LAI denotes a more rapid canopy development, favouring greater and faster soil coverage and thus less water loss from direct evaporation. In general, drought suppresses leaf initiation and growth and consequently affects LAI 29 . Therefore, a decrease in LAI is expected in plants exposed to water scarcity. As expected, MUNASQA and TJ2049 plants showed a reduction in LAI in response to drought, more noticeable after 8 d of stress and in the R5 stage.
Drought also alters the LAR: the leaf area development in relation to the total biomass produced 30 . Here, water de cit only affected LAR during the vegetative stage, displaying a signi cant reduction in both genotypes exposed to stress. The highest ratio between plant leaf area and biomass is reached at the beginning of the plant life cycle 31 , which explains the highest LAR values at the rst sampling day during the vegetative stage.
In MUNASQA and TJ2049, the relationship between leaf area expansion and the biomass produced over time (NAR) was also reduced due to water de cit at both phenological stages. Changes in photosynthetic e ciency were more signi cant at V3, where both genotypes exhibited greater NAR values. Noticeably, well-watered MUNASQA plants showed a signi cantly higher NAR than TJ2049 ones, while, in response to stress, MUNASQA's NAR increased in contrast to TJ2049 values, which were reduced. Although not as accentuated, a similar pattern was observed in plants exposed to water de cit in R5. Here, LAI, LAR and NAR results indicated that, in response to drought, MUNASQA plants regulated photosynthates allocation to leaves and the maintenance of photosynthetic e ciency more e ciently than TJ2049 ones.
In addition, the relative growth rate (RGR) or biomass produced over time was signi cantly reduced by water de cit in the vegetative stage (Table 3).
Interestingly, the reduction in RGR was more pronounced in TJ2049 plants. Regarding CGR, signi cant differences were observed in the V3 stage after 8 d of drought, while in R5 were detected after 4 d.
These results agree with the differences in yield and yield-DSI exhibited by MUNASQA and TJ2049 after water-de cit treatments in V3 and R5 (Fig. 2).
According to these ndings, when drought was applied in the V3 stage, TJ2049 showed a distinct but not signi cant yield penalty and a yield-DSI considerably higher than MUNASQA. Moreover, after a mild water de cit in R5, a highly moisture-sensitive phenological stage, TJ2049 exhibited the largest yield loss and a signi cantly higher yield-DSI than MUNASQA. Results that also agree with the ones reported by Pardo et al. (2015) 13 .

Water use markers
Maintaining tissue/cellular water content and/or metabolic activity at low water potentials are physiological strategies to survive drought 32 . Traits like pubescence, changes in stomatal density, slow wilting, leaf thickening, or tight control of stomatal closure, canopy temperature, RWC, and WUE are essential to determining drought tolerance in plants. Thus, all these water-use parameters were assessed in MUNASQA and TJ2049 plants in response to drought.
Morphological drought-adaptation traits, such as stomatal, trichrome density and leaf thickness 33 , showed signi cant differences after prolonged exposure to mild water de cit (21 d) ( Table 4). The stomatal and trichrome densities were substantially altered in response to drought in both genotypes. Interestingly, in MUNASQA, stomatal density increased in response to stress on both leaf sides; in contrast, TJ2049 plants exhibited a signi cant decrease. Regarding pubescence, no signi cant alterations were observed on the adaxial surface in response to drought; however, on the abaxial side, an increase of trichrome density was quanti ed in MUNASQA plants, which was higher than in TJ2049. Generally, plants that grow under water de cit develop denser stomata and trichomes 34 , mainly on the abaxial surface 35 , which helps prevent excessive water loss through transpiration. Our results suggest that MUNASQA morphoanatomic adaptations favour the genotype ability to cope with water de cit.
When evaluating the genotypes leaf thickness, no signi cant differences were detected under well-irrigated conditions. However, and unexpectedly, after 21 d of water de cit, MUNASQA leaves were signi cantly thinner, while TJ2049 ones were thicker. The knowledge about the links between leaf morpho-anatomy and its function under non-stressed/stressed conditions is relatively poor, especially for traits like leaf thickness 36 that is strongly related to transpiration 37 and reported by some authors as a drought-tolerance trait that maintains turgor pressure and enhances photosynthesis 38 . Yet, this feature, only apparent in TJ2049 plants under water de cit, could not be associated with adjustments in transpiration, nor with an increase in photosynthesis, or any other droughttolerant feature.
On the other hand, the regulation of stomatal aperture reinforces the drought-tolerant character of MUNASQA. Under non-stressed conditions, the susceptible TJ2049 showed more opened stomata ( 22% more than MUNASQA). Moreover, after 3 d of water de cit in V3, this difference was increased to almost 50% of stomatal aperture (Table 4). Although stomata represent a small percentage of the leaf lamina, large amounts of water evaporate through them 33 . Thus, the lack of stomatal control in TJ2049 might explain the fast-wilting phenotype displayed in response to air desiccation (Fig. 3). Measuring the water loss of detached leaves is a selection method for drought tolerance 39 . Here, R5 leaves of both genotypes were removed and air-dried. After 48 h, MUNASQA exhibited a greener, healthier and slow-wilting phenotype, previously linked to drought tolerance in studies with soybean cultivars 17 .
MUNASQA water-saving behaviour was also con rmed through parameters such as RWC and WUE, strongly regulated under drought ( Table 5). The RWC, a time-speci c measurement of the hydric status of a plant, is considered a physiological character recommended for drought-tolerance selection 40 . In V3 and R5 stages, well-irrigated MUNASQA plants showed higher RWC than TJ2049, even with a smaller canopy. As expected, under water de cit, the RWC values were reduced by 14% in MUNASQA and 20% in TJ2049, indirectly con rming the effectiveness of the stress imposed.
Water use e ciency (WUE), referring to the biomass produced per water unit, has been widely used as a breeding target in many rainfed crops, including soybean. Conservative water-use strategies are associated with high leaf WUE 34 . In agreement, and contrarily to TJ2049, V3 and R5 MUNASQA plants showed a gradual increase of WUE in response to water de cit that agrees with the tighter regulation of stomatal movements and the reduced water loss observed in the genotype (Table 5). Moreover, considering the discrete NAR reduction and the maintenance of a ~70% RWC under water de cit, we hypothesise that MUNASQA may display a stomatal control based on partial or total/partial closure intervals, therefore reducing transpiration and saving water through a smaller gas exchange (potential photosynthesis) penalty.
As a surrogate trait for stomatal conductance, the CTD, regarding the canopy temperature difference with the surrounding air, is a good indication of plant transpiration rate 41 . As expected, in response to drought, MUNASQA plants evidenced lower CTD values, a nding that supports the stomatal aperture results and strongly suggest the genotype water-saving behaviour. Plants with higher stomatal conductance transpire more and thus maintain a cooler canopy 42 . Thus, in TJ2049 stressed-plants, the high and positive CTD con rmed a higher stomatal aperture and transpiration rate that agrees with a water-spender behaviour. Moreover, TJ2049 also presented higher transpiration rates in unstressed conditions. Finding that could be evidence of a natural and predisposing difference between tolerant and susceptible genotypes.

Markers selection
Identifying and exploiting phenotyping markers will improve selection strategies for drought tolerance in legumes crops 17 . However, to successfully implement markers in a breeding program, it is imperative to validate their (i) accuracy, (ii) feasibility and (iii) strength of association with the desired trait. To further understand the markers contribution to drought tolerance and yield stabilization in MUNASQA and TJ2049, a Principal Component Analysis (PCA) was performed for each set of the parameters, previously grouped in (i) "stress response", (ii) "growth" and iii) "water use" categories. In addition, markers evaluated at V3 and R5 stages were measured by Pearson's correlation to determine their strength of association between the phenological stages (accuracy) and yield stabilization (desired trait) (Tables 2, 3 and 5).
Biochemical parameters such as enzymatic and non-enzymatic ROS scavengers, leaf pigments and MDA have been con rmed as adaptive responses to desiccation stress, frequently used for selecting plant genotypes under drought 43 . Clear discrimination of MUNASQA and TJ2049 drought responses were observed in PCA results (Fig. 4A). Here, the rst two principal components (PC) explained 96.6% of total variation (PC1 = 75.0% and PC2 = 21.6%). Data were clustered by irrigation treatment in PC1, suggesting that these markers are indicators of phenotypic plasticity. The CHL was associated with well-irrigated plants, while all enzymes, MDA, CAR and PRO, were related to drought stress (Fig.4A).
All these markers showed high accuracy between phenological stages due to signi cant (P < 0.001), strong (r 2 over 0.88) and positive correlations (Table 2). However, the correlation with yield showed inconsistent outcomes. PRO, MDA and CAR were signi cant and positive correlated with yield. Meanwhile, except SOD, all enzymes exhibited signi cant and positive correlations that were too variable in the strength of association and could not be linked with a speci c genotype or treatment. Thus, we considered these "stress-response" markers suitable for discriminating susceptible/resistant responses during early droughttolerance screenings in soybean. Still, due to the high environmental effect, their use as breeding traits is limited. Interestingly, a report 20 found that PRO and CHL were suitable markers for ranking soybean genotypes in response to drought in vegetative stages (5 days after emergence), while MDA could be useful during R5 as a sensitivity trait.
In legumes, features like LAI, smaller leaf area, leaf area maintenance and dry matter partitioning have been used to screen for drought tolerance 17 . Here, in the "growth" PCA ( Fig. 4B), the rst two PC explained 93.0% of total variation (PC1 = 73.0% and PC2 = 20.0%). In PC1, data were clustered by genotype, suggesting that these markers explained differences between MUNASQA and TJ2049 growth responses on their intrinsic genetic variability. Moreover, the markers analysed were only associated with MUNASQA, LAR and RGR associated with stressed plants. However, only LAI and LAR showed a weak correlation between phenological stage and yield (Table 3).
Water-saving features like denser leaf pubescence, a higher number of stomata, warmer canopies, RWC maintenance and increased WUE have been associated with drought tolerance in legumes and applied in drought-resistance breeding 17 . Here, "water-use" markers contributed the most to discriminating drought-tolerant and susceptible responses. In the two PCA performed, one for physiological markers and the other for morphological ones, data were clustered by genotype in PC1; thus, these markers are indicators of genetic variability.
In the PCA made with morphological "water-use" markers, the rst two PC explained 96.6% of total variation (PC1 = 73.8% and PC2 = 22.8%) (Fig. 4D). Here, pubescence and stomata abaxial density were strongly related to MUNASQA. Although the data were insu cient to execute good correlation analysis, these morphological markers have been demonstrated as clear indicators of water-saving strategies in legumes 34 .
A good drought-resistance marker linked to yield stabilization must also be accurate, cost-effective, if possible non-destructive and easily measurable. Hence, after evaluating marker accuracy (amid phenological stages) and assessing which ones better explained the phenotypic variability between genotypes and treatments, a nal selection was performed by cost-feasibility (CF) and statistical weight (SW) rankings. Based on the CF values, the degree of complexity and cost to assess each indicator, 9 markers were further selected encompassing the categories 1 (easy and cheap) and 2 (easy and expensive) ( Table 6).
Subsequently, after SW re-selection, 4 markers remained with "High" SW in both PC. The selected markers were (i) stomatal density on the adaxial and (ii) abaxial leaf surface, (iii) trichrome density on the abaxial side and (iv) CTD. These four traits were categorised as the most e cient phenotyping markers considering their high accuracy, strong association to water-saving strategies, feasible and non-destructive measurement and amenability to high throughput screening using high-resolution imaging.
Taken together these results: (i) con rmed the effectiveness of the phenotyping methodology and (ii) its successful application in early phenological stages, (iii) the usefulness of MUNASQA and TJ2049 as model genotypes for screening potential drought-tolerance markers and, (iv) suggested that pubescence, stomatal density and CTD are informative markers to discriminate soybean adaptations to withstand drought via dehydration avoidance.
Traditional crops phenotyping is labour-intensive, time-consuming and often destructive. Hence, a well-characterized drought-tolerant phenotype could be useful for developing new commercial genotypes or identifying drought-tolerance markers. In this context, our ndings con rm that MUNASQA and TJ2049 constitute suitable models for drought-tolerance screening in soybean. The deep-phenotyping results allowed us to ensure the methodology effectiveness and its possible implementation in early vegetative stages like V3, reducing time and money costs. Finally, identifying stomatal densities, abaxial pubescence and CTD as pro table and non-destructive markers, strongly associated with water-saving adaptation strategies, provides us with promising tools for early selection and prediction of drought tolerance in soybean.

Experimental approach
The response of MUNASQA and TJ2049 to mild water de cit applied in the R5 stage was assessed through transcriptional and leaf morphology analysis.
Subsequently, comparative studies were performed to determine the genotypes response to water de cit imposed in V3 and R5 stages. Next, all markers assessed were analysed according to their strength of association between phenological stages and yield, then were ranked by statistical weight and costfeasibility (Fig. S1). (2015) 13 , the water de cit was applied by maintaining the pots at 14% of VWC (Ψs =-0.65 MPa) for ten days. The desired Ψs was reached in 2-3 days. At the end of stress, plants were fully watered until harvest. The Ψs was daily monitored and recorded. Corrections for soil water status were made by weighing two plants per genotype and treatment every 3 days. The plant water status was monitored through the RWC 47 to ensure stress occurrence.

Plant material and growth conditions
All drought experiments were carried out for three consecutive years, always applying the previously described irrigation treatments. Water de cit in V3 and R5 phenological stages was imposed in independent plant sets; the sections below detailed the sampling process.

MACE-transcriptomic analysis and validation
Three biological replicates per treatment (Control and Stress) and genotype (MUNASQA and TJ2049) were collected from R5 plants after 72 h of water de cit (n=12), and RNA from fully mature expanded leaf between nodes 5 to 7 node was isolated for transcriptional analysis.
MACE-Seq libraries and sequencing were performed on an Illumina NextSeq500 machine (1x75 bp reads). The conversion was made with bcl2fastq2 software (version 2.19.1), and the cleaning of duplicate sequences was performed with "TrueQuant". In MACE-seq, the TrueQuant barcodes each DNA molecule before PCR ampli cation. As each barcode-template combination is statistically unique, PCR-duplicates can be identi ed and eliminated from the dataset to prevent PCR bias. Bases with low sequencing quality were clipped. Next, reads were mapped into genome version "Gmax_275_v2.0.fa" of soybean downloaded with standard parameters from Phytozome.net and Bowtie2 (version 2.2.4). Then, expression analysis was performed by in-house scripts and DESeq2 (R-package).
DEGs were de ned at an FDR of 0.1 and listed as either up or down-regulated. A heat map plot was generated using R software (version 3.4.1). Then, hierarchical clustering was applied by considering a cut-off threshold of 8 expression pro les (clusters). Venn diagrams were depicted using the VennDiagram R package.
A GO enrichment analysis was performed using the topGO R package 48 , while the GO annotation le was extracted from agriGO website 49 . Each GO term, containing at least two DEGs, was analysed by Fisher's exact test. The resulting P values were corrected by FDR multiple testing approach. GO terms with an FDR lower than 0.1 were considered for further analysis.
Ten DEGs were randomly selected and measured by qRT-PCR assays (Applied Biosystems) to validate MACE results. F-BOX gene was used as an internal reference to standarise the expression of target genes, and the ratio between treatments was calculated according to 50 . All primers used are listed in Suppl. Table 6. Data analysis and primer e ciencies were obtained using LinReg PCR software 51 . Relative expression ratios and statistical analysis were performed using fgStatistics software interface 52 . The cut-off for statistically signi cant differences was set as P < 0.05, indicated as *.

Antioxidant measurements
Additionally, antioxidant proteins encoded by DEGs detected in MACE were analysed. Five biological samples were collected per treatment and genotype Wilting air desiccation assay Response to air desiccation was evaluated in the R5 stage. Three whole leaves per genotype (n=6) were collected and exposed to air desiccation at 32 o C. After 0, 6, 24, 36 and 48 h of air exposure, plants were photographed with a Canon Power Shot SX520 HS (14 Mpx), and the wilting rate was assessed.

Comparative analysis of genotypes responses to water de cit in V3 and R5
The responses of MUNASQA and TJ2049 to water de cit, applied in V3 and R5 stages, were compared by measuring morphophysiological and biochemical parameters grouped by biological processes (BP) ( Table 1). Four treatments were de ned: (i) Control-V3, (ii) Control-R5, (iii) Stress-V3 and (iv) Stress-R5, and three sampling times were performed (0, 4 and 8 d of water de cit). Ten plants per genotype, treatment and time were collected and used for markers evaluation (n=240).
Additionally, for treatments (i), (iii) and (iv), 50 plants per genotype were harvested at physiological maturity (n=300) to quantify relative yield and calculate the relative yield DSI (Drought Susceptibility Index) according to Fischer and Maurer (1978) 61 .
The markers evaluated and their methodologies are detailed below.

Markers
As stress response markers, the activities of SOD, APX, POX and CAT proteins were assessed, together with the accumulation of free proline (PRO) 62 , malondialdehyde (MDA) 63 , total chlorophylls (CHL) 64 , and carotenoids (CAR) 65 .
As growth indicators, the plant total leaf area (TLA) and biomass (plant total dry weight) were determined. Then leaf area index (LAI) and leaf area ratio (LAR) 66 , the net assimilation rate (NAR) 67 , relative growth rate (RGR) 68 and crop growth rate (CGR) 69 were calculated.
Finally, as water-use parameters, plant RWC and WUE 70 were calculated. Moreover, the canopy temperature was monitored and recorded using a FLIR ONE-3 thermal camera (0.3456 Mpx) to calculate canopy temperature depression (CTD) 71 .

Univariate analysis
Data from stomatal apertures were subjected to a two-way ANOVA (Factor 1: genotype, Factor 2: treatment). The remaining data were analysed through ANOVA with post hoc contrast by Tukey's HSD test. Data were analysed with InfoStat statistical package 52 and presented as the arithmetic mean ±SE. Means were considered signi cantly different at P < 0.05.
All markers, grouped in (i) "stress response", (ii) "growth", and iii) "water use" sets, were subjected to PCA to determine which markers best explained the phenotypic variability between genotypes and treatments.
Additionally, the markers were ranked by CF and SW. The markers CF, in terms of their complexity and evaluation cost, was assigned according to 4 categories: easy and cheap (1), easy and expensive (2), complicated and cheap (3) or complicated and expensive (4). Meanwhile, the SW was obtained from PCA variables coe cients (autovectors e1 and e2) that were ranked and classi ed in "Low" (Low = [-2, 2]) and High (High = ℝ -[-2, 2]), according to their weigh on PC1 and PC2. Markers strongly correlated between phenological stages, if possible with yield, together with CF values of 1 or 2 and "High" SW in both PC, were selected as the most e cient phenotyping markers.    Tables   Table 1 Markers evaluated in MUNASQA and TJ2049, clustered by biological processes (BP).  Table 2 Effect of mild water de cit on stress-response enzymatic markers measured in MUNASQA and TJ2049. SOD, APX, POX and CAT activities, together with PRO, MDA, CHL and CAR contents, were obtained from plants submitted to water de cit (Ψs =-0.65 MPa) and well-watered treatments (Ψs =-0.05 MPa) applied in V3 and R5 stages. Two independent experiments (n= 10 per genotype/treatment) were conducted, assessing parameters at 0, 4, and 8 d after stress (DAS) imposition. Additionally, 50 plants per genotype and the following treatments: 1: Control, 2: V3-Stress and 3: R5-Stress, were harvested at physiological maturity to obtain relative yield. Average values followed by the same uppercase letter in the column and the same lowercase letter in the row do not differ statistically among them within each phenological stage, according to Tukey's HSD test at 5%. The strength of association between markers evaluated in V3 and R5 stages (n=240) and between markers and yield (n=300) was measured by Pearson's correlation analysis adjusted by Bonferroni (P > 0.05 indicated as ns; P < 0.05 indicated as *; P < 0.01 ** and P < 0.001 ***). Correlation coe cients (r 2 ) were classi ed as "S: Strong" (> ±0.60) and "W: weak" (below ±0.59).  Table 3 Effect of mild water de cit on growth markers measured in MUNASQA and TJ2049. LAI, LAR, NAR, RGR and CGR were assessed in plants submitted to water de cit (Ψs=-0.65 MPa) and well-watered treatments (Ψs= -0.05 MPa) in V3 and R5 stages. Two independent experiments (n= 10 per genotype/treatment) were conducted, assessing parameters at 0, 4, and 8 d after stress (DAS) imposition. Additionally, 50 plants per genotype and the following treatments: 1: Control, 2: V3-Stress and 3: R5-Stress, were harvested at physiological maturity to obtain relative yield. Average values followed by the same uppercase letter in the column and the same lowercase letter in the row do not differ statistically among them within each phenological stage, according to Tukey's HSD test at 5%. The strength of association between markers evaluated in V3 and R5 stages (n=240) and between markers and yield (n=300) was measured by Pearson's correlation analysis adjusted by Bonferroni (P > 0.05 indicated as ns; P < 0.05 indicated as *; P < 0.01 ** and P < 0.001 ***). Correlation coe cients (r 2 ) were classi ed as "S: Strong" (> ±0.60) and "W: weak" (below ±0.59).  Table 5 Effect of mild water de cit on water-use physiological markers measured in MUNASQA and TJ2049. RWC, WUE and CTD were assessed in plants submitted to water de cit (Ψs=-0.65 MPa) and well-watered treatments (Ψs=-0.05 MPa) in V3 and R5 stages. Two independent experiments (n= 10 per genotype/treatment) were conducted, assessing parameters at 0, 4, and 8 d after stress (DAS) imposition. Additionally, 50 plants per genotype and the following treatments: 1: Control, 2: V3-Stress and 3: R5-Stress, were harvested at physiological maturity to obtain relative yield. Average values followed by the same uppercase letter in the column and the same lowercase letter in the row do not differ statistically among them within each phenological stage, according to Tukey's HSD test at 5%. The strength of association between markers evaluated in V3 and R5 stages (n=240) and between markers and yield (n=300) was measured by Pearson's correlation analysis adjusted by Bonferroni (P > 0.05 indicated as ns; P < 0.05 indicated as *; P < 0.01 ** and P < 0.001 ***). Correlation coe cients (r 2 ) were classi ed as "S: Strong" (> ±0.60) and "W: weak" (below ±0.59).

RWC (%) Correlations
Genotype  PCA for determining markers interaction with genotypes and treatments. Markers Set I: stress-response. PCA was performed using the SOD, APX, POX, CAT, MDA, PRO, CHL and CAR data (a). Markers Set II: growth. PCA was performed using LAI, LAR, NAR, RGR and CGR data (b). Physiological markers of Set III: water use. PCA was performed using the RWC, WUE and CTD data (c). Morphological markers of Set III: water use. PCA was performed using the LT, TD_AB, TD_AD, ST_AB and SD_AD data (d). For (a to c), MUNASQA and TJ2049 were exposed to water de cit (Ψs=-0.65 MPa) and well-watered (Ψs=-0.05 MPa) treatments applied in V3 and R5. Markers were determined at 0, 4 and 8 d after stress treatment was initiated. For (d), both genotypes were exposed to the same water regimen in R5. Parameters were determined 21 d after stress imposition.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.