A natural population of 207 soybean accessions collected from 11 provinces was used to conduct association mapping analysis. All of the soybeans were obtained from the Soybean Research Institute, College of Pharmacy and Life Science, Jiujiang University (Jiujiang, China).

Twenty healthy and good-quality seeds were selected and sterilized with 1% NaClO₂ for 30 s and then rinsed with distilled water (ddH₂O) three times. Then the seeds were placed on two sheets of filter paper in sterilized 9 cm petri dishes. Drought stress was created by the addition of 15 ml polyethylene glycol 6000 (PEG-6000, 20% w/v) (Zhao et al., 2022b), whereas 15 ml ddH₂O was employed in control conditions. The plates were incubated in a germination chamber at 25 ± 1°C, and the experiment was conducted in the dark for 6 days. Each processing was repeated three times. Subsequently, the treated seeds were rinsed by the corresponding solution (ddH₂O or 20% (w/v) of PEG-6000), and 2 ml of the corresponding solution was added every day. The soybean seeds were considered to be germinated when the radicle length was beyond half of the longitudinal length of the seed. The count of germinated seeds was recorded every 24h for 6 days.

The germination traits were evaluated by the following formulas: Germination index (GI) = Σ (Gt/Dt), where Gt is the accumulated number of germinated seeds on day t, and Dt indicates the germination days corresponding to Gt in days; germination potential (GP) = (number of germinated seeds on the third day/total number of seeds) × 100%; germination rate (GR) = (number of germinated seeds on the sixth day/total number of seeds) × 100%. On the sixth day, for each cultivar, three normal and uniform growth seedlings were selected and used for phenotyping. The lengths of these radicles were measured by ruler as their seedling length (SL), and their fresh weights were measured by electronic balance and recorded as seedling fresh weight (SFW). The relative value of these traits was calculated by the formula: trait value under treatment/trait value of control × 100 (%), and then the relative germination index (RGI), relative germination potential (RGP), relative germination rate (RGR), relative seedling length (RSL), and relative seedling fresh weight (RSFW) were obtained.

The descriptive statistics, analysis of variance (ANOVA), and correlation analysis of the five germination-related traits were conducted using the programs MEANS, CORR, and PROC GLM by SAS 9.4 (SAS Institute, Cary, NC). The broad-sense heritability was estimated by the equation: h 2 = σ g 2 / ( σ g 2 + σ e 2 / r ) , where σ g 2 (g=1, 2, 3,…207) is the genotypic variance, σ e 2 represents random error variance, and r is the number of replicates (Nyquist and Baker, 1991).

The 207 soybean accessions were genotyped by ZDX1 (Sun et al., 2022), a high-throughput functional array, and a total of 123,540 SNPs were obtained (unpublished data). After filtering out SNPs with minor allele frequencies (MAF) ≤ 5%, a subset of 95,043 high-quality SNPs was used for further analyses.

The principal components analysis (PCA) and neighbor-joining tree were constructed using TASSEL 5.0 software (Bradbury et al., 2007). Linkage disequilibrium (LD) analysis was estimated using HAPLOVIEW 4.2 software. The 500-kb sliding window along each chromosome was used to calculate LD between all pairs of SNPs. The LD was estimated with the squared correlation coefficient r ². The LD decay distance was defined as the chromosomal distance where the r ² dropped to half of its maximum value (Barrett et al., 2005).

The RTM-GWAS procedure (https://github.com/njau-sri/rtm-gwas) was performed for QTL mapping (He et al., 2017). This association mapping was conducted in two stages. In the first stage, a single-locus model association analysis based on a simple linear model with a threshold of P = 0.05 was performed to preliminarily screen markers. In the second stage, the significance level of P = 0.05 under the multi-locus model was used for the preselected markers to identify genome-wide QTLs through stepwise regression. Manhattan and quantile–quantile plots were generated using the R software. The nomenclature (McCouch et al., 1997) with modifications was used to name the QTL in this study. Such as qRGI-3-1, q represents the QTL; RGI represents the relative germination index, -3 represents chromosome 3; -1 represents the first QTL on that chromosome. Genes located within ± 200 kb (the average LD decay distance across all 20 chromosomes) of the QTLs were identified as candidate genes. The Phytozome database (http://www.phytozome.org/) and the SoyBase database (http://www.soybase.org/) were used to annotate the candidate genes.

Two accessions, JJS98 (high drought-tolerant, DT) and JJS294 (high drought-sensitive, DS), were selected from the association panel for transcriptomic analysis. Briefly, the seeds of JJS98 and JJS294 were germinated under normal conditions, and after five days the germinating seeds were transferred into the ddH₂O or 20% (w/v) of PEG-6000. For each accession, the root tips were respectively collected at 0h, 12h, and 24h after drought treatment, with three biological repetitions. The total RNA was isolated using a Plant RNA Extract Kit (TianGen, Beijing, China) according to the manufacturer’s instructions. RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States) ( Supplementary Table S1 ). In total, 18 RNA samples were submitted to Biomarker Technologies Co., Ltd. (Beijing, China) for sequencing using the Illumina HiSeq 4000 sequencing platform.

The raw paired-end reads were first filtered with Fastp software to obtain high-quality clean data and then aligned to the soybean reference genome Williams 82 (Glycine max v2.1 genome) using Hisat2 v2.0.5 software. The FPKM (fragments per kilobase of transcript per million fragments mapped) method was utilized to normalize and estimate gene expression values. DEGs were identified using the false discovery rate (FDR) < 0.01 and |log₂ ^FoldChang| ≥2. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using BMKCloud (www.biocloud.net).

The co-expression network was constructed in R using the WGCNA (v1.47) package. For high reliability of the results, the FPKM greater than 1 of all DEGs was selected for WGCNA analysis. According to the correlation of gene expression, the gene clustering tree was constructed, and then the genes with similar expression patterns were classified into the same module. The modules were merged with the threshold value of module eigenvalue similarity > 0.25, and the minimum gene number of modules was set as 30 (Wang et al., 2024). The value of the soft threshold (β) was selected as 13, the hierarchical clustering tree was a dynamic hybrid tree cut algorithm, and the weighted network was unsigned. Modules that are significantly correlated with specific samples are identified according to the heatmap of the sample expression patterns, and the DEGs in these modules were used for KEGG enrichment analysis on the BMKCloud platform. In each module, the genes with eigengene-based connectivity value (|KME|) > 0.98 were regarded as hub genes.

Five hub genes were selected for qRT-PCR validation. Total RNA was extracted from the same samples that were used for sequencing. First-strand cDNA was synthesized using a HiScript IV RT SuperMix for qPCR (+gDNA wiper) (Vazyme Biotech Co., Ltd, Nanjing, China). The qRT-PCR experiment was performed on a CFX Connect (TM) Real-Time PCR Detection System (BIO-RAD, Singapore) using SupRealQ Purple Universal SYBR qPCR Master Mix (U+) (Vazyme Biotech Co., Ltd, Nanjing, China) according to the instructions. Soybean gene GmActin11 (Glyma.18G290800) was the internal control gene, and fold change was calculated using the 2^−ΔΔCT method (Livak and Schmittgen, 2001). The detail of primers was shown in Supplementary Table S2 .

Five germination-related traits, namely GI, GP, GR, SL, and SFW, were measured for germinating seeds of the 207 soybean accessions under ddH₂O (CK) or 20% (w/v) of PEG-6000 (PEG) ( Table 1 ). The mean values of GI, GP, GR, SL, and SFW for the seeds treated with CK were 17.72, 0.59, 0.67, 7.59, and 2.02, while the means for the PEG were 6.21, 0.12, 0.56, 2.40, and 1.41, respectively. Obviously, all of these traits were substantially higher under the control group than under drought-stressed conditions, indicating that after 20% PEG-6000 drought treatment, the germination of soybean germplasms was greatly inhibited. Moreover, a descriptive statistical analysis was also conducted on RGI, RGP, RGR, RSL, and RSFW of 207 soybean accessions. The average values of RGI, RGP, RGR, RSL, and RSFW were 0.36, 0.22, 0.83, 0.37, and 0.71, respectively. The coefficients of variation (CVs) of the RGI, RSL, RGR, and RSFW were 89.69%, 63.53%, 59.60%, and 26.04%, respectively. The CVs of RGP were 168.91%, surpassing the 100% threshold. These indicated that the population used in this study exhibits rich variation under drought conditions. Besides, the heritability of all five germination-related traits was high and ranged from 65.25% to 97.38% ( Table 1 ), indicating that most of the phenotypic variance of seed-droughting tolerance was mainly controlled by genetic variation.

An ANOVA indicated that there were significant (P < 0.001) differences in the five germination-related traits among genotypes, treatments, and genotype-by-treatment interactions. However, there was no significant variation in the five germination-related traits among treatments (replications) ( Table 2 ), indicating experimental data is relatively reliable.

To explore the correlation among RGI, RGP, RGR, RSL, and RSFW, Pearson’s correlation analysis was conducted based on the means for the 207 accessions ( Table 3 ). The result showed that there were significant positive correlations for all the five germination-related traits. In particular, RGI exhibited a significant positive correlation with RGP, RGR, and RSL (P < 0.001), with large correlation coefficients of 0.86, 0.82, and 0.62, respectively ( Table 3 ). These results demonstrated that there is a certain correlation between different drought-tolerant traits.

The neighbor-joining tree and PCA analyses were used to determine the population structure of 207 soybean accessions. Based on the neighbor-joining tree analysis, these 207 accessions were clearly grouped into three groups ( Supplementary Figure S1A ). For the three groups, the group 3 accessions were primarily from Southern China, the Huanghuai region, and local collections, while the Groups 2 and 1 accessions mainly originated from Hebei, Fujian, and Guangxi provinces. Furthermore, the result of the PCA analysis was consistent with the neighbor-joining tree analysis ( Supplementary Figure S1B ).

In the present study, the whole-genome LD decay rate of the population was estimated. Calculation results showed that the maximum r ² value was 0.71. Then, when the r ² value dropped to half of its maximum value, the LD decay distance reached approximately 200 kb ( Supplementary Figure S2 ).

The RTM-GWAS method was used to perform association analysis on five germination-related traits. A total of 58 QTL were detected with a significance level of P < 0.05 ( Figure 1 and Supplementary Table S3 ). The -log P values of these QTLs ranged from 2.03 to 23.22 and explained 1.11% to 40.63% of the genetic contribution (R ²).

Among these QTLs, there were 12 QTLs for RGI, and the phenotypic variation of these loci was 1.47% to 24.12%, which together explained a total of 84.36% of genetic contribution. Eight QTLs were identified for RGP, and the phenotypic variance varied from 2.20% to 44.36%, collectively accounting for 71.38% of the genetic contribution. There were 17 QTLs for RGR, and the phenotypic variation was 1.11%–19.28%, explaining a total of 88.35% of genetic contribution. Thirteen QTLs were identified for RSL, and the phenotypic variance varied from 2.31% to 12.58%, collectively accounting for 84.39% of the genetic contribution. For 8 RSFW QTLs, the range of phenotypic variation was 2.99% to 40.63% and explained a total of 78.22% of genetic contribution ( Supplementary Table S3 ). More importantly, the two pairs of QTLs, qRGI-17–1 and qRGP-17-1, and qRGI-19–1 and qRGP-19-1, were located on the same genomic region. qRGI-17–1 and qRGP-17–1 were mapped to the same SNP of Gm17_BLOCK_37041250_37044388, accounting for 24.12% and 44.36% of the phenotypic variation, respectively. qRGI-19–1 and qRGP-19–1 were co-associated with the SNP marker Gm19_BLOCK_30608933_30612180, explaining 15.95% and 5.77% phenotypic variation, respectively. Out of these QTLs, 16 QTLs were co-localized with the previously reported soybean drought tolerance QTL, and the others were identified in this study ( Supplementary Table S3 ).

In addition, loci with a higher phenotype variation (R ²>10%) were considered as large-contribution QTLs in this study, including qRGI-16-1, qRGI-17-1, qRGI-19-1, qRGP-17-1, qRGR-7-1, qRGR-12-1, qRGR-15-1, qRSL-3-1, qRSFW-5-1, and qRSFW-14-1. These large-contribution QTLs were used for the following candidate gene mining.

Via transcriptome sequencing for the two accessions with contrasting drought tolerance, a total of 137.74 Gb of clean data were obtained from 18 libraries. In different samples, the Q30 percentage was over 94.42%, and more than 90.56% of clean reads were mapped to the soybean genome, while the percentage of unique mapped reads ranged from 88.61 to 93.63% ( Supplementary Table S4 ), indicating that the RNA-seq libraries are of high quality. Based on three biological replicates at different time points after drought treatment, the PCA was conducted on all expressed genes in 18 samples. The first two principal components, PC1 and PC2, accounted for 24.35% and 15.72%, respectively ( Supplementary Figure S3A ). These results indicated that there were different gene expression patterns among the soybean roots under drought stress at different time points ( Supplementary Figure S3B ).

Subsequently, the analysis of differentially expressed genes (DEGs) was performed in two aspects. On the one hand, drought-responsive genes were determined by comparing drought-treated samples with controls for two accessions at two time points. Relative to the gene expression levels at 0h, 7,055 (3,196 upregulated and 3,859 downregulated) and 5,198 (2,509 upregulated and 2,689 downregulated) DEGs were detected in the DS accession at 12h and 24h after drought treatment, respectively. While 5,307 (2,414 upregulated and 2,893 downregulated) and 6,560 (2,792 upregulated and 3,768 downregulated) DEGs were detected in the DT accession at 12h and 24h after drought treatment, respectively ( Figure 2A ). There were 2,041 DEGs continuously responding to drought stress both in DS and DT accessions ( Supplementary Figure S4 ). On the other hand, DEGs between the two accessions at two time points were identified. There were 608 DEGs after 12h of drought treatments and 740 DEGs after 24h between the two accessions, which did not include DEGs that are not drought-responsive ( Figures 2B, C ). The number of DEGs identified after 24h of drought treatments more than 12h suggested that more gene networks are triggered by drought stress in the two accessions with time increasing. Of the 608 DEGs after 12h, 213 (35.0%) were shared by DS and DT, while the remaining 284 (46.7%) and 111 (18.3%) were DS- and DT-specific DEGs, respectively ( Figure 2B ). By contrast, the relative proportion of DEGs identified after 24h of drought treatment was equivalent to after 12h, as 226 of 740 DEGs (30.5%) were shared by DS and DT, and 225 (30.4%) and 289 (39.1%) uniquely belonged to DS and DT accessions, respectively. These results indicated that accession-specific responses to drought stress occurred ( Figure 2C ). In addition, 1,183, 1,354, and 1,581 DEGs between DS and DT after 0h, 12h, and 24h of drought treatment were separately obtained, and 478 genes showed significant expression differences between DS and DT at all time points ( Figure 2D ). Two DEGs were commonly detected in all pairwise comparison groups ( Supplementary Figure S5 ).

To know the functional information of the DEGs between DS and DT at two time points, the GO enrichment and KEGG analyses were performed. GO enrichment analysis indicated that 1,354 DEGs identified after 12h of drought treatment were mainly concentrated in hydrogen peroxide catabolic process (GO:0042744), response to oxidative stress (GO:0006979), root morphogenesis (GO:0010015), and signal transduction (GO:0007165) ( Supplementary Figure S6A ). Interestingly, the GO terms, such as hydrogen peroxide catabolic process (GO:0042744) and response to oxidative stress (GO:0006979), were also significantly enriched after 24h of drought treatment ( Supplementary Figure S6B ). While the GO term of glutathione metabolic process (GO:0006749) was only significantly enriched after 24h ( Supplementary Figure S6B ). KEGG analysis showed that three identical pathways were significantly enriched after 12h and 24h of drought treatment, including phenylpropanoid biosynthesis, starch and sucrose metabolism, and fatty acid degradation ( Figures 3A, B ). In the phenylpropanoid biosynthesis pathway, 14 of 32 DEGs after 12h of drought treatments and 19 of 42 DEGs after 24h were annotated as peroxidase activity, which corresponded to the results of GO enrichment analysis ( Supplementary Table S5 ). These observations indicated that the enrichment of antioxidant and secondary metabolite pathways was related to the ability of soybean drought tolerance at the germination stage.

Besides, 2,041 drought-responsive DEGs shared by DS and DT at two time points were also analyzed. GO analysis revealed a stark enrichment of drought-responsive DEGs in several functional categories, namely hydrogen peroxide catabolic process (GO:0042744), response to oxidative stress (GO:0006979), and UDP-glycosyltransferase activity (GO:0008194) ( Supplementary Figure S7 ). The KEGG annotations indicated that the pathways enriched with these drought-responsive DEGs were closely related to phenylpropanoid biosynthesis, glutathione metabolism, isoflavonoid biosynthesis, and the MAPK signaling pathway ( Supplementary Figure S8 ).

To further identify the hub genes or modules that are associated with the drought stress response in soybean roots, all DEGs between two accessions and between treatments at different times were applied for WGCNA analysis. Based on the co-expression patterns of these DEGs, the gene cluster dendrogram was constructed, and 13 distinct modules (labeled with different colors) were identified ( Figure 4A ). Next, the expression patterns of each module in different samples were analyzed based on the module eigenvalues ( Figure 4B ). Compared to other samples, the DT-12 and DT-24 samples were most significantly correlated with the pale turquoise and dark olive green modules, respectively. While both of these two samples were significantly correlated with the skyblue2 module. In contrast, the medium orchid module showed positive correlation with both DS-12 and DS-24 samples ( Figure 4B ). Therefore, the above-mentioned four co-expression modules were considered as key modules for further analysis. Analysis of the expression mode of DEGs showed that there was significant specificity in the expression of DEGs in each key module. For example, the pale turquoise and dark olivegreen modules respectively represent DEGs that showed high expression specifically at 12h and 24h after drought treatment in DT accession ( Figures 5A, B ). Additionally, DEGs in the skyblue2 and medium orchid modules were continuously highly expressed after drought treatment in DT and DS accessions, respectively ( Figures 5C, D ).

Furthermore, KEGG pathway enrichment analysis showed that the highly enriched pathways in the paleturquoise module included phenylpropanoid biosynthesis, starch and sucrose metabolism, and peroxisome ( Supplementary Figure S10A ). The highly enriched pathways in the dark olive green module included phenylpropanoid biosynthesis, fatty acid degradation, and the MAPK signaling pathway ( Supplementary Figure S10B ). Similarly, the skyblue2 module showed enrichment of KEGG pathways related to glutathione metabolism, phenylpropanoid biosynthesis, and plant hormone signal transduction ( Supplementary Figure S10C ). The medium orchid module is involved in glutathione metabolism, isoflavonoid biosynthesis, and zeatin biosynthesis pathways ( Supplementary Figure S10D ).

Additionally, a total of 40 hub genes were discovered in the four modules according to the kME value (|KME| > 0.98) and gene annotation, of which 28, 3, 4, and 5 existed in the dark olive green, medium orchid, pale turquoise, and skyblue2 modules, respectively ( Supplementary Table S6 ). Most strikingly, in the dark olive green module, Glyma.13G229300 encodes abscisic acid receptor PYL11-related protein, and Glyma.08G188300 encodes SNF1-related protein kinases 2 (SnRK2s), which are two key elements of ABA signaling that play an essential role in the abiotic stress response. Four genes (Glyma.05G226700, Glyma.16G181200, Glyma.03G066700, and Glyma.03G066600) were related to the oxidation-reduction process that is essential for basic life functions, including photosynthesis and respiration. Two transcription factors, BHLH57 and NF-YA8, respectively encoded by Glyma.12G178500 and Glyma.17G051400, were identified as hub genes. Furthermore, Glyma.04G231400, encoding an HD-ZIP family transcription factor; Glyma.09G052800, encoding an ethylene-responsive transcription factor; Glyma.12G091200 encoding an NAC transcription factor; and Glyma.13G124900, encoding a GRAS domain transcription factor, were identified as hub genes in the medium orchid, pale turquoise, and skyblue2 modules, respectively.

Among the 10 large-contribution QTLs, a total of 294 genes were retrieved within ± 200 kb (average LD decay distance across all 20 chromosomes) of these QTLs, according to the soybean genome Williams 82 (Glycine max v2.1 genome) ( Supplementary Table S7 ). Among the 294 genes detected by RTM-GWAS, three genes are differentially expressed between DS and DT accessions, and 19 genes belong to drought-responsive DEGs in both accessions, according to the transcriptome data ( Supplementary Table S8 ). Herein, these 22 genes were identified as putative candidate genes.

According to the functional annotations, the five candidate genes, Glyma.12G141700, Glyma.15G040000, Glyma.05G049300, Glyma.14G105900, and Glyma.15G041100, were further considered as priority candidate genes. Among them, three transcription factors, Glyma.12G141700, Glyma.05G049300, and Glyma.15G041100, belong to the bHLH gene family, GRAS gene family, and MYB gene family, respectively; Glyma.15G040000 and Glyma.14G105900 encode a cytosolic isoform of UDP-glucuronic acid decarboxylase and an E3 ligase, respectively. Among these genes, the transcriptome data uncovered that the Glyma.15G040000 and Glyma.05G049300 genes were downregulated both in DT and DS after 12 and 24h of drought treatment, whereas Glyma.12G141700 and Glyma.14G105900 were upregulated. Besides, Glyma.15G041100 was downregulated at 12h and upregulated at 24h after drought treatment in DT, while it was upregulated both at 12 and 24h after drought treatment in DS ( Figure 6 ).

To confirm the accuracy of the RNA sequencing results, five priority candidate genes were selected for qRT-PCR analysis. The qRT-PCR results of five genes were consistent with the expression pattern of RNA-seq data ( Figure 6 ). Genes significantly up-regulated in RNA-Seq data also exhibited an upregulation in qRT-PCR, and vice versa, suggesting the reliability and accuracy of our RNA-Seq data.