Draft genome sequence of a less-known wild Vigna: Beach pea (V. marina cv. ANBp-14-03)

Awninr Kumr Singh*, A. Velmurugn, Dejyoti Sen Gupt**, Jitenr KumrRvi Kesri, Arvin Kon, Nrenr Prtp Singh, Sinryn Dm Roy,Utpl Bisws, R.Rhul Kumr, Snjy SingheICAR-Centrl Isln Agriulturl Reserh Institute,Port Blir 744101,Anmn&Nior Islns,Ini

bDivision of Crop Improvement,ICAR-Indian Institute of Pulses Research,Kanpur 208024,Uttar Pradesh,India

cDepartment of Plant Breeding and Genetics,Bhola Paswan Shastri Agricultural College,Purnea,Bihar 854302,India

dDivision of Plant Biotechnology,ICAR-Indian Institute of Pulses Research,Kanpur 208024,Uttar Pradesh, India

eICAR-National Research Centre for Biotechnology,LBS Centre,Pusa Campus,New Delhi 110012,India

Keywords:Beach pea Vigna marina NextSeq 500 Whole-genome sequencing Salinity tolerance

ABSTRACT Beach pea or beach cowpea(Vigna marina(Burm.)Merr.)belongs to the family Fabaceae.It is a close relative of cultivated Vigna species such as adzuki bean (V. angularis), cowpea (V.unguiculata),mung bean(V.radiata),and blackgram(V.mungo),and is distributed throughout the tropics.With its ability to tolerate salt stress,beach pea has great potential to contribute salt-tolerance genes for developing salt-tolerant cultivars in cultivated Vigna species.However, it is still underutilized in Vigna breeding programs. A draft genome sequence of beach pea was generated using a high-throughput next-generation sequencing platform,yielding 23.7 Gb of sequence from 79,929,868 filtered reads. A de novo genome assembly containing 68,731 scaffolds gave an N50 length of 10,272 bp and the assembled sequences totaled 365.6 Mb.A total of 35,448 SSRs,including 3574 compound SSRs,were identified and primer pairs for most of these SSRs were designed.Genome analysis identified 50,670 genes with mean coding sequence length 1042 bp. Phylogenetic analysis revealed highest sequence similarity with V. angularis, followed by V. radiata. Comparison with the V.angularis genome revealed 16,699 SNPs and 2253 InDels and comparison with the V.radiata genome revealed 17,538 SNPs and 2300 InDels. To our knowledge this is the first draft genome sequence of beach pea derived from an accession (ANBp-14-03) adapted locally in the Andaman and Nicobar Islands of India. The draft genome sequence may facilitate the genetic enhancement in cultivated Vigna species.

1. Introduction

Vigna is a pantropical genus comprising 104 species distributed widely in tropical and subtropical regions [1] and containing several species of economic importance. Beach pea(Vigna marina Burm.Merr.),belongs to the family Fabaceae and is a close relative of cultivated Vigna species such as adzuki bean(V.angularis),cowpea(V.unguiculata),mung bean(V. radiata), and blackgram (V. mungo) [2]. It is also known as nanea, dune bean, notched cowpea, and beach cowpea, is distributed throughout the tropics, and is similar in floral structure to mung bean and urd bean [3,4]. Beach pea grows well on sandy beaches[5,6].It has a determinate growth habit and tolerates salt stress [2,7]. Because it can be grown in sunny areas as a ground cover, it is especially favored for beachfront properties and prevention of coastal erosion. The latter property, along with its nitrogen-fixing ability, favors the growth of neighboring plant communities.

The Andaman-Nicobar group of islands lies between two major biodiversity hotspots: the Indian subcontinent and the Malaysian-Indonesian region. The islands thus show high biodiversity within a limited geographical area. The flora of the Andaman group of islands show closer affinity to the Indo-Myanmarese-Thai flora, while those of the Nicobar islands are closer to those of Malaysia-Indonesia [8]. The islands host a diversity of wild forms of mung bean,blackgram, beach pea, red gram, and cowpea, and the diversity of cultivated pulse crops is presently maintained in home gardens by settlers from the mainland and adjoining countries[6,9].It is desirable to study the genetic relationship of beach pea with cultivated Vigna species for application in breeding programs. Naturally occurring resistance alleles for various biotic and abiotic stresses in wild relatives of Vigna can be identified and deployed in breeding programs[10]. The availability of genomic resources in a crop species helps to accelerate its genetic improvement. In recent years,progress has been made in developing genomic resources in other Vigna species,for example mung bean[11-13].However,few genomic studies of this underutilized legume-V. marina have been reported [5,14-16]. A draft genome sequence can be used to identify molecular markers such as SNPs and SSRs,for anchoring the entire genome and also serves as a reference for next-generation sequencing (NGS)-based genotyping methods, such as genotyping-by-sequencing(GBS)[17].

2. Material and methods

2.1. Plant material and DNA extraction and quantification

The beach pea accession ANBp-14-03, collected from Car Nicobar,Andaman and Nicobar Islands,India,was used.It is a creeper with a spreading growth habit, is determinate, and has ovate leaves. The 100-seed weight is 7.13 g. It is a diploid species with chromosome number 2n = 22 and is selfpollinating. DNA was isolated from 15-day-old seedling with a DNeasy Kit (Qiagen, New Delhi, India) per the manufacturer’s protocol.The quality of the DNA was checked on an 0.8%agarose gel and quantification was performed with a by Qubit 3.0 Fluorometer (Thermo Fisher Scientific, New Delhi,India).

2.2. Preparation of a 2 × 150 bp NextSeq 500 library

Library construction and sequencing were performed at Eurofins Genomics India Private Limited, Bangalore, India. A paired-end sequencing library was prepared with a NEB Next Ultra DNA Library Prep Kit, Illumina. Briefly, approximately 1 μg of gDNA was fragmented using a focused ultrasonicator(Covaris M220,New Delhi,India)to generate a mean fragment distribution of 400 bp. The fragments were then subjected to end repair. This process converts the overhangs resulting from fragmentation into blunt ends using End Repair Mix from Illumina. The 3′ to 5′ exonuclease activity of this mix removes the 3′overhangs and the 5′to 3′polymerase activity fills in the 5′ overhangs followed by adapter ligation to the fragments. This strategy ensures a low rate of chimera(concatenated template) formation. The ligated product was size-selected using AMPure XP beads (Beckman Coulter India Private Limited, New Delhi, India). The size-selected product was PCR amplified with the index primer as described in the kit protocol.

2.3. Quantity and quality check of library on automated electrophoresis tool

The PCR amplified library was analyzed in an Agilent 4200 TapeStation (Agilent Technologies, Bangalore, India) using a High Sensitivity D5000 ScreenTape assay kit as per the manufacturer’s instructions.

2.4. Cluster generation and sequencing

After obtaining the DNA concentration in Qubit Fluorometer and the mean peak size from an Agilent Tape-Station profile, the paired-end Illumina library was loaded onto a NextSeq 500 cartridge for cluster generation and sequencing. Paired-end sequencing allows template fragments to be sequenced in both forward and reverse directions. The adapter kit reagents were used to bind samples to complementary adapter oligos in a paired-end flow cell. The adapters were designed to allow selective cleavage in forward strands after resynthesis of the reverse strand during sequencing.The copied reverse strand was then used to sequence from the opposite end of the fragment.

2.5. Read generation, filtering and processing

The paired-end reads of 2 × 150 bp chemistry were generated on NextSeq 500. The raw reads generated were filtered using Trimmomatic (v 0.35) (http://www.usadellab.org/cms/index.php?page=trimmomatic) with quality value ＞30 and adapters were trimmed away. The filtered high-quality reads were assembled into a scaffold using Velvet (v.1.2.10) (https://omictools.com/velvet-tool) on optimized kmer value of 63 at which assembly with the highest continuity was found.

2.5.1.Identification of SSRs, SNPS, and InDels

The draft assembly was used to search for simple sequence repeats(SSRs)loci using MISA(http://pgrc.ipk-gatersleben.de/misa/). SSR flanking primer sequences were designed with Primer3 2.3.6 [18]. High-quality reads were mapped against Vigna angularis(GenBank accession number GCA_001190045.1)and Vigna radiata (GenBank accession number:GCA_000741045.2) genome using BWA mem (Version 0.7.12-r1039) with default parameters including minimum seed length 19 and penalty for mismatch 4. The alignment was obtained in BAM file format, which was used for SNP discovery. The alignment was obtained in BAM file format and used for SNP discovery. The standard script pipeline of Samtools mpileup was used with default parameters,but the minimum read depth was set to 15 to call SNPs and InDels.

2.5.2. Gene function prediction

Gene predictions were made with AUGUSTUS. Arabidopsis thaliana was used as a model to predict genes in the assembled scaffold. Functional annotation of the genes was performed using BLASTx searches of their sequences against NR(non-redundant protein database)in NCBI.

3. Results

3.1. Genome sequencing and assembly

Sequencing V. marina with Illumina NGS technology generated 92,144,442 reads. The raw reads generated were filtered using Trimmomatic (v 0.35); raw and filtered read statistics are presented in Table 1.Assembly statistics are presented in Table 2.A total of 68,731 scaffolds with N50 value of 10,272 bp and average size of 5320 bp were generated. The maximum and minimum sizes of the scaffolds were 73,319 bp and 700 bp respectively.

3.2. Detection of SSRs

SSRs, also known as microsatellites, are tandem repeated motifs of 2-6 bases and serve as the most important molecular markers in population and conservation genetics, molecular epidemiology, and gene mapping. SSRs were detected using Microsatellite Identification Tool (MISA v1.0). The statistics of the detected SSRs are presented in Table 3. A total of 35,448 SSRs were identified. The proportions of di-, tri- , tetra-, penta-, and hexanucleotide repeats were 47%, 49%, 3%, 7%, and 5%, respectively. Primers were designed for most of these SSRs (Table S1). Ten different dinucleotide repeats were identified, with repeats ranging from 10 to 30 (Table S2). The total number of dinucleotide repeats was 16,588 (46.79%)(Table 3). The most prevalent dinucleotide repeat was AT (51.56%), followed by TA (32.89%) and TC (3.45%) (Table S2). Sixty different trinucleotide repeats were found, with repeatnumbers ranging from 5 to 22 (Table 4). The total number of trinucleotide repeats was 17,283 (48.75%) (Table 3). The most prevalent trinucleotide repeat was AAT (13.42%), followed by TTA (11.77%) and ATA (10.25%) (Table S2). Ninety-nine different tetranucleotide repeats were identified,with repeats ranging from 5 to 15 (Table S2). The total number of tetranucleotide repeats was 1175 (3.31%) (Table 3). The most prevalent tetranucleotide repeat was AAAT(14.64%),followed by TTTA(8.60%)and AATA(6.89%)(Table S2).Pentanucleotide repeats were of 100 types with repeats ranging from 5 to 10(Table S2). The total number of pentanucleotide repeats was 244 (0.68%) (Table 3). The most prevalent pentanucleotide repeat was AAAAT (13.52%), followed by TTTAT (9.01%) and TTTTA (6.15%) (Table S2). There were 102 different hexanucleotide repeats, with repeat numbers ranging from 5 to 8 (Table S2). The total number of hexanucleotide repeats was 158 (0.44%)(Table 3). The most prevalent hexanucleotide repeat was TTATTT(20.25%),followed by AAAAAT(10.13)and AAAATA(3.80%)(Table S2).

Table 1-Raw and filtered data statistics of beach pea(Vigna marina)cv.ANBp-14-03.

Table 2-Assembly statistics for beach pea(Vigna marina)cv.ANBp-14-03.

3.3. Identification of SNPs and InDels

The statistics of the identified SNPs are presented in Fig. 1. Alignment with V. angularis revealed 16,699 SNPs and 2253 InDels and that with the V. radiata genome sequence revealed 17,538 SNPs and 2300 InDels (Fig. 1, Tables S3, S4).

3.4. Gene prediction and annotation

Gene predictions are presented in Fig. 2. In 50,670 genes, the mean coding sequence length was 1042 bp. The maximum and minimum gene length were 13,905 and 200 bp, respectively. The statistics of the predicted genes are presented in Table 4. A total of 14,470 genes could be functionally annotated (Table S4). In annotation, the greatest number of hits (number of similar or identical sequences in NR database) was with V. angularis. Large numbers of hits with mung bean and common bean (Phaseolus vulgaris) sequences were also found. In the GO category Biological Process, the first five processes were cellular process, metabolic process, biological regulation, regulation of biological process, and response to stimulus (Fig. 3). The first five functions in the GO category Molecular Function were binding, catalytic activity, molecular function regulator, transcription regulator activity, and translation regulator activity (Fig. 3). In the GO category Cellular Components, the first five components were membrane, membrane part, cell, cell part, and organelle (Fig. 3). Six genes encoding different pyrophosphatases were found with high sequence similarity with V. radiata (Table S5). The respective enzymes and their corresponding genes in V.radiata are as follows: ‘inorganic pyrophosphatase 2-like’(gene id:LOC106773663), ‘inorganic pyrophosphatase 2-like’ (gene id:LOC106774221), ‘dCTP 189 pyrophosphatase 1-like’ (gene id:LOC106767643), ‘dCTP pyrophosphatase 1’ (gene id:LOC106758825), ‘soluble inorganic pyrophosphatase 4’(gene id:LOC106757746), ‘soluble inorganic pyrophosphatase 6, chloroplastic’ (gene id:LOC106775936), and one gene matching top sequence similarity with gene (gene id:PHAVU_009G248400g) encoding plasma membrane H+192 transporting ATPase enzyme in P.vulgaris(Table S5).

Fig. 1-Numbers of SNPs and InDels in beach pea(Vigna marina)cv.ANBp-14-03.

3.5. Phylogenetic analysis

The phylogenetic relationship of the sample ANBp-14-03 with other plant species was analyzed and it was found that sample ANBp-14-03 is closely related to Vigna species (V.angularis and V.radiata)(Fig.4).

4. Discussion

Among 35,448 SSRs identified, the proportions of dinucleotide, tri-nucleotide, tetra-nucleotide, penta-nucleotide,and hexa-nucleotide unit SSRs were higher in the V. marina than in the V.radiata genome(25.83%,8.91%,0.90%,0.40%,and 0.08%), respectively [11]. However, the number of compound SSRs was lower (10.08%) than that in V. radiata (13%) [11].These findings indicate that the genome of V. marina is less complex than that of V.radiata.

Fig.2-BLAST top-hits results for prediction of genes in beach pea(Vigna marina)cv.ANBp-14-03 as compared to mung bean and other model crop species.

Fig.3- Go category:Biological process(BP),Molecular function(MF),and Cellular component(CC)summary of the annotated 14,472 genes.

To date, few efforts have been made to develop sufficient genomic resources in Vigna. This first genome sequencing effort in beach pea has generated SSRs, SNPs, InDels, and functional annotations for a huge set of genes. This information holds great promise for use in trait mapping, genomic selections,and diversity assessment.The quality of assembly generated will also advance comparative genomics in Vigna species, as whole genome sequences of prominent Vigna species including mung bean and adzuki bean are already available [11]. High-quality resequencing of various accessions within species is still needed for studying domestication and for performing marker-assisted breeding for cropimprovement [19]. In the present study, NGS has allowed the discovery of a large number of DNA polymorphisms, such as SNP and InDels markers,in a relatively short time at low cost[20,21].

Table 3-Number and distribution of 35,448 SSRs identified in the beach pea (Vigna marina) cv. ANBp-14-03 genome.

Given that of the 50,670 predicted genes only 14,470 could be functionally annotated, the remaining genes can be investigated for their role in RNAi or similar regulatory roles.In the V.radiata genome 22,427 genes could be annotated with high confidence [11]. Six genes encoding different pyrophosphatases and showing high sequence similarity with V. radiata were identified and can be used for candidate gene analysis of salt tolerance or for functional validation in V. marina. To our knowledge, the draft genome sequence described here is the first one for beach pea derived from the locally collected accession ANBp-14-03. An in-depth comparison of the new beach pea sequence with the recently published draft genomes of mung bean and adzuki bean would shed light on the mechanisms underlying salt anddrought tolerance in Vigna species. It would also allow marker-assisted breeding of high-yielding mung bean and urd bean cultivars with durable salt and drought tolerance.

Table 4-Detailed information about genes identified in beach pea(Vigna marina) cv.ANBp-14-03.

Fig.4- Diagrammatic representation of the phylogenetic relationship between an ANBp sample(V.marina)with other plant species based on BLAST hits.

The raw data is available at National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA445647) Bioproject with accession number PRJNA445647 and ID. 445647. Supplementary data to this article can be found online at https://doi.org/10.1016/j.cj.2019.05.007.

Acknowledgments

The authors are grateful to the Indian Council of Agricultural Research, Ministry of Agriculture & Farmers’ Welfare, Government of India,New Delhi,India and Director,ICAR-Central Island Agricultural Research Institute, Port Blair, Andaman &Nicobar Islands, India for financial support. The authors are also grateful to the Director, ICAR- Indian Institute of Pulses Research,Kanpur,Uttar Pradesh,India for providing technical support and guidance.