Detection and validation of single feature polymorphisms in cowpea (Vigna unguiculata L. Walp) using a soybean genome array

Cross species platform for array hybridization

SFP detection and validation

Robustified projection pursuit (RPP) [20] was used for SFP detection between two parental genotypes (CB46 and IT93K-503-1) of a RIL population segregating for several important agronomic traits. As stated in Methods, RPP provides a list of SFP probe sets and single probes positioned directly over genetic polymorphisms. Because the RPP method depends on the availability of "present" calls in both of the comparison genotypes, and because only about 20% of the probe sets are called "present" using cowpea RNA, only about 7,000 of the 37,500 soybean probe sets have the potential to detect SFPs in cowpea from typical cowpea RNA samples. Using a top 15% outlying score cut-off, we generated a list of 1058 putative SFP probes between genotypes CB46 and IT93K-503-1. A full list of these SFPs and their outlying scores is provided in Additional file 1. The question then was, "Do these statistically detected SFPs represent real genetic polymorphisms in cowpea?"

A schematic diagram of SFP detection and validation is given in Figure 1. Plots of the log intensities, affinity differences and individual outlying scores for a representative probe set (Gma.1863.1.S1_at) under both drought-stress and non-stress conditions are shown in Figure 2. The intensity differentiation is highest at probe 2 between CB46 and IT93K-503-1, indicating polymorphism at this probe position [20]. Probe 2 was selected as the best SFP probe by the RPP method based on having the highest outlying score.

PCR and gel electrophoresis

A list of primers used and the expected and observed PCR amplicon sizes are provided in Table 2. PCR amplicon sizes ranged from 400 – 450 bp. A representative gel image of PCR amplification is shown in Figure 3. For the nine primer pairs (GS1 to GS9), the amplicon sizes were as predicted by Primer3 (see Methods), the primer designing software. For primer pair GS 10, non specific bands were detected in both genotypes. For this primer pair the 'T_m' was increased during PCR but still no specific band of the predicted size was amplified. Primer pairs yielding such anomalies, or which yielded no amplicon from one or both cowpea genotypes (not shown), were not pursued with amplicon sequencing. For calculation of the validation rate, the numerator is the number of amplicons containing a polymorphism within the region spanned by a SFP probe, and the denominator is the number of amplicons from which sequences were generated.

Alignment of amplicon sequences for SFP validation

To estimate the SFP validation rate we selected 25 putative SFP probe sets. A total of 17 (68%) of these 25 SFP probe sets were validated by amplicon sequences. A representative alignment of genomic amplicon sequences with the target sequence of probe set Gma.1863.1.S1_at is shown in Figure 4. Probe position 2 spans a polymorphism (a SNP in this case), as predicted by the RPP method. Alignments of the 20 validated SFP probes from the 17 validated probe sets (three probe sets each contained two SFP probes) are shown in Figure 5. Among these 20 validated SFP probes, 15 SFP (75%) were positioned over a single SNP, 2 (10%) were positioned over only a single nucleotide insertion or deletion (INDEL), 1 (5%) spanned one SNP and one single nucleotide INDEL, 1 (5%) spanned two SNPs, and 1 (5%) spanned two SNPs and a dinucleotide INDEL.

Array-based genotyping

The availability of arrays designed to measure gene expression for a wide range of species has spawned considerable interest in identifying single feature polymorphisms (SFPs) from DNA hybridization and transcriptome data [27]. SFPs derived from transcript sequences provide a link to structural genes and are particularly attractive for species which lack sequence resources for marker development. The RPP method [20] is distinct from several commonly used statistical approaches to SFP identification [27] in that the RPP method does not accept expression level polymorphisms, which may be attributed to linked cis-acting regulatory regions or to unlinked trans-acting factors. Instead, RPP limits its SFP calls to those that are most likely to be positioned over a genetic polymorphism. This maximizes the relationship between the probe set functional annotation and the gene containing the detected polymorphism, which consequently maximizes the potential for cross-species synteny mapping using SFPs. A simplification of the RPP method has been used for 16,000 Mb hexaploid wheat to map markers termed "high variance probe sets" to delimit the positions of translocation breakpoints using wheat-rice synteny [28]. Here we report cross species application of the Affymetrix soybean genome array by RPP for polymorphism identification in parents of a cowpea recombinant inbred line mapping population.

The false positive rate in SFP detection is an important consideration for application of SFPs to genetic mapping. The false positive rate is essentially the percentage of statistically defined SFPs that are not valid at the DNA sequence level. Ideally, the false positive rate would be 100% minus the validation rate from amplicon sequencing. However, for cowpea and other species which lack complete genome sequence information, incomplete knowledge of genes, and gene families in particular, limits the validation rate to a theoretical maximum less than 100%. This is because PCR primers for amplicon production must be designed only from known sequences and therefore cannot always target the same gene detected by a SFP probe. Amplification of the wrong member of a gene family bearing no polymorphism can result in a false non-validation of a SFP (false negative). We had access to sequences of less than 70% of all cowpea expressed genes, and therefore also to less than 70% of the members of any given multigene family. This absence of sequence information would tend to cause an overestimation of the false positive rate. However, the selection of SFP probe sets for validation was somewhat skewed toward SFPs with higher outlying scores, which may cause an overestimation of the validation rate. Therefore, the estimated 68% validation rate and inferred 32% false positive are only rough estimates. We have observed that at least 50% of the non-validated cowpea SFPs are in genes that are multi-gene family members. An extreme case of multigene families comes from a study of Medicago truncatula nodules, which revealed the NCR (nodule-specific cysteine rich) gene family with more than 300 members [29].

A false positive rate of about 5% in SFP detection was reported in yeast [17] and 3% in Arabidopsis [18] using genomic DNA. In contrast, approximately a 25% false positive rate was reported in 440 Mb rice using genomic DNA [19], and 10–20% in 5300 Mb barley using RNA-based datasets and RPP [20].

Cross-species platform for orphan crops/organisms

Recently a number of interspecies comparisons of gene expression have been carried out including human versus monkeys [30, 31], between rodents [32], human versus mouse [33], within Xenopus [34] and within Drosophila [35]. Cross-species analysis of gene expression in non-model mammals was reported by Nieto-Díaz et al. [36]. The reproducibility of probe data obtained from hybridizing deer, Old-World primates, and human RNA samples to the Affymetrix human GeneChip^® U133 Plus 2.0 was compared. Cross-species hybridization affected neither the distribution of the hybridization reproducibility among different categories nor the reproducibility values of the individual probes. Studies such as these encouraged us to extend the cross-species concept to SFP detection for cowpea using the readily available Affymetrix soybean genome array.

In the SFP validation process we identified polymorphic sites including SNPs and INDELs in sequences neighbouring the SFP probe positions. Such polymorphisms, which generally represent haplotypes, also can be used for mapping and therefore are useful by-products of genomic amplicon sequencing for SFP validation. As shown in Figure 6, there are even instances of amplicon sequence polymorphisms when no polymorphism is detected in an SFP position. In general cowpea seems to have abundant genetic polymorphism.

Cost of array-based genotyping

The costs of marker development and application are important to consider. At our cost of $600 per soybean genome array for purchase, labelling, hybridization and data collection, the array-related costs of genotyping 100 RILs from two cowpea genotypes using one array per RIL would be $60,000 USD. If this would yield 1000 mapped SFP markers, then this would amount to $60 USD per mapped marker. Depending on the circumstances, this may or may not be attractive relative to other marker options, and array costs may vary at different locations. Nevertheless, $60 USD per mapped SFP in a 100-RIL population may be useful as an approximate point of reference. Alternatively, using RPP one could readily generate a short list of SFP-bearing probe sets using datasets derived just from parental genotypes and then develop smaller, less expensive arrays containing only the probe sets that contain SFPs, reducing the cost per mapped SFP accordingly.

Another possible application of cross-species SFP analysis would be simply to compare any two genotypes within a species in order to focus amplicon sequence-based marker development on a set of target sequences that have a high chance of revealing a polymorphism. A modest expenditure to produce two datasets from each genotype of interest could provide a sizeable cost savings in materials and labour by increasing the success rate of amplicon sequencing in search of polymorphisms.

Plant materials

The genotypes used in this study were California Blackeye No. 46 (CB46) and IT93K-503-1. CB46 was developed at the University of California, Davis, released in 1989 [37] and is currently the most widely grown blackeye-type cowpea cultivar in the United States. IT93K-503-1 is a breeding line developed by the International Institute of Tropical Agriculture (IITA) in Ibadan, Nigeria to produce high grain yield in the savanna agro-ecological zone of West Africa. The parental genotypes are inbred pure lines and differ in a number of important agronomic traits including grain size, photoperiod sensitivity, seedling drought tolerance, and resistance to Fusarium wilt race 4 and brown blotch fungal diseases. In order to determine the inheritance of these traits and to develop nearby markers, a RIL population with 135 lines was created.

Soybean genome array

Phylogenetic relationships based on the conserved sequences within Papilionoideae legumes imply that Vigna (cowpea) is closely related to soybean [38]. Since a cowpea genome array was not available, a soybean genome array was used to identify SFPs in cowpea. The soybean genome array contains 37,500 probe sets derived from soybean (Glycine max L.) unigenes. This represents 61% of the total probe sets on the chip, with the remainder targeting two pathogens important for soybean research, of which 15,800 (26%) probe sets target Phytophthora sojae (a water mold) and 7,500 (12%) probe sets target Heterodera glycines (soybean cyst nematode). This array uses probe sets composed of 11 probe pairs to measure the expression of each gene. Each probe pair consists of a perfect match (PM) probe and a mismatch probe.

RNA purification

Seeds were germinated in two sets of pots and grown in controlled glasshouse conditions during spring 2005 under natural photoperiod. The temperature was maintained at day/night cycle of 35/18°C. The growing axillary buds were harvested 14 days after planting (DAP) from one set of pots to provide tissues for control samples (flash frozen and stored at -80°C). Plants in the second set of pots were exposed to drought stress induced by withholding water after germination. The growing axillary buds were harvested from drought stressed seedlings at 14 DAP (flash frozen and stored at -80°C). With two genotypes, two treatments and two replicates, there were a total of 8 samples. Total RNA was isolated using TRIzol (Gibco BRL Life Technologies, Rockville, MD) reagent. RNA was further purified using an RNeasy spin column (Qiagen, Chatsworth, CA) and an on-column DNase treatment. RNA integrity was assessed prior to target preparation using RNA Lab-On-A-Chip (Caliper Technologies Corp., Mountain View, CA) evaluated on an Agilent Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA).

Labelling and hybridization of cowpea RNA on soybean genome array

Cowpea RNA samples were used to make biotin tagged cRNAs. These were hybridized to an Affymetrix soybean genome array as recommended by Affymetrix (Affymetrix GeneChip Expression Analysis Technical Manual; Affymetrix Inc., Santa Clara, CA) at the Institute for Integrative Genome Biology Microarray Facility at the University of California, Riverside. The hybridization data were scanned for visible defects and then extracted using default settings and tabulated as CEL files using Affymetrix GeneChip Operating Software (GCOS 1.2). A global scaling factor of 500, a normalization value of 1, and default parameter settings were used for the soybean genome array. The detection calls (present, absent, or marginal) for the probe sets were made by GCOS as follows (abridged from the Affymetrix GeneChip Expression Analysis Technical Manual; Affymetrix Inc., Santa Clara, CA). The detection algorithm uses probe pair intensities to generate a detection p-value and assign a "present", "marginal", or "absent" call. Each probe pair in a probe set has a potential vote in determining whether the measured transcript is or is not "present". The vote is described by the discrimination score (R). R is calculated for each probe pair and compared to a predefined threshold, Tau. Probe pairs with R higher than Tau vote "present" and the voting result is summarized as a p-value. The greater the number of discrimination scores (R) that are above Tau, the smaller the p-value and the more likely the given transcript is truly present in the sample.

Method for identifying SFPs

Expression data were generated by hybridizing cowpea cRNA to the soybean genome array. A statistical method called robustified projection pursuit (RPP) was used for SFP analysis [20]. Only the values from the PM probes were utilized. The use of RNA as a surrogate for genomic DNA eliminated interference from highly repetitive DNA as a technical impediment to SFP detection. An important aspect of the RPP method is that it first utilizes a probe set level analysis to identify SFP-containing probe sets and then chooses individual probes from within each SFP-containing probe set. The net result is the identification of probes that directly overlay polymorphic sequences.

Separate comparisons were made between two genotypes (with two replicates each) for unstressed and drought stressed treatments, resulting in two SFP lists. In the context of SFPs, there is no necessity to have separate stress and control lists; in fact it would be simpler and less costly to have only one SFP list from highly complex RNA made by blending stressed and unstressed RNA, as in Cui et al. [20]. In our case, two separate lists were available as a consequence of another study not described here which compared gene expression patterns in stressed and control plants (data not shown). At 15% outlying score cut-off, we detected 488 SFP probes in stressed and 661 SFP probes in unstressed treatments. The union of these two lists contained 1058 SFP probes and the intersection contained 91. A total of 37 primer pairs targeting 37 putative SFP probe sets were initially tested, of which 25 yielded single amplicons of the expected sizes from both parents. These 25 amplicons targeted 14 probe sets selected from the intersection of the two SFP probe set lists and 11 from the remaining SFP probe sets. As shown in Additional file 1, 9 of the 14 SFP probe sets (64%) from the intersection list were validated at the DNA sequence level and 8 of the other 11 (73%) were validated.

Genomic DNA isolation

Young leaves from one-week-old cowpea seedlings were ground with liquid nitrogen. Approximately ~0.1 g of ground tissue was used for genomic DNA isolation using a DNeasy Plant Mini Kit (Qiagen, USA) according to the manufacturer's recommendations. Eluted genomic DNA was examined by UV absorbance and with 1.0% agarose gel electrophoresis for quality and quantity evaluation.

Primer design and PCR

The soybean genome array unigene sequences were used to query (using blastx) Arabidopsis translated gene models (version 7.0) from The Arabidopsis Information Resource (TAIR). Annotations for the Affymetrix soybean probe sets were compiled into a browser called HarvEST:SoyChip which can be accessed online or downloaded for Windows installation [39]. Cowpea ESTs available from GenBank were assembled using the CAP3 program [40] and compiled into a browser called HarvEST:Cowpea which can be accessed online or downloaded for Windows installation [41]. Cowpea methyl-filtered sequences were obtained from Dr. Michael Timko, Department of Biology, University of Virginia [42]. The soybean unigenes corresponding to each SFP probe set were used to query cowpea unigenes in HarvEST:Cowpea and cowpea methyl-filtered sequences. The cowpea sequences corresponding to SFPs were then used for PCR primer design in the flanking regions of SFP position. Primers were designed using Primer3 web version software [43]. The strategy used for primer design is illustrated in Additional file 2. All primer pairs were designed with T_m~55°C.

PCR was performed in 20 μl reactions containing 20~25 ng of genomic DNA, 0.1 μM of primers, 0.2 mM dNTPs, and 1 unit of Taq DNA polymerase (New England Biolabs, USA). The reaction included an initial 5 min denaturation at 95°C, followed by 35 cycles of PCR (94°C, 30 sec; 55°C, 45 sec; 72°C, 1 min), and a final 5 min extension at 72°C. Aliquots (4 μl) of the PCR products were loaded and separated on 1.2% agarose gel by electrophoresis. A higher T_m was used when non-specific bands were amplified. PCR products were purified using QIAquick PCR purification Kit (Qiagen, USA) after confirming their uniformity on agarose gels.

DNA sequencing and analysis

DNA sequencing was performed using the dideoxynucleotide chain termination method [44]. Both strands of the amplified PCR products were sequenced with an ABI-PRISM 3730xl Autosequencer (Aplied Biosystems, USA) at the Core Instrumentation Facility of the UC Riverside Institute for Integrative Genome Biology. The forward and reverse sequence reads for each genotype for a particular primer pair were assembled using Contig Express (Invitrogen, USA) and the consensus sequence was used for alignment. These sequences were then compared with each soybean Sequence Information File (SIF, the target sequence which extends from the 5' end of the 5'-most probe to the 3' end of the 3'-most probe) and cowpea sequence using AlignX (Invitrogen, USA). Comparisons of nucleotide sequence similarity were displayed using GeneDoc [45]. Cowpea genomic amplicon sequences have been deposited into the dbGSS Data Library [GenBank: ET041523 to ET041556].

Data Availability

All expression data are available through the Gene Expression Omnibus (GEO) under platform GPL 4592, Series GSE 10284.