Transcriptome Sequencing of Mung Bean (Vigna radiate L.) Genes and the Identification of EST-SSR Markers

Abstract

Mung bean (Vigna radiate (L.) Wilczek) is an important traditional food legume crop, with high economic and nutritional value. It is widely grown in China and other Asian countries. Despite its importance, genomic information is currently unavailable for this crop plant species or some of its close relatives in the Vigna genus. In this study, more than 103 million high quality cDNA sequence reads were obtained from mung bean using Illumina paired-end sequencing technology. The processed reads were assembled into 48,693 unigenes with an average length of 874 bp. Of these unigenes, 25,820 (53.0%) and 23,235 (47.7%) showed significant similarity to proteins in the NCBI non-redundant protein and nucleotide sequence databases, respectively. Furthermore, 19,242 (39.5%) could be classified into gene ontology categories, 18,316 (37.6%) into Swiss-Prot categories and 10,918 (22.4%) into KOG database categories (E-value < 1.0E-5). A total of 6,585 (8.3%) were mapped onto 244 pathways using the Kyoto Encyclopedia of Genes and Genome (KEGG) pathway database. Among the unigenes, 10,053 sequences contained a unique simple sequence repeat (SSR), and 2,303 sequences contained more than one SSR together in the same expressed sequence tag (EST). A total of 13,134 EST-SSRs were identified as potential molecular markers, with mono-nucleotide A/T repeats being the most abundant motif class and G/C repeats being rare. In this SSR analysis, we found five main repeat motifs: AG/CT (30.8%), GAA/TTC (12.6%), AAAT/ATTT (6.8%), AAAAT/ATTTT (6.2%) and AAAAAT/ATTTTT (1.9%). A total of 200 SSR loci were randomly selected for validation by PCR amplification as EST-SSR markers. Of these, 66 marker primer pairs produced reproducible amplicons that were polymorphic among 31 mung bean accessions selected from diverse geographical locations. The large number of SSR-containing sequences found in this study will be valuable for the construction of a high-resolution genetic linkage maps, association or comparative mapping and genetic analyses of various Vigna species.

About the Authors:

Honglin Chen

Affiliation: The National Key Facility for Crop Gene Resources and Genetic Improvement, Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

Lixia Wang

Affiliation: The National Key Facility for Crop Gene Resources and Genetic Improvement, Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

Suhua Wang

Affiliation: The National Key Facility for Crop Gene Resources and Genetic Improvement, Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

Chunji Liu

Affiliation: CSIRO Plant Industry, Queensland Bioscience Precinct, Queensland, Australia

Matthew Wohlgemuth Blair

Affiliation: Department of Agricultural and Environmental Sciences, Tennessee State University, Nashville, Tennessee, United States of America

Xuzhen Cheng

* E-mail: [email protected]

Affiliation: The National Key Facility for Crop Gene Resources and Genetic Improvement, Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

Introduction

Mung bean belongs to the Vigna genus within the Phaseoleae tribe and is a diploid crop (2n = 2x = 22) with a genome size of approximately 560 Mb. It is widely grown in Asia as an important nutritional dry grain, food legume pulse that is complementary to rice for the balanced nutrition it provides to millions of people across China, Cambodia, Laos and Vietnam, to name a few of the countries where the crop is grown. Mung bean is thought to have health promoting and nutritional characteristics in the diet and can be used to improve soil fertility given high rates of nitrogen fixation [1]. In addition, mung bean vegetable sprouts are popular in Asian cuisine, and they are good sources of protein, fiber, vitamin C and minerals.

Studies in genetic diversity, map-based cloning and molecular breeding of mung bean have lagged behind other legume crops due to the lack of genomic information for this pulse crop species [2]. Previous efforts to develop molecular markers for mung bean have not generated sufficient markers for linkage saturated map construction because they were either monomorphic or not fully informative for bi-parental mapping populations. This has been the case for RFLP [3], RAPD [4], AFLP [5], CAPS [6] and SNP [7] markers. Therefore, microsatellites (also known as SSRs based on their Simple Sequence Repeat core) are a logical choice for broadening the scope of markers available to mung bean researchers.

A promising source of SSR markers for legumes is found in the EST sequences generated by traditional or full transcriptome evaluation techniques. Sanger sequence ESTs in legumes as in other crops have for the most part contained an abundance of SSR repeat types [8]. With up-to-date technology, next generation sequencing of the transcriptome, a method often called RNA-seq, provides deeper and broader coverage of the transcriptome than traditional Sanger sequencing [9]. In addition, RNA-seq provides a lower background to signal ratio, better coverage of adenylation signals and a larger dynamic range of gene expression levels for mRNA evaluation than previous sequencing methods [10,11].

RNA-seq technology has been successfully and ubiquitously applied to both model and non-model organisms [12–14]. In mung bean, transcriptome studies have been limited and to date only 454 FLX rather than Illumina sequencing has been used with this species. Therefore, the overall goal of this research was to conduct transcriptome analysis with RNA-seq and to obtain usable EST-SSR markers from sequencing with Illumina technology. There are few reports on the development of SSR markers in mung bean to date. Somta et al. (2011) designed 157 genic microsatellite markers in mung bean but these were of low polymorphism [15]. In follow up studies, Moe et al. (2011) identified 1,630 SSR loci from mung bean mRNAs of the genotype Jangan derived from 454 sequencing technology [16]; and Gupta et al (2014) designed 1,742 SSR markers from EST sequences of the same variety [17]. However, far fewer SSR markers are reported in mung bean than in common bean [8,18], chickpea [19–20], pigeon pea [21] and soybean [22].

The full objective of this study was to use RNA-seq technology and Illumina based transcriptome evaluation of two mung bean cultivars to develop EST-based SSRs for the crop given the low number of markers available for the species. We characterize the distribution of SSR motifs in the sequences generated and validate a group of EST-SSR for further use in diversity analysis. We discuss the utility of the microsatellite markers for comparative mapping.

Materials and Methods

Plant material

A total of 33 mung bean accessions were used in this study including the varieties, ZL1 and V6, for RNA-seq and an additional 31 accessions for genetic diversity analysis (S1 Dataset). Of the 31 accessions, 8 genotypes were used for initial screening and validation of marker polymorphism. These mung bean accessions were obtained from the National Center for Crop Germplasm Resources Preservation located in Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China and were grown for DNA and RNA extractions in a research field at the same location.

RNA extraction

Tissue samples of roots, stems and leaves were collected at 15 days after sowing and quickly frozen in liquid nitrogen for storage at -80°C. RNA from each of the samples was isolated using the Trizol Reagent with manufacturer’s instructions (Invitrogen, Life Technologies, Carlsbad, USA). Total RNA was then treated with RNase-free DNase I (Takara, Kyoto, Japan) for 30 min at 37°C to remove residual DNA. RNA quality was verified using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) and was also checked by RNase free agarose gel electrophoresis. The concentration of the total RNA was further quantified with a RNA NanoDrop (Thermo Fisher Scientific Inc., Waltham, MA, USA).

cDNA library construction

Aliquots of 20 μg each of total RNA from the two different mung bean cultivars were separately processed for cNDA library construction. In both cases, a concentration of ≥ 400 ng/μl, OD260/280 = 1.8~2.2, RNA 28S:18S ≥ 1.0, and RNA Integrity Number (RIN) ≥ 7.0 was used for the preparation of a cDNA libraries. Poly-T oligonucleotide labeled magnetic beads (Illumina Inc., San Diego USA) was used to isolate poly (A) mRNA from the total RNA. Subsequently, the isolated mRNA was purified and fragmented into smaller pieces (200–700 nt) using divalent cations at 94°C for 5 min. First strand cDNA was synthesized with SuperScript II reverse transcriptase and random primers using the small fragment RNAs as templates. Second-strand cDNA synthesis was carried out using GEX second strand buffer, dNTPs, RNase H and DNA polymerase I. The cDNA fragments were further processed with end repair and phosphorylation using T4 DNA polymerase, Klenow DNA polymerase, and T4 polynucleotide kinase. The repaired cDNA fragments were 3’ adenylated using Klenow enzyme (Exo-) before end-ligating with Illumina paired-end adapters. The products from this ligation reaction were electrophoresed on a 2% (w/v) TAE-agarose gel and purified to select templates of different sizes for downstream enrichment. Only cDNA fragments of 200 bp (±25 bp) were excised from the gel and subjected to PCR. Thermocycling enrichment consisted of 15 cycles of PCR amplification performed using PCR primers PE1.0 and PE2.0 with Phusion DNA Polymerase.

Illumina Sequencing, data filtering and de novo assembly

The new cDNA libraries of the two mung bean cultivars were sequenced with Illumina paired-end sequencing technology [23] and an Illumina Hiseq 2000 sequencer which automatically collected the data and generated FASTQ files (.fq) containing raw data for all the reads. The files for ZL1 and V6, based on cultivar, were submitted to the sequence read archive (SRA) database at Genbank (www.ncbi.nlm.nih.gov), where they were combined and given the accession number SRP043316. The raw data was stringently filtered for preliminary assembly. All reads with more than 10% of bases with a poor quality score (Q<20), or non-coding RNA (such as rRNA, tRNA and miRNA), as well as ambiguous sequences containing an excess of “N” nucleotide calls or adaptor contamination, were removed. We also discarded the reads that did not pass the Illumina failed-chastity filter according to the relation “failed-chastity ≤ 1”, with a chastity threshold of 0.6 on the first 25 cycles. After this, de novo transcriptome assembly was performed with the software Trinity [24] by uploading the high-quality reads onto a computer for further analysis to 1) reduce the graph complexity by resolving repetitive sequences shorter than the read length in the graph; 2) clip the short tips in the graphs with lengths less than 2K (58 bp); 3) filter the low-coverage links that appeared only once along with their related edges; 4) merge the detected bubbles into a single path if the sequences of the parallel paths had a difference of fewer than four base pairs with >90% identity. After all these steps, the connections on the simplified graphs were broken at any repeat boundaries. These bioinformatics processes resulted in sequences without redundancy that contained the least amount of “N” nucleotide calls un-extended on either end. Only these stringently compiled sequences were defined as unigenes.

Unigene annotation and classification

The annotation of unigenes was performed using various bioinformatics procedures. The unigenes were aligned with BLASTX to four protein databases (NCBI non-redundant or Nr proteins, Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes or KEGG and euKaryotic Ortholog Groups or KOG) and one nucleotide database (NCBI nucleotide or Nt sequences) with an E-value threshold of 1.0E-5 for all except KOG with a threshold of 1.0E-3 [25,26]. The proteins with highest sequence similarity were retrieved and annotated to each unigene. With nucleotide based annotation, Blast2GO [27] software was used to obtain GO annotation categories defined by molecular function, cellular component and biological process ontologies. The KOG database was used to predict possible functions while pathway assignments were determined with KEGG.

EST-SSR search and primer design

The MIcroSAtellite (MISA) search engine (http://pgrc.ipk-gatersleben.de/misa) was employed for SSR mining and identification. The minimum numbers of repeats used for selecting the EST-SSRs were ten for mono-nucleotide based loci, six for di-nucleotide loci, five for tri-nucleotide loci and three for all larger repeat types (tetra- to hexa-nucleotide motifs). SSR marker primer pairs were designed based on sequences flanking the selected microsatellite loci using the software package Premier 5.0 (PREMIER Biosoft International, Palo Alto, CA) with targeted sizes of PCR products in the range between 100 to 300 bp.

Marker validation and genomic DNA extraction

Validation of the EST-SSR markers was conducted with the 31 mung bean accessions mentioned previously. Genomic DNA was extracted from young leaves of these accessions using the Hexadecyl trimethyl ammonium Bromide (CTAB) extraction method [28]. DNA quality was evaluated on a 1.0% agarose gel electrophoresis. The working concentration of DNA was adjusted to 50 ng/ml for use in marker evaluations. Amplification was performed in 20 μl volume reactions containing 0.5 U of Taq DNA polymerase, 1 × PCR Buffer II, 1.5 mM MgCl2, 25 μM of dNTP, 0.4 μM of each primer, and 50 ng of genomic DNA. Microsatellite loci were amplified on a Heijingang Thermal Cycler (Eastwin, Beijing, China). PCR amplification was carried out with the following cycling conditions: one cycle of 4 min at 94°C, 30–35 cycles at 94°C for 30 s, 55–60°C for 30 s and 72°C for 30 s. The final extension was performed at 72°C for 10 min. The PCR products were analyzed by 8.0% non-denaturing PAGE (Polyacrylamide gel electrophoresis) using silver staining. Fragment sizes were estimated based on the 1 Kb size marker as a DNA ladder (Promega, Madison, WI, USA).

Genetic similarity analysis

A distance tree was built based on a genetic similarity matrix for the 31 mung bean accessions and branch support was estimated with 10,000 bootstraps. The number of alleles (Na), observed heterozygosities (Ho) and polymorphism information content (PIC) for each of the EST-SSR markers were calculated using the software POPGEN 1.32 [29]. The cluster analysis of genotypes was carried out based on Nei’s unbiased measures of genetic distance by using the unweighted pair-group method with arithmetic average (UPGMA) and coefficients of genetic similarity for the mung bean accessions calculated using the same program [29].

Results

Sequencing and de novo assembly of Illumina paired-end reads

A total of 52.7 and 51.7 million paired-end raw reads were generated in Illumina next generation sequencing runs for the ZL1 and V6 varieties, respectively. After removal of the low quality reads, 51.9 and 50.9 million clean reads remained, with GC content of 43.1% and 44.8% in ZL1 and V6, respectively. In terms of sequence quality, ZL1 and V6 had 98.4% and 98.2% of Q20 or above bases and 94.5% and 93.9% of Q30 or above bases, respectively.

The combined sequence length of the Illumina reads was 10.3 Gb and could be assembled de novo into 48,693 unigenes and 83,542 individual transcripts. The average length of the assembled transcripts was 1,194 bp (N50 = 1,936 bp), which was longer than the average length of the assembled unigenes (874 bp, N50 = 1,563 bp). The range in length of the assembled unigenes was from 200 bp to 20,214 bp. A total of 25,590 unigenes (52.6%) were short, with lengths no longer than 500 bp. The next two size classifications of unigenes were of similar frequency with 9,141 unigenes (18.8%) having 501 to 1,000 bp in length, and 8,643 unigenes (17.8%) with lengths ranging from 1,000 to 2,000 bp. Finally 5,319 unigenes (10.9%) were longer than 2,000 bp (Fig 1).

[Figure omitted. See PDF.]

Fig 1. Length distribution of assembled transcripts and unigenes.

A. Size distribution of transcripts. B. Size distribution of Unigenes. C. Length distribution of transcripts. D. Length distribution of Unigenes.

https://doi.org/10.1371/journal.pone.0120273.g001

Sequence annotation

For validation and annotation of the sequence assembly contigs and unique singletons, all unigenes were searched against the five databases as described earlier. A total of 48,693 unigenes provide a significant BLAST result, with 25,820 (53.0%) showing significant similarity to known proteins in the Nr sequence database, with 18,316 from Swiss-Prot (37.6%) and 17,652 from PFAM (36.3%), only 4,064 unigenes were annotated in all databases but 28,613 could be annotated in at least one while the rest (20,080 unigenes) were not annotated to the existing databases.

Assembled unigenes were classified in various ways (Fig 2). Based on Nr annotation, 19,242 unigenes (39.5%) were assigned gene ontology (GO) terms (Fig 2A). The sequences that belonged to the biological process, cellular component, and molecular function clusters were categorized into 55 functional groups. Binding (11,440, 59.5%), cellular processes (11,392, 59.2%), metabolic processes (11,049, 57.4%), catalytic activity (9,595, 49.9%) and cell part (6,524, 33.9%) were the dominant five groups respectively, however, only 1 unigene each was assigned to metallo-chaperone activity and symplast localization.

[Figure omitted. See PDF.]

Fig 2. Classification of assembled unigenes.

A. EuKaryotic Ortholog Groups (KOG) classification of assembled unigenes. B. Gene ontology (GO) classification of assembled unigenes. C. Kyoto Encyclopedia of Genes and Genomes (KEGG) classification of assembled unigenes.

https://doi.org/10.1371/journal.pone.0120273.g002

After functional prediction and classification, 10,918 unigenes (22.4%) could be classified by hits with the KOG database (Fig 2B). The KOG-annotated putative proteins were functionally classified into 26 molecular families. General function only (2,001, 19.3%), post-translational modification, protein turn over and chaperon activity (1,345, 13.0%), signal transduction (1,074, 10.4%), transcription (750, 7.2%) and intracellular trafficking (646, 6.2%) were the dominant five groups, whereas only 4 unigenes each were assigned to cell motility and un-named protein groupings. To further analyze the transcriptome of mung bean, all the unigenes were analyzed in the KEGG database (Fig 2C), where a total of 6,585 unigenes had significant matches and were assigned to 5 main categories in 244 pathways. Among these positive KEGG hit unigenes, metabolic pathways contained 3,068 unigenes, followed by genetic information processing (1,420, 21.6%), organismal systems (1,051, 16.0%), cellular processes (729, 11.1%) and environmental information processing (560, 8.5%).

Frequency and distribution of different types of EST-SSR markers

Of the 48,693 unigenes found in the current study, 10,053 (20.6%) contained one or more SSR sequences. Of these, 2,303 (4.9%) contained at least two separate SSR sequences and 978 (2.1%) contained compound SSRs of different motifs. The proportion of EST-SSR was not evenly distributed among EST-SSR unit sizes or groups. Mono-nucleotide motifs were the most abundant (4,751, 36.2%) with tetra- (2,813, 21.4%), tri- (1,915, 14.6%), di- (1,809, 13.8%), penta-nucleotide (969, 7.4%) and hexa- (877, 6.7%) nucleotide motif repeats being the next most common in consecutive order (Table 1).

[Figure omitted. See PDF.]

Table 1. Summary of the number of repeat units in mung bean EST-SSR loci.

https://doi.org/10.1371/journal.pone.0120273.t001

The number of SSR repeats per locus ranged from 3 to 23, and SSRs with three repeats were the most abundant, followed by those with ten, six and five random repeats. Motifs that showed more than 16 repeats were rare, with a frequency of only 4.4%. Among mono-nucleotide repeats, the (A/T)n repeats were far more abundant (99.7%) compared to the (G/C)n repeats. The six other main motif types were the (AG/CT)n di-nucleotide repeat (30.8%), the (GAA/TTC)n tri-nucleotide repeat (12.6%), the (AAAT/ATTT)n tetra-nucleotide repeat (6.8%), the (AAAAT/ATTTT)n penta-nucleotide repeat (6.2%), and then the (AAAAAT/ATTTTT)n hexa-nucleotide repeat (1.9%), consecutively (S2 Dataset).

Development of polymorphic EST-SSR markers in mung bean

A total of 13,134 in silico EST-SSR markers could be developed form the 10,053 SSR containing sequences using Primer3 (S3 Dataset). A subset of 200 markers was randomly chosen from these loci to validate EST-SSR marker usefulness in monitoring polymorphisms for eight mung bean accessions (S4 Dataset). Of the markers tested, 129 primer pairs (65.0%) produced clear PCR amplicons of the expected sizes, 36 markers amplified non-specific products, and 35 did not amplify any clear DNA bands. Of the successfully amplifying EST-SSR markers, 66 (or 51.2%) were polymorphic and consisted of 4 mono-, 2 di-, 33 tri-, 6 tetra-, 3 penta- and 18 hexa motif based marker (Table 2) while the other 97 markers were monomorphic. An average of 3.0, 2.5, 2.2, 2.2, 2.0 and 2.6 alleles were generated for the mono-, di-, tri-, tetra-, penta- and hexa motif markers, respectively.

[Figure omitted. See PDF.]

Table 2. The evaluation of microsatellite markers for different repeat classes.

https://doi.org/10.1371/journal.pone.0120273.t002

Gene functions of the unigene sequences containing polymorphic EST-SSRs

To determine the possible functions of the 66 validated EST-SSRs, they were subjected to BLASTn analysis with a non-redundant database of legume sequences. The results showed that most of the sequences were similar to known or hypothetical protein-encoding genes from common bean (Phaseolus vulgaris L.) and soybean (Glycine max L.) with a lesser proportion homologous to cowpea (Vigna unguiculata L. [Walp]) genes (S5 Dataset). Among the positive hits were genes for auxin efflux carrier component, dof zinc finger, F-box, gibberellin receptor, helicase, mitogen-activated and leucine-rich repeat extensin-like proteins as examples.

Phylogenetic analysis of the cultivated mung bean accessions

The 66 polymorphic EST-SSR markers developed in this study were used to assess the genetic diversity of 31 mung bean accessions from the complete geographic distribution of the crop for which a total of 154 alleles were detected and scored. The number of alleles per marker ranged from 2 to 5. Effective number of alleles per locus (Ne) varied from 1.074 (for marker MB64504) to 3.014 (MB27164) averaging 1.810, expected heterozygosity (He) ranged from 0.070 (MB64504) to 0.675 (MB9309) averaging 0.429, observed heterozygosity (Ho) varied from 0 (MB21076) to 0.897 (MB27164) averaging 0.100. Shannon's Information index (I) values ranged from 0.154 (MB64504) to 1.259 (MB27164) averaging 0.649 and PIC values ranged from 0.067 (MB17985) to 0.613 (MB25181) averaging 0.344 (Table 3).

[Figure omitted. See PDF.]

Table 3. Informativeness of EST-SSR loci following amplification from 31 geographically diverse accessions of mung bean.

https://doi.org/10.1371/journal.pone.0120273.t003

Phylogenetic relationships between the accessions grouped the 31 accessions into two main clusters in a dendogram (Fig 3). Cluster 1 was comprised of accessions from Southeast and South Asian countries such as Thailand, Vietnam, Philippines, Indonesia, Myanmar, Nepal and India. Cluster 2 was comprised of accessions from East Asian and Northeast Asian countries such as China, Japan, Korea and Russia. Results indicated that geographical distances between collection sites for the accessions were associated with Nei’s genetic distances between accessions.

[Figure omitted. See PDF.]

Fig 3. UPGMA dendrogram of 31 genotypes of mung bean.

https://doi.org/10.1371/journal.pone.0120273.g003

Discussion

Transcriptome sequencing and de novo assembly has proven to be an important tool for gene discovery in many organisms and an effective method for molecular marker development [30,31]. Our results also proved that the short reads from Illumina paired-end sequencing of mung bean cDNAs can be easily assembled and used for transcriptome analysis, marker development and gene identification even without a reference genome for the crop. The marker validation confirmed previous evaluations of SSRs in common bean [8,32] and other legume crops, where EST-SSR markers detect moderate polymorphism.

Our work complements previous analysis with 454 sequencing of two mung bean cDNA libraries which resulted in the discovery of 1,630 and 1,334 EST-SSR primer pairs from the leaves of Jangan and Sunhwa varieties, respectively [16]. Here we concentrated on the use of Illumina sequencing to develop a larger total number of in silico EST-SSR markers to increase the number of SSRs available for mung bean. We found that the EST-SSR marker validation rate was similar to the success rate of SSR development in a previous study in mung bean [33] but slightly lower than when using sequenced BAC end sequences or small-insert genomic libraries in common bean [34,35]. Comparing these legumes, the polymorphism ratio of EST-SSR markers in mung bean was slightly higher than for EST-SSR in common bean [18].

In terms of the types of motifs found in SSR loci other than the mono- and large sized repeats, we found similar results as in previous work with plant microsatellites [33]. For example, the proportions of di- and tri-nucleotide repeats were quite close (13.8% versus 14.6%) as was found in previous results [30] The relative abundance of di- and tri-nucleotide repeats in ESTs sequences has been observed in many other legumes including common bean [8], cowpea [36] and chickpea [37]. The most common tri-nucleotide repeats found in the mung bean varieties studied here were GAA/TTC followed by TCT/AGA and CTT/AAG, which are similar with previous reports in mung bean [17] and common bean [8,18], possibly indicating a shared origin among the Phaseoleae tribe. As in common bean, AG/CT motif was the most abundant repeat motif (30.8%), followed by AT/TA (28.5%) contrasting slightly with other legumes outside the Phaseoleae tribe [36–38] but similar to estimates in common bean [8].

A larger number of repeat units were generally correlated with greater allelic variability for an SSR locus. Therefore, the shorter motif loci such as those with mono- and di-nucleotides repeats usually had to possess more repeats to be of equivalent polymorphism to longer motif repeats such as those with tri-nucleotide motifs. Previous studies in legumes have mainly focused on di-, tri-, and tetra-nucleotide SSRs [18], whereas mono-nucleotide SSRs perhaps have not drawn enough attention for marker development. We found that mono-nucleotide SSRs had higher polymorphism rates than previously thought, followed by tetra-nucleotide SSRs, justifying their inclusion in future SSR evaluations.

To determine the level of polymorphism of our new EST-SSR markers, we validated 200 loci, of which 129 markers (65.0%) produced successful amplicons, which is in between previously reported success rates of 21.0% [15] and 78.0% [17] in mung bean. The failure of 35 primer pairs to generate amplicons may be due to long intervening introns which would not allow successful genomic amplification of markers based on transcribed mRNA based sequences. Alternatively the location of primers across splice sites or regions of poor sequence quality could explain non-amplification.

Despite these issues, about two thirds of the EST-SSR markers were successful, suggesting that the transcriptome sequencing was accurate and the assembled unigenes were of high quality. In terms of allele detection, only half of the successfully amplified SSRs produced more than one allele and most had no more than 4 alleles, which was in agreement with a previous study [39] for mung bean. PIC values in this study were in line with previously reported values for mung bean SSRs [17,40]. They were also similar to EST-SSRs from adzuki bean which can amplify products in mung bean [41]. One advantage of EST-SSR markers, is that they may detect valuable genetic diversity possibly associated with traits of interest for breeding because of their location in genes.

In summary, the accomplishments of our study were 1) the detection of a large number of unigenes for mung bean and 2) the discovery of over 10,000 SSR containing sequences in the transcriptome of the crop. We observed that the number and lengths of unigenes in mung bean compared favorably to previous analyses with the generation here of approximately 25 million paired-end reads for the transcriptome which assembled into over 48 thousand unigenes with an average length of over 850 bp. Sanger sequencing of cDNAs do not efficiently produce this number of unigenes or sufficient overall contig lengths because of a limitation in the depth of sequencing, even when full-length cDNA libraries are used [31,42,43]. The use of transcriptomic data for in silico microsatellite development was shown to be promising and we were able to increase the number of possible EST-SSRs tenfold compared to previous studies [33]. The newly developed SSR sequences and EST-SSR markers we made will be important resources for basic research and together with SNP resources can significantly enhance the ability to find closely linked markers for traits of interest in the molecular breeding of mung bean.

Supporting Information

[Figure omitted. See PDF.]

S1 Dataset. Germplasm accessions used in this study of mung bean.

https://doi.org/10.1371/journal.pone.0120273.s001

(DOC)

S2 Dataset. Frequencies of different repeat motifs in EST-SSRs from mung bean.

https://doi.org/10.1371/journal.pone.0120273.s002

(DOC)

S3 Dataset. Characteristics of 13,134 mung bean EST-SSR markers in this study.

https://doi.org/10.1371/journal.pone.0120273.s003

(XLS)

S4 Dataset. Primer sequences of a total of 200 SSR markers for validation.

https://doi.org/10.1371/journal.pone.0120273.s004

(XLS)

S5 Dataset. The putative proteins identified by BLASTX of 66 unigene sequences containing polymorphic EST-SSRs.

https://doi.org/10.1371/journal.pone.0120273.s005

(DOC)

Author Contributions

Conceived and designed the experiments: HLC XZC. Performed the experiments: HLC. Analyzed the data: HLC. Contributed reagents/materials/analysis tools: SHW LXW. Wrote the paper: HLC XZC. Contributed to manuscript revision: MB CJL.

Citation: Chen H, Wang L, Wang S, Liu C, Blair MW, Cheng X (2015) Transcriptome Sequencing of Mung Bean (Vigna radiate L.) Genes and the Identification of EST-SSR Markers. PLoS ONE 10(4): e0120273. https://doi.org/10.1371/journal.pone.0120273

References

1. Nair RM, Yang RY, Easdown WJ, Thavarajah D, Thavarajah P, et al. Biofortification of mungbean (Vigna radiata) as a whole food to enhance human health. J Sci Food Agr. 2013; 93: 1805–1813. pmid:23426879

2. Tangphatsornruang S, Sangsrakru D, Chanprasert J, Uthaipaisanwong P, Yoocha T, Jomchai N, et al. The chloroplast genome sequence of mungbean (Vigna radiata) determined by high-throughput pyrosequencing: structural organization and phylogenetic relationships. DNA Res. 2010; 17: 11–22. pmid:20007682

3. Young ND, Kumar L, Menancio-Hautea D, Danesh D, Talekar NS, Shanmugasundarum S. et al. RFLP mapping of a major bruchid resistance gene in mungbean (Vigna radiata, L. Wilczek). Theor Appl Genet. 1992; 84: 839–844. pmid:24201484

4. Kaga A, Ishimoto M. Genetic localization of a bruchid resistance gene and its relationship to insecticidal cyclopeptide alkaloids, the vignatic acids, in mungbean (Vigna radiata L. Wilczek). Mol Gen Genet. 1998; 258: 378–384. pmid:9648742

5. Moe KT, Gwag JG, Park YJ. Efficiency of POWERCORE in core set development using amplified fragment length polymorphic markers in mungbean. Plant Breeding. 2012; 131: 110–117.

6. Chen HM, Liu CA, George KC, Chien CM, Sun HC, Huang CC, et al. Development of a molecular marker for a bruchid (Callosobruchus chinensis L.) resistance gene in mungbean. Euphytica. 2007; 157: 113–122.

7. Van K, Kang YJ, Han KS, Lee YH, Gwag JG, Moon JK. et al. Genome-wide SNP discovery in mungbean by Illumina HiSeq. Theor Appl Genet. 2013; 126: 2017–2027. pmid:23674132

8. Blair MW, Hurtado N, Chavarro CM, Munoz-Torres MC, Giraldo MC, Pedraza F, et al. Gene-based SSR markers for common bean (Phaseolus vulgaris L.) derived from root and leaf tissue ESTs: an integration of the BMc series. BMC Plant Biol. 2011; 11: 50. pmid:21426554

9. Levin JZ, Yassour M, Adiconis X, Nusbaum C, Thompson DA, Friedman N, et al. Comprehensive comparative analysis of strand-specific RNA sequencing methods. Nat Methods. 2010; 7: 709–715. pmid:20711195

10. Nookaew I, Papini M, Pornputtapong N, Scalcinati G, Fagerberg L, Uhlen M, et al. A comprehensive comparison of RNA-Seq-based transcriptome analysis from reads to differential gene expression and cross-comparison with microarrays: a case study in Saccharomyces cerevisiae. Nucleic Acids Res. 2012; 40: 10084–10097. pmid:22965124

11. Agarwal A, Koppstein D, Rozowsky J, Sboner A, Habegger L, Hillier LW, et al. Comparison and calibration of transcriptome data from RNA-Seq and tiling arrays. BMC Genomics. 2010; 11: 383. pmid:20565764

12. Choudhary S, Sethy NK, Shokeen B, Bhatia S. Development of chickpea EST-SSR markers and analysis of allelic variation across related species. Theor Appl Genet. 2009; 118: 591–608. pmid:19020854

13. Garg R, Patel RK, Tyagi AK, Jain M. De novo assembly of chickpea transcriptome using short reads for gene discovery and marker identification. DNA Res. 2011; 18: 53–63. pmid:21217129

14. Jhanwar S, Priya P, Garg R, Parida SK, Tyagi AK, Jain M, et al. Transcriptome sequencing of wild chickpea as a rich resource for marker development. Plant Biotechnol. J. 2012; 10: 690–702. pmid:22672127

15. Somta P, Seehalak W, Srinives P. Development, characterization and cross-species amplification of mungbean (Vigna radiata) genic microsatellite markers. Conserv Genet. 2009; 10: 1939–1943.

16. Moe KT, Chung JW, Cho YI, Moon JK, Ku JH, Jung JK, et al. Sequence information on simple sequence repeats and single nucleotide polymorphisms through transcriptome analysis of mungbean. J Integr Plant Biol. 2011; 53(1): 63–73. pmid:21205180

17. Gupta SK, Bansal R, Gopalakrishna T. Development and characterization of genic SSR markers for mungbean (Vigna radiata (L.) Wilczek). Euphytica. 2014; 195(2): 245–258.

18. Blair MW, Torres MM, Giraldo MC, Pedraza F. Development and diversity of Andean-derived, gene-based microsatellites for common bean (Phaseolus vulgaris L.). BMC Plant Biol. 2009; 9: 100. pmid:19646251

19. Choudhary S, Sethy NK, Shokeen B, Bhatia S. Development of chickpea EST-SSR markers and analysis of allelic variation across related species. Theor Appl Genet. 2009; 118: 591–608. pmid:19020854

20. Garg R, Patel RK, Tyagi AK, Jain M. De Novo Assembly of Chickpea Transcriptome Using Short Reads for Gene Discovery and Marker Identification. DNA Res. 2011; 18: 53–63. pmid:21217129

21. Dutta S, Kumawat G, Singh BP, Gupta DK, Singh S, et al. Development of genic-SSR markers by deep transcriptome sequencing in pigeonpea [Cajanus cajan (L.) Millspaugh]. BMC Plant Biol. 2011; 11: 17. pmid:21251263

22. Xin D, Sun J, Wang J, Jiang H, Hu G, et al. Identification and characterization of SSRs from soybean (Glycine max) ESTs. Mol Biol Rep. 2012; 39: 9047–9057. pmid:22744420

23. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009; 10: 57–63. pmid:19015660

24. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013; 8: 1494–1512. pmid:23845962

25. Cameron M, Williams HE, Cannane A. Improved gapped alignment in BLAST. IEEE/ACM Trans Comput Biol Bioinform. 2004; 1: 116–129. pmid:17048387

26. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997; 25: 3389–3402. pmid:9254694

27. Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, et al. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005; 21: 3674–3676. pmid:16081474

28. Englen MD, Kelley LC. A rapid DNA isolation procedure for the identification of Campylobacter jejuni by the polymerase chain reaction. Lett Appl Microbiol. 2000; 31: 421–426. pmid:11123549

29. Krawczak M, Nikolaus S, von Eberstein H, Croucher PJ, El Mokhtari NE, et al. PopGen: population-based recruitment of patients and controls for the analysis of complex genotype-phenotype relationships. Community Genet. 2006; 9: 55–61. pmid:16490960

30. Huang D, Zhang Y, Jin M, Li H, Song Z, et al. Characterization and high cross-species transferability of microsatellite markers from the floral transcriptome of Aspidistra saxicola (Asparagaceae). Mol Ecol Resour. 2014; 14: 569–577. pmid:24286608

31. Zheng X, Pan C, Diao Y, You Y, Yang C, et al. Development of microsatellite markers by transcriptome sequencing in two species of Amorphophallus (Araceae). BMC Genomics. 2013; 14: 490. pmid:23870214

32. Zhang X, Blair MW, Wang S. Genetic diversity of Chinese common bean (Phaseolus vulgaris L.) landraces assessed with simple sequence repeat markers. Theor Appl Genet. 2008; 117: 629–640. pmid:18548226

33. La Rota M, Kantety RV, Yu JK, Sorrells ME. Nonrandom distribution and frequencies of genomic and EST-derived microsatellite markers in rice, wheat, and barley. BMC Genomics. 2005; 6: 23. pmid:15720707

34. Blair MW, Buendia HF, Giraldo MC, Metais I, Peltier D. Characterization of AT-rich microsatellites in common bean (Phaseolus vulgaris L.). Theor Appl Genet. 2008; 118: 91–103. pmid:18784914

35. Blair MW, Torres MM, Pedraza F, Giraldo MC, Buendia HF, et al. Development of microsatellite markers for common bean (Phaseolus vulgaris L.) based on screening of non-enriched, small-insert genomic libraries. Genome. 2009; 52: 772–782. pmid:19935925

36. Gupta SK, Gopalakrishna T. Development of unigene-derived SSR markers in cowpea (Vigna unguiculata) and their transferability to other Vigna species. Genome. 2010; 53: 508–523. pmid:20616873

37. Choudhary S, Sethy NK, Shokeen B, Bhatia S. Development of chickpea EST-SSR markers and analysis of allelic variation across related species. Theor Appl Genet. 2009; 118: 591–608. pmid:19020854

38. Tian AG, Wang J, Cui P, Han YJ, Xu H, et al. Characterization of soybean genomic features by analysis of its expressed sequence tags. Theor Appl Genet. 2004; 08: 903–913.

39. Gupta SK, Bansal R, Vaidya UJ, Gopalakrishna T. Development of EST-derived microsatellite markers in mungbean [Vigna radiata (L.) Wilczek] and their transferability to other Vigna species. Indian J Genet Pl Br. 2012; 72: 468–471.

40. Gwag JG, Dixit A, Park YJ, Ma KH, Kwon SJ, et al. Assessment of genetic diversity and population structure in mungbean. Genes & Genomics. 2010; 32: 299–308.

41. Chankaew S, Isemura T, Isobe S, Kaga A, Tomooka N, et al. Detection of genome donor species of neglected tetraploid crop Vigna reflexo-pilosa (creole bean), and genetic structure of diploid species based on newly developed EST-SSR markers from azuki bean (Vigna angularis). PLoS One. 2014; 9: e104990. pmid:25153330

42. Thiel T, Michalek W, Varshney RK, Graner A. Exploiting EST databases for the development and characterization of gene-derived SSR-markers in barley (Hordeum vulgare L.). Theor Appl Genet. 2003; 106: 411–422. pmid:12589540

43. Blair MW, Fernandez AC, Ishitani M, Moreta D, Seki M, et al. Construction and EST sequencing of full-length, drought stress cDNA libraries for common beans (Phaseolus vulgaris L.). BMC Plant Biol. 2011; 11: 171. pmid:22118559