QTL for yield per plant under water deficit and well-watered conditions and drought susceptibility index in soybean (Glycine max (L.) Merr.)

Abstract

Drought has historically represented the greatest abiotic stress to adversely affect the stability of soybean (Glycine max (L.) Merr.) yields in non-irrigated field conditions. To investigate the genetic basis underlying drought tolerance in soybean, we screened the seed yield per plant (YP) in a natural diversity panel of 149 accessions under both water deficit (drought-stressed; YP-S) and well-watered (full irrigation; YP-W) conditions in field and greenhouse trials and calculated a drought susceptibility index (DSI) based on the differences in yield between treatments within each accession. A total of 19 quantitative trait loci (QTLs) were identified: eight were significantly associated with YP-W, 10 were associated with YP-S, 2 were associated with DSI and one QTL was associated with both YP-S and DSI. We then compared QTLs identified here with previously reported markers and found that these loci were located in regions associated with yield-related and other agronomic traits. In particular, YP-S-associated single nucleotide polymorphism (SNP) ss246509299 in chromosome (Chr.) 8, and YP-W-associated ss249310678 on Chr. 17 were both previously correlated with canopy wilt. Eight significant QTLs were located within eight genes. Glyma.18g092200, contained ss249600231 and annotated as GmWRKY168, is reportedly responsive to cadmium. This study helps to resolve which loci contribute to yield under drought stress in soybean, and can potentially serve as markers for selection of elite, high-yield and drought tolerant varieties.

Keywords:

• Soybean
• drought tolerance
• drought susceptibility index (DSI)
• seed yield per plant (YP)
• genome-wide association study (GWAS)

Introduction

Drought represents the single most significant environmental stress adversely affecting agricultural production on a global scale [1–2]. In particular, water deficit has been predicted to increase in both frequency and severity. Improving crop tolerance to drought stress through breeding is currently considered the most economical strategy for increasing agricultural productivity and reducing water use [3–5].

Soybean (Glycine max (L.) Merr.; 2n = 40) is one of the most protein- and oil-rich (40 and 20%, respectively) grain legumes and is widely consumed by both humans and livestock [6]. The amount of water soybean plants require during the growth season is approximately 450–700 mm, depending on the hybrid characteristics [7]. However, the sustainability of soybean yields is under threat by predicted climatic changes, with persistent droughts reported in agricultural regions of the world [8–11]. During drought, yield losses in soybean reportedly reach about 40% [12], even amounting to 80% at certain developmental stages and particularly severe episodes [13,14]. However, irrigation is unviable in many soybean-growing areas because of high costs and/or water scarcity [15]. The development of drought-tolerant soybean lines is thus one of the best ways to enhance the stability and sustainability of soybean production and therefore serves as a major target for crop improvement [16].

Drought response in crop plants is genetically and physiologically complex, and as a result, it has remained poorly understood. To address this problem, it is necessary to establish accurate phenotyping methods and to target traits for evaluating drought tolerance that are stable and robust to environmental variation. Several quantitative traits have been reported for screening the response to drought conditions. Major drought stress-related traits in soybean include seed yield under water deficit [13], relative yield (determined by the difference between rain-fed and irrigated yields) [13], drought susceptibility index [9], leaf wilting [17–21], water use efficiency (WUE; the quantity of crop dry matter accumulated per unit of soil water transpired) [11, 22, 23], relative water content [19], excised leaf water loss [19], seed number per plant and stem diameter [24]. Plants use various adaptive physiological strategies to cope with drought stress. In particular, the capacity for high seed yield in drought-affected fields represents the ultimate indicator of drought tolerance [9, 25].

Strong genotype-by-environment (GE) interactions have been frequently described in the results of drought tolerance screens [26] which, along with low reproducibility of response to water deficit under field conditions, can result in inaccurate evaluation of drought tolerance [27,28]. In addition, the process of conventional breeding for drought-tolerant crop varieties is time-consuming and labor intensive owing to the need for extensive data collection for quantitative drought tolerance traits and difficulties in selection of tolerant accessions [29]. To develop improved cultivars via marker-assisted selection, researchers need to identify genomic regions associated with drought tolerance.

Previous works in quantitative trait loci (QTL) mapping of drought-related traits in soybean have cumulatively resulted in the publication of 236 QTLs related to drought in Soybase (http://www.soybase.org), which typically focus on traits such as WUE [22, 23, 30–32], canopy wilting [33–36], drought tolerance [37] and drought index [9].

A major difficulty in using QTLs for environmentally dependent traits is their instability throughout different experiments [19]. In order to identify alleles that exert an effect across a wide range of drought conditions, more experiments are necessary that characterize responses to different drought severities, different developmental stages and different populations. Moreover, few studies have identified QTLs specifically related to yield-associated traits under drought conditions.

In the present study, we aimed to identify the genetic mechanisms underlying drought tolerance using association analysis in a germplasm diversity panel of 149 commercially available soybean cultivars. The objectives were to (i) analyze the genetic variation in traits that contribute to drought tolerance; and (ii) to identify the QTLs associated with both drought tolerance and yield under full irrigation and water deficit conditions. The major loci detected in this study can provide a useful resource for molecular breeding of soybean to improve drought tolerance in soybean molecular breeding.

Materials and methods

Plant materials and phenotyping

We constructed an association mapping panel consisting of 149 cultivars, among which 59 accessions were from Heilongjiang province, 86 were from Jilin province, 2 from Liaoning province, 1 from Beijing, China and 1 from Japan. Detailed information for these 149 accessions is presented in Supplemental Table S1. Three independent experiments were conducted in 2016. Field sites were located in the experimental station at Dunhuang (N40°08′51.35″E94°38′30.36″, average annual rainfall = 39.9 mm, cumulated irrigated soil) in Gansu province, Sanya (N18°23′30.37″E109°11′29.45″, average annual rainfall = 32.7 mm during November to February, sandy loam soil) in Hainan province, and one field trial was conducted with a rain shelter at the experimental research station of Nankou (N40°13′31.06″E116°05′38.86″, fluvo-aquic soil), Beijing, China. All accessions were tested under two conditions: well-watered (full irrigation; WW) and water-stressed (drought-stressed; WS). Plots were irrigated at sowing time to ensure maximum seed germination. Then, to simulate conditions that induce drought stress under WS condition, irrigation was ceased for plots in Dunhuang and Sanya, and the rain shelter was closed during precipitation in Nankou. Under WW conditions, plots were irrigated as needed to ensure normal plant growth. Field management, including fertilization, pest and disease control, was the same as that used for the general maintenance of each field site.

Table 1. Analysis of variance (ANOVA) of soybean seed yield per plant (YP) under drought stressed and full irrigation conditions for 149 soybean accessions.

Source^a	DF^b	Sum of squares	Mean square	F value	Pr > F
Block (Env × Treatment)	12	613.1578	51.0965	11.19	<0.0001
Treatment	1	10259.0812	10259.0812	2246.78	<0.0001
Env	2	19780.1371	9890.0686	2165.97	<0.0001
Env × Treatment	2	1031.5607	515.7803	112.96	<0.0001
Geno	148	5221.2699	35.2789	7.73	<0.0001
Geno × Treatment	148	1210.2212	8.1772	1.79	<0.0001
Geno × Env	296	4956.8638	16.7462	3.67	<0.0001
Geno × Env × Treatment	295	1795.3308	6.0859	1.33	0.0004

^aBlock (Env × Treatment) means the block is nested within Environment × Treatment; Env means environment; Geno means genotype; and Geno × Env means genotype by environment interaction.

^bDegrees of freedom.

At the Dunhuang experimental station, accessions were grown in two-row plots 2 m in length, with 0.50 m spacing between rows and 0.10 m spacing between plants in a randomized block design with 3 replications. At the Nankou experimental station, all accessions were grown in one-row plots 2 m long, with 0.50 m spacing between rows and 0.10 m spacing between plants, in a randomized block design with 3 replications. At the Sanya experimental station, all accessions were grown in one-row plots 5 m in length, with 0.35 m spacing between rows and 0.06 m spacing between plants in a randomized block design with 3 replications.

Five plants in the middle of each row were selected to determine seed yield per plant under full irrigation (YP-W) and drought-stressed (YP-S) conditions at maturity.

Phenotypic data analysis

The best linear unbiased prediction (BLUP) values of YP under both conditions were calculated using the MIXED procedure in SAS software v.9.1.3 [38] and were used to calculate the drought susceptibility index (DSI) using the following formula: $$\text{DSI}=\left(1–Y_{ws}/Y_{ww}\right)/\text{DII}$$ where Y_ws and Y_ww are BLUP values of YP under WS and WW conditions, respectively, and DII = 1–${\overline{Y}}_{ws}/{\overline{Y}}_{ww}$, where ${\overline{Y}}_{ws}$ and ${\overline{Y}}_{ww}$ are the averaged values across all genotypes under WS and WW conditions, respectively [9, 39, 40]. Thus, the higher the yield under WS conditions, the smaller the DSI score, indicating higher drought tolerance for accessions of interest. Hereafter, YP-W, YP-S and DSI were used to evaluate seed yield under full irrigation, seed yield under water deficit and drought stress index, respectively, in the following genetic study. Analysis of variance (ANOVA) for phenotypic data was conducted using the MIXED procedure of SAS software v.9.1.3 [38], with the accessions as the fixed effect and the replications as random effects. The correlation coefficients between YP-W, YP-S and DSI were calculated using the CORR procedure of SAS software v.9.1.3 [38].

Genotyping and genotypic data analyses

In general, the genotyping protocol and genotypic data analyses, including population structure analysis, principal component analysis (PCA), unrooted neighbor-joining tree, analysis of molecular variance and F-statistics analysis, followed the same pipeline as described in Liu et al. [41]. All of the accessions in our diversity panel were genotyped with the Illumina SoySNP6k iSelect BeadChip (Illumina, USA), which contained 5361 single nucleotide polymorphisms (SNPs) [42]; the chromosomal distributions and quality control for these SNPs were demonstrated in Wen et al. [43]. A subset of 4995 SNPs (with missing rates lower than 0.25 and minor allele frequency higher than 0.05) were selected from the total 5361 SNPs and used in the following analyses.

Genome-wide association study

For marker–trait associations, a multi-locus random-SNP-effect mixed linear model controlling both population structure and kinship matrices was used, and implemented by mrMLM.GUI 3.2 software [44] containing six GWAS models [44–52]. The kinship coefficient between each accession pair was estimated in TASSEL 5.0 using the Loiselle algorithm [53]. Markers were defined as significantly associated with the trait of interest by LOD ≥ 3 [44]. To avoid false positives due to the strong LD for the nearby SNPs, the LD blocks were determined by LD decay distance (r² = 0.1) and used as the support intervals to determine the QTLs significantly associated with a given target trait.

Results

Phenotypic variation and correlation analysis

ANOVA revealed significant variation among the different sites, as well as between blocks within each field trial, between genotypes, and significant genotype-by-environment (GE) interactions (Table 1). Significant variation emerged among replicate blocks within each of the three environments, which indicated that effects on the ANOVA linear model could be introduced by distribution across blocks (Supplemental Tables S2 and S3). The results also showed smaller effects of GE interaction in comparison with the effects of genotype. This finding suggested that phenotypic variation was likely dependent on variations in genotype.

Table 2. Descriptive statistics of seed yield per plant (YP) under drought stressed and full irrigation conditions and DSI-YP values for 149 soybean cultivars.

Traits	Minimum	Maximum	Mean	SD^d	CV^e
YP-W^a	4.9127	13.1329	7.2484	0.9676	13.35
YP-S^b	2.9344	6.4222	4.0045	0.6324	15.79
DSI^c	0.5745	1.2719	0.9964	0.1439	14.43

^aYP-W is seed yield per plant under full irrigation conditions.

^bYP-S is seed yield per plant under drought stressed conditions.

^cDSI is drought susceptibility index.

^dSD is standard deviation.

^eCV is coefficient of variation.

Table 3. Phenotypic BLUP (best linear unbiased prediction) correlations between seed yield per plant (YP) under drought stressed and full irrigation conditions and DSI-YP in the 149 soybean cultivars.

Seed yield per plant	YP-S	DSI^c
YP-W^a	0.6950***	0.1865*
YP-S^b		−0.5711***

*** Significance level <0.0001 in Pearson correlation test

*Significance level <0.05 in Pearson correlation test.

^aYP-W is soybean yield under full irrigation conditions.

^bYP-S is soybean yield under drought stressed conditions.

^cDSI is drought susceptibility index.

The average YP-W was 7.248 (range: 4.9127 − 13.1329; coefficient of variation, CV: 13.35), which was significantly higher than that of YP-S, that is, 4.0045 (range: 2.9344 − 6.4222; CV: 15.79) (Table 2; Figure 1). Consequently, the average DSI was 0.9964 (range: 0.5745 − 1.2719; CV: 14.43). Pearson’s correlation analysis among treatments and DSI scores revealed a positive correlation between YP-W and YP-S (r = 0.6950; p < 0.001), while DSI shared a positive relationship with YP-W (r of 0.1865; p < 0.05), but displayed a negative correlation with YP-S (r = −0.5711; p < 0.001) (Table 3). These results thus indicated that cultivars with high yield under full irrigation conditions were apparently more likely to exhibit a lower reduction in yield under water deficit conditions, therefore suggesting that selection of high-yield cultivars was necessary for higher yield during water limitation. Furthermore, the higher the yield for a given cultivar under water deficit, the smaller the yield reduction and the lower their DSI.

Figure 1. Frequency distributions of YP-W (A), YP-S (B) and DSI (C) in the 149 soybean accessions.

Based on the YP-S, YP-W and DSI, elite accessions were screened. Some accessions, including Jiunong 33 (ZDD24511), Changnong 8 (ZDD23773) and Jiyu 54 (ZDD23720), presenting higher YP-S and YP-W and smaller DSI, were selected. Those accessions have potential to be used as the founder parents in soybean breeding with high drought tolerance.

Genetic diversity of the 149 accessions based on 4955 SNPs

Genetic diversity analyses based on 4955 SNPs showed average values of 0.2310 for MAF, 0.3064 for gene diversity, 0.0737 for heterozygosity and 0.2450 for PIC (Supplemental Table S4, Supplemental Figure S1). In addition, LD analysis revealed mean r² = 0.2678 for the whole genome, which decreased to half that of the maximum at ∼2100 kb decay distance. By contrast, r² was < 0.1 at a decay distance of 10,500 kb (Supplemental Figure S2). Taken together, these findings showed that LD decay was relatively slow throughout the diversity panel population. Moreover, the SNP marker density was 188.77 kb across the full test population, indicating that this set of SNPs was sufficiently robust for accurate detection of major drought tolerance QTLs in this population.

Figure 2. Genetic structure and relatedness of the 149 soybean accessions. (A) Evolution of the natural logarithm probability of the data against K and the magnitude of Δk for each k value. (B) Clustering using k = 2 for the entire panel of soybean accessions. Each accession is represented by a vertical bar, partitioned into colored segments with the length of each segment representing the proportion of each individual genome contributed by either group when k = 2. (C) Unrooted neighbor-joining tree constructed using SNP data; soybean accessions from Heilongjiang are shown in red; from Jilin are shown in black; from Beijing are shown in purple; and from Japan are shown in blue. (D) Principal component analysis for the entire panel of soybean accessions.

Table 4. SNPs significantly associated with YP under drought stressed and full irrigation conditions and DSI-YP.

Marker	Chr.	Position	Associated Traits (R^2 %)	Method	−log10(P)	Reported QTLs/genes
ss244375521	1	5369218	YP-W(4.71/6.60)	FASTmrMLM/ pLARmEB	4.7878/ 6.2699
ss244871057	2	47979951	DSI-YP(1.70)/YP-S(2.33)	pLARmEB/ pLARmEB	4.0728/ 4.2563	Seed coat hardness [56]; Lateral root density [57]; Root dry weight [58]; Seed weight [59]; Plant dry weight [58]; Seed protein [60]; Seed oil [60]
ss245133420	3	33143761	YP-S(4.00)	pLARmEB	4.2432
ss245630456	5	8589397	YP-W(26.65)	FASTmrEMMA	11.5589
ss245842048	6	8979504	DSI-YP(12.02/7.68)	mrMLM/ FASTmrEMMA	4.8821/ 4.0254	Seed isoflavone [61]
ss245889259	6	15057987	YP-W(5.78)	FASTmrMLM	4.2239	Seed weight [62]; Seed set [63]; Seed protein [64]; Pod maturity [65]; Sudden death syndrome [66]
ss246509299	8	10755506	YP-S(2.18)	pLARmEB	3.997	Canopy wilt [34]; Seed coat hardness [56]; Root dry weight [67]; Root fresh weight [67]; Seed weight [62]; Seed set [63]; Seed oil [60]; Soybean cyst nematode [68]
ss246573652	8	17927603	YP-S(4.58)	pLARmEB	4.4939
ss246821668	9	6300834	YP-W(5.67)	pLARmEB	7.2029
ss247081352	10	2715719	YP-W(8.49/13.73)	FASTmrMLM/ pLARmEB	4.487/ 8.3561	Seed yield [69–70]; Seed set [71]; Root to shoot weight ratio [72]
ss247927074	13	9019809	YP-S(3.13)	pLARmEB	3.904	Phytophthora root rot [73]; Seed stearic acid [74]
ss248616287	15	12547291	YP-W(4.41)	pLARmEB	6.6676
ss249310678	17	9428953	YP-W(2.61)	pLARmEB	3.9854	Canopy wilt [33, 35, 36]; Drought index [9]; Seed weight [75]; Seed yield [76]; Seed oil [77]
ss249576244	18	7314370	YP-W(4.26)	FASTmrMLM	4.9919
ss249600231	18	9265696	YP-S(8.15/5.07)	mrMLM/ FASTmrMLM	4.0758/ 4.6635
ss249889306	18	55508285	YP-S(3.51)	ISIS EM-BLASSO	4.0714	Seed yield [9]; Seed weight [78]; Plant height [79]; Soybean cyst nematode [80]
ss250063822	19	8965290	YP-S(4.83)	FASTmrMLM	4.6180	Plant height [81]; Seed oil [60]
ss250170677	19	35475642	YP-S(10.05/5.72)	mrMLM/ ISIS EM-BLASSO	5.6023/ 4.1864	Plant height [81]; Seed oil [82]
ss250223327	19	40292810	YP-S(4.94)	pLARmEB	4.0462	Seed weight [77]; Seed protein [60, 83]; Seed oil [60]

We next performed STRUCTURE analysis, and the resulting Δk values indicated an appropriate subdivision based on k = 2, that is, two main clusters emerged from the total association panel population (Figure 2A, B). Accessions largely from Jilin province (71 out of 78) comprised one of the groups, while the other main group consisted mostly of accessions from Heilongjiang province (54 out of 71) (Supplemental Table S5). In order to better understand the evolutionary relationships among these accessions, we also reconstructed a neighbor-joining phylogeny based on shared allele frequency. The results were in agreement with that of the STRUCTURE analysis and further supported the division of the panel of 149 accessions into two clades associated with geographical origin (Figure 2C). Subsequent PCA analysis of population structure using the 4995 SNPs revealed an obvious separation of the accessions into two clusters (Figure 2D), which further supported the results of STRUCTURE and phylogenetic analyses. We also performed analysis of molecular variance [54] and F-statistics analysis [55] to identify differentiation within the panel, and the resulting significant pairwise F_ST value of 0.17 (p < 0.001, Supplemental Table S6) indicated that the 149 accessions could be separated into two distinct sub-populations. Collectively, these results indicated that GWAS should use a Q matrix generated with k = 2 to screen for loci associated with seed yield under two treatments and drought susceptibility index.

QTLs for the drought tolerance traits

GWAS revealed that 19 QTLs were identified and distributed across 13 chromosomes (Figure 3), each of which explained 6.51% of the phenotypic variation on average (Table 4). Two QTLs were found on each of chromosomes 6 and 8, while chromosomes 18 and 19 each harbored three QTLs. Among the 19 QTLs detected by GWAS, 8 were associated with YP-W, 10 were associated with YP-S and 2 were associated with DSI, one of which was associated with both YP-S and DSI. It warrants mention that 11 of the QTLs reported here were located within 500 kb of QTLs identified using segregating populations in a previous study (Table 4). Among these previously reported QTLs, two were significantly associated with canopy-wilt, and one of these was also associated with drought index. Two of the other QTLs reported in this region were associated with seed coat hardness and root-related traits (Table 4).

Figure 3. Manhattan plots from the Q + K model across YP-W (A), YP-S (B) and DSI-YP (C). The names of candidate genes that are hit by significant SNPs are labeled on the top.

Candidate gene identification

We also looked for genes in close proximity or directly harboring these associated QTLs detected in this study (Supplemental Table S7) and found that eight QTLs were located within eight genes. In addition, two were found in exons, five within introns, and one was located in a 3’UTR. Specifically, Glyma.06g111000, encoding a pyridoxal phosphate (PLP)-dependent transferase superfamily protein, contained ss245842048 and was expressed with high level in leaf tissue, while Glyma.18g092200, encoding a WRKY family transcription factor, contained ss249600231. In addition, Glyma.18g092200, annotated as GmWRKY168, is reportedly responsive to cadmium [84], but exhibits no response to high salt concentrations [85]. Another gene, Glyma.03g120000, encoding alpha-l-arabinofuranosidase and highly expressed in roots, which contains a SNP in its 3′UTR, was reported to be down-regulated gene by shading [86]. We also found eleven QTLs located directly upstream or downstream of genes, among which four appeared likely to be located in the promoter regions, including ss250063822, ss246509299, ss245889259 and ss246573652. However, the function of these genes remains uncharacterized in soybean, and future work will clarify whether and how they could affect drought tolerance and yield in soybean.

Discussion

Decreased yield due to drought stress poses a long-standing problem for soybean breeders [87], given the role of drought as the most detrimental environmental factor limiting yield in cropping systems dependent on rain for irrigation. Currently, little is known about the underlying genetic basis of drought tolerance in soybean. While our results showed that YP-W was significantly higher than YP-S across accessions (Tables 1 and 2), the results of ANOVA revealed significant variation between replicates associated with field site (environment) and strong genotype-by-environment interactions for both YP-W and YP-S treatments (Supplemental Tables S2 and S3). These confounding factors thus increase the difficulty of determining effects due strictly to water deficit and can complicate field selection of tolerant accessions. To circumvent this issue in this work, we used the BLUP of genotypic value for each accession under each treatment as the phenotype.

In our study, YP-S and YP-W were significantly correlated (Table 3). The positive correlation between seed yield in water limited versus fully irrigated environments is in agreement with similar findings by Ramirez-Vallejo and Kelly [88] and Du et al. [9, 19]. Therefore, cultivars with the highest yield under water deficit conditions are most likely to be selected from among cultivars with the highest yield under full irrigation. Previous studies have proposed that cultivars with robust drought tolerance should exhibit stable yields, meaning that they have high yields in both normal and drought years, or at least display a limited reduction in yield under water deficit [89]. Ideal varieties can be bred through selection under both water-limited and full water environmental conditions, and the combination of mean seed yield in both drought stress-inducing and non-stress environments is recommended as the most effective criterion for selecting drought-resistant genotypes [39].

Seminal work by Fischer and Maurer [90] showed that the effects of variation on plant phenological characteristics and yield potential in tolerant genotypes can be minimized through the use of a yield-based drought susceptibility index (DSI) calculated by the relative loss in yield between drought and non-drought treatments within individual genotypes. Thus, the DSI may be the same for both high and low yielding varieties, dependent on their proportional reduction in yield under water restriction. The strategy of combining yield under both conditions, to calculate DSI, was further supported by Du et al. [9] for screening elite drought-tolerant lines.

In this study, we found eight QTLs that were significantly associated with YP-W, 10 associated with YP-S, but none that were significantly associated with seed yield under both treatments. Previous studies have proposed that instability in yields between drought-stressed and full irrigation conditions could be at least partially explained by regulatory contributions of several, unlinked (i.e., distributed across several chromosomes), small-effect loci [91]. Similarly, other early studies proposed that environmental variability can introduce GE interactions that lead to differential gene expression during drought stress compared to their expression profiles under full irrigation [9].

Since environmental conditions influence QTL stability, which can confound marker-assisted selection, QTLs that are robust to environmental variability may be therefore more informative than those detectable only under specific conditions. Thus, the QTLs identified here were compared with QTLs reported within a 500 kb vicinity using the Soybase (http://www.soybase.org) web-based toolkit. This analysis showed that 11 of the 19 QTLs we found were mapped to genomic regions previously associated with drought tolerance, seed yield and other traits. In addition, the overlap of QTLs associated with different traits could also be due to common physiological features shared between these traits or even a directly causal relationship between them [92]. However, these hypothetical interactions could potentially reduce the likelihood of finding false positive QTLs in overlapping regions. While this possibility is interesting, these overlapping QTLs require further detailed characterization and validation. The other eight detected QTLs that we identified in this study have not been reported in previous studies.

Notably, we found that ss244871057 (Chr.2) associated with both DSI-YP and YP-S, was located in the same region as QTLs responsible for seed coat hardness [56], lateral root density [57] and root dry weight [58]. Root characteristics such as deeper taproots or fibrous lateral roots are reportedly related to soybean drought tolerance [93–94]. In addition, drought tolerance in this crop is also purportedly associated with seed coat hardness, and cultivars with strong drought tolerance frequently also exhibit high percentage of seed hardness [95]. Our analysis also showed that the region containing the QTL for SNP marker ss246509299 (Chr.8) correlated with YP-S is 121,598 bp distant from Satt424, detected by Abdel-Haleem et al. [34], which is associated with delayed canopy wilting. Interestingly, one QTL related with seed coat hardness [56] and a QTL related with root weight [67] were also identified in the same region as ss246509299. It also deserves mention that the region containing the QTL for SNP marker ss249310678 (Chr.17), significantly associated with YP-W, is 147,691 bp away from Satt154, reported by both Charlson et al. [33] and Hwang et al. [35] as also significantly associated with delayed canopy wilting. This satt154 marker was also described by Hwang et al. [36] as significantly related to canopy wilting, and by Du et al. [9] as associated with DSI. Collectively, these results indicate that these QTLs may be valuable for improvement of drought tolerance in elite, high-yield lines.

This work also shows that seven QTLs, distributed across Chr. 2, 6, 8, 10, 17, 18 and 19, are genetically linked to markers associated with seed weight- and yield-related traits (Table 4). In particular, the YP-S-associated SNP marker ss250223327 (Chr. 19) is located 344,261 bp away from Satt156, which is reportedly correlated with seed weight [77] and seed protein content [83]. This SNP is also 27,174 bp from Satt481, reportedly associated with seed protein and seed oil [60]. These findings suggest that these QTL regions could also be used in the concurrent improvement of tolerance and yield.

Conclusions

In this study, 19 SNP–trait associations were identified by GWAS. Among them two SNPS were significantly associated with DSI, and 11 SNPs were located in or close to agronomically important QTL mapped by linkage analysis. Eight candidate genes were identified and Glyma.18g092200 annotated as GmWRKY168. The important drought-tolerant SNPs and potential genes in this study, will be useful and feasible to accelerate genetic gain in soybean breeding.

Authors’ contributions

Zhangxiong Liu, Lijuan Qiu and Xusheng Qi: Conceptualization and Writing – revision and editing. Zhangxiong Liu and Lijuan Qiu: Resources. Zhangxiong Liu, Xingrong Wang, Yanjun Zhang and Zuowang Gou: Investigation. Huihui Li, Honglei Ren, Zixiang Wen, Yinghui Li, Lili Yu, Huawei Gao, and Dechun Wang: Data curation. Zhangxiong Liu, Huihui Li, Xingrong Wang, Yanjun Zhang and Xingzhen Zhao: Writing - original draft.

Disclosure statement

The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability statement

The data that support the findings of this study are available from the corresponding author, [LJQ], upon reasonable request.

Additional information

Funding

This work was supported by the National Key R & D Program for Crop Breeding (Grant 2021YFD1201603), the Platform of National Crop Germplasm Resources of China (Grant NCGRC-2020-05), the Agricultural Science and Technology Innovation Program of the Chinese Academy of Agricultural Sciences (CAAS), and the Improvement of Soybean Abiotic Stress Tolerance to Address the Climate Change (Grant PJ0121092018).