Introduction

Soybean (Glycine max (L.) Merr.) is among the world's five most important crops, providing food, feed, protein, and edible oil. However, it is subjected to several abiotic stresses like drought, salinity, and heat stress1. Among these stresses, drought stress is going to be a major threat to crop plants. In the current scenario, increased temperature, extreme variability of seasonal rainfall, and prolonged dry season due to climate change have contributed to the water scarcity challenge2. Water scarcity is becoming a major constraint for crop production throughout the world3. In soybean, drought stress mostly occurs during the whole growing season which leads to considerable yield loss and quality deterioration mainly in arid and semi-arid zones. Some studies have reported that drought stress could cause 24–50% yield losses in soybean4,5 The research studies have also shown that during juvenile and reproductive stages soybean is sensitive to drought stress6,7 In most of the soybean growing regions, irrigation is not a possible option due to lack of water availability and the high cost of irrigation. Therefore, the most economical way of sustainable soybean cultivation is to develop soybean genotypes having the ability to germinate and survived under drought stress condition8.

The most apparent symptoms of drought stress in soybean are reduction in germination percentage, cell elongation and expansion, hypocotyl length, shoot and root fresh weight, increased root length, leaf flipping, leaf clamping and initiation of early flowering5. Drought stress also disturbs carbon assimilation and plant phenology9. The germination and seedling stages are the most sensitive to drought stress in soybean and other crops like rice10, maize11 and wheat12. To combat with the challenge of drought stress, plants have adopted various morphological, physiological and molecular mechanisms including the robust root system, reduced leaf area and stomatal conductance, maintenance of turgor pressure, early maturation, activation of drought responsive genes and transcriptional factors e.g. DREBs, WRKY, ZIP, ZFP, NAC, MYB and ERFs13,14,15,16. It is predicted that in future, plants will face drought stress more frequently due to high temperature extremes17. Therefore, efforts are necessary to explore the genetic base and genetic mechanism regulating the drought tolerance in soybean to combat climate change.

Drought tolerance is a polygenic trait influenced by genotypic and environmental factors18. Due to complex genetic architecture of drought tolerance, less genetic diversity, segregation distortion, linkage drag, and distortion and suppression of recombination conventional breeding fails to meet the requirement for stabilizing yield potential in plants under variable conditions19,20,21,22. Advances in genetic and genomic approaches have assisted in the dissection of genetic control of various complex traits which have facilitated mapping of numerous quantitative trait loci (QTLs) in different crops, including soybean23,24. However, QTL mapping provides inadequate genetic diversity as reported in Oryza sativa25. To far, numerous genetic maps have been created to investigate QTLs that affect drought tolerance in soybeans. But in practical breeding, confirmation and verifications of QTLs are recommended prior to their application in the breeding program26.

Advancement in next-generation sequencing (NGS) has facilitated exploiting marker assisted selection (MAS) to enhance the genetic gains in an economical way. Genome wide association studies (GWAS) by utilizing genome wide single nucleotide polymorphism (SNPs) have the potential to detect variants that could be used directly for crop plant improvement23,27. GWAS could overwhelm the constraints of traditional QTL mapping and can be utilized to investigate the minor gene effects controlling complex traits in crop plants28. GWAS has been successfully used in breeding of various crops like cotton, wheat, barely, rice and corn13,23,29,30.

Association mapping accuracy relies on marker density and mapping populations31. Among all the approaches utilized for GWAS, the mixed linear model (MLM) is applied widely for association mapping32. Furthermore, advanced MLM- based methods have also been established like compressed MLM (cMLM), which saves time for computation and enhances the statistical performance. Bonferroni correction is indeed for such tests since these approaches rely on single-locus genome-wide scanning32. Currently, a more refined method is purposed known as random-SNP impact MLM (mrMLM), which combines the advantages of both multiple locus model and single locus genome scanning. The mrMLM revealed significantly better accuracy and power in the detection of QTN and affected the estimation24.

Several GWAS have been conducted and reported in soybean33,34. Still, literature related to GWAS for drought tolerance at germination is limited35. Germination is a critical step contributing toward final soybean production. Cultivars with a better germination percentage during drought stress produced higher yields36. The information is scarce about molecular genetics controls drought tolerance during the germination stage. As a result, it's critical to identify the QTLs linked to drought tolerance in soybeans during the germination stage.

The present research includes marker-trait association depending upon phenotypic and genotypic information from 240 worldwide germplasm accessions with wide diversity. The 240 accessions were screened for drought tolerance at germination stage with popular method for phenotyping. We used both single and multi-locus models viz., MLM and mrMLM to detect the significant SNPs/QTNs, as well as for the prediction of candidate genes linked to drought tolerance during germination in soybean. Physiological mechanisms associated with drought tolerance at the germination stage may be akin to those that help plants to survive under drought stress during seedling and reproductive stages. The objective of present study is to unravel the genetic constitution of drought tolerance in 240 soybean accessions at the germination stage using GWAS.

Results

Phenotypic characterization of the 240 germplasm

The drought tolerance of 240 soybean accessions was assessed on the base of three germination linked traits i.e. RL, GR and WSL. The values of mean, standard deviation, range, ANOVA (F-value) and estimated h2 for all traits of 240 lines in two different years 2018 and 2019 are shown in Table 2-1. During 2018, the average value of GR, WSL and RL were 0.30, 0.32 and 0.43 and their range varied from 0–1.25, 0.11–1.75 to 0.0028–0.96 respectively. However, the mean comparison revealed that, during 2019, the average values of GR and WSL were higher than the mean values in 2018 whereas RL growth pattern was found contradictory to it (higher in 2018). These findings depicted that, in soybean, drought tolerance is influenced by environment. The results of ANOVA indicated highly significant difference between genotypes (G) and environments (E). Furthermore, significant genotype× environment (G×E) interaction was detected for all three studied traits (Table 1). Furthermore, h2 for GR, WSL and RL across both environments were high in the range of 0.76–0.99, showing that genetic factors play a vital role in drought tolerance as compared to environmental factors.

Table 1 Descriptive statistics, heritability (H2) and F values from ANOVA of three drought tolerance indices in 2 years (2018 and 2019).
Full size table

SNPs distribution, population structure and LD estimation

In this study, out of a total 42,509 SNPs, 17.98% SNPs having MAF ≤ 0.05 were extracted for analyses, while the leftover 34,817 SNPs with MAF > 0.05 were used to investigate the genetic variability in 240 soybean genotypes and other subsequent analyses. These targeted SNPs were randomly located on 20 chromosomes of G. max (Fig. 1A). On an average, each chromosome (Chr) carried 1762 SNPs. Maximum number of SNPs (2922) were mapped on Chr.18 while minimum number of SNPs (1223) were found on Chr. 12. Furthermore, the SNPs distribution density on 20 chromosomes is presented in Fig. 1B that varied from 25.87 SNPs/Mb on Chr.01 to 51.01 SNPs/Mb on Chr.13.

Fig. 1
figure 1

Distribution of 34,817 SNPs across 240 soybean PIs. (A) Number of SNPs per chromosome (B) SNP density distribution (SNPs/Mb).

Full size image

Besides this, three different methods, principal component analysis (PCA), phylogenetic analysis, and bayesian model-based methods were employed to build the population structure of 240 soybean accessions. STRUCTURE software was used to compute Delta K (∆K) value (Fig. 2A; K = 1–10, Fig. 2B), and the result revealed that 240 soybean PIs were grouped into three subpopulations i.e. K = 3. Moreover, PCA analysis and phylogenetic tree showed consistent results which agree with the findings of population structure analysis (Fig. 2C, D).

Fig. 2
figure 2

An examination of the population structure of 240 soybean lines. (A) ΔK plot calculated for K = 1–10 (B) structure software was used to determine the population structure. Three subpopulations are represented by three colors (red, green, and blue) (C) the tested 240 lines’ maximum likelihood neighbor-joining tree (D) principal component analysis plot.

Full size image

Furthermore, the r2 values were calculated for those SNPs pairs that resided within 10 Mb region and the LD decay rate was established based on the regression of negative natural logarithm. In G. max, LD decay was comparatively higher37 than other crop plants like A. thaliana38, O. sativa39, Z. mays40. In addition, our results showed that r2 values decreased gradually as the distance increases and At 1.76 Mb, the LD decay distance was measured, and r2 dropped to half of its highest value (Fig. 3).

Fig. 3
figure 3

Linkage disequilibrium (LD) decay across the genome in a population of 240 soybean lines. The squared correlation coefficient was used to get the mean LD decay rates.

Full size image

Association mapping for traits via MLM

The MLM model was utilised in GWAS research to identify important QTNs related to GR, WSL, and RL. Results highlighted that a total of 4, 9 and 10 QTNs were found to be associated with GR, WSL, and RL, respectively, at a threshold value of − log10P = 4.5 (Table 2 and Fig. 4, Fig. 5, Fig. 6 respectively). Four QTNs associated with GR were mapped on two chromosomes only: Chr.03 and Chr.08. Among these QTNs, two QTN qGR3-1 and qGR8-1 were distributed on Chr.3 and 8 respectively and covered 8.55 to 10.59% phenotypic variation (PV) and they were identified consistently in all three environments (i.e. 2018, 2019 and combined environment) (Table 2, Fig. 4). These QTNs showed significant association with SNP marker Gm03-4741942, Gm08-12368145 and exhibited largest phenotypic variation effect (10.59%) in combined environment. Besides this, qGR3-2 and qGR8-2 were detected in two environments only and they covered 8.99–9.14% of the PV (Table 2).

Table 2 Quantitative trait nucleotides (QTNs) associated with germination rate (GR), whole seedling length (WSL) and root length (RL) via MLM.
Full size table
Fig. 4
figure 4

Manhattan plots and quantile–quantile (QQ) plots were used in the population-based genome-wide association studies for (A) GR-2018. (B) GR-2019. (C) GR-Combined.

Full size image
Fig. 5
figure 5

Manhattan plots and quantile–quantile (QQ) plots were used in the population-based genome-wide association studies for (A) WSL-2018. (B) WSL-2019. (C) WSL-Combined.

Full size image
Fig. 6
figure 6

Manhattan plots and quantile–quantile (QQ) plots were used in the population-based genome-wide association studies for (A) RL-2018. (B) RL-2019. (C) RL-Combined.

Full size image

For WSL, a total of 9 QTNs were identified on six different chromosomes i.e. Chr.08, 09, 12, 13, 16 and 18 (Table 2, Fig. 5). Among these, qWSL-8 and qWSL13-1 were found in all the three environments and categorized as stable QTNs linked with WSL. These two QTNs accounted for 8.01%–8.63% phenotypic variations, respectively. Two QTNs linked to WSL viz., qWSL13-2 and qWSL13-3 were identified in two environments and explicated 8.84–9.77% of the PV respectively. However, qWSL-16–1 and qWSL-9–1 were detected during 2018 and 2019, respectively (Table 2).

We used the MLM approach to identify 10 QTNs in RL, and all these QTNs were located on 5 distinct chromosomes viz., Chr.06, 08, 15, 14 and 16 (Table 2). Out of these QTNs, only qRL-8 and qRL15-4 were continually detected in all environments (2018, 2019 and Combined), elucidating 8.12–9.83% of the PV and these QTNs were linked SNP marker Gm08-12368145 on Chr.08 and Gm15-4497855 on Chr.15, respectively. Other QTNs, qRL16-1 and qRL-14 were detected in 2018 and 2019 environments, explaining 8.14% and 8.49% of the phenotypic variance, respectively. The three remaining QTNs associated with RL were only found in one environment and accounted for 8.08 percent–8.12 percent of the PV (Table 2 and Fig. 6).

Association mapping for traits via mrMLM

Multi-locus GWAS analysis was performed by a mrMLM approach to ensure the QTNs dependability detected by the MLM model and to identify further QTNs related with GR, WSL and RL. As shown in Table 2-3, a total of 8, 8 and 11 QTNs were linked with GR, WSL and RL respectively (Table 3). For GR, one QTN qGR8-1 linked with the Gm08-12368145 marker was constantly found in all environments, and this QTN showed broad PV effects ranging from 8.40 to 9.95%. The remaining seven GR QTNs were only found in one environment, accounting for 4.04 to 9.95% of the PV.

Table 3 Quantitative trait nucleotides (QTNs) associated with germination rate (GR), whole seedling length (WSL) and root length (RL) via mrMLM.
Full size table

In case of WSL, 8 QTNs were identified, scattered on 6 different chromosomes. Among them, only single QTN qWSL-13 was mapped in all three environments with larger effect, while qWSL-18 was identified in two environments (2019 and combined). However, the rest of 6 QTNs were identified in single environment only. For RL, a total of eleven QTNs were detected, explaining 3.7–18.6% of the PV. Out of these, qRL-8-1 associated with RL in all environments (2018, 2019 and combined). Two QTNs qRL7-1 and qRL7-2 were observed in two environments (2019 and combined). However, the remaining eight QTNs were environment specific QTNs (Table 3). In short, a total of 23 and 27 QTNs linked with all studied traits of drought tolerance were observed via MLM and mrMLM model, respectively (Tables 2 and 3).

Comparison between MLM and mrMLM results

In comparison with the MLM and mrMLM, For MLM, six QTNs, viz., qGR3-1, qGR8-1, qWSL8, qWSL13-1, qRL-8 and qRL15-4 were detected in all three environments. As predicted, for mrMLM, three major QTN qGR8-1, qWSL13-1 and qRL-8 were linked with GR, WSL and RL were detected. These three critical QTNs were detected by both MLM and mrMLM. Nevertheless, the qGR3-1 and qWSL8-1 detected by MLM could not be confirmed by mrMLM. Furthermore, qRL6 and qRL14 were also detected by both MLM and mrMLM in 2018 and 2019, respectively (Table 3). Based on above mentioned results, two major and stable QTN i.e. qGR8-1(Gm08-12368145), qWSL13-1(Gm13-35394159) and qRL8 (Gm08-12368145) were identified in all studied environments via MLM and mrMLM methods. The QTNs, qRL8 and qGR8-1 were renamed as QTN8. QTN qRL15-4 (Gm15-4497855) was identified in all three environments in MLM but not in mrMLM. Thus, these major QTN on Chr.8 and 13 were potentially employed as the hotspot area for fine mapping and analysis of candidate genes linked with drought tolerance in soybean.

Candidate gene prediction

As the stable QTN found on Chr.08 (QTN8) and Chr.13 (qWSL13-1) were constantly detected in all studied traits i.e. GR, WSL and RL by both methods in all environments. For further confirmation, candidate gene prediction analysis has been conducted in the genomic region adjacent to the QTN8-1 and QTN13-1. Based on LD distance, 500 kb up as well as downstream region carrying peak SNPs markers was explored and within the interval of 1.0 Mb, a total of 100 and 139 genes were mapped on Chr.08 and Chr.13, respectively. The haplotype block analysis was carried out with Haploview 4.2 for all the SNPs located within ± 100 kb candidate region (Fig. 7). Based on gene annotation data derived from Soybase (http://www.soybase.org) and the Phytozome (https://phytozome.jgi.doe.gov), 15 genes were identified as possible candidate genes regulating the drought tolerance in soybean (Table2-4). Out of 15 candidate genes, Glyma08g156800, Glyma08g160000, Glyma08g162700 and Glyma13g249600 only showed higher expression in soybean root tissues.

Fig. 7
figure 7

Candidate regions and haplotype blocks of the genome significantly linked to identified major QTNs for GR and RL on Chr.08 and for WSL on Chr.13 in 240 soybean PIs.

Full size image
Table 4 In two distinct environments, phenotypic correlations (r) between GR, WSL, and RL.
Full size table

Discussion

In this study, significant diversity has been observed in 240 soybean PIs regarding the GR, WSL and RL values that are three important drought-linked parameters at germination stage. Results depicted that GWAS analysis using 34,817 SNPs detected QTLs related with drought tolerance (Table 2, Table 3). Furthermore, strong positive correlations of r = 0.79 was observed between GR and RL in the 240 soybean PIs (Table 4). Drought stress impairs germination by limiting water imbibition41. Therefore, for evaluation of soybean genotypes, we screened the population at germination stage by measuring GR, WSL and RL. On the basis of correlation analysis, it was cleared that WSL showed a positive correlation with GR and RL (Table 4).

Two stable QTNs were discovered among the QTNs identified, viz. QTN8, (qGR8-1, qRL-8) and qWSL13-1 were detected in all environments (2018, 2019, and Combined) via MLM, explaining 8.34–10.59% of the PV, whereas these QTNs (QTN8 and qWSL13-1) were also detected by mrMLM. Remarkably, Stable QTN (QTN8) located on Chr.08, was linked with SNP marker Gm08-12368145 for two traits (qGR8-1 and qRL-8), this QTN was discovered in all of the environments examined. Furthermore, mrMLM was used to effectively identify this stable QTN8 on Chr. 8. Another major QTN (qWSL13-1) associated with Gm13-35394159 was detected in all environments (Table 3). Thus, this stable QTN8 detected on Chr.08 in present study is a newly QTN linked with drought-tolerance in soybean.

However, the second stable QTNs viz., qWSL13-1 identified via MLM was detected in all environments, was only linked to a single trait. Therefore, on the base of above mentioned results, it is summarized that in this research all the identified stable QTNs were novel for drought tolerance in soybean and mostly carried larger phenotypic effects. Besides this, QTN8 and qWSL13-1 exerted larger phenotypic effects on drought tolerance. That’s why, we carried out candidate gene prediction within 1.0 Mb region adjacent to the QTN8/Gm08-12368145 and Gm13-35394159. Based on functional annotations, fifteen candidate genes which are directly or indirectly linked with drought tolerance were selected (Table 5).

Table 5 Functional annotation of candidate genes.
Full size table

Up till now, various studies have been carried out in soybean to explore the genetic mechanism that regulates drought tolerance. Several QTLs linked with drought tolerance have been identified across different environments42,43,44. However, most of the research related to QTLs mapping was carried out by employing linkage mapping analysis on the basis of low-density genetic maps. Drought tolerance is a complex character as controlled by several minor loci with significant environmental influence. To overcome the linkage mapping limitations, LD based GWAS analysis was conducted in this research. Several studies have confirmed that utilization of GWAS analysis could be an alternative and powerful approach for precise detection and defining the genomic regions associated with quantitative traits. Utilization of GWAS is an efficient alternative approach for precise detection and marking the genomic region of loci/genes linked with complex traits16,45.

Linkage disequilibrium refers to non-random association of SNP alleles that are located at different locus. Nonetheless, the amount of genome-wide LD is critical in GWAS research. As it is influenced by several forces viz., population stratification, linkage, genetic drift, and relatedness13,46. For the soybean panel, the LD decay was approximately assessed at 1.76 Mb, which was markedly higher as compared to other plant species i.e. rice (75 kb in indica)39, Arabidopsis thaliana (10 kb)39 and 1.1 kb in sunflower47. In the self-pollinated species like soybean comparatively higher LD value has been recorded as compared to cross-pollinated species48,49. Furthermore, resequencing of 31 soybean genomes, Lam, Xu, Liu, Chen, Yang, Wong, Li, He, Qin and Wang37 stated that both cultivated and wild soybean species showed high LD, while the comparison of both species showed that cultivated soybean demonstrated higher LD level than wild species. Therefore, for effective improvement of the genetic variability of GWAS population genotypes should have diverse genetic background (including wild species, cultivars, and landraces) collected from different geographical regions.

From previous functionally verified genes, transgenic lines with GmNAC8 overexpression show noticeably greater drought tolerance than those with GmNAC8 deletion, indicating that GmNAC8 is a positive regulator of drought tolerance50. GmWRKY54 enhanced drought tolerance by activating genes in the abscisic acid and Ca2+ signaling pathways from functionally confirmed genes. In Arabidopsis, overexpression of WRKY20 from a wild soybean greatly increased drought tolerance51,52. Using a constitutive expression promoter, it was possible to overexpress the LOS5/ABA3 gene in soybeans, encoding the molybdenum cofactor sulfurase, therefore improving drought tolerance and increasing seed output by at least 21% during drought stress53. In plants, U-box (PUB) proteins serve as E3 ligases. Under drought stress conditions, the transgenic lines of Arabidopsis showed reduced plant viability when GmPUB6 was overexpressed54.

Manavalan, Prince, Musket, Chaky, Deshmukh, Vuong, Song, Cregan, Nelson and Shannon43 has reported one QTN associated with drought tolerance linked with marker Satt315 on Chr.08 within a confidence interval of 15 cm. Depending upon the genomic sequence information, it was found that this QTL was linked with Satt315 from SoyBase (http://www.soybase.org), this QTL was present between 6751690 and 6751725 Mb of physical interval on Chr.08. In another study, Abdul Haleem et al. (2011) identified one QTL Q_root_Gm08 (flanked with Satt228 and Satt429 on Gm08) via MLM relevant to the fibrous root system. This QTL was located between 46366169 and 46366092 Mb of physical interval of Chr.08. (QTN8)55 also reported a cluster of 4 QTL associated with seed and root related traits (RFW, RDW, SFW, and SDW) on LG A2 (chromosome 8). This QTL is linked with marker Satt177 within a confidence interval of 15 cM. Charlson et al. (2009) also identified QTL (Sat_319) on MLG A2 contributed slow-wilting. This QTL was located between 4717934 and 4717907 Mb of physical interval of Chr.08. Wang et al., (2012) also identified QTL (I-Satt327 on G2-A2) controlling drought tolerant traits during germination stage. This QTL was located between 20511527 and 20511574 Mb of physical interval of Chr.08 but the physical position of Gm08-12368145 was 12.36 Mb of Chr.08.

Materials and methods

Plant materials and field evaluation

An association mapping population comprised 240 worldwide (Table S2-1) soybean plant introductions (PIs). The collections were obtained from 22 countries across Europe (EU), West Asia and North Africa, East Asia, and the Americas. The list of germplasm accessions and their origins are presented in Supplementary Table 1. Seeds of the 240 germplasm were collected from National Center for Soybean Improvement, Nanjing Agricultural University, Nanjing–China. The seeds obtained were sown in Jiangpu Experimental station of Nanjing Agricultural University in 2018 and 2019. The healthy and uniform seeds were selected for drought tolerance screening in the laboratory.

Germination test

To test the germination percentage, twenty seeds per entry were randomly selected for each of three replications. For evaluation of growth performance and phenotypic diversity related to seedling traits, seeds were kept on wetted filter paper in a 9 cm diameter petri dish. Drought stress was created by the addition of Polyethylene glycol 1000 (PEG-1000, 15% w/v) whereas distilled water (without PEG) was employed in control conditions. The Petri dishes were retained in an incubator at 25 ℃ in the darkness. The experiment was carried out in complete randomized design (CRD). The seeds were categorized as germinated when radicle attained the length of 2 mm. Germination was recorded at 24 h’ intervals for five consecutive days. To evaluate germination traits, we measured germination rate (GR, whole seedling length (WSL) and root length (RL) under drought and control conditions. Germination related traits were scored as per International Seed Testing Association rules (ISTA) (ISTA, 2003):

$${\text{Germination percentage }} = {\text{n}}/{\text{N }} \times { 1}00$$

where n and N represent number of seeds germinated and the total number of seeds sown, respectively.

For measurement of whole seedling length (WSL) and root length (RL), three plants of each accession were taken under control conditions and 15% PEG stress. After five days, RL and WSL (in cm) were measured manually by scaled ruler56. The data for these traits was utilised to create indices, which were calculated as the ratio of GR, WSL, and RL, germination traits, at drought conditions to the control conditions of same trait. The germination rate (GR) was scored by using the following equations:

$${\text{GR}} = \frac{{\text{Germination rate under drought condition}}}{{\text{Germination rate under normal condition}}}$$

Data analyses

Using SAS PROC UNIVARIATE tools, descriptive statistics such as mean, standard deviation (SD), and range were calculated for 240 soybean PIs in two environments: 2018 and 2019 respectively. The PROC GLM and PROC CORR programmes of SAS were used to calculate analysis of variance (ANOVA) and Pearson’s correlation for drought traits57. The broad sense heritability of studied traits under drought and normal condition was calculated for combined environment using the below formula

$${H}^{2}=\frac{{\upsigma }^{2}g}{({\upsigma }^{2}g+\frac{{\upsigma }^{2}ge}{n}+\frac{{\upsigma }^{2}e}{nr})}$$

For individual environment heritability was calculated by using the following formula:

$${H}^{2}=\frac{{\upsigma }^{2}g}{{(\upsigma }^{2}g+{\upsigma }^{2}e)}$$

where σ2g is the genotypic variance, σ2ge is interaction variance between genotypes and environment, σ2e is the error variance, n is the number of environments, r and e represent replications and error respectively58.

Genotyping and quality control

The SNP data for the 240 accessions is from SoySNP50K BeadChip available on SoyBase database (soybase.org/snps/index.php) with 42,449 SNPs59. Briefly, a high-throughput SNP assay, the SoySNP50K iSelect BeadChip containing 42,449 SNPs was utilized for genotyping following procedure reported by60. Chromosomal distributions of SNPs and quality control were estimated following the procedure of Wen et al. (2014). A subset of 34,817 SNPs (MAF > 0.05 and missing rates < 0.25) was selected out of 42,449 SNPs for further analyses.

Principal component, linkage disequilibrium, and population structure analysis

After screening using the PLINK in-depth-pairwise command 7358 SNPs were retrieved for this study61. Using the software STRUCTURE 2.2, a model-based cluster assay was performed to find the genetic structure and define the cluster (gene pools) in the data62. Potential populations (K), the number was from 1 to 10 with four repetitions, and each Q and P-value was estimated. The Delta K method was used to calculate the most likely number of subpopulations63. TASSEL 5.0 software was used to create a neighbor-joining tree (NJT)64. R studio was used to perform principal component analysis (PCA). TASSEL 5.0 was used to estimate Kinship, and further pairwise SNPs linkage disequilibrium (LD) was estimated as a squared correlation coefficient (r2) using RTM-GWASV1.165. The distance of chromosome at which r2 dropped to half of its highest value was used to determine the linkage disequilibrium decay rate.

Marker-trait association and haplotype block analysis (GWAS)

After eliminating SNPs with an MAF < 5%, 34,817 SNPs were remained for the 240 PIs. To reduce false-positives and enhance statistical performance, the principal component (PC) and kinship (K) matrix were calculated. For the SLM (Single locus model) viz., MLM (PCA + K) was carried out in TASSEL5.2 (Bradbury et al., 2007) while for multi-locus model, mrMLM was conducted in R using mrMLMGUI package (version 2.1)66. A threshold level (–log10 (P) value = 4.5; P = 1/34,817) was set to calculate a significant linkage of SNPs with drought tolerance parameters in MLM whereas a critical LOD value of 3 was adopted for mrMLM (Xu et al., 2018; Misra et al., 2018). The SNPs identified by both models and in all environments was regarded as comparatively stable SNPs. The QTNs were named after the nomenclature of McCouch et al. (1997) with minor changes. Furthermore, the abbreviation was started with ‘q’, followed by trait name (GR, germination rate; WSL, whole seedling length; RL, root length) and the chromosome name, succeeded by number of QTN identified on same chromosome. Thus, qGR8-1 designated a locus on chromosome 08 related with GR in all environments. Haplotype block analysis of the comparatively stable SNPs was carried out in Haploview software version 4.2 with default parameters (Barrett et al.37.

Prediction of candidate genes

Possible candidate genes were obtained within major QTN region found on Chr.08 and Chr.13 by either single or multi-locus models in the soybean (William 82) gene models 1.0 (Grant et al., 2010). The functional annotations of reference genes were obtained from SoyBase and phytozome and candidate genes were selected manually based on their expression pattern in roots.

Conclusions

A total of 23 and 27 QTNs linked with GR, WSL and RL were detected using MLM and mrMLM approaches in all three environments (2018, 2019, and Combined). Among these QTNs, six QTNs, viz., qGR3-1, qGR8-1, qWSL8, qWSL13-1, qRL-8 and qRL15-4 were identified in all three environments via MLM and three of them (qGR8-1, qWSL13-1 and qRL-8) were also detected using mrMLM method. On the base of MLM and mrMLM comparative analysis, the most stable QTN8 was linked with RL and GR and it was confirmed through MLM and mrMLM methods in all environments. Eventually, the 1.0 Mb region around QTN8/ Gm08-12368145 and Gm13-35394159 ± 500 kb was utilized for candidate gene prediction. Based on the expression data, Glyma08g156800, Glyma08g16000, Glyma08g162700 and Glyma13g249600 were regarded as most possible candidate genes which regulate drought tolerance in soybean. Although, functional validation is further needed to define its role in drought tolerance. In short, this study broadens our vision to understand the genetic mechanisms regulating drought tolerance in soybean and it will be helpful in future to breed soybean cultivars with increased drought tolerance via gene cloning and marker assisted selection.