The genome assembly of asparagus bean, Vigna unguiculata ssp. sesquipedialis

Asparagus bean (Vigna. unguiculata ssp. sesquipedialis), known for its very long and tender green pods, is an important vegetable crop broadly grown in the developing Asian countries. In this study, we reported a 632.8 Mb assembly (549.81 Mb non-N size) of asparagus bean based on the whole genome shotgun sequencing strategy. We also generated a linkage map for asparagus bean, which helped anchor 94.42% of the scaffolds into 11 pseudo-chromosomes. A total of 42,609 protein-coding genes and 3,579 non-protein-coding genes were predicted from the assembly. Taken together, these genomic resources of asparagus bean will help develop a pan-genome of V. unguiculata and facilitate the investigation of economically valuable traits in this species, so that the cultivation of this plant would help combat the protein and energy malnutrition in the developing world.


Background & Summary
Asparagus bean (Vigna unguiculata ssp. sesquipedialis, 2n = 2× = 22) is a warm-season and drought-tolerant subspecies of cowpea (Vigna unguiculata) with a wide cultivation area in East and Southeast Asia 1 . This plant is also known as yardlong bean because of its characteristic pod that grows up to 50-100 cm in length 2 . The long pod trait is believed to be the result of intensive local domestication after it was brought to Asia from sub-Saharan Africa 3 . Unlike the grain-type subspecies common cowpea (Vigna. unguiculata ssp. unguiculata, or black-eyed pea), asparagus bean is harvested while its pod is still tender, thereby providing a very good source of protein, minerals, vitamins, and dietary fiber 4 . Due to the low requirement for cultivation management and its high nutritional value, asparagus bean is one of the top crops that help combat malnutrition and food insecurity in most developing countries 5 .
As the DNA sequencing technologies became more advanced and affordable for the past decade, previous research had mainly focused on delineating the genome of common cowpea (estimated genome size of 620 Mb 6 ). The first study of cowpea genomics was reported in 2008, in which the gene-rich space of cowpea was sequenced and assembled into 52,149 assemblies (41,260 assemblies were annotated) and 70,679 singletons 7 . Then the common cowpea (variety IT97K-499-35) genomic resources including a partial 323 Mb whole-genome shotgun assembly 8 , a 497 Mb bacterial artificial chromosome physical map 8 , and consensus genetic maps based on either 10 K 9 or 50 K single nucleotide polymorphisms (SNPs) were available 8 . A more recent research reported two survey genomes of common cowpea (varieties IT97K-499-35 and IT86D-1010) with substantially improved assembly sizes (568 Mb and 609 Mb, respectively) 10 . In addition, a draft IT97K-499-35 variety reference genome was assembled by incorporating the single molecule real-time technology, yielding an assembly size of 519.4 Mb into 722 scaffolds and 11 pseudo-chromosomes 11 . Three genetic maps were derived from either simple sequence repeat markers 12,13 or restriction-site associated DNA sequencing for asparagus bean 14 . Most of these genetic resources are focus on the grain-type cowpea, but there are many differences between the two types of cowpea, such as morphology, growing environments and parts for use 12 .
In this study, we aimed to fill the knowledge gap with regard to the asparagus bean genome and provide new genetic resources for breeding cowpea and related legume species. A schematic workflow of the research is shown in Fig. 1. In brief, a series of short-insert and large-insert Illumina libraries were sequenced on an Illumina HiSeq 4000 platform, yielding a total of 222.9 Gb clean data (Table 1). Since the genome size of asparagus bean was estimated to be about 590 Mb using the K-mer distribution analysis (Table 2) (Fig. 2), the clean data used for genome assembly represented about 340× coverage. The software SOAPdenovo 15 was used to generate a draft contig assembly of 549.8 Mb with a contig N50 size of 15.2 kb (Table 3). After scaffolding and gap closing, the final asparagus genome was 632.8 Mb (549.81 Mb non-N size) in size with scaffold N50 size of 2.7 Mb (Table 3). We also obtained 536,824 high-confident SNPs from the whole-genome sequencing data of 97 asparagus bean F2 individuals and two parents from a well-controlled selfing population. These SNPs were used to construct a high-density genetic map for asparagus bean, in which 1,556 scaffolds were successfully anchored onto 11 pseudo-chromosomes (Table 4). Furthermore, the asparagus bean genome contained 294.95 Mb of transposable elements, accounting for 46.47% of the assembly (Tables 5 and 6). The gene prediction was performed on a combination of de novo, homologous, and RNA-Seq-based approaches. It resulted in 42,609 protein-coding genes and 3,579 non-protein-coding genes, respectively (Table 7).

Materials.
All plant accessions were provided by Hubei Natural Science Resource Center for Edible Legumes in Wuhan of China. A single plant of the widely cultivated asparagus bean variety 'Xiabao II' (Vigna unguiculata ssp. sesquipedialis var. 'Xiabao II') was used for de novo sequencing and genome assembly. A F2 sequencing population was obtained for making the genetic map according to the following procedure. First, the F1 population were obtained by crossing 'Xiabao II' (male, same plant used for de novo sequencing) with a cultivar from the other subspecies, 'Duanjiangdou' (Vigna unguiculata ssp. unguiculata var. 'Duanjiangdou'; female). This step yielded 17 seeds, from which only 12 seeds survived till flowering. These F1 individuals were bagged to promote selfing, which produced 561 seeds in total (the F2 generation). Only 367 of the F2 individuals were able to germinate and mature into full plants. We selected 97 of the 367 F2 individuals for genome sequencing and genetic map construction.
Whole-genome shotgun sequencing. Young leaves were collected from a single 'Xiabao II' plant and used for genomic DNA extraction by the CTAB method 16 . About 10 µg of genomic DNA were used for library construction. Four short-insert libraries (350 bp, 445 bp, 758 bp, and 912 bp) and five large-insert libraries (2 kb, 3 kb, To ensure high-quality reads for the subsequent de novo assembly step, we filtered out the low-quality data by the following criteria: (a) reads with >2% unidentified nucleotides (N) or with poly-A structure; (b) reads with ≥40% bases having low quality for short insert-size libraries and ≥60% for large insert-size libraries; (c) reads with adapters or PCR duplication; (d) reads with 20 bp in 5′ terminal and 5 bp in 3′ terminal. Subsequently, about 222.9 Gb clean data were retrieved 17,18 , covering 341.66-fold of the estimated genome ( Table 1). The genomic DNA was extracted with the same procedure for the parents and all 97 F2 individuals in the resequencing population. Each DNA was used to construct 500 bp insert size libraries, which were then sequenced on an Illumina HiSeq 4000 platform. Each individual was sequenced to at least 4× coverage. NGSQCToolkit_ v2.3.3 19 was used to filter low-quality reads (parameters: −l 70 −s 25) and trim the poor-quality terminal bases (parameters: −l 5 −r 5). A total of 882.67 Gb clean bases were kept, which represented 99% of the raw sequencing data 17,18 . estimation of the genome size. The genome size of asparagus bean was estimated by the k-mer analysis approach using 69.42 Gb filtered short-insert sequencing data. The number of effective k-mers and the peak depth of a series of k values (17, 19, 21, 23, 25, 27 29 and 31) were generated by Jeffyfish (v2.2.6) 20 with the C-setting and the genome size was estimated to be about 590 Mb, according to the formula Genome_Size = (Total k-mers -Erroneous k-mers)/Peak_depth (Table 2). It is worth noting that this number could be an underestimate, in that the GC rich regions and repetitive sequences could not be properly resolved by k-mer analysis. Nonetheless, our estimated genome size was within the range of previously reported sizes (560.3 Mb 11~6 20 Mb 6 ). The rate of genome heterozygosity was calculated with the k-mer frequency distribution by the GenomeScope (v1.0.0) 21 and the result was around 0.77% (Fig. 2).
De novo genome assembly. Clean data from short insert-size libraries were corrected with the Error Correction program in SOAPdenovo package 15 . Genome assembly was performed based on the de Bruijn graph algorithm using SOAPdenovo package 22 by the following steps: (1) the paired-end reads of all libraries were used to construct the contig sequences while the K-mer values were set as 95 and 85 at the pregraph step and map step, respectively; (2) mapped paired reads were used to construct scaffolds; (3) The GapCloser package was used to map reads to the flanking sequences of gaps and to close gaps between the scaffolds; (4) genome sequence was randomly broken to re-scaffold with SSPACE package. Gaps were then filled again by GapCloser to obtain the final assembly. In the end, there were 54,864 out of 80,696 contigs with sizes longer than 1 kb. The total length of the contig assembly was 549.81 Mb ( Table 3). The longest scaffold was 14,145,393 bp, and a total of 5,621 scaffolds were longer than 1,000 bp 17,23-25 . The total length of the scaffold assembly was 632.8 Mb (Table 3). www.nature.com/scientificdata www.nature.com/scientificdata/ High-density genetic map construction and genome assembly anchoring. All clean data obtained from the two parents and the 97 F2 individuals were mapped to the asparagus bean scaffold assembly using the Burrows-Wheeler-Alignment tool (BWA) 26 mem algorithm. The SAM files were converted to BAM files using SAMtools 27 . Then the bam files were used to call SNP by the GATK software package 19 with parameters "-T HaplotypeCaller -stand_call_conf 30.0 -stand_emit_conf 10.0" and "-T SelectVariants -selectType    www.nature.com/scientificdata www.nature.com/scientificdata/ SNP". The SNPs were filtered using GATK with parameters as the following:-filterExpression "QD <2.0 || ReadPosRankSum <−8.0 || FS >60.0 || MQ <40.0 || SOR >3.0 || MQRankSum <−10.0 || QUAL <30" -log-ging_level ERROR-missingValuesInExpressionsShouldEvaluateAsFailing. After genotyping, the raw SNPs were filtered with the following criteria: missing rate <0.3 and heterozygous genotypes <0.5, resulting in a total of 836,933 high-confidence SNPs 23 .
For the genetic map construction, 50 SNPs were selected to generate bin markers from the two termini and middle part of each scaffold. These bin markers were grouped into 11 linkage groups by JoinMap v4.1 28 with the regression mapping algorithm. The grouped bins were then sorted and genetic distance was calculated by MSTmap with the Kosambi model 29 . According to this linkage map, scaffolds were anchored onto 11 pseudo-chromosomes. The SNPs were then assigned chromosome positions and a sliding window method (window size of 50 SNPs; step size of one SNP) was adopted to identify recombination events for each individual. All the recombination sites were merged and sorted with 20 kb intervals 30 . In the end, the filtered 536,824 SNPs 23 were combined into 2,013 bins 23 . According to the distribution of bins, low-recombination regions were indicated (Fig. 3). These were used to construct 11 linkage maps, resulting in 2180.14 cM spanning the whole genome 23 , which was within the reported genetic map size ranged from 643 cM 31 to 2,670 cM 32 . It is worth noting that the map lengths could be inflated by genotyping errors, as well as some biological phenomena causing double    www.nature.com/scientificdata www.nature.com/scientificdata/ recombinations 33 . The sliding window method and bin markers were used to reduce genotype errors. Since the parents were not 100% homozygous, the F1 plants were not identical, which might result inflated genetic size. In addition, 1,556 scaffolds with 597.52 Mb were anchored 23 , accounting for 94.42% of the assembled genome (Table 4). transposable elements annotation. Transposable elements (TEs) annotation were performed by a combination of homology-based and de novo prediction approaches. Homology-based approach involved searching commonly used databases for known TEs at both DNA and protein level. With default parameters, RepeatMasker 3.3.0 34 was used to identify TEs against the Repbase TE library 18.07 35 and RepeatProteinMask 34 was used to identify TEs at the protein level in the genome assembly. For de novo prediction, RepeatModeler software (http:// www.repeatmasker.org/) was used in constructing the de novo repeat library. Tandem repeats were then predicted by TRF 36 with parameters set to "Match = 2, Mismatch = 7, Delta = 7, PM = 80, PI = 10, Minscore = 50 and MaxPeriod = 2000". In total, we identified 294.95 Mb of the transposable elements, accounting for 46.47% of the asparagus bean genome (Tables 5 and 6). Among all TEs, long terminal repeat (LTR), which are important determinants of angiosperm genome size variation, constituted 19.24% of the assembled genome. DNA TEs accounted for 7.2% of the total sequence. Gene annotation. We used de novo, homology and RNA-Seq-based prediction methods to annotate protein-coding genes in the asparagus bean genome. Three de novo prediction programs, Augustus 37 , Genscan 38 and GlimmerHMM 38 were used to annotate protein-coding genes while gene model parameters were trained from Arabidopsis thaliana. For homology-based prediction, protein sequences of all the protein-coding genes of eleven species including common bean (Phaseolus vulgaris), soybean (Glycine max), pigeonpea (Cajanus cajan), chickpea (Cicer arietinum), mungbean (Vigna radiata), adzuki bean (Vigna angularis), lotus (Lotus japonicus), medick (Medicago truncatula), Arabidopsis (Arabidopsis thaliana), grape (Vitis vinifera), and rice (Orzya sativa), were first mapped to the asparagus bean genome using TblastN with the parameter E-value = 10 −5 . GeneWise 39 was then used to predict gene structure within each protein-coding region. RNA-Seq data of root and stem tissues 40 were aligned to the asparagus bean genome using TopHat on default settings. Finally, the predicted genes were merged by EvidenceModeler (EVM) 41 to generate a consensus and non-redundant gene set. This process produced 42,609 protein-coding genes with an average length of 3,156 bp (Table 7).
With BLASTP (E-value ≤ 10 −5 ), gene functions were assigned according to the best hit of alignment to SwissProt 42 , TrEMBL 43 , and KEGG 44 database. Functional domains and motifs of asparagus bean genes were determined by InterProScan 45 , which analyzed peptide sequences against protein databases including SMART, ProDom, Pfam, PRINTS, PROSITE and PANTHER. Gene Ontology (GO) terms for each gene were extracted from the corresponding InterPro entries. The result showed that 75.40% (32,126) of the total genes were supported by TrEMBL, 56.22% (23,953) by Swiss-Prot, and 59.27% (25,254) by InterPro. In addition, 10,096 (23.69%) genes could not be functionally annotated with current databases ( Table 8).
The tRNA genes were identified by tRNAscan-SE software 46 with default parameters. The rRNA genes were identified based on homology search to previously published plant rRNA sequences using BLASTN with parameters of "E-value = 10 −5 ". The snRNA and miRNA genes were identified by INFERNAL v1.0 47 software against the Rfam database with default parameters. In all, 3,579 non-protein-coding genes were identified in the asparagus bean genome, including 1593 tRNAs, 1,076 rRNAs, 350 snRNAs, and 210 microRNAs (Table 9).

Data Records
The authors declare that all data reported here are fully and freely available from the date of publication. The data resulting from each experimental and analytic step are indicated in a table (Table 10). The assembly genome and annotation are available at CNSA 17 , figshare 23 , GenBank and have accessions CP039345 24 to CP039355 25 . Raw read files of genome sequencing are available at NCBI Sequence Read Archive 18 and CNSA 17 . The SNP sets of each pseudo chromosome, the anchored scaffolds information, the filtered SNPs set identified by GATK, the www.nature.com/scientificdata www.nature.com/scientificdata/ information of bin markers and the linkage map constructed by bin markers are deposited in figshare 23 . The RNA-seq data was deposited in figshare 40 . technical Validation DNa sample quality. DNA was quantified using 0.8% agarose gel electrophoresis and Qubit Fluorometer (Invitrogen, US). DNA concentrations were normalized to 100 ng/µl for subsequent library construction. assessment of the genome assembly and annotation. Completeness of the genome assembly was assessed with default settings using the Benchmarking Universal Single-Copy Orthologs (BUSCO) 48 approach with a total of 1440 orthologue groups from the Embryophyta Dataset. The results showed that 93.2% of the core orthologs could be found in the asparagus bean genome, indicating a high-integrity assembly superior to the other four legume genomes. We aligned the raw reads from short insert-size sequencing back to the assembly and showed that approximately 94.88% of short reads could be successfully mapped. Furthermore, a previously reported high-density linkage map ("ZZ" linkage map v.2) 5 was used to assess the quality of anchored scaffolds. The sequences of 7,964 SNPs markers were aligned onto the 11 pseudo chromosomes using BLAT with parameters of "-fine" 49 . High accordance was shown between the assembled genome and the linkage map (Fig. 4). Whole genome comparative analysis was also conducted between this assembly and the cowpea genome available from Phytozome 11 by the MUMmer with nucmer 50 , presenting high collinearities with four inconsistent areas, which are located in the low-recombination regions (Fig. 5).   www.nature.com/scientificdata www.nature.com/scientificdata/ Comparison of asparagus bean genome with published common cowpea genomes. A comparison was performed (Table 11) between the asparagus bean genome and previously published common cowpea assemblies 8,10,11 . The asparagus bean genome assembly (549.81 Mb, non-N) was significantly larger than the first published IT97K-499-35 b genome 8 . Its size was close to the other three common cowpea survey assemblies (IT97K-499-35 a11 ,IT97K-499-35 c and IT86D-1010) 10 . The scaffold N50 size of our asparagus bean genome was 2.7 Mb, longer than the other three genomes assembled by the next-generation sequencing technology. Moreover, the asparagus bean assembly had about 94% of the scaffolds anchored into 11 pseudo-chromosomes according to the high-density genetic map. In addition, a set of 42,287 common cowpea coding sequences (CDS) derived from the single molecule real-time technology 11 could be blasted back to our asparagus bean genome with 90% similarity. All these results showed that the asparagus bean genome was of high quality.