We used V. stipulacea (accession JP252948 in NARO Genebank) and its ethyl methanesulfonate-induced isi2 mutant (Takahashi et al., 2019). Plants were cultivated in 5-L pots with seven biological replicates in a single greenhouse maintained at 25°C or higher. Plants were scored for days to flowering, height 46 days after sowing, largest terminal leaflet width and length, shoot dry weight, seeds per pod, pod length, 100-seed weight, and water uptake in 20 seeds per individual following storage at 4°C for 1, 6, and 24 months after harvest. We counted the number of seeds taking up water at 1, 2, 3, 7, 10, 14, 21, and 28 days after watering. Data for one-month-old seeds and six-month-old seeds were published in our previous study (Takahashi et al., 2019).

We used the SPAD chlorophyll meter (SPAD-502, KONICA MINOLTA, INC., Tokyo, Japan) to measure early senescence. Fully expanded third terminal leaflets were detached from plants cultivated for one month in 5-L pots with 5 biological replicates in a single greenhouse maintained at 30°C or lower, then incubated on the water surface in 9-cm Petri dishes for 28 days under a cycle of 12 h of light at 25°C and 12 h of darkness at 22°C.

Water entry points were identified using the method of Soltani et al. (2021). Seeds were submerged in Lugol solution containing 33 g/L of I₂ (092-05422, FUJIFILM Wako Pure Chemical, Tokyo, Japan) and 66 g/L of KI (160-03952, FUJIFILM Wako Pure Chemical) as a contrasting agent for 6 h at 25°C, followed by rinsing with water and drying with a paper towel. Seeds were imaged using an inspeXio SMX-225CT FPD HR (Shimadzu Corporation, Kyoto, Japan). Each scan included 1200 projections, averaging two frames over 360° (pixel detector resolution 3000 × 3000), at 4 frames per second. Finally, 1024 × 1024-pixel slices were computed at a final spatial resolution of 15 μm. Tube voltages of 75 kV and tube currents of 100 μA were used. Metal filters were not used during imaging. Reconstructed CT volumes were visualized using VG Studio MAX 3.2 software (Volume Graphics, Heidelberg, Germany). Regions with strong signal intensity were assumed to be Lugol solution-permeated.

Seed surface morphologies were analyzed using a JEOL JSM-5610 instrument (LV, Tokyo, Japan). Samples were coated using a fully automated sputter coater (JFC-1600, JEOL Ltd., Tokyo, Japan) with fine-grained platinum for 100 s at 80 mA. The hilum side and the lens groove were photographed at 50X and 500X magnifications, respectively.

Before drying, mature seeds were sectioned at 200 μm using a Plant Microtome MTH-1 (Nippon Medical & Chemical Instruments Co., Ltd., Osaka, Japan). Sections were stained with 0.01% toluidine blue O and observed under an ECLIPSE Ci (OLYMPUS, Tokyo, Japan).

Genomic DNA from wild-type plants was extracted using the NucleoBond HMW DNA extraction kit (Macherey-Nagel, Düren, Germany) and sequenced using the PacBio Sequel II platform (Pacific Biosciences, Menlo Park, CA, USA). HiFi reads were created using the ccs and lima tools (https://github.com/PacificBiosciences/pbbioconda) with “–min-passes 5 –min-rq 0.995” and “–min-score 26 –min-length 10000” options, respectively. De novo assembly of the HiFi reads was performed using Canu ver. 2.0 (Koren et al., 2017), followed by the removal of haplotigs using Purge_Dups (Guan et al., 2020) with default settings.

To construct a chromosome-level assembly of the wild-type genome, we performed RAD-seq using 384 F2 plants derived from a cross between wild-type and another accession JP252958. Genomic DNA was extracted from young leaves using the CTAB method, digested with BglII and EcoRI, ligated with Y-shaped adaptors, amplified using polymerase chain reaction (PCR) with KAPA HiFi HS ReadyMix (Kapa Biosystems Inc., Wilmington, MA, USA), and size-selected using E-Gel size select (Life Technologies, Carlsbad, CA, USA). Approximately 450-bp library fragments were retrieved. Further details of the library preparation have been previously described (Sakaguchi et al., 2015). Sequencing was performed using the paired-end 151-bp mode of HiSeq X (Illumina Inc., San Diego, CA, USA). Quality filtering was performed using Trimmomatic ver. 0.39 (Bolger et al., 2014). Filtered sequence reads were mapped to contigs using BWA-MEM ver. 0.7.17 (Li and Durbin, 2009) with the “-q 20 -F 0×100” option. Single nucleotide polymorphisms (SNPs) were called using Stacks ver. 2.53 (Rochette et al., 2019) using default settings. Missing genotypes were imputed using ABHGenotypeR (Furuta et al., 2017). We removed SNPs for which genotypes were missing in more than 30% of the 384 samples or allele frequencies deviated from Hardy–Weinberg equilibria (p < 1.0E-10). The genetic map was constructed using OneMap software (Margarido et al., 2007) with the “max.rf=0.14” option (for linkage grouping) and the Kosambi map function. SNPs were ordered using three algorithms: “rec,” “rcd,” and “ug.” Chromosome sequences were reconstructed by aligning contigs on the genetic map according to the order of the SNPs.

Protein-coding genes were identified as previously described (Takahashi et al., 2019). For ab initio- and transcriptome-based gene prediction, we used the RNA-Seq data published by Takahashi et al. (2019). We combined protein sequence data deposited in UniProt (UniProt Consortium, 2021) with those of 10 Vigna species downloaded from the Vigna genome server (Sakai et al., 2016): V. angularis; Vigna exilis Tateishi and Maxted; Vigna indica T. M. Dixit, K. V. Bhat, and S. R. Yadav; Vigna marina (Burm. f.) Merr.; Vigna minima (Roxb.) Ohwi and Ohashi; Vigna mungo (L.) Hepper; Vigna riukiuensis (Ohwi) Ohwi and Ohashi; Vigna trilobata (L.) Verdc.; Vigna vexillata (L.) A. Rich.; and Vigna sp. (NI1135 in Meise Botanic Garden, Belgium). These data were used to create a non-redundant set of sequences using MMseqs2 (Steinegger and Söding, 2017) with the “–min-seq-id 0.99 -c 0.99 –cov-mode 0” options. We annotated transposable elements using the Extensive de novo Transposable Element Annotator with default settings (Ou et al., 2019).

V. stipulacea and V. angularis genomes were aligned one-to-one using LAST (Frith and Noé, 2014). Furthermore, unaligned regions (≥100 bp) in one species where flanking sequences (≥1000 bp for each end) were aligned with another species with no large gaps (>100 bp) (PVs candidates) were extracted. PV coordinates were determined by aligning each PV and flanking sequences of the two species using MAFFT (Katoh and Standley, 2013).

F2 plants were produced by crossing the wild-type with the isi2 mutant. In total, 288 F2 plants were cultivated in 1-L pots; DNA was extracted from their young leaves, and water uptake in 20 one-month-old seeds 10 days after watering and 20-seed weight (to estimate 100-seed weight) were measured. Five plants exhibiting the wild-type phenotype that did not produce sufficient seeds were excluded from the evaluation of seed weight. A total of 163 F3 plants derived from 11 F2 plants with recombination among the SNP positions 31,842,848, 32,352,997, and 32,970,540 nt on chromosome 4 were cultivated in 0.2-L pots; DNA was extracted from their young leaves, and water uptake in 10 2-week-old seeds 28 days after watering, 10-seed weights (to calculate 100-seed weights), and numbers of seeds per pod for three pods were measured.

We adopted the MutMap strategy (Abe et al., 2012) to identify candidate SNPs for the isi2 phenotype. Two DNA pools of 215 wild-type phenotypes and 73 isi2 phenotypes in the F2 plants from the wild-type and the isi2 mutant were sequenced using the paired-end 151-bp mode of HiSeq X (Illumina). Sequence reads were quality-filtered using Trimmomatic and mapped to the genome assembly using BWA-MEM. SNPs were identified using GenotypeGVCFs of GATK ver. 4.2.1.0 (Van der Auwera and O’Connor, 2020) and filtered with the “QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0” options of VariantFiltration implemented in GATK. The retained variants were used to calculate SNP indices. Each variant was functionally annotated with the SnpEff ver. 4.3t (Cingolani et al., 2012) using annotated gene structures.

The genotypes for the SNPs at 31,842,848 nt, 32,352,997 nt, and 32,970,540 nt on chromosome 4 were determined for 288 F2 and 163 F3 plants. Primers used for PCR and Sanger sequencing are shown in Table 1 . For Sanger sequencing, we amplified template DNA using AmpliTaq Gold 360 Master Mix (Thermo Fisher Scientific, Waltham, MA, USA) followed by sequencing with BigDye Terminator v3.1 (Thermo Fisher Scientific) and an ABI Genetic Analyzer 3130xl (Thermo Fisher Scientific) according to the manufacturer’s protocols.

Young pods of approximately 5 mm in length 3 days after flowering and unexpanded third leaves were harvested. Total RNA was extracted using the RNeasy Plant Mini Kit (QIAGEN, Hilden, Germany). Vigst.04G198600 cDNA was amplified using the PrimeScript One Step RT-PCR Kit Ver. 2 (TAKARA Bio Inc., Shiga, Japan) with the primers shown in Table 1 . DNA extracted from unexpanded third leaves was used as a negative control. Electrophoresis was performed using 2% agarose gels.

To identify Vigst.04G198600 orthologs in the genomes of closely related legume species, we obtained protein sequences from 11 Vigna species from the Vigna genome server as described above, and Phytozome (for V. unguiculata), as well as those of P. vulgaris, G. max, and Glycine soja Siebold and Zucc. from Phytozome (Goodstein et al., 2012). We ran Orthofinder ver. 2.5.4 (Emms and Kelly, 2019) using default settings to identify the single orthogroup including Vigst.04G198600. We measured collinearity between V. stipulacea and each of the legume species using MCScanX and confirmed that all genes in the orthogroup were syntenic with Vigst.04G198600. To construct a phylogenetic tree of orthologs, we created codon alignments and generated a neighbor-joining tree using the Nei–Gojobori method in MEGA ver. 10 (Kumar et al., 2018). Amino-acid sequences were aligned using the Clustal Omega multiple sequence alignment tool (http://www.ebi.ac.uk/Tools/msa/clustalo/).

Seeds preserved at 4°C for one, six, and twenty-four months after harvest were submerged in water. In wild-type V. stipulacea, water uptake at one month after watering was 1% in 1-month-old seeds, 11% in 6-month-old seeds, and 28% in 2-year-old seeds. These results indicate that water uptake was low in young seeds but higher in older wild-type seeds ( Figure 1A , Supplementary Table 1 ). In contrast, water uptake in isi2 mutant seeds was consistently higher at 96–98% a month after submersion regardless of the storage period but tended to be higher for older seeds until 3 days of submersion ( Figure 1A , Supplementary Table 1 ).

In addition, CT scanning identified water entry at the lens in the isi2 mutant ( Figure 1B, C ). Longitudinal sections showed entry into the cotyledon through the lens ( Figure 1D, E ). SEM revealed that the seed coat microstructure of isi2 was different from that of the wild-type. The mutant had decreased honeycomb-like epicuticular wax deposition in the lens groove compared to the wild-type; however, no differences were observed in the entire seed coat ( Figure 1F-I , Supplementary Figure 1 ). Similarly, optical microscopy revealed no reproducible differences in seed coat cross sections between the mutant and the wild-type ( Figure 1J, K ). Three biological replicates are shown in the Supplementary Figure 2 .

We sequenced the wild-type genome and obtained 20.62 Gbp of HiFi reads. De novo assembly generated 823 contigs with a total length of 413.48 Mbp and an N50 length of 11.87 Mbp ( Table 2 ). We constructed a genetic map of the 11 chromosomes based on 1,163 SNP markers derived from RAD-seq and anchored 44 contigs to it. These anchored contigs spanned 387.57 Mbp, covering 88% of the estimated genome size (441.54 Mbp). The repetitive fraction of the genome was 38% (156.62 Mbp), smaller than that of the V. angularis genome (45%). The most abundant class of repetitive elements was Gypsy type LTR-retrotransposons (7.5%), followed by unknown (6.8%) and Copia type (6.4%) LTR-retrotransposons.

We identified 30,963 protein-coding genes by combining ab initio-based, transcriptome-based, and protein-based approaches ( Table 2 ). BUSCO analysis of the embryophyta_odb10 dataset showed 98.7% completeness in the predicted proteins, slightly higher than our previous version (95.8%; Takahashi et al., 2019).

In addition, we assessed synteny between the genomes of V. stipulacea and azuki bean (Vigna angularis (Willd.) Ohwi et H. Ohashi) and between V. stipulacea and P. vulgaris. The genomes of V. stipulacea and V. angularis were highly syntenic across all 11 chromosomes. We detected 669 synteny blocks between V. stipulacea and V. angularis with 23,574 (76%) of the 30,963 genes overlapping. In contrast, the genomes of V. stipulacea and P. vulgaris had several structural differences, including translocations and inversions ( Figure 2 ).

Despite the high synteny between V. stipulacea and V. angularis, precise comparisons based on the whole-genome alignment revealed several structural variations. We detected 6,298 presence variations (PVs) in the V. stipulacea genome, with a total length of 5,094,208 bp, whereas the V. angularis genome harbored 8,090 PVs with a total length of 9,822,523 bp. Of the 5,094,208 bp of PVs in V. stipulacea, 3,154,013 (62%) were annotated as repetitive elements whereas 27 protein-coding genes were annotated in the PV regions.

Upon testing the genetic factor, we concluded that the isi2 phenotype was caused by a mutation in a single locus according to the law of independence. We defined the isi2 phenotype as one in which more than 80% of seeds exhibited water uptake within one to four weeks after watering. Among F2 (WT × isi2) plants, 215 showed the wild-type phenotype and 73 showed the isi2 phenotype, consistent with the expected 3:1 Mendelian ratio (216:72; p = 0.89, based on the Chi-square test; Supplementary Table 2 ).

For bulked segregant analysis, we pooled F2 (WT × isi2) plants into wild-type and isi2 phenotypes and sequenced each of them. We identified 61,664 SNPs across 11 chromosomes with an average depth of 64.6 for the wild-type pool and 70.8 for the mutant pool. Of these, 9 located between 31,842,848 nt and 35,015,640 nt were completely or nearly fixed with non-reference alleles ( Figure 3 , Supplementary Figure 3 , Table 3 ). Of the nine SNPs, two were located in the coding regions of Vigst.04G193900 and Vigst.04G214200 and were predicted to cause single-residue substitutions whereas one was located in the splice donor site of intron 10 of Vigst.04G198600 ( Figure 4A , Table 3 ). The other six were located in intergenic regions or introns and were not further characterized.

To determine which of the three SNPs is responsible for the isi2 phenotype, we genotyped 288 F2 (WT × isi2) and 163 F3 (WT × isi2) plants for the T>A SNP at 31,842,848 nt, the C>T SNP at 32,352,997 nt, and the C>T SNP at 32,970,540 nt on chromosome 4. Complete linkage was found only for the C>T SNP at 32,352,997 nt, whereas recombination was observed between this SNP and both of the other two SNPs ( Table 4 , Supplementary Table 2 ). Data from 163 F3 (WT × isi2) plants derived from 11 F2 (WT × isi2) plants with recombination among the three SNPs were consistent with the F2 data ( Table 5 , Supplementary Table 3 ). Thus, the C>T SNP at 32,352,997 nt, located in the splice donor site of the 10^th intron in Vigst.04G198600, was the single remaining candidate for the isi2 phenotype.

A BLAST search identified Vigst.04G198600 as a putative ortholog of the A. thaliana gene encoding phospholipid sterol acyltransferase 1 (AtPSAT1/ERP1, AT1G04010). Therefore, we renamed it VsPSAT1. Querying the V. stipulacea genome with the amino-acid sequence of AtPSAT1 identified only VsPSAT1. The predicted VsPSAT1 protein is 623 residues in length, with an overall identity of 76% with AtPSAT1. The conserved LCAT-like domain (Banas et al., 2005; Lara et al., 2018) has a total of 41 residues in six regions, of which 5 residues in 3 regions were replaced by strongly similar residues, and all 3 residues of the catalytic triad Ser-Asp-His were identical ( Figure 5 ).

We performed RT-PCR to determine whether the SNP at the splice donor site affected the mRNA sequence of VsPSAT1. The region between the exons 8 and 11 in the isi2 mutant was approximately 50 bp shorter than that in the wild-type ( Figure 4B ). Sequencing revealed that exon 10 had been skipped ( Supplementary Figure 4 ), leading to a frameshift and an ectopic stop codon that truncated the protein from 623 to 452 residues ( Figure 4C ).

We next searched for putative orthologs in species of the tribe Phaseoleae Bronn ex DC., specifically other Vigna species, P. vulgaris, and G. max. The diploids (Vigna species and P. vulgaris) had one copy and the palaeopolyploid soybean had two copies ( Figure 4D ), without any significant degradation in the conservation of the LCAT-like proteins, including the catalytic triad ( Supplementary Figure 5 ).

As the A. thaliana AtPSAT1 mutant exhibits early leaf senescence due to overaccumulation of sterol esters (Bouvier-Navé et al., 2010), we investigated whether this phenotype occurs in the isi2 mutant plants. Browning was observed at isi2 leaflet margins four weeks after detached leaves were floated in water ( Figure 6A, B ). The SPAD value, indicating chlorophyll content, was significantly (p < 0.05, t-test) different between the wild-type (27.3; SD = 1.37) and isi2 mutant (20.5; SD = 1.81); however, it was not significantly different between both intact leaves ( Supplementary Table 4 ).

We next evaluated the isi2 mutant for eight traits to elucidate whether it showed pleiotropic effects and found that days to flowering increased by 3 days, the largest leaflet width increased by 7%, the number of seeds per pod decreased by 50%, the pod length decreased by 23%, and the 100-seed weight increased by 37%, compared to the wild-type ( Figure 6C-E , Table 6 ). No significant differences were observed in height, largest leaf length, and shoot dry weight between the wild-type and isi2 mutant (p > 0.05, t-test).

We next determined if the isi2 VsPSAT1 mutation was genetically linked with the 100-seed weight and number of seeds per pod. The mean 100-seed weight in 73 isi2 F2 (WT × isi2) plants was 1.38 g, significantly higher than that of 210 wild-type F2 (WT × isi2) plants (1.19 g; p < 0.05, t-test; Figure 6 , Supplementary Table 2 ). In addition to seed size, the number of seeds per pod was tightly linked to the isi2 phenotype in 163 F3 (WT × isi2) plants derived from 11 F2 (WT × isi2) plants with recombination among the 3 SNPs ( Table 5 ). However, the linkage of 100-seed weight and number of seeds per pod to the isi2 phenotype was not complete on an individual basis, owing to high variability ( Supplementary Table 3 ).