Table 1 summarizes the materials tested in this study. Besides V. stipulacea, we used two or three accessions of three domesticated species, soybean [Glycine max (L.) Merr.], common bean (Phaseolus vulgaris L.), and cowpea [Vigna unguiculata (L.) Walp.], to evaluate pod shattering and seed dormancy. All the materials were provided by the NARO gene bank (https://www.gene.affrc.go.jp/index_en.php). We cultivated V. stipulacea in a field or in a bucket with gardening soil in a greenhouse of our institute, Tsukuba, Japan (36.030577, 140.098021). We cultivated three plants of the domesticated species in a bucket with gardening soil in a greenhouse of our institute, Tsukuba, Japan.

The process of mutant screening is shown in Figure 1 . To increase mutation density, we adopted two cycles of successive chemical mutagenesis as described by Tsuda et al. (2015). We started with 12,000 seeds of V. stipulacea (JP245503). We first scratched the seed coat with a knife because it is hard and waterproof, treated the scratched seeds with 0.35% EMS solution for 8 h, and thoroughly washed them with distilled water (M1 seeds). The M1 seeds were incubated in a wet condition for 3 days and those with germination were sown. The M1 plants were cultivated in our greenhouse, where the temperature was kept above >20°C. Three pods per plant were harvested from 3,000 fertile M1 plants (M2 seeds). Of the harvested seeds, four per line (12,000 in total) were scratched and treated with 0.35% EMS solution as described above (M2M1 seeds). Germinating M2M1 seeds were sown and cultivated in a field, located in Tsukuba, Japan. Three thousand fertile plants were selected and three pods per plant were harvested (M2M2 seeds). Six seeds per M2M2 line were sown in a 1 L plastic pot filled with gardening soil and cultivated in the greenhouse until three pods per line were harvested (M2M3 seeds). Of the M2M3 lines, four lines with mutant phenotypes were selected and seven seeds per line were sown in a 7 L bucket filled with the gardening soil and cultivated in the greenhouse.

To evaluate pod shattering, we harvested 20 pods of each domesticated accession and three pods per M2M2 plant and counted the number of shattered pods. Before evaluation, harvested pods were left at room temperature for a month and then completely dried at 40°C for 24 h in an incubator.

To evaluate seed dormancy, 20 seeds were soaked in distilled water and incubated in a dark incubator at 25°C, and imbibed seeds were counted daily until the 3^rd day, twice a week until the 14^th day, and weekly until the 28^th day. We replicated the measurements by seven times and performed Dunnett test to compare averages of each time point between the wild type and each mutant line.

We observed the imbibing seeds to detect the imbibition start point. We also observed pod sclerenchyma (tissue with dead cells with thickened secondary wall) with a stereoscopic microscope (ECLIPSE Ci-L, Nikon, Tokyo, Japan). Three pods were collected from each individual and were sliced by a microtome (MTH-1, Nippon Medial & Chemical Instruments Co., Ltd. Osaka, Japan), and stained with Phloroglucinol-HCL solution (1 g phloroglucin in 50 mL ethanol + 25 mL concentrated hydrochloric acid).

We sequenced the whole genome of V. stipulacea (JP245503) with RSII sequencer (Pacific Biosciences, Menlo Park, CA), as we have done previously for azuki bean (Sakai et al., 2015). DNA was isolated from 1 g of unexpanded leaves with CTAB method and purified with Genomic Tip 20/G (Qiagen K. K. Tokyo). The extracted DNA was sheared into 20 kb fragments using g-TUBE (Covaris, MA, USA) and converted into 20 kb SMRTbell template libraries. The library was size-selected for a lower cutoff of 10 kb with BluePippin (Sage Science, MA, USA). Sequencing was performed on the PacBio RS II using P5 polymerase binding and C3 sequencing kits with 360 min acquisition. In total, 52 SMRT cells were used to obtain ∼19.6 Gb of subreads.

In total, ∼3.3 million PacBio reads were used for de novo assembly with Celera Assembler 8.3rc1 (asmOvlErrorRate = 0.1, asmUtgErrorRate = 0.06, asmCgwErrorRate = 0.1, asmCnsErrorRate = 0.1, asmObtErrorRate = 0.08, utgGraphErrorRate = 0.05, utgMergeErrorRate = 0.05) (Berlin et al., 2015). About 25.2x of the longest error-corrected and trimmed reads were assembled to contigs. Redundant contigs were discarded by conducting all-to-all BLASTN searches. The non-redundant contigs were polished by PacBio subreads by using Quiver in SMRT Analysis v2.2.0 (Pacific Biosciences of California, Inc.) and then further polished by Illumina short-reads using BWA 0.7.9a (Li and Durbin, 2009), Samtools 0.1.19 (Li et al., 2009), Picard 1.94 (http://picard.sourceforge.net/), and GATK 3.3 (McKenna et al., 2010). The polished contigs were scaffolded by Reference-Assisted Chromosome Assembly (RACA) program v.0.9.1.1 (Kim et al., 2013b) using the genome sequences of Vigna angularis and P. vulgaris as the reference and outgroup species, respectively.

Repetitive sequences in the genome assembly were predicted using Censor (Kohany et al., 2006) with a composite library consisting of de novo created library constructed by RepeatModeler 1.0.8 (http://www.repeatmasker.org) and the MIPS Repeat Element Database ver. 9.3 (Nussbaumer et al., 2013).

Ab initio gene prediction was done by BRAKER version 1.6 (Hoff et al., 2016) with RNA-Seq data. Besides, gene structures were predicted by genome-guided and de novo RNA-Seq data assembly approaches using TopHat 2.1.0 (Kim et al., 2013a), Cufflinks 2.2.1 (Trapnell et al., 2010), Trinity 2.1.1 (Grabherr et al., 2011), and PASA pipeline 2.0.2 (Haas et al., 2003; Rhind et al., 2011). Open reading frames were predicted by Transdecoder 2.0.1 and Trinotate 2.0.2 (Bryant et al., 2017). Protein sequences of the G. max (Wm82.a2.v1), P. vulgaris (v.1.0), Medicago truncatula Gaertn. (Mt4.0v1), and V. angularis (Willd.) Ohwi & H.Ohashi (VANGULARIS_V1.A1) were downloaded from Phytozome (G. max, P. vulgaris, and M. truncatula) (Goodstein et al., 2012) and VigGS (V. angularis) (Sakai et al., 2016) and mapped to the genome assembly by Exonerate 2.2.0 (Slater and Birney, 2005). The ab initio gene models and transcript and protein alignments were then combined by EvidenceModeler 1.1.1 (Haas et al., 2008) and the predicted gene models were updated by PASA. Gene models with extremely long introns and those merged by PASA were manually curated. BUSCO v3 (Waterhouse et al., 2017) was used to evaluate protein sequences of annotated genes.

We performed MutMap (Abe et al., 2012) to isolate the candidate gene for the mutant phenotype. We crossed one of the mutants of increased seed imbibition (isi1) to the wild type (JP245503) and obtained 280 F2 seeds. We cultivated all of them, extracted DNA from an unexpanded leaf of each F2 plants by CTAB method, and obtained F3 seeds. To evaluate seed dormancy, 20 seeds of each F3 line was soaked with distilled water for 24 h to test imbibition. DNA of 59 F2 plants with mutant phenotype (100% seed imbibition) and DNA of 221 F2 plants with wild type phenotype (0% imbibition) were pooled to make the MT pool and the WT pool, respectively, and were sequenced with Illumina HiSeq X (Illumina Co. Ltd, San Diego, USA).

The obtained sequences were mapped to the draft genome sequences with bwa-0.7.17 (Li and Durbin 2009) and formatted with samtools-1.9 (Li et al., 2009) as described in the distributed manuals. SNPs were called with bcftools-1.9, and the DP4 values were used to calculate the SNP-index. If a SNP site is responsible for the mutant phenotype, the locus must be fixed with the alternative allele [a] in the MT pool ([aa] shares 100%). On the other hand, in the WT pool, the expected ratio of the reference allele [A] and the alternative [a] is 2:1, because the WT pool [AA]:[Aa]:[aa] should be 1:2:0. Thus, at the responsible locus, the expected values of SNP-index are 1.0 in the MT pool and 0.33 in the WT pool. Thus, after calculating SNP-index, we screened all the SNP sites for those that met the criterion above.

We genotyped 280 F2 plants by Sanger sequencing for two SNPs whose SNP-index was 1.0 in the MutMap analysis. Two primer pairs were designed to amplify a 232 bp genomic region including 6,804,429 nt (5′-gagggaatacgaagagtttaaggtt-3′ and 5′-ttgaaaaccaggtcttttctctcta-3′) and to amplify a 243 bp genomic region including 7,009,873 nt (5′-acagagcaaaagattaaacgagaga-3′ and 5′-aaagccgcttcctagtccttac-3′) on scf0015. For Sanger sequencing, we amplified the template DNA with AmpliTaq Gold 360 Master Mix (Thermo Fisher Scientific K. K., Tokyo), performed sequencing reaction with BigDye Terminator v3.1 (Thermo Fisher Scientific K. K., Tokyo), and sequenced with ABI Genetic Analyzer 3130xl (Thermo Fisher Scientific K. K., Tokyo), according to the provider’s protocol.

Though many studies have investigated seed dormancy in legume crops, there are still some arguments on the sites of water entry (Smýkal et al., 2014). Thus, we evaluated seed imbibition and pod shattering in soybean, common bean, and cowpea.

As a result, domesticated species showed variation in both imbibition start sites and time to imbibe ( Figure 2 ). Among the domesticated accessions tested, soybean imbibed within a few minutes and was the quickest. In soybean, the whole testa seemed permeable. The second quickest was the cowpea “JP244182,” where the testa was partially permeable and imbibition was observed in an hour. The common bean “JP41232” showed imbibition in 2 to 3 h, which started at the site of micropyle. It took 4 to 5 h for the common bean “JP41234” and the cowpea “JP239215,” where imbibition started at the site of lens. All the seeds were fully imbibed within 48 h.

We also observed abscission layers between the valves of seed pods and the sclerenchyma on the endocarp, because pod shattering is dependent on these tissues (Murgia et al., 2017). In soybean, as described by Dong et al. (2014), abscission layer was not completely formed at the fiber cap cells ( Figure 3 ). Abscission layer was even less completely formed in the common bean accession JP41963 whereas it was fully formed in all other accessions ( Figure 3 ). Pod sclerenchyma was thicker in the soybean accessions and the cowpea accession JP244182, whereas it was thinner in common bean accessions ( Figure 3 ).

To obtain mutants with reduced seed dormancy and reduced pod shattering, we treated 12,000 seeds of V. stipulacea with EMS and developed 3,000 M2M2 plants. Of the 3,000, we found one line which exhibited reduced pod shattering and designated it as reduced pod shattering1 (rps1) ( Figure 4 ). We also screened M2M3 seeds for seed imbibition and obtained three lines which exhibited increased seed imbibition and designated them as increased seed imbibition1 (isi1), increased seed imbibition2 (isi2), and increased seed imbibition3 (isi3), respectively ( Figure 5 ). Of them, the isi1 seeds had cracks in the hilum, and the isi3 seeds showed reduced pigmentation in the seed coat ( Figure 5 ).

The mode of imbibition also varied across the mutants ( Figure 5 ). In the isi1, the water entry was not necessarily through the cracks in the hilum but through the whole seed coat. Even though we sealed the hilum with glue, imbibition was quickly initiated [see the “isi1 (sealed hilum)” in Figure 5 ]. In isi2, the imbibition was initiated at the lens. In isi3, the imbibition started near the hilum, but not necessarily through the micropyle or the lens.

To elucidate the extent of seed dormancy in the mutant lines, we put the seeds on wet filter paper and evaluated the rate of imbibed seeds ( Figure 6 ). We also tested the seeds of 1 and 6 months old to elucidate the effect of duration time after harvesting. When we did the experiment with the wild type seeds, the imbibition rate was zero for at least 4 weeks in 1-month-old seeds and was 11% in the 6-month-old seeds ( Figure 6 ).

When we tested the mutant seeds of 1 month old, the imbibition was quick in isi1 and isi2, whereas it was slow in isi3. In isi1, all the seeds fully imbibed within 1 h ( Figure 6 ). In isi2, 85% of the seeds imbibed within a week, and almost 100% did so within 4 weeks. In contrast, the imbibition rate of the isi3 seeds was only 16% in 4 weeks, which was not significantly different from the wild type. The imbibition rate of the reduced shattering mutant rps1 was also higher but not significantly different compared to the wild type ( Figure 6 ).

When we tested the mutant seeds of 6 months old, the overall imbibition rate was increased compared to the experiment with 1-month-old seeds ( Figure 6 ). In isi1, the seeds also fully imbibed within 1 h. In isi2, the imbibition rate was already significantly different from the wild type in 1 day, more than 90% of the seeds imbibed within a week, and reached plateau (∼95%) in two weeks. In isi3, 20% of the seeds imbibed in 3 days, 36% in a week, and 71% in 4 weeks. However, this line exhibited a large standard deviation between the replicates and was not significantly different before 2 weeks ( Figure 6 ). In rps1, ∼20% of the seeds imbibed in 4 weeks, which was twice as high as the wild type but not significantly different from the wild type was not significant.

To evaluate pod shattering in the mutant lines, we calculated the rate of shattering of the harvested pods which were completely dried in the incubator. Whereas the shattering rate was 100% in the wild type, it was 0% in the rps1 mutant ( Table 2 ). The rps1 mutant also showed a reduced twisting of the seed pod. The number of twists/cm in the pods was 0.371 ± 0.018 in the rps1 mutant, which was less than half of the wild type (0.866 ± 0.022) ( Table 2 ).

Interestingly, the isi1, one of the mutants of seed imbibition, also exhibited slightly reduced shattering rate (73.99 ± 18.53%) and number of twists/cm (0.579 ± 0.093) ( Table 2 ). Other mutants were also slightly reduced in number of twists/cm, but their shattering rate was 100% ( Table 2 ).

We also observed cross-sections of seed pods and found the rps1 mutant did not form abscission layer between the valves at all ( Figure 4 ).

To build a reference sequence of V. stipulacea, we sequenced the genomic DNA with a PacBio sequencer, assembled, and annotated. We obtained 19.6 Gbp of subreads using 52 SMRT cells. Genome size of the V. stipulacea was estimated to be ∼445.1 Mbp based on the k-mer frequency distribution obtained from 10.3 Gbp of the Illumina short-reads. The assembled contigs showed N50 length of ∼1.9 Mbp and mean length of ∼169.8 kbp and covered ∼387.7 Mbp (87.9%) of the estimated genome size ( Table 3 ). Out of the 2,102 scaffolds, 52 scaffolds were comprised of 233 contigs scaffolded by RACA.

By combining the ab initio gene models and transcript and protein alignments, 26,038 protein coding genes were predicted on the assembled genome ( Table 3 ). We assessed the completeness of the gene set by BUSCO and found that 95.8% of the BUSCO gene models were completely detected ( Table 3 ). The BUSCO metrics was almost the same as that assessed for the recently published high quality genome sequence of cowpea (Lonardi et al., 2019), suggesting that our genome assembly covered nearly complete gene space of the V. stipulacea genome.

To identify the responsible genes for the mutant phenotype, we performed MutMap analysis. Of the 280 F2 plants derived from isi1 x WT, 221 showed the wild type phenotype and 59 showed the mutant phenotype, which was close to the expected 3:1 ratio (p = 0.13). We sequenced the genomic DNA of the wild type plants (WT pool) and the mutants (MT pool) and identified 33,936 SNPs across the whole genome. We then calculated SNP-index of each SNP and screened for those where it was 1.0 in MT pool and 0.3 ± 0.1 in WT pool.

As a result, we found only two on scf0015, a C to T substitution at 6,804,429 nt and a G to A at 7,009,873 nt. The C to T SNP was present in the 8^th exon of Vigst.0015s042600.01, which encoded CELLULOSE SYNTHASES A 7 (CesA7) protein, and turned the Glu codon into a STOP codon. The G to A SNP was present in Vigst.0015s044500.01, which encoded a nucleoporin-like protein, but was synonymous.

For further mapping, we looked for recombinants between the two candidate SNPs by direct sequencing of the SNP-containing regions in all the F2 plants. As a result, one F2 plant, which showed the wild type phenotype, was heterozygous at Vigst.0015s042600.01 but was fixed with the mutant A allele at Vigst.0015s044500.01 (see #130 in Supplementary Table 1 ). Another F2 plant, which showed the mutant phenotype, was fixed with the alternative T allele at Vigst.0015s042600.01 but was heterozygous at Vigst.0015s044500.01 (see #251 in Supplementary Table 1 ). The result was further confirmed by genotyping and phenotyping F3 plants derived from the F2 plants with recombination between the two candidate SNPs (#130, #251, and #275 in Supplementary Table 1 ). As shown in Supplementary Table 2 , the phenotype was completely linked with the SNP at scf0015_6804429 on CesA7 gene but not with the other.