|Home | About | Journals | Submit | Contact Us | Français|
A comprehensive transcriptome assembly for pigeonpea has been developed by analyzing 128.9 million short Illumina GA IIx single end reads, 2.19 million single end FLX/454 reads, and 18353 Sanger expressed sequenced tags from more than 16 genotypes. The resultant transcriptome assembly, referred to as CcTA v2, comprised 21434 transcript assembly contigs (TACs) with an N50 of 1510bp, the largest one being ~8kb. Of the 21434 TACs, 16622 (77.5%) could be mapped on to the soybean genome build 1.0.9 under fairly stringent alignment parameters. Based on knowledge of intron junctions, 10009 primer pairs were designed from 5033 TACs for amplifying intron spanning regions (ISRs). By using in silico mapping of BAC-end-derived SSR loci of pigeonpea on the soybean genome as a reference, putative mapping positions at the chromosome level were predicted for 6284 ISR markers, covering all 11 pigeonpea chromosomes. A subset of 128 ISR markers were analyzed on a set of eight genotypes. While 116 markers were validated, 70 markers showed one to three alleles, with an average of 0.16 polymorphism information content (PIC) value. In summary, the CcTA v2 transcript assembly and ISR markers will serve as a useful resource to accelerate genetic research and breeding applications in pigeonpea.
Pigeonpea (Cajanus cajan (L.) Millspaugh) is an important food legume crop of tropical and subtropical regions of the world. It is a diploid (2n=2x=22), with moderate genome size of 858Mbp (Greilhuber and Obermayer, 1998). The genus Cajanus comprises 32 species, most of which are found in India and Australia, and one is native to West Africa (Bohra et al., 2011a). Pigeonpea is grown in 4.67Mha across the world and India is the world’s largest producer (van der Maesen, 1990). It is an important crop in south Asia, the Caribbean, and parts of Africa and South America. Pigeonpea is a vital source of protein (with 20–22% protein by dry weight), especially in vegetarian diets (Duke, 1981).
Key limitations to sustainable pigeonpea production are several abiotic stresses (e.g. drought, salinity and water-logging) and biotic stresses (e.g. Fusarium wilt (FW), sterility mosaic disease (SMD), and pod borer insects). Addressing these limitations is critical to meeting the demands of resource-poor people where pigeonpea is grown. Although there are continuing efforts for pigeonpea improvement through conventional breeding (Reddy et al., 1978; Saxena et al., 1983; Wanjari and Patel, 2003; Saxena, 2008), and molecular breeding has a great potential to enhance crop productivity (Varshney et al., 2010), limited availability of genomic resources coupled with a narrow genetic diversity in the cultivated gene pool have been serious bottlenecks to successful molecular breeding in pigeonpea (Varshney et al., 2009a).
Genomic resources like molecular markers, genetic maps, transcriptomic, or genome sequence data are prerequisites for undertaking molecular breeding in any crop. In the case of pigeonpea, efforts have been made only recently to develop some genomic resources. These include 88860 bacterial artificial chromosome (BAC)-end sequences (BESs), 3072 BES-derived simple sequence repeat (SSR, or BES-SSR) markers, a 239 BES-SSR locus genetic map (Bohra et al., 2011a), 18353 Sanger ESTs (Raju et al., 2010; unpublished), and 1.696 million FLX/454 reads (Dutta et al., 2011). Additionally, 494353 FLX/454 reads and 128.9 million short Illumina GA IIx single end reads have also been generated (Dubey et al., 2011). The 494353 FLX/454 reads along with the 10817 Sanger ESTs available at the time were merged to generate a transcript assembly (CcTA v1) comprising 48726 contigs (Dubey et al., 2011).
To improve the resources for pigeonpea genetics research and breeding applications, the present study was undertaken to develop a comprehensive transcriptome assembly based on a hybrid approach consisting of Sanger ESTs and mRNA-Seq data from two different next-generation sequencing platforms (Illumina GA IIx and FLX/454). This new assembly (CcTA v2) is available through the Legume Information System (LIS) website at http://cajca.comparative-legumes.org/data/lista_cajca-201012.tgz. These transcript assembly contigs (TACs) were aligned to the genome sequence of soybean (Glycine max), also in the Phaseoleae tribe and separated from Cajanus by about 20Mya (Stefanovic et al., 2009). With the help of anchoring points between pigeonpea and soybean genomes, intron spanning region (ISR) markers were developed. Probable chromosomal assignments for more than 6000 ISR markers have been predicted for pigeonpea on the basis of location of the ISR markers in the soybean genome, and the locations of mapped BES-SSR loci in the pigeonpea genetic map that were used as anchor points between the pigeonpea and soybean genomes. Validation of a subset of ISR markers underlines the utility of these markers to enrich the existing pigeonpea genetic maps and identification of the quantitative trait loci (QTLs) for resistance/tolerance to biotic/abiotic stresses.
Four datasets, including 128.9 million Illumina GA IIx reads from 12 genotypes (Dataset I), 2.19 million FLX/454 reads from three genotypes (Datasets II and III), and 18353 Sanger ESTs from six genotypes (Dataset IV) were processed to generate the transcriptome assembly CcTA v2 (Table 1). The transcriptome assembly comprises 21434 TACs.
When the datasets above were analyzed individually in earlier studies, a wide range of TAC counts were reported: 4557 contigs from 9888 Sanger ESTs (Raju et al., 2010) and 43324 contigs from 1.696 million FLX/454 reads (Dutta et al., 2011). The CcTA v1 (Dubey et al., 2011), assembled from 494353 454/FLX reads and 10817 Sanger ESTs, produced an assembly of 48726 contigs. The transcriptome assembly in this study, referred to as CcTA v2, has numerous improved characteristics (Table 2). For instance, the CcTA v2 has a total of 21434 TACs, with N50 of 1510bp, while the CcTA v1 included 48726 contigs, with N50 length of only 285bp (Dubey et al., 2011). The largest TAC in CcTA v2 is 7909bp, which is almost four times larger than that in the CcTA v1, with 2067bp.
To check for microbial contamination, if any, all 21434 TACs of the transcriptome assembly (CcTA v2) were BLASTed against NCBI’s bacterial genomes database. Only 49 TACs had significant hits to bacterial genomes. These 49 TACs were further BLASTed against NCBI’s EST_others database (non-human, non-mouse). While 45 TACs had hit with pigeonpea cDNA libraries or other plant cDNAs, only four TACs (lista_cajca-201012 TACs# 9664, 18316, 21003, and 21371) still showed hit with microbial ESTs. However, three of these TACs (lista_cajca-201012 TACs# 9664, 18316, and 21371) could be mapped to the soybean genome at >90% sequence similarity and 80% coverage of query length. Therefore, these TACs also could be considered legitimate transcript contigs. The remaining one TAC (lista_cajca-201012 TAC# 21003) could also be mapped to soybean genome at a lower stringency, namely >85% sequence similarity and >50% coverage of query length. Therefore, the developed transcriptome assembly seems to be of high quality without any microbial contamination.
All TACs of CcTA v2 were aligned to the soybean genome (ftp://ftp.jgi-psf.org/pub/JGI_data/phytozome/v7.0/Gmax/) using the HMM-based alignment program Exonerate 2.2.0 (Slater and Birney, 2005), to investigate gene coverage and gene structures of pigeonpea. Using the alignment criteria of 80% identity and 50% coverage, of the 21434 TACs, 16622 (77.5%) could be aligned. These alignments can be visualized in the soybean genome browser, using the track ‘Cajanus cajan (pigeonpea) v.2’, at http://soybase.org/gb2/gbrowse/gmax1.01. It is encouraging to note that 27490 predicted soybean genes have matches, and that the overall distribution of soybean gene models covered by the CcTA v2 is relatively even across all 20 soybean chromosomes, ranging from 50 to 73% of predicted soybean genes per chromosome with pigeonpea matches (Table 3). Further, 9863 (59.3%) out of 16622 CcTA v2 TACs had two good matches with the soybean genome, consistent with the recent genome duplication in soybean with respect to pigeonpea, while an additional 20% had three or four good matches, as might be expected given the ancient whole-genome duplication event in pappilionoids ~59Mya (Schmutz et al., 2010). Details of the number of pigeonpea TACs mapped onto soybean at a given number of times are given in Supplemental Table 1. Synteny and commonalities between pigeonpea TACs and soybean genes can be visualized genome-wide in soybean Gbrowse (http://soybase.org/gb2/gbrowse/gmax1.01) (Figure 1).
The alignment of the CcTA v2 transcriptome assembly with the soybean genome predicted 10009 intron spanning regions (ISRs), for a total of 5033 TACs. The alignments and primer sets can be viewed on the SoyBase genome browser at http://bit.ly/Cajca_ISR3, and are available for download at http://cajca.comparative-legumes.org/data/contrib/cajanus_cajan_v2_primers.txt.gz. A minimum of one and a maximum of 19 ISRs were designed against each of the matched soybean gene (varying based on the number of introns in a gene and the ability of the primer prediction software to identify low-copy ISR markers across the introns). The largest number of genes with ISR markers (1516) contained two markers.
In addition to the identification of ISR markers in pigeonpea, approximate mapping positions for a portion of these markers were predicted based on syntenic regions between pigeonpea and soybean genomes, using mapped BES-SSR loci in the pigeonpea genetic map as anchor points. In this context, 239 pigeonpea BES-SSR loci (Bohra et al., 2011a), mapped onto the genome sequence of soybean, were used for identifying the syntenic regions between the pigeonpea and soybean genomes. Of these 239 BES-SSR loci, 93 showed probable synteny in soybean chromosomes, and were used to identify putative linkage groups for 6284 pigeonpea ISR markers.
This method produced putative linkage group assignments for all 11 of the pigeonpea linkage groups. The strongest associations were for Gm10 and CcLG02, and for Gm14 and CcLG10; for these associations, there were 12 and 10 markers in syntenic regions, respectively; and putative placements for 888 and 738 ISRs on CcLG02 and CcLG10, respectively. Details of correspondences of other ISR markers between soybean chromosome and expected pigeonpea CcLGs are given in Supplemental Table 2.
Primer pairs were designed and synthesized for a total of 128 ISR markers. All these primer pairs were screened for amplification of DNA from two pigeonpea genotypes, namely ICP 28 and the popular variety Asha (ICPL 87119). This analysis identified a set of 116 markers (90.6%) with scorable amplicons. These 128 ISR markers corresponded to 10 pigeonpea linkage groups, generally distributed evenly—as assessed relative to the syntenic regions on soybean chromosomes. Screening of these 116 ISRs on eight pigeonpea genotypes including seven cultivated (ICPL 332, ICPL 99050, ICPB 2049, ICP 28, ICPL 20096, ICPL 87091, and ICPL 87119) and one wild (ICPW 94) showed length polymorphism (two to three alleles) (Figure 2) with 70 (54.6%) markers (Table 4 and Supplemental Table 3). The polymorphism information content (PIC) value for the polymorphic markers ranged from 0.19 to 1.00, with an average of 0.16.
With an objective of identification of SNPs and indels with ISR markers, as an example, sequence data were generated for eight pigeonpea genotypes with six ISR markers (TAC12160, TAC08538, TAC10495, TAC04354, TAC18853, and TAC07243). After trimming low-quality sequences at both ends, sequence data were aligned and analyzed for occurrence of SNPs and indels in the eight genotypes. While six markers showed 27 SNPs and two ISR markers (TAC08538 and TAC07243) showed two indels between the wild species, ICPW 94 and seven cultivated species (ICPL 332, ICPL 99050, ICPB 2049, ICP 28, ICPL 20096, ICPL 87091, and ICPL 87119) (Supplemental Table 4 and Supplemental Figure 1). Indels observed in the cases of TAC08538 and TAC07243 markers were confirmed on mutation detection enhancement (MDE) gel as mentioned above (Figure 2).
In pigeonpea, conventional breeding to date has proceeded without the support of molecular methods. Limited use of germplasm over the course of pigeonpea domestication has also resulted in a very narrow genetic base (Yang et al., 2006). As a result, pigeonpea genetic improvement programs have made comparatively little progress and hence have faced problems in addressing key constrains to crop production, including a range of abiotic stresses (e.g. drought, salinity, water-logging) and biotic stresses (e.g. Fusarium wilt, sterility mosaic disease, Helicoverpa armigera). Only during last five years some efforts have been made to develop genomic resources such as SSRs, ESTs, genetic, maps and transcriptome assemblies in pigeonpea (Varshney et al., 2009a, 2010). The available genetic maps, especially based on intra-specific mapping populations, do not have a good marker density (Bohra et al., 2011b; Gnanesh et al., 2011).
This study reports a comprehensive transcriptome assembly that is based on 131 million sequence reads coming from a range of tissues of more than 16 different pigeonpea genotypes. We employed two assembly programs in order to take advantage of the characteristics of the constituent sequences. All Illumina GA IIx reads were assembled into contigs with the ABySS program, and FLX/454 reads were assembled with the miraEST program. Contigs generated for Illumina GA IIx and FLX/454 reads were assembled together with Sanger ESTs using the miraEST program. Therefore, the developed assembly is a hybrid assembly. Hybrid assemblies are considered superior over pure assembly based on sequence data coming from one sequencing platform, as weaknesses from single sequencing platforms may be compensated by different characteristics of sequences from other platforms (Schatz et al., 2010; Garg et al., 2011).
The assembly described in this study is coming from 31 tissues representing a range of development and growth stages as well as challenged by different stresses. Thus, this assembly (CcTA v2) can be considered the most comprehensive transcriptome assembly of pigeonpea. The completeness and quality of this assembly can be assessed by comparing it with other earlier assemblies (Raju et al., 2010; Dubey et al., 2011; Dutta et al., 2011). The average TAC length as well as N50 of TAC in the CcTA v2 is much better than the CcTA v1. It is important to mention here that earlier assemblies were developed based on CAP3 (Dubey et al., 2011; Raju et al., 2010) and Lasergene SeqMan Pro™ v8.0.12 (Dutta et al., 2011) programs, while the present assembly has been developed using two powerful assembly programs that can accommodate large amounts of next-generation sequence: ABySS and miraEST.
The transcriptome assembly developed here can be used for a variety of applications to advance genetics research and breeding applications in pigeonpea. For instance, this assembly can provide the information about gene content and function, for identification of candidate genes, for development of molecular markers such as SSRs, SNPs, etc. (Varshney et al., 2009b). The majority of these applications have already been explored and a preview of such applications for the pigeonpea transcriptome was presented in our earlier studies (Dubey et al., 2011; Dutta et al., 2011). Therefore, these topics are not being discussed in this study. We have demonstrated one major application of the transcriptome assembly in the development of genome-wide marker datasets for enriching the genetic map of pigeonpea, using a comparative genomics approach that employs the soybean genome sequence and the BES-SSR loci-based genetic map of pigeonpea (Bohra et al., 2011a). For instance, comparison of the CcTA v2 with soybean genome identified the homologs for 77.5% of the pigeonpea TACs, and covering 27490 genes in the soybean genome. Of these, the majority of TACs (9863/16622) mapped twice against soybean genome—as expected, because of the ~13-Mya genome duplication in soybean. Occurrence of more than two hits in the soybean genome for a given pigeonpea TAC is also not surprising, as ~75% of soybean genes were found in multiple copies, due to older genome duplications in soybean’s history (Schmutz et al., 2010). Based on exon–intron boundaries, 10009 primer pairs, designed for 5033 TACs, can be used for checking the length or sequence polymorphism between the parental genotypes of mapping populations. To shortlist a set of markers that can be mapped across the pigeonpea genome, 177 anchor points were identified between the soybean and pigeonpea genomes, with the genetically mapped pigeonpea BES-SSR loci being placed by sequence homology onto the soybean genome sequence. Based on this information, 6284 ISR markers were identified that have putative chromosomal placements in the pigeonpea genome. A subset (128) of these markers was further analyzed for length (indel) polymorphism in eight parental genotypes of mapping populations segregating for three important traits, FW, SMD, and pod borer that are significant for pigeonpea improvement.
While 90.6% (116) markers provided scorable amplicons, 54.6% (70) markers showed polymorphism with two to three alleles in the genotypes analyzed on MDE gel. Utility of single-strand confirmation polymorphism-based SNPs and indels has been shown earlier in several species like common bean (Galeano et al., 2009), pearl millet (Thudi et al., 2010), etc. Although the majority of the markers showed polymorphism between ICPW 94 (C. scrabaeoides) and other genotype(s) of cultivated species (C. cajan), only three markers showed length polymorphism between the parental combinations of cultivated species. These results are not unexpected, as very low levels of polymorphism in some other cultivated species have been observed in several earlier studies (Yang et al., 2006; Varshney et al., 2009a; Raju et al., 2010; Varshney et al., 2010; Bohra et al., 2011a, 2011b; Dubey et al., 2011; Dutta et al., 2011; Gnanesh et al., 2011). On the other hand, sequence analysis of amplicons generated for six ISR markers showed their utility for identification of SNPs and indels at sequence level. Two common markers used for length as well as sequence polymorphism also confirmed the length polymorphism on MDE gel at sequence level. In brief, the validation results with some markers for detection of polymorphisms on MDE gel and sequence level underline the importance of developed resource of ISR markers. These markers should be useful for genetic mapping and trait mapping in breeding programs to develop the superior pigeonpea varieties with enhanced crop productivity.
In conclusion, the present study demonstrated a high-quality comprehensive transcriptome assembly of the important legume crop pigeonpea using Sanger and second-generation sequencing (FLX/454 and Illumina GA IIx) technologies. The results deliver novel information for future genetic studies in pigeonpea and provide a robust transcriptome assembly. The identification of syntenic regions between the pigeonpea and the sequence of a related phaseolid legume, soybean, provides greater insight into the gene content of pigeonpea. For the ISR markers identified, their putative mapping positions and parental polymorphism information will be a useful resource for molecular breeding programs to develop elite pigeonpea cultivars.
The following four datasets were used for defining the transcriptome assembly: (1) 128.9 million Illumina GA IIx short single end reads (1x36-nt) generated from 12 genotypes at ICRISAT and NCGR (Dubey et al., 2011; unpublished), referred to as Dataset I; (2) 1.696 million FLX/454 reads generated from the genotypes ‘Asha’ and ‘UPAS 120’ at NRCPB (Dutta et al., 2011), referred to as Dataset II; (3) 494353 FLX/454 reads generated from ‘PusaAgeti’ at ICRISAT and J. Craig Venter Institute (JCVI) (Dubey et al., 2011), referred to as Dataset III; and (4) 18353 vector-trimmed Sanger ESTs downloaded from dbEST (www.ncbi.nlm.nih.gov/dbEST/dbEST_summary.html) (the majority of which were generated at ICRISAT (Raju et al., 2010) and NRCPB, referred to as Dataset IV) (see Table 1).
Sequence datasets, as mentioned above, were assembled using the programs ABySS (Simpson et al., 2009) and miraEST (Chevreux et al., 2004), using the following three steps. In the first step, all Illumina GA IIx reads (Dataset I) from 12 genotypes were pooled and assembled together using ABySS. In the second step, FLX/454 reads from three genotypes (Datasets II and III) were trimmed of adapter sequences and assembled individually using the miraEST assembler. Subsequently, the pooled Illumina GA IIx (step 1 by ABySS) and FLX/454 (step 2 by miraEST) assemblies were merged with vector-trimmed Sanger ESTs of Dataset IV using the miraEST program. Both programs were run with the default settings, except for the following parameters: for ABySS, scaffolding ‘on’ at the paired end stage; and for miraEST: number of threads was seven and these options specified as ‘no’: Load straindata, Enforce presence of qualities, Extra gap penalty, and Wants quality file. In order to decrease runtime, number of processors used was seven. Since we were interested in a consensus assembly, the ‘Load straindata’ option was turned off. During the second stage of the assembly in which FLX/454 and Sanger ESTs were merged, there were no quality scores for the synthetic ESTs. Therefore, ‘Enforce presence of qualities’ and ‘Wants quality file’ options were specified to ‘no’. By turning off ‘Extra gap penalty’, we avoided penalizing gaps during the smith waterman alignment, especially since FLX/454 data are known to have homopolymer errors.
In order to check for microbial contamination, BLAST of all 21434 CcTA v2 contigs was carried out against NCBI’s microbial genomes database (www.ncbi.nlm.nih.gov/genomes/MICROBES/microbial_taxtree.html). The TACs of CcTA v2 having significant hits to bacterial genomes have been later on run against NCBI’s EST_others database (non-human, non-mouse) to check whether they could be mapped to other plant genomes.
All TACs of the CcTA v2 assembly were aligned to soybean genome (Schmutz et al., 2010) build 1.0.9 using Exonerate 2.2.0 (Slater and Birney, 2005), with parameters and flags ‘percent 25’ (to report only alignments over 25% of the maximum score attainable by each query) and ‘refine region’ (to perform an exhaustive alignment over the region in which the heuristic alignment was found). Alignments were filtered to require at least 80% alignment identity and 50% query coverage. If this resulted in more than 12 matches for a given sequence, the sequence was considered repetitive, and all matches were discarded.
All genetically mapped BES-SSR loci onto the pigeonpea genetic map (Bohra et al., 2011a) were anchored to the soybean genome using BLASTN (Altschul et al., 1997) of the corresponding BESs (with maximum E-value 1e-8), followed by manual selection for best hits with matches up to two homologous soybean regions.
Alignment results of pigeonpea TACs with the soybean genome were analyzed for identification of flanking intron junctions. The Exonerate alignment of the TACs, in Exonerate ‘vulgar’ (Verbose Useful Labelled Gapped Alignment Report) output format, was used to identify intron junctions in the TAC fasta sequences. These junctions were used to design the primer pairs using Primer3 (Rosen and Skaletsky, 2000) and BatchPrimer3 (You et al., 2008). Primer pairs were remapped to the soybean genome (to evaluate for repetitive sequences) using e-PCR (Schuler, 1997), with parameters ‘-n3 -g1 -t3 -m400 -d50-1000’. These parameters have the following effects: ‘-n3’ allows up to three mismatches per primer; ‘-g1’ allows up to one gap per primer; ‘-t3’ specifies output in tabular format; ‘-m400’ specifies an allowable margin for the product of 400 bases; and ‘-d50-1000’ specifies the default PCR product size range. Primer pairs with more than two alignments at these parameters were discarded.
Putative approximate mapping positions for the identified ISR markers were imputed based on anchoring points between pigeonpea and soybean genetic maps using BES-SSR loci of pigeonpea. Where there are two or more pigeonpea SSR markers with proximity in both pigeonpea and soybean (i.e. with nearby cM values in pigeonpea and nearby nucleotide positions in soybean chromosome pseudomolecules), tentative pigeonpea linkage groups (CcLGs) were assigned for ISR candidate markers occurring between the neighboring pigeonpea SSR markers.
Polymerase chain reactions (PCRs) for amplification of ISR loci were performed on eight pigeonpea genotypes (seven cultivated and one wild species) in a 5-μl reaction volume as described by Gujaria et al. (2011). Amplified products were denatured and separation was undertaken on MDE gel electrophoresis as described earlier (Thudi et al., 2010).
For detection of SNP or indel polymorphism in the case of ISR markers, PCR products for eight pigeonpea genotypes using six ISR markers were sequenced in both directions using Sanger sequencing methodology. Sequence data analysis and SNP identification among the selected genotypes were carried out as described in our earlier study (Gujaria et al., 2011).
Supplementary Data are available at Molecular Plant Online.
The authors thank the CGIAR Generation Challenge Programme (GCP), Mexico (G.D.M., R.K.V., and N.K.S.) and the Indian Council of Agricultural Research (ICAR), India (N.K.S. and R.K.V.) for sponsoring this research. No conflict of interest declared.
The authors are thankful to Anuja Dubey and Rachit Saxena for their help extended throughout this study.