In this study, drought susceptible and tolerant soybean genotypes were obtained from different research institutions (Table 1), were evaluated for their responses to drought stress using physiological, biochemical, and molecular approaches. A two factorial experiment with three replications was conducted in a completely randomized block design (RCBD), with genotype (first factor) and treatment (well-watered and drought stress) as the second factor. The pot trial was set up at the research station of King Abdulaziz University, Jeddah, Saudi Arabia (21°32′36″N, 39°10′22″E; 12 meters above sea level).

Seeds of Glycine max (soybean) cultivars were surface-disinfected using 2% sodium hypochlorite for 3 min, followed by thorough rinsing with double distilled water. Sterilized seeds were sown in 25 cm diameter plastic pots filled with 2:1:1 (v/v) mixture of loam soil, sand, and well-decomposed farmyard manure to facilitate sufficient aeration and drainage. Pots were grown under controlled greenhouse conditions and irrigated regularly to facilitate uniform germination and seedling establishment. Two plants per pot were retained and for each treatment five pots were used. The glasshouse temperature was adjusted at 28 ± 2°C during the day, and 22 ± 2°C at night, with a relative humidity level of 60–70% and a 14 h light/10 h dark photoperiod following Makbul et al. (2011). The experiment was continued for 45 days after sowing, with the drought stress period lasting 15 days, which induced measurable stress responses without causing irreversible damage. Besides, soybean plants were maintained under well-watered conditions until the V₃ vegetative stage. Afterward, the drought stress was imposed by withholding irrigation or maintaining soil moisture at 30–40% field capacity (Board and Kahlon, 2011). Soil moisture at 30–40% field capacity was measured using the gravimeter. Field capacity at 100% was measured after saturating the soil and allowing free drainage and oven-drying a subsample at 105°C for 24 h. The required moisture level was maintained by weighing the pot. Furthermore, soil moisture at 30–40% of field capacity was calculated as 30–40% of the water content at 100% field capacity. It was maintained throughout the experiment and the measurements for FC were taken using the following equation,

For physiological, biochemical and genetic comparisons of genotypes the measurements were taken at once, at the end of the drought stress period, when stress responses were fully expressed before permanent wilting (Makbul et al., 2011). For each treatment the data were taken from three random pots, and averaged before analysis. Besides, the effect of positional error was minimized by changing the positions of pots,

Activities of catalase (CAT), peroxidase (POD), and superoxide dismutase (SOD) were calculated spectrophotometrically following the protocol of Alici and Arabaci (2016). Fresh soybean leaf samples were rinsed, frozen at −20°C, and homogenized in 100 mM sodium phosphate buffer (pH 7.0) having 1 mM ascorbic acid and 0.5% (w/v) PVP. Afterward, the homogenate was incubated at 4°C for 5 min, filtered through muslin cloth, and centrifuged at 5000 × g for 5 min to attain the supernatant as the enzyme extract. The CAT activity from extract was quantified by monitoring the decline in absorbance at 240 nm due to H₂O₂ decomposition (Chance and Maehly, 1955). Besides, POD activity was assessed via the increase in absorbance at 420 nm using 4-methylcatechol as substrate (Kar and Mishra, 1976). Moreover, the SOD activity was measured based on inhibition of NBT photoreduction, with absorbance read at 560 nm (Beauchamp and Fridovich, 1971). The enzyme activities were expressed in standard enzyme units. Proline content was quantified following the acid ninhydrin method described by Bates et al. (1973), in which fresh leaf tissues were homogenized in 3% sulfosalicylic acid, reacted with acid ninhydrin and glacial acetic acid. Subsequently heated at 100°C for 1 h, and the chromophore was extracted with toluene before measuring absorbance at 520 nm. Moreover, Glycine betaine (GB) was measured following the protocol of Grieve and Grattan (1983), where plant material dried in oven was extracted in deionized water. Afterward, treated with iodine–potassium iodide reagent and subjected to incubation at 4°C until precipitate formation. The precipitate was dissolved in 1,2-dichloroethane to record the absorbance at 365 nm.

Photosynthetic rate (Pn) was calculated by a portable Infrared Gas Analyzer (IRGA; ADC Bioscientific, UK). The readings for Pn were taken between 09:00 and 11:00 am to minimize diurnal variation. Fully expanded upper canopy leaves were clamped in the leaf chamber with photosynthetically active radiation (PAR) ranged approximately 800–1200 μmol m⁻² s⁻¹, leaf temperature maintained at 25–30°C, and ambient CO₂ concentration at approximately 380–400 μmol mol⁻¹. Besides, relative humidity followed natural conditions, and data were taken when gas-exchange parameters stabilized to ensure accurate CO₂ assimilation values. Furthermore, total chlorophyll was determined following Wellburn (1994) method. Fresh leaf tissue (0.2 g) was extracted in 10 mL of 80% acetone and centrifuged at 10,000 rpm for 10 mins. Besides, spectrophotometric absorbance was recorded at 645 and 663 nm using 80% acetone as blank. Besides, total chlorophyll (mg L⁻¹) was calculated as (20.2 × A645 + 8.02 × A663), and expressed as g kg⁻¹ of fresh weight using the formula: total chlorophyll (g kg⁻¹ FW) = (mg L⁻¹ × extract volume in liters)/fresh weight in kilograms. Cell membrane stability percentage (CMSP) was calculated using the electrolyte leakage method described by Ibrahim and Quick (2001). Leaf discs were incubated in deionized water and initial electrical conductivity E1 was measured, followed by autoclaving to release total electrolytes to record final electrical conductivity E2. CMSP was calculated based on the ratio of initial to final conductivity values using the formula, [CMSP = 1−(E1/E2) ×100]. Relative water content (RWC) was measured following the protocol of Barrs and Weatherley (1962). Fresh leaf samples were weighed immediately, rehydrated to full turgidity, and then oven-dried to constant weight. RWC was calculated using fresh (FW), turgid (TW), and dry weights (DW) as

Leaf samples for transcript analysis were collected at the onset of drought symptoms, including wilting, leaf rolling, turgidity loss, slight drooping of upper leaves, and a dull green to pale-coloration as compared with well-watered control plants. Besides, the sampling was performed at the initiation of early stress symptoms to capture initial transcriptional responses before the severe stress-induced damages. Collected leaf samples from soybean for transcript analysis were immediately frozen in liquid nitrogen to store at −80°C. Besides, total RNA was extracted using the RNeasy Plant Mini Kit (Qiagen, Germany) and reverse-transcribed from 2 µg RNA using the QuantiTect Reverse Transcription Kit (Qiagen, Germany). Furthermore, Quantitative real-time PCR (qRT-PCR) was conducted using SYBR Green PCR Master Mix (BioFACT, Korea) to quantify the expression of drought-responsive genes, with GmActin1 as the internal reference. All reactions were executed in triplicate (biological and technical), and expression levels were calculated using the 2^−ΔΔCt method (Chen et al., 2013). The list of primers used is indicated in Table 2.

The experiment was performed in RCBD with two factorial arrangement of treatments. Two-way ANOVA was performed using Statistix version 8.1 (McGraw-Hill, 2008) at a 5% probability level to assess the effects of treatments, water regimes (control and drought) and genotypes on the measured traits. Each treatment combination was replicated three times (blocks). When ANOVA indicated significant differences, means were separated using the Least Significant Difference (LSD) test at p≤ 0.05. Besides, prior to analysis, the assumptions of normality and homogeneity of variance were examined.

Moreover, RStudio version 1.1.456 (RStudio Team, 2020) was used to establish correlogram, heatmap clusters and principal component analysis (PCA) biplot. The R packages “factoextra” and “FactoMineR” were used to establish the PCA biplot. Besides, the heatmap cluster dendrogram was constructed using the R packages, “pheatmap” and “complex Heatmap”. The R packages “GGally” and “ggplot2” were used to execute Pearson correlation.

Drought stress significantly (p ≤ 0.05) reduced chlorophyll content, photosynthetic rate (Pn), and cell membrane stability (CMS) across all soybean genotypes, with pronounced genotypic variation as evident in Figure 1. Under drought conditions, chlorophyll content was the lowest in susceptible genotypes ‘Anjasmoro’, ‘Grobogan’, and ‘Dering-1’, remaining below 0.6 g kg⁻¹ (Figure 1). In contrast, tolerant genotypes maintained comparatively higher chlorophyll levels, with ‘Ajmeri’, ‘NARC-21’, ‘DMX4561’, and ‘Rawal’ showing drought values ranging from about 1.3 to 1.6 g kg⁻¹, and ‘Rawal’ recording the highest chlorophyll content under stress (Figure 1).

Photosynthetic rate (Pn) followed a corresponding pattern. Drought stress caused a noteworthy decline in Pn among all genotypes; however, susceptible genotypes such as ‘Anjasmoro’ and ‘Grobogan’ exhibited the lowest rates under drought (approximately 17–18 μmol m⁻² s⁻¹) as shown in Figure 1. Conversely, the tolerant genotypes kept higher Pn activity, with ‘Ajmeri’, ‘NARC-21’, ‘DMX4561’, and ‘Rawal’ maintaining Pn values approximately between 24 and 28 μmol m⁻² s⁻¹, and ‘Rawal’ showing the highest Pn under drought conditions (Figure 1).

Cell membrane stability was also significantly affected by drought stress. CMS values declined across all genotypes, but susceptible genotypes (‘Anjasmoro’, ‘Grobogan’, and ‘Dering-1’) exhibited the lowest CMS under drought, remaining close to 55–60%. In contrast, tolerant genotypes maintained higher CMS values, with ‘Ajmeri’, ‘NARC-21’, ‘DMX4561’, and ‘Rawal’ showing drought CMS values of approximately 68–75%, and ‘NARC-21’ and ‘Rawal’ ranking highest among the genotypes (Figure 1). Overall, Figure 1 demonstrates that genotypes with greater drought tolerance consistently maintained higher chlorophyll content, Pn, and CMS under water-limited conditions.

As shown in Figure 1, drought stress significantly (p ≤ 0.05) induced the accumulation of osmoprotectants, glycine betaine (GB) and proline alongside a variable decline in relative water content (RWC) across all soybean genotypes. GB content increased significantly (p≤ 0.05) under drought stress in all genotypes (Figure 1). Under drought conditions, GB content ranged approximately from 85 to 145 µg g⁻¹ FW. The lowest GB accumulation was observed in ‘Grobogan’ (approximately 80–85 µg g⁻¹ FW) and ‘Dering-1 (approximately 75–80 µg g⁻¹ FW), whereas ‘Anjasmoro’ showed a moderate increase (90 µg g⁻¹ FW) as shown in Figure 1. In contrast, drought-tolerant genotypes accumulated substantially higher GB levels, with ‘Ajmeri’ (125 µg g⁻¹ FW), ‘NARC-21’ (135 µg g⁻¹ FW), DMX-4561 (140 µg g⁻¹ FW), and ‘Rawal’ (145 µg g⁻¹ FW) showing the highest accumulation (Figure 1).

Proline content also increased significantly (p≤ 0.05) under drought stress among all genotypes (Figure 1). Under drought conditions, proline content ranged approximately from 30 to 47 µg g⁻¹ FW (Figure 1). The lowest proline accumulation was recorded in ‘Anjasmoro’ (0 µg g⁻¹ FW) and ‘Grobogan’ (approximately 28–30 µg g⁻¹ FW), followed by ‘Dering-1’ (32 µg g⁻¹ FW). Higher proline accumulation was observed in ‘NARC-21’ (42 µg g⁻¹ FW) and ‘Ajmeri’ (44 µg g⁻¹ FW), while the maximum proline content was recorded in ‘DMX-4561’ (approximately 46–47 µg g⁻¹ FW) and ‘Rawal’ (approximately 43–45 µg g⁻¹ FW) as indicated in Figure 1. The letter groupings indicate significant genotypic differences under drought stress (LSD, p≤ 0.05).

Besides, RWC was substantially maintained in tolerant genotypes, with ‘Rawal’, ‘Ajmeri’, and ‘NARC-21’ exhibiting values >80% under drought. In contrast, ‘Anjasmoro’ and ‘Grobogan’ showed significant water loss (RWC< 65%) (Figure 1). These results demonstrate that the synthesis of compatible solutes and water retention are critical components of drought adaptation in soybean.

The activities of antioxidant enzymes catalase (CAT), superoxide dismutase (SOD), and peroxidase (POD) showed significant (p ≤ 0.05) genotypic variation under drought stress (Figure 1). Drought-tolerant genotypes exhibited a stronger enzymatic response compared to susceptible ones. Catalase activity was markedly higher in ‘NARC-21’ (15.3 U/mg protein) and ‘DMX4561’ (13.1 U/mg protein), followed closely by ‘Rawal’ and ‘Ajmeri’, indicating effective ROS scavenging under stress (Figure 1). In contrast, susceptible genotypes ‘Anjasmoro’, ‘Grobogan’, and ‘Dering-1’ maintained significantly lower CAT activity (<9 U/mg protein) (Figure 1). A similar trend was observed for SOD, where ‘Rawal’ and ‘Ajmeri’ recorded values above 50 U/mg protein, whereas ‘Grobogan’ and ‘Dering-1’ showed values below 40 U/mg (Figure 1). The POD activity increased significantly (p≤ 0.05) in all genotypes under drought, but the most pronounced elevation was seen in ‘NARC-21’, ‘Rawal’, ‘Ajmeri’, and ‘DMX4561’ (0.70-0.80 U/mg protein), whereas ‘Anjasmoro’, ‘Dering-1’ and ‘Grobogan’ exhibited limited increases (0.40-0.45 U/mg protein) (Figure 1). These results demonstrate that enhanced antioxidant activity in tolerant genotypes contributes to their superior drought resilience by mitigating oxidative stress.

The expression profiles of drought-responsive genes revealed a consistent difference between tolerant (‘Ajmeri’, ‘NARC-21’, ‘DMX4561’, ‘Rawal’) and susceptible (‘Anjasmoro’, ‘Grobogan’, ‘Dering-1’) soybean genotypes as shown in Figure 2. The transcription factor gene, GmDREB2 showed higher expression in tolerant genotypes, with fold increases ranging from ~2.0 to 3.3-fold relative to control, compared to ~1.8 to 2.0-fold in susceptible genotypes, indicating stronger activation of downstream stress defense pathways (Figure 2). Moreover, GmLEA-D11, associated with dehydration tolerance, also displayed greater upregulation in tolerant genotypes (~2.8 to 3.6-fold) than in susceptible ones (~1.7 to 2.1-fold) as indicated in Figure 2. The proline biosynthesis gene GmP5CS exhibited markedly higher expression (Figure 2) in tolerant genotypes (~2.9 to 3.8-fold) than in susceptible genotypes (~1.8 to 2.2-fold), closely reflecting the higher proline accumulation observed in biochemical analyses (Figure 1). Genes encoding antioxidant enzymes, GmSOD1, GmPOD1, and GmCAT1, also followed corresponding trend, with tolerant genotypes showing fold increases of ~2.7 to 3.5-fold compared to ~1.6 to 2.0-fold in susceptible genotypes (Figure 2), consistent with the elevated enzymatic activities recorded in Figure 1. The osmolyte biosynthesis gene GmBADH2, responsible for GB production, was more strongly expressed in tolerant genotypes (~2.8 to 3.4-fold) than in susceptible genotypes (~1.7 to 2.0-fold) (Figure 2), aligning with their higher GB content observed in Figure 1. For water transport, GmPIP2;9 expression was maintained or slightly increased in tolerant genotypes (~1.3 to 1.6-fold) but declined in susceptible ones (~0.8 to 1.0-fold) as indicated in Figure 2, indicating reduced aquaporin-mediated water movement in sensitive backgrounds. The chlorophyll biosynthesis gene GmCHLG was downregulated in all genotypes under drought but remained higher in tolerant lines (~0.8 to 1.0-fold) compared to susceptible lines (~0.5 to 0.7-fold) (Figure 2), corresponding to their better chlorophyll retention (Figure 1). These consistent fold-change patterns across genes reinforce the earlier physiological, biochemical, and statistical findings (Figure 1), demonstrating that drought tolerance in soybean is regulated by coordinated transcriptional activation of genes involved in stress signaling, osmotic adjustment, antioxidant defense, water transport, and chlorophyll biosynthesis.

The multivariate statistical analysis provided further validation and deeper insights into the physiological and biochemical responses of the soybean genotypes under drought stress. The correlation matrix and pairwise scatter plots in Figure 3 revealed highly significant positive associations among drought-responsive traits. Under drought conditions, strong correlations were observed between proline and GB (r = 0.991***), as well as between proline and antioxidant enzymes such as SOD (r = 0.952***), CAT (r = 0.918***), and POD (r = 0.961***) (Figure 3). Similarly, chlorophyll content, Pn, and RWC also exhibited strong intercorrelations (r > 0.95), indicating a coordinated physiological and biochemical adjustment under water deficit condition (Figure 3). These interdependent relationships suggested that the accumulation of osmoprotectants and enhanced antioxidant activities directly contributed to CMS and photosynthetic performance during drought stress. Moreover, principal component analysis (PCA) further delineated the underlying variation among genotypes and traits (Figure 4). Under drought, tolerant genotypes such as ‘Ajmeri’, ‘DMX4561’, ‘Rawal’, and ‘NARC-21’ aligned positively with vectors corresponding to proline, GB, SOD, CAT, POD, and CMS, indicating strong trait associations and superior adaptive capacity (Figure 4). In contrast, susceptible genotypes like ‘Anjasmoro’, ‘Grobogan’, and ‘Dering-1’ grouped on the opposite side of the PCA space, reflecting poor association with these drought-responsive traits and reduced resilience under stress (Figure 4). The positioning of tolerant genotypes near the osmolyte and enzyme vectors implies their higher contribution to the total variation explained by principal components, especially under drought.

Furthermore, the treatment-based PCA biplot illustrated the clear divergence between control and drought conditions (Figure 5). Control plants clustered around traits such as chlorophyll, RWC, and Pn, indicating high baseline photosynthetic efficiency and hydration (Figure 5). In contrast, drought-stressed genotypes aligned closely with antioxidant enzymes and osmolytes, suggesting that these traits become the major determinants of stress tolerance under water-deficit conditions (Figure 5). Notably, the drought cluster was dominated by tolerant genotypes, confirming their shift towards biochemical defense mechanisms when subjected to stress. The shape and spread of ellipses in the biplot further indicate that drought imposes a distinct directional shift in the physiological profile of tolerant genotypes compared to susceptible ones.

Genotype-based PCA showed that tolerant genotypes were grouped in close proximity to key drought-responsive traits, including proline, GB, CAT, and SOD, demonstrating their consistent upregulation (Figure 6). On the other hand, susceptible genotypes such as ‘Anjasmoro’ and ‘Grobogan’ occupied peripheral positions in the PCA space and remained distant from stress-related vectors, highlighting their limited adaptability (Figure 6). ‘Dering-1’, although slightly closer to the center, still displayed weak associations with protective traits confirming this genotype as drought-sensitive. This genotypic separation provides strong multivariate evidence that drought tolerance in soybean is a function of both osmotic regulation and oxidative defense, with tolerant genotypes excelling in both domains.

The two-way hierarchical clustering presented in Figure 7 rectified the PCA findings through heatmap visualization of trait intensities across genotypes. Under drought stress, genotypes displayed clear differential responses. ‘NARC-21’, ‘Rawal’, ‘DMX4561’, and ‘Ajmeri’ generally exhibited higher relative values (orange to red gradients) for several drought-responsive traits, particularly proline, GB, antioxidant enzymes (CAT, POD, SOD), and CMS, although the magnitude of response varied among traits and genotypes (Figure 7). In contrast, ‘Anjasmoro’, ‘Dering-1’, and ‘Grobogan’ were characterized by comparatively lower values (green to yellow gradients) for many traits, most notably CMS, RWC, and antioxidant-related parameters (Figure 7). Hierarchical clustering under drought stress resulted in a clearer separation of genotypes, with tolerant genotypes grouping more closely based on their overall physiological and biochemical traits, while susceptible genotypes formed a separate cluster, reflecting contrasting drought response patterns (Figure 7).

Together, these statistical analyses provide strong, convergent evidence that drought tolerance in soybean is regulated by a combined network of antioxidant activity, osmoprotectant accumulation, and physiological stability. The close clustering of tolerant genotypes with these stress-mitigating traits across correlation, PCA, and heatmap analyses illustrates their superior resilience and highlights their suitability for drought-prone environments.