Genetic and environmental factors regulating soybean reproductive stages and their transitions

complex trait, reproductive maturity, genetic effect , photoperiod, thermal 8 accumulation, maturity model 9 10 11 12 13 14 15 16 17 18 19 20 21

22 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 2 The reproductive stage of soybean is influence d by the effect of genotype, 23 environment, and their interactions . While days to flowering and days to full maturity 24 have been studied, a systematic and comprehensive study that investigates the 25 variation in days to each stage and the role of maturity-related genes and environmental 26 variables is lacking. Therefore, we s tudied 508 unique accessions from the USDA 27 germplasm collection from maturity group 0-IV, and a set of 67 near-isogenic lines 28 differing for maturity-related genes. Field experiments and evaluations were conducted 29 in central Iowa , USA . T he days to each of t he reproductive stages, R1 -R8, were 30 recorded. We report considerable variation in the duration of reproductive growth stages 31 between flowering and maturity, which is largely explainable by known flowering and 32 maturity genes as well as environmental variabl es, day length , and growing degree 33 days. Besides the known maturity -related genes E1, E2, and Dt1, w e identified two 34 novel SNPs, such as Glyma.01G180600 and Glyma.10G221300, as potential targets for 35 genetic regulation of reproductive stages. We also captured two other loci, 36 Glyma.08G216800 and Glyma.04G088100 for day length and growing degree days, 37 respectively, that revealed dynamic regulation of environmental gradients on the 38 reproductive stage s. Furthermore, we developed a random forest -based gene tic 39 maturity model that can predict genetic and environmental effects across a wide range 40 of genotypes. This study broadens the understanding of the factors that contribute to 41 reproductive development, which will help to develop cultivars that combine the optimal 42 combinations of stage durations for a higher seed yield and enhanced resilience. 43 44 45 46 47 48 49 1 INTRODUCTION 50 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 3 As with many crop species, soybean phenology is commonly characterized 51 based on easily measurable morphological characteristics using the staging code 52 developed by Fehr and Caviness 1977 (1). During the vegetative stages, soybean is 53 characterized by the number of nodes with unfurled leaves, which develop from the first 54 node above the cotyledons and up. In the vegetative period, soybeans grow very rapidly 55 and devote their energy to leaf, stem, and root tissue development. The initiation of 56 flowering signals indicates the transition to the reproductive stage. Eight reproductive 57 stages (R1 to R8) are generally recognized: two stages of flower initiation, two stages of 58 pod formation and elongation, two stages of seed formation and expansion, and two 59 stages of physiological maturity (Figure 1A). The growth and development aspects are 60 directly related to the farmer ’s decisions for the choice of cultivar, time of planting , and 61 days to maturity for harvesting the crop. The decision -making process is further 62 complicated by the weather (2). Stresses and injuries during the reproductive phases 63 can cause significant seed yield loss, as time and resources for compensatory growth 64 are limited (3,4). 65 One of the most important traits in soybean breeding and production is seed yield (5). 66 Manipulation of the soybean developmental timeline holds potential for yield 67 improvement through different mechanisms. Ideotype breeding attempts to optimize the 68 combination of traits to achieve higher yield through physiological and morphological 69 changes (6). The timing of reproductive stages may be altered to achieve a more 70 optimal source -sink relationship and to minimize source limitation. For example, a 71 prolonged flowering period may result in greater pod set, as fewer flowers would be 72 competing for resources at any time (7). The developmental stages and their duration 73 can be altered to avoid abiotic and biotic stresses (8). For example, if biotic or abiotic 74 stresses are likely at a specific time of year, such as due to weather patterns or insect 75 reproduction cycles, one way to mitigate or lessen the negative impact is to alter the 76 timing of developmental stages to minimize the risk of these stresses coinciding with the 77 most susceptible stages. Another farming strategy is to manipulate the planting date, 78 which has cascading effects on soybean phenology. Among the reproductive stages, R1 79 (Flowering), R4 ( Pod initiation), and R8 (physiological maturity) are more commonly 80 studied (1,9); however, there is a lack of literature on comprehensive reproductive stage 81 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 4 studies. To optimize soybean maturity and production, it is essential to comprehend how 82 genetics and the environmental variables regulate reproductive development. 83 Understanding the mechanism underlying phenotypic variation requires determining the 84 extent to which populations exhibit underlying genetic variation, and the effect of the 85 environmental variables in natural field conditions (10,11). However, it is challenging to 86 predict responses to variable climates without knowing the particular environmental 87 cues that cause plastic responses, including in more extreme and less studied climates 88 (12). Early planting in the U.S. Midwest region is a recommended management strategy 89 to boost soybean yield (13), but personnel, equipment, logistics, and environmental 90 issues might cause planting to be delayed. Soybeans are more vulnerable to spring 91 killing frost, early -season insects, seedling disease, and damaging rainfall events that 92 could lead to a substandard emergence when early planting is feasible (14). In the MGII 93 and earlier maturities, the risk of frost rises and the time to emergence falls with delayed 94 planting. Diseases of the seedlings and roots are less likely to occur when planting and 95 emergence times are shorter (2). Therefore, studies on soybean maturities, the effect of 96 known maturity -related genes, the usefulness of identifying new genetic factors, the 97 relationship of genetics with the reproductive stages, and their transitions are crucial. 98 The understanding of the genetic underpinnings of development al stage timing is 99 further complicated by both temperature and photoperiod, as they are known to affect 100 the timing of reproductive stages (15–18). The findings of a simulation study by 101 Forecast and Assessment of Cropping sysTemS (FACTS) indicated that te mperature 102 affects early vegetative development (days from planting to blooming) more than 103 photoperiod, and this effect is especially noticeable in the maturity group (MGII). 104 Temperature was significant in the maturity group MGIV, but the effect of photoper iod 105 (influenced by latitude) was greater in the MGII (2). The genetic control of flowering (R1) 106 and maturity (R8) timing has been extensively studied in soybean due to its effect on 107 adaptation to a specific environment (19–22). At least 14 genetic loc i have been 108 identified that have a clear effect on the latitudinal adaptability of a given soybean 109 variety in the field ( Supplementary File 1 - Table S1). These include eleven E genes 110 (E1- E11) and three additional genes (Dt1, Dt2, J), which cause major shifts in latitudinal 111 adaptability (between 5-30 days to maturity per E locus) (23,24). E1, a B3 transcription 112 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 5 factor, epistatically controls time to maturity through its interactions with several other E 113 genes (25). In the presence of epistasis, the allelic effect of E1 to e1 delays maturity by 114 nearly three weeks. E2 is an ortholog of the Arabidopsis gene GIGANTEA (26). E3 and 115 E4 are phytochrome A genes (27,28). Dt1 is an ortholog of Arabidopsis TERMINAL 116 FLOWER1 (29); while Dt2 is a MADS -domain factor gene (30). These are the 117 commonly included genes in soybean maturity models, as these genes have been 118 mapped in multiple studies, allele -specific markers have been developed, and genes 119 have been cloned (26–31). Eight (E1-E4, E9, Dt1, Dt2, and J) of these 14 genes have 120 the causal gene identified, with the corresponding gene sequence available from the 121 Williams 82 reference genome of soybean (32). 122 While the role of these genes on days to flowering and days to maturity has been 123 studied (22,33,34), minimal information is available on the effect of these genes during 124 different reproductive stages. In a quest for further information on the genetic control of 125 developmental stages, several genome -wide association studies have previ ously 126 reported on flowering and/or maturity timing (35–37). As mentioned previously, there is 127 a lack of information on intermediate stages between R1 and R8 (38). Additionally, as 128 plant breeders have traditionally worked with a narrow genetic base (39,40), a smaller 129 repertoire of diversity within development-related genes in the germplasm pool , 130 insufficient information is available about the genetics of maturity -influencing genes, 131 growth stage durations , and genetic control. In this study, we studied a large panel of 132 diverse soybean genotypes and a panel of near-isogenic lines (NILs) that were grown in 133 central Iowa for two years (2018 and 2019) and investigated the genetic variation for 134 eight reproductive growth stages, and the effect of environmental variables on different 135 reproductive stages (R1 -R8). The objectives of this work were to: (1) investigate the 136 genetic variation for reproductive growth stages in soybean , (2) perform genome-wide 137 association analysis on reproductive growth stages, and study the effect of E, J, and Dt 138 loci on reproductive growth stages , and (3) model the effect of genetic loci and 139 environmental variables influencing reproductive growth stages. 140 2 MATERIALS AND METHODS 141 2.1 Plant material and experimental design 142 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 6 A diversity pane l including SoyNAM parents and landrace accessions from the 143 core collection w as grown in 2018 (42°00'35.5"N, 93°44'01.9"W) and 2019 144 (42°01’9.32”N, 93°46’11.5”W) near Ames, IA , and was used in this study. In 2018, 5 08 145 lines were included within this study. The accessions that were viny-type, agronomically 146 inviable, or did not produce sufficient seed in the winter nursery increase were removed 147 from the experiment grown in 2019. Subsequently, 450 genotypes were included in the 148 2019 field study. The full list o f genotypes grown in each year is provided in 149 Supplementary File 2 - Table S1. In 2018 and 2019, a concurrent experiment that 150 included 67 near-isogenic lines ( NILs) for 14 maturity -related genes was planted in the 151 same field (Supplementary File 2 - Table S2). These NILs were primarily designed in 152 Clark and Harosoy backgrounds, with allele swapping for one or more of the known 153 maturity or stem termination genes (23). 154 In 2018, the diversity panel (508 accessions) was planted on May 29th in an 155 RCB design with 6 reps, each plot consisting of a hill plot with three seeds planted per 156 hill in a 30” x 30” grid. A single, representative plant was harvested for each genotype at 157 full maturity, and hand -shelled for winter increase in Chile to minimize heterogeneity 158 within accessions. Two pods from each of these plants were stored for genotyping. In 159 2019, 450 genotypes were planted on June 4th in an RCBD with 5 reps, each plot 160 consisting of 7’ long paired rows on 30” rows. For the 67 near-isogenic line (NIL) panel, 161 the plots were organized in an RCB design blocked by planting date, with three planting 162 dates per year, staggered approximately 10-14 days apart to allow for dissection of 163 photoperiod and thermal effects on time to each reproductive stage. All plots (in both 164 years) consisted of hill plots with three seeds per plot on a 30” x 30” grid. Planting dates 165 in 2018 were: May 29, June 7, and June 21, while planting dates in 2019 were: June 4, 166 June 14, and June 24. In the later analysis in this study, each year was consi dered as 167 an environment. 168 2.2 Phenotyping and environmental variables 169 Throughout the growing season, the development stage data were recorded 170 when each plot reached R1, R2, R3, R4, R5, R6, R7, and R8 stages. We calculated the 171 duration of each reproductive stage transition, such as R1 to R8 (R1R8), R2 to R8 172 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 7 (R2R8), and so on (1). Each plot was phenotyped at least once every three days , and if 173 needed, backdating was applied to plots that appeared to reach a growth stage one or 174 two days prior to the rating date. There were four environmental variables used in this 175 study: day length (DL), growing degree days (GDD), diurnal temperature range (DTR), 176 and photothermal time (PTT). For estimating the environmental indexes, the Global 177 Historical Climatology Network (GHCN) database at the National Oceanic and 178 Atmospheric Administration (NOAA)’s National Centers for Environmental Information 179 (NCEI) (https://www.ncdc.noaa.gov/) was used for retrieving daily maximum (T max) and 180 minimum (T min) air temperature (Fahrenheit). DL was calculated by using the function 181 daylength in the R package geosphere. GD Ds were obtained by using the formula 182 ((Tmax + Tmin)/2 – Tbase), with T max greater than 100°F adjusted to 100°F, Tmin lower than 183 50°F adjusted to 50°F, and Tbase considered as 50°F. The PTT was estimated as GDU x 184 DL. The daily DTR was calculated by subtracting T min from T max (i.e., Tmax – Tmin), 185 without any adjustment. 186 2.3 Phenotypic variability and selection of environmental indices 187 For assess ing the variation and effect of genotype, replication, and year ( 2018 188 and 2019 ) on the phenotype, we performed analysis of variance (ANOVA) for each 189 reproductive stage, and their transitions , as well as thei r respective environmental 190 variables such as DL, GDD, DTR, and PTT for individual year s and combined . In the 191 combined analysis, each year was considered as an environment. The c oefficient of 192 variation (CV) and Tukey’s honest significant difference (HSD) te st were performed to 193 evaluate the phenotypic differentiation among the genotypes. Further, we conducted 194 principal component analysis (PCA) using the prcomp() function in R to examine 195 multivariate patterns among variables and evaluate how environment and ge netics 196 contribute to phenotypic structure. For the identification of environmental variables that 197 significantly contribute to the reproductive stages and their transitions, we used the 198 means of DL, GDD, DTR, and PTT from the years (environment) 2018 and 20 19. The 199 environmental variables for each reproductive stage and their transitions were 200 correlated with the respective reproductive stages and their transitions. The 201 environmental variables that showed the most significant correlation with the 202 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 8 reproductive stages and their transitions, as well as physiologically relevant, were used 203 in further statistical analysis as environmental variables. All analyses were performed 204 using R version 4.3.3 (41). 205 2.3.1 Effect of genotype, environment, and their interaction 206 The best linear unbiased estimate (BLUE) of the genotype effect was estimated 207 for each reproductive stage and their transitions using a linear model in the R package 208 lme4 (42) and modelled as: 209 Yijk =  + Ei + B(i)j + Gk + EGik + ijk 210 Where,  indicates the mean of the population; E i, B(i) j, and G k is the effect of 211 environment (i.e., year), replication nested in environment, and genotype, respectively; 212 EGik indicates the interaction between the environment and genotype; and e ijk is the 213 error. 214 The model considered the fixed effects of genotype, replication, environment, 215 and genotype x environment (G x E). The model calculates each genotype's mean trait 216 performance across environments and the impacts of genotype, environment, and their 217 interaction. After adjusting for replication and environmental effects, the genotype effect 218 is the best linear unbiased estimate of each genotype's performance across 219 environments. To understand how genotypes respond to environmental variations , a 220 joint regression analysis was performed (43,44). The environmental means w ere first 221 estimated for each environment (2018 and 2019), and then the first environmental mean 222 regression mod el (Finlay -Wilkinson Regression) (45) was fitted to each of the 223 genotypes to assess the relation between the traits and the environmental mean, where 224 the traits were modelled as a function of the environmental mean. The environmental 225 mean regression model was described as: 226 Yik = k + k.EMi + ik 227 Where, Y ik is the value of k in the environment i; k is the intercept (baseline 228 performance of genotype k; k indicates the response of the genotype k to environment; 229 EMi indicates the mean trait value of environment i; and ik indicates residual error 230 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 9 In the second joint regression, a reaction norm was performed to understand the 231 G x E interaction with the identified environmental variables. To investigate the genotype 232 plasticity in response to the environmental variables, fitted values from the first 233 regression were used, where the identified environmental variables replaced the 234 environmental mean and were used as the explanatory variable. The model yielded the 235 slope (plasticity) and intercept of the reaction norm. 236 Yijk = k + k.EVijk + ijk 237 Where, Yijk is the fitted trait value for genotype k, environment i, replicate j from 238 the BLUE model; k indicates genotype -specific intercept; k is the reaction norm 239 plasticity slope in response to environmental variable (EV); EV ijk is a scaled 240 environmental variable; and ijk is the residual error 241 Variance component analysis (VCA) was performed using a mixed -effect model 242 applying the VCA() function in R to partition the variance of each reproductive stage and 243 their transitions into component s attributable to genotype, environment, replication, and 244 their interaction (G x E). This analysis allows for a deeper understanding of the factors 245 that contribute to the variability in the trait and guide s breeding strategies. Furthermore, 246 the magnitude of G x E interaction for each genotype was quantified by estimating the 247 mean absolute residual from the BLUE model. This measure provide d an estimate of 248 the genotype's performance across the environments. Based on the reaction norm and 249 G x E magnitude of the genotypes to the specific environmental variables, they are 250 classified into stable, plastic, and environment-dependent genotypes. 251 2.4 Genome-wide association studies (GWAS) 252 Genome-wide marker data for the diversity panel w ere downloaded from 253 SoyBase using Wms82.a2.v2 as the reference genome (46). GWAS analysis was 254 performed considering the attribu table component for each of the reprod uctive stages 255 and their transitions. In the GWAS analysis, we used BLUE estimates for capturing the 256 SNP related to genotypic effect only ; we also used reaction norm parameters, such as 257 plasticity slope, to identify the SNPs related to the effect of the environmental variable, 258 and intercept for identifying SNPs related to the average performance of the population 259 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 10 across all the environments. GWAS was conducted using four models: FarmCPU 260 (Random Model Circulating Probability Unification), BLINK (Bayesian -information and 261 Linkage Disequilibrium Iteratively Nested Keyway), MLM ( Mixed Linear Model), and 262 MLMM (Multi -Locus Mixed Linear Model) (47–49), and run in R using the GAPIT v3 263 (Genome A ssociation and Prediction Integrated Tool) R package (50). The GAPIT's 264 Bonferroni correction (0.05/number of SNPs) and a significance threshold of P -value 265 0.05 were used for multiple hypothesis testing. We used these four models to utilize 266 their statistical and computational power for capturing SNPs related to known maturity -267 related genes and to identify novel candidate genes. In our study, the known maturity-268 related genes captured by any of these model s were included; however, for the novel 269 candidates, o nly the SNPs that were captured by more than one model were 270 considered. The identified SNPs were annotated using the SoyBase database using 271 Wms82.a2.v2 as the reference genome (32). We further performed a multi -step 272 haplotype analysis utilizing the pegas package (51) in a custom R script to examine 273 haplotype structure and linkage disequilibrium (LD) around important SNPs found 274 through GWAS. It enabled the finding of local SNP blocks in strong LD with each GWAS 275 SNP and described their haplotype patterns across the genotypes. To evaluate pairwise 276 LD, we calculated the Pearson correlation coefficient with every other SNP in the 277 dataset. SNPs that had an absolute correlation of at least 0. 8 were considered in LD 278 and added to a particular GWAS SNP's haplotype block. 279 2.5 De novo marker analysis and development of a genetic maturity model 280 For de novo marker analysis , we conducted Kompetitive Allele Specific PCR 281 (KASP) genotyping (LGC Genomics) based on previously reported variation within E1, 282 E2, E3, E4, E9, E10, Dt1, Dt2, and J. Assays for twenty -two markers were developed 283 and deployed across the near-isogenic lines (NILs) and diversity panel. Of these, ten 284 markers showed no variation within either panel. The remaining markers (E1, E2, E3, 285 E4, Dt1, Dt2, J, E1E3, E1E4, E2E3, E2E4, and E3E4), were included in the analyses. 286 The E1E3, E1E4, E2E3, E2E4, and E3E4 indicated the interactions between the loci 287 that were not captured by individual loci. At first, we used the NIL panel for the 288 prediction of the timing of the reproductive stages (R1 -R8) through implementing a 289 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 11 Random Forest regression framework using caret and randomForest packages in R. 290 The dataset included genetic data such as genotype of known maturity and stem 291 termination genes and interactions between E loci. The dataset also included days to 292 reach each of the reproductive stages (R1 -R8) and photothermal variables (DL and 293 GDD) for each reproductive stage. In the model, the timing of the reproductive stages 294 (R1-R8) was designated as the target v ariables. Using 10 -fold cross -validation to 295 ensure model robustness, we trained a distinct Random Forest regression model for 296 each step (R1 –R8). The caret package's createDataPartition() was used to divide the 297 data into training (80%) and testing (20%) subsets. The performance of each model was 298 assessed using root mean square error (RMSE) and coefficient of determination (R 2) of 299 the test set after 500 trees (ntree = 500) were used for training. Further, the varImp() 300 function was utilized to extract the feature importance from the training models. 301 3 RESULTS 302 3.1 Phenotypic variation of the genotypes for reproductive maturity 303 Phenotypic variation was higher for environment 2019 compared to 2018 within the 304 diversity panel (Table 1, Supplementary File 1 - Figure S1) . Earlier maturity groups 305 reached each reproductive stage sooner than their late -maturing counterparts when 306 examined as a whole (Table 1, Supplementary File 1 - Figure S1A & S1C ). For 307 transitioning to the reproductive stages, such as flowering to pod initiation (R1 -R3), pod 308 development to seed filling (R4-R6), and towards full maturity (R7-R8), the average time 309 followed the same trend (Supplementary File 1 - Figure S1B & S1D, and Table S2). 310 However, within maturit y groups, there was still considerable variation . ANOVA results 311 for individual and combined environments showed significant genotypic effects (p<0.05). 312 In 2018 and 2019, the coefficient of variation (CV) was 9.15% and 8.90%, respectively 313 (Supplementary File 2 - Table S3, S4 & S 5). The CV in the combined dataset was 314 9.03%. We used PCA on standardized trait values from both environments to 315 demonstrate multivariate variance. Of the overall variation, PC1 and PC2 accounted for 316 21.99% and 18.14%, respectively. A considerable divergence between the environment 317 in 2018 and 2019 was shown by the PCA scatter plot (Figure 1B , and Supplementary 318 File 2 - Table S6). Despite some overlap, the genotypes from 2018 and 2019 formed 319 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 12 separate clusters that show ed environmental difference s. The top contributors to PC1 320 were the PTT and GDD, which were associated with later transitions, such as R2R7 and 321 R1R8. On the other hand, PC2 was mostly driven by DL during earlier transitions , R3 to 322 R5 (Figure 1C). 323 Table 1. Days from planting to each reproductive development stage in both years R1 R2 R3 R4 R5 R6 R7 R8 2018 All (450)a 30-81b (47.95c, 46d) 36-81 (49.69, 47) 42-97 (60.86, 59) 46-101 (67.91, 67) 52-108 (75.02, 75) 63-116 (90.50, 90) 80-135 (108.32, 109) 87-155 (117.37, 116) MG 0 (66) 30-55 (40.14, 39) 36-55 (42.03, 41) 42-65 (50.79, 49) 47-76 (57.70, 58) 52-81 (64.28, 64) 63-99 (78.12, 78) 80-112 (96.17, 96) 87-116 (103.84, 105) MG 1 (95) 31-67 (42.79, 42) 36-69 (44.14, 43) 43-83 (55.16, 55) 48-86 (62.38, 62) 57-88 (69.45, 68) 69-110 (85.03, 85) 81-120 (103.29, 105) 88-130 (110.78, 111) MG 2 (137) 36-76 (48.41, 46) 38-77 (50.27, 48) 42-84 (61.98, 62) 46-90 (68.93, 68) 58-98 (75.92, 75) 69-111 (91.90, 91) 91-125 (109.42, 110) 105-147 (116.87, 116) MG 3 (150) 36-81 (53.56, 53) 38-81 (55.40, 54) 46-97 (67.10, 65) 52-101 (74.15, 73) 58-108 (81.58, 81) 70-116 (97.11, 97) 96-135 (114.82, 114) 103-155 (126.70, 124) MG 4 (2) 43-54 (46.85, 46) 44-53 (47.77, 47) 49-77 (63.62, 64) 66-84 (73.46, 69) 69-89 (80.54, 82) 91-105 (99.92, 98) 108-127 (120.23, 121) 114-149 (136.62, 139) 2019 All (450) 34-75 (48.72, 47) 36-78 (52.06, 50) 42-84 (61.59, 61) 51-94 (69.77, 70) 56-98 (76.17, 76) 71-110 (88.80, 88) 85-136 (106.07, 107) 92-143 (115.62, 116) MG 0 (66) 34-71 (39.66, 38) 36-73 (43.32, 38) 42-83 (52.32, 51) 52-89 (59.96, 59) 56-92 (66.34, 65) 71-106 (80.47, 80) 85-118 (94.97, 94) 92-138 (101.29, 100) MG 1 (95) 34-70 (44.32, 43) 37-72 (48.05, 47) 45-81 (57.79, 57) 54-88 (65.62, 64) 60-92 (72.33, 71) 72-102 (85.76, 86) 91-123 (101.88, 101) 96-142 (109.40, 107) MG 2 (137) 36-71 (49.95, 47) 38-73 (53.40, 51) 48-81 (63.03, 62) 54-88 (71.40, 71) 61-96 (77.76, 77) 75-104 (90.12, 90) 89-136 (107.14, 107) 97-143 (116.66, 117) MG 3 (150) 37-75 (54.41, 54) 40-78 (57.27, 57) 47-84 (66.75, 66) 51-94 (75.18, 75) 56-98 (81.40, 82) 72-110 (93.12, 92) 93-136 (112.44, 112) 93-143 (124.68, 124) MG 4 (2) 40-48 (45.82, 47) 44-51 (49.36, 50) 55-72 (64.27, 62) 64-81 (74.45, 76) 76-85 (82.27, 82) 88-100 (94.00, 94) 107-127 (119.91, 122) 124-142 (135.64, 138) anumber genotypes used for the maturity group bminimum and maximum number of days to reach the reproductive stage .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 13 caverage number of days to reach the reproductive stage dmedian number of days to reach the reproductive stage MG – maturity group, R1-R8 – soybean reproductive stages 324 325 Figure 1. Phenotypic variations during reproductive stages and their transitions in a 326 soybean diversity panel. (A) Depicting different reproductive stages (R1 -R8) and their 327 transitions (R1R8-R7R8), (B) PCA scatter plot with 95% confidence ellipses indicating 328 genotype distributions for 2018 (dark blue) and 2019 (yellow), and (C) Top contributors 329 to principal component 1 (PC1) and principal component 1 (PC2) are shown in bar plot 330 form. The variables PTT, GDD, and DL indicate photothermal time, growth degree days, 331 and day length, respectively. The contributions of each trait to the primary axis are 332 reflected in the loadings. Part of this figure was created in https://BioRender.com. 333 3.2 Effect of environmental variables and genotype-environment (G x E) interaction 334 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 14 We evaluated correlations between reproductive phases and four environmental 335 variables: DL, DTR, GDD, and PTT to investigate how environmental cues influence 336 reproductive development (Fig ure 2A, Supplementary File 1 - Figure S2 ). The 337 reproductive stages showed a strong correlation with DL and GDD (Figure 2B, 338 Supplementary File 2 - Table S7). Environmental factors also had a significant impact 339 on transition times, albeit the effects varied depending on the kind of environment 340 (Supplementary File 1 - Figure S3). Based on these results, DL and GDD were 341 identified as important regulators as the y represent two separate biological 342 mechanisms, such as photoperiod and thermal time , that control flo wering and 343 reproductive development. 344 345 Figure 2. Correlation between reproductive stages and environmental variables in a 346 soybean diversity panel. (A) Pearson correlation coefficients for environmental variables 347 such as day length (DL), and growth degree days (GDD) with reproductive stages (R1–348 R8), (B) Top ten correlations between environmental variables and reproductive stages. 349 Genotype and G x E interaction effects significantly contributed to the observed 350 phenotypic variability for reproductive stages and their transition traits. The DL and GDD 351 showed a wide range of environment -adjusted means for each reproductive stage (R1 –352 R8). The VCA revealed that genotype effects accounted for the lowest phenotypic 353 variation at R6 (66%) to the highest at R1 (81%), with the G x E interaction accounting 354 for an additional 6.81% at R1 to 8.08% at R6, while replication and residual variance 355 proportions were quite minor (Supplementary File 2 - Table S8). The developmental 356 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 15 trajectory of most genotypes for both DL and GDD models w as parallel or slightly 357 divergent, suggesting modest plasticity and G x E interaction (Supplementary File 1 - 358 Figure S4A and S4B). The average responses at R1 and R2 for DL (0.040 and 0.055, 359 respectively) were marginally stronger than those for GDD (0.049 and 0.056, 360 respectively). However, the GDD model revealed significantly higher plasticity in the 361 mid-reproductive stages (R3 -R5): mean slopes at R3 and R5 were 0.029 and 0.026, 362 respectively. For DL, mean slopes were higher: mean slopes at R3 and R 5 were 0.034 363 and 0.028, respectively. At R8, DL had similar slope as GDD (DL mean slope = 0.026; 364 GDD = 0.027) (Supplementary File 2 - Table S9). Genotype effects for reproduction 365 transitions accounted for an average of about 33.3% of the total variance , with a n 366 additional 12.4% explained by the G x E interaction (Supplementary File 2 - Table S10). 367 Reaction norm analysis of transition traits exhibited significant variance in the flexibility 368 of genotypes throughout developmental intervals. According to mea n absolute slope 369 values, the first transition (R1R2) had the highest flexibility , and DL -based models 370 captured somewhat larger average slopes (0.3225) than GD D-based models (0.3054). 371 The mid-stage transitions (such as R2R6, R4R6, and R5R6) responded better to GDD, 372 indicating that thermal accumulation is a key factor in controlling these transitions. Slope 373 values remained high for later transitions (R6R7 –R7R8), but variation increased, 374 especially for DL, suggesting greater sensitivity to coupled environmenta l signals or 375 genotypic variability (Supplementary File 2 - Table S11). 376 3.3 Effect of genomic loci during reproductive stages and their transitions 377 We conducted the GWAS using BLUE estimates from both DL and GDD models 378 for the reproductive stages and their transitions. Two SNPs encoding the Dt1 gene 379 (Glyma.19G194300) and E2 gene (Glyma.10G221500) were consistently detected at 380 chromosome 19 and chromosome 10, respectively, for both DL and GDD 381 (Supplementary File 1 - Figure S5A & S5B, and Supplementary File 2 - Table S12 & 382 S13). The Dt1 gene was detected for the mid reproductive stages, such as R4-R6, while 383 the E2 gene was found to be associated with later maturity stages, R7 and R8. 384 GWAS using two reaction norm parameter s, plasticity slope and intercept , 385 confirmed the importance of Dt1. Using slope, the Dt1 gene was detected for the R5 386 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 16 stage using the DL data across genotypes. E2 gene was not detected for any stage 387 using the DL data . However, we identified SNPs for the E1 gene (Glyma.06G207800) 388 on chromosome 6 for reproductive stage R2, R6, and R8 (Supplementary File 1 - Figure 389 S5C, Supplementary File 2 - Table S12 & S13). The E1 gene was also detected for the 390 R2 stage as well as when intercept was used as a response variable using the GDD 391 data. None of the known maturity -related genes was detected for both plasticity slope 392 and intercept. However, we found one QTL from DL data at position Chr8: 17603027 393 (Glyma.08G216800) associated with plasticity slope and intercept for R3 stage, as well 394 as BLUE and plas ticity slope for R1 stage. This QTL was identified as PPPDE putative 395 thiol peptidase family protein (Supplementary File 2 - Table S12). For GDD, we also 396 identified one SNP at position Chr4:7556514 associated with plasticity slope and 397 intercept for the R7 stage. Using Soybase genome browse, we found that th e nearest 398 (within 1Kbp) gene , annotated as Glyma.04G088100, which is an RNA-binding KH 399 domain-containing protein. In case of reproductive stage transitions using DL, we 400 identified the Dt1 and Dt2 genes on chromosomes 19 and 18, respectively, and the E1 401 gene on chromosome 6 (Supplementary File 2 - Table S14). For GDD, we also detected 402 the Dt1 and Dt2 genes on chromosomes 19 and 18, respectively, and the E1 and E2 403 genes were detected on chromosomes 6 and 10, respectively (Supplementary File 2 - 404 Table S15). 405 Besides the above -mentioned maturity -related genes, several high -confidence 406 SNPs were detected from multiple GWAS methods. Among them, most of the SNPs 407 were unknown or uncharacterized. Considering the DL environmental gradient, a few 408 top SNPs that were detected through three or four GWAS models were: 409 Chr01:51691145 ( Glyma.01G180600), Tetratricopeptide repeat (TPR) -like superf amily 410 protein (BLUE: R1, R2, R3, R4, R5, R6; Intercept: R2, R3, R4); Chr06:3396872 411 (Glyma.06G044600), NAD(P) -binding Rossmann -fold superfamily protein (Slope: R1, 412 R3, R5, R7); Chr08:1017668 ( Glyma.08G013000), Myosin heavy chain -related protein 413 (Slope: R4); and Chr10:45269968 (Glyma.10G221300), S-adenosylmethionine carrier 1 414 (BLUE: R7, R8; Intercept: R1, R5, R6, R7, R8). The SNP at position Chr01:51691145 415 was also detected for GDD through three or four GWAS models, and the SNP at 416 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 17 Chr10:45269968 was found to be strongly connected with the E2 gene at the haplotype 417 network. 418 Haplotype blocks were built around each GWAS -significant SNP based on SNPs 419 in high linkage disequilibrium (LD ≥ 0.8) to understand the genetic structure and linkage 420 patterns of detected matur ity genes. The analysis of six GWAS SNPs revealed 68 421 distinct haplotypes, indicating a significant level of allelic variation within genomic 422 segments characterized by LD (Supplementary File 2 - Table S16). LD correlation (r) 423 values ranged from 0.91 to 1.0 , and a mean r value of 0.96 . Ten LD blocks were found 424 for the top five GWAS hits, highlighting their potential as key hubs in the design of 425 haplotype networks (Supplementary File 1 - Figure S 6, and Supplementary File 2 - 426 Table S1 7). Compared to single SNPs , these haplotype blocks offer a finer -scale 427 resolution for association mapping and might be better at capturing causal variation. 428 Strong LD and significant haplotype diversity point to strong genomic areas that need 429 more functional investigation, especially considering the observed trait variability. 430 3.4 Random Forest-based maturity model dissecting genetic and environmental 431 contributions to soybean phenology 432 Random forest (RF) models trained on the NIL and diversity panels exhibited 433 high predictive power across soybean reproductive stages (R1 to R8), with some 434 variability observed between early and late stages. With RMSE values ranging from 435 0.01 to 1.98 and corresponding R 2 values between 0.96 and 0.99, the model with NIL 436 panel performed well for all reproductive phases except R8 (Supplementary File 1 - 437 Table S3). On the other hand, the same model using the diversity panel demonstrated 438 consistently better accuracy at every level, including R8. Different patterns of 439 environmental variable relevance over phases were found by the feature importance 440 analysis based on the percentage increase in mean squared error (MSE) (Figure 3A 441 and Figure 4A). Environmental predictors were the main factors affectin g stage 442 transitions in both panels, especially GDD and DL. In particular, GDD at R1 and R2 443 stages continuously received the greatest relevance ratings during early stages (R1 –444 R3) in both panels. As we moved toward a more complex environmental integration, 445 both GDD and DL from prior stages (e.g., GDD_R4, DL_R5) helped to explain variance 446 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 18 for intermediate stages (R4 –R6). The diversity panel indicated that DL_R8 and 447 GDD_R8 were predominantly dominant for late stages (R7 –R8), whereas the NIL panel 448 emphasized DL_R2 and DL_R7 (Figure 3A and Figure 4A). Throughout most phases, 449 maturity-related genes ( E1–E4, Dt1, Dt2) made a minor contribution to the R8 stage 450 prediction, particularly in the diversity panel. 451 The RF model also provided information regarding specific developmental 452 stages. For instance, at growth stage R8, the partial dependency plots showed that E1, 453 E2, and Dt1 genes had distinct stepwise effects in both panels (Figure 4B and Figure 454 4B). However, overall marginal effects on the anticipated R8 value wer e small, which is 455 consistent with their low feature selection relevance rankings. Figures 3C and 4C show 456 the overall interaction contributions of each variable used in the model as measured by 457 the interaction strength analysis. The NIL panel's most interactive features were DL at 458 the R2 and R7 stages, followed by DL at the R3 and R5 stages, indicative of NIL’s 459 sensitivity to cross-stage photoperiods. Particularly in the diversity panel, where DL and 460 GDD at the R8 stage were predominant drivers of interaction effects, genetic factors (E-461 loci) exhibited the least amount of interaction strength, but environmental predictors (DL 462 and GDD) contributed the most to interaction terms. 463 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 19 464 Figure 3. Random Forest -based maturity model of the reprodu ctive stages from NIL 465 population. (A) Feature importance by stage as a percentage increase in mean square 466 error (MSE) for reproductive stages R1 through R8, (B) One -dimensional partial 467 dependence plot (PDP) that illustrates how E1, E2, and Dt1 marginally affect the 468 expected R8 values, and (C) Analysis of feature interaction strength for stage R8. DL 469 and GDD indicate day length and growing degree days, respectively. 470 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 20 471 Figure 4. Random Forest -based maturity model of the reproductive stages from th e 472 diversity panel. (A) Feature importance of the variables used in the model for 473 reproductive stages as a percentage increase in mean square error (MSE), (B) One -474 dimensional partial dependence plot (PDP) showing the effect of E1, E2, and Dt1 gene 475 on R8 sta ge, and (C) Feature interaction strength for R8 stage. DL and GDD indicate 476 day length and growing degree days, respectively. 477 4. DISCUSSION 478 In this study, we presented results from a developmental perspective of 479 dissecting phenotypic variation in reaching e ach of the reproductive stages under 480 natural field conditions , considering environmental gradients, and providing insight into 481 the genetic loci regulated by the effect of genotype, environ mental variable, and their 482 interaction. We provide the first compreh ensive stage interval study from R1 to R8, 483 elucidating novel insights. Finally, we developed a random forest -based model that can 484 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 21 predict the effects of genotype and environmental variables in the reproductive stages 485 of soybean. 486 Using two years of data and a diversified soybean panel consisting of >450 PI 487 accessions and >65 NILs , our study showed significant phenotypic variation in 488 reproductive stages and their transitions. The change in the environment highlights how 489 day length has a significant impact on the early stages of reproductive development in 490 soybean. These kind s of variations due to genetic and non -genetic factors were also 491 reported previously (52,53). These findings highlight that, even for the reproductive 492 stages under strong genetic cont rol, their expression is context-dependent. 493 Environmental factors become the main determinants of reproductive development . The 494 temperature and photoperiod control flowering and maturity through mostly separate but 495 temporally synchronized processes (14,54,55). In this study, GDD accounted for more 496 variance in mid-to-late reproductive stages, which is consistent with the requirement for 497 thermal accumulation for pod filling and maturity, whereas DL had a significant influence 498 on early transitions (Figure 2A). Early transitions are driven by photoperiod, where short 499 DL triggers early flowering while long DL delayed flowering (56,57). G x E occurs 500 frequently, is linked to a variety of genetic factors and molecular events, and is 501 frequently brought about by changes in the strength of genetic effects in response to 502 environmental stimuli (58). In our study, we found stage -specific contributions of 503 genotype, environmental variables, and G x E interaction effects . For instance, early 504 stages like R1 dis played the highest genetic effect , which suggests that early floral 505 initiation is more genetically controlled than the start of pod filling. Early stages were 506 more receptive to DL, while mid stages are governed by GDD or heat accumulation that 507 is related t o pod setting and seed development , and later maturity is tandemly 508 influenced by both cues . These kinds of dynamic environmental influences on 509 reproductive trajectories were also reported by earlier studies (59,60), which have 510 shown that this sensitivity to temperature and photoperiod varies between different 511 maturity groups and genetic origins. Photoperiod sensitivity is often higher in later 512 maturity groups, and high temperatures can either increase or decrease these 513 sensitivities (61). The cumulative impacts of earlier and later developmental stimuli were 514 also incorporated in transition stages (e.g., R1R8, R2R8). When compared to individual 515 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 22 reproductive stages, the transition stages showed larger G x E and residual variance, 516 indicating that the impact of the environment increased when considering reproductive 517 transition stages . With the large genetic panel and reproductive stages and their 518 transitions, these findings further strengthen our understanding of the gen otypic effect 519 and its int eraction with the environment, both for reproductive stages and their 520 transitions. However, genetics is one of the main factors influencing growth stages and 521 their intervals. 522 GWAS from this study identified key loci, including Dt1, E1, and E2, for different 523 reproductive stages and their transitions (Supplementary File 1 - Figure S5), which were 524 also detected in a meta-GWAS study and known to regulate flowering and maturity (62). 525 Recent research showed that E1Lb, independent of E1, delays flowering in long -day 526 circumstances (63). The E2 gene, orthologous to the GIGANTEA (GI) gene in 527 Arabidopsis thaliana , was found to be involved in delaying flowering by inhibiting the 528 expression of GmFT2a during long days (26). The Dt1 gene is involved in the 529 determination of plant growth type , such as determinate or indeterminate type (64). 530 These genes' pleiotropic roles across reproductive phases were confirmed by the 531 detection of these genes not only for reproductive stages but also for sta ge transitions. 532 The findings of the study conducted by Miranda et al. (2020) enabled the comparison of 533 various variant alleles of those genes and showed notable genotype -based differences 534 between days to flowering and days to maturity (34). We report nov el loci, including 535 Glyma.08G216800 (PPPDE putative thiol peptidase family protein) and 536 Glyma.04G088100 (RNA-binding KH domain-containing protein) from GWAS with slope 537 and intercept of the reaction norm , which revealed more levels of regulatory intricacy 538 that may be connected to environmental responsiveness , which varies by reproductive 539 stage. The PPPDE putative thiol peptidase family protein was identified for early 540 reproductive stages, indicating that, irrespective of genotype and environmental effect, it 541 might have a regulatory role with photoperiod during the early reproductive stage . But 542 its functions are yet unclear and need additional research. It has been reported that t he 543 PPPDE superfamily is a deubiquitinase (DUB) that is conserved in eukaryotes, including 544 humans (65). Arabidopsis AtC3H59 controls cell division b y interacting with the PPPDE 545 family protein Desi1 via its WD40 domain, and the PPPDE domain of the deubiquitinase 546 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 23 PICI1 serves as an immune hub for pattern-triggered immunity and effector -triggered 547 immunity in rice (66). We also identified t he RNA -binding K homology ( KH) domain-548 containing protein associated with the late reproductive stage , suggesting that this 549 protein might be regulated by thermal o r heat accumulation of plants rather than by the 550 genotype effect. Guan et al. (2013) reported that RCF3 is a nuclear -localized putative 551 RNA-binding protein that contains a KH domain. The rcf3 Arabidopsis mutant was more 552 resistant to heat stress than the wi ld-type, which is consistent with the general higher 553 accumulation of heat -responsive genes (67). However, in our study, we did not 554 investigate the effect of temperature on the reproductive stages, and we did not 555 functionally validate the RNA -binding K homology (KH) domain -containing protein. The 556 novel SNP at position Chr01: 51691145 ( Glyma.01G180600), which is a 557 tetratricopeptide repeat (TPR) -like superfamily protein detected from BLUE and 558 intercept for both DL and GDD for early to mid -reproductive stages, indicated that this 559 protein could be a potential target for genetic regulation of reproductive st ages. This 560 TPR protein has been reported to be involved in hormonal regulation and drought stress 561 in plants (68,69). In the haplotype network , the E2 gene or the GIGANTEA protein was 562 closely connected with NOD26 -like intrinsic protein, Nucleotide -sugar transporter family 563 protein, Ribosomal protein S13/S18 family, S -adenosylmethionine carrier 1, Inositol 564 polyphosphate 5-phosphate-related protein, and five other uncharacterized or unknown 565 proteins. Among the uncharacterized or unknown proteins, Glyma.10G221900 had the 566 highest LD with the E2 gene (Supplementary File 1 - Figure S 6). According to the 567 reported function of the above genes that were in high LD with the E2 gene, they are 568 mainly involved in the growth and development of plants in several cellular pathways 569 (70–78). However, for better understanding of the mechanism of development during 570 maturity could be possible by discovering the function of unknow n or uncharacterized 571 proteins found in strong LD and varied haplotype structures surrounding GWAS -572 significant SNPs. 573 The Random Forest models corroborated experimental findings by attributing 574 dominant predictive power to environmental variables, especially DL and GDD, with 575 minor contributions from known maturity genes (Figure 3 and 4). The difference in 576 prediction accuracy between the diversity and NIL panels demonstrates how more 577 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 24 genetic variation improves the generalizability of the model. The improved mod el 578 performance of the diversity panel could be attributed to increased resolution and 579 variability in genotype and environmental responses, which would enable the model to 580 more effectively partition and learn intricate patterns influencing phenology. Analysis of 581 feature importance and interaction strength showed a temporal pattern of cue 582 dominance. It's interesting to note that even well -known maturity genes had negligible 583 effect sizes in prediction models, confirming the idea that the environment mostly 584 controls the timing of reproductive development in a diverse population. Since this study 585 didn’t include commercial varieties, it is difficult to generalize these findings to advanced 586 breeding stock. However, the large genetic collection and NILs were very us eful to 587 develop a prediction model that can be utilized t o predict stage-specific regulation of 588 genotypic and environmental factors with applications in breeding programs . These 589 discoveries open the door to more focused breeding tactics that take advantage of 590 developmental plasticity and environmental sensitivity. For example, breeders can 591 leverage this information for ideotype breeding to maximize crop productivity by 592 optimizing reproductive stages and their transitions in conjunction with physiological 593 characteristics in cultivation regions (6). The crop modelers can utilize the data and 594

presented in this study for more rigorous crop modeling applications. However, 595 a common approach for such studies is to perform multi -year studies at the same 596 geographical location to minimize the effect of photoperiod, and to some extent, 597 temperature or thermal accumulation. While less comprehensive and a more limited 598 scope of results interpretation and extrapolation than an extensive study across multiple 599 years, locations, and geographies, it has advantages due to better and more 600 standardized experimental factors. Also, central Iowa in the U.S. falls in primarily the 601 MGII maturity group, which has among the highest soybean production in the U.S. and 602 Canada (13,79). 603 5. CONCLUSION 604 From a comprehensive viewpoint based on the results of phenotypic and 605 genomic analysis and s tatistical modeling using a large panel of genotypes , the study 606 demonstrated that intricate interactions between genotype, environment, and their 607 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 25 interaction drive the reproductive development of soybeans. This study provides 608 practical insights for breeders looking to optimize soybean maturity and yield under 609 varying environmental grad ients for applying ideotype breeding concepts. Additionally, 610 farmers can make better decisions about variety selection by identifying stage -specific 611 environmental sensitivity, where photoperiod (DL) mostly influences early and late 612 reproductive stages and thermal time (GDD) primarily influences mid -stages. Growers 613 can lessen the chance of yield loss during crucial reproductive times due to heat, 614 drought, or frost stress by coordinating reproductive transitions with favorable 615 photothermal windows. The identification of novel loci governing plasticity and the major 616 maturity genes ( E1, E2, Dt1, and novel candidate genes identified in this study ) gives 617 plant breeders a genetic toolbox for creating varieties with specific adaption. Breeders 618 can find genotypes with advantageous environmental responsiveness under a variety of 619 situations by combining BLUE -based genetic value estimate s, reaction norm modeling, 620 GWAS, and a random forest -based predictive model . When taken as a whole, these 621

lend credence to breed ing plans that try to increase phenological plasticity and 622 give farmers the ability to match genotypes to certain environmental profiles, which 623 eventually improves output stability and resource efficiency. However, our study was 624 limited to central Iowa onl y and was conducted for two years. Multi-year data from 625 multiple geographical locations could provide more rigorous results of the effect of 626 environmental gradients and G x E interaction. To speed up breeding for environment -627 specific maturity optimization, future studies should concentrate on the functional 628 validation of novel loci and the mechanistic analysis of genotype -by-environment 629 responsiveness. 630 Data availability: 631 Upon the publication of this article, all data related to this study will be available at 632 https://github.com/SoylabSingh/SoyMaturity 633 Author Contributions: 634 Conceptualization: SC, JMS, AKS; Data curation: JMS; Methodology: SC, JMS, AKS; 635 Formal data analysis and visualization: SC, JMS; Writing original draft: SC; Review and 636 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 26 editing: SC, AKS; Project administration , funding acquisition , resources, and 637 supervision: AKS. 638

639 J. Shook was partially supported by funding from the NSF -NRT. Funding support came 640 from R. F. Baker Center for Plant Breeding, Iowa Soybean Association, USDA -NIFA, 641 and National Science Foundation. Research support from members of AKS group at 642 ISU - staff, post-doctoral fellows, and students - is sincerely acknowledged. 643 Conflict of interest statement 644 The authors declare no conflicts of interest. 645

646 1. Fehr W, Caviness C. Stages of soybean development. 1977; 647 2. FACTS. Soybean Planting Decision Tool, Iowa State University. 2014 [cited 2025 648 Aug 8]. Understanding interactions between soybean planting date and maturity 649 across environments. Available from: https://facts.extension.iastate.edu/soybean -650 planting-decision-tool 651 3. Klubertanz TH. Effects of drought stress and defoliation on soybean physiology 652 and yield. Iowa State University; 1993. 653 4. Foroud N, Mündel HH, Saindon G, Entz T. Effect of level and timing of moisture 654 stress on soybean plant development and yield components. Irrig Sci. 655 1993;13(4):149–55. 656 5. Singh DP, Singh AK, Singh A. Plant breeding: past, present, and future 657 perspectives. Plant breeding and cultivar development. 2021;1–24. 658 6. Singh DP, Singh AK, Singh A. Breeding of crop ideotypes. Plant Breeding and 659 Cultivar Development; Academic Press: Cambridge, MA, USA. 2021;497–516. 660 7. Nico M, Miralles DJ, Kantolic AG. Natural post -flowering photoperiod and 661 photoperiod sensitivity: Roles in yield -determining processes in soybean. Field 662 Crops Res. 2019;231:141–52. 663 8. Chen G, Wiatrak P. Soybean develop ment and yield are influenced by planting 664 date and environmental conditions in the southeastern coastal plain, United 665 States. Agron J. 2010;102(6):1731–7. 666 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 27 9. Nleya T, Sexton P, Gustafson K, Miller JM. Soybean growth stages. IGrow 667 soybean: Best Management Practices for soybean production. 2013;3. 668 10. Mu Q, Guo T, Li X, Yu J. Phenotypic plasticity in plant height shaped by interaction 669 between genetic loci and diurnal temperature range. New Phytologist. 670 2022;233(4):1768–79. 671 11. Napier JD, Heckman RW, Juenge r TE. Gene -by-environment interactions in 672 plants: Molecular mechanisms, environmental drivers, and adaptive plasticity. 673 Plant Cell. 2023;35(1):109–24. 674 12. Bonamour S, Chevin LM, Charmantier A, Teplitsky C. Phenotypic plasticity in 675 response to climate chan ge: the importance of cue variation. Philosophical 676 Transactions of the Royal Society B. 2019;374(1768):20180178. 677 13. De Bruin JL, Pedersen P. Soybean seed yield response to planting date and 678 seeding rate in the Upper Midwest. Agron J. 2008;100(3):696–703. 679 14. Mourtzinis S, Gaspar AP, Naeve SL, Conley SP. Planting date, maturity, and 680 temperature effects on soybean seed y ield and composition. Agron J. 681 2017;109(5):2040–9. 682 15. Garner WW, Allard HA. Photoperiodic response of soybeans in relation to 683 temperature and other environmental factors. United States Department of 684 Agriculture Washington, DC, USA; 1930. 685 16. Egli DB, Wa rdlaw IF. Temperature response of seed growth characteristics of 686 soybeans 1. Agron J. 1980;72(3):560–4. 687 17. Board JE, Hall W. Premature flowering in soybean yield reductions at nonoptimal 688 planting dates as influenced by temperature and photoperiod 1. Agro n J. 689 1984;76(4):700–4. 690 18. Seddigh M, Jolliff GD, Orf JH. Night temperature effects on soybean phenology. 691 Crop Sci. 1989;29(2):400–6. 692 19. Bernard RL. Two major genes for time of flowering and maturity in soybeans 1. 693 Crop Sci. 1971;11(2):242–4. 694 20. Yamanaka N, Ninomiya S, Hoshi M, Tsubokura Y , Yano M, Nagamura Y , et al. An 695 informative linkage map of soybean reveals QTLs for flowering time, leaflet 696 morphology and regions of segregation distortion. DNA research. 2001;8(2):61 –697 72. 698 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 28 21. Ellis RH, Asumadu H, Qi A, Summerfield RJ. Effects of photoperiod and maturity 699 genes on plant growth, partitioning, radiation use efficiency, and yield in soyabean 700 [Glycine max (L.) Merrill]‘Clark.’ Ann Bot. 2000;85(3):335–43. 701 22. McBlain BA, Bernard RL. A new gene affecting the time of flowering and maturity 702 in soybeans. Journal of Heredity. 1987;78(3):160–2. 703 23. Cober ER, Tanner JW, Voldeng HD. Soybean photoperiod‐sensitivity loci respond 704 differentially to light quality. Crop Sci. 1996;36(3):606–10. 705 24. Cober ER, Tanner JW, Voldeng HD. Genetic control of photoperiod response in 706 early‐maturing, near‐isogenic soybean lines. Crop Sci. 1996;36(3):601–5. 707 25. Xia Z, Watanabe S, Yamada T, Tsubokura Y , Nakashima H, Zhai H, et al. 708 Positional cl oning and characterization reveal the molecular basis for soybean 709 maturity locus E1 that regulates photoperiodic flowering. Proceedings of the 710 National Academy of Sciences. 2012;109(32):E2155–64. 711 26. Watanabe S, Xia Z, Hideshima R, Tsubokura Y , Sato S, Ya manaka N, et al. A 712 map-based cloning strategy employing a residual heterozygous line reveals that 713 the GIGANTEA gene is involved in soybean maturity and flowering. Genetics. 714 2011;188(2):395–407. 715 27. Watanabe S, Hideshima R, Xia Z, Tsubokura Y , Sato S, Nakamoto Y , et al. Map-716 based cloning of the gene associated with the soybean maturity locus E3. 717 Genetics. 2009;182(4):1251–62. 718 28. Liu B, Kanazawa A, Matsumura H, Takahashi R, Harada K, Abe J. Genetic 719 redundancy in soybean photoresponses associated with dupli cation of the 720 phytochrome A gene. Genetics. 2008;180(2):995–1007. 721 29. Liu B, Watanabe S, Uchiyama T, Kong F, Kanazawa A, Xia Z, et al. The soybean 722 stem growth habit gene Dt1 is an ortholog of Arabidopsis TERMINAL FLOWER1. 723 Plant Physiol. 2010;153(1):198–210. 724 30. Ping J, Liu Y , Sun L, Zhao M, Li Y , She M, et al. Dt2 is a gain -of-function MADS-725 domain factor gene that specifies semideterminacy in soybean. Plant Cell. 726 2014;26(7):2831–42. 727 31. Xia Z, Watanabe S, Yamada T, Tsubokura Y , Nakashima H, Zhai H, et al. 728 Positional cloning and characterization reveal the molecular basis for soybean 729 maturity locus E1 that regulates photoperiodic flowering. Proceedings of the 730 National Academy of Sciences. 2012;109(32):E2155–64. 731 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 29 32. Grant D, N elson RT, Cannon SB, Shoemaker RC. SoyBase, the USDA -ARS 732 soybean genetics and genomics database. Nucleic Acids Res. 733 2010;38(suppl_1):D843–6. 734 33. Liu L, Song W, Wang L, Sun X, Qi Y , Wu T, et al. Allele combinations of maturity 735 genes E1 -E4 affect adaptation of soybean to diverse geographic regions and 736 farming systems in China. PLoS One. 2020;15(7):e0235397. 737 34. Miranda C, Scaboo A, Cober E, Denwar N, Bilyeu K. The effects and interaction of 738 soybean maturity gene alleles controlling flowering time, maturity, and adaptation 739 in tropical environments. BMC Plant Biol. 2020;20(1):65. 740 35. Zhang J, Song Q, Cregan PB, Nelson RL, Wang X, Wu J, et al. Genome -wide 741 association study for flowering time, maturity dates and plant height in early 742 maturing soybean (Glycine max) germplasm. BMC Genomics. 2015;16(1):217. 743 36. Contreras-Soto RI, Mora F, Lazzari F, de Oliveira MAR, Scapim CA, Schuster I. 744 Genome-wide association mapping for flowering and maturity in tropical soybean: 745 implications for breeding strategies. Breed Sci. 2017;67(5):435–49. 746 37. Kim KH, Kim JY , Lim WJ, Jeong S, Lee HY , Cho Y , et al. Genome -wide 747 association and epistatic interactions of flowering time in soybean cultivar. PLoS 748 One. 2020;15(1):e0228114. 749 38. Li M, Liu Y , Tao Y , Xu C, Li X, Zhang X, et al. Id entification of genetic loci and 750 candidate genes related to soybean flowering through genome wide association 751 study. BMC Genomics. 2019;20(1):987. 752 39. Gizlice Z, Carter Jr TE, Burton JW. Genetic base for North American public 753 soybean cultivars released be tween 1947 and 1988. Crop Sci. 1994;34(5):1143 –754 51. 755 40. Wysmierski PT, Vello NA. The genetic base of Brazilian soybean cultivars: 756 evolution over time and breeding implications. Genet Mol Biol. 2013;36:547–55. 757 41. R Core Team R. R: A language and environment for statistical computing. 758 Citeseer; 2020. 759 42. Bates DM. lme4: Mixed-effects modeling with R. Springer New York; 2010. 760 43. Li X, Guo T, Mu Q, Li X, Yu J. Genomic and environmental determinants and their 761 interplay underlying phenotypic plasticity. Proceedings of the National Academy of 762 Sciences. 2018;115(26):6679–84. 763 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 30 44. Guo T, Mu Q, Wang J, Vanous AE, Onogi A, Iwata H, et al. Dynamic effects of 764 interacting genes underlying rice flowering -time phenotypic plasti city and global 765 adaptation. Genome Res. 2020;30(5):673–83. 766 45. Piepho H, Blancon J. Extending Finlay –Wilkinson regression with environmental 767 covariates. Plant Breeding. 2023;142(5):621–31. 768 46. Song Q, Hyten DL, Jia G, Quigley C V, Fickus EW, Nelson RL, et al. Development 769 and evaluation of SoySNP50K, a high -density genotyping array for soybean. 770 PLoS One. 2013;8(1):e54985. 771 47. Segura V, Vilhjálmsson BJ, Platt A, Korte A, Seren Ü, Long Q, et al. An efficient 772 multi-locus mixed -model approach for genome -wide association studies in 773 structured populations. Nat Genet. 2012;44(7):825–30. 774 48. Liu X, Huang M, Fan B, Buckler ES, Zhang Z. Iterative usage of fixed and random 775 effect models for powerful and efficient genome -wide association studies. PLoS 776 Genet. 2016;12(2):e1005767. 777 49. Huang M, Liu X, Zhou Y , Summers RM, Zhang Z. BLINK: a package for the next 778 level of genome-wide association studies with both individuals and markers in the 779 millions. Gigascience. 2019;8(2):giy154. 780 50. Wang J, Zhang Z. GAPIT version 3: boosting power and accuracy for genomic 781 association and prediction. Genomics Proteomics Bioinformatics. 782 2021;19(4):629–40. 783 51. Paradis E. pegas: an R package for population genetics with an integrated –784 modular approach. Bioinformatics. 2010;26(3):419–20. 785 52. Krause MD, Dias KOG, Singh AK, Beavis WD. Using soybean historical field trial 786 data to study genotype by environment variation and identify mega ‐environments 787 with the integration of genetic and non ‐genetic factors. Agr on J. 788 2025;117(1):e70023. 789 53. Bianchi MC, Bruzi AT, Soares IO, Ribeiro FDO, Gesteira GDS. Heritability and the 790 genotype× environment interaction in soybean. Agrosystems, Geosciences & 791 Environment. 2020;3(1):e20020. 792 54. Mathur S, Seo B, Reddy KR, Reddy VR . Morpho-phenological responses of early 793 and late soybean maturity groups to diverse day lengths and temperature 794 regimes. Ind Crops Prod. 2024;221:119362. 795 55. Jung C, Müller AE. Flowering time control and applications in plant breeding. 796 Trends Plant Sci. 2009;14(10):563–73. 797 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 31 56. Cober ER, Curtis DF, Stewart DW, Morrison MJ. Quantifying the effects of 798 photoperiod, temperature and daily irradiance on flowering time of soybean 799 isolines. Plants. 2014;3(4):476–97. 800 57. Mathur S, Seo B, Reddy KR, Reddy VR. Morph o-phenological responses of early 801 and late soybean maturity groups to diverse day lengths and temperature 802 regimes. Ind Crops Prod. 2024;221:119362. 803 58. Des Marais DL, Hernandez KM, Juenger TE. Genotype -by-environment 804 interaction and plasticity: exploring genomic responses of plants to the abiotic 805 environment. Annu Rev Ecol Evol Syst. 2013;44(1):5–29. 806 59. No DH, Baek D, Lee SH, Cheong MS, Chun HJ, Park MS, et al. High-temperature 807 conditions promote soybean flowering through the transcriptional reprograming of 808 flowering genes in the photoperiod pathway. Int J Mol Sci. 2021;22(3):1314. 809 60. Setiyono TD, Weiss A, Specht J, Bastidas AM, Cassman KG, Dobermann A . 810 Understanding and modeling the effect of temperature and daylength on soybean 811 phenology under high-yield conditions. Field Crops Res. 2007;100(2–3):257–71. 812 61. Ort NWW, Morrison MJ, Cober ER, Samanfar B, Lawley YE. Photoperiod affects 813 node appearance ra te and flowering in early maturing soybean. Plants. 814 2022;11(7):871. 815 62. Shook JM, Zhang J, Jones SE, Singh A, Diers BW, Singh AK. Meta -GWAS for 816 quantitative trait loci identification in soybean. G3. 2021;11(7):jkab117. 817 63. Zhu J, Takeshima R, Harigai K, Xu M, Kong F, Liu B, et al. Loss of function of the 818 E1-Like-b gene associates with early flowering under long -day conditions in 819 soybean. Front Plant Sci. 2019;9:1867. 820 64. Zimmer G, Miller MJ, Steketee CJ, Jackson SA, de Tunes LVM, Li Z. Genetic 821 control and allele variation among soybean maturity groups 000 through IX. Plant 822 Genome. 2021;14(3):e20146. 823 65. Xie X, Wang X, Jiang D, Wang J, Fei R, Cong X, et al. PPPDE1 is a novel 824 deubiquitinase belonging to a cysteine isopeptidase family. Biochem Biophys Res 825 Commun. 2017;488(2):291–6. 826 66. Zhai K, Liang DI, Li H, Jiao F, Yan B, Liu J, et al. NLRs guard metabolism to 827 coordinate pattern-and effector-triggered immunity. Nature. 2022;601(7892):245 –828 51. 829 67. Guan Q, Wen C, Zeng H, Zhu J. A KH domain -containing putative RNA -binding 830 protein is critical for heat stress -responsive gene regulation and thermotolerance 831 in Arabidopsis. Mol Plant. 2013;6(2):386–95. 832 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 32 68. Schapire AL, Valpuesta V, Botella MA. TPR protein s in plant hormone signaling. 833 Plant Signal Behav. 2006;1(5):229–30. 834 69. Liu X, Wang DR, Chen GL, Wang X, Hao SY , Qu MS, et al. MdTPR16, an apple 835 tetratricopeptide repeat (TPR) -like superfamily gene, positively regulates drought 836 stress in apple. Plant Physiology and Biochemistry. 2024;210:108572. 837 70. Kruse E, Uehlein N, Kaldenhoff R. The aquaporins. Genome Biol. 2006;7(2):206. 838 71. Venkataraghavan A, Schwerdt JG, Tyerman SD, Hrmova M. Barley Nodulin 26 -839 like intrinsic protein permeates water, metalloids, sa ccharides, and ion pairs due 840 to structural plasticity and diversification. Journal of Biological Chemistry. 841 2023;299(12):105410. 842 72. Park JW, Benatti TR, Marconi T, Yu Q, Solis -Gracia N, Mora V, et al. Cold 843 responsive gene expression profiling of sugarcan e and Saccharum spontaneum 844 with functional analysis of a cold inducible Saccharum homolog of NOD26 -like 845 intrinsic protein to salt and water stress. PLoS One. 2015;10(5):e0125810. 846 73. Guo H, Guan Z, Liu Y , Chao K, Zhu Q, Zhou Y , et al. Comprehensive identification 847 and expression analyses of sugar transporter genes reveal the role of GmSTP22 848 in salt stress resistance in soybean. Plant Physiology and Biochemistry. 849 2024;216:109095. 850 74. Van Lijsebettens M, Vanderhaeghen R, De Block M, Bauw G, Villarroel R, Van 851 Montagu M. An S18 ribosomal protein gene copy at the Arabidopsis PFL locus 852 affects plant development by its specific expression in meristems. EMBO J. 853 1994;13(14):3378–88. 854 75. Ito T, Kim G, Shinozaki K. Disruption of an Arabidopsis cytoplasmic ribosomal 855 protein S13 ‐homologous gene by transposon ‐mediated mutagenesis causes 856 aberrant growth and development. The Plant Journal. 2000;22(3):257–64. 857 76. Sekula B, Ruszkowski M, Dauter Z. S -adenosylmethionine synthases in plants: 858 Structural characterization of type I and II isoenzymes from Arabidopsis thaliana 859 and Medicago truncatula. Int J Biol Macromol. 2020;151:554–65. 860 77. Bouvier F, Linka N, Isner JC, Mutterer J, Weber APM, Camara B. Arabidopsis 861 SAMT1 defines a plastid transporter regulating plastid biogenesis and plant 862 development. Plant Cell. 2006;18(11):3088–105. 863 78. Chen X, Lin WH, Wang Y , Luan S, Xue HW. An inositol polyphosphate 5 -864 phosphatase functions in PHOTOTROPIN1 signali ng in Arabidopis by altering 865 cytosolic Ca2+. Plant Cell. 2008;20(2):353–66. 866 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint 33 79. Mourtzinis S, Conley SP. Delineating soybean maturity groups across the United 867 States. Agron J. 2017;109(4):1397–403. 868 869 870 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted September 16, 2025. ; https://doi.org/10.1101/2025.09.10.675456doi: bioRxiv preprint