GB13 RIL population, comprising 248 lines, constructed with Guizao1 as the female parent and B13 as the male parent was planted on three replications in July 2018, July 2019 and July 2020 at the Teaching and Research Base of Zengcheng (23°14´N, 113°38´E), and in July 2019 and July 2020 at the Experimental Teaching Base of Guangzhou (23°10´N, 113°22´E) (designated as 18ZC, 19ZC, 20ZC, 19GZ and 20GZ, respectively). In addition, one F₁ and 248 F_8:11 RILs of the GB13 RIL population were developed using the single seed descent method (Jiang et al., 2019). The experimental materials were planted in a completely randomized zonal design with one row per strain, a row length of 2.0 m, a row spacing of 0.5 m, and a plant spacing of 0.1 m. Field management was carried out according to conventional standards, and no pest or disease outbreaks occurred during growth.

Seed length (SL), seed width (SW) and seed thickness (ST) of soybean seeds were measured using an electronic vernier caliper ( Figure 1 ). Twenty seeds from each strain were randomly selected for seed size measurement, after which they were weighed on an electronic balance. After three repetitions, five times the average of the weights were used as phenotypic data for the seed weight trait (100-seed weight, HSW).

The mean of the phenotypic data from the GB13 RIL population in five natural environments was calculated as the phenotypic data of the sixth environment (combined environment, designated as CE). Descriptive statistics, linear regression analysis and calculation of the skewness and the kurtosis were performed on phenotypic data from six environments using GraphPad Prism software (GraphPad Software Inc. v.7.0.4, United States), and the corresponding tables or figures were generated. Analysis of variance (ANOVA) and estimation of broad-sense heritability (h ²) were performed on phenotypic data from the five natural environments using Genstat 12_th Edition software. The formulae (Nyquist and Baker, 1991; Elattar et al., 2021) were as follows:

where σ² _g denotes genotype variance, σ² _ge denotes genotype × environment interaction variance, and σ² _e denotes variance error. “n” denotes the number of environments, and “r” denotes the number of replicates per environment.

The composite interval mapping method (Zeng, 1994) was used to map QTLs with a high-density bin linkage map previously constructed in the laboratory, as described, the map contained 56,561 SNPs and 3715 bin markers covering a total genome length of 3049.2 cM, with the length of the linkage cluster ranging from 120.22-211.37 cM, The average distance between markers was 0.80 cM, and the number of bin markers on each cluster varied from 147 to 259. (Jiang et al., 2019). The additive effect (ADD) and PVE of each QTL were analyzed by the WinQTLCart 2.5 software. The threshold of LOD for different traits was set to 2.5 at the 5% significance level, and LOD thresholds greater than 2.5 were considered as the basis for determining the presence of QTLs. QTLs located on the same chromosome with adjacent marker intervals were combined into one QTL if their confidence intervals overlapped. QTLs were named according to “q + trait + chromosome + sequence number”, with the symbol “-” in between (Mccouch et al., 1997).

Generally, QTLs detected in at least two environments were considered stable QTLs (Qi et al., 2017), but there were six experimental environments in this study, so QTLs detected in at least three environments were considered as stable QTLs in this study to identify the stable QTLs.

All QTLs were sorted by chromosome as the primary condition and physical position as the secondary condition. QTLs with identical chromosomes in overlapping or adjacent physical position (less than 5cM) were grouped into a cluster and identified as QTL cluster if it was associated with at least two traits (Liu et al., 2017; Wang et al., 2022). To mine the main effect and stable candidate genes, the identified QTL clusters contain at least one stable QTL in this study.

First, GO (Gene Ontology) term enrichment analysis was performed for the genes within the interval of the identified QTL clusters. This is useful for selecting candidate genes based on the sequence variation of genes between the parents. The GO annotation numbers were downloaded from the soybean data website SoyBase (https://www.soybase.org/), and then the WEGO2.0 online gene enrichment analysis mapping web page (http://wego.genomics.org.cn/) was used to create the gene enrichment cellular composition, molecular function and biological processes visualization diagrams. Next, genes that were highly or specifically expressed in soybean seeds were to be further screened. Gene expression data were obtained from the soybean data website SoyBase. Finally, the selected genes were functionally annotated by the Phytozome database (https://phytozome-next.jgi.doe.gov/), and then the candidate genes were identified based on the relevant literature and related gene functions. Meanwhile, the heat map of candidate genes expression was employed to analyze tissue-specific expression (Jia et al., 2022) by the online site (https://www.omicshare.com/tools/Home/Soft/heatmap), and the structures of candidate genes were mapped using the online site of the Gene Structure Display Server (http://gsds.gao-lab.org/index.php). The resequencing data of the parental lines were referred to compare the variation of candidate genes between parental lines. The whole genome of parental lines were sequenced using the Illumina HiSeq X Ten p, with an average sequencing depth of 8× (Xian et al., 2022). High-quality sequencing data were assessed to predict the gene structural variations.

The regression linear analysis of seed size with seed weight shows that the SL, SW and ST were all highly significantly and positively correlated with HSW (p<0.01), with R-Squared (R²) of regression linear analysis ranging from 0.440 to 0.804 ( Figure 2 ). In 18ZC, 19ZC, 19GZ, 20ZC and CE, the R² between SW and HSW were the highest at 0.804, 0.721, 0.753, 0.552 and 0.742, respectively, while in 20GZ, the R² between ST and HSW was the highest at 0.754. Thus, seed size is indeed closely related to seed weight. 18ZC, 19ZC and CE had the next largest R² between SL and HSW at 0.685, 0.680 and 0.614, respectively and the smallest R² between ST and HSW at 0.679, 0.626 and 0.588, respectively. 19GZ and 20ZC had the next largest R² between ST and HSW at 0.713 and 0.467, and the smallest R² between SL and HSW at 0.664 and 0.440. 20GZ had the next largest R² between SL and HSW at 0.721, and the smallest R² between SW and HSW at 0.717. Also, the above results suggest that the correlation between seed size and seed weight seems to be applied in field selection breeding as a reference for selecting high-yielding varieties.

The phenotypic data were analyzed by descriptive statistics, and the results showed that the SL, SW, ST ( Supplementary Figure 1 ) and HSW of the female parent Guizao1 were higher than those of the male parent B13 with significant differences ( Table 1 ), which laid the foundation for QTL mapping. The GB13 RIL population differed remarkably insize and weight traits, both of which showed transgressive segregation. The coefficient of variation (CV) for SL ranged from 4.33 to 5.25%, for SW from 3.96 to 4.78%, for ST from 4.71 to 5.18%, and for HSW from 10.68 to 12.98%. The CV for HSW was greater than that for SL, SW and ST, but they were relatively stable, which indicated that both parents contained genes that were acting in the phenotypic variation.

The absolute values of the skewness and the kurtosis of SL, SW, ST and HSW of the GB13 RIL population were less than 1 ( Table 1 ). In addition, the frequency distribution graph ( Figure 3 ) showed that the phenotypic data of seed size and weight traits displayed continuous variation. The above results indicated that seed size and weight traits of the GB13 RIL population obeyed normal distribution, which was consistent with the characteristics of the RIL population and belonged to quantitative genetic traits.

The results of ANOVA for SL, SW, ST and HSW of the GB13 RIL population in five natural environments ( Table 2 ) show that the genotype, environment and interaction between genotype and environment had significant effects on seed size and weight traits of the GB13 RIL population (p<0.0001). The h ² for seed size and weight traits of the GB13 RIL population were relatively high ranging from 0.74 to 0.83 ( Table 2 ). The traits are suitable for further QTL localization analysis.

In terms of the proportion of total variation accounted for variation of different sources ( Table 2 ), the largest source of variation in SL and SW was genotype, while the largest source of variation in ST and HSW was genotype × environment interactions. The most influenced by the environment was HSW, accounted for 30.71% of the total variation. It was followed by SL, SW and ST, which accounting for 19.51%, 18.19% and 17.53% of the total variation, respectively. The above results indicated that the seed weight trait was more influenced by the environment compared to the seed size trait.

A total of 85 QTL were detected for seed size and weight in six environments with PVE in the range of 3.14% - 21.93%, of which 19 QTLs were identified for SL, 23 for SW, 18 for ST, and 25 for HSW ( Supplementary Tables 1 – 4 ). By collation and analysis, a total of 18 stable QTLs were identified ( Table 3 ), including 3 for SL, 6 for SW, 4 for ST, and 5 for HSW ( Figure 4 ). One QTL, qSL-3-1, associated with SL was detected in all environments, and it explained more than 10% of PV in three environments in the range of 11.93% - 15.91%. Three QTLs (qHSW-11-3, qHSW-20-3 and qSW-20-3) were detected in four environments. Among them, qSW-20-3 explained 9.22%, 10.79%, 20.89% and 21.93% of PV in four environments, respectively. Another 14 QTLs distributed on chromosomes 03, 04, 05, 07, 11, 12, 13, 18 and 20 were detected in three environments. In summary, integrating the PVE for QTL, the number of detections in various environments, and the relevant references, qSL-3-1 and qSW-20-3 can be considered as novel stable QTL with the main effects. Meanwhile, these stable QTLs above were important references for the excavation of genes controlling seed size and weight traits.

According to the classification criteria of QTL clusters, we obtained five QTL clusters ( Table 4 ) from the QTL mapping results ( Supplementary Tables 1 – 4 ). These five QTL clusters were distributed on chromosomes 03, 04, 13 and 20. All of them were associated with HSW, which verified that the seed size was closely related to HSW from another aspect. The four QTL clusters related to SL were qLH3-1, qLH3-2, qLWTH4 and qLH13; the two QTL clusters related to SW were qLWTP4 and qWTH20; the two QTL clusters related to ST were qLWTP4 and qWTP20. In terms of the number of traits controlled, one QTL cluster for four traits was qLWTP4, one QTL cluster for three traits was qWTH20, and three QTL clusters for two traits were qLH3-1, qLH3-2, and qLH13. Among them, qLP3-1 and qLP3-2 containing major and stable QTL(s) were located on the forearm of chromosome 03 at the physical positions between 4387438 and 6456915 bp and between 8784766 and 17199876 bp, respectively. Similarly, qWTH20 also contained major and stable QTLs and was located at the physical position of chromosome 20 between 12056958 and 33222868 bp.

In summary, qLP3-1, qLP3-2 and qWTH20 may be QTL clusters containing major effect genes that regulate soybean seed size and weight traits; qLWTP4 and qLP13 may be QTL clusters containing micro-effector genes that harmonize the control of soybean seed size and weight traits. Therefore, these five QTL clusters can provide a reference for mining the target genes controlling seed size and weight traits, and were used as the target intervals for gene GO enrichment analysis.

By Gene GO enrichment analysis, we found that most of the genes within these five QTL cluster intervals were involved in cellular processes and metabolic processes ( Figure 5 ). Most of the genes in the qLP3-1 were involved in intracellular metabolic processes, bioregulation and reaction to stimuli, and were mostly involved in binding and catalytic activities in terms of molecular functions. While most of the genes in the qLP3-2 were similar in function to those in the qLP3-1. In comparison, most of the genes in the qLWTP4 were also involved in intracellular metabolic processes, but a small number of genes were involved in regulating growth and cell proliferation. The qLP13 contained only 11 annotated genes, most of which were involved in growth and cellular processes. Similarly, most of the genes in the qWTP20-2 were involved in intracellular metabolic processes.

Through the above gene GO enrichment analysis, followed by gene expression screening, gene function annotation, and sequence variation analysis, we obtained a total of 15 candidate genes ( Table 5 ) that might be participating in the regulation of seed size and weight traits. Among them, six genes were located on chromosome 03, three on chromosome 04, two on chromosome 13, and four on chromosome 20. Among the 15 genes screened, Glyma.03g040400 was the gene with the same function as the gene (LOC_Os07g19000) associated with seed size in rice; Glyma.03g045400 and Glyma.20g084000 were related to DNA methylation; Glyma.04g100400, Glyma.20g081600 and Glyma.20g087000 were hormone-related and might be involved in the regulation of seed size and weight traits by hormone signaling pathway; Glyma.04g100100, Glyma.04g107100, Glyma.13g159500 and Glyma.20g084500 might negatively regulate seed size and weight through the ubiquitin-proteasome pathway. Glyma.03g036600 and Glyma.03g065700 were associated with transcription factor regulation and might be involved in the regulation of seed size and weight traits through the transcription factor regulation pathway; Glyma.03g044200, Glyma.03g065900 and Glyma.13g160400 might be related to the metabolic synthesis of protein and oil, and indirectly regulated the seed size and weight traits.

The expression heat map ( Figure 6A ) of these 15 candidate genes showed that two genes each were highly expressed in young leaf and nodule of soybean; three genes were specifically expressed at the pod shell development stage in soybean; eight genes including Glyma.03g036600, Glyma.03g040400, Glyma.03g045400, Glyma.03g065700, Glyma.03g065900, Glyma.04g100400, Glyma.04g107100 and Glyma.20g087000 were specifically and progressively expressed at the seed maturity stage in soybean. Meanwhile, the structural maps for most of the candidate genes ( Figure 6B ) were matched with their sequence in the parental lines, but ten candidate genes varied in the sequence of the parental lines, with the variant region on introns, 3’UTR (untranslated region) and 5’UTR ( Table 6 ). They are Glyma.03g036600, Glyma.03g040400, Glyma.03g044200, Glyma.03g045400, Glyma.03g065700, Glyma.04g100100, Glyma.04g107100, Glyma.13g160400, Glyma.20g084000 and Glyma.20g084500, and deserve to be investigated in depth.