Phenological Development, Environmental Effects, and Selection of Soybean Lines from Crosses with Contrasting Relative Maturity Groups

1. Introduction

Soybean (Glycine max (L.) Merrill) is the main oilseed crop cultivated worldwide, with a major impact on protein and vegetable oil supply, being important for human and animal nutrition as well as bioenergy production [1]. In the 2025/2026 growing season in Brazil, soybean production is estimated at 180.13 million tons [2], maintaining the country as the world’s largest producer. Because of its economic importance and wide cultivation across diverse environments, agronomic traits related to the phenological cycle, such as number of days to flowering and maturity, have gained prominence in plant breeding programs. Genetic advances in soybean have been accompanied by changes in cycle duration, resulting in the development of more productive cultivars adapted to environmental conditions [3]. For this purpose, it is essential to classify cultivars according to their relative maturity group (RMG), allowing their proper allocation to different growing environments, since phenological maturity is strongly influenced by environmental factors such as latitude, altitude, and climate [4] in each region.

RMG is a widely used classification to describe the phenological cycle of soybean cultivars, based on plant responses to photoperiod and temperature, which are key regulatory factors for flowering and maturity [4]. The RMG system is an adaptation of the North American system, grouping genotypes into 13 classes (000, 00, 0, I, II, III, IV, V, VI, VII, VIII, IX, and X) representing adaptation ranges, from extremely early materials (group 000, high latitudes) to the latest ones (group X, latitude 0°) [5]. In Brazil, this classification was converted into a continuous numerical scale ranging from 5.0 (early cultivars) to 9.0 (late cultivars), facilitating regional cultivar placement and breeding cross planning [4]. In this context, RMG is also essential for defining sowing windows and suitable growing environments for each cultivar, since each adaptation zone is expected to meet the developmental requirements of specific maturity groups. However, when cultivars are grown outside the optimal sowing window or adaptation zone, their phenological cycle may be shortened or extended beyond the expected pattern [6].

From a biological perspective, this classification is associated with genetic differences underlying soybean phenological development. According to previous studies [7,8,9], genetic advances have enabled the identification of important loci and regulatory pathways associated with soybean flowering and maturity, including E-series genes, the J gene, circadian clock components, FT-like genes, and several QTLs. These loci influence the transition from vegetative to reproductive development, the duration of the reproductive period, and the time required to reach maturity. Thus, differences among genotypes in RMG represent not only agronomic variation, but also differences in the genetic control of phenology under different environmental conditions.

In addition to genetic factors inherent to each RMG, phenotypic expression is strongly influenced by environmental variables. In soybean, the main environmental effects related to changes in phenological behavior are temperature, photoperiod, and precipitation [10], as well as soil type and agricultural practices [11], also affecting yield and seed composition-related traits [12]. Considering these complex interactions, appropriate statistical methodologies are essential. Mixed models via REML/BLUP are particularly useful in plant breeding because they can handle unbalanced data, estimate variance components, predict genetic values with high accuracy, and incorporate environmental covariates into the analysis [13].

Despite these advances, there are still gaps in the understanding of soybean phenological behavior regarding the interaction between relative maturity groups and environmental variables. In this context, the objective of the present study was to: (i) characterize the response of soybean populations in terms of the duration of phenological periods across different growing seasons derived from crosses contrasting in RMG; (ii) select lines with high performance for yield and earliness using REML/BLUP in each population; and (iii) evaluate the effect of precipitation and temperature on the phenological cycle of these populations.

2. Materials and Methods

2.1. Experimental Area Characterization

The experiments were conducted at the Teaching, Research and Extension Farm (FEPE) of FCAV – UNESP, Jaboticabal Campus. The municipality is located in the northwestern region of the state of São Paulo, at latitude 21º 14’ 59” S and longitude 48º 17’ 8” W, at an altitude of 575 meters. The soil in the area is classified as a eutrophic dark-red Latosol, with a very clayey texture and gently undulating relief. According to the Köppen-Geiger classification, the climate of the region is type Cwa [14].

2.2. Genetic Material

Two populations were synthesized from narrow biparental crosses in relation to the Relative Maturity Group (RMG) of the cultivars used as parents. The combinations used to obtain the populations were as follows: BMX Potência RR (RMG 6.7) (♀) × BMX Energia RR (RMG 5.3) (♂), forming population E2, derived from parents with early and intermediate maturity cycles; and BRS 245 RR (RMG 7.3) (♀) × BRS 278 RR (RMG 9.4) (♂), forming population E3 [15], derived from parents with intermediate and late maturity cycles.

2.3. Experimental Design and Treatments

The experimental design used was an augmented block design [16], in which the lines (treatments) were arranged in single 5-meter rows, with 0.5 m spacing between rows (2.5 m² of useful area) and were included only once in the experiment without replication. Three check cultivars were interspersed within the blocks, each repeated twice per block, totaling 30 check plots across the five blocks. At least one of these checks consisted of one of the parental lines.

The checks used were: BMX Potência RR (6.7), AS 3680 IPRO (6.8), and M 6410 (6.4) for population E2; and BRS 245 RR (7.3), HO APORE IPRO (7.5), and BS 2606 IPRO (6.0) for population E3. The selection of checks was based on their maturity group proximity to that of the parental lines of the populations.

The populations were evaluated during the growing seasons 2020/2021, 2021/2022, 2022/2023, and 2023/2024. Population E2 (F7:10) consisted of 150 lines, while population E3 (F7:10) consisted of 80 lines, all at advanced generations of inbreeding.

2.4. Management of the Experimental Area

The management practices were carried out according to Embrapa recommendations [17]. The crops were established under a no-tillage system (NTS), with a sowing density of 18 plants per meter. Soil fertility was evaluated annually to ensure appropriate application of amendments and fertilizers. Seeds were inoculated with bacteria of the genus Bradyrhizobium to supply nitrogen. Pest, disease, and weed management were conducted as needed.

2.5. Evaluated Traits

2.5.1. Agronomic Traits

The evaluated traits were the number of days to flowering (NDF), number of reproductive days (NRD), number of days to maturity (NDM), and grain yield (GY). The NDF trait was assessed at the R1 phenological stage when plants transition from the vegetative to the reproductive phase. The NRD trait was calculated as the number of days from R1 to R8, and NDM was assessed at the R8 developmental stage, which corresponds to the period from sowing to the point when at least 50% of the plants have 95% or more mature pods, according to the classification of Fehr and Caviness [18] and is expressed in days. GY was obtained from the grain weight of the plot useful area (2.5 m²) after harvest and threshing. The data obtained in grams per plot were converted to kilograms per hectare (kg ha⁻¹) and adjusted to 13% moisture content.

2.5.2. Environmental Variables

In addition to the phenotypic evaluation of the traits, environmental variables associated with the lines were also assessed, namely accumulated precipitation (mm) and mean temperature (°C) during the cropping cycle. These variables were calculated individually for each line, considering the duration of each studied phenological period, to represent the climatic conditions experienced by each line throughout its development. This approach was adopted because fixed calendar windows could include periods in which earlier lines had already completed the evaluated phenological stage, potentially misrepresenting their environmental exposure. The meteorological data used in this study were provided by the Agroclimatological Station of the Department of Exact Sciences at FCAV/UNESP, Jaboticabal Campus.

2.6. Statistical Analyses

To visualize the behavior of NDF, NRD, NDM, and GY across growing seasons, boxplot graphs were used, as they allow comparison among cropping years, detection of patterns, and identification of differences related to line performance.

The linear mixed model used follows the structure presented by Bates et al. [19], defined by the conditional distribution of the response variable Y given the vector of random effects B:

$\left(Y\right|B=b)~N(X\beta +Zb,{\sigma }^2{W}^{-1})$where Xβ represents fixed effects, Zb random effects, σ² is the scale parameter (residual variance), and W is a diagonal matrix of known weights (assumed to be the identity).

For each experiment, the data were fitted using linear mixed models considering the structure of the Augmented Block Design (ABD), with genetic effects treated as random, using the lme4breeding package [20]. The model used was as follows:

$y_{ijkl}=\mu +a_i+c+g_j+{\left(gxa\right)}_{ij}+b_{k\left(i\right)}+{\epsilon }_{ijkl}$where:

$y_{ijkl}$: observed value (NDM or GY)

$\mu$: overall mean of the trait

$a_i$: fixed effect of growing season $i$$c$: fixed effect of checks

$g_j$: random effect of genotype $j$${\left(gxa\right)}_{ij}$: random effect of genotype x growing season interaction

$b_{k\left(i\right)}$: random effect of block

${\epsilon }_{ijkl}$: experimental error.

The quality of the models was evaluated using the coefficient of determination (R²), with a decomposition of these coefficients into fixed and random effects, estimated using the rsq package. To assess the significance of fixed effects in the mixed model, a Wald-based ANOVA (Type III) was performed.

Residual diagnostics were performed for all fitted mixed models to assess normality and homogeneity of variance. For grain yield (GY), number of days to flowering (NDF), and number of reproductive days (NRD), residuals did not show relevant deviations from normality or homoscedasticity. For number of days to maturity (NDM), some deviation from these assumptions was observed; however, since variance components and fixed-effect estimates were obtained via restricted maximum likelihood (REML) within a linear mixed-model framework, which is comparatively robust to moderate departures from normality and homogeneity of variance, this deviation was not considered to materially compromise the reliability of the genetic parameter estimates or the significance tests reported for this trait.

The variance components (σ_g², σ_gxe², σ_e²) were estimated by restricted maximum likelihood (REML). The BLUPs (Best Linear Unbiased Predictors) were extracted from the estimated random effects in the mixed models. Two-dimensional scatter plots were generated between NDM and GY BLUPs to identify dual-purpose lines for earliness and higher productivity in each population.

The broad-sense heritability (H²) was estimated from the variance components for each population, considering the different years. According to the equation:

$H^2={\displaystyle \frac{{\sigma }_g^2}{{\sigma }_g^2+\frac{{\sigma }_{gxa}^2}{a}+\frac{{\sigma }_e^2}{ar}}}$where a is the number of environments, corresponding to the evaluated growing seasons, and r is the number of within-environment replications. In the augmented block design used in this study, test lines were unreplicated within each growing season, therefore, r was considered equal to one for these lines. The replicated checks and blocks were used to adjust local environmental variation and improve experimental precision. The selective accuracy of the BLUPs was estimated as $r_{\widehat{g}g}=\sqrt{H²}$. To assess experimental precision, the experimental coefficient of variation ${(CV}_e)$ was calculated, as well as the ratio between the genetic coefficient of variation and the experimental coefficient of variation $\left({\displaystyle \frac{{CV}_g}{{CV}_e}}\right)$.

Additionally, mixed models were fitted for NDF, NRD, and NDM, using environmental covariates as fixed effects to evaluate their influence on the phenotypic response of the populations, using the lme4 package [19]. Preliminary analyses indicated collinearity between accumulated precipitation and mean temperature, so these covariates were fitted in separate single-covariate models to avoid unstable coefficient estimates. In cases of low variation, the random block effect was not considered. The model used was as follows:

$y_{ijk}=\mu +{\beta }_1{x}_i+c+{\left(gxa\right)}_{ij}+b_{k\left(i\right)}+{\epsilon }_{ijk}$where:

$y_{ijk}$: observed value (NDF, NRD and NDM)

$\mu$: overall mean of the variable

${\beta }_1{x}_i$: fixed effect of the environmental covariate $x_i$ (accumulated precipitation or mean temperature)

$c$: fixed effect of checks

${\left(gxa\right)}_{ij}$: random effect of the genotype x growing season interaction

$b_{k\left(i\right)}$: random effect of block

${\epsilon }_{ijk}$: experimental error.

Marginal prediction plots were generated, allowing visualization of the effect of environmental covariates on NDF, NRD, and NDM, using the ggeffects package. Model quality was evaluated using the coefficient of determination (R²), both marginal and conditional, estimated using the performance package, where marginal R² refers to the variance explained only by fixed effects and conditional R² represents the total variance explained by the model. The β coefficients were extracted from the fixed effects estimated by the mixed model. To assess the significance of fixed effects in the mixed model, a Wald-based ANOVA (Type III) was performed. All analyses and plots were generated using R software [21].

3. Results

The distributions of accumulated precipitation (mm) and temperatures (°C) differed among the evaluated growing seasons. The 2020–2021 and 2021–2022 seasons showed similar temperature and precipitation patterns, whereas the 2022–2023 season was the coolest and had the highest rainfall intensity. In contrast, the 2023–2024 season presented the highest temperatures and the lowest precipitation levels (Figure A1).

The traits NDF, NRD, NDM, and GY of Population E2 exhibited different behaviors across the growing seasons (Figure 1). For NDF in Population E2, similar patterns were observed in the 2020-2021 and 2022-2023 seasons, with plants flowering within a similar range of days. In the 2021-2022 season, the plants exhibited delayed flowering, while in the 2023-2024 season, flowering occurred earlier. For NRD and NDM in Population E2, similar patterns to NDF were observed across all growing seasons, with the duration of the reproductive period and time to maturity showing consistent behavior across the different years. Regarding GY, a progressive decrease was observed across the growing seasons (Figure 1).

Population E3 exhibits later flowering compared to Population E2, in which the influence of the growing seasons on this trait is also observed. In the first three seasons, similar NDF patterns were noted, while in the 2023-2024 season, an earlier flowering can be observed (Figure 2). Regarding NRD, it can be observed that in the 2023-2024 season, when flowering occurred earlier, the number of reproductive days of the plants was higher than in previous seasons. The same NRD pattern can be observed in NDM. As in Population E2, GY in Population E3 showed a progressive decrease across the seasons (Figure 2).

For the population’s behavior, the growing season had a significant effect in all evaluated years, in some cases, the checks did not show a significant effect (Table 1). The coefficients of determination of each model were generally high, except for GY, with alternating importance between fixed and random effects across models. For NDF in populations E2 and E3, the variation was mostly accounted for by the fixed effects (0.65 and 0.76, respectively). For NRD, in population E2, the variation was also primarily due to the fixed effects (0.43), whereas in population E3, it was mainly attributable to the random effects (0.68). For NDM in populations E2 and E3, the random factors represented the main component of the variation (0.69 and 0.56, respectively). For GY, in the models for populations E2 and E3, the explanatory power was relatively low, but most of the variation was still due to the random effects (0.32 and 0.36, respectively) (Table 1).

In Population E2, the genotypic variance for NDF, NRD, and NDM was higher than the other variance components (2.65, 4.80, and 9.41, respectively), except for GY, which showed a high residual variance. The same pattern was observed in Population E3 for NDF, NRD, and NDM (10.21, 15.84, and 21.23), as well as a high residual variance for GY. The genotype × year interaction variance showed changes depending on the population and the trait analyzed, while the block variance was relatively low (Table 2).

NDF and NDM exhibited high heritabilities, NRD showed intermediate heritability, and GY presented lower heritability, especially in Population E3. The CVg/CVe ratio was greater than 1 for NDF, NRD, and NDM in both populations, indicating potential for selection on these traits, as well as the selective accuracy of the BLUPs. Regarding the means, Population E2 was characterized as earlier than Population E3 for both NDF and NDM. For NRD, Population E2 had a higher mean than Population E3. Concerning GY, Population E2 (3,553.15 kg ha⁻¹) showed better performance than Population E3 (2,953.71 kg ha⁻¹) (Table 2).

${\sigma }_g^2$: genotypic variance; ${\sigma }_{gxy}^2$: genotype × year interaction variance; ${\sigma }_{b\left(y\right)}^2$: variance of blocks within years; ${\sigma }_e^2$: residual variance; $H²$: broad-sense heritability; $r_{\widehat{g}g}$: selective accuracy of BLUPs; ${CV}_e$: experimental coefficient of variation; ${CV}_g$: genotypic coefficient of variation.

Figure 3 and Figure 4 show a dispersion of the BLUPs estimated for NDM and GY. The zero (0) point on the axes represents the overall mean for each trait. In this context, positive BLUPs for NDM indicate late lines, while negative BLUPs indicate early lines, based on the population mean. The same applies to GY, where lines with positive BLUPs are more productive and lines with negative BLUPs are less productive.

In Population E2 (Figure 3), BLUPs for NDM ranged from –10 to +5 days, representing a difference of 15 days between the earliest and latest lines. For GY, BLUPs ranged from –1089 to +700 kg ha⁻¹. A complete list of inbreed lines with their BLUP values for NDM and GY is provided in the Appendix A (Table A1).

For Population E3 (Figure 4), BLUPs for NDM ranged from –9 to +9 days, representing a difference of 18 days between the earliest and latest lines in the cycle. For GY, BLUPs ranged from –550 to +814 kg ha⁻¹. Full list of inbred lines is in Appendix A (Table A2).

To visualize the influence of the environment on the variation of NDF, NRD, and NDM, mixed models were fitted for each population using the environmental covariates accumulated precipitation and mean temperature. The models exhibited high explanatory power, with conditional R² ranging from 0.79 to 0.97 (Figure 5 and Figure 6).

The environmental covariates were significant (p < 0.001) in all models fitted for Population E2 (Table A3). Lower precipitation levels were associated with early flowering, while increased precipitation was associated with delayed flowering. The model exhibited high explanatory power (R² = 0.91), with a considerable contribution from accumulated precipitation (R² = 0.56), presenting a positive β coefficient, that is, for each 10 mm increase in precipitation, NDF increased by 0.42 days (Figure 5A, Table A3).

Lower temperatures are associated with delayed flowering, while higher temperatures are associated with earlier flowering. The model showed an explanatory power of 0.87, with mean temperature as a fixed effect explaining 0.30. The β coefficient was negative; that is, for each 1 °C increase in mean temperature, NDF decreases by 2.29 days (Figure 5B, Table A3).

Regarding NRD, the pattern is repeated: increased precipitation is associated with a higher number of reproductive days, while increased mean temperature is associated with a reduction in NRD. The models showed explanatory powers of 0.81 and 0.79, respectively, and a low contribution from the fixed effects (0.13 and 0.22, respectively), with a β coefficient of 0.052 for accumulated precipitation and –3.50 for mean temperature (Figure 5C–D, Table A3).

For NDM, that is, considering the entire cycle of the lines, the environmental covariates had significant effects, although their contributions were smaller than those of the random effects, with R² values of 0.18 for accumulated precipitation and 0.32 for mean temperature. The full models exhibited high explanatory power (0.94 and 0.93, respectively). For NDM, each 10 mm increase in precipitation resulted in an increase of 0.101 days in the total cycle, while each 1 °C increase led to a reduction of 5.56 days in the total cycle of the lines (Figure 5E–F, Table A3).

In Population E3, for NDF, increased precipitation delayed flowering of the lines, while increased temperature advanced flowering of the population. The models showed good fit, with conditional R² of 0.97 for accumulated precipitation and 0.96 for mean temperature, such that the fixed effects (0.65 and 0.68, respectively) were highly important in explaining the variation in flowering. The β coefficient was 0.169 for accumulated precipitation and –2.587 for mean temperature (Figure 6A–B, Table A3).

In contrast to what was observed for Population E2, regarding NRD, in Population E3 increased precipitation was associated with a shortening of the reproductive days, while increased temperature was associated with a prolongation of these days. The models showed high R² values (0.89 and 0.91, respectively), but with a low contribution from the fixed effects (0.11 and 0.21, respectively). The β coefficient was negative for precipitation (–0.068) and positive for mean temperature (4.867) (Figure 6C–D, Table A3).

For NDM, both models showed a conditional R² of 0.95. However, the fixed effects were not significant (p = 0.783 for accumulated precipitation and p = 0.726 for mean temperature), having no influence on the shortening or prolongation of the total cycle of the lines (Figure 6E–F, Table 3A). The decomposition of the coefficients of determination into fixed and random effects is available in the Appendix A (Table A4).

4. Discussion

The variability observed (Figure 1 and Figure 2) for NDF, NRD, NDM, and GY in the populations can be attributed to a combination of genetic and environmental factors: genetic factors associated with the maturity group (RMG) of each line and, consequently, adaptation to the location where the experiments were conducted; and environmental factors related to the climatic conditions of each growing season (Figure A1). Due to being a less contrasting cross in terms of RMG, population E2 presented lower genotypic variance for NDF, NRD, and NDM, characterizing a more early-maturing population, inheriting this trait from its parents. In contrast, population E3 showed higher genotypic variance for NDF, NRD, and NDM, largely due to the greater genetic distance between its parents (Table 2). This relationship between parental genetic distance and progeny variability was demonstrated by Jean et al. [22], as crosses between genetically distant parents generate a wide distribution of variance among the resulting progenies, whereas crosses between closely related parents generate lower variance among progenies for the trait in question.

Regarding the check cultivars, a significant effect (p < 0.001) was detected for NDF and NDM in population E3 (Table 1), whereas no significant check effect was observed in population E2 for any trait. This result likely reflects genetic differences among the three check cultivars used in population E3 (BRS 245 RR, HO APORE IPRO, and BS 2606 IPRO), which span a wider range of relative maturity groups (6.0 to 7.5) compared with the checks used in population E2. Because check performance was used to adjust for local environmental variation within blocks, this significant effect indicates that the checks captured genuine differences in phenological behavior rather than purely environmental noise, reinforcing the validity of using them as a reference for adjusting the BLUPs of the unreplicated test lines in this population.

The growing seasons had a strong influence on the phenological variables, as well as on GY of the populations, since in each growing season the same genotypes showed different performances (Figure 1 and Figure 2, Table 1, and Figure A1). This can be explained by the environmental conditions inherent to each growing season, which can affect plant phenological development as well as yield performance. Gong et al. [23] observed that changes in precipitation, temperature, and sunlight hours can shorten or extending the soybean phenological cycle. Soybean may exhibit positive or negative phenotypic plasticity depending on genotype × environment interaction. Generally, under stress conditions, as observed across seasons with reduced precipitation and increased temperature (Figure A1), plants may exhibit negative plasticity, reducing some traits such as yield components [24].

High broad-sense heritability estimates for NDF, ranging from 0.77 to 0.90 in both populations, demonstrate the strong genetic control of this trait, as also reported in the literature: 0.82 [25], 0.88 [26], and 0.94 [27]. In addition, flowering has a strong genetic component, although environmental heterogeneity is a crucial factor contributing to variation in flowering time in soybean [28], as also observed in our results. NRD in both populations showed intermediate heritability (0.68 and 0.66) and higher variance compared to NDF, demonstrating that the environment influences the duration of the reproductive phase, i.e., redistributing time among phases [29]. NDM (Table 2), ranging from 0.82 to 0.84, showed that despite environmental influence on the crop cycle, the genetic component plays a predominant role in the phenotypic expression of the trait, a fact widely reported in the literature. Previous studies have reported similar or higher heritability values, with estimates of 0.72 [26], 0.85 [30], 0.91 [31], 0.93 [32], and up to 0.94 [27].

For GY, the largest variance component observed was the residual, indicating that this trait is strongly influenced by environmental conditions [33] and other factors not included in the model, resulting in lower heritability values. Heritability estimates for GY were lower compared to phenological traits, ranging from 0.43 to 0.62 in both populations, which is expected for complex traits such as GY (Table 2). These results are consistent with values reported in other studies, which estimated heritability for yield at 0.59 [26] and 0.62 [30].

The use of mixed models through BLUPs in several breeding programs has proven to be an efficient tool for the prediction of genetic values and selection support, especially for quantitative traits influenced by the environment [34,35,36,37]. The distribution of BLUPs (Figure 3 and Figure 4) highlights the existence of genetic variability for NDM and GY within each population, which was also observed by Bianchi et al. [27]. This variability allows the identification and selection of lines with superior performance for earliness and productivity based on their predicted genetic value. Thus, as demonstrated in the present study, the BLUP-based approach enables the simultaneous selection of dual-purpose genotypes, considering phenological and productive traits in an integrated manner [38].

In population E2 (Figure 5A–F), accumulated precipitation showed a positive relationship, while mean temperature showed a negative relationship with the duration of phenological stages in both growing periods. In years with higher precipitation, an extension of the phenological cycle was observed, whereas drier years were associated with a reduction in cycle duration. Similarly, higher temperatures were associated with a shortening of the cycle, while cooler conditions resulted in longer duration. These results are consistent with the literature, which indicates that increased precipitation tends to prolong crop cycles, while increased temperature accelerates development and reduces the growth period [10,24,39]. This behavior is attributed to the fact that higher temperatures increase plant development rates, shortening phenological stages [40].

Liu and Dai [39], when evaluating the phenological response of soybean under different precipitation and temperature conditions, observed distinct effects depending on the developmental stage. In the vegetative period, variations of −0.01 to −0.11 days/10 mm of precipitation and −0.63 to −2.51 days/°C of mean temperature were reported. In the reproductive period, the effects ranged from −0.02 to 0.11 days/10 mm and −0.98 to 0.56 days/°C, indicating greater variability in the climatic response during this phase. For the total crop cycle, precipitation showed effects ranging from −0.06 to 0.09 days/10 mm, while temperature reduced the cycle duration by approximately 2.51 to 5.19 days/°C, depending on the region and growing season.

The β coefficients estimated in the present study for precipitation showed positive relationships with flowering (0.42 days/10 mm), reproductive period (0.052 days/10 mm), and maturity (0.101 days/10 mm), suggesting that increased precipitation tends to slightly prolong the duration of phenological stages. In contrast, mean temperature showed a consistent negative effect throughout the crop cycle, with a less pronounced reduction at flowering (−2.289 days/°C), followed by the reproductive period (−3.502 days/°C) and maturity (−5.565 days/°C), indicating that increased temperature progressively accelerates crop development (Figure 5A–F).

For population E3, the same pattern observed in population E2 for NDF was maintained regarding the effects of precipitation and temperature. For NRD, a shortening of the period was observed with increasing precipitation, and an extension of the period with increasing temperature, although the contribution of these variables was low (R² = 0.11 and 0.21, respectively; Figure 6C–D). This pattern is consistent with a possible self-compensation mechanism between phenological phases within the genotype itself; however, since no physiological measurements (e.g., photoassimilate partitioning or developmental rate at the sub-phase level) were taken in this study, this interpretation remains a hypothesis and warrants direct investigation in future research. For NDM, environmental covariates did not show a significant effect, suggesting that later-maturing genotypes are less sensitive to precipitation and temperature effects, maintaining a more stable phenological pattern across years (Figure 6A–F). This result is consistent with that reported by Sobko et al. [12], who found that early-maturing soybean genotypes are more sensitive to precipitation and temperature than later-maturing genotypes.

Both results obtained for populations E2 and E3 corroborate the study by Xin et al. [41], which demonstrates that soybean phenology is influenced by climatic conditions, with temperature generally exerting a negative effect on the duration of phenological stages, indicating an acceleration of crop development. And the precipitation shows a predominantly positive effect, being associated with the prolongation of phenological stages and the overall development cycle.

These results indicate that the inclusion of environmental variables in models becomes important for understanding how certain factors influence agronomic traits, allowing the quantification of specific effects on plant development. However, careful attention must be given to the selection of variables, as they should have a causal or physiological relationship with the process being modeled, and additionally environmental variables are generally collinear, which can make it difficult to distinguish the true contributions of different climatic factors [10].

Although this study provides relevant evidence on the effect of the environment on the phenological and productive modulation of biparental populations with contrasting maturity patterns, some limitations should be acknowledged. First, the experiments were conducted at a single location, although over four growing seasons. Therefore, although seasonal climatic variation was captured, the extrapolation of the results to other soybean-producing regions with different latitude, soil, altitude, and management conditions should be made with caution. Second, the study was based on only two biparental populations, which limits the generalization of the observed responses to the broader genetic diversity of soybean regarding relative maturity groups. Third, considering the large number of lines evaluated, the use of an augmented block design was an appropriate choice; however, this design did not allow replication of the lines within each growing season, which could have provided greater experimental precision.

An additional limitation concerns the plot structure used for grain yield estimation. Test lines were evaluated in single, unreplicated 5-m rows without border rows separating adjacent plots. While this configuration is commonly accepted for phenological traits such as NDF, NDM, and NRD, which are less susceptible to inter-plot interference and grain yield is known to be more sensitive to competition and border effects between neighboring plots, particularly in single-row, non-bordered designs. Shading, differential resource competition, and edge effects from adjacent genotypes with contrasting growth habits or maturity could have contributed to the high residual variance and comparatively low heritability observed for GY (Table 2), independently of genuine environmental effects. Therefore, part of the unexplained variance attributed to environmental influence may also reflect this methodological source of noise rather than climatic variation alone. Future evaluations aimed specifically at yield estimation would benefit from replicated plots with border rows or alternative spatial correction methods to better isolate genetic and environmental effects on this trait.

In addition, the environmental-covariate models focused on accumulated precipitation and mean temperature, although other factors relevant to the context of the study may also influence population responses, such as solar radiation, soil water availability, growing degree days, and biotic stresses. Environmental variables are often collinear, which may make it difficult to distinguish the true contributions of different climatic factors [10]. Because precipitation and temperature were evaluated separately due to collinearity, the estimated effects should be interpreted as associations with phenological responses, rather than as fully independent causal effects. Therefore, future studies including additional locations, greater genotype diversity with a broader range of RMG, and the use of replicated experimental designs would be useful to confirm the applicability of the results observed in this study.

5. Conclusions

Population E2 showed lower variability among lines and a predominance of early-maturing genotypes, characterized by uniformity, while population E3 also exhibited uniformity, but with predominantly late-maturing lines, resulting from parents that were not highly contrasting in terms of relative maturity group.

The use of mixed models via REML/BLUP was efficient for predicting genetic values and enabled the identification of lines with potentially higher productivity, as well as the identification of early-maturing lines.

Population E2 was more sensitive to environmental variations in precipitation and temperature, whereas population E3 showed a lower influence of these covariates on the lines cycle, suggesting that early-maturing genotypes exhibit greater phenotypic plasticity, while late-maturing genotypes tend to present greater phenological stability, compensating between vegetative and reproductive periods and resulting in a more stable total cycle duration.