*Article* **A Simple and Accurate Method Based on a Water-Consumption Model for Phenotyping Soybean Genotypes under Hydric Deficit Conditions**

**Sebastián Simondi 1,†, Esteban Casaretto 2,†, Gastón Quero <sup>2</sup> , Sergio Ceretta <sup>3</sup> , Victoria Bonnecarrère <sup>4</sup> and Omar Borsani 2,\***


**Abstract:** Drought limits crop productivity and reduces yield stability. Drought tolerance as a selection criterion in breeding programs requires the development of high-throughput, precise, and low-cost phenotyping strategies. We developed a mathematical model, based on biological approaches, for evaluating soybean plants' response to drought under controlled growth conditions. The model describes the kinetics of water consumption of a plant pot substrate system (PPS) with low sampling requirements. The model generated two parameters, *t*0.5 (time necessary for the PPS to reach half of the maximum amount of evapotranspirable water) and *Gw*(*t*0.5) (stomatal conductance [*Gw*] at *t*0.5), which determined the water- consumption curve of each genotype. An analysis of the kinetics of water consumption in response to a progressive water deficit in a biparental and breeding population was performed as a preliminary test of the model. A correspondence analysis between the *t*0.5 and *Gw*(*t*0.5) parameters with the genetic structure of the populations shows a genetic association. The phenotyping methodology presented in this work and drought susceptibility in field conditions are discussed based on previous results. This work could be useful for improving the selection of soybean genotypes in relation to their performance under drought conditions.

**Keywords:** drought; stomatal conductance; mathematical modeling; crop breeding

## **1. Introduction**

Despite increases in soybean (*Glycine max* L.) crop yields achieved by breeding and better agricultural practices, its productivity and yield stability are especially susceptible to drought events [1,2]. Drought stress affects both the vegetative and reproductive stages in soybean, reducing leaf area, increasing flower and pod abortion, and diminishing pod and seed size [3]. Several works relate the response in the vegetative stage to drought tolerance during reproductive stages [4–6]. According to Kron et al. (2008) [4], there is a "developmental window" in the V4 stage, in which plants subjected to water stress improve the subsequent drought-stress tolerance in the reproductive stage. Furthermore, Sinclair et al. (2010) [6], using a simulation model, determined that water conservation through an early decrease in stomatal conductance and reduced transpiration rate explains the increase in soybean yield throughout 70 years with drought events.

**Citation:** Simondi, S.; Casaretto, E.; Quero, G.; Ceretta, S.; Bonnecarrère, V.; Borsani, O. A Simple and Accurate Method Based on a Water-Consumption Model for Phenotyping Soybean Genotypes under Hydric Deficit Conditions. *Agronomy* **2022**, *12*, 575. https://doi.org/10.3390/ agronomy12030575

Academic Editor: Yuan-Ming Zhang

Received: 14 December 2021 Accepted: 24 February 2022 Published: 25 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Soybean breeding programs are mostly focused on increasing total yield and yield stability, which are mainly affected by water shortage during the crop cycle. These programs, particularly public breeding programs, do not use drought response as selection criteria due to the high costs of the massive phenotyping equipment currently available. Most of the high-throughput phenotyping platforms are based on the diagnosis of changes in the plant's physiology, such as leaf conductance, leaf area, and root system development, using a complex and expensive system of analysis [7]. Hence, including drought tolerance as a selection criterion in breeding programs requires the development of high-throughput, precise, and low-cost phenotyping strategies [8]. However, as Valdez et al. (2013) [9] indicate, strategies could include the quantification of two crucial aspects of the plant's water budget: (i) the ability to capture more water, and (ii) the ability to conserve and use captured water more efficiently. Nowadays, it is well documented that water extraction during the key crop stages greatly informs crop performance in water restricted scenarios [9–11].

Water balance in a land crop is greatly dependent on evapotranspiration phenomena. This is the combination of two independent processes involved in water losses from the soil. One is the evaporation of the water content in the soil surface, and the second is the loss of water contained in the plant tissues by transpiration [12]. In crops, the transpiration increases as the plants grow due to the increase of leaf area. At the same time, this causes a contrary effect on the evaporation, which decreases progressively. Evapotranspiration can be determined by measuring several components of water balance in the soil. Specifically, in a controlled close system with no run-off or percolation, the total water content of the system within a certain period can be quantified as weight.

At the whole plant level and under constant water demand, water uptake or water use depends mainly on root-system development, leaf area, transpiration rate, and leaf conductance [9]. Under conditions with increasing water restriction, each of the abovementioned parameters has a different role on water use. However, there is a consensus within the scientific community that transpiration changes are a critical component in contexts where the available water content in the soil is changing [13–16]. Transpiration is determined by the demand and controlled by stomata; during soil desiccation, the water extractable by the roots is continuously reduced, which determines that full transpiration demand cannot be supported. In this situation, the plants respond with stomatal closure to avoid shoot desiccation. In addition, stomatal opening is quite sensitive to the evaporative demand, and a high vapor-pressure deficit (VPD) reduces stomatal opening to restrict water losses. Soil water content and atmospheric demand ultimately determine the dynamic of water consumption [17]. Therefore, leaf conductance measurements using gravimetric methods have been explored as robust parameters for breeding programs [18–21]. Moreover, the variation of stomatal conductance in soybeans in response to VPD and soil water content has been shown to be dependent on the genotype [22,23].

The increase in the number of new genotypes evaluated by breeding programs calls for a high capacity for data acquisition and processing, yield prediction, and responses under different environmental conditions. To this end, the use of tools such as sensors, data analysis programs, algorithms, and models is essential [24]. However, the evaluation of the models with experimental data is necessary for their improvement and adjustment [25,26]. Under optimal development conditions, mathematical models can predict plants' responses more accurately. Under stress conditions, however, a greater range of responses are generated and the prediction accuracy of the models becomes weaker [24]. The results of models focused on the water absorption of crops [27–29] show a great variation among themselves as a consequence of the many factors involved. For example, the substrate information must be improved, as well as the plant-soil interaction [24]. It is necessary to have a better theoretical and empirical basis and an appropriate knowledge of the environmental conditions where the model will be applied. Another aspect that stands out that is not always taken into account when determining the performance of the model is the competence of its users. Therefore, it is necessary to simplify and facilitate the input data models. Most of the relevant mathematical models that we know of are designed to explain crop behavior

in response to field environmental changes. Moreover, the time variable is not included in these models, which would not account for the physiological changes at the plant level.

In this study, we developed a simple and informative model based on biological and mathematical approaches for evaluating, under controlled growth conditions, the response of soybean plants to water restriction through the quantification of water consumption in the vegetative stage. To find a physiological explanation of the variables generated by the model, they were correlated with the stomatal conductance dynamic. Also, a preliminary test of the model was performed separately in a biparental segregating population and in a breeding population. Results show that the approach proposed is a valid option to be included in plant phenotyping protocols with limited manpower and infrastructure.

#### **2. Materials and Methods**

#### *2.1. Phenotyping Method Based on Water Consumption under Controlled Environmental Conditions*

For the development of the mathematical model, we performed an experiment with a soybean genotype (GENESIS 5601). Plants were grown in a 0.5 L plastic bottle (pot) filled with a mix of sand:vermiculite (1:1). This combination of plant, pot, and substrate was defined as a Plant Pot Substrate system (PPS) (Figure S1). Plants were grown in an environment defined by day/night cycle temperatures of 30/20 ◦C, respectively, and a light/darkness photoperiod of 16/8 h, respectively. Relative humidity (RH) was controlled at 35/40% during the entire growth period. Three seeds per pot were sown, and only one seedling remained after the cotyledon expanded. The homogeneity of the plants was carefully analyzed to avoid any interference related to developmental phenotype. During the first 16 days after sowing (developmental stage V2–3), soybean seedlings were grown without water restriction, and substrate was kept at field capacity with Rigaud and Puppo (1975) [30] medium supplemented with KNO<sup>3</sup> (10 mM final concentration). Since day 17, watering was suspended, and water substrate content was measured daily by gravimetry (water gravimetric content) during the next 10 days of water deficit (dwd) (Figure 1). Stomatal conductance was measured simultaneously with PPS on the abaxial leaf surface with a Porometer Model SC-1 (Decagon Device), as instructed by the manufacturer, since the suspension of watering. Five biological replicates (*n* = 5) were used for determinations, consisting of five independent PPSs. *Agronomy* **2022**, *12*, x FOR PEER REVIEW 4 of 17

**Figure 1.** Schematic illustration of the experiment for determining water consumption**.** Plants were maintained at field-capacity condition until day 16, after which the water supply was suspended. The Plant Pot Substrate system (PPS) was weighted daily for 10 days (*t* =0 until *t* =10 ). Evaporation: *E* ; transpiration: *T* ; days of water deficit: dwd. **Figure 1.** Schematic illustration of the experiment for determining water consumption. Plants were maintained at field-capacity condition until day 16, after which the water supply was suspended. The Plant Pot Substrate system (PPS) was weighted daily for 10 days (*t* = 0 until *t* = 10). Evaporation: *E*; transpiration: *T*; days of water deficit: dwd.

lines SO7.6557 x DM6.8 were phenotyped using the methodology and the mathematical model developed in this study. The phenotyping experiment was laid out in a randomized incomplete block design, with three replications in each experiment. Growth conditions were the same as described in the Section 2.1. PPS weight and stomatal conductance

A local breeding population composed of a set of 89 genotypes [31] was also phenotyped using the methodology described above and the mathematical model developed in this study. Five well-known commercial varieties were included as checks in all phenotyping experiments. In this case, plants were grown in 2.851 L PVC tubes (11 cm in diameter and 30 cm long) with a mix of sand:vermiculite (1:1) under the same environmental conditions as described previously. Plants were grown without any watering restriction for 30 days (developmental stage V4–5), after which watering was suspended and the PPS weight was registered at days 0, 4, and 8 after suspending the water supply. PPS weight

The F3 segregating population of 177 genotypes was genotyped using a SoySNP6k chip [32] in an Iscan system (Illumina, San Diego, CA, USA). SNP calling was done using the GenomeStudio software (Illumina, San Diego, CA, USA). Genotyping-by-sequencing (GBS) data were obtained for the breeding population according to the methodology pro-

*2.3. Phenotyping Strategy Applied to a Breeding Population* 

and stomatal conductance were measured simultaneously.

*2.4. Genotyping by Sequencing and SNP Calling* 

posed by Quero et al. (2021) [31].

were measured simultaneously.

#### *2.2. Phenotyping Strategy Applied to an F3 Segregating Population*

An F3 segregating population of 177 genotypes derived from the crossing of parental lines SO7.6557 × DM6.8 were phenotyped using the methodology and the mathematical model developed in this study. The phenotyping experiment was laid out in a randomized incomplete block design, with three replications in each experiment. Growth conditions were the same as described in the Section 2.1. PPS weight and stomatal conductance were measured simultaneously.

#### *2.3. Phenotyping Strategy Applied to a Breeding Population*

A local breeding population composed of a set of 89 genotypes [31] was also phenotyped using the methodology described above and the mathematical model developed in this study. Five well-known commercial varieties were included as checks in all phenotyping experiments. In this case, plants were grown in 2.851 L PVC tubes (11 cm in diameter and 30 cm long) with a mix of sand:vermiculite (1:1) under the same environmental conditions as described previously. Plants were grown without any watering restriction for 30 days (developmental stage V4–5), after which watering was suspended and the PPS weight was registered at days 0, 4, and 8 after suspending the water supply. PPS weight and stomatal conductance were measured simultaneously.

#### *2.4. Genotyping by Sequencing and SNP Calling*

The F3 segregating population of 177 genotypes was genotyped using a SoySNP6k chip [32] in an Iscan system (Illumina, San Diego, CA, USA). SNP calling was done using the GenomeStudio software (Illumina, San Diego, CA, USA). Genotyping-by-sequencing (GBS) data were obtained for the breeding population according to the methodology proposed by Quero et al. (2021) [31].

#### *2.5. Data Analysis*

#### 2.5.1. Nonlinear Models

To model the weight of the plant pot substrate system and its evapotranspiration, a onephase exponential decay and one-phase association function were used, respectively. In both cases, the initial values and parameters estimation is described by Equations (11) and (12) of Section 3.6. It is worth noting that the nonlinear models were fitted using the *nls* function of the *stats* [33] package of the statistical R software [33].

#### 2.5.2. Multivariate Characterization

The multivariate characterization of the soybean populations was accomplished through some exploratory analyses and visualization tools. First, genotypes were clustered in groups based on the PCA of the genotypic variability using a hierarchical clustering algorithm (HCPC). Genotypes were grouped based on similarity from the Euclidean distance using the Ward method, and the number of groups was determined by the highest relative loss in inertia using the function *HCPC* of the *FactoMineR* [34] package in the statistical software R [33]. Second, a correspondence analysis (CA) was performed based on the contingency tables of the HCPC analysis from the genotypic characterization, and the clustering was performed according to *t*0.5 and *Gw*(*t*0.5). The correspondence analysis was performed to visualize the relationship between the two grouping strategies (i.e., based on genetic and phenotypic variables). The CA was conducted using the *CA* function of the *FactoMineR* [34] package in the R software [33].

#### 2.5.3. Statistical Model and Adjustment of Phenotypic Means

Best linear unbiased estimators (BLUEs) for each advanced inbred line were obtained with mixed models to include experimental design components using the following model:

$$y\_{ij} = \mu + \alpha\_i + \beta\_j + \varepsilon\_{ij}$$

where *yij* is the response variable, *µ* is the overall mean or intercept, *α<sup>i</sup>* is a random variable associated with the *i*th assay with *αi*~N (0, *σ* 2 *<sup>A</sup>*), *β<sup>j</sup>* is the effect of the *j*th line, and *εij* is the residual error with *εij*~N (0, *σ* 2 *ε*). The model was adjusted in the R software [33] using the *lme4* package [35], while the BLUEs were estimated using the *emmeans* package [36], also in R [33]. In all cases, BLUEs were further used as genotypic values for the model, PCA, and CA analyses.

#### **3. Results**

In this work, we developed an empirical mathematical model for describing the kinetics of water consumption of a PPS (Figure 2) using an important plant such as soybean. This model was applied for developing a phenotyping methodology for water-deficit response in soybean plants under controlled environmental conditions. *Agronomy* **2022**, *12*, x FOR PEER REVIEW 6 of 17

**Figure 2.** (**a**) Empirical model representing evapotranspiration over time and the empirical model adjustment. Weight (*W* ) of the PPS along the days of water deficit (dwd). Sum of dry weight substrate, pot, and plant weight ( *S* ), residual water ( *AR*), evapotranspirable water ( *AET* ), and potential evapotranspiration ( *B* ) of plants at field-capacity conditions. (**b**) Mathematical analysis for modeling *ET t*( ) as a function of time. The parameter 0.5 *t* is the time required for the PPS to halve **Figure 2.** (**a**) Empirical model representing evapotranspiration over time and the empirical model adjustment. Weight (*W*) of the PPS along the days of water deficit (dwd). Sum of dry weight substrate, pot, and plant weight (*S*), residual water (*AR*), evapotranspirable water (*AET*), and potential evapotranspiration (*B*) of plants at field-capacity conditions. (**b**) Mathematical analysis for modeling *ET*(*t*) as a function of time. The parameter *t*0.5 is the time required for the PPS to halve the potentially evapotranspirable water.

#### the potentially evapotranspirable water. *3.1. Mathematical Model Development*

*3.1. Mathematical Model Development*  In the experiment, the PPS system weight (*W* ) is defined as the weight of water ( *A* ) plus the rest of the components of the system ( *S* ). *S* is the sum of substrate, pot, and plant weight. The latter was considered constant during the assay because, although it could vary throughout the days, it is insignificant in comparison to the rest of components of *S* . *A* is the sum of the transpirable water ( *AET* ) plus a percentage of that which cannot be evapotranspired by the PPS throughout the whole assay. The non-evapotranspirable water is defined as residual water ( *AR* ). The values of *AR* depend on the matric In the experiment, the PPS system weight (*W*) is defined as the weight of water (*A*) plus the rest of the components of the system (*S*). *S* is the sum of substrate, pot, and plant weight. The latter was considered constant during the assay because, although it could vary throughout the days, it is insignificant in comparison to the rest of components of *S*. *A* is the sum of the transpirable water (*AET*) plus a percentage of that which cannot be evapotranspired by the PPS throughout the whole assay. The non-evapotranspirable water is defined as residual water (*AR*). The values of *A<sup>R</sup>* depend on the matric potential of the substrate and the transpiratory capacity of the plants (Figure 2a). The amount of potentially evapotranspirable water of the PPS is a function dependent on time (*t*), and it is named *AET*(*t*). Therefore, *A*(*t*) = *AET*(*t*) + *AR*.

potential of the substrate and the transpiratory capacity of the plants (Figure 2a). The

In order to find an algebraic expression for the *W t*( ) function, we assume the following: "The velocity with which the PPS weight varies is directly proportional to the amount of water that can be evapotranspired for the PPS in this fraction of time". This hypothesis can be mathematically expressed through the following differential equation:

*kA t*

*dt* = − (2)

( ) *ET* ( ) *dW t*

#### amount of potentially evapotranspirable water of the PPS is a function dependent on time *3.2. PPS Weight Modelling over Time*

(*t)*, and it is named *AET* ( )*t .* Therefore, () () *A ET R t At A* = + . If *W*(*t*) is a function of the PPS weight over time (*t*), then

If *W t*( ) is a function of the PPS weight over time (*t*), then

$$\mathcal{W}(t) = A(t) + \mathcal{S} = A\_{ET}(t) + A\_{R} + \mathcal{S} \tag{1}$$

In order to find an algebraic expression for the *W*(*t*) function, we assume the following: "The velocity with which the PPS weight varies is directly proportional to the amount of water that can be evapotranspired for the PPS in this fraction of time". This hypothesis can be mathematically expressed through the following differential equation:

$$\frac{dW(t)}{dt} = -kA\_{ET}(t)\tag{2}$$

where *k* > 0 is a constant of proportionality. The negative sign of the equation is due to the decrease of the PPS weight over time, as the assay was performed withdrawing watering at *t* = 0.

Since *A<sup>R</sup>* and *S* are constant magnitudes over time, the derivation of the equality (1) results in *dW*(*t*) *dt* = *dAET*(*t*) *dt* , thus, the Equation (2) can be rewritten as follows:

$$\frac{dA\_{ET}(t)}{dt} = -kA\_{ET}(t)\tag{3}$$

obtaining the first order homogeneous linear differential equation for the function *AET*(*t*). Solving the Equation (3) using separation variables method (Figure S2), the solution can be written as:

$$A\_{ET}(t) = Be^{-kt} \tag{4}$$

where *B* is a positive constant.

By combining the Equations (1) and (4), we obtain the following equation:

$$W(t) = Be^{-kt} + A\_R + S \tag{5}$$

The graphic representation and the experimental data adjustment of Equation (5) are shown in Figure 2a.

#### *3.3. Evapotranspiration Modelling as a Function of Time*

To quantify the evapotranspiration (Figure 2b) (*ET*(*t*)) of the PPS from the precise moment when watering was suspended (*t* = 0) to time *t*, is the difference between the PPS weight in both times, that is

$$ET(t) = W(0) - W(t) \tag{6}$$

Combining the Equations (5) and (6), we have that:

$$\begin{array}{l} ET(t) &= W(0) - W(t) = A\_{ET}(0) + A\_R + S - (A\_{ET}(t) + A\_R + S) \\ &= A\_{ET}(0) + A\_R + S - A\_{ET}(t) - A\_R - S \\ &= A\_{ET}(0) - A\_{ET}(t) = B - Be^{-kt} \\ &= B \left(1 - e^{-kt}\right) \end{array}$$

Therefore,

$$ET(t) = B\left(1 - e^{-kt}\right) \tag{7}$$

The graphic representation and the experimental data adjusting of the Equation (7) are shown in Figure 2b.

#### *3.4. Potential Evapotranspiration Estimated by the Model*

At the moment watering was suspended (*t* = 0), the PPS had the maximum quantity of evapotranspirable water. By definition, this quantity is *AET*(0) = *B*, as observed in Figure 2a.

On the other hand, since the constant of proportionality (*k*) of Equation (7) is positive, we have that

$$\lim\_{t \to \infty} ET(t) = \lim\_{t \to \infty} B\left(1 - e^{-kt}\right) = B \tag{8}$$

That is, the parameter *B* is the horizontal asymptote of the function *ET*(*t*), and represents the potential evapotranspiration of PPS, as observed in Figure 2b.

#### *3.5. Half-Time of ET*

Half-time is, by definition, the time required for the PPS to reduce the potential evapotranspiration in half. As the potential evapotranspiration is *B*, the half-time, *t*0.5, is expressed as follows:

$$ET(t\_{0.5}) = \frac{B}{2} \tag{9}$$

Equation (7) is used to calculate *t*0.5, therefore,

$$t\_{0.5} = \frac{1}{k} \ln(2) \tag{10}$$

Thus, the half-time is inversely proportional to the constant *k*, and independent of *B*. The graphic representation of *t*0.5 is shown in Figure 2b.

#### *3.6. Parameters of the ET Model Calculated from the Experiment Data*

If *W<sup>j</sup>* is the PPS weight registered on the *<sup>j</sup>*th day since suspending watering, *j*, *W<sup>j</sup> <sup>N</sup> j*=0 is the set of data obtained from the assay (Figure 2a).

To determine the *ET* parameters as a function of time (Figure 2b), with this data set, we can use the least-squares method that requires computational resources or an analytical method from Equation (7). For this last option, we took three specific determinations, *W*0, *Wn*, and *W*2*n*, and performed them on days *0*, *n*, and *2n*, respectively. Replacing this data on the Equation (7), we can directly find the parameters, obtaining the following:

$$B = \frac{\left(\mathcal{W}\_{\text{n}} - \mathcal{W}\_{0}\right)^{2}}{\mathcal{W}\_{2\text{n}} - \mathcal{W}\_{\text{n}} + \mathcal{W}\_{0}}\tag{11}$$

$$k = -\frac{1}{n} \ln \left( \frac{\mathcal{W}\_{2n} - \mathcal{W}\_0}{\mathcal{W}\_n - \mathcal{W}\_0} \right) \tag{12}$$

To develop an empirical model, the PPS weight measurements (*W*(*t*)) during the whole water-deficit period were used to determine the curve fitting according to the experimental methods described in Figure 2. As shown in Figure 1, 10 days of water restriction determined a curve of PPS weight loss, determined by the water transpired by the plants. As indicated in Figure 2a, *W* at time *t* is defined by the parameter *B* determined by the weight when *t* is 0 and it represents the potentially evapotranspirable water, *A<sup>R</sup>* indicates the water not extractable, and *S* indicates the weight of the support and dry substrate.

The evapotranspiration of the PPS over time from the day when watering was withdrawn was determined by Equation (7) and defined by the parameters *B* and *k*. It is important to note that *B* represents the potential evapotranspiration of the PPS. By definition, half-time (*t*0.5) is the time required for the PPS to reach half of the potential evapotranspiration (Figure 2b). These parameters show the kinetics of water consumption of the genotypes and could help in the characterization of soybean genotypes.

#### *3.7. The Model Minimizes Sampling Requirements in Phenotyping Protocols*

To simplify the data collection procedure and increase the high throughput capacity, the model parameters were estimated using the minimum sampling. Figure 3a shows the values of *K* in Equation (7) estimated by the Gauss Newton numeric methods versus the values of *K* estimated by analytic methods using a sampling of the PPS weight every two days. As observed in Figure 3a, an optimum adjustment was obtained when the values of the PPS weight were sampled on days 4 and 8. The parameters estimated with the weight of these two days are closer to the best model fit. The same results were obtained with

parameter *B* of Equation (7) (Figure 3a), but with an adjustment of 0.96. Both parameters of the model are critical to evaluate the water-consumption curve of a specific genotype. Hence, the PPS weight on days 4 and 8 appears to be enough to describe the kinetics of water consumption throughout the water restriction period. Figure 3b shows that the adjustment curve of Equation (7) with parameters *B* and *K* is the same as the one obtained using all the PPS weight data in the Gauss Newton estimation method. Hence, the PPS weight on days 4 and 8 appears to be enough to describe the kinetics of water consumption throughout the water restriction period. Figure 3b shows that the adjustment curve of Equation (7) with parameters *B* and *K* is the same as the one obtained using all the PPS weight data in the Gauss Newton estimation method.

*Agronomy* **2022**, *12*, x FOR PEER REVIEW 9 of 17

inition, half-time ( 0.5 *t* ) is the time required for the PPS to reach half of the potential evap-

otranspiration (Figure 2b). These parameters show the kinetics of water consumption of

To simplify the data collection procedure and increase the high throughput capacity,

the model parameters were estimated using the minimum sampling. Figure 3a shows the values of *K* in Equation (7) estimated by the Gauss Newton numeric methods versus the

values of *K* estimated by analytic methods using a sampling of the PPS weight every two days. As observed in Figure 3a, an optimum adjustment was obtained when the values of the PPS weight were sampled on days 4 and 8. The parameters estimated with the weight of these two days are closer to the best model fit. The same results were obtained with parameter *B* of Equation (7) (Figure 3a), but with an adjustment of 0.96. Both parameters of the model are critical to evaluate the water-consumption curve of a specific genotype.

the genotypes and could help in the characterization of soybean genotypes.

*3.7. The Model Minimizes Sampling Requirements in Phenotyping Protocols* 

**Figure 3.** Minimum sampling requirements analysis and adjustment of model. (**a**) Relation between *K* or *B* estimated by the Gauss Newton method and the predicted *K* or *B* generated by analytic methods with two-time sampling, respectively. Axis indicates the parameters and the days of the sampling. (**b**) *ET* over time. Graphic representation of Equation (7) adjusted according to the data **Figure 3.** Minimum sampling requirements analysis and adjustment of model. (**a**) Relation between *K* or *B* estimated by the Gauss Newton method and the predicted *K* or *B* generated by analytic methods with two-time sampling, respectively. Axis indicates the parameters and the days of the sampling. (**b**) *ET* over time. Graphic representation of Equation (7) adjusted according to the data of the numeric (solid curve) or analytic methods (dashed colored curve). Colors indicate the two specific sampling days used for modeling the curve. Days of water deficit: dwd.

of the numeric (solid curve) or analytic methods (dashed colored curve). Colors indicate the two

#### *3.8. Stomatal Conductance as a Function of Time* chamber, or the water vapor exiting through the stomata pore of the leaf. Since most of

*3.8. Stomatal Conductance as a Function of Time* 

*Agronomy* **2022**, *12*, x FOR PEER REVIEW 10 of 17

By definition, stomatal conductance (*Gw*) is the rate of CO<sup>2</sup> entering the substomatal chamber, or the water vapor exiting through the stomata pore of the leaf. Since most of the water lost by the PPS is from stomatal transpiration, it is possible to suppose that the evaporation is negligible. Thus, the variation in weight can be attributed exclusively to transpiration, that is: the water lost by the PPS is from stomatal transpiration, it is possible to suppose that the evaporation is negligible. Thus, the variation in weight can be attributed exclusively to transpiration, that is: ( ) *<sup>w</sup>* ( ) *dW t G t dt* = −τ(13)

By definition, stomatal conductance (*Gw* ) is the rate of CO2 entering the substomatal

$$\frac{d\mathcal{W}(t)}{dt} = -\tau \mathcal{G}\_{\mathcal{W}}(t) \tag{13}$$

where in *τ* represents a constant of proportion. By deriving the Equation (5) and replacing in (13), we obtain the following:

By deriving the Equation (5) and replacing in (13), we obtain the following: *<sup>k</sup> G t Be*

$$\mathcal{G}\_{\mathcal{W}}(t) = \frac{k}{\pi} B e^{-kt} \tag{14}$$

The graphic representation and the experimental data adjustment of the Equation (14) are shown in Figure 4a. (14) are shown in Figure 4a.

The graphic representation and the experimental data adjustment of the Equation

**Figure 4.** Modeling of stomatal conductance in response to water availability. (**a**) Empirical model representing the conductance over time and the empirical model adjustment. (**b**) Lineal regression between conductance and substrate weight. Gw: stomatal conductance; W: PPS weight; dwd: days of water deficit. **Figure 4.** Modeling of stomatal conductance in response to water availability. (**a**) Empirical model representing the conductance over time and the empirical model adjustment. (**b**) Lineal regression between conductance and substrate weight. *Gw*: stomatal conductance; *W*: PPS weight; dwd: days of water deficit.

*3.9. Conductance as a Function of PPS Weight*

ττ

*3.9. Conductance as a Function of PPS Weight*  By combining the Equation (5) with (14), we obtain the following:

$$G\_W(t) = \frac{k}{\tau} B e^{-kt} = \frac{k}{\tau} (\mathcal{W}(t) - A\_R - S) = \frac{k}{\tau} \mathcal{W}(t) + \frac{k}{\tau} (-A\_R - S)$$
 
$$\text{Thus}$$

*w RR*

Thus,

Thus,

$$G\_w = \frac{k}{\tau} \mathcal{W} + a \tag{15}$$

where in *a* = *<sup>k</sup> τ* (−*A<sup>R</sup>* − *S*).

*w <sup>k</sup> G Wa* τ = + (15) *<sup>k</sup> a AS* Thus, stomata conductance is a linear function of the PPS weight, as shown in Figure 4b. To understand the physiological component of the kinetics of water consumption dissected by the proposed model, *G<sup>w</sup>* over time was included in the mathematical modeling.

where in ( ) *<sup>R</sup>* τ = −− . Thus, stomata conductance is a linear function of the PPS weight, as shown in Figure 4b. A modeled curve of *G<sup>w</sup>* over the time was performed from Equations (5) and (13) (described in Sections 3.2 and 3.9). Figure 4a shows how the model generated by theoretical tools is adjusted with the experimental data. Using function (4), it is possible to estimate

the half-time for *Gw*(*t*0.5) from the moment watering is suspended. This parameter is relevant as a useful indicator to identify contrasting responses of different genotypes to the hydric deficit. Figure 4b shows the direct theoretical relation between the *G<sup>w</sup>* obtained from Equation (15) and the PPS weight (red curve), and the high correlation with the experimental data (*R* <sup>2</sup> = 0.92).

#### *3.10. Application of the Phenotyping Methodology in Two Breeding Populations at Different Plant Developmental Stages*

The methodology was tested in a phenotyping approach in two breeding populations. The parameters of the water-consumption curve defined by the model (*t*0.5 and *Gw*(*t*0.5)) were analyzed and related to the genetic variability in order to improve the understanding of the plant's water-consumption behavior in soybean.

#### 3.10.1. Half Time and Stomatal Conductance

An analysis of the kinetics of water consumption in response to a progressive water deficit was performed in a biparental population. Two parameters, *t*0.5 and *Gw*(*t*0.5), were selected to characterize the variability of the recombinant genotypes (Figure 5a,b). Values of *t*0.5 between 2 and 6 days showed high variability in the kinetics of water consumption. However, 25% of the genotypes had a *t*0.5 lower than 3 days, and the remaining 75% had a *t*0.5 higher than 3 days. The normal distribution observed confirms that the values of the parameters had a biological behavior inside the population. On the other hand, the distribution of *Gw*(*t*0.5) ranked between 101.8 and 202.8 mmol H2O m<sup>2</sup> s −1 shows normal distribution that accomplishes the behavior of a water-consumption response (Figure 5b).

When the phenotyping protocol was evaluated in V5 plants of a breeding population, a similar range of *t*0.5 values were obtained in comparison to those obtained in the biparental population when evaluating in V3 plants (Figure 5a,c). However, when *Gw*(*t*0.5) was analyzed, a wider range of values was obtained (6–300), showing that this parameter is more affected by the developmental stage and genotypic variability (Figure 5b,d).

It is important to point out that *Gw*(*t*0.5) is the conductance reached at the time when the plant has consumed half of the potentially evapotranspirable water, and not a simple calculation of half *G<sup>w</sup>* at the initial time (*t*0). This confirms the idea that parameters identified by the model are biologically relevant and, at least in soybean, that could help in the analysis of plant response to water deficit.

#### 3.10.2. Genetic Structure and Correspondence Analysis

Based on the genetic structure of both breeding populations, we identified three main groups in the biparental population and six groups in the breeding population (Figure 6a,b). The latter shows a higher genetic diversity, which would explain the results of the analysis in the distribution of values for the parameters of the model (Figure 5). The correspondence analysis between HCPC genotypic and HPC phenotypic groups showed the following correspondence: for the *t*0.5 parameter and in the biparental population; tq3 and g2; and tq1 and g3. The relation between tq2 and g1 is not clear. In the case of the breeding population, the correspondence indicated the relation between tq1 and g3; tq2 and g2; and g6, tq3, and g4. It was not possible to find a correspondence for g1 and g5 (Figure 6c,d). When the same correspondence analysis was performed between the genetic groups and the *G<sup>w</sup>* parameters, the results showed correspondences between Gwq3 and g2, and between Gwq2 and g1 in the biparental population. It was not possible to find a relation between Gwq1 and g3 (Figure 6e). In the case of the breeding population, correspondences between Gwq3, g3 and g6; Gwq2 and g1; and g2 and g5 were found (Figure 6f). It was not possible to find a relation between Gwq1 and g4. This type of analysis indicates that the changes in water consumption under hydric conditions, quantified by the parameters of the model, are related to the genetic components in the genotype.

**Figure 5.** Graphic distribution of t0.5 and *Gw*(*t*0.5) values obtained from two soybean populations. t0.5 is expressed in days and *Gw*(*t*0.5) in mmol H2O m2 s−1. A frequency analysis was performed in 177 and 89 genotypes of a biparental and breeding population, respectively. The darker regions represents the 25% of the genotypes with superior and inferior *t*0.5 and *Gw*(*t*0.5) values of the populations. **Figure 5.** Graphic distribution of *t*0.5 and *Gw*(*t*0.5) values obtained from two soybean populations. *t*0.5 is expressed in days and *Gw*(*t*0.5) in mmol H2O m<sup>2</sup> s −1 . A frequency analysis was performed in 177 and 89 genotypes of a biparental and breeding population, respectively. The darker regions represents the 25% of the genotypes with superior and inferior *t*0.5 and *Gw*(*t*0.5) values of the populations. (**a**,**b**) biparental population; (**c**,**d**) breeding population.

When the same correspondence analysis was performed between the genetic groups and the Gw parameters, the results showed correspondences between Gwq3 and g2, and between Gwq2 and g1 in the biparental population. It was not possible to find a relation between Gwq1 and g3 (Figure 6e). In the case of the breeding population, correspondences between Gwq3, g3 and g6; Gwq2 and g1; and g2 and g5 were found (Figure 6f). It was not possible to find a relation between Gwq1 and g4. This type of analysis indicates that the

Based on the genetic structure of both breeding populations, we identified three main

groups in the biparental population and six groups in the breeding population (Figure 6a,b). The latter shows a higher genetic diversity, which would explain the results of the analysis in the distribution of values for the parameters of the model (Figure 5). The correspondence analysis between HCPC genotypic and HPC phenotypic groups showed the following correspondence: for the *t*0.5 parameter and in the biparental population; tq3 and g2; and tq1 and g3. The relation between tq2 and g1 is not clear. In the case of the breeding

(**a**,**b**) biparental population; (**c**,**d**) breeding population.

3.10.2. Genetic Structure and Correspondence Analysis

**4. Discussion** 

the model, are related to the genetic components in the genotype.

**Figure 6.** Population structure calculated from genotypic data (molecular markers), phenotypic data (parameters of the model), and the correspondence analysis. The groups were obtained by hierarchical clustering on the principal component (HCPCgenotypic). (**a**) Biparental population: g1, g2, and g3 groups. (**b**) Breeding population: g1, g2, g3, g4, g5, and g6 groups. Correspondence analysis of the contingency table between groups based on the parameters of the models and groups according to genotype. (**c**,**e**) Biparental population. (**d**,**f**) Breeding population. Population structure calculated from parameters *t*0.5 and *Gw* of the model. tqn and Gwqn were obtained by hierarchical clustering on the principal component (HCPCphenotypic). **Figure 6.** Population structure calculated from genotypic data (molecular markers), phenotypic data (parameters of the model), and the correspondence analysis. The groups were obtained by hierarchical clustering on the principal component (HCPCgenotypic). (**a**) Biparental population: g1, g2, and g3 groups. (**b**) Breeding population: g1, g2, g3, g4, g5, and g6 groups. Correspondence analysis of the contingency table between groups based on the parameters of the models and groups according to genotype. (**c**,**e**) Biparental population. (**d**,**f**) Breeding population. Population structure calculated from parameters *t*0.5 and *G<sup>w</sup>* of the model. tqn and Gwqn were obtained by hierarchical clustering on the principal component (HCPCphenotypic).

to develop varieties with a better ability to cope with periods of water shortage. The com-

#### **4. Discussion**

In the current climate change scenario, soybean breeders must use effective strategies to develop varieties with a better ability to cope with periods of water shortage. The complexity of drought-tolerance traits has prevented the development of successful and accessible phenotyping strategies for selection, especially in small-scale breeding programs. This situation motivated our work to develop an effective phenotyping strategy. The simplification of phenotyping strategies becomes necessary when a high number of plants must be evaluated at the same time. In this regard, considerable efforts were made to reach this objective. In our case, we set the focus on the development of a model able to characterize and predict the water-consumption curve with low sampling requirements. Since the plant growth system minimizes water losses by evaporation, the water-consumption curve could be related to the transpiration curve, while also including the stomatal response as a parameter of the model. As demonstrated in several studies, water consumption using gravimetric methods correlates with the measurement of the transpiration rate under specific VPD conditions [37], so a fairly high throughput analysis based on that concept could be applied.

The data confirms the strong relation between the kinetics of water consumption and stomatal conductance, which means that the water transpired from the PPS defined in the study is regulated by stomatal conductance, and this last variable is regulated by the water availability in the pot. The model generated in this work contemplates the time variable, therefore it considers the changes in the parameters of water consumption. These types of models included dynamic variables which are scarce and would allow more precise characterization of complex biological phenomena such as plant adaptation to changes in the levels of available water. In addition, the data input of the model is based on a simple variable and easily quantifiable as the weight of PPS.

Moreover, the model can predict the value of water when the conductance is 0 and define the limit of water extraction. This trait has been identified as an important factor, because genotypes with high values of water thresholds begin to partially close their stomata at a relatively high water content, therefore saving soil water [9]. A study using data from different regions and years of the USA has shown, through simulation tools, that this trait would lead to a significant increase in soybean yield, especially in crop seasons classified as dry [6]. An early and accurate screening of the genotypes with specific responses in the water-consumption curve under water deficit seems to be an interesting advantage to be explored in plant breeding programs.

How drought episodes are established in field conditions varies depending on the agro-climatic region, the rain patterns, the soil characteristics, and the atmospheric demand. Under alternating drought conditions, in which there is a frequent period of stress alleviation, genotypes with high evapotranspiration capacity and water extraction (higher *B* and lower *AR*) could be more interesting than those with a contrary response (lower *B* and high *AR*). However, in a drought situation with a low probability of soil water recovery, the selection of genotypes with low *B* and high *A<sup>R</sup>* could be the main objective for breeding.

The gravimetric measurements of the transpiration-curves trait performed under different VPD could lead to an increase in the adaptation capacity of several crop genotypes to different environments, using the proposed method with a high scale-up potential. Phenotyping protocols have been used in different ways to classify genotypes in response to drought. In this work, we propose the assignment of some parameters of the model as a specific trait of genotypes, so these parameters should be included in drought phenotyping methods. Moreover, associations among the variables of the model increase the possibility of identifying other more informative parameters related to plant-drought tolerance. For example, the model demonstrated that the stomatal conductance could be included as an explicative variable of plant response to progressive water deficit. It is clear that regulating the speed of stomatal closure is a key element in the drought response [9]. In this context, *Gw*(*t*0.5) appears to be an effective measurement to explain the responsiveness of genotypes to changes in water availability.

Moreover, the phenotype of soybean plants in response to the hydric deficit from two breeding populations could be characterized by the proposed model. Since the phenotype of each genotype is determined by the genome-environment interaction [38,39], the variables generated by the model *t*0.5 and *Gw*(*t*0.5) were subjected to a correspondence analysis. The analysis confirmed a possible genetic association between the response of genotypes to water restriction and the parameter of the model. As expected, a clearer grouping association was observed in the biparental population than in the breeding population, where the genetic diversity of the germplasm was greater. However, the different number of genotypes included in both analyses could explain the results in the correspondence analysis. Moreover, a genome-wide association analysis could contribute to identifying the specific genomic region and genetic marker associated with the parameters defined by the model.

A recent report has classified the same genotype collection using indexes such as *drought susceptibility index* (DSI) and *yield stability index* (YSI) in relation to the crop cycle group [31], and the authors were able to identify QTLs associated with those traits. A preliminary but not confirmatory correspondence analysis between the parameter of the model *t*0.5 and DSI showed a grouping (Figure S3), indicating that the phenotyping strategy discussed in this work could be useful for improving genotype selection in relation to the performance at field conditions. However, a specific validation assay of waterdeficit response in field conditions of a set of genotypes previously characterized by the model represents the next challenge. This point is critical for proposing the phenotyping methodology as a tool to be included in crop breeding programs, especially in those with low-income support.

#### **5. Conclusions**

The developed model in the study characterizes and predicts, using gravimetric measurements and with low sampling requirements, the water-consumption curve of soybean plants when watering is withdrawn. The model confirms the strong relationship between the kinetics of water consumption and stomatal conductance. The correspondence analysis between model parameters and the response of the genotypes to water restriction in two different soybean breeding populations confirmed a possible genetic association between parameters of the model and genotypic identity. A preliminary approximation shows that the phenotyping methodology presented in this work could be included in crop breeding programs, especially in those with low-income support.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/agronomy12030575/s1, Figure S1: Experimental setup and plant growth system. Figure S2: Solution of the first-order linear homogeneous differential Equation (3). Figure S3: Cluster analysis based on drought susceptibility index (DSI) and the parameter of the model *t*0.5.

**Author Contributions:** S.S. performed the theoretical mathematical model. O.B., V.B., E.C., G.Q. and S.C. were involved in the planning of the work. G.Q. conducted all data analyses of the experiments. E.C. and G.Q. performed the phenotypic evaluation. S.C. generated the breeding populations. O.B. and V.B. wrote the manuscript. All authors corrected the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the following projects: CSIC-UdelaR Grupo 418 Estrés abiótico en plantas; Innovagro FSA\_1\_2013\_1\_12924, funded by Agencia Nacional de Investigación e Innovación; and Red Nacional de Biotecnología Agrícola, RTS\_1\_2014\_1, funded by Agencia Nacional de Investigación e Innovación, Instituto Nacional de Investigación Agropecuaria, Barraca Erro S.A., Lebu SRL, Fadisol SA, CALMER and COPAGRAN.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing is not applicable to this article as all new created data is already contained within this article.

**Acknowledgments:** E.C. would like to thank Comisión Academica de Posgrado (Universidad de la República) for granting him a Doctoral fellowship.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**

