All the accessions used in this study were provided by NARO Genebank (Tsukuba, Japan) (Figure 1 and Supplementary Figure 1).

For azuki bean, we started from the BC₁F₂ plants derived from a cross between domesticated azuki bean (JP81481, the recurrent parent) and a wild relative, Vigna nepalensis Tateishi & Maxted (JP107881) (Isemura et al., 2007). Of them, we selected one BC₁F₂ plant where the locus for pod shattering is fixed with the wild-type allele but most of other loci are fixed with domestication-type alleles and crossed again to the recurrent parent. We further backcrossed the obtained BC₂F₁ plants to obtain BC₃F₁ plants. We selfed them and obtained BC₃F₂ seeds from those with pod shattering phenotypes. For BC₃F₃, BC₃F₄, and BC₃F₅ populations, we kept selecting and selfing those with recombination within the candidate region. We also included plants that were heterozygous throughout the candidate region for obtaining more recombination and some that are homozygous in the same region to reconfirm the relationship of genotype and phenotype.

For yard-long bean, we started from BC₁F₂ plants derived from a cross between yard-long bean (JP89083, the recurrent parent) and a wild cowpea, Vigna unguiculata subsp. dekindtiana (Harms) Verdc. (JP81610) (Kongjaimun et al., 2013). Of them, we selected one BC₁F₂ plant that are fixed with wild-type allele at the locus for pod tenderness but are fixed with domestication-type at most other loci. The selected BC₁F₂ plants were backcrossed to the recurrent parent and further backcrossed to obtain BC₃F₁ plants. We selfed them and obtained BC₃F₂ seeds from those with hard pods phenotype. We kept selecting and selfing up to BC₃F₅ populations as done in azuki bean.

All the plants were grown from July through November in a field located in Tsukuba city, Japan, except BC₃F₅ plants of yard-long bean that were grown in a greenhouse from September through January.

We observed sclerenchymal tissues in seed pods by staining lignin with the Wiesner (phloroglucinol–HCl) reaction (Pormar et al., 2002). We harvested three undried mature seed pods per plant, sliced them to 200-μm sections with a Plant Microtome MTH-1 (Nippon Medical & Chemical Instruments Co., Ltd., Osaka, Japan), incubated in phloroglucinol–HCl solution for 1 min and observed with stereoscopic microscope SZX7 (OLYMPUS, Tokyo) for BC₃F₅ plants of azuki bean and yard-long bean. Phloroglucinol–HCl solution was prepared by dissolving 1 g of phloroglucin in 50 ml ethanol and adding 25 ml of concentrated hydrochloric acid.

For azuki bean populations, we evaluated pod shattering by counting the number of twists in a pod, as described in Isemura et al. (2007). For BC₃F₂ and BC₃F₃, we harvested five pods per plant before shattering and measured the length. We then incubated the pods at 60°C to let them totally dry. Unless the pods dehisced, we manually opened the pods and counted the number of twists and calculated the number of twists/cm. For BC₃F₄ and BC₃F₅, we harvested 10 pods per plant, let them dry, and evaluated the rate of shattering pods before counting the number of twists.

For yard-long bean, we manually evaluated the tenderness of five young pods per plant by binarizing hard or soft.

We also measured 100 seed weight on BC₃F₂ and BC₃F₅ plants of the azuki bean population and BC₃F₄ and BC₃F₅ plants of the yard-long bean population.

We sequenced the whole genome of wild cowpea (JP81610) with RSII sequencer (Pacific Biosciences, Menlo Park, CA, United States), as we have done previously for azuki bean (Sakai et al., 2015). DNA was isolated from 1 g of unexpanded leaves with the cetyl trimethylammonium bromide (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, United States) 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, United States). Sequencing was performed on the PacBio RS II using P6 polymerase binding and C4 sequencing kits with 360 min acquisition. In total, 21 SMRT cells were used to obtain ∼26.4 Gb of raw reads.

In total, 4 million PacBio reads were used for de novo assembly with Canu v1.6 under the default settings (corOvlErrorRate = 0.2400, obtOvlErrorRate = 0.0450, utgOvlErrorRate = 0.0450, corErrorRate = 0.3000, obtErrorRate = 0.0450, utgErrorRate = 0.0450, cnsErrorRate = 0.0750). About 23.5x error corrected and trimmed reads longer than 1,000 bp were assembled to contigs. Assembled contigs were polished by PacBio raw reads by using the arrow in GenomicConsensus v2.3.2 (Pacific Biosciences of California, Inc.).

Repeat detection was conducted by RepeatMasker ver. 4.0 (Smit et al., 2013)¹. A de novo repeats library of wild cowpea genome constructed by RepeatModeler ver. 1.0.11 (see footnote 1) and the Fabaceae repeats library in RepBase24.02 (Smit and Hubley, 2008)² were used for the prediction.

Ab initio gene prediction was done with AUGUSTUS (version 3.3.2) (Stanke et al., 2008). A set of gene annotation information of recently published cowpea genome was used for training AUGUSTUS (Lonardi et al., 2019). We trained a new model twice using 1,000 high-confidential genes selected based on abundance of annotations of domains, pathway networks, and gene ontology information. BUSCO v3 (Waterhouse et al., 2017) was used to evaluate protein sequences of annotated genes.

Single-nucleotide polymorphisms (SNPs) and short insertions/deletions (indels) were detected by using MUMmer v3.23 (Kurtz et al., 2004). Genome alignment of our wild cowpea assembly and the reference cowpea genome (Lonardi et al., 2019) was conducted by the nucmer command with the following options: –maxgap = 500 –mincluster = 100. Then one-to-one alignment was extracted with delta-filter command with an option: −1. Based on the alignments, SNPs and indels were reported by the show-snps command.

Raw sequence data, the assembled sequences, and gene annotations are all available from DNA Data Bank of Japan (DDBJ)³ under the BioProject ID PRJDB8129.

We extracted DNA from the seeds as described by Kamiya and Kiguchi (2003) and genotyped them by fragment analysis with capillary electrophoresis for simple sequence repeats (SSRs) and INDELs as described by Isemura et al. (2007) or by directly sequencing SNP sites (see section “Direct Sequencing” below). The information of the markers we used is summarized in Supplementary Table 1.

We designed primers by Primer3 (Untergasser et al., 2012) to amplify polymorphic sites between the domesticated azuki bean and V. nepalensis that are available in Vigna Genome Server (VigGS⁴; Sakai et al., 2016) and those between the domesticated cowpea genome in Legume Information System⁵ and the wild cowpea genome sequenced in this study. Parameters for designing primers on Primer3⁶ were 20–30 bp in length, 55°C–65°C in annealing temperature, and 100–350 bp and >700 bp in expected length of amplified fragments for fragment analysis and direct sequencing, respectively.

For azuki bean, the markers we used were as follows:

For yard-long bean, the markers we used were as follows:

To sequence the (potentially) polymorphic sites, 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.

To sequence both the domesticated and the wild alleles of Vigan.07g034400.01 and Vigun05g273500 genes, we sequenced the transcribed sequences of both loci. We extracted total RNA from 100 mg of flower buds right before flowering of all the parental accessions using RNeasy Plant Mini Kit (QIAGEN) with RNase-free DNase I (Invitrogen). Total RNA of 1 μg was converted into first-strand cDNA with Super Script II Reverse Transcriptase (Invitrogen) and Oligo(dT)12-18 Primer (Invitrogen) following the manufacturer’s instructions. The cDNA sequences of Vigan.07g034400.01 and Vigun05g273300.01 were then amplified and sequenced with the primers (Supplementary Table 1) and then transferred to the direct-sequencing protocol described above.

To confirm that azuki bean and yard-long bean have lost or reduced sclerenchyma in seed pods, we observed cross sections of seed pods stained with phloroglucinol–HCl and found a clear-cut difference between the wild and domesticated accessions. In the wild azuki bean and cowpea, which have shattering phenotypes, thick layers of sclerenchyma (∼0.3 mm) were formed on the endocarps of seed pods (Figure 1). However, in the domesticated azuki bean and yard-long bean, which are both non-shattering, the sclerenchyma layer was slightly formed (<0.1 mm) or not formed at all, respectively (Figure 1).

We also evaluated the correlation between the thickness of pod sclerenchyma and helical shape change in seed pods (number of twists/cm) in the BC₃F₅ plants of azuki bean and yard-long bean. As a result, all the shattering plants in azuki bean population formed thicker sclerenchyma (0.15–0.20 mm) and showed a stronger helical shape change (0.30–0.43) than the non-shattering plants did (0.00–0.08 mm and 0.00–0.05, respectively) (Figure 1 and Supplementary Figure 1). The correlation coefficient between the sclerenchyma thickness and the helical shape change was 0.92. In the yard-long bean population, all the plants with hard pod phenotype formed sclerenchyma (0.16–0.22 mm) and showed little helical shape change (0.09–0.12), whereas those with soft pod phenotype did not form sclerenchyma or show helical shape change at all (Figure 1 and Supplementary Figure 1). The correlation coefficient between the sclerenchyma thickness and the helical shape change was 0.99.

In addition, we observed cross sections of seed pods in a wild soybean and three domesticated soybean cultivars and found thick sclerenchyma layers in all the accessions (Supplementary Figure 2).

We previously mapped the QTL for pod shattering in between CEDG182 and CEDG174 on LG7, which is around 1.4–4.0 Mbp in Chr7 of the reference azuki bean genome (Sakai et al., 2015) (Figure 2). To more finely locate the genetic factor for pod shattering, we genotyped 1,049 BC₃F₂ plants and selected 238 for phenotyping. The obtained phenotype and genotype data revealed the pod shattering factor was completely linked with CEDG064 and located in between SPD03 and SPD04 (Figure 2 and Supplementary Table 2). We selected 46 plants out of the 238 (Supplementary Table 2), selfed them, and obtained 4,222 BC₃F₃ seeds.

Of the BC₃F₃ seeds, we selected 53 that had recombination between SPD03 and SPD04. Phenotyping and genotyping on these 53 revealed that pod shattering factor was still completely linked with CEDG064 and was in between SPD08 and SPD09 (Figure 2 and Supplementary Table 3), which is about 19 kb long and contained three genes (Figure 2).

Since this region seemed to have a relatively higher recombination rate, we expected it would be possible to further locate the pod shattering factor. Thus, of the 53 BC₃F₃ plants, we selected 24 plants (Supplementary Table 3) and selfed them to obtain BC₃F₄ seeds. We then cultivated 1–7 plants/line (81 in total; Supplementary Table 4) for further genotyping and phenotyping (Figure 2).

As a result, we again found CEDG064 was completely linked with the phenotype, but we also obtained 29 recombinants between SPD10 and CED064 and five between CEDG064 and SPD11 (Supplementary Table 4). Thus, the pod shattering factor was located in a 9-kb region between SPD10 and SPD11, where there is only one gene Vigan.07g034400 (Figure 2), which encoded VaMYB26, an ortholog of AtMYB26.

To see if there were any loss-of-function mutations in VaMYB26 in azuki bean, we searched the polymorphism data between wild and domesticated azuki bean provided by VigGS. The database showed that there were six SNPs (SNP01-SNP06) and one INDEL (INDEL01) within its open reading frame (ORF). Of these, all the SNPs were in the introns or in the 3′-untranslated region (UTR), but the INDEL was within the coding DNA sequence (CDS) (a T insertion 4 bp before the stop codon in the domesticated azuki bean). Since the wild azuki bean did not have this T insertion, the CDS could be 405 bp (125 aa) longer (Supplementary Figure 3). Interestingly, a BLAST search of MYB26 gene revealed that the longer version is widely conserved across plant taxa (Figure 2). Thus, the VaMYB26 in the domesticated azuki bean seemed to have an immature stop codon.

To confirm the polymorphisms identified in the database, we determined the genomic and the transcribed sequences of VaMYB26 locus of both parents and found all the polymorphisms in the ORF were true (Supplementary Figure 3). We also sequenced the 17 lines of BC₃F₅, which we obtained from the selected BC₃F₄ plants (those with fixed genotypes between SPD03 and SPD09; Supplementary Table 5), to further confirm the relationship between the genotypes and the phenotype. Of the 17, one had recombination between SNP03 and INDEL01, but none had recombination between INDEL01 and CEDG064 (Supplementary Table 5). Thus, the genotypes at INDEL01, SNP04–SNP06, and CEDG064 were completely linked with the pod shattering phenotype, locating the pod shattering factor within the 4-kb region between SNP03 and SPD11 (Figure 2 and Supplementary Table 5).

Though the reference genome of cowpea had been released (Lonardi et al., 2019), it was yet not easy to obtain polymorphic markers because the genome sequence of the wild cowpea was not available. Thus, to facilitate our fine-mapping the factor of pod tenderness, we sequenced and assembled the genome of V. unguiculata subsp. dekindtiana (JP81610), which is a wild cowpea accession used to develop the mapping population of yard-long bean. We used PacBio RSII and obtained 4,413,480 raw reads with an average read length of 6.0 kbp and N50 length of 10.6 kbp, covering 44.9X of the estimated genome size (585.8 Mbp) (Lonardi et al., 2019) (Supplementary Figure 4 and Supplementary Table 6). The assembly produced 4,285 contigs covering 488.5 Mbp (83.4%) of the cowpea genome, with the contig NG50 of 438.6 kbp and the maximum contig length of 2.6 Mbp, respectively (Supplementary Table 7). Of the assembled sequences, 45.7% composed of transposable elements, of which LTR retrotransposons share the largest content (20.5% of the assembly) (Supplementary Table 8). Our ab initio gene prediction detected 36,061 protein-coding genes, and 27,345 of those were non-repeat related (Supplementary Table 9). The annotated genes of the wild cowpea contained 93.2% (85.6% complete and 7.6% partial) of the 1,440 plant genes in BUSCO v3 (Waterhouse et al., 2017).

We aligned the genome sequences of the wild cowpea to those of the reference cowpea (Lonardi et al., 2019) and identified 5,661,319 SNPs and 1,626,169 INDELs. As expected, we detected fewer SNPS and INDELs around centromeric and pericentromeric regions but many in chromosome arms, which were presumably gene-rich regions (Supplementary Figure 5).

The QTL for pod tenderness was previously located between cp06388 and VR294 on LG7 (Kongjaimun et al., 2013) (Figure 3), which corresponds to about 47.5–45.5 Mbp region in Chr5 of the cowpea genome (Vigna unguiculata v1.0, NSF, UCR, USAID, DOE-JGI)⁷. To more finely locate the pod tenderness factor, we designed more markers based on the SNPs and INDELs we identified above, genotyped 2,304 BC₃F₄ seeds, and selected 195 for phenotyping (Supplementary Table 10).

The obtained data of genotype and phenotype located the pod tenderness factor between VuPT03 and VuPT04, which was about 47 kbp (Figure 3 and Supplementary Table 10). Although no marker was completely linked with the phenotype, eight plants had recombination between this region (Figure 3 and Supplementary Table 10).

Of the 195 BC₃F₄ plants we tested, we selfed the eight plants and obtained BC₃F₅ for further mapping. As a result, we located the candidate regions in between VuPT08 and VuPT11, which was 13 kbp long and contained three genes, Vigun05g27300, Vigun05g273400, and Vigun05g273500. Of note, Vigun05g273500 encoded VuMYB26 (Figure 3 and Supplementary Table 11).

According to the SNPs and INDELs data, Vigun05g273500 had an A to G substitution which might disrupt the junction site of the first intron and the second exon (Supplementary Figure 6). On the other hand, Vigun05g27300 seemed to have only synonymous SNPs, and Vigun05g273400 was a non-coding gene. All the SNPs were confirmed by directly sequencing the ORFs of these genes. Interestingly, the ORF sequences of the three genes were exactly the same between the reference cowpea and yard-long bean, which means yard-long bean and the reference cowpea share the domestication-type allele in this locus (Supplementary Figure 6).

To test whether the A to G substitution disrupted splicing, we sequenced cDNAs of this locus prepared from both parents. As a result, we found, in the domestication-type allele, the junction site of the first intron and the second exon had shifted downstream by 8 bp, which could result in a frameshift and an immature stop codon in the middle of the second exon (Supplementary Figure 6). This product would encode a protein of 60 aa, which would be 305 aa shorter than that of the wild-type allele.

Because Murgia et al. (2017) suggested the pod shattering phenotypes load an energy cost on plants and limit seed size, we measured 100 seed weight of the mapping population (BC₃F₂ of azuki bean and BC₃F₄ of yard-long bean).

As expected, those with domestication-type trait produced larger seeds compared to those with the wild-type trait (Figure 4). In the azuki bean population, the seeds of non-shattering plants were 15.4 ± 2.4 g/100 seeds, whereas the seeds of shattering plants were 13.3 ± 1.9 g/100 seeds. In the yard-long bean population, the seeds of tender-pod plants were 13.6 ± 1.9 g/100 seeds, whereas the seeds of hard-pod plants were 12.9 ± 1.7 g/100 seeds. The following t-test revealed that in both cases, the seeds of the plants with the domestication-type phenotypes produced significantly larger seeds than the others (p = 4.6 × 10^–7 for azuki bean and p = 0.019 for yard-long bean).

We also measured 100 seed weight of BC₃F₅ plants and observed the same trend, though the differences were not significant (Supplementary Figure 7). In azuki bean, the seeds of non-shattering plants were 9.7 ± 1.4 g/100 seeds, whereas the seeds of shattering plants were 8.7 ± 1.2 g/100 seeds. In yard-long bean, the seeds of tender-pod plants were 14.2 ± 1.5 g/100 seeds, whereas the seeds of hard-pod plants were 12.2 ± 1.2 g/100 seeds.