Development of ESTs from chickpea roots and their use in diversity analysis of the Cicergenus

Analysis of ESTs and assembly of consensus sequences

Gene annotations were based on similarities to either known or putative ESTs in the public databases. All annotations were based on Blast searches, with a score threshold of ≥ 200 for BLASTn. For tBLASTx a score threshold of >100 was set, as these generally had e-values <10^-5 with a minimum of 50% identity over at least 30% of the length of the protein, which are the commonly used thresholds for reliable sequence annotation [39, 40]. Tentative functional annotations of EST sequences based on tBLASTx were grouped under 12 general categories based on the biochemical functions of the predicted proteins (Figure 1). The 2858 sequences were assembled into 210 contigs with 267 ESTs remaining unassociated (singletons). Putative identifications could be assigned to 2179 of the 2858 sequences (76%), representing 73% of the unigenes (178 singletons and 162 contigs), and all of these proteins were represented in other plant species. The 2858 sequences were assembled into 210 contigs with 267 ESTs remaining unassociated (singletons). Assuming that in most cases each contig represents one gene [41] a maximum number of 477 genes are represented in this dataset (available at NCBI: accession numbers CK148643–CK149150). The total length of contigs ranged from 240 to 1291 bases (with an average of 648 bases) while the largest number of ESTs assembled into one contig was 255. EST assembly leading to the generation of tentative consensus (TC) sequences has dramatically reduced the redundancy in this database, as reported elsewhere [42] (see Additional file 1).

Over three quarters of the 210 consensus sequences (77% equivalent to 164 contigs) were found to have significant similarity (tBLASTx >100) to sequences in public databases. Many of these annotations (47 consensus sequences) were validated through comparison with the CDD database using RPS-BLAST which compares a protein query sequence to a position-specific score matrix prepared from the underlying conserved protein domain alignment [43]. Conserved motifs consistent with tBLASTx annotations were identified in 47 consensus sequences providing an additional level of confidence for these assigned gene functions (see Additional file 2).

As expected, we observed a substantial number of the EST contigs and singletons to have no significant homology (NSH: 48 contigs and 89 singletons; 29%) to either nucleotide or protein sequences in public databases at the time of analysis. This compares well to the higher proportion of sequences (44%) reported with little or no homology from chickpea leaf tissue [44]. In addition, a substantial number of ESTs (approximately 19% of the unigenes) were found to have similarity with unannotated hypothetical, putative or unknown proteins contained in the public databases and were, therefore, placed in the unidentified function (UF) class.

PCR optimisation and unit costs

Optimization of PCR conditions is a critical precursor for accurate and robust large-scale marker screening and offers additional benefits for reducing the unit cost of genotyping [45]. As optimization relies on the sequential investigation of each reaction variable, the process can take a long time, thus the conditions of optimization are rarely identified empirically and individual components are rarely all tested simultaneously. However, Cobb and Clarkson [46] devised a strategy in which the components of a PCR reaction could be tested together using the least number ofexperiments based on the Taguchi principle [47]. We adapted the Cobb and Clarkson [46], 9 reactions of 4 PCR components to a 5 reactions of 5 PCR components, the additional component being enzyme concentration. This further reduced the time and cost of optimization, whilst also minimizing the cost of the PCR screening by also considering enzyme concentration. Final primer (0.15 - 0.5 pmol) and Taq DNA polymerase (0.2 – 0.5 U) concentrations resulted in a reagent unit cost of US$0.08 per PCR sample. This compares very favorable with other reports of PCR unit costs [48, 49]. Clearly DNA extraction costs are an equally important component for molecular breeding programs and can also be greatly reduced as has been reported elsewhere [50]. The primer sequences and optimization parameters of all the polymorphic markers are given in Table 2. Optimization also resulted in fewer spurious amplification products and clearer polymorphic bands amongst the accessions studied (Figure 2).

EST marker diversity analysis

Fourteen of the 106 EST's contained SSR motifs, of which 10 contained perfect trinucleotide repeats. A relatively large proportion of these EST-SSR markers (8) were derived from the no significant homology (NSH) and unknown protein functional (UF) classes, based on sequence comparisons (BLASTn tBLASTx and PSI-BLAST). Markers from these EST classes have readily detected polymorphism in the diverse germplasm tested (PIC 0.20–0.82). Since over 48% of the EST's generated from the chickpea root library fall into NSH and UF classes, these are clearly useful markers with potentially important, yet unknown functions.

Eleven of the polymorphic EST's were from stress annotated transcripts (AGLC2, AGLC16, AGLC 20, AGLC29, AGLC 34, AGLC45, AGLC52, AGLC53, AGLC55, AGLC66, AGLC68 and AGLC93); of these four contained SSR motifs (AGLC 34, AGLC55, AGLC67 and AGLC84). These stress related EST markers may be useful for allele-mining of germplasm collections for the identification of candidate accessions with new sources of agronomically important traits and for candidate gene mapping [51].

The EST data-set generated in this study was used to compute pair-wise genetic distances between species according to Band similarity coefficient and UPGMA clustering. In order to display, with minimal distortion, the genetic relationships between species, a Principle Coordinate Analysis (PCoA) was carried out. The first two dimensions of the PCoA plot indicate the presence of 5 clusters, which account for 20% and 17% of the total variation (Figure 3). Dollo and polymorphism parsimony analysis (DOLLOP) conducted on this data-set found one most parsimonious tree with 301 steps (Figure 4) corresponding well with the PCoA analysis.

Both PCoA and parsimony analysis provide clear separation of most species. Both analysis grouped C. echinospermum, C. reticulatum and C. arietinum, together. These three species are classified together as the primary gene pool based on crossability studies. Similarly, all previous diversity studies based on allozymes [52], RAPD [16, 53], microsatellite [54], ISSR [15, 55] and AFLP [17] markers have also grouped C. echinospermum, C. reticulatum and C. arietinum together. The EST-based diversity analysis presented here shows C. arietinum (believed to be the wild annual progenitor [56] of cultivated chickpea) to be slightly distant to C. echinospermum and C. reticulatum but still within the same cluster (best reflected by the PCoA analysis Figure 3). All previous molecular diversity studies have clustered C. bijugum, C. pinnatifidum and C. judaicum together, although sometimes not well differentiated from the primary gene-pool species based on molecular diversity analysis. These three species are commonly considered to represent the secondary gene-pool as they require embryo rescue to recover viable progeny when crossed with the cultigen [6]. Both parsimony and PCoA analyses of the EST data-set clustered C. judaicum, C. pinnatifidum and C. bijugum together. This is in agreement with previous studies which demonstrated the close association of these three species. The current analysis suggests that C. pinnatifidum and C. judaicum are more closely related. This is in agreement with previous reports based on morphology analysis [56], seed storage protein analysis [57] and isozyme analysis [12]. In contrast, ecogeographical studies of the wild Cicer germplasm have suggested that C. judaicum to be quite diverged [7], which may explain the outlying position of one C. judaicum accession studied here (IG 69986) (Figure 4). Most other species generally fall into one of two clusters generally collectively referred to as the tertiary gene pool (Figure 3 and Figure 4) from which no viable hybrids have been reported from crosses with the cultigen.

Parsimony analysis grouped C. chorassanicum with C. cuneatum, but placed the C. yamashitae accessions as outliers (Figure 4). This data supports the tertiary group proposed by Croser et al. [6] to consist of; C. chorassanicum which has never been successfully crossed with cultivated chickpea and C. yamashitae and C. cuneatum which have been crossed with the cultigen but the resultant progeny have failed to flower or proven sterile. However, this is not reflected in the PCoA (Figure 3).

The distribution of accessions from the same species within secondary and tertiary clusters suggests that these two gene-pools may not be as distinct as previously expected. This may indicate that there are some tertiary gene-pool species that might be more readily crossed with the cultigen than their taxonomic classification would otherwise indicate. However, since the current study is based on a limited numbers of accessions for both C. cuneatum and C. yamashitae further research is required to fully define the relationships between these species and the scope of new opportunities for plant breeders.

Plant material and growth conditions for SSH

Fifteen seeds from each chickpea accession (ICC 4958 and Annigeri) were sterilized with 25% chlorex (v/v) for 10 mins, rinsed twice with sterile distilled water, before placing in pre-sterilized plant pots containing sterilized vertisol. The soil was collected from the field, passed through a 2 mm sieve and autoclaved in plastic bags using a 60 min cycle at 120°C and repeated three times on consecutive days. Seeds were sown immediately after the last autoclave cycle.

Three seeds were sown per pot; the plants were maintained in a Conviron growth chamber (Conviron, Winnipeg, Man.) under optimal physiological conditions for chickpea growth: 25°C day (11 h) and 12°C night (13 h), 455 μmoles m^-2 sec^-1 illumination intensity, 80% day and 40% night relative humidity, and irrigation as required.

Roots were harvested 28 days after germination (with flowering expected after 30 days). The soil was washed off the roots in running water, the intact whole plant was lifted out of the soil and soil particles were gently rubbed off in running water with gloved hands. The root tissue was subsequently treated with (0.1% v/v) DEPC solution and immersed in sterile Milli Q water, excess water from the roots was blotted onto filter paper and roots were weighed (approx. 2 g), flash frozen in liquid nitrogen and stored at -80°C.

RNA isolation and suppressed subtractive hybridization

The subtraction of the two genomes was achieved through a process of PCR-based suppression of genes common to both genotypes, and genes that are differentially expressed are enriched. As the process is based on PCR, low copy number cDNA can be detected from the tester (ICC 4958). RNA isolation and library construction was carried out by Avestha Gengraine Technologies (Bangalore, India). Total RNA was isolated from the root tissue of ICC 4958 and Annigeri using Trizol reagents (GIBCO-BRL) for RNA isolation, mRNA was purified using Oligotex (Qiagen GmbH, Hilden, Germany).

The suppression subtractive hybridization (SSH) process was carried out using a Clontech PCR-Select cDNA subtraction kit (Clontech INC. USA). Tester (ICC 4958) and driver (Annigeri) cDNA was restricted with RsaI for 90 mins and cDNA from the ICC 4958 was ligated with the required adaptors. Adaptor ligation was confirmed by PCR, followed by one round of over-night hybridisation with Annigeri as the driver and ICC 4958 as the tester. A second round of hybridisation was followed by primary PCR (27 cycles) and secondary PCR (12 cycles) according to the manufacturer's instructions for the PCR-Select. The PCR products of the subtraction were analyzed by gel electrophoresis and cloned using pCR4-TOPO cloning kit (Invitrogen Inc., USA). The cDNA clones were amplified in E. coli DH10B cells by electroporation and insert containing clones were selected by plating on LB agar containing 100 ug/ml ampicillin. Large colony forms were picked and used to regenerate single clone cultures in 96-well microtitre plates. After growth over-night at 37°C, glycerol was added to a final concentration of 15% and cultures were stored at -80°C.

EST sequencing

Cultures of transformed E. coli were grown over-night in 100 ul LB media containing 100 ug/ml ampicillin, the cells were centrifuged and suspended in sterile distilled water, heat denatured for 10 mins at 95°C, 1 ul of supernatant was used for insert amplification with M13 forward and reverse primers. Inserts (5 ul) were separated by electrophoresis on 1.2% agarose gel containing ethidium bromide to confirm amplification of a single insert and for quantification. After SAP/exonuclease I digest, 4 ul of the insert was sequenced using T₇ primer and 1/8 fraction of the recommended volume of dye terminator reagent (Big Dye terminator cycle sequencing kit v2, ABI-Foster City) for cycle sequencing. After removal of dye blobs by ethanol precipitation, DNA sequences were analyzed by capillary-electrophoresis using an Applied Biosystems genetic analyzer model ABI 3700 or ABI 3100.

Sequence analysis for predicting gene function

Base calling was performed using the ABI DNA Sequence Analysis Software (v3.7) or by Chromas v2.2 (Technelysium Pty Ltd. Australia). All sequences were scanned visually for quality of the peak shapes and corresponding base call. Vector sequences were trimmed using Sequencher v.4 (Gene Codes, Ann Arbor, MI) software and low quality bases (quality score <20) were trimmed from both ends of sequences. A total of 4000 sequences were analyzed, only sequences of more than 170 bp with less than 5% ambiguity were processed. All sequences containing interspersed or simple repeats were masked with the RepeatMasker software [61] using Arabidoposis as reference. The masked sequences were screened against the following DNA databases: Arabidopsis thaliana, Medicago truncatula, Glycine max, Zea mays and Oryza sativa, in the TIGR Gene Indices using default settings for BLASTn and tBLASTx searches [62]. The same sequences were also used to search the NCBI EST database (dbEST) using Viridiplantae as a limiter [63].

The comparison of chickpea sequences with those in public databases was carried out during 2002 – 2003, and blast searches with unique sequences were repeated in November 2004. An in-house script was used to retrieve best matches from the outputs of BLASTn and tBLASTx text files, retaining the following information: gene identifier, description, BLAST score, percent identity, alignment length, reading frame, and database name from the top-scoring alignments. Alignments with scores of ≥200 for BLASTn and ≥100 for tBLASTx were extracted into an MS-Access table. The individual ESTs were assembled into tentative consensus sequences (TCs) using Sequencher, the assembly parameters used for 'dirty data' were (minimum percentage match of 90 and minimum overlap 40 bp). A sequence was retained in a contig if it matched at least one sequence already in the contig based on pair-wise BLAST [63], with a percentage identity threshold of 90%. A consensus sequence was then derived for each cluster. Consensus sequences were translated into six possible reading frames with Transeq software [64], using the standard genetic code. The translated sequences were then submitted to Reverse Position Specific (RPS)-BLAST [63] to search for matches to the conserved domain database (CDD; all 4540-PSSMs) using default parameters. The RPS-BLAST outputs were similarly processed using an in-house script, retaining the following information from the top-scoring alignments: gene identifier, description, score and e-value. Each result was manually interpreted and classified according to the biochemical function of the predicted protein.

Plant material and DNA isolation

Primer design of EST markers

A total of 106 EST sequences were chosen that represented diverse functional/stress annotations, a selection of which contained SSR motifs. In addition, primers were also designed from a random selection of sequences categorized as having 'no significant homology' (NSH) and 'unknown protein function' (UF) when compared with either nucleotide or protein sequences in public databases.

Primer sequences were designed using Genefisher 1.1 software [65] using the following criteria 57 – 63°C melting temperature (Tm), 40–50% GC content, 17–24 bp primer length and the difference in Tm between forward and reverse primer pairs was limited to 2°C. For those EST sequences containing repeat motifs, primers were designed to flank the microsatellite region. Optimal primer pairs were selected using Netprimer, the primer design algorithm which analyzes all possible secondary structures including hairpins, self-dimers and cross-dimers [66]. The Tm and product length of 8 forward and reverse pairs per sequence from Genefisher software were evaluated by Netprimer software for minimal or no secondary structure. All primer and EST sequences are accessible through the ICRISAT Chickpea EST Database [38].

PCR optimization and evaluation

PCRs were performed in a total volume of 5 μl containing 5–10 ng of DNA and optimized concentration of the following: primer (0.15–0.3 pmoles), dNTP (0.1–0.2 mM), MgCl₂ (1.5–2 mM), Qiagen Taq (0.2–0.5 u) and 1 x buffer. PCR amplifications were carried out using a GeneAmp model 9700 thermocycler (Perkin Elmer-Applied Biosystems). As the annealing temperatures of the markers ranged from 47–65°C we opted to use 3 catagories of "touchdown" temperature cycles (see Table 4).

PCR amplified products with a size range of 200–600 bp were resolved by non-denaturing polyacrylamide gels electrophoresis (6–10% acrylamide/bisacrylamide, 20:1 in TBE). Bands were visualized through a modified silver staining protocol [67] as follows: gels were immersed in water for 3 mins, followed by 20 mins in 0.1% CTAB solution and 0.3% ammonia solution for 15 mins. Silver staining solution was freshly prepared each day, consisting of 0.1 % (w/v) AgNO₃ in 4 mM NaOH solution, 0.5 to 0.6 ml of 25% ammonia was titrated until the cloudy suspension became clear. Gels were gently agitated in the silver nitrate solution for 30 mins, and developed in 1.5% (w/v) sodium carbonate and 0.02% (v/v) formamide solution until bands appeared. The gels were rinsed in water, fixed in 1.5% glycerol solution and documented by scanning. PCR products above 600 bp were resolved by agarose gels electrophoresis (0.8–1%) containing ethidium bromide.

EST marker data collection and analysis

The amplification profile from each EST marker was checked for reproducibility and subsequently scored visually (and independently rescored) for the presence (1) or absence (0) of polymorphic bands. The degree of polymorphism was quantified using the polymorphic information content (PIC) calculator [68] based on the following formula:

PIC = 1 - ∑ Pi² - ∑ ∑ Pi²Pj²

i = 1 i = 1 j = i+1

where Pi is the frequency of an individual genotype.

Pair-wise genetic distances were calculated with NTSYS-pc software version 2.0 [69] as S_xy = 1 - (2n_xy/n_x + n_y) as S_xy = 1 - (2n_xy/n_x + n_y) based on band coefficient [70], where n_xy are shared bands amongst the individuals, n_x and n_y are the number of fragments exhibited by each individual. Distance based on Band coefficient is similar to the genetic distance (GD) formula by Nei and Li [71]: GD = 1 - S_xy. The 15 × 15 similarity matrix (based on Band coefficient) was subjected to sequential agglomerative hierarchical nested (SAHN) clustering using UPGMA (unweighted pair-group method analysis). Principle coordinate analysis based on these distance estimates was performed using NTSYS, using consecutive commands of 'DCENTRE' and 'EIGEN' in order to generate a PCoA plot which is more informative in highlighting differences between major taxonomic groups than dendrogram representations. Maximum parsimony analysis was conducted using the PHYLIP program Dollo and polymorphism parsimony analysis version 3.63 [72] and the consensus tree viewed in tree view.