Development and Validation of Single Nucleotide Polymorphisms (SNPs) Markers from Two Transcriptome 454-Runs of Turbot (Scophthalmus maximus) Using High-Throughput Genotyping

The turbot (Scophthalmus maximus) is a commercially valuable flatfish and one of the most promising aquaculture species in Europe. Two transcriptome 454-pyrosequencing runs were used in order to detect Single Nucleotide Polymorphisms (SNPs) in genes related to immune response and gonad differentiation. A total of 866 true SNPs were detected in 140 different contigs representing 262,093 bp as a whole. Only one true SNP was analyzed in each contig. One hundred and thirteen SNPs out of the 140 analyzed were feasible (genotyped), while III were polymorphic in a wild population. Transition/transversion ratio (1.354) was similar to that observed in other fish studies. Unbiased gene diversity (He) estimates ranged from 0.060 to 0.510 (mean = 0.351), minimum allele frequency (MAF) from 0.030 to 0.500 (mean = 0.259) and all loci were in Hardy-Weinberg equilibrium after Bonferroni correction. A large number of SNPs (49) were located in the coding region, 33 representing synonymous and 16 non-synonymous changes. Most SNP-containing genes were related to immune response and gonad differentiation processes, and could be candidates for functional changes leading to phenotypic changes. These markers will be useful for population screening to look for adaptive variation in wild and domestic turbot.


Introduction
The turbot (Scophthalmus maximus; Scophthalmidae, Pleuronectiformes) is a commercially valuable flatfish that has been intensively cultured since the 1980s. Its production has steadily increased up to the present figure of 8549 tons in 2011 (91.2% European production from Spain; [1]) and it appears to be one of the most promising aquaculture species in Europe. In response to turbot industry demands, genetic markers have been developed in this species in order to evaluate genetic resources in both wild and hatchery populations and perform parentage analysis to support genetic breeding programs [2][3][4]. These markers have also been applied to develop genomic tools to identify genomic regions associated with productive characters [5][6][7] and to detect selection footprints in wild populations [8]. Increasing growth rate, controlling sex ratio (females largely outgrow males) and enhancing disease resistance currently constitute the main goals of genetic breeding programs in this species.
The necessity of understanding the immune response to pathogens of industrial relevance and to identify genes involved in the sex differentiation pathway led us to increase genomic resources in turbot. As a consequence thereof, an Expressed Sequence Tag (EST) database from cDNA libraries of the main immune tissues was constructed using Sanger sequencing [9]. Recently, this database has been amplified with two 454 FLX runs [10,11] (454-Life Sciences, Brandford, CT, USA; for 454-technique methodology see [12,13]). Next Generation Sequencing (NGS) technologies offer the ability to produce an enormous volume of data with a very low sequencing cost per base [12]. Thus, this turbot EST database is currently composed of ~70,000 unique sequences (~20,000 contigs and ~50,000 singletons). ESTs are essential to ascertain the gene [14,15], but also to identify polymorphic gene-associated markers, such as microsatellites and single nucleotide polymorphisms (SNPs) (type I markers; [9,[16][17][18]). Type I markers are very useful for constructing genetic or physical maps, and for comparative mapping [7,19,20].
SNPs have several advantages over other markers when it comes to mapping genes or inferring population structure [21]. They can be easily evaluated in silico off public databases and their genotypes quickly assessed by mini-sequencing reactions [9,22] or by high-throughput technologies [23,24]. SNP alleles are almost exclusively identical-by-descent (IBD) and thus they prevent scoring errors associated to homoplasy [25]. They are extremely stable, due to low mutation rates [26], and occur more often in the genome than other markers. In the human genome, for instance, there is on average 1 SNP per 300 bp [27], and their frequency in non-model species has been estimated at ~1 in 200-500 bases for non-coding DNA and ~1 in 500-1000 bases for coding DNA [28]. In turbot, Vera et al. [29] estimated 1 true SNP every ~100 bp from the EST database composed only of Sanger sequences, suggesting the existence of large SNP resources in this species. During the last decade, SNP discovery pipelines have been developed for aquaculture species including fish [18,[30][31][32][33][34][35], shellfish [36][37][38] and crustaceans [39,40]. In turbot, a SNP calling tool was included in the turbot database [9] and it has been refined in the updated version [11]. In this study, we screened genomic resources available in an updated version of the turbot EST database using contigs containing NGS 454-sequences to identify and characterize SNPs associated to immune-and reproduction-related genes. These markers will be used for further structural genomic analysis focused on quantitative trait loci (QTLs) linked to productive traits, as well as for population screening to look for adaptive variation in wild and domestic turbot.

Database Exploitation and SNP Detection
The main characteristics of the turbot 454-transcriptome sequencing runs have been described in previous studies [10,11]. The used database (version 4.0 September 2011) was constituted by 71,033 unique sequences, 18,880 contigs and 52,153 singletons including 454-sequences and Sanger sequences [9] with a total length of 52,402,177 base pairs (bp, ~52 Mb). However, in order to avoid duplicates with the previous SNPs developed from sequences obtained with Sanger methodology [29], and since we were mainly interested in validating SNPs at new immune-and reproduction-related genes, only contigs composed exclusively of at least four 454-sequences were used for SNP detection. Thus, 140 contigs from the turbot database, which met these requirements, were taken into account for the SNP development. The total length analyzed was 262,093 bp and contig length ranged from 728 bp to 4885 bp, with a mean length value of 1872.09 ± 746.69 bp. The total number of true SNPs detected using the program QualitySNP (for true SNP definition see the experimental section) was 866, SNP number per contig ranged from 1 to 58, with a mean value of 6.18 ± 8.34. Thus, the expected frequency of SNP appearance in the analyzed sequences would be 1 SNP every 302 bp. This value is lower than that previously reported in S. maximus (1 SNP each ~100 bp; [29]), but similar to those described in non-model species [28]. The success of any genotyping method is reflected in what is referred to as the conversion rate and the global success rate. The former only considers the polymorphic markers, whereas the latter considers all the markers (monomorphic and polymorphic) that were successfully typed within the analyzed samples [41]. Of the 140 true SNPs tested, 27 (19.3%) could not be genotyped, and thus they were considered to be genotyping failures due to technical and/or genotyping problems. Only 2 out of the 113 feasible SNPs (see definition in the experimental section) were monomorphic. Therefore, the global success rate and conversion rate were 80.7% and 79.3%, respectively. Global success rate was very similar to that previously described in the species (78.4%), but conversion rate was much higher than previously reported using sequences from cDNA libraries (37.7%; see Vera et al. [29]), likely due to the different library construction methods and bioinformatic pipeline approaches followed in 454 and Sanger contigs (see experimental section).

SNP Performance
A total of 65 transitions (A/G and C/T) and 48 transversions (A/C, A/T, C/G and G/T) were detected among feasible SNPs, A/G being the most common (35) and A/C the least common (6) substitutions observed (Figure 1). This represented a transition/transversion (ts/tv) ratio of 1.354. This ratio was lower than that observed by Vera et al. [29] (1.885) and in silico (1.456) by Pardo et al. [9], but it was very similar to that described for common carp (Cyprinus carpio) (1.310) [42] and gilthead seabream (Sparus aurata) (1.375) [31]. Also, the most frequent transitions and transversions differed from previous reports: C/T and G/T, respectively [29], and A/G and A/C [9]. These discrepancies could be due to the opposite sequencing directions, as all sequences by Vera et al. [29] and Pardo et al. [9] were obtained from the 3' end using cDNA libraries, while those from the 454-run were randomly obtained by fragmentation of the whole cDNA according to the cDNA rapid library preparation method (Roche Farma, S. A. [43]). Moreover, the longer coding region portion analyzed in 454-runs regarding Sanger sequencing in our study may determine differences because of the different selective constraints of UTR regarding coding regions. No differences were detected among distribution of the variants between tested SNPs and feasible SNPs (χ 2 = 0.3115; p = 0.9974). All polymorphic SNP loci showed two alleles and all of them agreed with those expected from the database information.

SNP Diversity
Only two loci among the 113 feasible SNPs were monomorphic (SmaSNP_287 and SmaSNP_334). Among polymorphic SNPs, unbiased gene diversity (He) estimates ranged from 0.060 at SmaSNP_237, SmaSNP_245 and SmaSNP_305 to 0.510 at SmaSNP_225 with a mean value of 0.344 ± 0.149. The minimum allele frequency (MAF) in the polymorphic markers ranged from 0.030 (SmaSNP_237, SmaSNP_245 and SmaSNP_305) to 0.500 in SmaSNP_249 with a mean value of 0.259 ± 0.140. Departures from Hardy-Weinberg equilibrium (HWE) were detected in five markers (SmaSNP_253, SmaSNP_271, SmaSNP279, SmaSNP_289, SmaSNP_326; Table 1), although all markers were at equilibrium after Bonferroni correction (p = 0.0004). The samples from the Cantabrian turbot population were globally in accordance with HWE expectations when tested simultaneously for all loci (p = 0.9999). These polymorphic values were in the range to those previously described in the species [29], and they were also similar to those reported in other fish species [42,44]. No Linkage disequilibrium (LD) was detected among the 6328 loci pairs after Bonferroni correction (p = 0.0004).

SNP Position within Genes: Synonymous vs. Non-Synonymous Substitutions
Consensus sequences of contigs containing polymorphic SNPs were compared using NCBI BLAST with public databases, namely UniRef90, NCBI's nr, KEGG, COG, PFAM, LSU and SSU. The subsequent BLAST output was then parsed with Auto FACT [45]. All contigs containing feasible SNPs were annotated (except SmaSNP_320, Table 1). The informative strand, reading frame, and stop codon at each contig were recorded using homology with the highest homologous annotated sequence in public databases. Nine feasible SNPs (8.0%) could not be positioned, because no consistent reading frames were detected (indicated as "unknown" location on Table 2). Fifty-five SNPs (48.7%) were located in untranslated regions (UTR), either in the 5' UTR (17, 15.0%) or 3' UTR (38, 33.6%), which is in accordance with the approximately double length of 3' compared to 5' UTR [9]. On the other hand, 49 SNPs (43.4%) were localized in the coding region (Table 2), a percentage of SNPs higher than previously reported in the species (24.7%, [29]) and in other aquaculture fish species (e.g., Atlantic salmon 24%, [32]; Atlantic cod 17.4%, [34]). All these studies followed a 3' UTR Sanger sequencing strategy, and therefore the coding region was less represented than in the case of the 454 Roche runs after a cDNA rapid library preparation protocol, which accounts for the differences observed. This result shows the utility of the NGS methodologies for SNP detection in the coding region. Thirty-three (29.2%) of these 49 SNPs were synonymous, and 16 (14.2%) were non-synonymous. On the other hand, the relationship between synonymous vs. non-synonymous changes (2:1) was lower than in other species [46,47]. Evolutionary constraints should preferentially eliminate non-synonymous variation because it is usually associated with deleterious mutations [35].
Non-synonymous SNPs in coding regions represent alternative allelic variants of a gene, which can determine functional changes in the corresponding proteins and lead to phenotypic changes. Among these genes there can be found a retinol dehydrogenase (SmaSNP_264), three zona pellucida proteins (SmaSNP_212, SmaSNP_217, SmaSNP_282) related to reproduction processes, and a lipocalin (SmaSNP_325) involved in tear secretion ( Table 2).
In the present study, we used sequences obtained from two transcriptome 454-pyrosequencing runs, one related to immune system [10] and another one from the hypothalamic pituitary-gonad axis [11]. Thus, GO terms were mainly related to immune response and reproduction processes ( Table 2). The non-synonymous variation was associated with genes involving either immune response or sex differentiation processes. A large number of SNP linked to annotated genes were identified and validated. This set of markers are being used for population genomic studies and turbot genetic map enrichment, both approaches providing useful information for evolutionary and turbot industry applied studies.

EST Database, SNP Detection and Annotation
Sequences were obtained from two transcriptome 454-pyrosequencing runs of turbot cDNA libraries, one belonging to the immune transcriptome [10] and another one from the hypothalamic pituitary-gonad axis [11]. A brief description of both runs is shown in Table 3. All the 454-reads were assembled with MIRA [48], and they make up the 454-sequences incorporated into the turbot database. In order to create contigs and detect SNPs, these 454-sequences were assembled alongside Sanger sequences available [9] in the database with CAP3 [49] using default parameters. This is a common strategy when dealing with hybrid Sanger-454 assemblies [50]. The resulting ACE format assembly file was fed into QualitySNP [51] in conformity with the bioinformatic pipeline described by Vera et al. [29]. Briefly, QualitySNP uses three filters for the identification of reliable SNPs: Filter 1 screens for all potential SNPs with the requirement that every allele is represented in more than one sequence (each contig has to have at least a depth of 4 sequences); filter 2 uses a haplotype-based strategy to detect reliable SNPs after reconstructing confident haplotypes with an algorithm that minimizes false haplotypes due to the occurrence of sequencing errors; and filter three screen SNPs by calculating a confidence score based on sequence redundancy and quality (only sequences with PHRED >20 were used). SNPs that pass filters 1 and 2 are called real SNPs and those passing all filters are called true SNPs [51]. Resulting files were processed with our own custom Perl programs to extract relevant information. The obtained data were imported into a mySQL server [52]. A user-friendly web access interface was designed so that contig graphs are clickable and the output visually refined with color-coded nucleotide views [53]. A graphical interface allowing for SNP database search by alleles, contig depth, and annotation was set up. EST annotation of these contigs was performed using BLASTx, which searches proteins using a translated nucleotide query [54]. Only E-values lower than 10 −5 were considered for gene annotation ( Table 1, Table S1).

SNP Genotyping and Validation
DNA of all individuals analyzed was extracted from a piece of caudal fin using standard phenol-chloroform procedures [55].
SNPs identified were validated and genotyped with the MassARRAY platform (Sequenom, San Diego, CA, USA) following the protocols and recommendations provided by the manufacturer. Briefly, the technique consists of an initial locus-specific polymerase chain reaction (PCR), followed by single-base extension using mass-modified dideoxynucleotide terminators of an oligonucleotide primer that anneals immediately upstream of the polymorphic site (SNP) of interest (see [56,57] for more technical information). The distinct mass of the extended primer identifies the SNP allele. Primer sequences, SNP position, expected variants and annotation for the 140 tested SNPs are shown on Supplementary Table 1. MALDI-TOF mass spectrometry analysis in an Autoflex spectrometer was used for allele scoring.
Assays were designed for 140 true SNPs always located in different sequences and were combined in 7 multiplex reactions including 24 SNPs each except for multiplex 5 (23 SNPs), 6 (18 SNPs) and 7 (3 SNPs) (see Supplementary Table 1 for multiplex information). SNP multiplexes were designed in silico and tested on a panel of 8 turbot individuals from a wild Cantabrian (northern Spain) population. SNPs were classified based on manual inspection as "failed assays" (in case that the majority of genotypes could not be scored and/or the samples did not cluster well according to genotype), and feasible SNPs (markers with proper and reliable genotypes), these being either monomorphic or polymorphic.

Gene Diversity and Population Analysis
In order to estimate genetic diversity parameters, all SNPs were genotyped for polymorphism evaluation in a sample of 33 individuals (including the 8 individuals used for marker performance) from the wild Cantabrian population previously used.
Estimates of genetic diversity (unbiased expected heterozygosity (He) and minimum allele frequency (MAF)) were estimated using FSTAT 2.9.3 [58]. The conformance to Hardy-Weinberg (HW) and genotypic equilibria were obtained using GENEPOP 4.0 [59,60]. Conformance to HWE was checked using the complete enumeration method [61] because only two alleles were detected at each locus. Bonferroni correction was applied when multiple tests were performed [62].

Detection of Synonymous/Non-Synonymous SNPs
All the six possible reading frames of the consensus sequence of each containing SNP functionally annotated contig were obtained using ORF (Open Reading Frame) Finder application [63]. The best candidate frame (usually the longest one) was compared against the NCBI protein database using BLASTp and BLASTx, and the protein with highest E-value was downloaded and aligned with the selected frame for SNP location using Clustal W [64] implemented in BioEdit v. 7.1. [65]. This approach enabled us to locate SNPs by looking at the coding region. For those SNPs in the coding region, the resulting amino acid sequences of both variants were translated to determine whether SNP variants were synonymous or non-synonymous. Gene onthology (GO) terms were searched using QUICKGO [66] and AmiGO [67] utilities.

Conclusions
A total of 140 contigs (total length 262,093 bp) formed exclusively by 454-pyrosequencing reads were used to identify new putative SNPs in S. maximus. One hundred and thirteen SNPs of the 140 tested were amplified and genotyped, 111 being polymorphic in a wild Cantabrian population, showing the utility of the new NGS techniques for true SNP detection (conversion rate = 79.3%). Diversity levels at the population were similar to previous studies [29,42,44] and were in accordance with HWE expectations. An important number of these polymorphic SNPs were located in the coding region and 16 of them (14.4%) represented non-synonymous changes at genes related to immune response and gonad differentiation processes as shown by the detected GO terms. Therefore, these SNPs are valuable resources for future population genetics, high-resolution genetic maps, quantitative trait loci (QTL) identification, association studies and marker assisted selection (MAS) breeding in turbot.