Gene-based single nucleotide polymorphism discovery in bovine muscle using next-generation transcriptomic sequencing

Background Genetic information based on molecular markers has increasingly being used in cattle breeding improvement programmes, as a mean to improve conventionally phenotypic selection. Advances in molecular genetics have led to the identification of several genetic markers associated with genes affecting economic traits. Until recently, the identification of the causative genetic variants involved in the phenotypes of interest has remained a difficult task. The advent of novel sequencing technologies now offers a new opportunity for the identification of such variants. Despite sequencing costs plummeting, sequencing whole-genomes or large targeted regions is still too expensive for most laboratories. A transcriptomic-based sequencing approach offers a cheaper alternative to identify a large number of polymorphisms and possibly to discover causative variants. In the present study, we performed a gene-based single nucleotide polymorphism (SNP) discovery analysis in bovine Longissimus thoraci, using RNA-Seq. To our knowledge, this represents the first study done in bovine muscle. Results Messenger RNAs from Longissimus thoraci from three Limousin bull calves were subjected to high-throughput sequencing. Approximately 36–46 million paired-end reads were obtained per library. A total of 19,752 transcripts were identified and 34,376 different SNPs were detected. Fifty-five percent of the SNPs were found in coding regions and ~22% resulted in an amino acid change. Applying a very stringent SNP quality threshold, we detected 8,407 different high-confidence SNPs, 18% of which are non synonymous coding SNPs. To analyse the accuracy of RNA-Seq technology for SNP detection, 48 SNPs were selected for validation by genotyping. No discrepancies were observed when using the highest SNP probability threshold. To test the usefulness of the identified SNPs, the 48 selected SNPs were assessed by genotyping 93 bovine samples, representing mostly the nine major breeds used in France. Principal component analysis indicates a clear separation between the nine populations. Conclusions The RNA-Seq data and the collection of newly discovered coding SNPs improve the genomic resources available for cattle, especially for beef breeds. The large amount of variation present in genes expressed in Limousin Longissimus thoracis, especially the large number of non synonymous coding SNPs, may prove useful to study the mechanisms underlying the genetic variability of meat quality traits.


Background
Cattle (Bos taurus) are considered to have been one of the first animals domesticated by man for agricultural purposes. Approximately 10,000 years ago, cattle ancestors (aurochs) were tamed to provide milk, meat and hides and for draft purposes [1]. Bos taurus was also one of the first animal species to enter the genomics era. In the past few years, genetic information based on molecular markers has increasingly been used in cattle breeding improvement programmes, as a mean to improve conventionally phenotypic selection, particularly for traits with low heritability or for which measurement of phenotype is difficult, expensive, only possible late in life, sex-limited or not possible on selection candidates [2]. Advances in molecular genetics have led to the identification of several genes or genetic markers associated with genes that affect economic traits [3][4][5][6][7][8][9][10]. For example, the non conservative K232A substitution in the acylCoA:diacylglycerol acyltransferase (DGAT1) gene has a major effect on milk yield and composition [5]. Several of these genetic markers are now available and used in industry marker-assisted selection programmes [11,12].
Because of its economical importance Bos taurus was one of the first mammals to have its genome sequenced. In August 2006, the sequence of the cattle genome was released by the Human Genome Sequencing Center at Baylor College of Medicine [13]. During the sequencing more than 2.2 million putative single nucleotide polymorphisms (SNPs) were identified and deposited in public databases [14]. The Bovine Genome Sequencing Consortium has since discovered approximately 62,000 extra highquality SNPs [15]. These SNPs have been used to develop a whole-genome cattle SNP genotyping microarray [16]. More recently, a novel higher-density whole-genome bovine SNP BeadChip, containing~770,000 SNPs has being developed by Illumina [17].
With the availability of genome-wide dense marker maps and cost-effective genotyping methods, a novel genetic improvement method, called genomic selection, has been developed and is already revolutionising the cattle breeding industry. Genomic selection is a form of marker-assisted selection in which genetic markers covering the whole genome are used to estimate breeding values (genomic breeding values) [18]. However, since most of the SNPs present on the whole-genome cattle SNP genotyping microarrays commonly used, are not in genes and also because of the extent of linkage disequilibrium, SNPs associated with economically important traits, will most likely, not be involved directly in these traits. The identification of the causative genetic variants involved in the phenotypes of interest, remain a difficult task. It is therefore, crucial to develop strategies to pinpoint more rapidly causative genetic variants underlying phenotypes of interest.
The identification of these causative genetic variants, also known as quantitative trait nucleotides (QTNs) involves the mapping of quantitative trait loci (QTLs), the discovery of novel genetic markers in the QTL regions, the fine-mapping of QTLs and then the sequencing of candidate genes. This iterative process until recently was very time-consuming, but thanks to the availability of a large number of SNPs and to the relatively low-cost of whole-genome genotyping methodologies, the finemapping of QTL regions has now been expedited. In addition, the advent of novel sequencing technologies [19][20][21][22][23] offers now a new opportunity for the identification of QTNs, with the ability to partially or completely re-sequence mammalian genomes, in a relatively costeffective manner, and to identify polymorphisms responsible for the traits of interest.
The genome of animals from many species has now been sequenced, including the genomes of several bulls [24][25][26][27][28][29][30]. For example, Eck et al. (2009) generated the first single cattle genome sequence by a next-generation sequencing method [24]. By sequencing the wholegenome sequence of one Fleckvieh bull, they discovered more than 2 million novel cattle SNPs. Even though sequencing costs plummeting, sequencing whole-genomes or large targeted regions is still too expensive for most laboratories.
A whole-transcriptome RNA sequencing (RNA-Seq) method has recently been developed to identify and quantify genes expressed in different tissues [31,32]. This method has also been used to identify polymorphisms in transcribed regions, in different species, including in cattle [33,34]. A transcriptomic-based sequencing approach offers a cheaper alternative to identify a large number of polymorphisms and possibly to discover QTNs.
In the present study, we performed a gene-based SNP discovery analysis in bovine Longissimus thoraci, using a whole-transcriptome sequencing approach. To our knowledge, this represents the first study done in bovine muscle. For this purpose, muscle samples from three different Limousin bulls were analysed. We have identified more than 34,000 putative SNPs, including more than 60% novel polymorphisms. To evaluate the accuracy of the SNPs detected, 48 putative SNPs were genotyped. One-hundred percent concordance was observed when a stringent SNP quality criterion was chosen. The RNA-Seq data and the collection of newly discovered coding SNPs improve the genomic resources available for cattle, especially for beef breeds. The large amount of variation present in genes expressed in Limousin Longissimus thoracis, especially the large number of non synonymous coding SNPs, may prove useful to study the mechanisms underlying the genetic variability of meat quality traits.

RNA sequencing
To obtain a global view of the bovine Longissimus thoracis transcriptome at single-nucleotide resolution, poly(A)-enriched mRNA from three Limousin bull calves were retrotranscribed and subjected to highthroughput sequencing. The three RNA-Seq libraries were barcode-tagged and sequenced on one lane of an Illumina HiSeq2000 sequencer. Sequencing of cDNA libraries generated a total of 125,781,357 raw paired-end reads with a length of 100 bases, resulting in a total of 25 gigabases. The reads were de-multiplexed to assign reads to each sequenced sample according to its barcode index. Approximately 36 to 46 million paired-end reads were obtained for each library. Reads from each sample were then mapped back to the bovine reference transcriptome. We used the set of Bos taurus Ensembl transcripts v61 RefSeq genes as the reference transcriptome. This set contains transcripts for 22,915 known or novel genes but also pseudogenes. Based on mappings done using the Burrows-Wheeler Aligner (BWA) programme, 63% to 67% of the mapped reads were aligned properly paired (Table 1). Transcriptome contamination was negligible (0.19%-0.24%). A total of 19,752 transcripts (16,287 genes) were identified, with at least one paired-end read in all samples analysed. Similar RNA-Seq read mapping rate and the number of genes identified were obtained in other RNA-Seq bovine studies [33][34][35][36][37][38]. For example, Wickramasinghe et al. (2012) found that~65% of the RNA-Seq reads they generated while sequencing the milk transcriptome mapped uniquely onto the bovine genome. They also found that~17,000-19,000 genes were expressed in milk [35]. Baldwin and collaborators found, this time, by sequencing the rumen epithelium that~71% of the reads mapped onto~17,000 different genes [36].
Gene expression was normalised as paired-end reads mapped per million total uniquely mapped paired-end reads (FPKM). Amongst these transcripts, 14,298 (72%) were identified with more than 1 read per million in at least one library. Some transcripts were represented by many reads. Moreover, 50% of the reads mapped to only 77 transcript sequences and 90% mapped to 2,878 transcripts. The top twenty of these transcripts are shown in Table 2. Amongst these transcripts, several are associated with energy metabolism (cytochrome c oxidase subunit I, II and III, cytochrome b, ATP synthase subunit alpha, NADH dehydrogenase subunit I and NADHubiquinone oxidoreductase chain 3) or locomotion (alpha skeletal muscle actin, troponin T, myosin regulatory light chain 2, tropomyosin beta chain, myoglobin, myotilin, myosin 1 and myosin 7). These results were consistent with the physiological role of genes expected in the surveyed tissue.
To assess the consistency of gene expression profile measurements, the pairwise individual-to-individual Pearson correlation coefficient of the gene expression levels was calculated. The correlations were very high between individuals (r > 0.92) (Additional file 1: Table  S1). The shared and unique presence of transcripts is shown in Figure 1. 17,172 (87%) of the transcripts were shared among the three samples. However, approximately 2% of the transcripts are only expressed in one sample.

SNP discovery and annotation
For SNP calling, BWA was used to map the pairedreads from each sample to the bovine reference genome sequence. The SAM tools package was used for SNP discovery using stringent parameters (e.g. minimum coverage of 8 reads and mapping quality of 20). SAMtools can identify single base substitutions as well as small insertions and deletions; however, only SNPs were considered in the current analysis. In total 34,376 different SNP positions were detected with the RNA-Seq reads. Amongst these SNPs, 8,974 (26%) were homozygous in all three sequenced samples, corresponding presumably to differences between Limousin and the Hereford bovine whole-genome reference sequence [13]. A comparable number of SNPs were discovered by Canovas et al. (2010) using a similar total number of RNA-Seq reads (~118 millions reads). They identified~100,000 SNPs located in genes expressed in milk samples from Holstein cows. However, only 33,045 SNPs (32%) were polymorphic within their seven Holstein cows [33].
In our study, we found that there were 30,998 biallelic SNPs mapping to coding regions, 38.6% of which were previously found and recorded in dbSNP. This high percentage of novel SNPs, even though there are currently more than 9 millions SNPs in the public SNP database dbSNP (version 133), suggests that a large fraction of the genetic variability present in Limousin cattle still remains to be discovered. The proportion of transition substitutions were A/G, 36%, and C/T, 37%, compared to transversions A/C, 7%, G/T, 7%, A/T, 4% and C/G, 9%. This corresponds to a transition:transversion ratio of 2.65:1. The observed transition:transversion ratio is closed to the expected ratio (2:1) if all substitutions were equally likely.
Amongst these bi-allelic SNPs, 17,011 (55%) were found using Ensembl's Variant Effect Predictor in a predicted coding region. 3,791 (22.23%) resulted in an amino acid change (nonsynonymous coding SNP; nscSNP) found in 2,438 different genes. The percentage of nonsynonymous changes in the coding region found in our study was lower compared to whole-genome [24][25][26][27] [25]. They found 11,713 nscSNPs in 4,643 different genes. However, our results were similar to the rate found in another  transcriptome-based study [34]. Huang and collaborators (2012) found 1,779 nscSNPs (in 1,369 genes) out of 6,941 coding SNPs (~25%) identified by sequencing the transcriptomes of leukocytes from three animals from three different breeds [34]. The broader gene coverage when sequencing DNA versus RNA might contribute to the discrepancy in the rate of nscSNPs found between whole-genome and transcriptome-based studies. The deleterious effect of non-synonymous SNPs were analysed using the SIFT and PolyPhen algorithms. In order to use these programmes, sequences flanking the bovine nscSNPs were mapped onto the human genome and custom scripts were used to extract the human position orthologous to each bovine SNP position. We selected only bovines nscSNPs for which the two bases before and the two bases after the SNP exactly matched the human sequence. The human chromosomal position and the bovine alleles were combined to produce "pseudo human" variant positions and then used to query SIFT and PolyPhen. Using this conservative approach, we could retrieve the human "orthologous" position for 206 different bovine nscSNPs.
Using SIFT, we found that 90 different "pseudo human" coding variants were damaging. The three Limousin animals used were homozygous or heterozygous for 41 and 68 of these damaging SNPs, respectively. The difference between the number of SNPs found homozygous and heretozygous, reflects the fact that deleterious alleles are less likely to be homozygous. All three Limousin animals were homozygous for 17 damaging ncSNPs, including 13 SNPs with a genotype probability score above 20 (in all 3 samples) and 8 SNPs with a genotype probability score of 99 (in at least one sample).
Using PolyPhen-2, we found 69 different damaging "pseudo human" coding variants. 29 SNPs were homozygous and 52 SNPs heterozygous in at least one of the three Limousin samples. All Limousin animals were homozygous for 12 damaging nscSNPs, including 10 SNPs with a genotype probability score above 20 (in all 3 samples) and 6 SNPs with a genotype probability score of 99 (in at least one sample).
Fifty damaging nscSNPs were found by both SIFT and PolyPhen-2 algorithms, including 5 high-confidence nscSNPs for which all three Limousin animals are homozygous (Additional file 2: Table S2).
Gene Ontology analysis was performed with all genes containing nscSNPs. Out of the 2,438 genes, 1,092 (45%) were assigned to one or more GO annotations. In total 3,589, 2,892 and 8,172 GO terms were obtained for biological processes, cellular components and molecular functions, respectively. GO term analysis showed a significant enrichment of specific GO terms when comparing the annotations of SNP-containing genes against all unique transcripts from the bovine reference transcriptome. A summary of the classification of these genes into major biological process, cell component and molecular function categories is presented in Additional file 3: Table S3. Genes encoding proteins from the cytoskeleton and the extracellular matrix, or involved in cell cycle and cellular response are significantly over-represented. This finding might be explained by the high level of expression of these genes, that likely translates into greater sequence coverage and ultimately in a larger proportion of SNPs being identified in specific functional groups of genes. No significant enrichment in KEGG terms/pathways was found.
The positions of the 34,376 different SNPs predicted with the RNA-Seq reads were compared to the position on the UMD3.1 bovine genome assembly of know quantitative trait loci (QTLs) deposited in the public database AnimalQTLdb [39]. 32,631 SNPs were located in 3,855 different QTL regions (Additional file 4: Table S4). For example, 2,116 different SNPs are found in 16 QTL regions for meat tenderness score; whereas 14,560 SNPs are within 121 QTL regions for marbling score. QTLs were sorted into two groups (meat quality/muscle-related QTLs versus other QTLs) and the number of SNPs found in these two groups were counted. We then performed a Chi-squared test and found a significant difference (P = 0) in the number of SNPs between the two groups (Additional file 5: Table S5), suggesting an enrichment of SNPs in meat/muscle related QTLs. The high number of predicted SNPs located within known QTL regions, particularly in chromosomal regions harbouring QTLs for meat quality-related traits, indicates that the collection of SNPs found in the Longissimus thoraci transcriptome should allow the detection of candidate quantitative trait nucleotides responsible for the genetic variability of some of these traits.

Selection of candidate SNPs and validation
To analyse the accuracy of RNA-Seq technology for SNP detection, a set of SNPs were selected for validation by genotyping. Non-synonymous SNPs are of particular interest because they are more likely to alter the structure and biological function of a protein, and therefore could be the causative mutations underlying important phenotypes. We therefore selected nscSNPs for validation. All suitable putative bi-allelic nscSNPs were evaluated with the Illumina ADT software. 2,452 nscSNPs (65%) with ADT score >0.6 passed the filtering step. In order to increase the probability of an in silico detected SNP being a truly polymorphic site, we selected nscSNPs already found in dbSNP. Finally, 48 putative nscSNPs detected in 38 genes were selected ( Table 3).
The 48 selected SNPs were genotyped on the three original Limousin bull calves used for the RNA-Seq work, using llumina's GoldenGate BeadXpress system. From the 48 SNPs that were genotyped, 11 SNP assays failed to work (23%), equivalent to a conversion rate of~77%. We had 100% call rate for all remaining 37 SNPs with these three DNA samples (Table 4). A similarly low assay conversion rate was obtained in a recent SNP genotyping project using Illumina's GoldenGate BeadXpress system and was due to failure in the synthesis of some of the oligonucleotides (unpublished data).
A comparison between genotypes obtained by direct genotyping and predicted from the RNA-Seq data show 23 discrepancies (20%) ( Table 4). A quick survey shows that discordant genotyping calls occur when genotypes have been predicted from the RNA-Seq data with a low probability (score below 20). Only two discrepancies (1.8%) remained when RNA-Seq-based genotypes having at least a probability score of 20 were selected, and no discrepancies were observed when using the highest probability threshold (score of 99). It is important to point out that the RNA-Seq-based genotypes were derived from cDNA sequences whereas the genotypes produced by genotyping were obtained from DNA samples. The two discrepancies seen after filtering with a probability score above 20 (SNP26 AG versus AA and SNP31 GG versus AG; RNA-Seq-based genotype versus BeadXPress-based genotype) could therefore possibly be true differences between RNA and corresponding DNA samples, due to A-to-I (G) RNA editing (e.g. [40] and allele-specific expression [41], respectively. The SNP discovery analysis was performed initially without filtering the individual genotypes derived from the RNA-Seq data. Following on our validation study, we further filtered the identified SNPs, using this time the highest genotype probability score. We selected SNPs for which at least one individual had a heterozygous or the alternative homozygous genotype, with a probability score equal to 99. We detected 8,407 different high-confidence SNPs among 3,867 transcripts. Amongst these SNPs, 1,966 (23%) were homozygous in all three sequenced samples; 8,199 (97%) were biallelic SNPs; 3,123 (37%) were previously found in dbSNP; 6,158 (73%) were found in coding regions and 1,242 (18%) resulted in an amino acid change (in 948 different genes). A list of the high-confidence SNPs is available, as an additional file to this manuscript (Additional file 6: Table S6).

Population genetics screens
To test the usefulness of the identified SNPs, the 48 selected nscSNPs were assessed by genotyping a total of 90 bovine samples (including the three Limousin samples used for the RNA-Seq work) representing the 9 major breeds used in France, an African taurine breed (Watusi), and two other Bovinae species (European bison and Greater Koudou).
As reported above, 8 SNP assays failed to work in all samples. SNP call rate ranged from 55% (rs42555633) to 100%, whereas the call rate for bovine DNA samples ranged from 93% to 98%.
The majority (95%) of the selected SNPs with working assays, generated data with the European bison and the Greater Koudou samples (35/37 and 27/37 SNPs, respectively) ( Table 5). This could be expected since the markers were developed from (conserved) intra-genic regions. Only 3 SNPs exhibited polymorphisms in these two outcross species (2 SNPs in European bison and 2 SNPs in Greater Koudou). However, due to the small sample size (n = 1), this number is likely to be downwardly biased and a higher proportion of SNPs may in fact be polymorphic and therefore prove useful in these species. As expected from the phylogenetics of these species, the proportions of working SNPs were lower in the Greater Koudou than in the European bison.
The observed allele frequencies for the all autosomal SNPs with a SNP call rate above 92% are shown in Table 5, for each cattle population. All autosomal SNPs had a minor allele frequency (MAF) >= 0.04 in all populations, with the exception of 13 SNPs which had a fixed allele in at least one population. The highest SNP MAF observed was 0.50. The mean MAF for all autosomal markers ranged from 0.19 (HOL) to 0.27 (LIM).
The observed heterozygosities, expected heterozygosities under HWE for the observed population allele frequencies, and significance level for the test for departures from HWE for each autosomal SNP, are shown in Additional file 7: Table S7. All these markers were in agreement with HWE (P = 0.001). The mean observed heterozygosity estimated for all autosomal markers, for each population ranged from 0.259 (+/− 0.176) to 0.386 (+/− 0.230). The mean observed heterozygosities in our populations were similar to values estimated in previous studies, including a study that used a whole-genome SNP    panel to characterise the genetic diversity of several French cattle breeds [42,43]. The overall genetic differentiation among breeds was moderate (F ST = 10.9% and G ST = 9.86%) but highly significant from zero (unpublished data). This genetic differentiation among breeds implies that approximately 90% of the total genetic variation was explained by individual variability. A similar genetic differentiation was previously reported in a study carried out on French breeds, using microsatellite markers [44].
The exact test for population differentiation based on allele frequency variations shows that all breeds tested were significantly different from each other (P < 0.0001, unpublished data). Genetic distances between breeds were measured by pair-wise F ST as shown in Additional file 8: Table S8. The HOL breed was the most differentiated one. The largest similarity was detected between BLA and SAL animals (F ST = 0.0011). These results were in agreement with a previous study that analysed the genetic relationships between BLA, HOL, LIM and SAL populations [44]. Gautier and collaborators found in their study that HOL is also the most differentiated breed; however they found that AUB and LIM animals shared the smallest F ST (F ST = 0.0353) [42]. This discrepancy with our findings might mostly be due to the LIM population they surveyed. Since their study included US LIM animals, it is possible that these LIM animals were not pure-bred animals, unlike the LIM animals we used.
The degree of genetic differentiation among the breeds studied and the high levels of significance for the between-population F ST estimates indicate a relatively low gene flow between these breeds.
Principal component analysis was performed including all animals and all autosomal loci using allele frequencies to summarise breed relationships. The analysis indicates a clear separation between the nine populations ( Figure 2), but also some variability within each breed (Additional file 9: Figure S1). A total of approximately 69% of the variance accounted for the first three dimensions of the PCA.

Functional candidate mutations
The discovered coding SNPs, especially the 8,407 highconfidence SNPs may have a direct functional effet and some of them may be involved in the genetic variability of meat quality traits.
Among the high-confidence non synonymous coding SNPs, we have identified a single polymorphism resulting in a premature stop codon. SNP rs135279925 (ENSBTAT00000007104:c.1093C>T) is located within the 10 th and last exon of CD46, a membrane cofactor protein. This variant leads to a three amino acid shortened protein. None of the sampled animals were homozygous for this mutation. The corresponding bovine gene (ENSBTAG00000005397) has three known different transcripts encoding 343, 361 and 367 amino acid long proteins. The nscSNP modifies the longest bovine protein version; however, as the last three amino acids are not conserved within the bovine proteins or between species, the polymorphism is unlikely to have a functional impact.
We also found among the high-confidence nscSNPs, the previously reported F94L mutation (rs110065568: BTA2 g.6213980C>A) in the growth differentiation factor 8 (GDF8). GDF8 is a known muscle growth factor inhibitor commonly known as myostatin (MSTN). This gene has been identified as the gene responsible for the double-muscling phenotype in cattle [45][46][47]. Numerous mutations in MSTN have been described in many breeds that cause muscle hypertrophy [45][46][47][48][49][50][51], including a non synonymous amino acid substitution (F94L) in a region of the protein known to be the inhibitory domain of the MSTN propeptide [52]. Limousin cattle are not considered a double-muscled breed, however genotyping of the SNP rs110065568 has shown that the A allele is present at high frequency [48][49][50]53]. Interestingly, the three sampled animals were homozygous for this mutation. Several studies have shown that the F94L mutation is associated with increased muscle mass, carcass yield, meat tenderness and with a reduction of collagen content in Limousin and Limousin-cross cattle [54][55][56]. The high frequency of the mutant allele in Limousin most likely reflects the effects of selection for increased muscle mass.
We found among the high-confidence polymorphisms a nscSNP in another bovine gene known to be involved in meat quality traits: the mutation A127S (rs109995479: BTA2 g.107515456C>A) in the protein kinase adenosine monophosphate-activated α3-subunit (PRKAG3). Studies have shown that mutations in the porcine PRKAG3 affect the glycogen content in muscle, and consequently, ultimate pH, meat colour, water-holding capacity, drip loss, tenderness and cooking loss [57,58]. Because of the association of this gene with meat quality traits, polymorphism screens in the bovine PRKAG3 have also been performed and several non synonymous SNPs have been identified, including SNP rs109995479 [59][60][61]. Associations between another polymorphism within PRKAG3 and meat colour traits and cooking loss have been found in cattle [62]. It will be therefore interesting to test the effects of SNP rs109995479. This nscSNP is located within a region of the gene highly conserved in mammals; however, it is not located within any of the cystathione βsynthetase domains, where the two mutations with the highest phenotypic effects (I199V and R200Q) have been found, in pig.
In addition, we identified several polymorphisms in new candidate genes for several meat quality-related traits. For example, we found a high-confidence non synonymous coding SNP (rs109813896: BTA1 g.134130474G>C) in the gene encoding the mitochondrial propionyl-coA carbolylase beta subunit (PCCB), which is involved in the catabolism of propanoate, an important intermediate in the metabolism of several amino acids. Yang and collaborators [63] have shown that a polymorphism in PCCB is associated with fat weight, in pig. Interestingly, the bovine PCCB gene lies within a QTL region for fat thickness at the 12 th rib [64]. PCCB could therefore be a good candidate gene for this trait.
We also found seven high-confidence nscSNPs (including previously discovered SNPs: rs136458240, rs211315064 and rs209586352) in the gene encoding the heparin sulfate proteoglycan 2 (HSPG2, ENSBTAG00000017122). This gene encodes a large proteoglycan that is a component of the extracellular matrix. Choi and collaborators [65] found an association between a polymorphism within this gene and marbling score, in pig. The bovine HSPG2 gene is located within a marbling score QTL [66] and could therefore be a good candidate for this phenotype.

Conclusions
Our results represent the first study of gene-based SNPs discovered using RNA-seq in bovine muscle. Our results show that RNA-Seq is a fast and efficient method to identify SNPs in coding regions and we identified more than 34,000 putative SNPs (including more than 8,000 highconfidence SNPs). More than 60% of these SNPs are completely novel. The high percentage of validation confirms the utility of the SNP-mining process and the stringent quality criteria for distinguishing sequence variations from sequencing errors or artifacts introduced during the preparation of the cDNA libraries. The RNA-Seq data and the collection of newly discovered coding SNPs improve the genomic resources available for cattle, especially for beef breeds. The large amount of variation present in genes expressed in Limousin Longissimus thoracis, especially the large number of non synonymous coding SNPs, may prove useful to study the mechanisms underlying the genetic variability of meat quality traits. The coding SNPs could also be used to study allele-specific gene expression.
Our approach could be further improved in order to reduce the cost of SNP discovery and validation. Higher multiplexing of cDNA libraries prior to sequencing, would reduce sequencing cost while still allowing SNP discovery and genotype assignment. With continued improvements in next-generation DNA sequencing technologies, throughput will increase while sequencing costs are expected to decrease. When relevant tissue samples are available, it will soon be reasonable to directly perform association studies using a genotyping RNA-Seq-based approach.

Animal ethics
All animal experimentation complied with the French Veterinary Authorities' rules. No ethics approval was required by a specific committee, as the selected animals were not animals bred for experimental reasons.
Animals and tissue samples The study was conducted with three Limousin bull calves from a large study on the genetic determinism of beef and meat quality traits [67]. The three bull calves were not closely related to one another (for at least 4 generations) were fattened in a single feedlot and fed ad libidum with wet corn silage. They were humanely slaughtered in an accredited commercial slaughterhouse when they reached 16 months. Longissimus thoracis (LT) muscle samples were dissected immediately after death and tissue samples were snap frozen in liquid nitrogen and stored at −80°C until analysis.

RNA isolation and sequencing
After transfer to ice-cold RNeasy RLT lysis buffer (Qiagen, Courtaboeuf, France), LT tissue samples were homogenised using a Precellys tissue homogeniser (Bertin Technologie, Montigny-le-Bretonneux, France). Total RNA was isolated using RNeasy Midi columns (Qiagen) and then treated with RNAse-free DNase I (Qiagen) for 15 min at room temperature according to the manufacturer's protocols. The concentration of total RNA was measured with a Nanodrop ND-100 instrument (Thermo Scientific, Ilkirch, France) and the quality was assessed with an RNA 6000 Nano Labchip kit using an Agilent 2100 Bioanalyzer (Agilent Technologies, Massy, France). All three samples had an RNA Integrity Number (RIN) value greater than eight.
The mRNA-Seq libraries were prepared using the TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA) according to the manufacturer's instructions. Briefly, Poly-A containing mRNA molecules were purified from 4 μg total RNA of each sample using oligo(dT) magnetic beads and fragmented into 150-400 bp pieces using divalent cations at 94°C for 8 min. The cleaved mRNA fragments were converted to double-stranded cDNA using SuperScript II reverse transcriptase (Life Technologies, Saint Aubin, France) and primed by random primers. The resulting cDNA was purified using Agencourt AMPure W XP beads (Beckman Coulter, Villepinte, France). Then, cDNA was subjected to endrepair and phosphorylation and subsequent purification was performed using Agencourt AMPure W XP beads (Beckman Coulter). These repaired cDNA fragments were 3′-adenylated producing cDNA fragments with a single ' A' base overhung at their 3′-ends for subsequent adapter-ligation. Illumina adapters containing indexing tags were ligated to the ends of these 3′-adenylated cDNA fragments followed by two purification steps using Agencourt AMPure W XP beads (Beckman Coulter). Ten rounds of PCR amplification were performed to enrich the adapter-modified cDNA library using primers complementary to the ends of the adapters. The PCR products were purified using Agencourt AMPure W XP beads (Beckman Coulter) and size-selected (200 ± 25 bp) on a 2% agarose Invitrogen E-Gel (Thermo Scientific). Libraries were then checked on an Agilent Technologies 2100 Bioanalyzer using the Agilent High Sensitivity DNA Kit and quantified by quantitative PCR with the QPCR NGS Library Quantification kit (Agilent Technologies). After quantification, tagged cDNA libraries were pooled in equal ratios and a final qPCR check was performed post-pooling. The pooled libraries were used for 2×100 bp pairedend sequencing on one lane of the Illumina HiSeq2000 with a TruSeq SBS v3-HS Kit (Illumina). After sequencing, the samples were demultiplexed and the indexed adapter sequences were trimmed using the CASAVA v1.8.2 software (Illumina).

Mapping reads to reference transcriptome and gene expression counts
The Bos taurus reference transcriptome was downloaded from Ensembl (version 63, Bos_taurus.Btau_4.0.63.cdna. all.fa). To align the reads back to the assembled reference transcriptome the BWA programme (version 0.5.9-r16) was used [68]. Reads were mapped for each sample separately to the assembled transcriptome. The BWA default values were used for mapping. Properly paired reads with a mapping quality of at least 30 (−q = 30) were extracted from the resulting BAM file using SAMtools [69] for further analyses. Properly paired is defined as both left and right reads mapped in opposite directions on the same transcript at a distance compatible with the expected mean size of the fragments (<500-bp). Custom scripts were developed to identify paired-reads mapping to single locations and with the expected distance. Read pairs mapping to separate chromosomes were discarded for the present study. Transcriptome contamination was assessed by mapping with BWA reads on a sequence library, containing E. coli, phiX and yeast genome sequences. The number of pairedreads uniquely aligning to transcribed regions of each transcript was calculated for all genes in the annotated transcriptome. The transcript paired-read count was calculated as the number of unique paired-reads that aligned within the exons of each transcript, based on the coordinates of mapped reads. The expression level of each gene was calculated in FPKM (fragments per kilobase per million sequenced reads) using a custom script based on Tapnel et al. (2010) [70].

Polymorphism identification
BWA was also used to map reads onto the bovine genome reference sequence (version UMD3.1, [71]. Only reliable properly paired BWA mapped reads were considered for Single Nucleotide Polymorphism (SNP) calling. Indels were not considered because alternative splicing impedes reliable indel discovery. SNPs were called using the SAMtools software package. Genotype likelihoods were computed using the SAMtools utilities and variable positions in the aligned reads compared to the reference were called with the BCFtools utilities [72]. SNPs were called only for positions with a minimal mapping quality (−Q) of 30, a minimum coverage (−d) of 4 and a maximum read depth (−D) of 10,000,000.

Functional annotation of detected SNPs
The functional effect of the newly discovered SNPs on known transcripts were analysed using Ensembl's Variant Effect Predictor v2.5, following local installation [73].
The deleterious effect of non-synonymous SNPs were analysed using the SIFT (Sorting Intolerant From Tolerant; http://sift.bii.a-star.edu.sg; [74] and PolyPhen-2 (Polymorphism Phenotyping 2; http://genetics.bwh.harvard.edu/pph2/; [75] programmes. In order to use these two programmes, sequences flanking the bovine nscSNPs were mapped onto the human genome (version GRCh37/ hg19) using MegaBLAST [76] and custom scripts were used to extract the human position orthologous to each bovine SNP position. The human chromosomal position and the bovine alleles were then used to query SIFT and PolyPhen. Default settings were used for both programmes. We refered to damaging SNPs, SNPs that were identified as damaging and not tolerated, using PolyPhen-2 and SIFT, respectively.
In order to evaluate whether SNP-containing genes were significantly enriched for specific gene ontology (GO) terms and KEGG pathways compared to all annotated bovine genes, gene enrichment analyses were conducted using the FATIGO tool of the online software suite Babelomics (http://babelomics.bioinfo.cipf.es; [77]. Genes were assigned their Ensembl identities as input for Babelomics. Only one copy of each gene was used. Default parameter settings were used for the analysis. Statistical assessment of annotation differences between the two sets of sequences (SNP-containing genes versus all the other bovine genes) was carried out for each FATIGO analysis, using the Fisher Exact Test with correction for multiple testing.
Selection of candidate SNPs for genotyping assay After SNP detection, in silico evaluation of candidate SNPs was carried out to select a panel of candidate SNPs for validation. SNP selection was based on the results from the Illumina Assay Design Tool. The SNP score from the Illumina Assay Design Tool (referred to as the Assay Design Score/ADS) utilises factors including template GC content, melting temperature, sequence uniqueness, and self-complementarity to filter the candidate SNPs prior to further inspection. The Assay Design Score (assigned between 0 and 1) is indicative of the ability to design suitable oligos within the 60 bp up/downstream flanking regions, and the expected success of the assay when genotyped with the Illumina GoldenGate chemistry. Following the Illumina guidelines, all SNPs with a score below 0.4 should be discarded; SNPs with a score above 0.4 accepted, with SNPs scoring above 0.6 being used preferentially. SNP flanking sequences were retrieved and only SNP sequences with unambiguous 121 bases (60 bases up/down-stream of each SNP position) were submitted to Illumina to assess the design quality. SNPs with ADS showing a quality score above of 0.6 were retained for analysis.

SNP validation by high-throughput genotyping
Ninety bovine DNA samples were genotyped for each selected SNP using Illumina's GoldenGate assay. These samples include 11 Aubrac (AUB), 11 Blonde d' Aquitaine (BLA), 11 Charolais (CHA), 11 Holstein (HOL), 11 Limousin (LIM), 11 Montbéliard (MON), 11 Salers (SAL), 6 Maine-Anjou (MAN), 6 Normande (NOR) and 1 Watusi (WAT) animals. These animals were not closely related to one another (for at least 4 generations) according to genealogical records from the French Centre de Traitement de l'Information Génétique (INRA, Jouy-en-Josas, France). To assess the utility of developed markers in related species, two Bovinae species; the European bison (BIS, Bison bonasus) and a more distantly related species; the Greater Koudou (KOU, Tragelaphus strepsiceros) were also genotyped. Blood samples were collected at the Parc du Rénou Zoo (Le Vigen, France). Genomic DNA was extracted from whole-blood or semen samples using the Qiasymphony SP robotic system and DNA Midi kit (Qiagen). Quality of DNA was checked using a Nanodrop ND-100 spectrophotometer (Thermo Scientific) and quantity was estimated with Quant-iT Picogreen dsDNA kit (Life Technologies) on an ABI 7900HT (Life Technologies). All DNA samples were standardised to 50 ng/μL. All animal manipulations were done according to good animal practice as defined by the French Veterinary Authorities.
High-throughput genotyping reactions were performed using Illumina's GoldenGate BeadXpress system, according to the manufacturer's protocol. Oligonucleotides were designed, synthesised, and assembled into a custom oligo pooled assay (OPA) by Illumina. Automatic allele calling for each SNP was accomplished with the GenomeStudio software (Illumina). All genotypes were manually checked and re-scored if any errors in calling homozygous or heterozygous clusters were evident. Genotype calls were exported in spreadsheets from the GenomeStudio data analysis software for further analysis.

Population genetics analyses
Genetic diversity parameters within each population were calculated using the GENETIX 4.05.2 software package [78]. Tests for deviation from Hardy-Weinberg equilibrium were performed by the GENEPOP 3.4 software [79], using the exact test of Guo and Thompson (1992) [80]. Genetic differentiation among and within the populations was estimated based on F-statistics (F ST ) according to Weir and Cockerham (1984) [81]) using the GENEPOP and GENETIX software packages. Test for population differentiation was performed as implemented in GENEPOP. The Reynolds genetic distance (D R ) was calculated for each pair of populations based on allele frequencies [82] using the GENETIX software. Principal component analysis (PCA) was performed using the GENETIX programme from allele doses for each individual.

Data availability
The sequencing data have been submitted to the European Nucleotide Archive (accession number ERP002220).