Sequence diversity in three tomato species: SNPs, markers, and molecular evolution

EST assembly

The unigene database version 200607 build 1 from SGN consists of 239,172 ESTs grouped into 34,829 unigenes, each one containing between 1 and 1087 ESTs, with an average of 6.9 ESTs per unigene (data not shown). Unigenes can be formed by EST sequences from any tomato race or species, although most unigenes (87%) contain sequences from a single species (Figure 1).

We mined sequence assemblies of S. lycopersicum (L) ESTs for intraspecific polymorphisms and assemblies containing sequences from both S. lycopersicum and S. habrochaites (LxH) or S. lycopersicum and S. pennellii (LxP) for interspecific SNPs. To reduce the complexity of each analysis we reassembled the unigenes using only ESTs from the species relevant to the analysis. The consensus sequence for each reassembled contig was then aligned to the original unigene sequence and, when available (88.89% of the unigenes, data not shown), to the predicted coding sequence (CDS). In the cases where a predicted CDS was available, each nucleotide in the assembly was annotated as 5' UTR, CDS or 3' UTR. We assigned quality scores to each position in the assemblies by calculating separately for each species the sum of the qualities of all nucleotides at that particular position.

Several filters were applied to diminish the number of false positive SNPs predicted in our analyses. Assembly regions that did not align with unigene sequences were removed. We also discarded positions where only a single sequence was found and, in the interspecific analyses, assembly regions where sequences from only one of the species assayed were detected. Using the remaining portions of the assemblies we determined an optimum quality threshold for each analysis at which the average sequence quality was maximized (see Methods, Figure 2). Assembly positions that did not pass the quality threshold were removed, leaving for SNP mining a total of 4,712 unigenes spanning 1,736 Kb of interspecific assemblies and 19,159 unigenes including 11,058 Kb of intraspecific S. lycopersicum data (Table 1).

Analysis	LxP			LxH			L
	CDS	UTR	NA	CDS	UTR	NA	CDS	UTR	NA
Unigenes analyzed		2249			2463			19159
Kb analyzed	504.9	116.8	21.5	860.1	203.6	29.7	8446.3	1962.1	649.6
	3728/	1104/		8457/	3501/	349/	14455/	8229/	919/
SNPs/indels	206	534	92/30	614	1611	64	1980	3006	221
Unigenes with SNPs		1341			1934			5084
SNP rate		0.8852			1.3349			0.2605
SNP rate per region	0.779	1.4025	0.5682	1.0546	2.5107	1.3925	0.1946	0.5726	0.1755

Only nonconsecutive SNPs are displayed. NA- gene region is not available for unigenes in which a CDS could not be predicted. SNP rates are calculated as the number of SNPs per 100 bp.

In every analysis the majority of nucleotides assembled were located in predicted coding regions. Between 2.8 and 6.2% of the positions could not be assigned to a gene region due to the lack of predicted CDS (Table 1). The average quality of the nucleotides considered in the assemblies was always highest in the 5' regions and lowest in the 3' regions (Figure 3). This could be explained by the bias towards sequencing the ESTs from the 5', as an average of 95.87% of the sequences from each library were readings from the top strand (data not shown). Interestingly, assemblies from unigenes for which there is no predicted CDS had overall average qualities similar to the 5' region (Figure 3). This raises the possibility that those unigenes contain mostly 5' UTR sequence, in turn explaining why ESTscan [29, 30] was not able to predict a CDS.

SNP mining

We mined the resulting sequence assemblies for nonconsecutive SNPs taking into account the position of the SNP with respect to the predicted gene regions. In total 40,834 substitutions and 8,266 insertions/deletions (indels) were detected (summarized in Table 1; see Additional files 1, 2, 3 &4 for full data). Several observations confirm the performance of our SNP detection algorithms. As shown in Table 1, indels appear between 2.6 and 8.7 times more frequently in UTRs than in the CDS. Insertion and deletion of nucleotides are likely to produce changes in the open reading frames disrupting correct translation, therefore selection against them is expected in coding regions. Similarly, predicted noncoding regions yielded between 1.8 and 2.9 times more SNPs per kilobase analyzed than coding regions (Table 1), as reported in previous tomato analyses [15, 19, 23, 25]. We next used the estimated CDS for each unigene to calculate the codon position for each SNP. In translated regions the majority of the mutations are located in the third position of the codon (Figure 4). This result was expected since selective pressure in coding regions reduces the number of non-synonymous substitutions and the redundancy of the genetic code is mainly in the third base of the codons. We then calculated hypothetical codons in the UTRs by extending the ORF from the origin of transcription in the predicted CDS. The third position bias disappeared in these hypothetical noncoding regions (Figure 4) suggesting that on average the SNP, CDS, and codon position predictions are accurate. In addition, when looking at all SNPs in coding regions together, the percentage of non-synonymous SNPs is 46.37%, similar to genome-wide analyses in Arabidopsis (45.34%) [31] and humans (46.46%) [32].

Regarding the differences between intra and interspecific analyses, we found only 26.5% of the analyzed loci to be polymorphic when mining a single species versus the 69.5% when analyzing pairs (Table 1). For interspecific analysis the SNP rates were comparable to those published before between S. pennellii and S. lycopersicum (1.02 to 1.61 SNPs/100 bp) making no distinction between coding and non coding regions [14, 33]. The same holds true for S. lycopersicum intraspecific analyses, where reported SNP rates range between 0.0117 and 0.585 SNPs/100 bp on all EST sequences analyzed regardless of their location in the gene [14–16, 25].

SNP representation in Gene Ontology classes

We reasoned that certain gene families might be more variable among tomato cultivars and species than others. For example, since tomato species are adapted to diverse environments, genes involved in environmental response might accumulate a higher number of non-synonymous SNPs due to selection. We used Gene Ontology (GO) categories [34], which group genes into functionally related classes, to assess the differential polymorphic rates of unigenes encoding specific classes of proteins. We assigned GO categories to each nucleotide position in the assemblies based on the closest Arabidopsis thaliana homolog, and looked for categories with over- or under-representation of non-synonymous SNPs. To achieve higher power we grouped together the data from the interspecific analyses. The larger number of sequences available in the intraspecific analysis allowed greater statistical power, although the most over- and underrepresented gene classes are in agreement in both analyses (Figure 5). As expected, GO categories related to environment interaction, such as responses to stress and abiotic stimulus, had an over-representation of non-synonymous polymorphisms both between and within species (p < 0.001, Figure 5). On the other hand, categories involved in basic biological processes and transcription regulation showed a relative lack of polymorphisms, as had been found in Arabidopsis thaliana accessions [31]. Surprisingly, genes encoding ribosomal proteins also presented more SNPs that expected, perhaps due to the challenge of distinguishing homologs from paralogs in this gene family, which shows complex patterns of copy number variants in Arabidopsis [35].

McDonald-Kreitman test

The neutral evolution theory holds that most within species polymorphisms and between species differences have little fitness consequence [36]. As a result, the ratio of non-synonymous to synonymous substitutions should be similar for a particular gene both within and between species. Deviations from this can be a sign of non-neutral evolution [26]. To examine this prediction, we constructed alignments of the wild and cultivated consensus sequences containing only the good quality positions and the SNPs predicted by our algorithms, and surveyed each alignment for signatures of positive selection using the McDonald-Kreitman test [26]. With this method we tested 1425 unigenes from the S. lycopersicum and S. habrochaites analysis and 924 unigenes from the S. lycopersicum and S. pennellii analysis. Before correction for multiple testing we found significant excess of non-synonymous mutations (p < 0.05) in 16 unigenes in the LxH analysis and 3 unigenes in the LxP analysis (data not shown). However, none of those unigenes survived the Benjamini and Hochberg correction for multiple testing [37].

Maximum likelihood codon-substitution models

Molecular marker design

Using the information gathered from the assemblies we developed a set of molecular markers for detecting the predicted SNPs. Successful design of PCR based molecular markers from EST sequences requires the ability to avoid amplifying introns. Since many Arabidopsis intron positions are conserved in tomato [23], we used this information (see Methods) to inform the design of molecular markers for the polymorphisms identified in the interspecific analyses. First, we tested our intron predictions by designing three primer pairs surrounding predicted introns. One of those produced a band 700 bp greater than expected without introns and the other two did not amplify, suggesting the existence of introns where predicted (data not shown). We did not carry out similar analysis for the intraspecific SNPs due to the lack of information within SGN regarding the S. lycopersicum accessions used. Although it is well known that S. lycopersicum ESTs in SGN come from several cultivars, the information for each individual EST library is not available. For each high-quality, interspecific SNP we developed a database containing suggested primers to amplify the SNP region, restriction enzymes to detect the polymorphism, and predicted fragment sizes before and after digestion of each allele. This information can be accessed at http://www.plb.ucdavis.edu/labs/maloof/TomatoSNP/index.html.

We validated a subset of markers by using 491 primer pairs designed to amplify fragments from S. lycopersicum and from S. pennellii, S. habrochaites, or both. We calculated the size differences between the amplified products cut with the appropriate restriction enzymes when needed or uncut in the case of indels. From the 491 primer pairs, 6 pairs amplified fragments that were smaller than expected and were discarded. Another 48 pairs failed to amplify fragments in at least one species, leaving us with 437 primer pairs that would test 281 polymorphisms between S. lycopersicum and S. habrochaites and 228 polymorphisms between S. lycopersicum and S. pennellii. Among these, 30 primers pairs yielded bands that were bigger than expected, probably due to the existence of non-conserved introns, nevertheless 23 of these were polymorphic as expected after restriction enzyme digestion. For the remaining amplifications we found 87% (220 of 261 LxH and 193 of 210 LxP) of the molecular markers to work as predicted. It is worth noting that 10% of the successful markers showed bands corresponding to both alleles in at least one of the species tested, suggesting heterozygosity in the lines used or the amplification of fragments from a family of genes sharing that sequence. Markers that failed could be due to undetected SNPs in the sequence that modify the predicted restriction sites, or errors in SNP prediction. All tested markers are available in the Additional file 5.

Sequence data

Fasta files containing version 200607 build 1 of the unigene sequences, EST sequences, EST qualities and estimated coding region for each unigene were downloaded from SGN. ESTs in this database originated from the sequencing of at least 43 S. lycopersicum cDNA libraries belonging to at least two different accessions, 2 S. pennellii and 2 S. habrochaites cDNA libraries plus some individual cDNAs. Unigenes are the consensus sequences of these ESTs assembled with cap3 software as described in [17]. CDSs for the unigenes were calculated by SGN with ESTscan software [29, 30]. This software returns an optimum open reading frame based on Markov models and the nucleotide usage found in coding regions. Lukas Mueller at SGN kindly generated a custom list of all the unigenes and their constituent ESTs. The highest BLAST hit of every unigene versus the Arabidopsis thaliana genome was bulk queried and downloaded from SGN. Arabidopsis thaliana sequences were downloaded from TAIR (version 20080412, http://www.arabidopsis.org)

Assembly construction

We mined intraspecific SNPs in S. lycopersicum (L) ESTs and interspecific SNPs between S. lycopersicum and S. habrochaites (LxH) and between S. lycopersicum and S. pennellii (LxP) ESTs. Intraspecific analysis of S. pennellii and S. habrochaites ESTs were not performed as the EST libraries for those species were developed from a single accession.

For each unigene we used cap3 with relaxed parameters (-p 66 -b 99 -e 200) to assemble the EST sequences from the species participating in each of the analyses into analysis-specific contigs. Using these assemblies we calculated the sum and average of the qualities for each nucleotide call at every position and obtained analysis-specific consensus sequences. These consensus sequences were aligned to the original unigene sequence and the predicted CDS to estimate the beginning and end of the coding region. Every step in this process was carried out and verified using custom Perl scripts.

Perl and R [40] scripts were developed to estimate the quality score thresholds for sequence inclusion in SNP discovery. Once a threshold was determined we exclusively considered the positions in the assemblies above the threshold as follows. For interspecific analyses we required the sum of all the qualities of the nucleotides from each species to be over the threshold. For the intraspecific analysis we developed an algorithm that partitions all qualities at a given position into two groups with the minimum difference of the sums. We then considered only those positions in the assemblies at which the sum of the qualities of both groups was over the threshold. The threshold was defined as the quality score that maximizes the average sequence quality calculated separately for each species (interspecific analyses), or for each group of qualities (intraspecific analysis) taking into account only those positions over the threshold (Figure 2). We determined a quality threshold score of 50 for the intraspecific S. lycopersicum analysis and of 40 for both interspecific analyses (Figure 2).

SNP discovery

We defined SNPs in the intraspecific analysis as any assembly position where two and only two different bases were registered and for which the sum of qualities for each nucleotide call was over the quality threshold imposed for the analysis. For the interspecific analysis we considered SNPs whose positions within the assemblies presented a single and different nucleotide call for each species, and whose sum of qualities was greater than the imposed quality threshold. For quantification, SNPs that were consecutive in the assemblies were counted as a single polymorphisms. SNP rates were calculated as the number of non-contiguous SNPs per 100 bases.

Amino acid translation, SNP codon position and transitions/transversion ratios were evaluated with custom Perl and R scripts that analyzed the EST assemblies, the unigene sequence and the predicted CDS sequences.

Differential representation of SNPs in GO terms

Each nucleotide in the EST assemblies was assigned one or more GO categories based on the terms from the homologous Arabidopsis thaliana locus. To increase the accuracy of the test, we used only the parts of the assemblies that corresponded to coding regions and only those SNPs that has been predicted to produce amino acid changes. GO categories text files (version 20080712) were downloaded from TAIR [34]. For the interspecific analyses, we pooled the nucleotide calls for both (LxH and LxP) analyses. For duplicate loci we removed the locus with the shortest sequence. R scripts were developed to calculate over and under-representation of SNPs in nucleotide pools of each GO category versus all SNPs/nucleotides detected in the analysis using Fisher's exact test.

Tests for selection

We performed McDonald-Kreitman tests and estimated Maximum likelihood from codon-substitution models on those unigenes that contained S. lycopersicum together with S. pennellii or S. habrochaites sequences. For the McDonald-Kreitman test we built fasta files for each unigene with the estimated CDS for the wild species alleles and two S. lycopersicum alleles differing in the SNPs identified in the L analysis. To maintain the fidelity of the analysis, those positions that did not pass the quality threshold using the methods described above were substituted with the 'unknown' character and not considered in the subsequent analysis. The number of synonymous and non-synonymous substitutions and p-values were calculated using the previously described MK.pl Perl script [41].

Two maximum likelihood codon-substitution models were fit to test the hypothesis of existence of positive selection in each unigene [27, 28]. First, we fit the null model with a single d_n/d_s ratio with equal ratios in every branch. The second model allowed for two d_n/d_s ratios: one for the S. lycopersicum lineage and one for the rest of the tree. Then, a likelihood ratio test of the hypothesis of two branch rates was calculated by comparing the likelihood values from both models as in [42]. To do this we created Perl scripts that used ClustalW [43] to align the cultivated and wild species predicted CDS and the homologous Arabidopsis thaliana coding sequences. Alignments whose sum of qualities were not over the imposed quality threshold were substituted with the 'unknown' character. We also removed the parts of the alignments that lacked sequences from any of the three species. We constructed a phylogenetic tree of the three species and used PAML (v 3.14 [44]) to calculate the maximum likelihood of the models.

The resulting p-values from these experiments were corrected for multiple testing using the Benjamini and Hochberg algorithm in the Bioconductor multtest package [37, 45].

Molecular marker design

For each polymorphic unigene in the interspecific analysis we aligned its predicted protein sequence with the protein sequence of its Arabidopsis best BLAST hit using standalone BLAST and custom Perl scripts. We calculated intron positions in the unigene based on intron positions in the Arabidopsis CDS. For each SNP we designed primers using Primer3 [46] with the unigene sequences as input, adjusting the program to design the primers within the predicted exon where the SNP was located.

We used Bioperl to generate virtual PCR fragments for each SNP allele based on the primers designed, and to find restriction endonucleases that would differentially cut the fragments, thus creating molecular markers. A set of these molecular markers was tested on genomic DNA from S. lycopersicum VF36, S. pennellii LA716 and S. habrochaites LA1347. Touchdown PCR was performed in a MJ Research PTC-200 Thermocycler with a starting annealing temperature of 58°C, which decreased 0.5°C per cycle for 15 cycles and stayed constant at 55°C for 30 cycles. Extension time was 40 seconds and denaturizing steps were performed for 30 seconds at 96°C. PCR products were digested with the appropriate restriction enzymes to detect the polymorphisms. We developed a database and a website holding the molecular marker information for each interspecific SNP.

All R and Perl scripts are available by request.