Development and application of a 6.5 million feature Affymetrix Genechip® for massively parallel discovery of single position polymorphisms in lettuce (Lactuca spp.)
- Kevin Stoffel1,
- Hans van Leeuwen†1, 6,
- Alexander Kozik†2,
- David Caldwell†1, 7,
- Hamid Ashrafi†1,
- Xinping Cui4, 5,
- Xiaoping Tan1,
- Theresa Hill1,
- Sebastian Reyes-Chin-Wo1,
- Maria-Jose Truco2,
- Richard W Michelmore2, 3Email author and
- Allen Van Deynze1, 3Email author
© Stoffel et al.; licensee BioMed Central Ltd. 2012
Received: 14 November 2011
Accepted: 27 February 2012
Published: 14 May 2012
High-resolution genetic maps are needed in many crops to help characterize the genetic diversity that determines agriculturally important traits. Hybridization to microarrays to detect single feature polymorphisms is a powerful technique for marker discovery and genotyping because of its highly parallel nature. However, microarrays designed for gene expression analysis rarely provide sufficient gene coverage for optimal detection of nucleotide polymorphisms, which limits utility in species with low rates of polymorphism such as lettuce (Lactuca sativa).
We developed a 6.5 million feature Affymetrix GeneChip® for efficient polymorphism discovery and genotyping, as well as for analysis of gene expression in lettuce. Probes on the microarray were designed from 26,809 unigenes from cultivated lettuce and an additional 8,819 unigenes from four related species (L. serriola, L. saligna, L. virosa and L. perennis). Where possible, probes were tiled with a 2 bp stagger, alternating on each DNA strand; providing an average of 187 probes covering approximately 600 bp for each of over 35,000 unigenes; resulting in up to 13 fold redundancy in coverage per nucleotide. We developed protocols for hybridization of genomic DNA to the GeneChip® and refined custom algorithms that utilized coverage from multiple, high quality probes to detect single position polymorphisms in 2 bp sliding windows across each unigene. This allowed us to detect greater than 18,000 polymorphisms between the parental lines of our core mapping population, as well as numerous polymorphisms between cultivated lettuce and wild species in the lettuce genepool. Using marker data from our diversity panel comprised of 52 accessions from the five species listed above, we were able to separate accessions by species using both phylogenetic and principal component analyses. Additionally, we estimated the diversity between different types of cultivated lettuce and distinguished morphological types.
By hybridizing genomic DNA to a custom oligonucleotide array designed for maximum gene coverage, we were able to identify polymorphisms using two approaches for pair-wise comparisons, as well as a highly parallel method that compared all 52 genotypes simultaneously.
Various types of microarrays have been used extensively for gene expression studies and, more recently, for genotyping and marker discovery [1–5]. Affymetrix microarrays in particular offer a massively-parallel approach to genotyping. The basis of identifying polymorphisms, termed single feature polymorphisms (SFPs), is differential hybridization of template RNA or DNA onto 25 bp oligonucleotide probes on the array due to the presence of single nucleotide polymorphisms (SNPs) or small insertion/deletions (InDels). Using this approach, thousands of genes can be queried and simultaneously analyzed allowing whole genome approaches to mapping genes and quantitative trait loci (QTL) discovery , as well as determining linkage disequilibrium (LD)  and population structure [8, 9]. When arrays represent coding sequences, they can also be used for genotyping closely related species .
Although Affymetrix expression arrays can be used for genotyping, their drawback is that all except the most recently produced microarrays have been designed with approximately 11 perfect match probes per unigene giving only limited coverage of each gene. Gresham et al. showed that an array designed with 25 bp oligos in a 5 bp overlapping tile format had greater sensitivity (ability to detect true SFPs) in yeast and increased specificity (reduced rate of false positives) in calling SFPs. This overlapping tile design offers technical reproducibility and extensive genome coverage if the number of features on the microarray is sufficient.
Genotyping by microarray hybridization has proven to be challenging in species with complex genomes. Microarrays have been successfully used to detect SFPs in small genomes, for example the 13.5 Mb genome of yeast , the 145 Mb genome of Arabidopsis [1, 8] and, more recently, the 389 Mb genome of rice . Although different algorithms were used for each of these three species, 87% of the SFPs in yeast  and 75% of those in rice and Arabidopsis[1, 3] were independently validated. To identify polymorphisms in the barley genome, complexity was circumvented by using RNA to hybridize to the microarray and 67%  to 80%  validation of SFPs was achieved. When DNA from barley (5,200 Mb genome) was hybridized directly to the Barley1 GeneChip® a significant overlap between SFPs identified using genomic DNA and those identified and validated using RNA was reported . The increased efficiency reported by Rostoks et al. from using RNA is, however, complicated by incomplete and variable transcriptome representation due to tissue- and environment-specific gene expression and false SFP discovery due to alternative splicing or adenylation  associated with sampling RNA versus DNA.
Several types of analyses have been implemented for SFP detection from microarray data. Generally, the data have been processed using expression analysis software to correct for background signal variation using Robust Multi-array Analysis (RMA)  followed by correction for overall signal variation by quantile normalization across arrays . To call SFPs, a modified T-test , Robustified Projection Pursuit (RPP)  and SFP de-viation  have been developed to first estimate the normalized hybridization of a reference set of probes and then test with appropriate statistics or ratios for differential hybridization of specific probes across genotypes. In addition, a maximum likelihood procedure using the source of sequence on the chip as a reference was deve-loped by Gresham et al. to take advantage of overlapping tile data to call SFPs. As each microarray and experiment design tends to be different, new methods for analysis have been developed in attempts to gain greater specificity and sensitivity.
Cultivated lettuce, L. sativa, has substantial genetic and genomic resources including approximately 76,000 ESTs and another 20,000 to 50,000 ESTs in each of four related species (http://compgenomics.ucdavis.edu/). Furthermore, several mapping populations have been developed including a reference mapping population of 214F7 recombinant inbred lines (RILs) between L. sativa cv. Salinas and L. serriola acc. US96UC23 (; RW Michelmore et al., unpublished). This population segregates for multiple agronomic, disease resistance and quality traits. It has approximately 1,500 mapped DNA markers comprised of approximately 700 mapped unigenes with the remainder amplified fragment length polymorphisms (AFLPs), restriction fragment length polymorphisms (RFLPs) and simple sequence repeats (SSRs) . The large number of available sequences and recombinant inbred lines provide ideal resources for further marker discovery and high density mapping. Considerable genetic resources are also available within germplasm collections of L. sativa. Furthermore, several Lactuca species have variable cross-compatibility with L. sativa and represent a diverse genetic resource for investigations of novel alleles and population structure.
In this paper, we describe the development of a microarray designed to provide extensive gene coverage and maximize detection of SFPs for marker discovery and genotyping in lettuce. We analyzed the parents of the reference L. sativa x L. serriola mapping population to demonstrate that DNA representing complex genomes (2,639 Mb)  can be effectively hybridized onto microarrays. Parameters affecting DNA hybridization and accurate detection of polymorphism were optimized. Algorithms from West et al. and Borevitz et al. were modified to take advantage of the overlapping tile design to detect polymorphisms. The use of the multiple probes covering a single position led to the identification of single position polymorphisms (SPPs). Additionally, we assessed SPPs in a diverse panel of Lactuca species concentrating on the cultivated L. sativa.
Identification of a non-redundant consolidated unigene set from Lactuca spp. for design of an oligonucleotide array
The number of ESTs and Unigenes for each species before and after filtering
Unigenes after filtering
Design of microarray with overlapping probes
In collaboration with Affymetrix, probes from both sense and anti-sense strands were selected to create a 2 bp overlapping tiling path (See Additional file 1: Figure S1). The resulting 11.4 million candidate probes designed from the 35,788 submitted sequences were triaged down to ~6.5 million that could be accommodated on the chip through a series of steps based on: 1) Affymetrix probe quality score (> 0.25) except for a select set of unigenes with putative polymorphisms where probes with a quality score above 0.1 were retained. 2) probes matching mitochondrial or chloroplast genomes were discarded. 3) Probes that matched to more than one target were synthesized only once on the chip and computationally associated to corresponding unigene for analysis.
In addition to lettuce probes, background and standard Affymetrix control probes  were included on the microarray. In order to determine background hybridization signal, 13,567 anti-genomic (AG) background probes were synthesized on the microarray, with an average of 1,355 probes representing each G/C bin (probes with the same number of guanines and/or cytosines in the 25 bp probe). These AG probes represent sequences that had no BLAST hits in GenBank at the time of chip design. The use of AG probes reduced the number of probes included on the chip for background correction by 99% compared to the use of mismatch probes (allocated to half of the array positions in pre-vious designs), without substantially compromising the ability to perform accurate background correction . From 8,000 visually interrogated EST contigs, ~2,000 putatively polymorphic regions (50 to 100 bp) were represented from 1,184 contigs.
In total, 6,410,923 lettuce probes representing 35,788 unigenes were synthesized on the microarray. The average and median number of probes representing a unigene were both 187, with ~80% of the unigenes being represented by 50 to 275 probes per sequence (See Additional file 2: Figure S2). Each unigene had an average of 591 bps and median of 603 bps covered by probes; the average number of contiguous stretches of overlapping probes per unigene is 3.1. Due to the selection parameters described above, a contiguous overlapping tile across the unigenes was not possible. Consequently, the average and median lengths of the contiguous stretches of overlapping probes within probe sets are 190 and 120 bps, respectively. Regions of high or low G/C content were sparsely covered by probes (See Additional file 3: Figure S3) likely due to the removal of probes with low Affymetrix probe quality scores. The total number of probes on the array was 6,482,479.
Preparation and hybridization of genomic DNA to the lettuce GeneChip
Improvement of the algorithm for detecting polymorphisms
We also modified the algorithm described by Borevitz et al. for detection of SFPs to exploit the tiling array design. This modified SFP algorithm (MSA) interrogated all weighted probes above background at every 2 bp position, and calculated a D-stat similar to that described by Borevitz et al.. This, however, was done for each position rather than each probe. A false disco-very rate was calculated for each threshold cutoff value for both pairwise analyses using permutation analysis as described by Borevitz et al..
Detection of SPPs between L. sativa cv. Salinas and L. serriola acc. US96UC23
Genomic DNA from L. sativa cv. Salinas (SAL) and wild L. serriola acc. US96UC23 (SER) were hybridized to the GeneChip® in three technical replicates. The Affymetrix .CEL file data were background-corrected by RMA and quantile-normalized across all chips. The data were analyzed for SPPs using both our modified SFPdev and MSA algorithms. SPPs were filtered to require a minimum 4 bp range and two informative probes covering the interrogated position to increase confidence in the SPPs called. SPPs were defined by the range of positions (bp) that met a selected FDR value of 0.1 and a Delta value of 0.2 for the SFPdev and MSA methods, respectively. Furthermore, only the SFPdev values with a ratio (SAL/ SER) less than 1 were considered, as values above one had an actual FDR of 79%. With these requirements the SFPdev method identified 40,462 SPPs in 19,345 contigs; while 40,960 SPPs in 18,290 contigs were identified with the MSA method. The coincidence of reported SPPs between the two methods showed that 73.6% of SPPs detected by the SFPdev method were found by MSA, and 81.1% of SPPs reported by MSA were found by SFPdev.
Observed FDR for each cut-off value for the three SPP prediction methods
Analysis of a Diversity Panel (DP)
Individuals in the diversity panel and their relative species and horticultural class within L. sativa
Number of SPPs and Contigs containing SPPs identified for each class/species by DP analysis
L. Sativa types
Highly parallel genotyping has become an important component of genomics. Hybridization of genomic DNA and RNA to microarrays has been used in the past for detection of polymorphisms between genotypes [1, 4, 5, 10, 22]. However, the previously available arrays for complex genomes only provided limited transcriptome coverage. We developed an array designed to maximize transcriptome coverage while maintaining the possibility of performing other analyses. Our custom designed Lettuce GeneChip® combined the benefits of overlapping probes across unigenes, similar to that demonstrated by Gresham et al. for yeast, with the use of anti-genomic probes to maximize the possible coverage of unigenes while maintaining the sensitivity to detect polymorphisms and retaining appropriate controls to normalize and correct for background noise. The tiling path design allows for multiple assessments of hybridization differences between lines for single positions rather than single assessments of a few positions as obtained with most expression arrays. We developed custom scripts for analysis of our hybridization data taking into account the multiple probes covering a single position as well as filtering out poorly performing probes. We used recent advances in high throughput sequencing technology to validate our SPP calls as well as filter out potentially unreliable data.
Genomic DNA and cDNA are two options for hybridization to an array for SFP detection. The decision of which to use becomes more difficult as genome size and complexity increases. DNA as well as cDNA are both viable targets for species with smaller genomes such as Arabidopsis [1, 5] and rice [3, 23]. However, with larger and more complex genomes such as barley, cDNA was indicated as a more reliable option for hybridization even with the added difficulty of subtracting out expression effects . The genome of lettuce is nearly 17x larger than Arabidopsis although it is half the size of barley. Given the difficulty of accounting for spatio-temporal expression effects as seen in cDNA, we focused on developing methods to use genomic DNA. Rostoks et al. suggested that genomic DNA may be a feasible target in larger genomes with added replication. With the redundancy of the overlapping probes in the lettuce array, the need for additional replication was reduced because they provide technical replicates within a chip. The need for additional replication may also have been reduced by using elevated amounts of genomic DNA and the use of end-labeling rather than BioPrime may have increased the reliability of calls. The protocol used for hybridization of lettuce genomic DNA was also subsequently highly effective for pepper (genome size = 3,000 MB) and other Solanaceae . Furthermore, the use of genomic DNA is a desirable target because SFPs identified using cDNA may be a result of alternate splicing or gene expression differences . Rostoks et al.  indicated that 40% of the SFPs they identified may have been falsely called and partially explained them as being mRNA structural variants. They also reported a high predicted false positive rate of 22% (their mid-stringency cutoff value) for SFPs detected using genomic DNA. We concluded that fragmented, end-labeled genomic DNA provided a suitable target for detection of polymorphisms while reducing false positive sequence polymorphism (2 to 24% in our experiments, Table 2).
The overlapping tile design increases the likelihood of detecting polymorphisms due to redundancy at individual positions, coverage along the contigs and optimal position of the SNP within a probe [1, 11]. Furthermore, the number of probes and hence the possible genome coverage was increased by substituting mismatch probes with AG probes for background correction and normalization of data. Because the peripheral 1 to 6 bases of a 25 bp oligonucleotide are less sensitive than the central bases, in terms of detecting sequence polymorphisms [1, 4], Figure 4, , the tiling strategy reduces the loss of coverage due to probe position. The number and reliability of SPP calls in our experiments demonstrates that the overlapping tiling array design has improved coverage, sensitivity and specificity to detect polymorphisms.
SPP calls were validated using several approaches. The data from the two pair-wise comparisons (MSA and SFPdev) yielded 20 to 41 thousand and 27 to 40 thousand SPPs respectively, depending on the criteria used for specificity and sensitivity. When SPPs from MSA and SFPdev were compared to the 51,552 SNPs detected between RNAseq reads of Salinas and US96UC23, 61.5% and 57.8% were found in or within at least 8 bp of the SPP range respectively, similar to that described by Gresham et al. . However, because of the high FDR associated with duplicated sequences, SPPs that were found to have a duplicated locus within the chip assembly, the gene space assembly or the genome assembly were removed from consideration; one third of the SPPs called that had duplicated loci did not contain a SNP in any of our validation tests. These identified SPPs likely were due to differences between paralogs rather than alleles at a single locus. Due to the increased redundancy (up to 107 genetic replicates) provided by the mapping population of 213 RILs compared to the pair-wise comparison of the parents, SPPs in the SFPdev and MSA pair-wise comparisons that coincided with SPP mapped by Truco et al.  but were absent of a SNP were considered real. Removal of duplicated loci and inclusion of mapped SPPs provided a balance between false positive and false negative rates and allowed us to optimize FDR while still discovering a high number of SPPs. Taking into consideration the lower observed FDR (Table 2) we concluded that the MSA method performed best as a pairwise comparison; however using multiple detection methods would yield a higher confidence in the subset of SPPs identified by both methods.
The SPPs identified in the diversity panel (DP) that were polymorphic between L. sativa cv. Salinas (SAL) and L. serriola acc. US96UC23 (SER) showed a low FDR. However, as a result of the filtering, sensitivity of this analysis (detection of true positives) was reduced compared to the two-genotype analyses by MSA and SFPdev. Specific analysis of the DP data for regions containing known SNPs showed that SFPdev values would have been significant in a pair-wise comparison, between SAL and SER but due to inclusion of data from all genotypes in the DP, the two were not called as polymorphic (data not shown). The lack of some called SPPs in the DP may be due to larger genetic differences between L. perennis, L. virosa, or L. saligna relative to L. serriola and L. sativa. As a result of smaller hybridization differences between the more closely related genotypes, genotypes differing at a locus may have been grouped together reducing the number of SPPs called between the two genotypes. Consequently, the DP analysis showed a lower false positive rate, but a higher false negative rate when comparing SAL and SER to sequence and mapping data.
As part of our goal was to investigate the diversity and relationships of the genotypes in the DP, SPPs identified by the DP analysis were evaluated. Removal of SPPs in duplicated regions with inconsistent data or missing data (see Methods) was a reasonable method of removing unreliable data as these data may be from poorly performing probes in one or all replicates, heterozygous loci, paralogous genes or deleted genes. There was not a large difference in the observed FDRs for the three SPFdev cutoff values (1.2, 1.5, & 2.0) for the DP analysis; so in order to maximize the number of markers used in our phylogenetic analysis and principal component analysis, we used the least stringent cutoff value of 1.2. As the assumptions for analysis with the PHYLIP  package were not violated with the large number of markers, they were left as independent. To meet the constraints of the PC analysis software, markers were limited to those that were mapped.
The markers discovered in our DP analysis were used to generate a phylogenetic tree showing species separation with 100% boot strap support. L. virosa and L. saligna are sexually incompatible species with L. sativa and appear to be more closely related to each other than to other species in the DP. Our data supports the conclusion by Kesseli et al., that these two species are not progenitors of L. sativa. By limiting markers to those polymorphic within cultivated lettuce we are able to separate most genotypes into classes representing each of the plant types. The butterhead type formed a distinct clade from the iceberg and cos types with 100% bootstrap support. However, the leafy type and the Batavia type both showed a wide distribution across the L. sativa clade. This is not unexpected and may reflect admixture between types during breeding programs. Alternatively, this distribution may indicate that these types are artificial polyphyletic groups based on loose morphological criteria . The leafy types are non-heading with a broad range of leaf morphology [28, 29]. Batavia types vary from heading to non-heading phenotypes. Batavia and iceberg cultivars are both considered crisphead types ; however our phylogenetic and PC analyses showed that the two did not cluster together and are significantly different from each other (See Additional file 8: Table S1).
Rapid advancements in sequencing technology today are changing the methods for genetic analyses. Microarray technology presented in this paper yielded an in depth analysis of diversity for lettuce germplasm separating even closely related lines such as the crisphead class. It also has potentially several other uses including: detection of copy number variants, splice site identification, expression analysis or use with other species within the Compositae. The SPPs identified in this study were highly reproducible and showed similar false positive results to current sequencing methods in the literature [30, 31]. This technology has also been used to create an ultra-dense, inter-specific genetic map between L. sativa cv. Salinas and L. serriola acc. US96UC23 to dissect phenotypic traits as well as validate and align genomic assemblies of lettuce into chromosomal linkage groups .
We designed and exploited a custom lettuce microarray using an overlapping tiling path and by using anti-genomic probes rather than mismatch probes to provide comprehensive unigene coverage. Our protocols for genomic DNA preparation and labeling, assisted by positional vs. feature-based analyses reliably identified DNA polymorphisms using both pair-wise genotype comparisons as well as a highly parallel comparison within a diverse panel of genotypes including five species and focused on the cultivated L. sativa. The phylogenetic and principal component analyses clearly distinguished species while the analysis of L. sativa supports previous analyses of cultivated lettuce and revealed differences among the more heterogeneous horticultural types as well as polymorphisms within the most genetically narrow type.
EST assemblies of Lactuca spp.
ESTs for five Lactuca species (L. sativa, L. serriola, L. saligna, L. virosa and L. perennis) generated as part of the Compositae Genome Project (http://compgenomics.ucdavis.edu) were trimmed to phred scores of 15 and vector sequences were removed with custom Perl scripts (https://code.google.com/p/atgc-xyz/source/browse/#svn%2Ftrunk). Trimmed ESTs for each of the five species of Lactuca were assembled separately using CAP3 . To minimize assembly of paralogous sequences, stringent assembly parameters were used for CAP3 (95% identity, 100 nt minimum overlap). The assemblies were combined iteratively using L. sativa as a base and adding unique sequences from individual species in order of genetic distance from L. sativa.
Genomic DNA was extracted following the protocol described by Kozik  with minor modifications. Frozen leaf tissue from three-week old, greenhouse-grown plants (2.5 g fresh weight) was ground in liquid nitrogen. Fifteen ml of 2X extraction buffer (100 mM Tris–HCl pH 8.0, 1.4 M NaCl, 20 mM EDTA, 2% w/v CTAB, 10 μl/ml β-mercaptoethanol) was added to the tissue, mixed and incubated at 65°C. One volume of chloroform:isoamyl alcohol (24:1) (ChIA) was added and mixed, and the sample centrifuged at 3,500 rpm for 20 min twice. The aqueous phase was then transferred and 3 to 3.5 volumes of precipitation buffer (50 mM Tris–HCl pH 8.0, 10 mM EDTA, 1% w/v CTAB) were added. The sample was incubated overnight at room temperature to precipitate the DNA and then centrifuged at 3500 rpm for 15 min. The DNA pellet was washed with dH2O and centrifuged for 10 min; 5 ml of 1.5 M NaCl and 6 μl of 10 mg/ml RNaseA was added to the pellet and incubated at 37°C until completely re-suspended. A chloroform extraction was performed as above to remove RNaseA and any additional contaminants. The aqueous phase was collected and DNA was precipitated and washed with ethanol. Samples were centrifuged again for 10 min at 3,500 rpm. The pellet was allowed to dry then re-suspended in 100 μl ddH2O. The sample was diluted 1:100 and run on a 1% agarose gel with known amounts of uncut lambda DNA to estimate the concentration of the genomic DNA. The 2.5 g of starting material yielded approximately 200 μg of genomic DNA, which was sufficient for hybridization to six chips at 30 μg per chip.
Whole genome amplification
The GenomePlex® Complete WGA kit (WGA, Sigma-Aldrich Corp, St. Louis, MO, USA) was evaluated for whole genome amplification. The maximum amount of genomic DNA that could be accommodated by the WGA kit (20 ng) was amplified and then re-amplified using the GenomePlex WGA re-amplification kit following the manufacturer’s instructions. Ten parallel re-amplification reactions were required to obtain sufficient amounts of amplified DNA for hybridization. The samples were then purified using Qiaquick PCR purification columns (QIAGEN, Valencia, CA, USA), pooled and concentrated using a Millipore Centricon YM-10 (Millipore, Temecula, CA, USA). Concentrated samples were then re-suspended in 37 μl ddH2O and quantified using the ND-1000 (Nanodrop, Wilmington, DE, USA).
DNA preparation, labeling and hybridization
In a 50 μl reaction, fragmentation of DNA was achieved by incubation of 30 μg genomic DNA with 0.0015 U RQ1 DNase I (Promega, Madison, WI, USA) per μg of DNA in the presence of: 0.175 mg/ml BSA, 10 mM Tris-acetate pH 7.5, 10 mM magnesium acetate, and 50 mM potassium acetate. Reactions were carried out for 30 minutes at 37°C followed 15 min at 99°C. Fragmentation profiles were examined by electrophoresis of 1.5 μl of fragmented product on a 2% agarose gel alongside a 50 bp DNA ladder (Fermentas, Glen Burnie, MD, USA) (Figure 3). Fragmented DNAs ranging between 25 to 250 bp were considered acceptable for hybridization.
The 3’ termini were end-labeled by addition of a biotinylated oligonucleotide (Affymetrix 7.5 mM Labeling Reagent) using terminal deoxynucleotidyl transferase (TdT) (Promega, Madison, WI, USA) as per a modified Affymetrix Labeling of Fragmented Double-Stranded DNA instructions . The reaction was scaled up to accommodate 30 μg of fragmented DNA rather than the 7.5 μg fragmented cDNA. The BioPrime DNA Labeling System (Invitrogen, Carlsbad, CA, USA) was also tested for labeling DNA as per the manufacturer’s instructions.
Hybridizations were carried out following the protocol described in Affymetrix GeneChip® Whole Transcript Double-Stranded Target Assay (WTDSTA) Manual (pp. 51–52; Affymetrix 2005) with minor modification. Samples were comprised of 70 μl of fragmented and labeled DNA containing the elevated amount of DNA (30 μg total), the remaining components of the GeneChip® Hybridization Kit (Affymetrix, Santa Clara, CA, USA), adjusted to 1.1x standard volumes, and the 20x hybridization control from the Hybridization Control Kit (Affymetrix, Santa Clara, CA, USA) were added as per instructions for a total volume of 220 μl. After centrifugation, the microarray cassette was injected with 200 μl of the cocktail. Post hybridization, washing and staining were carried out as described in the Affymetrix GeneChip® WTDSTA Manual using program FS450_0001 (pp. 55–61 Affymetrix 2005). Hybridization intensities were measured using an Affymetrix GeneChip® Scanner 3000 with 7 G up-grade.
Data processing and management
The raw .CEL files generated by the GCOS software (Affymetrix, Santa Clara, CA, USA) were subjected to background correction by using the Robust Multiarray Analysis (RMA) method and quantile normalization  of the Affymetrix package  with R and Bioconductor Affymetrix package . A MySQL 5.0 database was developed to curate all the information regarding experiments, probe characteristics, hybridization data, and quality control data (https://pgfmars.ucdavis.edu/phpmyadmin). A web interface was designed using PHP scripts for retrieval of the data. All scripts of our “SPPscan Package for Lettuce GeneChip® Data Management, Processing and Analysis” and detailed descriptions can be obtained from http://chiplett.ucdavis.edu/.
SPPdev detection algorithm
Where R(Y i ) is the average value of all the high quality probes on the same hybridized lettuce chip with the same GC content as probe Y i . The high quality probes were defined as those probes whose signal intensities exceeded those of a defined percentile of anti-genomic probes. The threshold value selected for our analysis was 90% (See Results).
Where w k is the weighting factor for the high quality probe k that covers position s in the unigene. The weighting factor takes into account the distance from the central nucleotide of each probe to the evaluated position (Figure 5). The weighting factor was determined empirically by analyzing the hybridization behavior of probes covering a training set of 1,000 SNPs obtained from the alignment of multiple ESTs from L. sativa cv. Salinas with L. serriola acc. US96UC23 (See Results). All scripts are available from http://chiplett.ucdavis.edu/.
Statistical evaluation of putative SPPs
Arrange rows for all the probes of the unigene and columns for all the replicates of the two genotypes.
For the b th permutation, b = 1, 2,…B:
Permute the columns of the data matrix obtained in S1
Compute for each position s, s = 1, 2,…m.
After carrying out the B permutations, sort all the values obtained from Step 2 in increasing order. The permutation p-value for testing the null hypothesis (that no polymorphism exists at position s) is determined by:
Obtain a FDR-adjusted p-value , P * s for each ps.
A position s is claimed to be a putative SPP if P * s < 0.01.
Note that the null distribution of R s for each position s is assumed to be the same. Therefore, the null distribution of R s is estimated by pooling the permutation distributions of R s at all positions.
Alternative SPP detection algorithm (Modified SFP Algorithm; MSA)
Mean of parent A or B at a given 2 bp position is natural logarithm (L n ) of weighted average intensity of all probes that interrogating a position,
S 0 is a constant for all positions, which is the median of all S i values for the entire 2 bp positions on the chip (>10,000,000 positions).
RNA was isolated from multiple tissues and treatment conditions (http://cgpdb.ucdavis.edu/cgpdb2/data_info_files/CGP_Library_Construction_V02.html) and pooled in equimolar concentration for each of the two genotypes. cDNA was synthesized and prepared for sequencing following protocols as described . L. sativa cv. Salinas and L. serriola acc. US96UC23 were sheared and selected for fragments around 260 bp and 200 bp respectively. Fragments were normalized using duplex specific nuclease (DSN, Evrogen, Moscow, Russia) after denaturation and re-association for five hours to reduce high copy sequences digesting double stranded DNA. The L. sativa cv. Salinas library was sequenced in one direction while the L. serriola acc. US96UC23 library was sequenced as paired-ends. Sequences are available in the Sequence Read Archive at NCBI (Study numbers SRP004854 and SRP008310).
Using BWA  and SAMtools  to identify a set of high quality SNPs, heterozygous positions, and InDels in the IGA set. Custom scripts were used to identify false positive SPPs identified by each of the SFPdev and MSA methods by comparison to the IGA set, the SNPs mined from the EST assembly used to develop the array, and SPP markers mapped in the core mapping population . A SPP was considered a true positive if SNPs were within 8 bp of the SPP ranges. An 8 bp range on either side of detected SPPs was chosen to account for detection of SPPs with the overlapping 25 bp oligo design and empirically determined sensitivity of oligos (Figure 5). Additionally, if the SPP range in the pair-wise comparison overlapped the SPP range mapped within the core RIL population, it was considered real. SPPs with no IGA sequence were excluded as we were not able to identify if they were true or false SPPs.
In order to identify if the apparent SPP was present in multiple contigs, each SPP plus and minus 8 bases on either side of the SPP range was compared using BLAST to three separate sets of sequences: a lettuce whole genome assembly of shotgun Illumina reads (Michelmore et al., unpublished data), a gene space assembly (Kozik, unpublished data) and the EST assembly used for the array design (http://compgenomics.ucdavis.edu) using BLASTn . To maintain confidence in the BLAST hits there had to be a minimum of 95% of the subject in the query with no more than two mismatches. SPPs identified as duplicated loci were removed from the analysis. SPPs detected in the SFPdev pair-wise comparison method that had a SFPdev ratio greater than one were removed as we empirically determined a 79% false positive rate through our validation procedures.
Diversity panel analysis
As described in Truco et al. the RIL algorithm calculates the SPPdev value for each array and plots that value for each position. The distribution of data points for each position is then evaluated for a bimodal distribution while maximizing the number of individual arrays under each peak. Those that fall between the identified populations of the bimodal distribution, or are completely missing, were treated as a missing data point (−). Arrays with hybridization signals in the higher part of the distribution were designated the A allele while those with lower hybridization values (because they deviated from the probe sequence), were designated the B allele. Data points for each array were then summarized based for each individual or cultivar as follows. Marker calls based on three replicates were designated as A, B, C, D, I (inconsistent), or – (missing). An A marker call resulted from an A in all three replicates. A B marker call resulted from a B in all three replicates. C calls were not A and were missing in one of the three reps i.e. (B/B/-) and D was not B and missing one of the three reps i.e. (A/A/-). An I was any combination of the three chips that contained two missing and an A or B or, any combination that included an A and B call from one of the three reps. – was missing data from all three reps.
Utilizing the package PHYLIP-3.69  the SPP markers were treated as restriction site markers and scored as 1 or 0. The SeqBoot module was used to create 100 re-samples data sets for bootstrap calculation. The resulting 100 replicates were used to create distance matrices subsequently used for construction of phylogenetic trees with the Fitch module. Global rearrangement and randomized input order of species with 10 jumbles for each of the multiple data sets was used to construct 100 trees for use in building a consensus tree with bootstrap values. The Consense module was then used with species L. perennis used as an out-group root. The resulting consensus tree was then visualized with Mega4  and branches were rotated for legibility. See Additional file 9: Table S2.
Principal component analysis
For principal component analysis, markers within contigs (9,513) that were common to the L. sativa cv. Salinas by L. serriola acc. US96UC23 map  were converted to haplotype frequencies using custom Perl scripts. Frequencies were calculated as the number of times the haplotype occurs in the panel for a given contig divided by the total number of genotypes in the panel. A limitation on the number of markers to use with the SAS/STAT software required us to reduce our set of markers to 10,000 markers or less. Creating haplotypes for all markers within a contig allowed us to reduce the number of input data points to 8,381 while retaining more genetic information than a randomly selected marker for each mapped contig. The markers used for calculation of haplotype frequencies in the entire DP were then filtered to include those polymorphic within L. sativa lines only and haplotype frequencies were calculated for 3,311 contigs to resolve L. sativa varieties. Frequencies of haplotypes in the DP panel were used as input for principal component analysis using the SAS/STAT® software’s PRINCOMP procedure (Version 9.1.3, SAS Institute, Cary, NC). Eigenvalues for each of the first three principal components were calculated for each genotype and an analysis of variance (ANOVA) and means separation using Students t-test were performed using JMP (JMP, Version 7. SAS Institute Inc., Cary, NC, 1989–2007) to determine if species and/or classes were significantly different. The two data sets described above were used to differentiate lettuce species and lines. See Additional file 10: Table S3 and Additional file 11: Table S4.
We thank Smitha Mathrakott, Gene Wong, and Charlyn Suarez for technical assistance, Katherine Stapleton, Eric Schell, Gene Tanimoto (Affymetrix) for technical advice in array design, and Belinda Martineau for editorial assistance. This project was supported by an award from the University of California UC Discovery Program Bio04-1047 with matching funds from Rijk Zwaan BV, Enza Zaden BV, and Vilmorin Co. In addition, no personal funds were taken from these entities from any authors. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
- Borevitz JO, Liang D, Plouffe D, Chang HS, Zhu T, Weigel D, Berry CC, Winzeler E, Chory J: Large-scale identification of single-feature polymorphisms in complex genomes. Genome Res. 2003, 13: 513-523. 10.1101/gr.541303.PubMed CentralView ArticlePubMedGoogle Scholar
- Das S, Bhat PR, Sudhakar C, Ehlers JD, Wanamaker S, Roberts PA, Cui X, Close TJ: Detection and validation of single feature polymorphisms in cowpea (vigna unguiculata l. Walp) using a soybean genome array. BMC Genomics. 2008, 9: 107-10.1186/1471-2164-9-107.PubMed CentralView ArticlePubMedGoogle Scholar
- Kumar R, Qiu J, Joshi T, Valliyodan B, Xu D, Nguyen HT: Single feature polymorphism discovery in rice. PLoS One. 2007, 2: e284-10.1371/journal.pone.0000284.PubMed CentralView ArticlePubMedGoogle Scholar
- Rostoks N, Borevitz JO, Hedley PE, Russell J, Mudie S, Morris J, Cardle L, Marshall DF, Waugh R: Single-feature polymorphism discovery in the barley transcriptome. Genome Biol. 2005, 6: R54-10.1186/gb-2005-6-6-r54.PubMed CentralView ArticlePubMedGoogle Scholar
- West MA, van Leeuwen H, Kozik A, Kliebenstein DJ, Doerge RW, Clair DA, Michelmore RW: High-density haplotyping with microarray-based expression and single feature polymorphism markers in arabidopsis. Genome Res. 2006, 16: 787-795. 10.1101/gr.5011206.PubMed CentralView ArticlePubMedGoogle Scholar
- Werner JD, Borevitz JO, Warthmann N, Trainer GT, Ecker JR, Chory J, Weigel D: Quantitative trait locus mapping and DNA array hybridization identify an flm deletion as a cause for natural flowering-time variation. Proc Natl Acad Sci U S A. 2005, 102: 2460-2465. 10.1073/pnas.0409474102.PubMed CentralView ArticlePubMedGoogle Scholar
- Kim S, Zhao K, Jiang R, Molitor J, Borevitz JO, Nordborg M, Marjoram P: Association mapping with single-feature polymorphisms. Genetics. 2006, 173: 1125-1133. 10.1534/genetics.105.052720.PubMed CentralView ArticlePubMedGoogle Scholar
- Borevitz JO, Hazen SP, Michael TP, Morris GP, Baxter IR, Hu TT, Chen H, Werner JD, Nordborg M, Salt DE: Genome-wide patterns of single-feature polymorphism in arabidopsis thaliana. Proc Natl Acad Sci U S A. 2007, 104: 12057-12062. 10.1073/pnas.0705323104.PubMed CentralView ArticlePubMedGoogle Scholar
- Jiang R, Marjoram P, Borevitz JO, Tavare S: Inferring population parameters from single-feature polymorphism data. Genetics. 2006, 173: 2257-2267. 10.1534/genetics.105.047472.PubMed CentralView ArticlePubMedGoogle Scholar
- Gresham D, Ruderfer DM, Pratt SC, Schacherer J, Dunham MJ, Botstein D, Kruglyak L: Genome-wide detection of polymorphisms at nucleotide resolution with a single DNA microarray. Science. 2006, 311: 1932-1936. 10.1126/science.1123726.View ArticlePubMedGoogle Scholar
- Cui X, Xu J, Asghar R, Condamine P, Svensson JT, Wanamaker S, Stein N, Roose M, Close TJ: Detecting single-feature polymorphisms using oligonucleotide arrays and robustified projection pursuit. Bioinformatics. 2005, 21: 3852-3858. 10.1093/bioinformatics/bti640.View ArticlePubMedGoogle Scholar
- Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP: Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003, 4: 249-264. 10.1093/biostatistics/4.2.249.View ArticlePubMedGoogle Scholar
- Bolstad BM, Irizarry RA, Astrand M, Speed TP: A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003, 19: 185-193. 10.1093/bioinformatics/19.2.185.View ArticlePubMedGoogle Scholar
- Truco MJ, Antonise R, Lavelle D, Ochoa O, Kozik A, Witsenboer H, Fort SB, Jeuken MJW, Kesseli RV, Lindhout P: A high-density, integrated genetic linkage map of lettuce (lactuca spp.). Theor Appl Genet. 2007, 115: 735-746. 10.1007/s00122-007-0599-9.View ArticlePubMedGoogle Scholar
- Zohary D: The wild genetic-resources of cultivated lettuce (lactuca-sativa l). Euphytica. 1991, 53: 31-35. 10.1007/BF00032029.View ArticleGoogle Scholar
- Arumuganathan K, Earle ED: Nuclear DNA content of some important plant species. Plant Mol Biol Rep. 1991, 9: 11-View ArticleGoogle Scholar
- Huang XQ, Madan A: Cap3: a DNA sequence assembly program. Genome Res. 1999, 9: 868-877. 10.1101/gr.9.9.868.PubMed CentralView ArticlePubMedGoogle Scholar
- Affymetrix: Quality assessment of exon and gene arrays. 2007Google Scholar
- Affymetrix: Exon array background correction. 2005Google Scholar
- Truco MJ, Ashrafi H, Kozik A, van Leeuwen H, Bowers J, Reyes Chin Wo S, Stoffel K, Lavelle D, Hill T, Caldwell D: An ultra-high-density, transcript-based, genetic map of lettuce. Prep. 2012Google Scholar
- Felstein J: Phylip - phylogeny inference package (version 3.2). Cladistics. 1989, 5: 164-166.Google Scholar
- Bernardo AN, Bradbury PJ, Ma HX, Hu SW, Bowden RL, Buckler ES, Bai GH: Discovery and mapping of single feature polymorphisms in wheat using affymetrix arrays. BMC Genomics. 2009, 10: 19-10.1186/1471-2164-10-19.View ArticleGoogle Scholar
- Kim S-H, Bhat PR, Cui X, Walia H, Xu J, Wanamaker S, Ismail AM, Wilson C, Close TJ: Detection and validation of single feature polymorphisms using rna expression data from a rice genome array. Bmc Plant Biology. 2009, 9: 65-10.1186/1471-2229-9-65.PubMed CentralView ArticlePubMedGoogle Scholar
- Hill T, Hill TA, Ashrafi H, Yao J, De Jong W, Francis D, Kozik A, Van Deynze A: The application of a whole genome pepper array to solanaceae crops. Solanaceae Genomics. 2008, Cologne, GermanyGoogle Scholar
- Luo ZW, Potokina E, Druka A, Wise R, Waugh R, Kearsey MJ: SFP genotyping from affymetrix arrays is robust but largely detects cis-acting expression regulators. Genetics. 2007, 176: 789-800.PubMed CentralView ArticlePubMedGoogle Scholar
- Devries IM: Crossing experiments of lettuce cultivars and species (lactuca sect lactuca, compositae). Plant Systematics and Evolution. 1990, 171: 233-248. 10.1007/BF00940608.View ArticleGoogle Scholar
- Kesseli R, Ochoa O, Michelmore R: Variation at RFLP loci in lactuca spp and origin of cultivated lettuce (l-sativa). Genome. 1991, 34: 430-436. 10.1139/g91-065.View ArticleGoogle Scholar
- Kristkova E, Dolezalova I, Lebeda A, Vinter V, Novotna A: Description of morphological characters of lettuce (lactuca sativa l.) genetic resources. Hortic Sci. 2008, 35: 113-129.Google Scholar
- deVries IM: Origin and domestication of lactuca sativa l. Genet Resour Crop Evol. 1997, 44: 165-174. 10.1023/A:1008611200727.View ArticleGoogle Scholar
- Hamilton JP, Hansey CN, Whitty BR, Stoffel K, Massa AN, Van Deynze A, De Jong WS, Douches DS, Buell CR: Single nucleotide polymorphism discovery in elite north american potato germplasm. Bmc Genomics. 2011, 12: 302-10.1186/1471-2164-12-302.PubMed CentralView ArticlePubMedGoogle Scholar
- Iorizzo M, Senalik DA, Grzebelus D, Bowman M, Cavagnaro PF, Matvienko M, Ashrafi H, Van Deynze A, Simon PW: De novo assembly and characterization of the carrot transcriptome reveals novel genes, new markers, and genetic diversity. Bmc Genomics. 2011, 12: 389-10.1186/1471-2164-12-389.PubMed CentralView ArticlePubMedGoogle Scholar
- Kozik A: Fine mapping of the sym2 locus of pea linkage group 1. 1996, Kozik: Proefschrift WageningenGoogle Scholar
- Affymetrix: Genechip® whole transcript (wt) double-stranded target assay manual. 2005Google Scholar
- Irizarry RA, Bolstad BM, Collin F, Cope LM, Hobbs B, Speed TP: Summaries of affymetrix genechip probe level data. Nucleic Acids Res. 2003, 31: e15-10.1093/nar/gng015.PubMed CentralView ArticlePubMedGoogle Scholar
- Gautier L, Cope L, Bolstad BM, Irizarry RA: Affy–analysis of affymetrix genechip data at the probe level. Bioinformatics. 2004, 20: 307-315. 10.1093/bioinformatics/btg405.View ArticlePubMedGoogle Scholar
- Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J: Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004, 5: R80-10.1186/gb-2004-5-10-r80.PubMed CentralView ArticlePubMedGoogle Scholar
- Good P, Good P: Springer series in statistics: Permutation tests: A practical guide to resampling methods for testing hypotheses. 1994, Springer-Verlag: Springer-Verlag New York IncGoogle Scholar
- Benjamini Y, Hochberg Y: Controlling the false discovery rate - a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B-Methodological. 1995, 57: 289-300.Google Scholar
- Illumina: Mrna sequencing sample preparation guide. 2009Google Scholar
- Li H, Durbin R: Fast and accurate long-read alignment with burrows-wheeler transform. Bioinformatics. 2010, 26: 589-595. 10.1093/bioinformatics/btp698.PubMed CentralView ArticlePubMedGoogle Scholar
- Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, Genome Project Data P: The sequence alignment/map format and samtools. Bioinformatics. 2009, 25: 2078-2079. 10.1093/bioinformatics/btp352.PubMed CentralView ArticlePubMedGoogle Scholar
- Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: Basic local alignment search tool. J Mol Biol. 1990, 215: 403-410.View ArticlePubMedGoogle Scholar
- Tamura K, Dudley J, Nei M, Kumar S, Mega4: Molecular evolutionary genetics analysis (mega) software version 4.0. Mol Biol Evol. 2007, 24: 1596-1599. 10.1093/molbev/msm092.View ArticlePubMedGoogle Scholar