Establishment of a pipeline to analyse non-synonymous SNPs in Bos taurus

Single nucleotide polymorphisms (SNPs) are an abundant form of genetic variation in every species. SNPs are more stable than other genetic markers, such as simple sequence repeats although they are not as informative due to their predominantly biallelic nature. They are distributed throughout the genome and in Homo sapiens SNPs have been estimated to occur at a frequency of approximately one in every 1,000 bp [1] making them the most common form of genetic variation. SNPs are extremely useful for association studies, gene mapping and phylogenetic studies. They provide genomic landmarks of past events shaping individuals, breeds and species. In species of agricultural importance, SNPs can be used to identify regions influencing important economic traits. Single nucleotide polymorphisms have also been shown to be responsible for a number of economically desirable phenotypes in sheep, such as the muscle hypertrophy phenotype callipyge and the fertility phenotype Inverdale [2, 3]. A public SNP database (dbSNP) has been established by the National Center for Biotechnology Information to store sequence variation information from a variety of species [4]. There are currently 19,395 Bos taurus SNPs deposited in this database.

SNPs can be identified in coding regions (cSNPs) or non-coding regions. A subset of SNPs in coding regions give rise to variation in the amino acid sequence of the encoded protein and are known as non-synonymous SNPs (nsSNPs). Such SNPs have been reported to occur less frequently than synonymous SNPs, presumably due to evolutionary constraints as selection eliminates deleterious substitutions from the population [5]. Non-synonymous SNPs are therefore of interest as they are more likely to affect the function of the encoded protein and so may influence a phenotype of interest. Indeed, it has been estimated that 20–30% of nsSNPs affect protein function [6, 7].

SNPs are very amenable to high-throughput analysis but in vitro SNP discovery can be a lengthy and expensive process. Identifying SNPs in silico and then validating the polymorphisms in vitro provides a cost effective method of SNP discovery for mapping or association studies.

We developed a pipeline for detecting cSNPs and validating these SNPs cost-effectively. Available cattle expressed sequence tagged (EST) sequences were mined to identify cSNPs yielding 15,353 putative transcribed SNPs, of which 6,325 were predicted to be in coding regions. 2,868 of the cSNPs were predicted to be non-synonymous. A subset of the non-synonymous SNPs was genotyped through pools containing DNA from a number of cattle of different breeds and 29% of these SNPs were confirmed to be polymorphic in the animals tested (minor allele frequency >10%). Finally association was calculated between a subset of the validated non-synonymous SNPs and a number of economically important phenotypes.

In order to detect non-synonymous SNPs 528,218 public bovine expressed sequences were assembled into contigs and polymorphisms in these contigs detected using conservative criteria. In total the region examined encompassed over 11,000 kbp of bovine expressed sequence and resulted in the detection of 15,353 putative SNPs in 7,145 contigs. The sequence of the contigs is given in Additional file 1 and the SNP allele frequencies are listed in Additional file 2. 6,325 of the SNPs were predicted to be coding SNPs and 2,868 of the cSNPs were predicted to result in a change in the amino acid encoded at that position. The severity of the change of amino acid was calculated using an index that combined the scores from a Blosum90 and PAM250 substitution matrix. The putative synonymous and non-synonymous cSNPs are given in Additional file 3. If all base substitutions were equally likely, 76% of the cSNPs would be expected to result in an amino acid change. We found 45.3% of the cSNPs were predicted to be non-synonymous indicating that there is strong selection pressure against amino acid changes in the sequences studied. The ratio of non-synonymous SNPs to synonymous SNPs reported here is in good agreement with the findings of groups who have calculated this ratio from human sequences [5, 8].

The frequency of base substitutions at each of the codon positions for non-synonymous and synonymous SNPs was calculated and is shown in Figure 1. As expected the majority of the substitutions (93.5%) occurred in the third position for synonymous SNPs and at the first and second positions (87.1%) for non-synonymous SNPs. Candidate SNPs were also categorized according to nucleotide substitution as either transitions (A/G or C/T) or transversions (A/C, A/T, C/G, or G/T). Of the 15,353 putative SNPs detected, 15,335 SNPs could be categorized in this manner with the remaining 18 showing 3 alleles at the SNP position. 66.2% of the SNPs detected were the result of transitions, even though transitions represent only one third of the potential types of mutation. A high frequency of transitions has been observed in other SNP discovery programs and the transition frequency observed here is in good agreement with human data where approximately 70% of SNPs were found to be the result of transitions [9–11]. We found the frequency of transitions was significantly higher among synonymous SNPs (81.9%) than non-synonymous SNPs (51.6%; χ² = 665.8; P = 8.0 × 10^-147). The frequency of the various substitutions is illustrated in Figure 2. Within the non-synonymous SNPs, the frequency of transitions between conservative substitutions (index score ≥ 2) and radical substitutions (index score ≤ -2) was also examined. In this case no significant difference was found between the two groups (43.9% and 43.8% respectively; χ² = 0.33; P = 0.97). A high frequency of transitions has been observed in other SNP discovery programs [11] and reflects the high frequency of C to T mutation after methylation.

For each contig, the proportion of nsSNPs was plotted against SNP density and a polynomial regression fitted, weighting by the total number of SNPs the contig contained (Figure 3). Only the linear term was found to explain any effect of SNP density, and this accounted for less than 1% of the variation in the proportion of nsSNPs. Therefore, although the proportion of SNPs that are non-synonymous increases with SNP density, SNP density is a poor predictor of this proportion. The nsSNPs detected in this study appear to come from contigs with a range of SNP densities and are not substantially more frequent in contigs with high sequence variation.

In order to determine the percentage of in silico derived non-synonymous SNPs that were polymorphic, 920 nsSNPs were genotyped through pools of cattle DNA using chip based MALDI-TOF mass spectrometry. The pools used consisted of either a single pool containing DNA from 57 animals of a variety of breeds (see methods) or four pools each containing DNA from eight animals from a total of seven breeds plus a back-cross containing two breeds (see Table 1). The pooling strategy provided power to reliably identify SNPs with a minor allele frequency (MAF) greater than 10% and may detect SNPs with a MAF greater than 3%. In total, 795 of the 920 putative nsSNPs tested could be reliably scored giving a strike rate of 86%, while 125 of the assays failed. In total 227 of the putative nsSNPs were confirmed as polymorphic in the animals tested, while 568 of the putative nsSNPs were not confirmed in the animals tested indicating that they are not polymorphic, are rare (MAF <10%), or polymorphism is specific to breeds not tested in this study. The number of nsSNPs that were polymorphic in the pools tested is illustrated in Figure 4. In total 19 SNPs were polymorphic in 1 pool, 25 were polymorphic in 2 pools. 44 were polymorphic in 3 pools and 139 were polymorphic in all 4 pools tested. These 139 SNPs are found in 131 contigs, of which 128 could be annotated with a human RefSeq, giving 112 unique genes. These SNPs represent a subset of highly polymorphic non-synonymous SNPs and include a number of SNPs in immune response genes, such as Major Histocompatibility Complex genes. These genes are known to be extremely polymorphic. The full list of nsSNPs tested and their allele calls is available in Additional file 4.

Pool	Breed	Number of animals	Total number of animals
Pool A	Belgian Blue	8	16
	Piedmontese	8
Pool B	Murray Grey	8	16
	Shorthorn	8
Pool C	Hereford	8	16
	Simmental	8
Pool D	Friesian	8	16
	Jersey × Limousin back-cross	8

As there were a large number of nsSNPs that were polymorphic in all pools tested we needed to eliminate the possibility that this was due to spurious amplification of a closely related gene with a different nucleotide at the extension position. We therefore genotyped 55 nsSNPs, which we scored as polymorphic in all four breed pools, across four individual animals. Only one of these SNPs was heterozygous in all four individual animals tested indicating that amplification of closely related genes is not a serious contaminant in our cSNP discovery program. Approximately 0.8% of the candidate nsSNPs detected in this study were previously deposited in dbSNP. Of the 227 polymorphic nsSNPs validated in this study, 3 had been previously deposited in dbSNP. We have deposited all 15,353 putative SNPs detected in this study in dbSNP.

It was found that the allele frequencies of some of the genotyped SNPs were not in Hardy-Weinberg equilibrium. These SNPs are highlighted in bold in Table 2. The reason these SNPs do not obey Hardy-Weinberg equilibrium is unknown, however, it may reflect the fact that the effective population sizes are relatively small, the animals were not randomly mated and the animals genotyped were not a random sample but were chosen as they were outliers in the population.

Using public bovine expressed sequence data we have identified 15,353 candidate SNPs in the genome of Bos taurus. One putative SNP was detected approximately every 716 bp for the sequences studied. 6,325 of the putative SNPs are predicted to be coding sequence SNPs and 2,868 of the cSNPs are predicted to be non-synonymous. The low frequency of non-synonymous SNPs detected, with respect to the frequency expected if all base substitutions in coding regions were equally likely, implies that the genes studied are under strong selection pressure. Recently, the deleterious mutation rate was estimated by comparing the coding sequences of many genes between Bos taurus and Homo sapiens. This analysis concluded that 87% of amino-acid changing mutations were deleterious [13]. Therefore, any non-synonymous SNPs detected in this study may affect the function of the encoded protein. The non-synonymous SNPs we detected were in a variety of genes with differing biological functions and on a number of different bovine chromosomes and may influence traits of economic interest. Additional programs such as SIFT or SNPs3D may be used to predict whether the amino acid substitution will affect protein function [14, 15]. However, it is likely that direct association or linkage studies will be needed to determine the real effect of allele substitution on phenotype in an animal or population.

Polymorphisms that affect biological function can also occur outside coding sequences [16], particularly in the proximity of genes such as in splice junctions and promoters. As all of the SNPs detected in this study are in the vicinity of genes, any of these SNPs could potentially influence phenotype or be in linkage disequilibrium with a phenotype-altering mutation. Therefore all SNPs detected in this study are useful for association studies although non-synonymous SNPs are of particular interest.

Approximately 66% of the polymorphisms detected in this study were the result of transitions. Transitions represent only one third of the possible types of substitutions yet they are reported to account for approximately 70% of all substitutions. Our data is in good agreement with previous studies that have calculated the frequency of transitions in human genes [9, 10]. The fact that our transition frequency is so close to the expected value implies that our SNP detection strategy is accurate. A high false positive SNP discovery rate would be expected to result in a high frequency of transversions.

The percentage of nsSNPs that were polymorphic was determined by genotyping 920 putative SNPs through pooled DNA from cattle of different breeds. 29% of the SNPs were found to be polymorphic in the samples tested. However, the pooling strategy only allowed reliable detection of SNPs with a minor allele frequency >10%. Therefore many more of the SNPs may be polymorphic, but with a lower MAF. The SNPs were also genotyped through a limited number of cattle breeds, therefore, some SNPs may be present or polymorphic in breeds not tested in this study. Also, many of the ESTs used for SNP detection were derived from cDNA from a limited number of animals, therefore some of the SNPs detected may be restricted to particular animals. Additionally, high-throughput EST sequencing is prone to sequencing errors and this may result in some false positive SNPs. A previous study identifying bovine SNPs from EST data reported a SNP validation rate of 38% for non-synonymous SNPs which had been detected more than once in EST data, although they validated only a small number of SNPs [17]. Therefore, our results are consistent with these findings.

Approximately 0.8% of the non-synonymous SNPs predicted from this study were previously identified and deposited in NCBI dbSNP with the remainder representing novel bovine SNPs. Many of these SNPs could potentially be used as genetic markers in known genes. All SNPs detected in this study and their allele frequencies have been deposited in dbSNP.

Forty one validated polymorphic SNPs in 33 genes of interest were genotyped through animal resources which had been measured for age of puberty, retail beef yield or resistance to facial eczema. The number of nominally significant associations and the strength of these associations with facial eczema resistance was higher than expected by chance although no single SNP reached a genome wide level of significance. The most significant SNP was in the gene B cell lymphoma 10 (BCL10) and showed association with facial eczema resistance in a Jersey population. BCL10 encodes a protein that appears to both activate and be activated by NFκB and overexpression of BCL10 causes weak apoptosis mediated via an amino terminal caspase recruitment domain (CARD) [18, 19]. The BCL10 protein is a positive regulator of lymphocyte activation and proliferation [20]. It is more likely that this gene is in linkage disequilibrium with a resistance conferring polymorphism in Jerseys, than responsible for facial eczema resistance per se, as the SNP in this gene was not associated with facial eczema resistance in Friesians.

We have developed a pipeline for detecting, validating and genotyping large numbers of single nucleotide polymorphisms using Bos taurus as a model organism. However, this strategy could be successfully applied to any organism for which sufficient EST sequence data exists. We have detected 15,353 putative SNPs in the vicinity of genes in this study and 6,325 candidate cSNPs, of which 2,868 were predicted to be non-synonymous. Approximately 0.8% of the non-synonymous SNPs detected in this study were previously deposited in dbSNP while 99% appear to be novel bovine SNPs. These novel SNPs have the potential to add greatly to the previously discovered bovine SNPs although some additional SNP validation work is required. These SNPs can be used for association studies or gene mapping and we have found SNP associations for resistance to facial eczema in Jersey or Friesian populations at a higher rate than expected by chance, although no single SNP reached a genome wide level of significance. Whether these genes are involved with facial eczema resistance in specific breeds, false positives, or are merely in linkage disequilibrium with a resistance enhancing polymorphisms remains to be determined.