A natural population of 496 soybean accessions was used in this study, consisting of 280 landraces and 216 improved cultivars with a wide geographical distribution in China, and identified 9,357,842 single-nucleotide polymorphisms (SNPs). Principal component analysis of above SNPs showed a continuous distribution without distinct clusters. The linkage disequilibrium (LD) decreased to half of its maximum value at 455 kb for all samples (Zhang et al., 2021). The 148 historical soybean varieties were selected from 216 improved soybean lines, and origin, soybean-growing region and releasing years used in this study were collected on the basis of public registration information or personal communication with soybean breeders (detailed information was listed in Supplementary Table 1 ).

This experiment was carried out simultaneously at the Dishang experimental station of the Institute of Cereal and Oil Crops, Hebei Academy of Agricultural and Forestry Sciences (Shijiazhuang) and the crop experimental base of Huazhong Agricultural University (Wuhan) from 2016 to 2018 (detailed information was listed in Supplementary Table 2 ). Normal full seeds were selected for sowing in both environments, with three rows in each plot, each row being 1 meter long, with 0.4 meters between rows, and 0.1 meter between plants. Three repeats were set in each experimental environment using a randomized block design. At three weeks after sowing, the overall planting density was manually thinned to ~ 250,000 plants per hectare.

The same measurement methods and standards (as listed in Supplementary Table 3 ) were used to investigate agronomic traits under all experimental conditions. Twelve agronomic traits were investigated, comprising nine plant architecture-related traits and three seed-related traits. The plant architecture traits were measured after plant maturity, and seed traits were measured after harvest. Seed length, seed width, and hundred-seed weight were measured using mixed samples harvested in each plot. Ten healthy seeds were randomly selected and measured with a vernier caliper. Hundred-seed weight was measured by weighing 100 healthy seeds harvested in each plot using a balance. The 1^st branch angle, 2^nd branch angle, 3^rd branch angle, seed length, seed width, and hundred-seed weight were measured at both Shijiazhuang and Wuhan experimental stations only in 2016. Plant height, node number, lowest pod height, and branch number were measured in 2016, 2017, and 2018 at both experimental stations. Pod range was calculated by dividing plant height minus lowest pod height by plant height, and internode length was calculated by dividing plant height by node number. The R package “lme4” to was used to calculate the best linear unbiased estimator (BLUE) values and broad-sense heritability for phenotypic data from multiple environments for GWAS.

We estimated the variance of genotype, environment, and the interactions between genotype and environment for plant height (PH), node number (NN), internode length (IL), pod range (PR), lowest pod height (LPH), the 1^st branch angle (1-BA), the 2^nd branch angle (2-BA), the 3^rd branch angle (3-BA), and branch number (BN) by analysis of variance (ANOVA). We only estimated the variance of genotype and environment for hundred-seed weight (HSW), seed width (SW), and seed length (SL) for these traits using mixed samples grown in Shijiazhuang in 2016.

GWAS was performed based on phenotypic data for agronomic traits in plants from multiple environments (E1 to E6) and BLUE data for each trait, using the genotypic dataset as previously reported (Zhang et al., 2021). A multi-locus mixed-model method (MLMM) model in the R package “GAPIT3” was used as the GWAS model to detect minor-effect genomic regions, the kinship matrix, and the first ten principal components as co-variants (Segura et al., 2012; Wang and Zhang, 2021).

Significantly associated variants were binned together into peaks using a clumping method using the PLINK (Version 1.90) command plink –clump –clump-p1 1E-5 –clump-p2 1E-3 –clump-kb 455 for clump analysis as previously reported (Crowell et al., 2016). A clump was regarded as a singleton if fewer than five single nucleotide polymorphisms (SNPs) with –log10(p)<2.5 × 10^–4 were present within the clump interval, and the genome coordinate interval of each clump was defined as a GWAS QTL.

All accessible GWAS QTL and mapping QTL information were collected from SoyBase (https://www.soybase.org/). All QTL pairs with overlapping genome intervals were considered to be co-located QTL pairs.

To evaluate the direct genetic correlations of QTLs, all QTL pairs located on the same chromosome in the GWAS results with all phenotypic BLUE data were identified; the Inter-LD between QTLs was calculated using the following formula as previously reported (Fang et al., 2017):

,where aveLD (QTL1, QTL2) is the mean value of linkage disequilibrium (LD) between all variants in the interval of QTL1 and all variants in the interval of QTL2, and PmaxLD (QTL1) is the maximum value of the mean value of LD between any two variants in QTL1. PmaxLD (QTL2) is the maximum value of the average LD values of any two variants within QTL2. LD is the r² value between the specific SNPs calculated using PLINK. All QTL pairs with Inter-LD≥0.4 were screened. Cytoscape (Version 3.7.2) was used to construct the genetic association network in which agronomic traits and QTLs as nodes and QTL LD as edges are represented by dashed lines. Node diameter represents the minimum p-value of each variant in the QTL.

All genes in the overlapping region of GWAS results for branch angle related traits were extracted from the Williams 82 reference genome (https://phytozome-next.jgi.doe.gov/info/Gmax_Wm82_a2_v1). R function modified from R package IntAssoPlot (He et al., 2020) was used to visualize the regional GWAS results and the regional LD heatmap.

The physical positions of all genes in each clump interval in GWAS results were extracted from the Williams 82 reference genome. The genotypes of all variants within the interval of each gene of interest, including the 2000-bp upstream regions, were extracted from the soybean population genotype dataset which were generated by variant calling, imputation and annotation pipelines in a previous study, and all cultivars in this population were divided into independent haplotypes based on the different combinations of all variants in the interest region (Zhang et al., 2021). The differences of haplotypes were estimated by ANOVA or Wilcoxon test and visualized by R package ggpubr.

To functionally annotate these interest genes, the encoded protein sequences of these genes were extracted from Williams 82 reference genome, and these sequences were then aligned to Arabidopsis thaliana reference genome by BLAST to identified the homologous genes and extract the function of these genes.

In this study, we planted 496 soybean accessions in two regions (Wuhan and Shijiazhuang) of China during the summer over three consecutive years from 2016 to 2018 ( Supplementary Table 2 ). We recorded nine plant architecture traits: plant height (PH), node number (NN), pod range (PR), internode length (IL), lowest pod height (LPH), branch number (BN), the 1^st branch angle (1-BA), the 2^nd branch angle (2-BA), and the 3^rd branch angle (3-BA). We also recorded three seed-related traits: hundred-seed weight (HSW), seed width (SW), and seed length (SL) ( Supplementary Table 3 ). The coefficient of variation CV (%) for each of the 12 traits ranged from 5.9% (PR) to 35.4% (PH), highlighting their broad variation among this soybean core collection ( Supplementary Table 4 ). Compared to the three seed-related phenotypes, soybean architecture traits were more sensitive to environmental differences ( Supplementary Figure 1 ). For example, the CV for 3-BA was greater in Shijiazhuang (28.05%) than in Wuhan (16.80%), while the CV of HSW was 31.03% and 34.74% in Shijiazhuang and Wuhan, respectively ( Supplementary Table 4 ). The broad-sense heritability (H² ) of these traits varied from 0.27 (2-BA) to 0.94 (PH and NN). The H² values of the three seed-related traits were all greater than 0.8, while the three branch angle-related traits exhibited lower heritability values of 0.32, 0.27, and 0.32, respectively ( Supplementary Table 4 ). ANOVA indicated that both genotype and environment have significant effects on the 12 traits, and genotype-environment interactions were also highly significant for eight architecture traits with the PR being the sole exception ( Supplementary Table 4 ).

Analyzing the changes in agronomic traits of improved cultivars released in different years of soybean breeding history is important for providing directions for future breeding. We thus analyzed the changes for the 12 traits among 148 soybean varieties that were released from 1952 to 2015 ( Figure 1 ), 63% of which are grown in the Huang-Huai-Hai plain region in China ( Supplementary Table 1 ). In agreement with an earlier report, we determined that HSW has gradually been increasing (Qin et al., 2017). SW has also gradually been increasing, whereas IL and PR have gradually been decreasing, and the eight remaining traits were generally unchanged. These results indicated that seed-related traits have been the main focus of soybean improvement over the past ~ 60 years in China and that improving architecture traits has the potential for increasing yields in soybean breeding programs.

To identify major genetic loci that regulate these 12 traits, we conducted GWAS using the 496-soybean collection. Our analysis was based on the 5,294,054 markers in this natural population and the phenotypes of the 12 traits using seven types of data, comprising phenotypic data from two locations over three years (E1 to E6) and the corresponding best linear unbiased estimator (BLUE) data. Of the GWAS signals, 10,338 (–log₁₀ P ≥ 5) were significantly associated with the 12 traits (E1:1710, E2:1767, E3:1283, E4:740, E5:1042, E6:1224, BLUE:2572), and 188 signals were repeatedly observed in Wuhan (E1, E3, E5) and the Shijiazhuang (E2, E4, E6) ( Supplementary Table 5 ). 786 and 1786 were associated with seed and architecture traits in the BLUE dataset, respectively ( Figure 2 and Supplementary Figure 2 ). To identify the genomic regions defined for all these traits, we used a clumping-based method as previously described (Crowell et al., 2016). We successfully identified 169 genomic regions in the BLUE dataset, of which 129 were associated with architecture traits and 40 were associated with seed-related traits ( Table 1 and Supplementary Table 6 ). Furthermore, we compared all detected genomic regions with previously reported architecture- or yield-related genomic regions or QTLs detected by GWAS or bi-parental QTL mapping, respectively (https://www.soybase.org). We observed that 61 genomic regions co-localize with previously identified QTLs, with the remaining 108 genomic regions being newly identified. Among these genomic regions co-localizing with QTLs, 4 GWAS regions detected for plant height-related traits (PH, NN, IL and PR) overlapped with the previously characterized gene Dt1, and the GWAS regions harboring Dt2 locus were detected for 3 agronomic traits (PH, NN, and BN) (Ping et al., 2014; Fang et al., 2017; Liang et al., 2022). The identification of these overlapping regions (36% of total) indicated that our data collection and analysis were robust. In summary, we identified 108 new genomic regions associated with 12 soybean agronomic traits.

Correlations between traits is a common phenomenon. Thus, breeders must consider the correlation among traits in variety development, especially when the goal is to increase the values of two negatively correlated traits. We evaluated the phenotypic relationships among the 12 investigated agronomic traits and established that all 12 traits can be grouped into different categories ( Supplementary Figure 3 ). Seed traits were strongly and positively correlated (with correlation coefficients ranging from 0.56 to 0.96). For architecture traits, the three branch angle-related traits (1-BA, 2-BA, and 3-BA) were positively correlated (the correlation coefficients ranged from 0.66 to 0.78) and independent from the other traits (with correlation coefficients ranging from 0.01 to 0.08). The four plant height-related traits, PH, NN, IL and PR, were highly and positively correlated, while the same traits were negatively correlated with HSW and SW (with correlation coefficients ranging from -0.29 to -0.50) ( Supplementary Figure 3 ). These results suggested that regulating plant height-related traits might have negative effects on seed width and seed weight. By contrast, optimizing branch angle-related traits might have minor effects on seed-related traits in soybean.

The above results indicated that phenotypic correlations were common for soybean agronomic traits. In particular, Dt1 was associated not only with PH and NN, but also with IL and PR, and the Dt2 locus was simultaneously detected for PH, NN, and BN, suggesting that these related traits might be genetically co-regulated. Association networks allow regions of the genome associated with multiple phenotypes to be quantitatively identified. Therefore, we built a genetic association network by integrating the above GWAS results to explore the genetic relationship of the target traits. In total, we included 169 genomic regions defined by GWAS in this network, 90 of which were associated with a single agronomic trait, while the remaining 79 were linked to other regions by Inter-linkage disequilibrium (LD) or different traits. Of these 79 regions, 26 were associated with two traits, and 7 hub-nodes that linked to at least three traits were present in the multiple trait association network ( Figure 3 ). Notably, 12 and 6 GWAS regions were solely associated with IL and NN, respectively, suggesting that several relatively independent loci involved in regulating plant height-related traits were identified in this study ( Figure 3 ). qHub_2 is a newly identified hub-node linking three seed-related traits (HSW, SW and SL) and two plant height-related traits (PH and IL) ( Figure 3 and Supplementary Table 7 ). This genetic association might partially explain the significant phenotypic correlation between seed-related traits and plant height-related traits.

Dt1 and Dt2 are two key genes regulating plant height and node number in soybean, and they were presented in two hub-nodes, qHub_5 and qHub_3, respectively ( Figure 3 ). In agreement with previous reports, in this network, the Dt2 locus overlapped with a pleiotropic locus (qHub_3) associated with PH, NN, and BN (Ping et al., 2014; Fang et al., 2017; Liang et al., 2022). Dt1 was reported to affect plant height, node number, branch density, stem pod density, the number of three-seeded pods, and total seed number (Fang et al., 2017). In the association network, qHub_5 (containing Dt1) was linked to four plant height-related traits: PH, NN, IL and PR ( Figure 3 ). To investigate whether Dt1 might affect IL and PR, we conducted haplotype analysis for above plant height-related traits based on previously reported Dt1 haplotypes information (H1, H2, H3, H4) that were detected in 496 soybean accessions (Liu et al., 2010; Tian et al., 2010; Yue et al., 2021). Cultivars with the Dt1-H1 haplotype carrying the functional Dt1 allele were taller and produced more nodes than those carrying the dt1 haplotypes (H2, H3 and H4), as reported previously (Yue et al., 2021). Moreover, cultivars with the functional Dt1 allele had longer internodes and wider pod ranges than those with the dt1 alleles (H2, H3) ( Supplementary Figure 4 ). In addition, according to the released transcriptome data retrieved from the Plant Public RNA-seq Database (https://plantrnadb.com/), Dt1 is expressed in axillary meristems and internodes ( Supplementary Figure 5 ), suggesting that Dt1 might also affect soybean internode length in addition to node number and plant height.

Optimizing genes that regulate tiller/leaf angles has greatly contributed to boosting maize and rice yields under high-density planting (Li et al., 2007; Tian et al., 2019; Huang et al., 2021). To date, only a few studies have focused on genetically modifying soybean branch angle traits (Clark et al., 2021). Moreover, correlation analysis showed that branch angle traits were independent from other traits ( Supplementary Figure 3 ). In the genetic association network, three branch angle traits are associated with 16 genomic regions, and only one genomic region associated with 2-BA was linked to internode length by Inter-LD (value of Inter-LD is 0.6) ( Figure 3 and Figure 4A ). We also discovered a hub-node region (qHub_7) for branch angle traits. Consistent with our correlation analysis ( Supplementary Figure 3 ), the genetic association network indicated that branch angle traits were relatively independent and were strongly and positively correlated to each other as compared with the nine other agronomic traits ( Figure 3 ).

Since these three branching phenotypes were similar and closely and positively correlated ( Figure 4B ), we reasoned that candidate genes regulating these three traits are most likely present in qHub_7. The interval lengths of the GWAS regions for 1-BA, 2-BA, and 3-BA were 728 kb, 588 kb, and 706 kb, respectively ( Figures 4C–F ). We extracted the sequences of all genes in each region and performed haplotype analysis for the corresponding branch angles. We identified 16, 26, and 19 genes that were significantly associated with phenotypic differences among the haplotypes of these genes for 1-BA, 2-BA, and 3-BA, respectively. We considered the overlapping genes among these gene sets (eight genes) to be interesting candidates for branch angle-related traits ( Supplementary Table 8 and Supplementary Figure 6 ). We then searched a public gene expression database for the expression patterns of these candidate genes in different tissues and developmental stages (Shen et al., 2014). Only four genes (Glyma.09G227100, Glyma.09G228200, Glyma.09G228300, and Glyma.09G229000) were expressed, whereas the four other genes, whose expression was not detected, were excluded from further analysis ( Supplementary Figure 7 ).

To narrow down the list of candidate genes regulating branch angle, we examined data from functional studies of their homologs in Arabidopsis thaliana. The homologs of Glyma.09G228300 in Arabidopsis have no known functions and are annotated as encoding a protein with unknown function (DUF3411) and an uncharacterized protein. Glyma.09G227100 encodes a UvrABC system C protein involved in DNA damage repair. Glyma.09G229000 encodes a homolog of Arabidopsis actin-related protein 7 (ARP7) according to TAIR (https://www.arabidopsis.org/). ARP7 plays vital roles in embryogenesis and plant development, as knocking down this gene in Arabidopsis led to dwarfism, small rosette leaves, retarded root growth, altered flower development, delayed perianth abscission, and reduced fertility (Kandasamy et al., 2005). Intriguingly, Glyma.09G228200 encodes a protein homologous to the FK506-binding protein 42 (FKBP42) protein TWISTED DWARF 1 (ATWD1). TWD1 is involved in the auxin signal transduction pathway as well as brassinosteroid (BR) signaling (Wu et al., 2010; Zhao et al., 2016; Zhu et al., 2016). Both auxin and BRs play vital roles in modulating tiller or leaf angle in rice or maize (Sun et al., 2015; Tian et al., 2019; Wang et al., 2020b). The Arabidopsis twd1 mutant shows drastically reduced cell elongation combined with a disoriented growth behavior (Geisler et al., 2003). Thus, these results suggested that Glyma.09G228200 and Glyma.09G229000 are most likely involved in branch angle regulation. It will be necessary to genetically manipulate these two genes in soybean in order to verify their functions.