The mapping population consisted of 168 F_2:3 progenies from a cross between two cowpea subspecies accessions: (V. unguiculata ssp. unguiculata) cultivar “Green pod cowpea” (female) and yardlong bean (V. unguiculata ssp. sesquipedalis) cultivar “Xiabao II” (male) (Figure 1, Supplementary Figure 1). The two parents showed distinct morphological traits, including seed coat color, pod length, carpopodium length, grain number per pod and grain size.

Field trials were performed in the Caidian district (114.03E, 30.58N), Wuhan city, Hubei province, P.R. China from July to October 2014. Each sample was planted 1 m apart in rows 1.5 m in length.

The authors declare that the experiments comply with the current laws of the country in which they were performed.

For each F_2:3 progeny and the two parents, trait evaluation included pod length (PL, cm), carpopodium length (CL, cm), grain number per pod (GN), and thousand-grain weight (TGW, g). PL was measured as the average length from the top to the end of ten pods. CL was measured as the average length from the upper part to the end of carpopodium of ten pods. GN was measured as the average number of seeds of 10 uniform capsules from each plant. To measure the TGW, the mean weight of three independent samples of fifty grains was calculated first, and then was converted to TGW. All samples of F_2:3 population were tested for phenotyping at the harvest stage.

Young leaves were collected from each sample and then frozen at −80°C. Total genomic DNA was extracted from each sample using the modified CTAB method as described by Doyle and Doyle (1987). Genomic DNA from each sample was quantified using spectrophotometry (Qubit 2.0 Fluorometer, Invitrogen) and checked for genomic integrity by 0.8% agarose gel electrophoresis.

Restriction-site associated DNA (RAD) strategy combined with Illumina DNA sequencing was used for the fast and effective identification of SNP markers. RAD libraries were constructed following a protocol from Baird et al. (2008). The general protocol was as follows: first, each sample of genomic DNA (~0.5 μg per sample) from the 168 individuals and their two parents was digested with 20 units (U) of EcoRI or Sbf I (New England Biolabs, NEB) for 15 min at 37°C in a 50-μL (microlitre) reaction and then heat denatured at 65°C. P1 Adapters, each with a unique 4–8 bp (base pair) molecular-identifying sequence (MID), were ligated to the genomic DNA of each sample. The adapter sequences were as follows: EcoRI digestion, top: 5′-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCTxxxxx-3′ [x = barcode], bottom: 5′-Phos-AATTxxxxxAGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT-3′, for Sbf I digestion, top: 5′-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCTxxxxTGCA-3′, bottom: 5′-Phos-xxxxAGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT-3′. The ligation reaction was incubated at room temperature (RT) for 20 min in a final volume of 10 μL containing 2.5 μL P1Adapter (100 nM), 1 μL rATP (100 mM) (Promega), 1 μL 10 × EcoRI buffer, 0.5 μL T4 DNA Ligase (1,000 U) (NEB), and 5 μL H₂O. Samples were treated with heat denaturation at 65°C for 20 min again. The pooled samples of each library were randomly cut to DNA fragments of 400–700 bp. The size of sheared DNA was isolated through 1% agarose gel electrophoresis using a DNA gel extraction kit (Tiangen, China). The ends of the DNA were blunted by using the Quick Blunting Kit (NEB). After purifying the samples using a DNA purification kit (Tiangen), 15 U of Klenow exo⁻ (NEB) was added for adenine (Fermentas) overhangs on the 3′end of the DNA at 37°C. Then, 1 μL P2 Adapter (10 nM) was ligated to the DNA fragments at RT. Samples were again purified and eluted in 50 μL. Five microlitres of this product was used in a PCR amplification along with 50 μL Phusion Master Mix (NEB), 40 μL H₂O and 5 μL primer mix (10 μM) (P1-forward primer: 5′-AATGATACGGCGACCACCGA-3′; P2-reverse primer: 5′-CAAGCAGAAGACGGCATACGA-3′). The PCR enrichment was carried out following instructions (NEB) for 11 cycles in total. Samples were gel purified, excising DNA 400–700 bp, and diluted to 10 nM.

Each library of the RAD products from the 170 individuals was sequenced on an individual lane of the Illumina HiSeq2000 next-generation sequencing platform (BGI-Shenzhen, Shenzhen, P.R. China). Sequencing data for each individual were then obtained in terms of the specific MID. Raw sequences are available in the National Centre for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database under the accession number SRR5097134.

Illumina raw sequence reads were filtered out at first if lacking expected restriction enzyme or ambiguous barcodes. Thus, the raw reads were ensured to more than 90% of the nucleotides having a quality value above Q30 (equals 0.1% sequencing error) and more than 95% above Q20 (equals 1% sequencing error). SNPs were discovered using a two-step process. In step 1, four individuals having the maximum sequencing reads were used to generate “reference tags” using Ustacks and Cstacks (Catchen et al., 2011). Then, the SNPs of the two parents and their offspring were identified by aligning the clean reads of the offspring to the reference tags. SNPs were regarded as true polymorphisms when each allele was represented by a minimum coverage of three reads or alleles were observed in at least five reads with a score >20 (p > 0.05). Genotypes of offspring individuals were determined by comparing to the parental genotypes. SNP markers with <10% missing data among the 170 individuals and meeting the above criteria were used in the following linkage mapping.

The construction of genetic linkage map was performed by using JoinMap4 (Van Ooijen, 2006) and MSTmap (Wu et al., 2008). Linkage groups (LGs) were first defined using a minimum LOD threshold of 15 and a maximum recombination of 45% with the software JoinMap4. The Kosambi mapping function in the software MSTmap was used to convert recombination frequencies into map distance for each linkage group. The final marker order of each linkage group was verified by the software program RECORD (Van et al., 2005). The final marker order was verified in each linkage group by using RECORD (Van et al., 2005). The linkage map was graphically visualized with MapChart2.2 (Voorrips, 2002). The corrected length of the linkage map was assessed by multiplying the length of each linkage group by (m + 1)/(m−1), where m is the number of markers in the linkage group (Chakravarti et al., 1991). The coverage of the genome by the linkage map was then estimated by calculating c = 1 −e^−2dn/L, where d is the average interval of markers, n is the number of markers, and L is the length of the linkage map (McDaniel et al., 2007).

To verify the map quality, we compared the new SNP map and recently reported cowpea maps (Muñoz-Amatriaín et al., 2017). We identified identical and proximal markers between the linkage maps using a reference genome of a related species as an intermediary (e.g., common bean) because the cowpea genome shares a high degree of collinearity with the common bean (Vasconcelos et al., 2015).

To assess the quality and universality of genetic maps, comparison of marker names, sequences, LGs and cM positions were carried out between our RAD-based SNP map and a recent published genetic map of cowpea (Muñoz-Amatriaín et al., 2017). Comparisons of linkage group composition and order were investigated for P. vulgaris. The basic workflow was as follows. First, all marker sequences were combined into a single FASTA file and mapped to a reference genome (e.g., common bean genome, NCBI accession number: ANNZ00000000.1) and published scaffolds (https://www.ncbi.nlm.nih.gov/genome/genomes/380). When two markers mapped to the same reference genome scaffold or contig, the two closest markers were considered as a marker pair. Each marker can only be paired once, and any other marker that was seconds closest (or further) to the now-paired marker is discarded. The paired markers were visualized by MapChart (Voorrips, 2002).

The mean phenotypic data from all the 170 individuals (two parents and 168 F_2:3 progenies) were calculated for frequency distributions using function formula “NORMDIST” in Excel 2010.

Quantitative trait loci (QTL) were detected for each of the traits using the CIM (composite interval mapping) method implemented in Windows QTL Cartographer 2.5 (Wang S. et al., 2012). The LOD significance thresholds (P < 0.05) were analyzed by running 1,000 permutation tests. QTL regions were visualized by the program Matplotlib (version 1.4.3) in software Python version 2.7.10 (Sanner, 1999).

To compare the genome structures of cowpea and a closely related species, a comparative analysis was carried out between the cowpea linkage map and the whole genome sequences of adzuki bean V. angularis. BLASTN was used to perform similarity searches (using the nucleotide sequences from which the SNP markers were mapped in the current linkage map of cowpea) against each pseudomolecule of the genome sequences of V. angularis, with a threshold e-value of 1e-20. The graphical comparative maps were drawn using the Circos program (Ubi et al., 2000).

In addition, sequences of cowpea RAD tags located in QTL regions were chosen to reveal potential genetic information by using BLASTN. These sequences of the mapped markers of cowpea RAD-tags were searched against genomes of legumes deposited in the NCBI database using BLASTN with an e-value cut-off of 1e-10 (www.ncbi.nlm.nih.gov/blast/).

A total of 170 RAD-seq libraries from two parents and 168 offspring were constructed and sequenced on an Illumina HiSeq2000 platform. After data filtering, 768,592,488 raw reads, consisting of ~69.315 Gb of sequencing data, were individually assigned to RAD tags.

A total of 383,077 RAD-tags were obtained. Using all the tags as reference sequences, SNP loci derived from paired end reads of 170 individuals were identified. Overall, 34, 868 polymorphic SNPs between the offspring genotypes and the two parent genotypes were identified on basis of strict SNP selection and criteria.

The set of 34,868 SNPs were used for genetic map construction. A total of 17,996 SNPs were successfully mapped onto 11 different LGs, covering 1,194.25 cM of the cowpea genome with an average distance of 0.066 cM between adjacent markers (Table 1 and Supplementary Figure 2). The length of the individual LGs ranged between 76.19 cM (LG11) and 198.85 cM (LG1), with an average inter-locus distance of 0.03 (LG2) to 0.13 cM (LG10 and LG11). LG1 was the densest, having 3,343 SNP loci, with an average density of 0.42 cM, whereas LG11 had the least number of SNP loci (591). On average, each LG consisted of 1,636 SNP loci spanning 108.57 cM. The corrected length of the linkage map was estimated at 1,195.97 cM, which was converted to genome coverage of 86.3%. Loci names and genetic distance of SNP positions on the 11 LGs of the genetic map are listed in Supplementary File 1.

To reveal the characterization of SNPs in cowpea, we investigated distribution and percentages of different SNP types. We found that the 17,996 SNP markers were genome-wide distributed across the 11 LGs, implying 1,636 SNPs per linkage group (Table 1). The maximum number of SNPs (3,343) was found in LG1, while the minimum number of SNPs (591) was found in LG11. Among the SNPs, the dominant type was transition, and the C/T and G/A types contributed to 34.8 and 35.8% of the SNPs, respectively. The other four SNP types were transversion including C/G, G/T, C/A, and A/T. Their proportion varied from 6.2 to 7.8% accounting for 29.4% of all SNPs (Supplementary Table 1). The SNP-flanking sequences and the polymorphic sites are summarized in Supplementary File 2.

To compare the maker position of our linkage map and the newly published consensus cowpea linkage maps (Muñoz-Amatriaín et al., 2017), we employed the partial reference genome available of Phaseolus vulgaris (Schmutz et al., 2014). To assign LGs to chromosomes, both our linkage map and the newly published consensus cowpea linkage maps (Muñoz-Amatriaín et al., 2017) were compared by mapping marker sequences against the common bean genome. Markers from both cowpea linkage maps were identified and aligned if they hit against the same contig/scaffold. Descriptive statistics for the identified markers on the LGs are in Table 2. Then, based on these markers, extensive synteny was observed using Circos (Krzywinski et al., 2009). A full comparison was completed across the two existing maps and the common bean genome. A set of 3,695 SNPs from our linkage map was identified, while 24,643 SNPs from the consensus cowpea linkage map were detected (Muñoz-Amatriaín et al., 2017). Thus, a total number of 28,338 SNPs were assigned to the 11 chromosomes of the common bean genome (Figure 2, Table 2, Supplementary File 3).

To capture identical markers between the two cowpea linkage maps, markers from each cowpea linkage map were paired if they hit against the same contig/scaffold of the common bean genome by taking the closest two markers together as each pair. We obtained 171 paired SNP marker loci. These paired markers were distributed unevenly across the common bean genome (Supplementary Figure 3, Supplementary File 4). This demonstrates that the current cowpea map accounts for 185 of the 37,372 SNPs mapped by Muñoz-Amatriaín et al. (2017) in cowpea, and the two maps are syntenic across all LGs to a degree.

All traits showed a continuous distribution and transgressive segregation in the F_2:3 population (Figure 1) in this study, indicating that the traits were controlled by multiple genes. The near-normal curve distribution of PL, CL, GN, and TGW suggested a polygene mode for genetic control in these traits.

Overall, 11 significant QTL for cowpea pod yield were detected on seven LGs (LG4, 5, 6, 7, 9, 10, and 11) using the CIM method (Figure 3, Table 3, and Supplementary Figure 4). Most of these QTL were assigned on their respective LGs. On LG11, three QTL (Qcpl-4, Qcgn-2, and Qctgw-4) were located between the regions of 48.5–75.6 cM. On LG4, two QTL (Qcpl-1 and Qctgw-1) were found between the regions of 56.1–128.9 cM. On LG9, one QTL (Qcpl-3) was detected on LG9 at the position of 57.9 cM. Noticeably, Qcpl-3 located at 57.9 cM had the highest LOD value of 18.36 and correspondingly had the highest contribution to phenotypic variation at 17.32%. On LG10, two QTL (Qccl-1 and Qctgw-3) were detected at positions 4.4–68.2 cM. The other three LGs, (LG 5, 6, and 7,) contained one QTL locus, including Qcgn-1 on LG5 at 59.8 cM, Qcpl-2 on LG6 at 96.3 cM, and Qctgw-2 on LG7 at 77.4 cM. Most of the 11 QTL showed negative additive effects except Qcpl-2 (1.70) and Qccl-1 (0.08).

From these 11 QTL, four QTL for cowpea pod length (Qcpl-1≈Qcpl-4) explaining 42.64% of the total phenotypic variation and two major loci (explaining > 15% of the total variation) were detected. Moreover, we detected two major-effect QTL from the above mentioned four QTL that had a phenotypic effect (R²) of more than 15%, including Qcpl-3 (17.32%), and Qcpl-4 (15.45%).

The mapped SNP flanking sequences of V. unguiculata showed significant alignment with similar sequences in the adzuki bean (Vigna angularis) (Figure 4). Between V. unguiculata and V. angularis, 5,616 of 17,996 (31.2%) mapped SNP tags (94 bp per tag) were found, which had a total length of 527.9 Kb (Supplementary Table 2).

To uncover the synteny information between cowpea and the closely related species adzuki bean, all tags of the mapped SNPs of V. unguiculata were compared with the V. angularis genome. Based on the maximum number tags between the two species, all chromosomes of V. angularis were assigned to at least one linkage group of V. unguiculata, except chromosome 6 of V. angularis, which had two assigned LGs of cowpea (Supplementary Table 2), possibly due to translocation on the chromosome. Eleven LGs of V. unguiculata aligned to eleven chromosomes of V. angularis (Supplementary File 5), caused by two LGs (LG 2 and 6) sharing one single chromosome of V. angularis (chromosome 6). These conserved regions should represent homologous regions between the two genomes. Our result implies that, although both V. unguiculata and V. angularis contain the same number of chromosomes (2X = 2N = 22), the chromosomes of cowpea evolved independently and exhibit different rearrangements from the chromosomes of V. angularis.

In addition, 93 RAD-tags located in QTL regions of the four traits (PL, TGW, GN, and CL) were analyzed by BLASTN in the NCBI databases. Of these, eight RAD tags were found to have homology domains in two species (V. radiata and V. angularis) in genus Vigna (Supplementary File 6). These tags were mainly related to two traits (pod length and TGW). Noticeably, two RAD-tag (M26270 and M26271) sequences were revealed to show similarity with the sequences of E3 ubiquitin protein ligase DRIP2-like (LOC106770573), which targets transcription factors for degradation and thus regulates gene expression in response to environmental or hormonal signals (Qin et al., 2008).

Molecular markers and genetic maps play a significant role in the genetic mapping of important traits in crops (Zhou et al., 2014). With recent significant advancements, RAD-seq technology as well as other technologies, such as GBS based DArTseq markers has been employed as a useful tool for high-density genetic mapping and 30 QTL mapping (Baird et al., 2008; Raman et al., 2014, 2016; Zou et al., 2016). This technology can discover thousands of markers in the model or other organisms regardless of genome size or reference genomes (Davey and Blaxter, 2010). Currently, the RAD-seq technique has been gradually applied to several legume species, such as chickpea (Gaur et al., 2012) and peanut (Gupta et al., 2015). However, to the best of our knowledge, no study has ever been reported on RAD-seq-based SNP discovery technique in cowpea. In the present study, we constructed and characterized a high-density SNP linkage map in cowpea using the RAD-seq method. A total of 383,077 RAD tags were obtained from the two cowpea parents and 168 F_2:3 offspring individuals. Then, we identified 34,868 novel SNPs distributed throughout the whole cowpea genome, providing genetic variation resources for the genome-wide association investigations in the future. Of the discovered SNPs, 17,996 SNPs were successfully assigned to 11 LGs corresponding to 11 chromosome pairs of the cowpea genome. The map spanned 1194.25 cM, showing 0.066 cM/SNP locus. Therefore, the current version of the population-specific linkage maps for V. unguiculata was a considerable improvement compared with the formerly published linkage maps of cowpea (Table 4). The analysis of synteny between our linkage map and a recent map of cowpea (Muñoz-Amatriaín et al., 2017) revealed that the two maps show relatively high-level synteny across all cowpea LGs. This study generated a high density RAD-based linkage map of cowpea to date. The inter-marker distance of 0.066 cM indicated it would be beneficial for locating sequence scaffold on the physical map of the cowpea genome sequence (Deokar et al., 2014). Moreover, this high-density linkage map would be an ideal genetic map for QTL mapping, gene locating, and comparative genomics analysis for agriculturally important traits of cowpea. However, there are limitations for the utilization of this map because it was based on the SNPs between only two cowpea varieties, with a SNP flanking sequence of only 94 bp.

The QTL related to domestication traits are non-randomness distribution across the crop genome (Gepts, 2004). This can be attributed to its association with “cultivation magnetism” and the “protracted transition paradigm” of crop domestication (Allaby, 2010). The increasing yield is one of the most significant domestication targets in crops, but it is a tedious process influenced by complex environmental components (Wu et al., 2014). Yield-related traits dominate QTL studies in wheat, maize, rice and soybeans, reflecting the agronomic and economic importance of these crops (Salvi and Tuberosa, 2015). Based on path coefficient analysis, it was revealed that single pod weight and number of pods per plant had the major positive contribution toward the pod yield in Phaseolus vulgaris (Ankit et al., 2017).

Quantitative trait loci (QTL) of yield-related traits have been detected in Brassica napus (Cai et al., 2016; Zhao et al., 2016), Oryza nivara (Ma et al., 2016), and Hordeum vulgare (Mikołajczak et al., 2016). In B. napus, 14 crucial candidate genes, which might be involved in developmental processes (such as, the regulation of flowering and vegetative phase change) and biomass accumulation (such as, response to secondsary wall biosynthesis), were revealed to be associated with the yield establishment (Lu et al., 2017). In several legume crops, agricultural traits have been investigated, including adzuki bean (Kaga et al., 2008), rice bean (Isemura et al., 2010), chickpea (Das et al., 2015), and cultivated peanut (Luo et al., 2017). In Arachis hypogaea, Luo et al. (2017) identified three co-localized major QTL for pod size and weight to a 3.7 cM interval on chromosome A05. These traits were controlled by a few large QTL or by a major QTL and several minor QTL. The QTL governing domestication traits in cowpea were reported in previous studies, which included pod length, plant height, seed size, flowering time, seed coat color, flower color and other traits (Ubi et al., 2000; Peksen, 2004; Kongjaimun et al., 2012a; Xu et al., 2013, 2016).

In this study, the excavation and identification of yield-associated QTL would pave the way for cowpea breeding, particularly for marker assisted selection (MAS) in cowpea. The current high-density genetic map is useful for QTL analysis in V. unguiculata.

The level of synteny between cowpea and the closely related legume species V. angularis was assessed. The Circos plot of the cowpea linkage map vs. V. angularis genome domestrated that some syntenic SNP locations were uncovered between the cowpea LGs and the V. angularis genome. When comparing the sequences of these QTL regions with reported gene sequences in NCBI databases by BLASTN, we found that two tags (M26270 and M26271) in Qctgw-2 regions showed a homologous domain with sequences of E3 ubiquitin protein ligase DRIP2-like (LOC106770573) in V. radiata. In a previous study, DRIP2 was reported to be a C3HC4 RING–domain protein interacting with the key transcription factor DREB2A and control concentration of DREB2A protein by regulating DREB2A degradation (Qin et al., 2008). In further analysis, it was documented that the DREB2A protein was involved in the seed development in Arabidopsis thaliana, because the seed yield of 35S:DREB2A CA-b mutant was lowered to 40% of the wild type of Arabidopsis thaliana (Sakuma et al., 2006). Therefore, the region of Qctgw-2, which was associated with the trait of thousand grain weight, may have candidate genes that control the seed weight of cowpea.

The SNP type was detected to be different between the two parents in the two tags (M26270 and M26271). In the M26270-SNP tag, the SNP nucleotide is “G” in the female parent, while it is “C” in the male parent. In the M26271-SNP tag, the SNP nucleotide is “C” in the female parent, while it is “A” in the male parent. Considering the phenotypic difference of the seed between the two cowpea parents, the different SNP markers of the parents suggested the feasibility of utilizing marker-assisted selection in breeding programs targeting seed variation in cowpea.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.