To find the KIX domain in 59 species ( Figure 1 ), KIXBASE (http://www.nipgr.res.in/kixbase/home.php) database was used (Yadav et al., 2017). To identify the KIX domain in soybean (Glycine max Wm82.a2.v1 genome version), we used the sequences of Arabidopsis KIX domains as query sequences for performing BLASTP searches on the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/), Phytozome website (https://phytozome.jgi.doe.gov/pz/portal.html), and SoyBase databases (https://soybase.org). To confirm the KIX domain in the selected GmKIX proteins, we used the KIX_prediction tool in KIXBASE and the Searching Protein Sequence Motifs (MOTIF, https://www.genome.jp/tools/motif/) for validation.

Phylogenetic analysis was performed to identify the evolutionary and functional relationship among the species. To better understand the evolutionary relationship among the KIX domains in the genome of various plants, we constructed a phylogenetic tree using the amino acid full-length sequence of 591 KIX domains from 59 species representing major plant groups ( Figure 1A ). The phylogenetic tree was constructed using molecular genetics analysis (Wickett et al., 2014) X with the neighbor-joining (NJ) method and bootstrap analysis was conducted using 1,000 replicates. The consensus tree and unrooted tree were redrawn using Figtree (http://tree.bio.ed.ac.uk/software/figtree/).

To predict the function and investigate the structural characteristics of KIX domains, we collected protein and nucleotide sequences of 48 KIX domains from Arabidopsis (11), soybean (13), Cicer arietinum (5), Medicago truncatula (11), and Phaseolus vulgar (8) through Phytozome. The sequence information was obtained from Phytozome for Arabidopsis thaliana TAIR10, Cicer arietinum v1.0, Medicago truncatula Mt4.0v1, and Phaseolus vulgaris v2.1.The exon/intron structures of the collected KIX genes were visualized using Gene Structure Display Server (GSDS, http://gsds.gao-lab.org/) (Guo et al., 2007). The functional domains of all protein sequences encoded by the candidate KIX genes were predicted using Simple Modular Architecture Research Tool (SMART, http://smart.embl-heidelberg.de/) and Searching Protein Sequence Motifs (MOTIF, https://www.genome.jp/tools/motif/). To confirm the consistency of the KIX domains between Arabidopsis and soybean, Weblogo (https://weblogo.berkeley.edu/logo.cgi) (Crooks et al., 2004) and MEGA-X (Kumar et al., 2018a) were used.

Haplotype analyses of GmKIX genes were performed using whole-genome re-sequencing data from a soybean core collection of 422 accessions(Kim et al., 2021) ( Supplementary Table 8 ). The whole-genome re-sequencing data were utilized using the SRA accession: PRJNA555366, which has been made publicly available by Kim et al., 2021. To filter the re-sequencing data, monomorphic and low-coverage site SNP markers were removed, and those with a minor allele frequency (MAF) less than 0.05 were excluded to minimize the potential influence of rare alleles on the analysis. Additionally, SNPs with missing data for more than 10% of the accessions were removed to reduce the impact of incomplete genotyping information. These filtering steps and genetic admixture analysis were performed using the QTLmax 3.0 program (https://www.qtlmax.com). The soybean core collection was cultivated in an experimental field at the National Institute of Crop Science in 2017 and 2018. After the soybean seeds were harvested and naturally dried to achieve a stable seed weight, phenotypic measurements were conducted indoors. The measured seed agronomic traits were 100-seed weight (100-SW), area, thickness, and major and minor axes ( Supplementary Table 8 ).

We performed filtering of the soybean core collection re-sequencing data using QTLmaxV3.0. The filtering process included a MAF threshold of < 5%, a limit of 0.1% for missing SNPs, and a stringent threshold for Hardy-Weinberg equilibrium (P-value < 10e-6). Using the filtered set of high-quality SNPs (542,422), we conducted a population structure analysis on the core soybean group, incrementally increasing the K value from 2 to 5, in order to identify an appropriate cluster.

To examine the soybean core collection distribution classified by population structure in haplotypes, we generated a haplotype network using PopART v1.7 (Leigh and Bryant, 2015).

Four cultivars of soybean, namely Hoseo, PI86490, KLS88035, and Soheung-2, were grown in greenhouses to analyze the expression of the GmKIX gene during seed development. To analyze the RNA expression levels during seed development, the process was divided into three stages. Stage 1 (S1) included small-seeded cultivars Hoseo and PI86490 with seed sizes less than 3 mm and large-seeded cultivars KLS88035 and Soheung2 with sizes less than 5 mm. Stage 2 (S2) included small-seeded cultivars with seed sizes ranging from 3 to 6 mm and large-seeded cultivars ranging from 5 to 10 mm. Stage 3 (S3) included small-seeded cultivars with sizes greater than 6 mm and large-seeded cultivars with sizes greater than 10 mm ( Supplementary Figure 2 ). For each developmental stage, a sufficient number of seeds were promptly collected, immediately frozen in liquid nitrogen, and stored at -80°C for subsequent analysis. Total RNA was extracted from organoid cells using the RNeasy^® Plant Mini Kit (Qiagen), and RNA-seq libraries were prepared according to the manufacturer’s protocol. Paired-end RNA-seq reads were generated on the Illumina Genome Analyzer platform, and the quality of the trimmed reads was assessed using FastQC. The expression levels were determined by calculating the reads per kilo-base of the exon per million mapped reads (RPKM). The gene expression profiles were visualized using Pheatmap software (Kolde, 2012). The RNA-seq data reported in this article has been deposited in NCBI under SRA accession: PRJNA1003551.

Statistical analyses were conducted using R package. Mean differences among the genotypic groups were analyzed using Fisher’s least significant difference test at a p value of 0.05 using PROC GLM.

To understand the distribution and evolutionary relationship of the KIX domain in various plant species, we extracted and analyzed amino acid full-length sequence of 591 KIX domains from 59 plant species, including primitive plants, non-seed plants, monocotyledons, and dicotyledons ( Figure 1A and Supplementary Table 1 ). In primitive plant groups classified as Chlorophyta and Streptophyta, a small number of (1 to 3) KIX domains were identified. Selaginella moellendorffii, a more evolved species belonging to Tracheophyta, is considered to have a closer relationship with higher plants among non-seed plants (Banks et al., 2011). In Selaginella moellendorffii, four KIX domains were found, whereas eight were detected in Amborella trichopoda, currently known as the most basal angiosperm (Islam et al., 2012). With an average of 8.6 KIX domains identified in monocots and approximately 14 in dicots, it was observed that KIX domain-containing proteins were conserved and increased in number throughout the evolutionary process. In monocotyledonous plants, the highest number of KIX domains were identified in Oryza rufipogon and Oryza sativa. Following them, Oryza meridionalis and Oryza nivara, closely related to Oryza sativa’s ancestors, also showed 11 KIX domains each. Similarly, barley and wheat, two major staple crops, encoded 11 KIX domains each. The abundance of KIX domains identified in plants primarily used as human food or subjected to domestication is intriguing. In dicotyledonous plants, particularly in the genus Brassica, a significant number of KIX domains were identified. Brassica oleracea and Brassica rapa, diploid species, displayed 24 and 26 KIX domains, respectively. In the tetraploid species Brassica napus, approximately twice the number of KIX domains, around 47, were identified. In tetraploid Brassica napus, there was a significant increase in the number of KIX domains, whereas in the hexaploid monocot Triticum aestivum, a dicotyledonous plant, only 11 KIX domains were identified, making it difficult to determine the variation of KIX domains based on ploidy. Furthermore, no specific trend has been observed concerning chromosome numbers. The interpretation of these results should consider the scope and accuracy of chromosome deciphering in plants while also being mindful of potential biases.

A phylogenetic tree was constructed using the full-length sequences of 591 KIX domains identified from 59 different plant species to investigate their evolutionary relationships and distribution. The phylogenetic analysis classified the proteins into four main groups ( Figure 1B and Supplementary Table 1 ). Previously identified KIX domains, including KIX8/9, HAC, MED15, and RECQL5, were distinctly divided into these four clusters. The groups were designated as Group-I, Group-II, Group-III, and Group-IV, containing KIX8/9, HAC, MED15, and RECQL5 proteins, respectively. Group-IV, which included RECQL5, also contains WPP proteins and uncharacterized proteins. Therefore, Group-IV was further divided into Group-IV-A, which included RECQL and tryptophan-proline-proline (WPP), and Group-IV-B and Group-IV-C, which include uncharacterized proteins ( Supplementary Table 1 ). Group-I contained AtKIX8/9, which interestingly showed even distribution in diverse plants but was rarely detected in mosses and plants in the Poaceae family ( Figure 1B ). Our phylogenetic analysis suggests that, similar to previous studies (Thakur et al., 2013), KIX8/9 proteins belonging to Group-I have evolved from a common ancestral sequence with the KIX domain of HAC proteins classified into Group-II. Interestingly, HAC proteins were detected in all of the 59 plant species except seven. Group-III was classified based on the MED15-related proteins, which means that the proteins belonging to the AtMED15 family were assigned to this group. In the case of the MED15 family, a limited distribution was observed in green algae, while a relatively even distribution was found in monocots and dicots. Approximately 40% of the 591 KIX domains belonged to the MED15 family. This highlights the significance of the MED15 family, as they make up a significant proportion of the KIX domain. Group-IV, a distinct cluster, primarily consisted of RECQL, tryptophan-proline-proline (WPP), and uncharacterized proteins. An interesting aspect of Group-IV was that it primarily comprised plant proteins with little to no research or functional predictions. Group-IV could be further subdivided into Group-IV-A, consisting of RECQL and WPP proteins, Group-IV-B, consisting of uncharacterized proteins closely clustered with MED15, and Group-IV-C, consisting of a total of 52 uncharacterized proteins primarily found in Poaceae plants.

Using publicly available complete genome sequences of soybean and Arabidopsis, we identified all possible genes that encode the KIX domain. Ultimately, we identified 13 orthologs corresponding to 11 KIX domains of Arabidopsis ( Figure 2 and Table 1 ). KIX domains from soybean were found to contain one or two KIX domains. To increase the reliability of our results and analyze the trends, we conducted a comprehensive analysis by including 24 KIX domains identified in Phaseolus vulgaris, Cicer arietinum, and Medicago truncatula, in addition to soybean and Arabidopsis. This analysis revealed the presence of four conserved groups ( Figure 2 ). Importantly, these four groups corresponded to the divisions made using 591 KIX domain sequences obtained from 59 different plant species, further supporting our findings. Group-I contained AtKIX8 and AtKIX9, which show high similarity to GmKIX8-1, GmKIX8-2, GmKIX8-3, GmKIX9-1, and GmKIX9-2. The differentiation between GmKIX8 and GmKIX9 was based on the similarity to AtKIX8 and AtKIX9, respectively ( Supplementary Table 2 ). Group-II contained AtKIX3, AtKIX4, AtKIX5, and AtKIX7, reported as p300/CBP related gene, AtHAC12, AtHAC1, and AtHAC5 respectively. The two GmKIX genes clustered with HAC were named GmHAC-1 and GmHAC-2, respectively. Group-III contained AtKIX1, AtKIX2, and AtKIX6, reported as AtMED15-1, MED15-like protein, and AtMED15-2, respectively. In soybean, four highly similar genes were identified and named GmMED15-1, GmMED15-2, GmMED15-3, and GmMED15-4. In Group-IV, AtKIX10 and AtKIX11 were classified and reported as AtRECQL3 and AtWPP1, respectively. In soybean, one similar protein each for AtRECQL3 and AtWPP1 was selected and named GmRECQL and GmWPP, respectively. In three different legume plants, a set of genes corresponding to each group were found to be distributed, with the presence of KIX domains.

To understand the structural and functional diversity of the KIX domains, a comparative analysis of protein architecture was performed ( Figure 2B ). Additionally, 24 KIX domains from Fabaceae species, including Phaseolus vulgaris, Cicer arietinum, and Medicago truncatula, were included in the study [62]. Using the KIX domains of Arabidopsis as a reference, other similar proteins were clustered into four groups ( Figure 2B ). We observed a remarkable similarity in the KIX domain structures within each group, including the presence and location of the KIX domains, as well as the protein size and arrangement of other domains. The KIX8/9 and MED15 proteins were clustered in Group-I and Group-III, respectively. Apart from two CHASE (Cyclases/Histidine kinases Associated Sensory Extracellular) domains identified in Medtr2g104430, no other known domains beyond the KIX domains were detected in these two groups. However, Group-II consisted of HAC proteins, which had conserved domains, including ZF-TAZ, PHD, Hat_Kat11, and ZZ, with the KIX domain. A notable feature in Group-IV is the conservation of RECQL proteins. Distinct from other KIX domains, RECQL proteins preserve the KIX domain at the C-terminal. Additionally, they contain helicase and DEAD domains. The WPP and uncharacterized proteins in Group-IV exhibit diverse sizes and structures, indicating that they are not uniform.

In order to determine the numbers and positions of exon/intron within each GmKIX gene, we compared the full-length gDNA sequences with the corresponding Arabidopsis KIX gene sequences ( Supplementary Figure 3 ). KIX genes possess multiple exons and introns, yet their structure and arrangement were largely conserved within the group, typically exhibiting similar patterns. The length of each exon is described in detail in Supplementary Tables 3, 4 . The genes of Group-I (KIX8/9) had four exons and each exon had a similar length. The KIX genes of Group-II (HAC) had the most exons with 16 to 18. Group-III (MED15) showed 11 to 12 exons Furthermore, it was confirmed that RECQL, which was included in Group-IV, preserved 19 exons. The observed exon lengths of GmKIX genes were very similar to those reported in Arabidopsis and Rice (Thakur et al., 2013). Additionally, we confirmed that the structure of KIX genes in four other Fabaceae species was also similar to that of Arabidopsis, indicating their conservation. On the other hand, it was observed that the UTRs and introns of KIX genes exhibited significant variations in length. Specifically, GmMED15-2 possessed a third intron of approximately 13 kb, whereas the remaining genes within the same group had introns with sizes up to 3 kb. The third intron of KIX8/9 also exhibited significant variation, with the third intron of AtKIX being 79 bp in length, while GmKIX9-2 had a much larger intron size of 2,301 bp. In summary, there was a tendency of exon conservation among KIX genes based on their respective groups, while introns exhibited significant variations.

Through gene structure and motif analysis, we had previously identified differences among groups. Further, we sought to understand functional diversity by analyzing the variations and patterns in the amino acid sequence of the primary KIX domains. The KIX domain is characterized by a conserved structural fold consisting of three helix bundles that mediate the interaction with binding proteins. Hydrophobic interactions between helices contribute to the formation of a robust fold in the domain and aid in stabilizing the binding with interacting partners (Radhakrishnan et al., 1997; Brzovic et al., 2011). Despite the conservation of this fold, KIX domain sequences exhibit significant diversity, contributing to their functional flexibility (Yadav et al., 2017). To confirm the preservation of the 3-helix structure in the selected KIX domains from soybean and the diversity of KIX domain sequences, we analyzed the KIX domain sequences of soybean and Arabidopsis ( Supplementary Table 5 , and Supplementary Figure 2, 4 ). The KIX domain of the selected GmKIX genes also maintained the 3-helix bundle structure, and the amino acid residues critical for structural stability were more highly conserved than other residues. Additionally, we investigated the conservation of domain sequences within each group of GmKIX proteins. As a result, among the four groups, Group-I, which included KIX8/9, was found to have the highest sequence conservation in the KIX domains. Query coverage scores were above 97%, and the positive scores were approximately 90%, indicating a match with Arabidopsis. Groups containing proteins such as MED15, HATs, and RECQL exhibited positive scores of around 55-62% for the KIX domains, suggesting their potential to contribute to functional flexibility. An interesting observation is the rarity of fully conserved amino acid sequences in the KIX domain sequences between soybean and Arabidopsis. In the domain sequence alignment, only the 22^nd and 49^th amino acids were perfectly conserved as Arg and Glu, respectively ( Supplementary Figure 4 ). Other sequences appeared to have diverged and undergone variations based on their respective genes and groups.

We reconstructed the evolutionary tree using only KIX domain sequences diversified according to their functionalities ( Supplementary Figure 4 ). The tree constructed using only KIX domain sequences exhibited remarkably high similarity to that constructed using full amino acid sequences. This result suggests that the sequence variations in the KIX domain have occurred in conjunction with the functional diversification of KIX domain-containing proteins. It is anticipated that the function of proteins containing the KIX domain can be predicted solely based on the sequences. We identified six key amino acid sequences that contribute to the classification of KIX domains into four groups. The selection of these six amino acid sequences was primarily based on their differential conservation across the four groups of KIX domains. In particular, the 66^th amino acid exhibited distinct characteristics in each group: Glu, a polar and negatively charged amino acid, in Group-I; Lys, a polar and positively charged amino acid, in Group-II; Gln, a polar and uncharged amino acid, in Group-III; and Gly, an amino acid classified under special cases group, in Group-IV. Hence, the selected six amino acids can serve as benchmarks for differentiating the functions of KIX domains through one or multiple combinations. Therefore, these six amino acids have undergone increased diversification during evolution, allowing for functional diversification of the KIX domain in polyploid plants, including soybean.

Gene structure and domain sequence analyses revealed that GmKIX genes within the same group exhibit conserved structures and domain sequences. Furthermore, we focused on variations in exon sequences to conduct a detailed analysis of sequence variations within the conserved coding region. In order to investigate the variations and diversity within GmKIX genes that have undergone functional differentiation, we performed haplotype analysis using re-sequencing data from the soybean core collection consisting of 422 accessions.

We removed heterozygous variants and sequencing errors from the filtered mutations in the re-sequencing data of 422 soybean varieties. We focused on haplotypes that exhibited non-synonymous substitutions and functional InDels among the filtered mutations. The analysis revealed an average of 3 haplotypes in Group-I, 8 in Group-II, 7 in Group-III, and 2 in Group-IV ( Figure 3 ). The haplotype distribution of each GmKIX gene was predominantly characterized by a single dominant haplotype, except for GmMED15-1, GmHAC-1, and GmRECQL, which exhibited two or more equally distributed haplotypes. Interestingly, GmKIX genes belonging to the same group, namely GmKIX8-3, GmKIX9-1, and GmKIX9-2, have approximately 900 bp exons, and among them, no significant variations were observed in the coding region. On the other hand, GmKIX8-1 and GmKIX8-2 harbored a 1,200 bp exon and exhibited six distinct haplotypes. These findings suggest that in cases where variations are limited, such as in GmKIX8-3, GmKIX9-1, and GmKIX9-2, sequence changes in the coding region may potentially have significant implications in the plant system.

We performed the haplotype-based association analysis of the 13 GmKIX genes related to seed-related agronomic traits using 422 soybean accessions. As a result, we confirmed the correlation between the genetic variations of GmKIX8-1 and GmMED15-4 and the phenotypes in seed-related agronomic traits ( Figures 4 , 5 ). GmKIX8-1, with four main haplotypes, had a sufficient number of resources for statistical analysis ( Figure 4A ). Variations in the haplotypes of GmKIX8-1 were observed only in the fourth exon, where a total of 8 non-synonymous mutations were present. Upon examination of seed-related agronomic traits according to each haplotype, Haplotypes Hap-1 and Hap-2 exhibited a relatively larger and heavier seed shape, while Hap-3 and Hap-4 showed a relatively smaller distribution in seed-related agronomic traits. The distribution of 100-SW revealed that the mean values for Hap-1, Hap-2, Hap-3, and Hap-4 were 24.6 ± 6.8 g, 27.3 ± 7.5 g, 15.5 ± 4.9 g, and 19.8 ± 6.9 g, respectively ( Figures 4B , 6C ). While there was no significant difference between Hap-1 and Hap-2, there was a significant difference between these two haplotypes and the other haplotype groups. The SNP that distinguishes Hap-1 & 2 from Hap-3 & 4 is the G to A change at position 8,909,112, resulting in the conversion of Gly²¹⁹ to Asp²¹⁹, which enables differentiation between the two haplotypes. Based on these results, it can be understood that the one-base substitution at position 8,909,112 may impact the function of GmKIX8-1 and thus affect the development related to seed-related agronomic traits. SNPs of GmMED15-4 were distinguished as two haplotypes. Specifically, a variation in the nucleotide sequence at position 2,201,699 of the KIX domain of GmMED15-4 was observed with A and C alleles, which encoded Thr⁷⁵ and Pro⁷⁵, respectively ( Figure 5A ). The correlation between haplotypes containing the variation in the amino acid sequence at position 75 and seed agronomic traits showed mostly small tendencies in Hap-1 compared to Hap-2. The average 100-SW ( Figure 5B ) for Hap-1 and Hap-2 were 24.48 ± 8.3 g and 16.3 ± 7.9 g, respectively, with a p-value <0.01 indicating a significant difference between the average 100-SW for Hap-1 and Hap-2. The two genes, GmKIX8-1 and GmMED15-4, showed differences not only in 100-SW but also in area, thickness, and minor axis, while no correlation was found between the major axis and haplotype. Both GmKIX8-1 and GmMED15-4 have been reported to have an impact on plant size and seed size through genome editing using CRISPR-Cas9 in soybean, as well as significant differences in SNP and seed development and morphology in rice (Thakur et al., 2013; Nguyen et al., 2021). These findings further enhance the credibility of our association analysis between haplotypes and seed-related agronomic traits. No significant difference was observed between the haplotypes of other GmKIX genes and seed-related agronomic traits in this study.

We further conducted population structure analysis and investigated the distribution of seed-related traits based on the genetic data of the core soybean collection in Korea. After filtering, a total of 542,422 high-quality SNPs were obtained from the re-sequencing data of the soybean core collection. Population structure analysis was performed using the high-quality SNPs, and based on the reference results of the 180K SNP analysis of the core soybean collection by (Lee et al., 2022), the optimal value for K was determined to be 4 ( Figure 6A ). The resulting clusters were labeled as Subpopulation (SP)-I, SP -II, SP-III, and SP-IV, comprising 105, 158, 101, and 51 resources, respectively. To investigate the association between the four SP and 100-SW, we initially performed linear regression analysis with 100-SW as the dependent variable. In this analysis, we used q-values obtained from the four subpopulations as covariates. Due to the high correlation among the q-values, which raised concerns of multicollinearity, we systematically omitted one q-value at a time and conducted the regression analysis with the remaining three q-values ( Supplementary Table 6 ). Consistently, the majority of q-value covariates displayed a statistically significant relationship with 100-SW values (p < 0.01). The R² value was determined as 0.46. Furthermore, in order to assess the association between the candidate gene’s haplotype and 100-SW, we conducted a linear regression analysis with population structure effects as covariates, revealing a statistically significant relationship with p < 0.01 ( Supplementary Table 7 ). Upon investigating the distribution of 100-SW among the resources within each cluster, it was observed that SP-II exhibited significantly higher 100-SW values than the resources in the other clusters ( Figure 6B ). Furthermore, significant differences were observed among the clusters regarding the admixture results. We further investigated the relationship between subpopulations within the haplotype distribution of the GmKIX8-1 and GmMED15-4 genes and their association with 100-SW ( Figures 6C, D ). In the case of GmKIX8-1, among its four haplotypes, Hap-1 and Hap-2 were characterized by relatively larger seed sizes, and these haplotypes predominantly constituted the resources of SP-II, which had the largest average seed sizes. Subsequently, the resources within SP-I were found to be distributed next. In contrast, Hap-3 and Hap-4, characterized by relatively smaller seed sizes, exhibited a higher distribution of SNPs in resources within SP-III and SP-IV, where smaller resources were predominant. Regarding the haplotypes of GmMED15-4, Hap-1, which includes relatively larger seeds, the highest number of SNPs in resources was seen within SP-II. On the other hand, Hap-2, characterized by relatively smaller sizes, had the highest proportion of resources classified under SP-IV, consisting of smaller resources. The results of the haplotype network and population structure analysis strengthen the confidence in the association analysis between GmKIX gene haplotypes and seed-related agronomic traits using re-sequencing data from the core soybean collection.

Transcriptome sequencing (RNA-seq) data from three different developmental stages of four soybean accessions were used. RPKM values were standardized as Z-scores to compare the expression according to the seed development stage for each GmKIX gene. The expression of GmKIX genes was downregulated according to the seed development process and showed differences between small (Hoseo and PI86490) and big seeds (KLS88035 and Soheung-2) ( Figure 7 ). Initial seed Stage (S1) had a high expression of GmKIX genes in both small and big seed breeds compared to the other stages. There was a difference in the expression of the GmKIX gene between the small seed and the big seed in the expansion stage (S2). Most of the GmKIX genes still showed relatively high expression in small seeds, while the amount of expression in the big seeds was significantly reduced compared to S1. In the filling stage (Sterner and Berger, 2000), it can be seen that the amount of expression decreased in both small and big seed varieties. Interestingly, it was confirmed that the small and big seeds showed similar expression patterns in S1 and S2, while the expression was maintained in the two accessions with small seeds in the vigorous seed development stage, such as S2. In addition, normalized RPKMs (log₂ scale) were checked to detect the number of genes actively expressed in the GmKIX gene ( Supplementary Figure 5 ). The GmKIX9-1 and GmKIX9-2 genes showed lower expression than other GmKIX genes, and the GmMED15 and GmHAC genes showed relatively higher expression. Upon investigating the expression patterns during the early stages of seed development, the expansion phase, and the seed filling phase, we concluded that the expression of these GmKIX genes is involved in the initial stages of seed development and may influence size determination.