Phenotyping and Quantitative Trait Locus Analysis for the Limited Transpiration Trait in an Upper-Mid South Soybean Recombinant Inbred Line Population (“Jackson” × “KS4895”): High Throughput Aquaporin Inhibitor Screening

Soybean is most often grown under rainfed conditions and negatively impacted by drought stress in the upper mid-south of the United States. Therefore, identification of drought-tolerance traits and their corresponding genetic components are required to minimize drought impacts on productivity. Limited transpiration (TRlim) under high vapor pressure deficit (VPD) is one trait that can help conserve soybean water-use during late-season drought. The main research objective was to evaluate a recombinant inbred line (RIL) population, from crossing two mid-south soybean lines (“Jackson” × “KS4895”), using a high-throughput technique with an aquaporin inhibitor, AgNO3, for the TRlim trait. A secondary objective was to undertake a genetic marker/quantitative trait locus (QTL) genetic analysis using the AgNO3 phenotyping results. A set of 122 soybean genotypes (120-RILs and parents) were grown in controlled environments (32/25-d/n °C). The transpiration rate (TR) responses of derooted soybean shoots before and after application of AgNO3 were measured under 37°C and >3.0 kPa VPD. Then, the decrease in transpiration rate (DTR) for each genotype was determined. Based on DTR rate, a diverse group (slow, moderate, and high wilting) of 26 RILs were selected and tested for the whole plant TRs under varying levels of VPD (0.0–4.0 kPa) at 32 and 37°C. The phenotyping results showed that 88% of slow, 50% of moderate, and 11% of high wilting genotypes expressed the TRlim trait at 32°C and 43, 10, and 0% at 37°C, respectively. Genetic mapping with the phenotypic data we collected revealed three QTL across two chromosomes, two associated with TRlim traits and one associated with leaf temperature. Analysis of Gene Ontologies of genes within QTL regions identified several intriguing candidate genes, including one gene that when overexpressed had previously been shown to confer enhanced tolerance to abiotic stress. Collectively these results will inform and guide ongoing efforts to understand how to deploy genetic tolerance for drought stress.


INTRODUCTION
Plant photosynthesis and transpirational rates are highly coupled (Sinclair, 2017). Most of the water lost by plants under stress is due to transpiration, which is linked to stomatal opening to allow CO 2 diffusion (Devi and Reddy, 2018a;Shekoofa and Sinclair, 2020). Hence, water loss due to transpiration is linked to plant biomass growth and yield (Tanner and Sinclair, 1983;Blum, 2009). Low water availability reduces assimilate partitioning to reproductive sinks, and on a cellular level impairs cell growth and division (Blum, 2011;Sarkar, 2020;Bennett et al., 2021;Sarkar et al., 2021a). Water related stresses, such as direct drought stress and associated factors including high temperature and evaporative demand, reduce both transpiration and photosynthesis, resulting in reduced crop yield (Sinclair, 2017;Balota et al., 2021).
Decreases in transpiration rate and stomatal conductance are affected by leaf temperature (LT) (Gates, 1964;Blum, 2011). Gates (1964) argued convincingly that increase in LT can be lethal for plants and that transpiration is required to cool the leaves. The study also argues that even a slow rate of transpiration can dissipate enough heat from leaves to have a huge impact on photosynthesis and growth of plants. This implies that plants with slower rate of transpiration (with subsequently warmer leaves than plants with higher rates of transpiration) could be better drought-tolerant by conserving moisture during low water stress, which would be available later in the season when Abbreviations: TR lim , limited transpiration; VPD, vapor pressure deficit; RIL, recombinant inbred line; QTL, quantitative trait locus; TR, transpiration rate; DTR, decrease in transpiration rate; LT, leaf temperature; AQPs, aquaporins; TRH 2 O, transpiration in water; TRAgNO 3 , transpiration rate following exposure to silver nitrate; PCA, principle component analysis; H 2 , broad sense heritability; DAP, days after planting; BP, breakpoint. evaporative demand is higher. Under field conditions, droughttolerant peanut genotypes (with partial stomata closure trait) displayed a downward LT slope with an increase in VPD during midday (1100 to 1400 EST) (Balota and Sarkar, 2020). Bai and Purcell (2018) observed an interaction between LT and slow and fast-wilting soyabean genotypes, where slow-wilting genotypes had a lower LT during water-deficit stress.
It is essential that any putative water saving trait should result in monetarily relevant yield benefit and should have genetic variability within the target crop species (Sadok and Sinclair, 2011). In fact, soybean genotype PI 416937 has been identified as expressing the TR lim trait (i.e., slow-wilting) phenotype in the field under high VPD (>2.5 kPa), and this has been traced to low leaf hydraulic conductance (Sinclair et al., 2008). Sadok and Sinclair (2012) suggested that the low leaf hydraulic conductance in PI 416937 might be related to a unique population of aquaporins (AQPs) in its leaves. Further studies have shown that soybean genotypes expressing TR lim trait are insensitive to aquaporin inhibitors, such as silver and zinc ions Sinclair, 2010, 2012;Devi et al., 2016).
In plants, AQPs occur in multiple isoforms in both plasmalemma and tonoplast membranes resulting in regulation of water flow in and out of cells. Physiological and molecular studies have identified AQPs as playing key roles in regulating hydraulic conductance in leaves and roots (Beaudette et al., 2007;Shekoofa and Sinclair, 2018;Sung et al., 2021). Therefore, the high-throughput phenotyping of structured soybean populations, along with molecular genotyping is needed for marker development. This type of research is essential for making the soybean genome sequence useful for breeding purposes, in particular for screening drought-tolerant soybean lines based on their sensitivity to aquaporin inhibitors.
In this work, we used a two-tiered screening method to identify soybean genotypes with the TR lim trait. Our method is similar to the three-tiered approach by Sinclair et al. (2000). Direct measures of the transpiration response rate to increasing VPD are low-throughput, tedious, and require specialized equipment. Furthermore, the number of genotypes that can be directly phenotyped for the TR lim trait is very limited. Thus, an indirect measurement of TR that may be less accurate but can allow characterization of a large number of genotypes would be helpful in breeding for the TR lim trait (Choudhary and Sinclair, 2014). Choudhary and Sinclair (2014) suggested that one possibility is that the response in TR to feeding of chemical inhibitors to leaves or roots of plants might allow discrimination among genotypes under varying VPD levels. Appropriate parents have to be chosen for crossing when it comes to screening for drought tolerance traits. Previous studies have reported that soybean genotypes "KS4895" and "Jackson" differ in drought tolerance traits such as delayed wilting, ureide, and nitrogen concentration Charlson et al., 2009;Hwang et al., 2013Hwang et al., , 2015aKing et al., 2014). Therefore, the objectives of this study were to (a) evaluate 120 recombinant inbred lines (RILs; F3 and F5 derived from "KS4895" × "Jackson") through a high-throughput phenotyping technique while measuring the RILs transpiration sensitivity rate to aquaporin inhibitor (AgNO 3 ) under high VPD condition (initial/indirect screening for TR lim ) (b) directly measure the TR of 26 selected RILs from objective "a" and categorize them for the TR lim trait under high VPD and temperature (direct screening for TR lim ), and finally (c) undertake genetics/quantitative trait loci (QTLs) analysis for the RIL population in objective "a" and associate the potential QTLs with the drought tolerance trait, i.e., TR lim .

Plant Materials and Study Site
The experiments were conducted at the University of Tennessee West Tennessee Research and Education Center (WTREC), in Jackson, TN between 2018 and 2020. A group of 122 soybean genotypes including 120 RILs and their parents "Jackson" and "KS4895" that were tested during the study ( Table 1) were provided by Dr. Larry Purcell (Hwang et al., 2015a). Genotype "KS4895" (PI 595081) is a maturity group (MG) IV developed in Kansas (Schapaugh and Dille, 1998), and "Jackson" (PI 548657) is an MG VII genotype developed by the USDA-ARS in North Carolina (Johnson, 1958). Due to greenhouse space limitations, and to increase the accuracy during data collection since a large number of bottles were involved in each weighing process, only a subset of genotypes was tested with the silver nitrate test. Thus, the RILs were grown in six sets under controlled environments in a greenhouse at 33-35 • C and 35-50% humidity during day, and 20-22 • C and 45-60% humidity during the night. Each set was grown for 4 weeks and included both parents. Five soybean seeds were planted in 3-L pots filled with soil (commercially available Miracle-Gro potting mix) and inoculated with Bradyrhizobium japonicum (Verdesian Life Sciences, Cary, NC, United States). They were thinned to three plants per pot after 1 week in three-pot replicates per genotype. Plants were maintained in a well-watered condition, in 250 mL pots (i.e., pot capacity) during the initial growth pretreatment period. After approximately 4 weeks, the plants were ready for AgNO 3 test. For the measurement of TR response to AQPs inhibitor, AgNO 3 , the technique explained by Sadok and Sinclair (2010), Shekoofa et al. (2013), and Devi et al. (2016) was followed.

Phenotyping and Phenotypic Traits
The evening before the application of AgNO 3 , three replicate plants per soybean genotype were removed from the soil and derooted by cutting the base of the plant stem underwater using a sharp blade ( Table 1). Then derooted soybean shoots were immediately placed in 150-mL Erlenmeyer flasks filled with deionized water. The shoots were kept in a dark laboratory room overnight for about 14 to 15 h with temperature maintained at 24 • C. The following morning, the derooted shoots were transferred to another 150-mL Erlenmeyer flask containing fresh deionized water and sealed with Parafilm-"M" (Pechiney Plastic Packaging, Chicago, IL, United States) to avoid direct water evaporation from the flask. Then, in each set of experiment, the derooted shoots were moved to a greenhouse in which both temperature and VPD were high ( Table 1).
Water-cooled lamps provided the photosynthetic photon flux density, which was 600 µmol m −2 s −1 at plant level (Mullen and Koller, 1988;Devi et al., 2016). The derooted shoots were allowed to acclimatize for 60 min and then the flasks plus derooted shoots were weighed (TR H 2 O initial). After another 60 min, the flasks plus derooted shoots were reweighed (TR H 2 O final), and the difference between the two weights divided by the time interval was used to calculate the TR in water (TR H2O ). Following the second weighing of the derooted shoots in water, the individual shoots were quickly transferred to dark-brown, 30-mL glass bottles for exposure to the solution of 200 µM, AgNO 3 aquaporin inhibitor (Sadok and Sinclair, 2010). The AgNO 3 solution was freshly prepared on the day before each set of experiment. The AgNO 3 solution concentration of 200 µM was selected because this concentration was the lowest concentration resulting in maximum transpiration decrease (Sadok and Sinclair, 2010). The derooted shoots were allowed to take up the AgNO 3 solution for 60 min by which time TR of the shoots had again reached a constant value (Sadok and Sinclair, 2010).
After the 60 min period, all bottles were weighed to get an initial weight for the transpiration measurement following exposure to AgNO 3 (TR AgNO 3 initial). After 180 min, each bottle was again weighed (TR AgNO 3 final). Measurements of TR for each shoot were generally completed in about 240 min after AgNO 3 treatment. The TR following exposure to silver (TR AgNO3 ) was calculated based on the difference between these two weights divided by the time interval. The difference between TR H2O and TR AgNO3 were used to quantify decrease in TR using various arithmetic combinations ( Table 2).

Leaf Temperature
In each set, LT was measured for all tested derooted shoots using a FLIR C3 Thermal Camera with WiFi (Teledyne FLIR LLC, Wilsonville, OR, United States) between 1300 and 1500 CST. In each set of experiment, the corresponding thermal images were taken 1 m above the plant at a resolution of 0.9 cm per pixel for three replications of derooted shoots. The temperature at the central point of a fully expanded trifoliate leaf was obtained for each of the three plants from a thermal image using the camera software. The thermal data were collected after the last stage of exposure to AgNO 3 when the final weighing was done. These data sets were used for calculating the average LT for each genotype.

Statistical Analysis
A mixed linear model was applied to TR response data using the Fit Model function of JMP 14 (SAS Institute) using an input of 120 RIL lines (and two parental lines), for which we had both genotypic and phenotypic information. For mixed linear models, three factors were included: (1) genotype; (2) experiment; and (3) replication within each experiment. Only genotype was considered fixed, whereas other factors were considered random, and replication was nested within experiment. Principle component analysis (PCA) was run to create a PCA biplot of phenotypic measurements. Pearson's correlation matrix was used to measure correlation among phenotypes. Polynomial regression was used to explain the variation of relative decrease in TR due to LT.

Broad Sense Heritability
Broad sense heritability (H 2 ) of all traits was calculated as the ratio of genotypic variance (σ 2 G ) by phenotypic variance (σ 2 P ) (Phansak et al., 2016). Variance was calculated as the ratio of total sum of squares (TSS) to population size (n). The heritability data have been included in Supplementary File 1.

Experiment II Direct Measurement of Transpiration Under High Vapor Pressure Deficit
Twenty-six soybean genotypes (24 RILs and both parents), based on RDTR K and RNDTR K (from the indirect measurement of transpiration, "Experiment I"), were selected to be tested for the direct TR measurement under varying levels of VPD (0.5 to 4.0 kPa). Soybean genotypes that had RDTR K and RNDTR K values from 0.00 to 0.50 were considered slow wilting, 0.51 to 1.00 were considered moderate wilting, and above 1.00 were considered as fast wilting. The plants were grown in pots constructed from polyvinyl chloride pipe (100-mm diameter and 200-mm long). The bottom of each pot was fitted with a flat end cap, in which a small hole was drilled to allow drainage of excess water. A toilet flange was attached to the top of the pot to allow easy attachment of a VPD chamber during measurements (Fletcher et al., 2007;Shekoofa et al., 2013). The pots were filled with commercially available Miracle-Gro potting mix. Five seeds per pot were sown and inoculated with B. japonicum (Verdesian Life Sciences, Cary, NC, United States). After 1 week, the plants were thinned to one plant per pot. Plants were fertilized with 200 mL of 0.075% V/V liquid fertilizer (0-10-10, N-P 2 O 5 -K 2 O, GH Inc., Sebastopol, CA, United States) at planting and again at 14 and 24 days after planting (DAP). Plants were grown for 28 days under well-watered conditions, with greenhouse temperatures regulated at 33 • C day/26 • C night. After approximately 4 weeks, four replicate plants of each genotype were selected and moved into a walk-in growth chamber and enclosed in individual humidity-regulated chambers (21 L). As there were 12 individual chambers, this meant that the genotypes had to be measured in batches that included three genotypes in each batch. Measurements of transpiration response to VPD were carried out at 32 and 37 • C. Each individual VPD chamber (21 L) was fitted with a 12-V, 80-mm-diameter cooling fan (Masscool) to continuously stir the air inside the chamber. Stirring of the air helped to maintain plant temperature near ambient air temperature within the VPD chamber. A humidity/temperature data logger (Lascar Electronics) was mounted through the sidewall of each container to monitor the environmental

RDTR K
Relative decrease in transpiration rate by "KS4895" DTR DTR of"KS4895

RRTR J
Relative ratio of transpiration rate by "Jackson" RTR RTR of"Jackson

RRTR K
Relative ratio of transpiration rate by "KS4895" RTR RTR of"KS4895

RNDTR J
Relative normalized decrease in transpiration rate by "Jackson" NDTR NDTR of"Jackson

RNDTR K
Relative normalized decrease in transpiration rate by "KS4895" conditions of the chamber. The plants were subjected to three VPD levels: low (0-1.5 kPa), medium (1.5-2.5 kPa), and high (2.5-4 kPa), on two consecutive days for each temperature (i.e., 32 and 37 • C). The humidity in the VPD chambers was obtained by adjusting the airflow rate through each chamber; in the case of the highest VPD treatment (2.5-4 kPa), the air was also initially flowed through a column of silica gel to dry the input air (Fletcher et al., 2007).
The observations of all plants for each genotype were combined for a two-segment linear regression analysis (Prism 8.0, GraphPad, Software Inc., San Diego, CA, United States) of TR versus varying VPD at both levels 32 and 37 • C temperatures.

Statistical Analysis
The outputs of a successful regression fit to the two-segment model were the coefficients defining two intersecting linear regressions and the VPD value at the breakpoint (BP) between the two linear segments: If VPD < BP, TR = Slope 1 (VPD) + Intercept 1.
The slopes of the two linear regressions (Slopes 1 and 2) were statistically compared to determine whether they differed significantly (p < 0.05). If the slopes differed, the double-linear regression was retained. When the slopes were not significantly different, a simple linear regression was applied to all the data (Devi et al., 2010;.

Genotyping and Quantitative Trait Locus Discovery
Genetic Map Construction DNA was isolated from lyophilized leaf tissue for 118 out of 120 RIL lines, as well as the parental lines "KS4895" and "Jackson" using a Promega AS1600 kit (Promega, Madison, WI, United States) using a Maxwell RSC instrument (Promega). Two genotypes did not germinate and the leaf tissue could not be collected. Genotypes for genetic map creation and QTL detection were obtained using the SoySNP6k iSelect BeadChip (Song et al., 2020), an Illumina Infinium R HD array. Genotypes were called using the iScan and Genome Studio software (Illumina, San Diego, CA, United States) by Dr. Qijian Song and Chuck Quigley of the USDA-ARS. Genotypic results were imported into TASSEL 5 (Bradbury et al., 2007), where nonsegregating markers were removed, missing genotypes were imputed using LinkImpute (Money et al., 2015), functionality incorporated into TASSEL 5 (settings = 30 High LD sites, 10 nearest, 10,000,000 maximum distance), and imputed genotypes converted to ABH format (AA = "Jackson" allele homozygote, AB = heterozygote, BB = "KS4895" allele homozygote) using the ABH Genotype add-in. Data were then error-corrected using the R package ABHGenotypeR (Reuscher and Furuta, 2016) with settings maxHapLength = 5 used throughout. This resulted in 2,184 total genetic markers in the final genetic map. Genotypic data were then merged with phenotypic information and imported in the QTL package in R (Broman and Sen, 2009). Distance between genetic markers was estimated using the est.map function with settings typical of an F3 population (error.prob = 0.01, overall genotypes ratios were AA = 33.0%, AB = 33.8%, and BB = 33.3%). Markers which introduced large gaps or increased genetic map length were identified using the droponemarker 1 command and manually removed if they reduced the size of the overall genetic mapping significantly.
Summary information on the genetic map is located Supplementary File 2 and the complete R/QTL file containing both genotypic and least square means (LSM) phenotypic data is located in Supplementary File 3.

Quantitative Trait Locus Analysis
Genetic mapping was performed within the R/QTL program (Broman and Sen, 2009) using two different QTL detection methods: standard interval mapping (SIM) and composite interval mapping (CIM) using the computer package R/QTL (Broman and Sen, 2009). Lod thresholds were determined from 1,000 permutation testing for each trait. Allelic effects were estimated after using "sim.geno" function (16 draws and an error probability of 0.01). QTL were used to build an additive model, which was fitted, refined, and then refitted using  the functions "makeqtl, " "refineqtl, " and "fitqtl." Results were extremely concordant between SIM and CIM; for brevity and clarity only CIM results are presented in this study.
considered as "slow wilting, " 0.51 to 1.00 as "moderate wilting, " and above 1.00 "fast wilting." Genotypes exhibited a decrease in TR slope over increasing VPD, indicating the expression of the TR lim trait, and those with constant linear slope exhibited no decrease in TR under high VPD, indicating the lack of the TR lim trait. For example, between two parents "KS4895" expressed the TR lim trait and "Jackson" did not (Figure 4). For those genotypes that expressed the TR lim trait, the point on x-axis where the slope FIGURE 3 | Regression curves of relative difference in transpiration rate response by "KS4895" (RDTR K ) and relative normalized difference in transpiration rate by "KS4895" (RNDTR K ) over the leaf temperature for 122 soybean genotypes.
FIGURE 4 | Transpiration rate response of two soybean genotypes ("Jackson" and "KS4895") over different levels of vapor pressure deficit at 32 • C.
changed was considered as the VPD breakpoint. Within slow wilting genotypes, 6 out of 7 were found to express the TR lim trait (i.e., VPD breakpoint) at high VPD under 32 • C ( Table 5). Among genotypes in moderate wilting category, 5 out of 10, and in high wilting only 1 out of 9 expressed the TR lim trait at 32 • C. Whereas, at 37 • C, the expression of TR lim trait was found to be limited only to 3 out of 7 in slow wilting, 1 out of 10 in moderate wilting, and none among fast wilting genotypes.

Genotyping and Quantitative Trait Locus Discovery
A genetic map was created from genotypic data obtained using the SoySNP6K Illumina array (Akond et al., 2013). The final genetic map was composed of 2,181 genetic markers (average 109/chromosome) encompassing a total of 2,856.4 cm (each chromosome was 142.8 cm on average) (Supplementary File 2).
Quantitative trait locus analysis was performed for four drought-related traits using LSM calculated from phenotypic measurements taken during our greenhouse experiments ( Table 6). Three QTLs were identified in total (Tables 6, 7), two QTLs for the TR traits (qTR_Gm10_1 and qTR_Gm12_1), and one associated with LT (qLT_Gm12_1) (Figure 5).
Furthermore, the effect of each QTL was estimated ( Table 6). Alleles from "KS4895" for the qTR_Gm10_1 QTL were associated with reduced transpiration relative to "Jackson" (−0.63 mg H 2 O m2 −1 s −1 ). In contrast, alleles from "Jackson" for the qTR_Gm12_1 QTL were associated with reduced transpiration relative to "KS4895" (−0.27 mg H 2 O m2 −1 s −1 ). The lowest TR (and presumably most water saving) lines had a combination of QTL from two different parents: qTR_Gm10_1 for "KS4895" allele and qTR_Gm12_1 for the "Jackson" allele ( Table 7). Lines which inherited alleles from "KS4895" for qLT_Gm12_1, which was distinct from the TR QTL, had lower leaf temperatures (0.84 • C lower on average).

Identification of Candidate Genes From Quantitative Trait Locus Regions
For genes within three detected QTL windows in this work (Table 6) using the "Williams 82" genome assembly 2 (Glyma.W82.a2.1), we examined gene ontologies (biological process, cellular component, and molecular function). Genes were termed as candidate genes ( Table 7 and Supplementary  Table 1) for the TR lim traits if they had biological process GO terms associated with any of the following terms: abscisic acid, water transport, root development, leaf senescence, jasmonic acid, heat acclimation, stomata, and/or salicylic acid. It should be noted that very few genes in soybean have direct experimental or empirical evidence for gene function, and in large part gene annotations are inferred from homologs identified by BLAST searches using genomic assemblies from other species (typically Arabidopsis). These genotypes have been selected based on "experiment 1" results. The wilting categories (slow, moderate, and fast) were determined based on the RDTR K and RNDTR K values. Values from 0.0 to 0.50 were considered slow wilting, 0.51 to 1.0 were considered moderate wilting, and above 1.0 were considered as fast wilting. Both parents were included with the RILs.
Within the qTR_Gm12_1 QTL region (Gm12:35946078-36388059; Table 7 and Supplementary Table 1), we identified eight candidate genes based on GO annotation (Supplementary Table 1). Unfortunately, none of the candidate genes has any direct evidence for gene function, with all annotations inferred based on BLAST identification of presumably homologous genes from other species. Four of the candidate genes in this region are putative transcription factors (Glyma.12G199100, Glyma.12G199200, Glyma.12G199600, and Glyma.12G203100) whose annotations are associated with abiotic stress responses. Another interesting candidate gene is a homolog of a stomatal patterning gene (Glyma.12G202700). The qTR_Gm12_1 QTL displayed significant allelic differences for TR using H 2 O ( Table 6), but no significant differences between parental alleles in the presence of the aquaporin inhibitor (TRAgNO 3 ). Although, no obvious aquaporin gene is present within the region identified in the "Williams 82" reference genome, it remains possible that a nonobvious aquaporin-encoding gene (but not identified as such either due to problems with the assembly or diverged between "Williams 82" and the parental lines) is present within the qTR_Gm12_1 QTL region, which was affected by the silver inhibitor and forms the underlying genetic basis for qTR_Gm12_1 QTL. Within the qLT_Gm12_1 QTL region associated with differential LT (Gm12:6971475.11867391, Table 7 and  Supplementary Table 1), we identified a total of 312 genes, of which 29 were termed candidate genes based on GO term and KOG annotations. Only one gene was noted to have any functional characterization; Glyma.12G104800 was identified through RNA-Seq/RT-PCR analysis as associated with response to salt stress (Sun et al., 2019). Other potential candidate genes for the QTL include an aquaporin (Glyma.12G097800); eight genes are involved in abscisic acid biosynthesis/signaling (Glyma.12G087200, Glyma.12G089200, Glyma.12G094800, Glyma.12G096100, Glyma.12G098900, Glyma.12G103100, Glyma.12G106400, and Glyma.12G108900), and a there were large number of transcription factors with annotations, which suggest involvement in abiotic stress responses (Glyma.12G087000, Glyma.12G094500, Glyma.12G094800, Glyma.12G096100, Glyma.12G098800, Glyma.12G100100, Glyma.12G100600, Glyma.12G103100, Glyma.12G104500, Glyma.12G104600, Glyma.12G104800, Glyma.12G105400, Glyma.12G105600, Glyma.12G106400, and Glyma.12G110400).
The results also indicated that AgNO 3 could be an effective way for high-throughput screening of drought tolerant soybean genotypes. Previous studies have suggested that silver nitrate (AgNO 3 ) may be an effective initial screen for the expression of the TR lim trait in a multitier screening system (Sadok and Sinclair, 2010;Choudhary and Sinclair, 2014;Shekoofa and Sinclair, 2018). In this work, the second tier (first being  AgNO 3 treatment) would involve selection of fewer genotypes for intensive and direct phenotyping for the TR response under increasing VPD conditions. Therefore, it allowed for considering not only the DTR trait, but also several other secondary traits to represent the reduction in TR response. A simple but less accurate screen that allows a large number of genotypes to be examined is a first-tier screen, followed by tiers of more sophisticated screens of decreasing numbers of genotypes (Sinclair, 2011).
Of all primary and secondary traits, those subjected to PCA, RDTR K , and RNDTR K had the best correlation with LT (Figure 2). The regression curve of RDTRK and RNDTR K with the LT shows that the slope is negligible for temperature rates below 32 • C, but it drastically changes with increasing the temperature above 32 • C (Figure 3). At optimum LT (24-32 • C), the TR response was almost constant. These findings were confirmed through the direct measurement of TR responses under high VPD conditions at 32 and 37 • C (i.e., the second-tier screening). Almost 88% of slow wilting genotypes (i.e., the least sensitive genotypes to AgNO 3 exposure), 50% of moderate wilting genotypes, and 11% of high wilting genotypes (i.e., highly sensitive genotypes to AgNO 3 exposure) expressed the TR lim trait at 32 • C. This was consistent with our hypothesis and previous studies done on other legume crops (Devi et al., 2010;Shekoofa et al., 2017). The genotypes with low sensitivity to aquaporin inhibitor, AgNO 3 , had the TR lim trait with VPD breakpoint (BP) at about 1.3-2.9 kPa ( Table 5). The TR lim at lower VPD rates (i.e., 1.3-2.9 kPa) allows those genotypes to conserve moisture so that there would be more water available in the soil during late-season drought to sustain soybeans physiological activities and complete seed fill. The plants with TR lim trait have also shown to delay wilting under rainfed conditions (Devi et al., 2010;. Although the TR response results observed at 37 • C (i.e., Experiment II) show that observations of the TR lim trait at 32 • C do not necessarily predict expression of the trait at higher temperatures, the expression of TR lim trait reduced to about 43, 10, and 0% for slow, moderate, and high wilting genotypes from Experiment I, respectively. This could be because the rest of genotypes which expressed the TR lim (i.e., VPD BP) at 32 • C started to increase TR under 37 • C to disperse excess temperature and heat. Shekoofa et al. (2016) reported that if very high temperatures are frequent, even to the point of being high enough to threaten plants with heat stress, it may be advantageous to consider corn hybrids that lose the TR lim trait at 38 • C rather than at 32 • C or other temperatures below 38 • C. However, soybean genotypes #16, #23, #73, and #87 allowed full expression of the TR lim trait under both temperatures (i.e., 32 and 37 • C) and displayed the maximum water conservation, no matter what the temperature regime was and performed better than both parents. Further evaluation of these genotypes should be done in field by simulating drought conditions.

Genotyping and Quantitative Trait Locus Discovery
Quantitative trait loci analysis for the RILs population was done to associate the potential QTLs with the drought tolerance traits. Two QTLs (qTR_Gm10_1 and qTR_Gm12_1; Table 6) overlapped with two previously identified QTLs (qSV_Gm10 and qSV_Gm12) associated with drought responses from PI 416937 (Carpentieri-Pipolo et al., 2012). QTL conditioning the sensitivity of plants to silver nitrate were located in the same chromosomal regions reported by Carpentieri-Pipolo et al. (2012). Although, from investigations of pedigrees, it does not appear that the population ("Jackson" × "KS4895") studied in the current research has any direct pedigree connection to the previously studied PI416937, but these QTL may be evidence that variation in the same underlying genes is causative.
Furthermore, it has been hypothesized that insensitivity to AgNO 3 could be related to QTLs and the trait could be heritable (Devi et al., 2016). Our study also found that all primary and secondary traits had high broad sense heritability (H 2 from 0.68 to 0.90) ( Table 8). Therefore, these drought-tolerant traits can be used for phenotypic selection of soybean genotypes in future breeding efforts.
We identified a relatively small number of candidate genes from within detected QTL regions (8, 8, and 29 for qTR_Gm10_1, qTR_Gm12_1, and qLT_Gm12_1, respectively) based on GO terms associated with abscisic acid, water transport, root development, leaf senescence, jasmonic acid, stomata, and/or salicylic acid. These candidate genes will be useful for future efforts to determine the causative genetic basis for the three QTLs we identified.

Limitations of Our Study and Future Directions
Although one of the RIL parental lines ("Jackson") has resequencing data publicly available, 3 the other parental line ("KS4895") to our knowledge has not been resequenced. Moreover, no resources for fine-level genetic mapping (e.g., near-isogenic lines) have previously been developed for detected QTL regions identified for the TR lim traits. As such, finemapping, identification of gene polymorphisms, correlation of polymorphisms with water use efficiency, and TR lim traits, and the ultimate cloning of causative polymorphisms from within QTL regions remains to future work.
Nevertheless, our results strongly suggest that advanced mapping populations can reveal QTLs for drought-tolerance traits under complicated genetic control to enhance the TR lim trait (i.e., ability to tolerate late-season drought) in a RIL soybean population from crossing "KS4895" × "Jackson." Validation of identified QTLs will be useful in molecular breeding of these favorable and informative QTL alleles for a superior cultivar with the ability to produce stable yield under water limiting conditions. It will be an interesting approach to associate these alleles with other agronomic traits, which are tightly linked to desirable drought tolerance traits to increase soybean yield and production under dryland conditions.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material. responsible for DNA isolation and quality control. Also, we thank Larry Purcell and Andy King at the University of Arkansas for providing the seeds for this research. We would furthermore like to thank Kendall Sheldon and Sara Safikhan for their assistance during data collection.