The 1210 soybean accessions from China were collected to obtain accurate phenotypic data for the GWAS. In 2020–2022, the plant materials were planted in three regions of China (Gongzhuling, Yantai, and Guangzhou), which represent typical climates (subtropical, tropical, and temperate zones). Seeds of different soybean accessions were dried after harvest to measure the SW and fatty acid content. The NIRSTM DS 2500 (FOSS Inc., Hillerod, Denmark) was used to measure the oil content of the soybean seeds. Five FAs (linolenic acid, stearic acid, oleic acid, palmitic acid, and linoleic acid) were determined by matrix-assisted laser desorption/ionization time-of-flight imaging mass spectrometry (MALDI-TOF IMS) (Bruker Daltonics, Bremen, Germany). Each calculation was repeated thrice. The correlation coefficients were calculated using the SPSS software 20.0 (SPSS Inc., Chicago, IL, USA).

Genomic DNA was extracted from young leaves of the natural population consisting of 1210 soybean accessions using the cetyltrimethylammonium bromide (CTAB) method (Porebski et al., 1997). For each accession, 6 μg of DNA was used to construct a sequencing library according to the manufacturer’s instructions. The same genotyping method was performed as previously described (Qin et al., 2021). The restriction endonuclease combination was RsaI-HaeIII. Single nucleotide polymorphism (SNP) molecular markers were used for the GWAS analysis.

Based on the phenotypic data of 1210 soybean germplasms, the SNP markers with a minor allele frequency (MAF) > 0.06 obtained from genotyping were included in the GWAS to analyze the association between genotype and phenotype datasets. Haploview 4.2 software was used for Manhattan map construction. We used plink² software to calculate the linkage disequilibrium (LD) decay distances. The candidate genes were within the LD distance (8.9 Kb). We used 3,212,756 high-quality SNPs (MAF > 0.06) to perform GWAS for SW in the 1210 accessions. The fastlmmc model were used for association analyses. The threshold value for GWAS was set to –log(p) > 6.20. The significantly associated regions were verified from the aligned sequencing reads against the Williams 82 genome (Glycine max Wm82.a2.v1) with SAMtools ver. 0.1.18 (Li and Durbin, 2009). For each statistically significant SNP of interest, the annotated genes within 150 Kb around the SNP were considered as candidate genes (based on Williams 82).

The GmWRI14-like gene from the soybean cultivar 010a (approval number: 2012010) was cloned into the BamHI-SacI site of the plasmid pTF101 named pTF101-GmWRI14-3×Flag, which was driven by the CaMV35S promoter. The target gene was terminated by the NOS terminator (Supplementary Figure A1). The recombinant plasmid was transformed into calli using the Agrobacterium tumefaciens strain (LBA4404). Transgenic soybean plants were generated as previously described (Holsters et al., 1978). Individual T₀ plant lines were established in a greenhouse and three independent transgenic soybean lines with the same plasmid were harvested: GmWRI14-1, GmWRI14-2, and GmWRI14-3. Following a previous protocol (Bai et al., 2020), the novel CRISPR/Cas9 vector pGES201 was used for soybean genome editing. The Cas9 gene driven by the pM4 promoter (Supplementary Figure A2) and the Cas9 vector were gifted by Fujian Agriculture and Forestry University. The exon sequence of the GmWRI14-like gene in the soybean genome was determined using the CRISPR direct website (http://crispr.dbcls.jp/). We designed three target site sequences with a length of 20 bp (Supplementary Figure A2): WRI14-pM4-T1-F1/WRI14-pM4-T1-R1, WRI14-pM4-T2-F1/WRI14-pM4-T2-R1, and WRI14-pM4-T3-F1/WRI14-pM4-T3-R1. We simultaneously designed the CRISPR primers WRI14-CasF-Test/WRI14-CasR-Test. Three independent gmwri14-like soybean mutant lines containing the same plasmid were harvested: gmwri14-like-1, gmwri14-like-2, and gmwri14-like-3.

To explore the changes in differentially expressed genes (DEGs) and differential metabolic pathways in GmWRI14-like transgenic soybeans, RNA sequencing was performed on transgenic soybeans using a high-throughput Illumina HiSeq sequencing platform. Total RNA was extracted using the Eastep^® Super total RNA extraction kit (TaKaRa, USA). RNA quality was assessed using the Nanodrop 2000c (Thermo Scientific, Hudson, NH, US). RNA-seq library preparation and sequencing were performed as previously described (Wang et al., 2014). The 2^-ΔΔCt method was used to assess the fold change in gene expression (Livak and Schmittgen, 2001). To analyze the transduction pathways induced by the GmWRI14-like gene, the total RNA of each sample (GmWRI14-like transgenic soybeans, CK, gmwri14-like mutants) was isolated from the seeds for the RNA-sequencing; 2 g of RNA was obtained from each plant sample. The RNA-sequencing depth was 10 ×. The bioinformatics software (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) was used to improve data quality. In total, we obtained 41.66 Gb of clean data in this project. The statistical power of this experimental design, calculated in cluster Profiler Power, was equal to a q value ≤ 0.01. Gene Ontology (http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genome (KEGG) statistical analysis (http://www.genome.jp/kegg/) were performed to assess the function of the DEGs.

After 10 days of podding, RNA was extracted from the seeds of the soybean lines (GmWRI14-like transgenic soybeans, gmwri14-like mutants, and control soybean lines). We then conducted 24 h sampling daily at 10, 11, 12, 13, and 14 days after podding (GT0, GT4, GT8, GT12, GT16, GT20, and GT24 soybean seeds, respectively) to extract total RNA. qRT-PCR was used to compare the overexpression of GmWRI14-like soybeans on different days and at different times to explore the molecular regulatory network related to the GmWRI14-like, gmwri14-like mutants, and control soybean plants involving the DEGs. cDNA synthesis was performed using a reverse transcription kit EZBiosciences (Omega. USA). Real-time fluorescence analysis was performed using a LightCycler (Roche, Rotkreuz, Switzerland) according to the manufacturer’s instructions (Livak and Schmittgen, 2001). All primers used are listed in the additional files (Supplementary Table A1).

We performed a Y2H assay by screening the cDNA library constructed from the seed aged 15 days. The full-length GmWRI14-like protein was used as a bait, and we identified its putative partner using Y2H library screening. Y2H assays were performed as described previously (Xiaofei et al., 2017). The BiFC assay was performed as previously described (Mou et al., 2021).

Over 3 years of measurements, the variation in SW ranged from 5.43 g to 27.11 g in 2020 (standard deviation (SD): 7.6), 17.88 g to 26.89 g in 2021 (SD: 3.7), and 2.21 g to 23.59 g in 2022 (SD: 6.7). The variation in oil content ranged from 3.35% to 25.12% in 2020 (SD: 7.9), 15.52% to 25.19% in 2021 (SD: 2.4), and 17.33% to 27.83% in 2022 (SD: 2.6). The distribution of the SW was continuous, and both the SW and oil content conformed to normal distributions (Figures 2A, B).

GWAS using the fastlmmc model was used to identify potential candidate genes associated with SW. The high-quality 3,212,756 SNP markers obtained by SLAF-Seq technology were used for GWAS for SW in the 1210 accessions. The SNP marker distributions are shown in Figure 2C. Manhattan maps were used to show the oil content (Figure 2D) and SW (Figure 2E) across different years; the threshold value was set at –log(p) > 6.20 (red). We identified one significant association loci for SW on chromosome 4 (Figure 2E), and we estimated the candidate region to be from 1.600 to 1.750 Mb (150 Kb) around the leading SNP in chromosome 4. This region contains 14 genes (Supplementary Table A2). According to Swissport annotation, these genes were annotated, such as gibberellin 20 oxidase, phosphatidylethanolamine-binding protein, and dof zinc finger protein, and had already been reported to affect SW (Wang et al., 2007; Lu et al., 2016a); therefore, we focused on Glyma.04G116500.1, which we labeled GmWRI14-like. It is speculated that GmWRI14-like belongs to the WRI1 family with two AP2/EREB domains (Supplementary Figure A3). The soybean line with the lowest SW (KH1088) harbored one SNP predicted to cause a gain in a stop codon, resulting in the premature termination of translation after 241 amino acids in the 710-amino acid GmWRI14-like protein.

Using a stable soybean genetic transformation system, transgenic soybean plants overexpressing GmWRI14-like were successfully obtained (Figure 3A). Overexpression of GmWRI14-like in soybeans significantly increased seed size and SW compared to the wild-type soyabeans (P<0.01), resulting in severe bending of soybean branches (Figure 3A). The gmwri14-like soybean mutant was obtained using CRISPR/Cas9 genome editing technology (Figure 3A). Compared with the non-transgenic soybean control, the seed size of the gmwri14-like soybean mutant was significantly reduced (P<0.01), and the seed epidermis of the gmwri14-like soybean mutant exhibited wrinkling (Figure 3B). Overexpression of GmWRI14-like resulted in a 15.32% increase in seed length (P<0.01), whereas the seed length of the gmwri14-like mutant was reduced by 23% (Figure 4F). Furthermore, the width of the GmWRI14-like seeds was 14.29% greater than that of the CK. GmWRI14-like-3 had the largest width, reaching 7.14 ± 0.03 mm, whereas the gmwri14-like mutant seeds showed a 18.23–23.57% decrease in width. The width of GmWRI14-like-3 reached 6.34 ± 0.16 mm (Figure 4G). The phenotype of GmWRI14-like transgenic soybeans showed a significant difference in the color of the seed coat from the different soybean lines (P<0.01). The non-transgenic soybean seed coat was yellow, whereas the seed coat color of GmWRI14-like soybeans tended to be white. These results suggested that the GmWRI14-like gene regulates the cytochrome gene expression related to the seed coat.

The oil content and distribution in gmwri14-like soybean mutants and overexpressing GmWRI14-like soybean seeds were detected by using MALDI-TOF IMS (Figure 4A). The GmWRI14-like overexpressing soybeans seeds showed a 16.14–19.92% increase in 100 seed weights (Figure 4B). The MALDI-TOF IMS results showed that, compared to control, the relative content of oleic acid in the gmwri1-like mutant significantly decreased from 22.12% to 15.14% (P<0.01), whereas the oleic acid content in the GmWRI14-like overexpressing soybeans increased from 22.12% to 74.31% (Figure 4C). The relative content of linoleic acid in the gmwri1-like mutant increased by 8.2% (P<0.01), whereas the linoleic acid content in GmWRI14-like soybeans decreased from 49.31% to 18.42% (Figure 4D). There was no significant difference in protein content among the GmWRI14-like soybean, gmwri14-like mutant, and non-transgenic soybean (Figure 4E). These results indicate that the GmWRI14-like gene does not affect the protein content in soybean seeds.

RNA-seq was performed on the transgenic and control groups (three replicates), resulting in more than 5.1 × 10⁷ clean reads; 2101 DEGs were screened through transcriptome analysis. The differential expression gene analysis (DEGVP) showed that two genes related to cytochrome metabolism were differentially expressed (more than 2.5 times), namely, GmCYP78A50 (Glyma.07G052300) and GmCYP78A69 (Glyma.16G021200) (Figure 5A). We analyzed the functions of the target genes and potential target genes using the KEGG pathway (Figure 5B), which showed that the DEGs were enriched in three signaling pathways related to seed development, including fatty acid biosynthesis (50.83%), fatty acid metabolism (36.37%), and cytochrome biosynthesis (19.17%) (Figure 5B).

To explore the role of GmWRI14-like in more detail, we measured the expression level of the GmWRI14-like gene in different soybean tissues (root, stem, leaf, and seed) using qRT-PCR and found that the GmWRI14-like gene had the highest expression level in soybean seeds (Figure 5C). The qRT-PCR results showed that GmWRI14-like gene expression was highly tissue-specific in soybean plants (Figure 5C). We used the lectin gene (Gene ID: 100775957) as a reference gene and detected the expression levels of GmCYP78A50 and GmCYP78A69 in different transgenic soybean tissues using qRT-PCR, which suggested that the expression levels of GmCYP78A50 and GmCYP78A69 were the highest in the epidermis of transgenic soybean seeds (Figure 5D). GmCYP78A50 and GmCYP78A69 showed an increased expression level in the transgenic plants than in the control plants (Figure 5E). Based on the RNA-seq and qRT-PCR results, we speculate that the GmWRI14-like gene regulates the expression of GmCYP78A50 and GmCYP78A69 in soybean seeds.

We found 140 unique cDNA clones using Y2H. Interestingly, among the interacting proteins, GmCYP78A50 (Glyma.07G052300) and GmCYP78A69 (Glyma.16G021200) were identified multiple times. We further confirmed the interaction between GmWRI14-like and GmCYP78A50/GmCYP78A69 by point-to-point Y2H. First, we used GmWRI14-like bait and GmCYP78A50/GmCYP78A69 as prey (i.e., fused with GAL4-BD and GAL4-AD, respectively) and confirmed the interaction through Y2H. As an additional control, we fused another coiled coil domain protein (structural maintenance of chromosome 5) with GAL4-AD and another AP2 transcription factor (dehydration response element binding protein 2A) with GAL4-8D and tested their interactions with WRI1-BD and GmCYP78A50/GmCYP78A69-AD, respectively. Only when GmWRI14-like-BD binds to GmCYP78A69/GmCYP78A50-AD clones can they grow on selective media, indicating that the interaction is specific (Figures 6A, B). In addition, the interaction between GmWRI14-like and GmCYP78A50/GmCYP78A69 was further confirmed by BiFC (Figure 6C). Therefore, these results show that the GmWRI14-like protein can interact with CYP78A50/GmCYP78A69.