Here, a set of 340 diverse accessions of soybean derived from the different growing regions of China were used as plant material (Supplementary Table S1). These soybean accessions were planted under natural conditions at two separate locations of the Jilin province viz., Changchun (CC) and Jiutai (JT). In the CC location, the accessions were grown in the experimental field of the Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences (44°0’N, 125°24’E), on the following dates: April 29^th, 2022 (CC22); April 28^th, 2023 (CC23); and April 30^th, 2024 (CC24), respectively; whereas these accessions were grown on May 10^th, 2024 in the experimental field of Changchun Polytechnic (44°4’N, 125°40’E, Jiutai District, Changchun), representing the environment JT24. Hence, the diverse soybean germplasm resources were accessed across the four different environments viz., CC22, CC23, CC24 & JT24. The randomized complete block design (RCBD) involving three replications was used for planting this germplasm. The plot of each soybean accession was consisted of three rows with the spacing between the rows maintained at 50 cm and the length of row kept at 200 cm. At the maturity stage, the seeds of each accession were harvested in the main season of all four environments. The soybean germplasm was cultivated following the normal agronomic practices.

The seed samples collected were placed in a cup with a bottom of 5 cm in diameter, and a height of 3 cm for near-infrared spectroscopy analysis; and the seeds were analyzed without any treatment and in good condition. The cup was filled with the equivalent of about 60 ml of soybean seed (about 20 g or 80 seeds). Spectra acquisition was carried out using a Fourier transform near-infrared (FTNIR) analyzer (Antaris II spectrometer: Thermofisher Scientific, France). The SOC and SPC were expressed in percentage by weight of seeds. Phenotypic data of SOC and SPC were estimated for each accession in three replications for all four environments. The SOC and SPC were estimated in each replication, and the average of three replications was used for further analysis. The combined environment (CE) data were derived from the data of individual environments by utilizing the “lme4” package in R environment (Bates et al., 2015). The predicted means (best linear unbiased predictions, BLUPs) for the CE were calculated by fitting the following model:

where Yijk is the phenotypic value of the i-th genotype in the k-th replication of the j-th environment, μ is the overall mean, Gi is the random effect of the i-th genotype, Ej is the random effect of the j-th environment, GEij is the random effect of the genotype-by-environment interaction, Repk(Ej) is the random effect of the block nested within the j-th environment, and ϵijk is the random error. The analysis of variance (ANOVA) was calculated by utilizing the R language. The broad-sense heritability (H) was estimated as σG2/(σG2+σGE2/n+σϵ2/nr), in which the σG2, σGE2 and σϵ2 represent the genotypic variance, genotype × environment interaction variance, and residual variance, respectively, n is the number of environments, and r is the number of replications within each environment.

The resequencing data for the soybean accessions were obtained from our previously published study (Zhang et al., 2024b). Briefly, the genomic DNA was extracted from fresh and healthy leaf tissues of three-week-old seedlings using the cetyltrimethylammonium bromide (CTAB) method (Murray and Thompson, 1980). The genomic DNA was fragmented to a target size of 350 bp. Library preparation was performed by Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China) sequenced on the DNBSEQ-T7 platform (BGI, Shenzhen, China) by Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China). The average depth of resequencing was 25 × and the volume of the total data produced was 8.57 Tb. The fastp software (Chen et al., 2018) was utilized to eliminate the reads with low quality, reads with > 10% N content and adapter sequences; and we retained only the clean reads. The Burrows-Wheeler Aligner software (Li, 2013) was utilized to map the sequence reads with the reference genome of Wm82.a2.v1. The Genome Analysis Toolkit (McKenna et al., 2010) was used to identify the sequence variants such as SNPs, insertions/deletions (InDels) and other larger variants; and in this way we identified a total of 25,442,248 genetic variants. Furthermore, we utilized the stringent measures of quality control viz., missing rate greater than 10% and minor allele frequency less than 0.05, and finally detected high-quality SNPs of 3,343,372, with the mean heterozygous rate of 0.094.

The PCA result and kinship of GWAS panel were generated by GAPIT package (Wang and Zhang, 2021) in R environment. The PopLDdecay software v3.43 (Zhang et al., 2018a) was used for LD decay calculation with the following parameters: -MaxDist 300, -MAF 0.05, and -Miss 0.1.

For association mapping analysis, the study utilized the GAPIT package (Wang and Zhang, 2021). Seven different GWAS models were used viz., GLM (General Linear Model), MLM (Mixed Linear Model), CMLM (Compressed Mixed Linear Model), MLMM (Multi-Locus Mixed Linear Model), SUPER (Settlement of MLM Under Progressively Exclusive Relationship), FarmCPU (Fixed and random model Circulating Probability Unification) and BLINK (Bayesian-information and Linkage-disequilibrium Iteratively Nested Keyway). The multiple models were used because various models of GWAS are built on divergent hypotheses involving different distribution characteristics of QTL effects. Among these, MLMM, BLINK, SUPER, and FarmCPU are considered multi-locus models, whereas GLM, MLM and CMLM are categorized as single-locus GWAS models. The default parameters of GAPIT were employed to detect significant SNP associations. The GAPIT package in the R program was employed for conducting principal component analysis (PCA) and kinship analysis to assess population structure (Lipka et al., 2012). To visualize population structure in R environment, we used the “plotly” package. A suggestive association threshold (i.e., 1/N_test) was set at P < 2.99 × 10^-7, which reflects the level at which, under the null hypothesis, one false positive result is expected across the entire genome scan (Lander and Kruglyak, 1995; Duggal et al., 2008; Wang et al., 2012; Li et al., 2013; Guo et al., 2017). The UpSetR package was used to draw the UpSet plots (Conway et al., 2017).

To determine the level of LD between pairs of SNPs, we employed the Haploview 4.2 software (Barrett et al., 2004). A haplotype block was determined when two or more than 2 significant SNPs were present in the same chromosome, and these SNPs were present within the LD decay distance in the upstream and downstream of chromosome such as ± 71.7 kb in the present study. Pairwise means were analyzed using Fisher’s Least Significant Difference (LSD) Test, and visualization of the results was performed in the R environment.

Gene models that were present across the genomic interval of 17 QTLs along with their function annotation were collected from SoyBase (https://www.soybase.org/); and soybean genome version of Wm82.a2.v1 was utilized to collect the gene model information. The information about gene ontology (GO) was obtained from SoyBase (https://www.soybase.org/). The analysis of GO enrichment was carried out by utilizing the agriGO V2.0 (https://systemsbiology.cau.edu.cn/agriGOv2/) tool (Tian et al., 2017). The publicly available RNA-seq data at SoyMD (https://yanglab.hzau.edu.cn/SoyMD/#/) for different tissues/organs as well as development stages were collected for the model genes (Yang et al., 2023). The heatmap of transcriptomic data was produced by utilizing the ComplexHeatmap package of R software (Gu et al., 2016). Initial screening of the putative candidate genes was carried out using the in-silico analysis and a literature survey. The genes selected based on the initial step were further screened using the variant annotation and haplotype analysis; and the genes fulfilling the criteria of variant annotation and haplotype differences were prioritized as potential candidate genes regulating the SOC and SPC in soybean. The variant annotation was carried out by utilizing the SnpEff software (Cingolani et al., 2012).

The sequencing data for the haplotype selection analysis of soybean accessions adapted to different latitudes of China has been collected from SoyMD (https://yanglab.hzau.edu.cn/SoyMD/#/). The different regions of China have been divided into China-I, China-II, China-III, China-IV, China-V and China-VI (Liu et al., 2023); China-I has the highest latitude including the Northeast China; China-II has the mid-latitude that includes Huanghuaihai regions; and China-III to China-VI include the Southern China and have low-latitude. Pairwise F_ST values were estimated by using the window size of 5 kb and a step size of 200 bp using VCFtools (Danecek et al., 2011). Nucleotide diversity (π) values within various populations were also evaluated using VCFtools (Danecek et al., 2011) with a fixed window of 5 kb and a step size of 200 bp. Tajima’s D values were evaluated using VCF-kit (Cook and Andersen, 2017) with a fixed window of 5 kb and a step size of 200 bp.

The mean values of the SOC and SPC across the four different environments, viz. CC22, CC23, CC24 and JT24 are revealed in the violin plots (Figures 1A, B; Supplementary Figure 1). The average values of SOC showed significant differences among the four environments; however, the average values of SOC exhibited non-significant differences between CC22 and CC23, while the remaining comparison groups of environments showed significant differences in the average values of SOC (Figure 1A). The average values of SPC revealed significant differences among the four environments, but the non-significant differences were observed between CC22 and JT24 for the SPC; while the significant differences were observed for the average values of SPC in the remaining comparison groups (Figure 1B). The frequency distributions of SOC and SPC traits in the set of 340 diverse soybean accessions across four individual environments, viz. CC22, CC23, CC24 & JT24, are shown in Supplementary Figure 1. The scatter plots of SOC and SPC traits in CC22, CC23, CC24, JT24 environments and the CE shown in Figure 1C and the Pearson correlation coefficients between these traits in different individual environments and the CE (-0.41, -0.46, -0.41, -0.36 and -0.34, P < 0.0001) underscored a consistent moderate negative correlation between SOC and SPC across different environments in our soybean population.

The parameters of descriptive statistics for the SOC and SPC traits in 340 soybean cultivars across four individual environments, viz. CC22, CC23, CC24 & JT24 plus the CE, are presented in Supplementary Table S2. The average values of the SOC ranged from 19.21 to 20.56% in the environments of CC23 and JT24, respectively; and the minimum and maximum values for the SOC were 13.59% and 24.67% in the environments of CC22 and JT24, respectively. The coefficients of variation (CVs) of SOC varied from 7.62 to 8.37% in the environments of JT24 and CC23, respectively; and the CV of SOC in the CE was 7.11%. The skewness of the SOC varied from -0.79 to -0.21 in the JT24 and CC23 environments, respectively; and the kurtosis of SOC ranged from 0.52 to 1.37 in the environments of CC23 and CE, respectively. Broad-sense heritability (H) for SOC was high equal to 0.86 (Supplementary Table S2). The results of ANOVA analysis showed highly significant (P<0.0001) variances of genotype (G), environment (E) and genotype × environment (G × E) interaction for the SOC (Table 1).

The average values of the SPC ranged from 39.75 to 42.23% in the environments of CC22 and CC23, respectively; and the minimum and maximum values of SPC were 31.72% and 50.83% in the environments of CC22 and CC23, respectively. The CVs for the SPC varied from 4.58 to 5.59% in the environments of CC24 and CC22, respectively; and the CV of SPC in the CE was 4.36%. The skewness for the SPC ranged from -0.11 to 0.02 in the CC22 and CC23 environments, respectively; and the kurtosis of SPC varied from 0.02 to 0.50 in the environments of CC24 and CC23, respectively. Broad-sense heritability (H) for SPC was high, and was equal to 0.79 (Supplementary Table S2). The results of ANOVA analysis showed highly significant (P<0.0001) variances of G, E and G × E interaction for the SPC (Table 1).

Here, we detected 3,343,372 high-quality SNPs by applying the stringent measures of quality control (Supplementary Table S3). The distribution of these SNPs was located across the whole genome of the soybean; the lowest and highest SNPs were present on Chr.11 (63,224) and Chr.15 (256,000), respectively, with an average number of 167,168.6 SNPs per chromosome (Figure 2B; Supplementary Table S3). Substantial variation was observed for the marker density among the various soybean chromosomes; for example, the lowest marker density of 1,818.51 SNPs/Mb was observed on Chr.11, and the highest marker density of 4,946.25 SNPs/Mb was observed on Chr.15, and the mean value of the marker density was 3,522.37 SNPs/Mb (Figure 2A; Supplementary Table S3). The kinship matrix showed no clear-cut clustering in the panel of soybean accessions (Figure 2E). Besides, no distinct structure was observed, and a continuous distribution of the soybean accessions has been shown by the analysis of population structure (Figure 2D; Supplementary Table S3).

We studied the features of LD across the panel of 340 diverse soybean accessions, and the outcomes are revealed in the graph (Figure 2C). The average value of r² was 0.32, with LD decay started at 0.68; the half decay reached at 0.34. The intersection occurs between the LD decay curve and half decay at the 71.7 kb; and this distance decides the linkage at the genome-wide level. The genomic region of ± 71.7 kb across the stable significant SNPs was designated as the QTL interval.

Here, we detected 105 and 59 significant SNPs linked with the SOC and SPC traits, respectively across the four individual environments (CC22, CC23, CC24 & JT24) plus CE and seven GWAS models (Figures 3, 4; Supplementary Table S4). The 105 significant SNPs of SOC were distributed across all the soybean chromosomes except Chr.07 and Chr.15; and the allocation of these SNPs in the various models and environments were presented in UpSet plot (Figures 5A, B). A maximum of 33 significant SNPs of SOC were present on Chr.08, and a minimum of one significant SNP was identified on Chr.16 and Chr. 20. The remaining 15 chromosomes possessed significant SNPs in the range of 2 to 15 (Supplementary Table S4). By considering the detection of SOC related significant SNPs by various GWAS models, the one SNP viz., Chr08_9071263 was persistently identified by all the seven models viz., GLM, MLM, CMLM, MLMM, SUPER, FarmCPU and BLINK; the rest of the 104 SOC related significant SNPs were identified by one to five GWAS models (Supplementary Table S4). Furthermore, the SNP Chr08_9071263 was persistently identified by three single environments plus CE. The SNPs viz., Chr08_9054741, Chr08_9056324 and Chr08_9071236 were persistently detected in two single environments plus CE; and Chr09_6929070, Chr10_44584365, Chr10_44584422 and Chr15_10147967 were detected by two single environments or one single environment plus CE. The rest of the 97 SOC related significant SNPs were identified in one single environment or CE (Supplementary Table S4).

The 59 significant SNPs detected for the SPC were distributed across the 17 out of total 20 soybean chromosomes viz., Chr.01, Chr.02, Chr.03, Chr.05, Chr.06, Chr.07, Chr.08, Chr.09, Chr.10, Chr.12, Chr.13, Chr.14, Chr.15, Chr.16, Chr.17, Chr.18 and Chr.20; no significant SNPs for the SPC were detected on the Chr.04, Chr.11 and Chr.19 (Supplementary Table S4). Distribution of these 59 SNPs across the different GWAS models and environments is presented in the UpSet plot (Figures 5C, D). We detected a maximum of eight SNPs on the Chr.18, and a minimum of one significant SNP on the chromosomes viz., Chr.07, Chr.10 and Chr.14; and the remaining 14 chromosomes possess the significant SNPs from two to six (Supplementary Table S4). The significant SNPs for the SPC were distributed across the different models as, the SNPs viz., Chr08_9054741 and Chr18_54408341 were detected by all the seven GWAS models, this is followed by SNP Chr18_54376853 that was detected by five models viz., GLM, MLMM, super, FarmCPU and BLINK; the remaining 56 SNPs were detected by one to three GWAS models (Supplementary Table S4). Moreover, the SNPs viz., Chr18_54408341 was detected by two single environments (CC22 & CC23) plus CE; and the SNPs viz., Chr01_11010539, Chr06_6762512 and Chr08_9054741 were detected by one single environment plus CE. The rest of the 55 SPC related SNPs were identified in one single environment or CE (Supplementary Table S4).

However, among the total of 105 and 59 SNPs detected to be associated with SOC and SPC, respectively; only three significant SNPs viz., Chr08_9054741, Chr08_9071236 and Chr08_9071263 were associated with both SOC and SPC (Supplementary Table S4). The effect of Chr08_9054741 was negative on both SOC (-1.34) and SPC (-2.04) traits, whereas the effect of Chr08_9071236 was positive on both SOC (1.35) and SPC (1.76); moreover, the effect of Chr08_9071263 was positive on SOC (1.32) as well as SPC (1.94). These results indicate that these three common SNPs govern the SOC and SPC traits in the same direction. The remaining SNPs were associated with either SOC or SPC (Supplementary Table S4).

The significant SNPs present in the same chromosome that were associated with the same trait or related traits, and were fitting into the LD decay range of ± 71.7 kb formed the haplotype blocks. Besides, the SNPs identified through different GWAS models but fell within the physical interval of ± 71.7 kb were also used for haplotype block construction. In this regard, our findings revealed that two, three, two, two, three, three, twenty-six, three, eleven, two, two, two, five, two and two significant SNPs associated with the SOC or/and SPC were present on the chromosomes viz., Chr.01, Chr.06, Chr.06, Chr.08, Chr.08, Chr.08, Chr.08, Chr.09, Chr.10, Chr.12, Chr.15, Chr.15, Chr.18, Chr.18 and Chr.20, respectively and are fitting into the half LD decay distance; therefore, they constituted 15 haplotype blocks viz., Hap-1, Hap-6.1, Hap-6.2, Hap-8.1, Hap-8.2, Hap-8.3, Hap-8.4, Hap-9, Hap-10, Hap-12, Hap-15.1, Hap-15.2, Hap-18.1, Hap-18.2 and Hap-20, respectively (Figures 6, 7; Supplementary Table S5).

The number of haplotype alleles of these 15 haplotype blocks varied from two (Hap-1) to four (Hap-6.1) (Figures 6, 7; Supplementary Table S5). Out of 15 haplotype blocks, the haplotype alleles of 14 blocks viz., Hap-1, Hap-6.1, Hap-6.2, Hap-8.1, Hap-8.2, Hap-8.3, Hap-8.4, Hap-9, Hap-10, Hap-15.1, Hap-15.2, Hap-18.1, Hap-18.2 and Hap-20 revealed differences at the significant level in the governing of SOC that varied from 16.68-21.15% (Figure 6; Supplementary Table S5). In contrast, the haplotype alleles of 11 out of total 16 haplotype blocks viz., Hap-1, Hap-6.1, Hap-6.2, Hap-8.1, Hap-8.2, Hap-8.3, Hap-8.4, Hap-12, Hap-15.2, Hap-18.1, and Hap-20 revealed differences at the significant levels in the governing of SPC that varied from 38.63-42.69%; whereas the haplotype alleles of the remaining four haplotype blocks viz., Hap-9, Hap-10, Hap-15.2, and Hap-18.1 showed differences at the non-significant levels in the governing of SPC (Figure 7; Supplementary Table S5).

In conclusion, the haplotype alleles of 14 and 11 haplotype blocks associated with SOC and SPC, respectively, regulate the different phenotypes of these two traits ranging from low-intermediate-high.

Multiple significant SNPs were detected in different environments and models across the LD decay (± 71.7 kb) distance of 15 haplotype blocks viz., Hap-1, Hap-6.1, Hap-6.2, Hap-8.1, Hap-8.2, Hap-8.3, Hap-8.4, Hap-9, Hap-10, Hap-12, Hap-15.1, Hap-15.2, Hap-18.1, Hap-18.2 and Hap-20, hence they also represented the 15 stable QTLs named as qSOC1, qSOC_SPC6.1, qSOC_SPC6.2, qSOC8.1, qSOC8.2, qSOC8.3, qSOC_SPC8, qSOC9, qSOC10, qSPC12, qSOC_SPC15.1, qSOC_SPC15.2, qSPC18, qSOC18 and qSPC20 (Supplementary Table S4; Supplementary Figure 2). Besides, the genomic intervals of ± half LD-decay distance across the SNPs detected by various models and environments such as Chr01_11010539 and Chr15_10147967, which were present on Chr.01 and Chr.15, respectively, were defined as stable QTLs. They represent two QTLs referred to as qSPC1 and qSOC15, respectively (Supplementary Table S4). Hence, together we identified a total of 17 QTLs; among these QTLs, the eight viz., qSOC1, qSOC8.1, qSOC8.2, qSOC8.3, qSOC9, qSOC10, qSOC15 and qSOC18 were associated with the SOC trait only, whereas the four QTLs viz., qSPC1, qSPC12, qSPC18, and qSPC20 were associated with only SPC trait. However, the five QTLs viz., qSOC_SPC6.1, qSOC_SPC6.2, qSOC_SPC8, qSOC_SPC15.1 and qSOC_SPC15.2 were linked with both SOC and SPC traits (Supplementary Table S4; Supplementary Figure 2). The above genomic loci/QTLs were identified through various environments and GWAS models, that in turn indicates the persistency and stability of these genomic loci in the governing of SOC or/and SPC.

To quantify the genetic basis underlying the observed phenotypic correlation between SOC and SPC, we performed a linear regression analysis to evaluate the phenotypic variance explained (PVE) by the identified stable QTLs based on the phenotypic data from the CE (Supplementary Table S6). The five major pleiotropic QTL regions (qSOC_SPC6.1, qSOC_SPC6.2, qSOC_SPC8, qSOC_SPC15.1 and qSOC_SPC15.2) were found to collectively account for 26.19% of the total variance for SOC and 22.26% for SPC. This indicated that around a quarter portion of the variation in both traits was governed by shared genomic “hotspots”. Furthermore, the other six stable QTLs of which the haplotype blocks exhibited significant individual effects on both traits (qSOC1, qSOC8.1, qSOC8.2, qSOC8.3, qSPC18 and qSPC20) explained a total PVE of 22.17% for SOC and 29.49% for SPC. These 11 QTLs explained 48.36% of the phenotypic variance of SOC and 51.74% of SPC, respectively. The remaining phenotypic residuals (~31%) likely reflected the influence of numerous sub-threshold minor-effect loci and environmental factors. These results provide a genetic explanation for the significant correlation between SOC and SPC observed in our germplasm.

Here, we identified 353 model genes across the genomic intervals of 17 QTLs viz., qSPC1, qSOC1, qSOC_SPC6.1, qSOC_SPC6.2, qSOC8.1, qSOC8.2, qSOC8.3, qSOC_SPC8, qSOC9, qSOC10, qSPC12, qSOC15, qSOC_SPC15.1, qSOC_SPC15.2, qSPC18, qSOC18 and qSPC20, that include 5, 15, 28, 24, 21, 19, 25, 41, 9, 30, 37, 13, 11, 22, 25, 14 and 14 genes, respectively. The gene ontology (GO) analysis was carried out by utilizing the online available web-based toolkit viz., agriGO V2.0; and this analysis showed the enrichment of various terms among the major categories viz., molecular function, biological process and cellular component within the 17 genomic loci linked with SOC and SPC traits (Figure 8A). The terms with the highest enrichment across the 17 QTLs are shown in Figure 8A. The gene expression levels of model genes residing in the physical intervals of 17 QTLs across the various tissues/organs and seed development stages of soybean are presented in the form of heatmap (Figure 8B). By using the in-silico analysis viz., function annotations, transcriptomic analysis and literature survey, we sorted out 121 genes across the genomic interval of 17 QTLs as possible candidates governing SOC and SPC traits in soybean (Supplementary Table S7). These include 3, 5, 11, 9, 11, 7, 7, 14, 3, 6, 5, 6, 4, 4, 8, 10, and 8 candidate genes across the physical interval of 17 QTLs viz., qSPC1, qSOC1, qSOC_SPC6.1, qSOC_SPC6.2, qSOC8.1, qSOC8.2, qSOC8.3, qSOC_SPC8, qSOC9, qSOC10, qSPC12, qSOC15, qSOC_SPC15.1, qSOC_SPC15.2, qSPC18, qSOC18 and qSPC20, respectively (Supplementary Table S7). Gene function annotations such as amino acid biosynthesis, protein biosynthesis, amino acid and protein transport, fatty acid biosynthesis, lipid biosynthesis, fatty acid and lipid transport as well as processes related to the protein and lipid biosynthesis and their accumulation, and seed development (Supplementary Table S7), and genes showing higher expression (as showed by RNA-seq dataset) during the seed development stages were the main criteria to select the candidate genes (Figure 8B). These 121 genes were regarded as putative candidates for regulating the SOC and SPC in soybean.

To further screen out the genes regulating the SOC and SPC traits, we performed the sequence variation analysis, variant annotation analysis and haplotype analysis. Out of these 121 genes, only 104 genes showed sequence variations (Supplementary Table S8). Among these 104 genes, the haplotype alleles of 80 genes revealed differences at significant levels in governing the SOC or/and SPC (Supplementary Tables S9). In these 80 genes, the haplotype allele number varied from a minimum of two to the maximum of 8, regulating the SOC and SPC in the ranges of 15.98-21.23% and 37.69%-43.30%, respectively. Among these 80 genes, the haplotype alleles of 11 genes only showed differences at the significant levels in governing the SOC, whereas the haplotype alleles of 6 genes only showed differences at the significant levels in governing the SPC. However, in the remaining 63 genes (out of the total 80) the haplotype alleles revealed differences at the significant levels in governing of both SOC and SPC. The haplotype alleles of these 80 genes regulate the different phenotypes of SOC or SPC or both traits ranging from minimum to maximum levels through intermediate phenotypes. The variant annotation analysis showed that 45 of the 80 genes showed non-synonymous mutations, whereas the remaining 35 genes showed synonymous mutations (Supplementary Table S8). Hence, these 80 genes were prioritized as candidate genes governing the SOC and SPC traits in soybean (Supplementary Table S9).

Furthermore, we performed haplotype selection analysis of the soybean accessions adapted to different latitudes i.e., from low-latitude to high-latitude for these 80 genes using the public data from the SoyMD (https://yanglab.hzau.edu.cn/SoyMD/#/; Supplementary Figures 3-8). Our haplotype results showed that 12 genes viz., Glyma.01G175500, Glyma.01G176000, Glyma.08G101100, Glyma.08G101900, Glyma.08G103100, Glyma.08G111700, Glyma.08G111800, Glyma.08G117000, Glyma.15G232200, Glyma.18G258000, Glyma.18G258100, Glyma.18G259100 underwent divergent selection for the specific haplotypes across the soybean population adapted to different latitudes (Table 2; Figure 9). For example, in Glyma.01G175500, the Hap_001 was predominantly selected in soybeans in Northeast China, whereas Hap_002 was predominantly selected in soybeans from southern parts of China. In the same way, the divergent specific haplotypes have been selected for the low- and high-latitude soybeans in the remaining 11 genes (Table 2). Moreover, the F_ST, π and Tajima’s D analysis also revealed the differential/divergent selection of these 12 genes among the soybean populations adapted to different latitudes (Supplementary Figures 3-8). Hence, these 12 genes are the most potential candidates regulating the difference in the SOC and SPC traits in the soybeans adapted to the environment of Northeast China.

Here, a set of 340 diverse accessions of soybean was assessed for SOC and SPC traits at two locations viz., Changchun and Jiutai research farms, corresponding to the entire of four environments viz., CC22, CC23, CC24 and JT24 plus CE. The effects of G, E and G × E interactions on the SOC and SPC traits were highly significant, which is in accordance with the findings documented by earlier authors (Li et al., 2019; Jin et al., 2023). These results in turn indicated the great potential of this soybean germplasm for the improvement of these traits; besides, the differences in the environmental variables among various environments have resulted in the different trait values of the SOC and SPC in the same set of soybean accessions. The higher value of H for the SOC and SPC traits indicated the same accessions of GWAS panel will perform alike if the same environmental conditions were provided. In agreement, many previous studies have reported similar outcomes for the SOC and SPC traits in soybean (Li et al., 2019; Diers et al., 2023; Kim et al., 2023). Moreover, highly significant G, E and G × E interactions shown by SOC and SPC revealed their complex pattern of inheritance; and similar outcomes are reported by earlier researchers (Jin et al., 2023).

In recent years, GWAS has been routinely used in plant research for the elucidation of genetic basis of different traits; and the special thanks go to the progress done in genome sequencing and high-throughput genotyping (Alqudah et al., 2020). In soybean, GWAS is utilized to explore the genetic elements of multiple traits including the SOC and SPC (Lee et al., 2019; Li et al., 2019; Jin et al., 2023). Hence, the current study used the advanced SNP genotyping in the GWAS approach to unravel the genomic loci and genes for SOC and SPC traits in soybean. Here, we detected the significant association of 105 and 59 SNPs with the SOC and SPC, respectively across the multiple environments (CC22, CC23, CC24 & JT24) plus CE and seven GWAS models. These 105 and 59 significant SNPs of SOC and SPC were located on the 18 and 17 chromosomes of soybean, respectively, indicating multigenic inheritance of both traits, i.e. regulated by multiple genes. In accordance with our findings, the previous research studies have also documented the similar results of SOC and SPC traits in soybean (Jin et al., 2023; Doszhanova et al., 2024). Among the 105 and 59 SNPs detected significantly associated with SOC and SPC traits, respectively; only three common significant SNPs were linked with both SOC and SPC viz., Chr08_9054741, Chr08_9071236 and Chr08_9071263, while the remaining significant SNPs were either associated with the SOC or SPC. The three common SNPs associated with both SOC and SPC traits indicated their pleiotropic effect on the two traits, that may be because these SNPs were associated with genes encoding the enzyme involved in the common metabolic pathways (Gao et al., 2024; Van et al., 2025).

Furthermore, our outcomes revealed substantial variation in the identification of significant SNPs by seven different models for both traits. For instance, the Chr08_9071263 linked with SOC was detected by all seven models; in comparison, only one GWAS model detected the 30 significant SNPs. Similarly, the significant SNPs viz., Chr08_9054741 and Chr18_54408341 of SPC were detected by all seven GWAS models, whereas 27 SNPs of SPC were detected by a single GWAS model. Similar results have been reported by the earlier studies for various crop plants including soybean (Leventhal et al., 2025), maize (Kaler et al., 2020) and wheat (Merrick et al., 2022). The reason to describe this phenomenon is because the various models of GWAS are built on divergent hypotheses involving different distribution characteristics of QTL effects; hence, this leads to the difference in SNP detection among various GWAS models (Odell et al., 2022).

The moderate negative phenotypic correlation between SOC and SPC might be explained at the genetic level because some of the SNPs, QTLs and haplotype blocks affecting both traits simultaneously. Besides, the effects of environment as well as G × E were highly significant in the GWAS population that might have further resulted in the negative correlation. The correlation between SOC and SPC was underpinned by a substantial shared genetic basis, as evidenced by the collective PVE of approximately 26% and 22% by the major pleiotropic QTLs, respectively (Supplementary Table S6). Considering the other six stable QTLs of which the haplotype blocks exhibited significant individual effects on both traits, the 11 detected QTLs explained 56.23% of the heritability of SOC and 65.50% of SPC, respectively. Our regression analysis showed that the detected stable QTLs have already explained more than half of the genetic commonality. This implies that the observed phenotypic linkage is not merely a localized phenomenon but the aggregate result of multiple stable genetic regions.

Besides, substantial variation has been observed for the SNPs detected for both traits across the different environments. In some genomic regions, a locus may surpass the significance threshold for one trait/environment while having a significance level below the threshold for the other one. These suggestive signals may show consistent peak patterns across both traits/environments, supporting shared genetic control, but might be masked by statistical corrections and G × E interactions. These outcomes revealed that SOC and SPC traits are substantially influenced by the environment, hence supporting the highly significant G × E reported in this study.

In the GWAS analysis, soybean researchers have mostly used biallelic SNP markers; and these biallelic markers are not able to detect the multiple alleles present across the diverse GWAS panel that might represent the superior and rare alleles, hence making this valuable genetic diversity unavailable for the crop improvement (He and Gai, 2023). Therefore, there is an important need for multiallelic markers in the GWAS analysis to trap these superior/rare alleles as well as non-linear variation across the diverse natural crop germplasm. In this context, multiallelic markers of haplotypes have promised to fix these challenges of biallelic SNP markers in the GWAS (Bhat et al., 2021; Voorrips and Tumino, 2022). Hence, these markers will make the considerable genetic diversity accessible for crop improvement (Voorrips and Tumino, 2022). In the last few years, some studies have identified the superior haplotypes for various traits in soybean such as traits related to yield (Shao et al., 2022), plant height (Qin et al., 2023) and salinity tolerance (Patil et al., 2016). Besides, superior haplotypes are also detected in other crop species such as in pigeonpea for water deficit tolerance (Sinha et al., 2020) and in rice for quality traits (Abbai et al., 2019). There are multiple mechanisms including selection, recombination and mutation that govern the haplotype variation across the diverse crop germplasm (Hellsten et al., 2013). Hence, identifying and utilizing the haplotypes in molecular breeding has a substantial impact in crop improvement (Sinha et al., 2020). Here, we identified a total of 15 haplotype blocks, and the haplotype alleles of 14 out of total 15 haplotype blocks showed differences at the significant levels in the governing of SOC, whereas among the 15 haplotype blocks, the haplotype alleles of 11 blocks revealed differences at the significant levels in the governing of SPC. The haplotype alleles of these 15 blocks vary from two to four regulating the SOC or SPC from the low-intermediate-high levels. Among the 15 haplotype blocks, the haplotype alleles of Hap-1, Hap-8.2, Hap-8.3, Hap-8.4, showed positive correlation in the regulation of SOC and SPC; whereas the haplotype alleles of Hap-6.2, Hap-15.2 and Hap-20 showed negative correlation in the regulation of the SOC and SPC traits. However, the haplotype alleles of Hap-9, Hap-10, Hap-15.1, Hap-18.2 showed non-significant differences on the regulation of SPC, but revealed differences at the significant level in the governing of SOC; whereas the haplotype alleles of Hap-12 revealed differences at the significant level in the governing of SPC but revealed differences at the non-significant level in the governing of SOC. These findings suggest the SOC and SPC traits of soybean possess complicated genetic structure, where some genetic loci have positive effects on both SOC and SPC traits, and some genetic loci showed negative correlations in the regulation of SOC and SPC; and some loci have independent effects on either SOC or SPC. Our study showed that alleles of each haplotype block regulate the varied phenotypes of SOC and SPC, thus delivering an opportunity to optimize the desired quantity of SOC and SPC in the soybean cultivar. These haplotypes can be used to produce cultivars of soybean with optimized as well as improved SOC and SPC, hence will substantially impact the quality and yield of soybeans.

By considering the genomic regions (± half LD-decay distance) flanking the two significant SNPs detected through various models and environments as well as the LD decay distance across the 15 identified haplotype blocks, we delineated a total of 17 stable QTLs associated with SOC or SPC or both traits. Among these 17 QTLs, 11 of them fell within the genomic interval of previously reported QTLs, hence not being considered as novel QTLs (Supplementary Figure 2). For example, Zhang et al. (2025) identified the significant SNP associated with SPC on Chr.06, and that SNP fell within the genomic interval of qSOC_SPC6.1 (5133255–5365962 bp). The qSOC_SPC6.2 present on Chr.06 was present in the genomic regions of earlier reported QTLs viz., cqSeed protein-015 (5449370–7326665 bp) and cqSeed oil-016 (5504147–8260733 bp) (Pathan et al., 2013). Hence, the qSOC_SPC6.2 might belong to the same cqSeed protein-015 and cqSeed oil-016 QTLs, as previously reported by Pathan et al. (2013). The qSOC8.1 fell within the genomic region (7525742- 9439751) of previously reported QTL viz., Seed oil 11-g2 (Yao et al., 2020). Similarly, qSOC8.2 fell within the genomic regions (7966391- 11892640) of earlier reported QTL viz., Seed oil 1-1 (Mansur et al., 1993). The physical interval of qSOC8.3 was located within the genomic region of earlier reported QTL viz., Seed oil 43-1 (Mao et al., 2013). The qSOC_SPC8 present on the Chr.08 fell within the physical intervals of earlier reported QTLs viz., Seed oil 8-g13 (Zhang et al., 2018b) and Seed protein 34-5 (Lu et al., 2013) QTLs, hence qSOC_SPC8 might be the same QTLs as that of Seed oil 8-g13 and Seed protein 34-5. The qSOC10 present on the Chr.10 was located within the physical interval (44580197- 44718071) of earlier reported QTL viz., Seed oil 39-16, hence representing the same QTL (Wang et al., 2014). Moreover, the qSPC12 fell within the genomic region (4520131- 6391062) of previously reported QTL viz., Seed protein 28-3 (Liang et al., 2010). qSOC15 located on the Chr.15 might be the same QTLs as previously reported mqSeed Oil-014 (Qi et al., 2011), as the qSOC15 fell within the genomic interval (10392903-10393032) of mqSeed Oil-014. The qSPC18 lied across the genomic region (1736692-56797602) of earlier documented QTL viz., Seed protein 36-25 (Mao et al., 2013). The qSPC20 fell within the genomic interval (2615587–41000124 bp) of previously reported QTL viz., Seed protein 36-26 (Mao et al., 2013). Zhao et al (2024) found that the SNP named rs20739 present on Chr.18 was associated with SOC, which fell within the genomic interval of qSOC18 (57712803-57856803). However, no QTL was identified to date that lies across the genomic interval of qSOC1, qSPC1, qSOC9, qSOC_SPC15.1 and qSOC_SPC15.2 (Supplementary Figure 2). Therefore, these five QTLs were newly reported/novel QTLs for SOC or/and SPC identified in the current study. Here, we report that GWAS analysis substantially diminished the genomic interval of the earlier detected QTLs, thus providing the proof for the higher resolution of GWAS analysis in QTL identification. The classical mapping method possesses reduced resolution because it utilizes the biparental population in the genetic mapping, that in turn results in a larger confidence interval of earlier documented QTLs (Yu et al., 2024). Hence, high-resolution stable QTLs identified by the GWAS in the current study will be potentially harnessed in the molecular breeding programs for the improvement of SOC and SPC in soybean.

The central objective of the researcher is to find-out true candidate genes governing the trait of interest that are underlying the major QTLs. Successful validation of function for the candidate genes underlying the QTLs will dictate their potential use in breeding programs. The true genes regulating the SOC and SPC traits have remained largely unexplored, and only a few genes regulating these traits in soybean have been characterized (as reviewed by Mo et al. (2024)). Based on the in-silico, sequence variation, variant annotation and haplotype analysis, we prioritized 80 genes as the candidates regulating the SOC and SPC in Northeast China soybeans. Among them, 12 genes undergo divergent selection between low-latitudes and high-latitudes; hence, these 12 genes will be considered as the potential candidates regulating SOC and SPC traits in soybean, as well as regulating the difference in the quality traits among the soybeans adapted to different environments. For example, the ortholog of Glyma.01G175500 in Arabidopsis thaliana is CYP90B1 (Zhang et al., 2024a). The CYP90B1 encodes the steroid C-22 hydroxylase that plays an important role in the sterol biosynthesis (Zhang et al., 2024a). The ortholog of Glyma.01G176000 is EMB2754, and the EMB2754 encodes the Vacuolar sorting protein 39, and has been reported to transport the proteins into vacuoles during seed development (Tsao et al., 2022), and lipid transfer in Arabidopsis (Michaud et al., 2017). The Glyma.08G101100 ortholog in Arabidopsis is AT1G67680 that encodes SRP72; the SRP protein mediates the transport of secretory proteins to the endoplasmic reticulum (ER) in plants, and it has been documented that storage proteins reach the protein storage vacuoles (PSV) through a vesicle-mediated biosynthetic trafficking pathway that includes the ER, the Golgi apparatus, the tubulovesicular network (TGN) and the endosomes/prevacuoles (Zeng et al., 2023). Glyma.08G101900 is an ortholog of AT5G40780 which encodes lysine histidine transporter 1 (LHT1) in Arabidopsis (Chen and Bush, 1997). The LHT1 is a high-affinity transporter for cellular amino acids in both roots and shoots, hence is important in the biosynthesis of proteins (Rabby et al., 2022). In soybean the Glyma.08G103100 is an ortholog of AT1G24330 which encodes the ARM repeat superfamily protein; and this protein belongs to E3 ubiquitin ligase. In Brassica species the ARM repeat protein viz., ARC1 acts downstream of the S receptor kinase to promote ubiquitination and protein degradation (Stone et al., 2003). The Glyma.08G111700 and Glyma.08G111800 are the ortholog of AT5G55160/SUMO2/ATSUMO2 which encodes for small ubiquitin-like modifier 2 in Arabidopsis (Table 2). The ubiquitin-like modifier stabilizes the proteins to which it is conjugated; and removal of the ubiquitin-like modifier results in degradation of the proteins (Miura et al., 2009). Gene viz., Glyma.08G117000 is an ortholog of AT5G45800/MEE62 that encodes for leucine-rich repeat protein kinase family protein, and these proteins have been documented to play important roles in plant growth, development and stress responses (Song et al., 2022). The Glyma.15G232200 is an ortholog of AT3G56710/SIB1 gene that encodes for sigma factor binding protein 1; and this gene has been documented to play an important role under nitrogen starvation, and nitrogen is the primary component in the protein synthesis (Ramel et al., 2007). The Glyma.18G258000 and Glyma.18G258100 are the orthologs of AT3G29590/AT5MAT and AT3G49190, respectively that encode for acyl-transferase family proteins, and these proteins are important in the de-novo fatty-acid biosynthesis (Yang et al., 2021). Moreover, the ortholog of Glyma.18G259100 in Arabidopsis thaliana is AT3G28910/MYB30 (Liu et al., 2014; Yan et al., 2020). The MYB30 belongs to subgroup I, that governs the biosynthesis of very-long-chain fatty acids (VLCFAs) by regulating genes encoding the four enzymes of the acyl-CoA elongase complex (Zhukov and Popov, 2022). Hence, the divergent selection of these 12 genes across the soybean populations adapted to different regions might be the reason of the unique characteristics of the SOC and SPC in Northeast China.

In conclusion, we identified the haplotype alleles within the 15 haplotype blocks as well as within the 80 candidate genes across the soybean germplasm of 340 accessions that regulate the SOC or SPC or both traits from the minimum to maximum levels through different intermediate phenotypes. Hence, these haplotypes can be used to produce cultivars of soybean with desired SOC and SPC.

Besides, we identified key genes that were divergently selected between the soybeans adapted to different regions of China, and they might be the important genes regulating the difference in the SOC and SPC across the various latitude regions.