The experiment was carried out in 2022–2023 at the Sanping Experimental Farm (116°41′E, 39°91′N; Urumqi, Xinjiang, China). Meteorological data (temperature and precipitation) for the soybean growing seasons are shown in Figure 1. The soil is classified as sandy loam, with initial properties as follows: organic matter content 13.8 g kg^-1, total N 0.82 g kg^-1, mineral N 56.65 mg kg^-1, available P 14.1 mg kg^-1, and available K 200.6 mg kg^-1.

The field experiment followed a split-plot arrangement over two consecutive growing seasons. The main plots consisted of two soybean cultivars: Xindadou 27 (a low-yielding cultivar) and Xinnongdou 2 (a high-yielding cultivar). The subplots comprised four nitrogen application levels: 0 (N₀), 120 (N₁₂₀), 180 (N₁₈₀), and 240 kg ha^-1 (N₂₄₀). Among these, 0 kg N ha^-1 served as the nitrogen-free control; 240 kg N ha^-1 represented the conventional high-yield application rate used by local farmers under plastic mulching with drip irrigation; 180 kg N ha^-1 was identified as the optimized threshold based on preliminary experiments, which maintained yield potential while significantly improving nitrogen use efficiency; and 120 kg N ha^-1 represented a moderate nitrogen reduction treatment. All nitrogen fertilizer was applied in a single dose at the beginning pod stage via drip irrigation. Each plot measured 48 m² (4.8 m × 10 m). Each treatment was replicated three times, with a 2 m buffer zone between adjacent plots. Before tillage, superphosphate (containing 19% P₂O₅) was applied at 300 kg ha^-¹. Soybeans were sown on April 27 and harvested from September 6 to 26. The experiment was conducted over two consecutive years. A uniform row spacing of 40 cm was used with a planting density of 33×10⁴ plants ha^-¹. The field was covered with 140 cm wide black plastic mulch, and drip irrigation tapes were laid beneath the mulch at 40 cm intervals (Φ16 mm with a discharge rate of 2.5–3.5 L h^-1). Irrigation was conducted annually from June 17–27 to August 5–15, once every 10–13 days, totaling six events with each irrigation applying 600–700 m³ ha^-1. The total irrigation volume was 3900 m³ ha^-1. Potassium dihydrogen phosphate (containing 34% K₂O and 52% P₂O₅) was applied through fertigation. The total nutrient application rates were 51 kg K₂O ha^-1 and 78 kg P₂O₅ ha^-1, applied in a 1:2 ratio at the initial flowering stage and the beginning of the pod stage, respectively. Two manual weeding were performed during the soybean growing season. All other management practices were consistent with conventional field production.

Root sampling was carried out at five key growth stages: full flowering stage⁠ (R₂), full pod stage (R₄), beginning pod stage (R₅), full seed stage (R₆), and full maturity stage (R₈). For each treatment, the profile trench method was used for sampling. The experiment employed a planting pattern with four rows under a single plastic film. To systematically control the influence of root overlap between adjacent plants, all root samples were collected from the two middle rows. The specific sampling location was centered on a single plant selected from the middle row of plants. From this center point, sampling extended 10 cm along the row direction in both directions (total length 20 cm) and 20 cm perpendicular to the row direction in both directions (total width 40 cm). The excavation depth was 20 cm, resulting in a single-layer soil sampling volume of 16,000 cm³. This volume was determined based on the local sandy loam soil texture and the characteristic that soybean roots are primarily distributed in the 0–60 cm soil layer, ensuring both representativeness of the root samples and practicality for field operation and subsequent root separation. Root samples were collected from the 0–20 cm, 20–40 cm, and 40–60 cm soil layers separately. They were placed on a 2 mm mesh sieve and rinsed with gently flowing water. All root residues retained on the sieve were collected (pre-experiment verification indicated root loss during this operation was minimal, not affecting data accuracy). The roots were then sorted, re-screened, and washed, and the taproots and lateral roots were separated. The root samples were placed in a transparent glass box filled with water. Each root was carefully extended and flattened with forceps to avoid overlap and crossing. Before scanning, the WinRHIZO-2004a system was calibrated using its standard root images to ensure accurate measurement of root length and surface area. The root samples were then scanned using the root scanning analyzer to obtain the relevant parameters. After scanning, the root samples were placed in a 105°C oven for 30 minutes to deactivate enzyme activity, then dried at 80°C to a constant weight, and the root dry weight was recorded. Based on the sampling volume of each soil layer, the root dry weight density (RDD, Equation 1), root length density (RLD, Equation 2), and root surface area density (RSD, Equation 3) were calculated.

During the full flowering stage (R₂), full pod stage (R₄), beginning pod stage (R₅), full seed stage (R₆), and full maturity stage (R₈), ten representative plants were collected from each treatment. Root systems were carefully excavated and rinsed with tap water. Five plants were then separated into constituent organs and dried in an oven at 60°C until constant weight was achieved for biomass assessment. The dried plant material was ground to pass through a 1 mm sieve and stored under dry conditions. Roots from the remaining five plants were preserved at –80°C for later biochemical analysis.

Fresh roots were extracted with an appropriate buffer, and the supernatant was collected for the determination of the following enzyme activities: nitrate reductase (NR), glutamine synthetase (GS), and glutamate synthase (GOGAT). NR activity was measured using the method of Ogawa et al (Ogawa et al. (1999)). with slight modifications: the reaction mixture contained 100 mM KNO₃, 5 mM NaHCO₃, 25 mM PBS, and 0.2 mM nicotinamide adenine dinucleotide (NADH). After the mixture had reacted thoroughly for 30 minutes, 250 μL of 1% sulfanilamide reagent and 250 μL of 1% α-naphthylamine reagent were sequentially added. A standard curve was prepared using sodium nitrite. GS activity was determined according to the method of O’Neal and Joy (O’Neal and Joy (1974)): the reaction system contained 50 mM Tris-HCl (pH 7.5), 4 mM ATP, 80 mM sodium glutamate, 30 mM MgSO₄, 10 mM NH₂OH, and 30 mM cysteine, with γ-glutamyl hydroxamate used to prepare the standard curve. GOGAT activity was assayed by the method of Singh and Srivastava (Singh and Srivastava (1986)): the reaction medium consisted of 100 mM potassium phosphate buffer (pH 7.6), 0.1% (v/v) 2-mercaptoethanol, 100 μM NADH, 2.5 mM 2-ketoglutarate, and 100 mM glutamine.

Plant nitrogen content was determined using the Kjeldahl method: plant organs were digested with concentrated H₂SO₄ and H₂O₂, and nitrogen content was measured using a K9840 automatic Kjeldahl nitrogen analyzer. The total N accumulation in the samples was calculated using the following formula (Equation 4):

Where DMA represents dry matter accumulation, and Nc represents the total nitrogen concentration.

The dynamics of total nitrogen accumulation in soybean were fitted using the Logistic equation as follows (Gao et al., 2021):

Where t is the days after emergence (DAE, d), Y is the plant biomass or total nitrogen amount (kg) at time t, K is the maximum biomass or total nitrogen accumulation (kg), and a and b are constants.

Using differential calculus on Equation 5 yields the following:

Where Δt is the fast accumulation period (FAP) of total nitrogen, t₁ is the start time, t₂ is the end time, and Vt is the mean accumulation rate during this FAP (Equation 6).

With reference to the methodology described by Mao et al. (2018) for quantifying biomass accumulation rates, the beta growth function (Equation 7) was employed to simulate the unimodal curve dynamics of total nitrogen accumulation in both roots and shoots.

Where V (kg ha^-1 d^-1) represents the total nitrogen accumulation rate at days after emergence (DAE) t (d); Te (d) denotes the termination time of accumulation, i.e., the accumulation duration; Tm (d) represents the occurrence time of the maximum accumulation rate Vm (kg ha^-1 d^-1). Tm and Te can be obtained by fitting the total nitrogen accumulation data using Equation 8, while Vm can be calculated using Equation 9.

Where W (kg ha^-1) and Wmax (kg ha^-1) represent the total nitrogen accumulation at days after emergence t and the maximum total nitrogen accumulation, respectively.

At the full maturity stage, ten uniform, disease-free plants per treatment were sampled in three replicates to determine yield components, including pods per plant, grains per plant, and 100-grain weight. Grain yield was measured by harvesting a central area of 6.4 m² (1.6 m × 4 m) per plot, excluding border rows and the first meter at both ends, also with three replicates. The harvested grains were air-dried, weighed, and the yield was adjusted to a standard moisture content of 13%.

The nitrogen use efficiency (NUE) is calculated as follows (Equation 10):

Statistical analyses were performed using SPSS 20.0 (SPSS Inc., Chicago, IL, USA). One-way ANOVA was used to evaluate the effects of year, cultivar, growth stage, and nitrogen treatment on root architecture and enzyme activities. Two-way ANOVA was applied to examine the interactions between cultivar and nitrogen application rate on root architecture, nitrogen accumulation parameters, yield, and nitrogen use efficiency (NUE). Data are presented as means. Structural equation modeling (SEM) was implemented in IBM SPSS AMOS 26 to quantify the relationships among root architecture, nitrogen accumulation, and yield. Chi-square (χ²/df), root mean square error of approximation (RMSEA), standardized root mean square residual (SRMR), and comparative fit index (CFI) were used to evaluate the applicability of SEM (Noll et al., 2022). If an index did not meet the SEM evaluation criteria, the model was adjusted by adding significant relationships or removing non-significant relationships until all indices met the evaluation criteria.

Analysis of variance (Table 1) revealed that cultivar, nitrogen application rate, and their interaction significantly influenced root dry weight density (RDD) in the 0–60 cm soil layer during the R₄~R₆ stages in both growing seasons (p < 0.05). The RDD in the 0–60 cm and 0–20 cm layers exhibited a consistent pattern, initially increasing and then decreasing with growth progression, peaking at the R₆ stage (Figure 2). Significant differences were observed among cultivars and nitrogen treatments (p < 0.05), with the N₁₈₀ treatment consistently exhibiting the highest RDD. Compared to other treatments, the N₁₈₀ treatment significantly increased RDD for both cultivars in both years. The trend in the 0–20 cm soil layer was consistent with that in the 0–60 cm layer.

Analysis of variance (Table 1) indicated that cultivar, nitrogen application rate, and their interaction significantly affected soybean root length density (RLD) in the 0–60 cm soil layer during the R₄~R₆ stages in both growing seasons (p < 0.05). The RLD in both the 0–60 cm and 0–20 cm soil layers showed an initial increase followed by a decrease during the growth period, reaching maximum values at the R₆ stage (Figure 3). Significant differences were observed among different cultivars and nitrogen treatments (p < 0.05), with the N₁₈₀ treatment consistently exhibiting significantly higher RLD than other treatments across both cultivars and years. The changing pattern of RLD in the 0–20 cm soil layer was generally consistent with that in the 0–60 cm layer.

Analysis of variance revealed that cultivar, nitrogen application rate, and their interaction significantly affected soybean root surface area density (RSD) in the 0–60 cm soil layer during the R₄~R₆ stages across both growing seasons (p < 0.05). The RSD in both the 0–60 cm and 0–20 cm soil layers exhibited a pattern of initial increase followed by a decrease during the reproductive growth period, reaching maximum values at the R₆ stage (Figure 4). Significant differences were observed among different cultivars and nitrogen treatments (p < 0.05), with the N180 treatment consistently showing significantly higher RSD than other treatments for both cultivars in both years. The RSD trend in the 0–20 cm soil layer was essentially consistent with that in the 0–60 cm layer.

As shown in Figure 5, during both the 2022 and 2023 growing seasons, the activities of NR, GS, and GOGAT in roots of Xindadou 27 and Xinnongdou 2 showed an initial increase followed by a decrease during the reproductive growth period, peaking at the R₄ stage. Nitrogen application significantly increased the activities of root NR, GS, and GOGAT during the R₄ and R₅ stages (p < 0.05). The trends of enzyme activities in 2023 were generally consistent with those observed in 2022.

During both the 2022 and 2023 growing seasons, the root nitrogen accumulation of both soybean cultivars exhibited a slowfast-slow pattern (Figure 6). The root nitrogen accumulation rate showed a typical single-peak curve, rising first and then falling (Figure 7). In 2022, as the nitrogen application rate increased, the total root nitrogen accumulation (K) and average accumulation rate (Vt) of both cultivars initially increased, and then decreased (Table 2). Nitrogen application significantly prolonged the fast accumulation period (FAP, Δt) of total root nitrogen. Compared to the N₀ treatment, the K value and Vt values of Xindadou 27 in the N₁₂₀, N₁₈₀, and N₂₄₀ treatments increased significantly by 42.55%~75.32% and 90.48%~214.29%, respectively; while those of Xinnongdou 2 increased by 45.37%~150.24% and 81.58%~213.16%, respectively.

Over the 2022 and 2023 growing seasons, aboveground nitrogen accumulation in both soybean cultivars followed a typical sigmoidal curve (Figure 8). The aboveground nitrogen accumulation rate showed a typical single-peak curve, rising first and then falling (Figure 9). As shown in Table 3, key accumulation parameters increased with the application of nitrogen. Notably, the N₁₈₀ treatment significantly enhanced the final accumulation potential (K), prolonging the rapid accumulation duration (Δt) by 23.97%~28.58% and substantially increasing the average accumulation rate (Vt) compared to the N0 control. This trend was consistent across both years, with the N₁₈₀ rate consistently promoting the most favorable nitrogen uptake dynamics for each cultivar.

The ANOVA results (Table 4) indicated that both cultivar and nitrogen application rate significantly affected pods per plant, seeds per plant, 100-seed weight, yield, and nitrogen use efficiency (NUE). Yield consistently reached its peak under the N₁₈₀ treatment for both cultivars across both years. Compared to N₀, the yield increases for Xindadou 27 were 19.32% in 2022 and 25.57% in 2023, and for Xinnongdou 2, they were 24.92% in 2022 and 21.81% in 2023. This yield advantage at N₁₈₀ was associated with a significant increase in both the number of pods per plant and the number of seeds per plant. Notably, while 100-seed weight often reached its maximum under the N₂₄₀ treatment, this increase in individual seed mass did not lead to higher overall yield beyond the optimum achieved at N₁₈₀. For NUE, both cultivars reached their maximum under N₁₈₀ in 2022. In 2023, the highest values for Xindadou 27 and Xinnongdou 2 were 4.13 kg kg^-1 under N₁₈₀ and 4.88 kg kg^-1 under N₁₂₀, respectively.

As shown in Table 5, all fit indices were within the acceptable range, with χ²/df=0.649, RMSEA = 0.029, SRMR = 0.031, and CFI = 1.000. All fit parameters of the initial model indicated a good fit. Structural equation modeling (SEM) was applied to evaluate the effects of root architecture and nitrogen accumulation on yield (Figure 10). Root dry weight density (RDD), root length density (RLD), and root surface area density (RSD) significantly and positively influenced root nitrogen accumulation, with path coefficients of 0.887, 0.205, and 0.198, respectively. Root nitrogen accumulation demonstrated a highly significant positive effect on aboveground nitrogen accumulation. Furthermore, aboveground nitrogen accumulation exhibited a significant positive effect on yield, with a path coefficient of 0.509.