Single nucleotide polymorphism discovery in rainbow trout by deep sequencing of a reduced representation library

Construction of Reduced Representation Libraries

Restriction enzymes Hpy CH4V and Bst UI did not entirely digest the rainbow trout DNA, therefore, they were not considered for further procedures. Trout DNA of an individual animal was completely digested with Alu I, Bst UI and Hae III. In contrast with bovine DNA[9], and similar to swine[12], repetitive DNA bands were not observed when separating the fragmented DNA in polyacrylamide gels with any of the three enzymes. The restriction enzyme Hae III (Figure 1) was selected for the library construction based on analysis of sequence data of 960 clones containing fragments from the approximately 440 bp size range: 94% of the sequenced fragments were unique and only 12% of the unique sequences matched other sequences when blasted against the NCBI salmonids nucleotide DNA database, indicating acceptably low content of repetitive DNA in the fragment pool.

After digesting the pooled DNA from 96 unrelated individuals with Hae III, the 440 bp fragments were selected for library construction (3.2 μg of fragment DNA was recovered from the preparative gel). An aliquot of the recovered fragments was cloned and Sanger sequenced before proceeding with the library construction protocols: 89% of 1,269 sequenced fragments were unique and 33% of the analyzed bases masked to the cGRASP salmon repeat masker database http://lucy.ceh.uvic.ca/repeatmasker/cbr_repeatmasker.py.

Sequencing results and SNP discovery

In order to discover SNPs present only in non-repetitive regions, the set of masked reads were mapped onto the two assemblies, with 75% mapping to the Newbler assembly and 63% mapping to the Velvet assembly. Both assemblies had similar distributions of coverage depth: ~40% less than 4×, ~10% greater than 10× and ~50% 4-10×. The average read depth at SNPs was skewed to the high end, ~9×. Despite the similarities, the two methods discussed thus far produced widely different predictions for the number of SNPs: 187,916 for Newbler and 52,942 for Velvet. The distribution of the SNPs across the contigs reveals that the vast majority of SNPs in the Newbler set came from contigs with 5 or more predicted SNPs. Filtering out such SNPs from both sets leaves 47,128 for Newbler and 33.624 for Velvet. Previous experience with SNP discovery in swine using a reduced representation library showed a much more reasonable distribution of putative SNPs across the assembled contigs [12]. One of the two attempts to lower the high rate of false positives was to use Velvet to create a different assembly. The other was to add an additional layering of filtering at the mapping stage (method N2 as described in the Methods section). N2 predicted 23,923 SNPs before filtering out the contigs with 5 or more SNPs and 22,022 after applying the filter. Results of the three analyses are summarized in Table 1.

A comparison of the SNPs from the three methods was done on a contig by contig basis, allowing for imperfect matches between velvet contigs and newbler contigs, as well as the occasionally varying numbers of predicted SNPs for the contigs. While compiling the lists of testable SNPs, some moderate filtering was performed to ensure sufficient quality and quantity of flanking sequence. As stated above, the two assemblies were more alike than different, yet the largest subgroup of SNPs is V1, indicating that the contigs unique to the Velvet assembly contain either a rich source of false positives or the best source of true SNPs. The next largest subgroup is N1, despite the fact that most of its potential SNPs were filtered out due to their presence in contigs with 5+ SNPs. The false positives seem to be infiltrating the unfiltered contigs as well. The third largest subgroup contains the SNPs predicted by all three methods, giving some assurance that the three approaches are similar enough to warrant comparison of validation rates as a means to incrementally improve the computational methods for SNP discovery. Even though the reads mapped in N2 were a subset (48%) of those mapped in N1, N2 contained SNPs not predicted in N1. These were from contigs that had 5 or more SNPs predicted by N1 but fewer than 5 predicted in N2. SNPs from each of the subgroups were tested and the validation results are discussed in the next section.

SNP validation

Call frequencies; minor allele frequencies and expected/observed heterozygosis of the validated SNPs are reported in Additional file 1. Call frequencies were higher than 90% for all markers and minor allele frequencies ranged from 1.3% to 49.7%. The double haploid samples were of great assistance for the validation process. Both Whale Rock samples were monomorphic for all validated markers; Skookumchuck, Hot Creek and Swanson were polymorphic for only one marker (OMS00033, OMS00163 and OMS00108, respectively) and Klamath for three markers (OMS00018, OMS00015 and OMS00108). The WGA DNA sample was successfully genotyped, 100% genotypes being equal to the original sample. Using SNaPshot, four SNPs were cross-validated. Genotyping results of the double-haploid fish were equal in both methods and 96% of the genotypes from the discovery samples were also equivalent.

Validated markers were also genotyped using the NCCCWA mapping families and 167 markers were polymorphic and mapped to the rainbow trout second generation map[21]. In Additional file 2, the position of the SNP markers in the linkage map is indicated. The positions are based on the two-point distance of the SNP markers to the closest linked marker in the map with the highest confidence in map location. The newly developed markers were distributed along all linkage groups of the rainbow trout map; there are at least two markers in each linkage group.

Construction and sequencing of Reduced Representation Libraries

Five different restriction enzymes (Alu I; Bst UI; Hae III; Hpy CH4V and Rsa I) having four-base recognition target sequences and producing fragments with blunt ends were tested for the library construction. Digestions of 1 mg of DNA from a single fish were performed for each enzyme as suggested by the manufacturer (New England Bio Labs). The digested DNA was separated on 3 mm thick 5% nondenaturing polyacrylamide gels. DNA fragments between 400 and 500 bp were excised from the gels (slices were cut every 2 mm or approximately 20 bp), fragments were eluted using a modified [31] "crush and soak" procedure and DNA was purified by alcohol precipitation. For trial Sanger sequencing to identify the best enzyme for reduced representation library construction, the fragments were ligated into the pSTBlue-1 vector of the Perfectly Blunt Giga Cloning kit (Novagen). DNA sequencing was carried out with an ABI Prism 3100 Genetic Analyzer. The ability to completely digest the trout genomic DNA and the number of DNA repetitive elements present in the selected fragment size range were evaluated for each enzyme.

Approximately equal amounts of DNA from 96 unrelated individuals from the NCCCWA 2005 and 2006 brood years [22] was pooled prior digestion with Hae III. Complete digestions of 1 mg of DNA from the pool (3 U enzyme/μg) were performed overnight at 37°C. DNA fragments were separated and eluted from the gel following the procedures previously described. Excised fragments approximately 440 bp in length were selected to construct the library; 600 ng of DNA obtained from the preparative gel was used for construction of the library and ligated to the sequencing adapters provided in the GS FLX library construction kit (Roche Applied Science). From this point, the library was prepared and sequenced on the GS FLX as per the manufacturer's protocols.

Sequence analysis and SNP discovery

Three independent SNP discovery approaches were used. The first two, which were called N1 and N2, were initiated by creating a reference sequence set by using the Newbler algorithm (version 1.1.03, provided with the GS FLX sequencer) to assemble all the unmasked sequence data into contigs. The individual reads were then screened for known repetitive elements and repeat masked using the cGRASP RepeatMasker database http://lucy.ceh.uvic.ca/repeatmasker/cbr_repeatmasker.py and sequences with at least 50 consecutive non-repetitive bases were mapped onto the assembly using ssaha2[32]. N2 differs from N1 in that the mappings for N2 were filtered so that only those with at least 90% sequence identity over at least 80% of the full length of the read were mapped onto the reference assembly. This step was done to decrease the likelihood of reads being mapped to duplicated regions of the genome. The third approach (V1) used the Velvet algorithm [33] for the assembly; in this case, all the sequence data was repeat-masked first, and the reads that contained at least 50 consecutive non-repetitive bases were used as input data. This same set of masked and filtered reads were then mapped back onto the Velvet assembly using ssaha2. The next step in all three approaches was to use ssaha_pileup [32] to detect sequence variation among the mapped reads. To be included as putative SNPs, both the major and the minor allele had to be detected at least twice each, without a third allele being observed. All SNPs from contigs that presented more than 5 putative SNPs (1 putative SNP/100 bp) were discarded, as a compromise balancing number of SNPs discovered and false discovery rate, since contigs with multiple apparent SNPs are more likely to represent paralogous loci (discarding loci with only two putative SNPs greatly decreases number of apparent SNPs but likely increases validation rate; allowing higher numbers adds additional SNPs but also increases false discover rate). Additionally, flanking sequence of at least 50 bp in one side (with a limit of 10% unknown bases) and 10 bp on the other side was required, in order to provide sufficient flanking sequence data to support assay design in the Illumina system to be used for genotyping.

SNP Validation

Illumina GoldenGate assays [34] were used to evaluate a fraction of the discovered SNPs. A set of 384 putative SNPs was randomly chosen (from among those with the highest scoring flanking sequences) to represent each one of the three sequence analysis methods (N1: 47; N2: 46 and V1: 100) and putative SNPs found in more than one analysis were also tested (N1 and N2: 47; V1 and N1: 47; V1 and N2: 47; V1, N1 and N2: 50).

Genomic DNA from 85 of the discovery panel fish and 98 samples from five NCCCWA mapping families [21] were used for the validation assays. In addition, seven double haploid samples from the Arlee, Hot creek, Klamath, Whale rock female, Whale rock male, Skookumchuck and Swanson clonal lines [35, 36] were genotyped and used as a control to eliminate the paralogous sequences that can erroneously be identified as SNPs. In some genotyping projects, the amount of DNA available for the assays is limited, and it is necessary to use Whole-Genome Amplified (WGA) DNA, as has already been used in GoldenGate assays [37]. A protocol based on the method described by [38] was used to prepare one WGA sample, which was added to the SNP validation panel, its respective template was also genotyped for comparison.

The BeadStudio software version 3.2 (Illumina) was used to analyze the Goldengate assays SNP data. Genotyping data quality was assessed following the guidelines suggested by Illumina [39]. The no-call threshold used was 0.25 and all the samples used for validation had at least a 0.85 Call rate score. A cross-validation assay was done using the ABI SNapshot method; four SNPs were randomly selected and the 96 discovery samples and the seven double-haploid individuals were genotyped.

A local BLAST search using the sequences of the validated SNPs against the rainbow trout transcriptome data from TIGR (BioEdit software at default settings http://www.mbio.ncsu.edu/BioEdit/BioEdit.html) was performed to identify coding regions.

Linkage analysis

Segregating data from the NCCCWA mapping families was used to anchor the polymorphic validated markers into the second generation rainbow trout linkage map [21]. Genotype data combined for both sexes were formatted using MAKEPED of the LINKAGE [40] program and checked for inconsistencies with Mendelian inheritance using PEDCHECK [41]. RECODE [42] and LNKTOCRI [43] were used to assemble the data into CRIMAP [44] format. SNP genotype data were added to that of Rexroad et al [21] and MULTIMAP [45] was used to conduct two-point linkage analyses to identify the closest markers from the published map which were ordered at the highest level of significance.