QTL analysis of domestication syndrome in zombi pea (Vigna vexillata), an underutilized legume crop

Zombi pea (Vigna vexillata (L.) A. Rich) is an underutilized crop belonging to the genus Vigna. Two domesticated forms of zombi pea are cultivated as crop plants; seed and tuber forms. The cultivated seed form is present in Africa, while the cultivated tuber form is present in a very limited part of Asia. Genetics of domestication have been investigated in most of cultivated Vigna crops by means of quantitative trait locus (QTL) mapping. In this study, we investigated genetics of domestication in zombi pea by QTL analysis using an F2 population of 139 plants derived from a cross between cultivated tuber form of V. vexillata (JP235863) and wild V. vexillata (AusTRCF66514). A linkage map with 11 linkage groups (LGs) was constructed from this F2 population using 145 SSR, 117 RAD-seq and 2 morphological markers. Many highly segregation distorted markers were found on LGs 5, 6, 7, 8, 10 and 11. Most of the distorted markers were clustered together and all the markers on LG8 were highly distorted markers. Comparing this V. vexillata linkage map with linkage maps of other four Vigna species demonstrated several genome rearrangements in V. vexillata. QTL analysis for 22 domestication-related traits was investigated by inclusive composite interval mapping in which 37 QTLs were identified for 18 traits; no QTL was detected for 4 traits. Number of QTLs detected in each trait ranged from 1 to 5 with an average of only 2.3. Five QTLs for tuber width and three QTLs for tuber weight. Interestingly, 2 QTLs each for tuber width and tuber weight detected on LG2 and LG4 were located at similar position and wild allele increased tuber width and weight. This indicated wild germplasm having small tuber have potential to increase yield of large tuber cultivated type. Large-effect QTLs (PVE > 20%) were on LG4 (pod length), LG5 (leaf size and seed thickness), and LG7 (for seed-related traits). Comparison of domestication-related QTLs of the zombi pea with those of cowpea (Vigna unguiculata), azuki bean (Vigna angularis), mungbean (Vigna radiata) and rice bean (Vigna umbellata) revealed that there was conservation of some QTLs for seed size, pod size and leaf size between zombi pea and cowpea and that QTLs associated with seed size (weight, length, width and thickness) in each species were clustered on same linkage.


Introduction
The genus Vigna is an important taxon of leguminous plants. This genus comprises more than 100 species that distribute in all major continents including Africa, America, Asia, Australia and Europe. Among those species, as high as nine Vigna species are domesticated/cultivated. These domesticated Vigna species include azuki bean (Vigna angularis (Wild.) Ohwi and Ohashi), black gram (Vigna mungo (L.) Hepper), créole bean (Vigna reflexo-pilosa Hayata), mungbean (Vigna radiata (L.) Wilczek), moth bean (Vigna aconitifolia (Jacq.) Maréchal), rice bean (Vigna umbellata (Thunb.) Ohwi and Ohashi), cowpea (Vigna unguiculata (L.) Walps.), Bambara groundnut (Vigna subterranea Verdc.) and zombi pea (Vigna vexillata (L.) A. Rich) [1,2]. The former six species are of Asian origin and belong to the same subgenus Ceratotropis, while the latter three species are of African origin. Cowpea and Bambara groundnut belong to the subgenus Vigna, while zombi pea belongs to the subgenus Plectrotropis. These crops are generally cultivated by resource-poor farmers as a sole crop or a component in various cropping systems.
Among the domesticated Vigna species, zombi pea is the least known crop. Two forms of zombi pea exists; seed type and tuber (storage root) types. The seed type is believed to be domesticated in Sudan (Africa) [3,4], while the tuber type is believed to be domesticated in Bali and Timor, Indonesia (Asia) [5] and India [6,7]. Genetic diversity analysis in a large set of V. vexillata germplasm using simple sequence repeat (SSR; also known as microsatellite) markers revealed clear difference between the two types and suggested that these two types were domesticated from different wild gene pools of V. vexillata [4]. The seed type zombi pea is classified as variety macrosperma. Possibly, this type is domesticated from wild zombi pea of East Africa [4]. It has erect growth habit and large seed size, and is insensitive to day length and no seed dormancy with some degree of pod dehiscence. The variety macrosperma is highly crosscompatible with wild V. vexillata [8]. The tuber type zombi pea has not yet been classified as a variety of V. vexillata. It has viny growth habit and large seed size, and is sensitive to day length and no seed dormancy with some degree of pod shattering. The tuber type is grown principally for fresh edible tuber roots which contain as high as 15% of proteins [5,6]. This cultivated type has been reported to be cross-incompatible with wild V. vexillata and the variety macrosperma [8]. The phenotypic difference between the two cultivated types of V. vexillata is a result of divergent selection during the evolution of the crop. Morphological and physiological differences that occur during the change from wild plants to cultivated crops is known as "domestication syndrome" [9]. Crop domestication is an accelerated evolutionary process that is the result from human intentional and unintentional selection and natural selections [10]. Domestication-related traits in legume crops include plant architecture (determinate growth habit), gigantism in the consumed plant organs (seed size and/or pod size), reduced seed dispersal (indehiscent pod), loss of seed dormancy and no response to photoperiod [11][12][13].
Knowledge on genetic basis of crop domestication is useful for identification of beneficial gene(s) in the wild relatives of the cultivated crops that can be used in plant breeding. Moreover, such knowledge can also provide insight into crop evolution and agriculture, such that the case of genetic analysis of indehiscent pod in soybean by Funastuki [14] which showed crucial role of the trait in the global expansion of soybean cultivation. Generally, genes controlling crop domestication are identified by quantitative trait loci (QTL) analysis. Among crops of the genus Vigna, QTL mappings of domestication syndrome traits have been carried out in azuki bean [13,15], mungbean [16], rice bean [17], cowpea/yardlong bean [18][19][20][21]. The genetics of domestication of zombi pea is of particular interest because two cultivated forms of this legume have experienced different domestications from wild zombi pea. Seed type zombi pea was domesticated from wild zombi pea of Africa, while tuber type zombi pea was domesticated from wild zombi pea of Asia [4,22]. Information on the genetics of domestication-related traits would be useful for zombi pea improvement programs, and for comparative genome studies among members of the genus Vigna. The objectives of this study were to identify QTLs for domestication-related traits in tuber type zombi pea and to compare them with previously reported QTLs for domestication in other Vigna species.

Materials and methods
An F 2 population of 139 plants was used in this study. The population was developed from a cross between accession JP235863 and accession AusTRCF66514 (Fig 1). JP235863 is a cultivated zombi pea collected from Bali Island, Indonesia which it is cultivated principally for edible tuber root, while AusTRCF66514 is a wild zombi pea from India. AusTRCF66514 was crossed with JP235863 as female and male parents, respectively, to produce F 1 hybrid. Only one F 1 seed was obtained and was self-pollinated to generate F 2 population. One hundred and eighty seven F 2 seeds together with parental seeds were sown in pots (1 seed per 1 pot) during June to November 2015 under a greenhouse condition at National Agriculture and Food Research Organization, Tsukuba, Japan. One hundred and fifty eight seeds germinated and developed into plants, but 139 of them grew vigorously and set flowers, and were used for DNA analysis and phenotypic data collection.

DNA extraction
Total genomic DNA was extracted from young leaves of the parental and the F 2 plants using a CTAB method described by Lodhi et al. [23]. DNA quality and quantity was checked by comparing with a known concentration ʎ-DNA (50 and 100 ng/μl) on 1% agarose gel. DNA concentration of each plant was adjusted to 5 ng/μl for PCR analysis. The adjusted DNA was also checked for quality using NanoDrop 8000 (Thermo Fisher Scientific, Wilmington, U.S.A.).

Trait measurements
Twenty-two traits related to domestication in Vigna crops [13,[15][16][17]20] were measured/evaluated (Table 1). Among those traits, seed color (SDC) and presence of black mottles (SDCBM) on seeds were considered as qualitative trait, while the others were considered as quantitative trait. Seedling traits, primary leaf width (LFPW) and primary leaf length (LFPL) were recorded when the first trifoliate was appeared. Stem thickness (STT), 10 th node length (STL10), maximum leaf length (LFML) and maximum leaf width (LFMW) were collected when eleventh trifoliate leaf was fully developed. Seed-related traits including number of seeds per pod (SDNPPD) was an average from ten pods, seed width (SDW), seed length (SDL) and seed thickness (SDT) was the average from ten seeds while 100-seed weight (SD100WT) was recorded from 100 seeds of each plant. The number of days from planting to first flowering (FLD) was recorded. Number of days to first mature pod (PDDM) chosen from first flower that had been developed to mature pod. After harvesting all pods, ten pods from each plant were used to measure/recorded for pod width (PDW), pod length (PDL), pod dehiscence (PDT) or number of twist along the length of pod (NTWP). PDT was recorded after the pods were kept in hot air oven at 40˚C for 24 hours. Tuber traits (NTB, TBWT, TBW and TBL) were collected at 260 days after planting. Broad-sense heritability of each traits was calculated by using the following equation; h 2 ¼ ðs 2 F2 À ðs 2 P1 À s 2 P2 Þ=2Þ=s 2 P , where h 2 is broad-sense heritability, s 2 F2 is the variance of F 2 population, s 2 P1 is the variance of JP235863 and s 2 P2 is the variance of AusTRCF66514.

RAD-seq analysis
RAD-seq analysis was conducted following a modified protocol described by Peterson et al. [35]. Genomic DNA of parents and each F 2 plant was digested by double-digest RAD library method. In brief, 10 ng of genomic DNA was digested with EcoRI and BgIII (New England Biolabs, Ipswichs, MA, USA), and purified. Each adaptor with unique 4-8 bp index was ligated to digested DNA samples. The adaptor sequence was as follows: TruSeq_EcoRI_adaptor

Extraction of RAD-tag and bi-allelic RAD-marker detection
RAD-tag sequences were extracted and bi-allelic RAD-markers were detected using Stacks ver. 1.30 [36] following the procedures described by Marubodee et al. [33]. With the Stacks pipeline, the low quality sequences were filtered to obtain 51 bp RAD-tag sequence reads, and classified the RAD-tags sequences into the parental accessions and F 2 individuals. Then ran denovo_map.pl program where RAD-tags with less than two sequence mismatches were grouped as a stack, parental stacks with less than three mismatches were estimated as those derived from homologous loci and lists of RAD-tag sequences and their count were constructed for each sample. We called genotypes of F 2 individuals only when stacks have more than five RAD-tags (minimum depth of 5).

Linkage map construction using SSR ad RAD-seq markers
A genetic linkage map of the F 2 population was constructed using software Join Map 4.0 [37]. Polymorphic markers showing more than 10% missing genotypic data were excluded from linkage analysis. Chi-square analysis was used to test segregation of the marker loci (1:2:1 for co-dominant markers and 3:1 for dominant markers). The markers were grouped using minimum logarithm of the odds (LOD) of 3 and recombination frequency of 0.40. Genetic distance between markers in centimorgan (cM) unit was calculated using Kosambi's mapping function [38].

QTL analysis
QTLs controlling domestication-related traits were located onto the linkage map by inclusive composite interval mapping (ICIM) method [39] using the software QTL IciMapping 4.1 [40]. Probability value for entering variables in stepwise regression of phenotype on marker variables (PIN) was 0.001. ICIM was performed at every 1 cM. Significant log of odds (LOD) threshold for QTL analysis of each trait was determined by 3,000 permutation test at P = 0.05.

RAD-seq markers
Illumina sequencing with HiSeq 2000 yielded a total of 29,568,103,697 RAD-tag 51-base reads from 578,901,731 raw reads. The number of RAD-tags of parents (AusTRCF66514 and JP235863) and F 2 individuals was 3,124,320, 2,700,159 and 573,077,252, respectively. The average number of RAD-tags per F 2 individual was 4,122,857. RAD-tags were aligned and clustered into 6,576 stacks. Average fill gaps rated for this 139 F 2 individuals was 0.747 (max = 0.949, min = 0). Among these stacks, 479 RAD markers showed homozygote polymorphic genotype between parents. However, among these 479 RAD markers, 362 RAD markers were discarded because they had percentage of missing data more than 10%.

Linkage map analysis
Finally, 264 markers (145 SSR, 117 RAD and 2 morphological markers) were used for the linkage map construction after discarding markers showing identical F 2 genotypes, missing data more than 10%. Among these markers, as high as 128 (48.5%) showed segregation distortion. The 264 markers were clustered into 11 linkage groups (LGs 1-11) ( Table 2 and Fig  2). The map spanned 704.8 cM in total, with a mean distance between markers of 2.87 cM. The lengths of linkage groups ranged from 30.9 cM (LG10) to 112.7 cM (LG2). A gap between markers more than 15 cM presented on only LG11 (between markers DMBSSR192 and CEDG066).
LGs 5 and 8 possessed many markers showing severe segregation distortion (P< 0.001). All the distorted markers on LG5 showed excessive cultivated parent genotype. Majority of the distorted markers on LGs 6 and 11 and about half of the distorted markers on LG10 showed excessive heterozygous genotype. All the distorted markers on LGs 7 and 8 and about half of the distorted markers on LG10 showed excessive homozygous wild genotype.
Primers CEDG065 and CEDG244 each amplified two different loci (markers). One of the loci from CEDG065 (named as CEDG065a) and one of the loci from CEDG244 (named as CEDG244a) were mapped near each other on LG2. The other loci from CEDG065 (named as CEDG065b) and the other loci from CEDG244 (named as CEDG244b) were mapped adjacent to each other on LG7.This indicated partial genome duplication and orthologous block of LG2 to LG7 or vice versa.

Comparative linkage analysis
The linkage map of zombi pea developed in this study were compared with previous linkage maps developed for Vigna crops including yardlong bean [19], azuki bean [41], mungbean [16] and rice bean [17] by using common SSR markers (Fig 3). The comparison revealed that in general the linkage groups and orders of the markers among the Vigna species were the same or highly similar. Nonetheless, remarkable macro genome rearrangements were found on LG5 and LG7 of our linkage map for V. vexillata. About a half of LG5 (middle to bottom) appeared to be an orthologous blocks from LG4, while a large portion of LG7 appeared to be orthologous blocks from LG4 and LG10. In addition, a small portion of the LG7 appeared to be duplication from LG2, and vice versa. The length of the zombi pea linkage map developed in this study was shorter that the linkage maps of others Vigna species.

Variation of domestication-related traits
The parents were clearly different in all the traits recorded. The cultivated zombi pea showed higher value of trait value than the wild zombi pea in all the traits except pod shattering, seeds per pod, tubers per plant, and tuber length. For qualitative traits, the cultivated parent had yellow seed, while the wild zombi pea had brown seed with black mottle (Fig 1). The mean, range and standard deviation, and broad-sense heritability of the quantitative traits are shown in Table 3. Mean of the traits of the F 2 population was between mean of parents for all traits except pod length, tuber length, 10 th node length and number of seeds per pod. The measured traits in the F 2 population showed nearly a normal distribution (Table 3). PDT, SD100WT, SDL, SDW, SDT, PDW, STT, LFPW, LFML, LFMW, TBWT, STL10 and PDDM showed moderate to high heritability (>60%), whereas LFPL, TBW, TBL and FLD showed medium to low heritability (<60%). PDL, NTB and SDNPPD showed heritability lower than 30%.
There was significant and positive correlation (P<0.05) between related traits in F 2 population, such as between seed-related traits (SD100WT, SDL, SDW and SDT), between pod width and seed-related traits (SD100WT, SDL, SDW and SDT), between SDNPPD and PDL, between LFPW and LFMW, between TBW and TBWT (S2 Table). Negative correlation was found between TBW, TBWT and SDNPPD, STL10 and LFML and NTB, TBL and FLD (S2 Table). In general, correlation between related traits was moderate to high (>0.50), while correlation between non-related traits was low or none.

QTLs for domestication-related traits
The results of the QTL analysis for each trait in F 2 population are shown in Table 4. Out of 22 traits, ICIM did not find any QTL for four traits. Among those 18 traits which QTLs were identified, six traits each had only one QTL.
Pod dehiscence (PDT). Both cultivated and wild parents showed small difference in number of twists along the pod (Table 3). In the F 2 population, there was little variation of the trait. As a result, no QTL was detected for this trait.
Increase in organ size. Seed size, pod size, stem thickness, leaf size and tuber-related traits were considered.
Seed size (SD100WT, SDL, SDW and SDT). One to three QTLs for traits related to seed size were located on LGs 2, 5, 7 and 8. At all QTLs detected, alleles from the cultivated parent increase the seed size. The QTLs with the largest phenotypic effects for all four traits were located on LG7 (28.19%, 32.60% and 24.02% for seed weight, seed length and seed width, respectively).
Pod size (PDL and PDW). The cultivated zombi pea had larger pod size than the wild zombi pea. Three QTLs for each PDL and PDW traits were detected on LGs 1, 2, 4, 5 and 7. The QTLs with the highest contributions to these traits (PVE = 20.65% and 18.11%) were found on LGs 4 and 2. QTLs for both PDL and PDW on LGs 5 and 7 were located close to QTLs for seed size traits (Fig 4).
Stem thickness (STT). Only one minor QTL (PVE = 11.25%) was detected on LG7 for stem thickness which was influenced by alleles from the cultivated parent that has thicker stem than the wild zombi pea.  Leaf size (LFPL, LFPW, LFML and LFMW). Only one QTL was detected on LG8 for LFPL. Three QTLs for LFPW were detected on LGs 1, 3 and 5. One QTL was detected on LG5 for LFML. Two QTLs for LFMW were detected on LGs 3 and 5. The QTLs on LG5 for LFPW, LFML and LFMW showed highest PVE and were located very close to each other and likely to be the same locus. Similarly, the QTLs located on LG3 for LFPW and LFMW were located not far from each other. Alleles from the cultivated zombi pea increased leaf width, while alleles from the wild zombi pea increased leaf length.
Plant type. Stem length was considered. Stem length (STL10). No QTL identified for this trait, although the trait showed high heritability (83.4%).
Earliness. Days to first flowering and days to first maturity were considered. Days to first flowering (FLD). The cultivated and the wild parents were very different in FLD (102.6 and 72.4 days, respectively). The F 2 population showed a high variation for the trait (Table 3). However, only one QTL was detected for this trait. The QTL was on LG1 and accounted for only 17.37% of the trait variation. The alleles from the cultivated parent prolong FLD.
Yield potential. Or Number of seed per pod (SDDNPPD). The cultivated parent produced lower SDDNPPD than the wild parent. Five QTLs for SDDNPPD were detected LGs 5, 7, 8 and 9. Two QTLs on LG7 had high effects (PVE = 25.14% and 11.89%). At QTLs Sddnppd5.1-and Sddnppd7.2-, the allele from the cultivated parents increased SDDNPPD, while at the other three QTLs the alleles from the wild parent increased SDDNPPD.
Pigmentation. Seed coat color (SDC) and black mottle on seed coat (SDCBM) Seed coat color (SDC). SDC was characterized as a qualitative trait. The cultivated parent had yellow seed coat but wild parent had brown seed coat color. F 2 plants segregated for this trait at ratio of 81 (brown) to 24 (yellow), fitting a 3:1 ratio (χ 2 = 0.26, P = 0.61). This indicated that SDC is controlled by a single locus, named as Sdc. This locus was mapped as a morphological marker onto LG4 (Fig 2).
Black mottle on seed coat (SDCBM). SDCBM was also characterized as a qualitative trait. The cultivated Bali had no black mottle, while the wild parent had black mottle on seed coat. F 2 progenies segregated for seed coat color at ratio of 76 (black mottle presence) to 29 (black mottle absence), fitting a 3:1 ratio (χ 2 = 0.38, P = 0.54). This indicated that SDCBM is controlled by a single locus. The resulted indicated that the presence of SDCBM is controlled by a single dominant locus, named as Sdcbm. Sdcbm was mapped next to the locus Sdc at marker interval VR045 and C_13986_405837 (Fig 2).
LG4. Large-effect QTLs for pod length (Pdl4.2+) was on this linkage group. In addition, seed coat color (SDC) and black mottle (SDCBM) were mapped as morphological makers onto this linkage group.

Comparison of QTLs for domestication in V. vexillata with those of other Vigna species
The QTLs for domestication-related traits for zombi pea identified in this study were compared with those for yardlong bean [20], azuki bean [13], mungbean [16], and rice bean [17] based on common SSR markers (Fig 4). The comparison showed that most of the QTLs detected for zombi pea were also found in other Vigna species; they were mapped to the same linkage groups. For example, most of the QTLs for seed-related traits with high PVE in zombi pea were mapped onto LG7 in which several QTLs for seed-related traits were detected for yardlong bean, azuki bean, mungbean and rice bean. Three QTLs, Pddm2.1-, Sdt8.1+ and Sdnppd9.1-were mapped to similar regions with other Vigna species. Nonetheless, none or a few of QTLs for seed-related traits were detected on LG6 and LG10 of zombi pea which is the same case as in other Vigna species except yardlong bean.

Fertility of progenies of cultivated zombi pea from Bali x wild zombi pea
A previous study on reproductive compatibility among cultivated zombi pea from Bali (tuber type zombi pea),cultivated zombi pea from Africa (seed type zombi pea), and wild zombi pea from Africa and Australia revealed high compatibility between the African cultivated zombi pea (seed type zombi pea) and the wild form from both Africa and Australia, but revealed various pre-and post-zygotic incompatibility between the Bali cultivated zombi pea (tuber type zombi pea) and the wild form from Africa and Australia [8]. In our present study, we successfully obtained a fertile F 2 population for the first time from a cross between the cultivated zombi pea from Bali and Indian wild zombi pea using the former as female parent. This suggested that the Bali cultivated zombi pea constitutes a primary gene pool with the Indian wild zombi pea. The difference between result in our study and that in Damayanti et al. [8] is possibly due to environmental factors. Damayanti et al. [8] noted that even self-pollination of the cultivated zombi pea from Bali results in low pod setting. In our study, when we conducted the hybridization during November 2014 to February 2015, the Bali cultivated zombi pea showed very low pod setting (<3%). When it was grown again during November 2015 to February 2016 and during November 2016 to February 2017, it set more flowers and pods (~15% and 70%, respectively; Somta and Dachapak, personal observation). In addition to environmental factors, genetic relatedness of the parents appeared to affect the success of hybridization. The wild zombi pea (JP235863 from India) and the Bali cultivated zombi pea were genetically closely related [4], while the wild zombi pea used in the study of Damayanti et al. [8] were from Africa and Australia and they were genetically highly differentiated from the Bali cultivated zombi pea [4].

Transferability and polymorphism of SSR markers
Although as high as 1,876 SSR primers pairs from various Vigna species were screened for polymorphism between the cultivated and wild parents, only 36.6% of them were amplifiable and only 10.7% of them were polymorphic. The low polymorphism of SSR markers was also observed in the mapping parents of V. vexillata used by Marubodee et al. [33] in which only 6.2% of 1,336 SSR markers were polymorphic. In addition, Dachapak et al. [4] reported that only 21.2% of 1,024 SSR markers screened in 6 accessions of V. vexillata from Asia, Africa, Australia and America were polymorphic. However, the percentage of amplifiable SSR markers in Marubodee et al. [33] and Dachapak et al. [4] (65.4% and 58.1%, respectively) was higher than in our study (36.6%). It is worth noting that almost all of the markers used by Marubodee et al. [33] and Dachapak et al. [4] were screened in our study. Nonetheless, these results indicate moderate transferability but low polymorphism of SSR markers from other Vigna species in V. vexillata. High transferability rate of SSRs (>80%) among several other Vigna species have been reported earlier [19,25,29,34,42]. The low transferability rate of SSRs from other Vigna species to V. vexillata is likely due to the fact that those Vigna species and V. vexillata are genetically highly differentiated [43]. Due to the low transferability and polymorphism of SSR markers in V. vexillata, other types of markers (RAD-seq in this study) should be used for genome analysis of V. vexillata.

Linkage map construction and orthologous blocks in zombi pea
Previously, there was only two genetic linkage maps developed for zombi pea [33,44]. The first map contained 14 linkage groups with all were dominant markers (70 RAPD and 47 AFLP) except one co-dominant SSR marker. The map reported by Marubodee et al. [33] comprised 11 linkage groups from 84 SSR and 475 RAD-seq markers. The number of linkage map of zombi pea constructed in our study comprised 11 linkage groups of 145 SSR markers and 117 RADseq markers. The number of linkage groups of the maps developed by Marubodee et al. [33] and by this present study corresponded with the chromosome number of V. vexillata. (n = 11).
Previous comparative genome studies in the genus Vigna by means of comparison of linkage maps revealed high genome conservation among the species [17,19,25,42,45]. We compared our V. vexillata linkage map with other Vigna linkage maps using common SSR markers (Fig 3). We found that middle to lower part LG5 of our V. vexillata linkage map was very likely an orthologous block from LG4 and that LG7 of our V. vexillata linkage map was orthologous blocks from LG4 and LG10. In addition, a part of the LG7 appeared to be duplication of LG2 or vice versa (Fig 3). Intraspecific macro-orthologous block rearrangement has been demonstrated in azuki bean by means of comparative linkage mapping and fluorescence in situ hybridization [46,47]. The rearrangement of orthologous block was reciprocal where LG4 and LG6 were intermingled and was found in many accessions of wild azuki bean and it possibly affected fitness of those wild accessions to some environments and seed size [47].

Segregation distortion
Segregation distortion is a normal phenomenon in interspecific hybridization and in even intraspecific hybridization between two highly genetically different genotypes such as between wild and cultivated form of the same species. Segregation distortion of DNA markers has been reported for intraspecific hybridization of several Vigna species including V. radiata [16], V. angularis [41,46], V. umbellata [17], V. mungo [42] and V. unguiculata [19,20]. The percentage of marker distortion in these reports ranged from 3.9% in V. angularis to 48.7% in V. unguiculata. In this study, as high as 48.5% of the markers showed segregation distortion, the value is considered very high as compared to that other reports in Vigna species. The distorted markers were clustered on LGs 5, 6, 7, 8, 10 and 11 (Fig 2). The presence of a gene(s) controlling sterility and/or compatibility may cause segregation distortion of nearby loci and clustering of distorted markers [48]. The highly distorted markers (P< 0.001) were on LGs 5, 7 and 8. The distortions found on LG5 and LG7 possibly stem from the chromosomal rearrangement in these two linkage groups (Fig 2; see also above discussion on orthologous block). The cause of the distortion on LG8 is not known, however, it is possible that this linkage group represents the most genome divergence between the wild and cultivated parents used in this study. It is worth noting that major QTLs for tuber-related traits (Tbw8.1-and Tbwt8.1-) were located on the LG8 (Fig 4 and Table 4). Tuber root trait in the cultivated V. vexillata from Bali used in this study is the most strikingly different adaptive/domestication trait distinguishing the cultivated V. vexillata Bali from other forms of V. vexillata. In yardlong bean (V. unguiculata), major QTLs for pod length and other pod-related traits (the most important domestication trait for this crop) and seed traits were found on LG7 where highly segregation distortions existed [19]. Similarly, in mungbean, major QTLs for seed dormancy, seed productivity and day length sensitivity were located on LG4 that showed high level of segregation distortion [16]. Most of markers mapped on LG11 of intra-and inter-specific crosses of Vigna species always showed segregation distortion [17,19,45]. This is also the case in our study where 58.8% of the markers mapped to LG11 showed segregation distortion (Fig 2). Therefore, LG11 appeared to play an important role in genetic differentiation within and among Vigna species.
Although pod length (PDL) and number of seed per pod (SDNPPD) were highly correlated (r = 0.83) (S2 Table), only one of the QTLs for these two traits were co-located (Table 4 and Fig 4). Since three of the five QTLs, including the largest effect QTL, detected for SDNPPD were mapped to genome regions that showed high segregation distortion, those QTLs may be associated with genetic compatibility and fertility, and thus seed production. In Vigna species, early generation (F 1 and F 2 ) hybrid progenies derived from hybridization of distantly-related genotypes always possesses pods with incomplete filled (empty seeds in some locules) [49,50]

QTLs for tuber traits in zombi pea
Most of tubers crops are domesticated from non-legume species such as potato, sweet potato, yam and cassava. In Vigna species, V. vexillata, V. lobatifolia and V. marina produced tuber roots [51]. Tubers produced by legume crops have much more nutrition than non-legume crops especially nitrogen or protein due to nitrogen fixation ability. The protein content in tuber roots of V. vexillata is about 15% which is five-fold higher than that in sweet potato. In potato, tuber weight was controlled by a few QTLs with PVE of 20% or less [52]. Similarly, only three QTLs were identified for tuber weight (TBWT) in V. vexillata with largest PVE of about 15%. Tuber weight of potato has been shown to be negatively correlated with leaf area (leaf size) [53]. On the contrary, cultivated alleles that increase tuber width (Tbw3.1-), primary leaf width (Lfpw3.1-) and maximum leaflet width (Lfmw3.1-) were detected at similar position on LG3, which suggested preiotropic effect of the same gene (Table 4 and Fig 4). It should be noted that alleles from wild V. vexillata increased the tuber weight (Tbwt2.1+, Tbwt4.1+) and tuber width (Tbw2.1+, Tbw4.1+), in spite of the fact that wild V. vexillata has smaller tubers ( Table 4 and Fig 1). Map positions of these QTLs are close and therefore tuber weight and tuber width increase by wild allele could be preiotropic effect of the same gene. These reverse effect alleles of wild V. vexillata may account for transgressive segregation found in the F 2 population (S1 Fig). In fact, there is a F 2 plant (F 2 -18) with 141-160 g range tuber weight, while tuber weight of cultivated parent is 61-80 g range (S1 Fig). Although QTL could not be detected, tuber number increase by wild Indian parent also contribute tuber weight increase. Therefore, these alleles from the wild V. vexillata will be useful for improvement of tuber yield of the Bali cultivated V. vexillata.

Genomic region and distribution of QTLs for other domestication traits
Many domestication-related traits were studied in Vigna crops including mungbean [16], yardlong bean [19,20], azuki bean [13] and rice bean [17]. In general, these studies found that (i) domestication-related traits were controlled by one or two major QTLs together with some minor QTLs in a narrow genome region, (ii) major QTLs for different traits were clustered and (iii) QTL clusters were not randomly distributed along linkage groups. For examples, major QTLs for seed size, pod size, and seed dormancy were clustered on LG4 of mungbean, and major QTLs for pod length, seed size, pod shattering, leaf size, and stem thickness were clustered on LG7 of yardlong bean. In this study, we did not find clusters of major QTLs for different traits in zombi pea (Fig 4 and Table 4).

Comparison of domestication QTLs between V. vexillata and yardlong bean
Since V. vexillata and is genetically closely related with cowpea/yardlong bean (V. unguiculata) and both of them originate from African, we compared domestication-related QTLs for V. vexillata found in this study with those for yardlong bean reported by [19,20]. Only QTLs having phenotypic effect more than 20% were considered ( Table 5).
The cultivated zombi pea used in this study and yardlong bean are believed to have been domesticated in Asia, although their origin of species is believed to be in Africa [4,19]. Major QTLs of yardlong bean were detected on four linkage groups (LGs 1, 3, 7 and 11), while largeeffect QTLs of zombi pea were found on three linkage groups (LGs 4, 5 and 7). Although some major QTLs were located on LG7 of both species, distributions of QTLs on this linkage group were different.
Seed weight. The 100-seed weight of the cultivated zombi pea was 7.3 g, while the wild zombi pea was only 2.1 g. Two QTLs for this trait were located on LG2 and LG7 (the largest effect (PVE = 28.19%) was on LG7). In yardlong bean, ten QTLs of 100-seed weight were detected on nine linkage groups with the highest QTL effect on LG7 (24.6%). Although the majors QTL for seed weight of zombi pea and yardlong bean were both located on the LG7, QTLs of domestication syndrome in zombi pea they appeared to be different because the QTL region in zombi pea was an orthologous block from LG10 as compared to yardlong bean (Fig 3). Pod size. Pod size of the cultivated zombi pea is not much longer and thicker than that of the wild zombi pea (94.3 vs. 85.6 mm. for length and 5.7 vs. 3.6 mm. for width, respectively). Six QTLs for pod size in zombi pea was found on LGs 1, 2, 4, 5 and 7 with the largest effect QTL (PVE = 20.7%) on LG4. Seventeen QTLs for pod size in yardlong bean was found on every LGs except LG10 with the largest effect QTL on LG7 for both PDL and PDW (PVE = 31.0% and 31.2%, respectively). One of the QTL for pod size on LG4 might be shared between zombi pea and yardlong bean (Fig 4) was found between markers CEDC055 and CEDG185.
Yield potential. Number of seeds per pod (SDNPPD) of the cultivated zombi pea was lower than that of the wild zombi pea (7.6 and 12.7, respectively). Five QTLs on LG5, LG7, LG8 and LG9 were detected with QTL showing highest PVE on LG7 (25.14%). In yardlong bean, two QTLs were detected on LGs 7 and 11. The QTL on LG11 showed highest PVE (70.1%).
Seed size and weight. Seven QTLs of seed size and weight were detected for zombi pea in which QTLs with high effect were on LG7 (Sd100wt7.1-, Sdl7.1-and Sdw7.1-). In azuki bean, rice bean and mungbean, large-effect QTLs for seed size and weight were co-located on the same region; LGs 1, 2, 3 and 9 for azuki bean, LG4 for rice bean, and LG8 for mungbean. This suggested that major genes controlling seed size with different biological mechanisms exist among domesticated Vigna species.
Pod size. Among the QTLs for pod length and pod width in zombi pea, the QTL Pdl4.2+ on LG4 showed the highest PVE (20.65%). QTL with the highest effect (PVE = 23.1%) for pod length in rice bean was also detected on LG4, while the QTL with the highest effect for pod length in azuki bean and mungbean were on LG7.
Leaf size. Among seven QTLs for leaf size in zombi pea, two QTLs (Lfpw5.1-and Lfmw5.1-) showed PVE higher than 20%; both of which were on LG5. In rice bean, a major QTL for leaf size was located on LG4. No QTL was detected for leaf size in azuki bean and mungbean due to small difference in the mapping parents.
Yield potential. Five QTLs for SDNPPD in zombi pea were located on LGs 5, 7, 8 and 9 with the QTL on LG7 showed highest PVE (25.14%), while three QTLs for SDNPPD in rice bean were detected on LG4 and LG9, and two QTLs for SDNPPD in mungbean were identified on LG1 and LG9. The QTLs on LG7 of zombi pea, rice bean, and mungbean were identified in a similar location and may be the same locus. However, all the QTLs identified in rice bean and mungbean possessed PVE lower than 20%.
Pigmentation. Seed coat color and presence of black mottle on seed coat were treated as morphological markers which were mapped next to each other on the LG4. In azuki bean and mungbean, gene controlling the presence of black mottle on seed coat was also mapped on LG4. This indicated that the gene controlling this trait is highly conserved among species in the genus Vigna.

Conclusions
An F 2 of V. vexillata population derived from a partially fertile hybrid between tuber-root-type cultivated form and wild type form were successfully developed. A genetic linkage map was constructed for this population utilizing 145 SSR, 117 RAD-seq and 2 morphological markers and to locate QTL for domestication syndrome traits. In total, 37 QTLs were detected for 18 domestication traits. Eight QTLs with large effect (>20%) were located on 4 out of 11 linkage groups. Domestication traits including seed size, pod size, leaf size, yield potential and seed pigmentation were controlled by one or two major QTLs and a few minor QTLs. QTLs for tubers were found on five different linkage group. Zombi pea had highest number of shared common QTLs with azuki bean, followed by yardlong bean, rice bean and mungbean. QTLs for seed size (SD100WT, SDL, SDW and SDT) were clustered together on different linkage groups of each species, suggesting specific and different seed size increase genes were used independently for each species. This study provides a genetic map and information for marker-assisted selection in zombi pea breeding. QTL analysis revealed potential of wild genetic resources to improve tuber size and yield. The comparative genomic analysis provides more understanding on genome evolution of Vigna species.