Skip to main content

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

Abstract

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 [310]. 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 high-quality 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 fine-mapping of QTL regions has now been expedited. In addition, the advent of novel sequencing technologies [1923] offers now a new opportunity for the identification of QTNs, with the ability to partially or completely re-sequence mammalian genomes, in a relatively cost-effective 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 [2430]. For example, Eck et al. (2009) generated the first single cattle genome sequence by a next-generation sequencing method [24]. By sequencing the whole-genome 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.

Results and discussion

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 high-throughput 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 [3338]. 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].

Table 1 Summary of reads mapping to the bovine transcriptome

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 NADH-ubiquinone 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.

Table 2 Top twenty transcripts with most assigned reads

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.

Figure 1
figure 1

Unique and shared transcripts within the three muscle samples (Venn diagram).

SNP discovery and annotation

For SNP calling, BWA was used to map the paired-reads 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 bi-allelic 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 [2427] studies performed previously in cattle. For example, Kawahara-Miki et al. (2011) have reported up to 57.3% of nscSNPs in coding regions in the whole-genome of a single individual of the Japanese Kuchinoshima-Ushi native cattle breed [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 extra-cellular 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).

Table 3 List of selected SNPs

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).

Table 4 Genotype comparison

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 bi-allelic 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.

Table 5 Details and allele frequencies of SNPs in the nine French cattle breeds, and genotypes in the three other samples

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 (FST = 10.9% and GST = 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 FST 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 (FST = 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 FST (FST = 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 FST 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.

Figure 2
figure 2

Principal component analysis. Per cent value in each axis indicates contribution to the total genetic variation.

Functional candidate mutations

The discovered coding SNPs, especially the 8,407 high-confidence 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 10th 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 [4547]. Numerous mutations in MSTN have been described in many breeds that cause muscle hypertrophy [4551], 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 [4850, 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 [5456]. 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 [5961]. 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 12th 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 high-confidence 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.

Methods

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® XP beads (Beckman Coulter, Villepinte, France). Then, cDNA was subjected to end-repair and phosphorylation and subsequent purification was performed using Agencourt AMPure® 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® 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® 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 paired-end 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 paired-reads 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/down-stream 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 (FST) 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 (DR) 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).

References

  1. Edwards CJ, Bradley DG, MacHugh DE, Dobney K, Martin L, Russell N: Ancient DNA analysis of 101 cattle remains: limits and prospects. J Archaeol Sci. 2004, 31: 695-710. 10.1016/j.jas.2003.11.001.

    Article  Google Scholar 

  2. Davis GP, DeNise SK: The impact of genetic markers on selection. J Anim Sci. 1998, 76: 2331-2339.

    CAS  PubMed  Google Scholar 

  3. Fujii J, Otsu K, Zorzato F, de Leon S, Khanna VK, Weiler JE: Identification of a mutation in porcine ryanodine receptor associated with malignant hyperthermia. Science. 1991, 253: 448-451. 10.1126/science.1862346.

    Article  CAS  PubMed  Google Scholar 

  4. Milan D, Jeon JT, Looft C, Amarger V, Robic A, Thelander M: A mutation in PRKAG3 associated with excess glycogen content in pig skeletal muscle. Science. 2000, 288: 1248-1251. 10.1126/science.288.5469.1248.

    Article  CAS  PubMed  Google Scholar 

  5. Grisart B, Coppieters W, Farnir F, Karim L, Ford C, Berzi P: Positional candidate cloning of a QTL in dairy cattle: identification of a missense mutation in the bovine DGAT1 gene with major effect on milk yield and composition. Genome Res. 2002, 12: 222-231. 10.1101/gr.224202.

    Article  CAS  PubMed  Google Scholar 

  6. Blott S, Kim JJ, Moisio S, Schmidt-Küntzel A, Cornet A, Berzi P: Molecular dissection of a quantitative trait locus: a phenylalanine-to-tyrosine substitution in the transmembrane domain of the bovine growth hormone receptor is associated with a major effect on milk yield and composition. Genetics. 2003, 163: 253-266.

    PubMed Central  CAS  PubMed  Google Scholar 

  7. Van Laere AS, Nguyen M, Braunschweig M, Nezer C, Collette C, Moreau L: A regulatory mutation in IGF2 causes a major QTL effect on muscle growth in the pig. Nature. 2003, 425: 832-836. 10.1038/nature02064.

    Article  CAS  PubMed  Google Scholar 

  8. Cohen-Zinder M, Seroussi E, Larkin DM, Loor JJ, Everts-van der Wind A, Lee JH: Identification of a missense mutation in the bovine ABCG2 gene with a major effect on the QTL on chromosome 6 affecting milk yield and composition in Holstein cattle. Genome Res. 2005, 15: 936-944. 10.1101/gr.3806705.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  9. Murphy SK, Nolan CM, Huang Z, Kucera KS, Freking BA, Smith TP: Callipyge mutation affects gene expression in cis: a potential role for chromatin structure. Genome Res. 2006, 16: 340-346. 10.1101/gr.4389306.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  10. Clop A, Marcq F, Takeda H, Pirottin D, Tordoir X, Bibé B: A mutation creating a potential illegitimate microRNA target site in the myostatin gene affects muscularity in sheep. Nat Genet. 2006, 38: 813-818. 10.1038/ng1810.

    Article  CAS  PubMed  Google Scholar 

  11. Dekkers JC: Commercial application of marker- and gene-assisted selection in livestock: strategies and lessons. J Anim Sci. 2004, 82: E313-328.

    PubMed  Google Scholar 

  12. Andersson L, Georges M: Domestic-animal genomics: deciphering the genetics of complex traits. Nat Rev Genet. 2004, 5: 202-212. 10.1038/nrg1294.

    Article  CAS  PubMed  Google Scholar 

  13. Bovine Genome Sequencing and Analysis Consortium: The genome sequence of taurine cattle: a window to ruminant biology and evolution. Science. 2009, 324: 522-528.

    Article  Google Scholar 

  14. Bovine HapMap Consortium: Genome-wide survey of SNP variation uncovers the genetic structure of cattle breeds. Science. 2009, 324: 528-532.

    Article  Google Scholar 

  15. Van Tassell CP, Smith TP, Matukumalli LK, Taylor JF, Schnabel RD, Lawley CT: SNP discovery and allele frequency estimation by deep sequencing of reduced representation libraries. Nat Methods. 2008, 5: 247-252. 10.1038/nmeth.1185.

    Article  CAS  PubMed  Google Scholar 

  16. Matukumalli LK, Lawley CT, Schnabel RD, Taylor JF, Allan MF, Heaton MP: Development and characterization of a high density SNP genotyping assay for cattle. PLoS One. 2009, 4: e5350-10.1371/journal.pone.0005350.

    Article  PubMed Central  PubMed  Google Scholar 

  17. Illumina’s BovineHD Genotyping BeadChip. http://www.illumina.com/documents/products/datasheet/datasheet_bovineHD.pdf,

  18. Meuwissen TH, Hayes BJ, Goddard ME: Prediction of total genetic value using genome-wide dense marker maps. Genetics. 2001, 157: 1819-1829.

    PubMed Central  CAS  PubMed  Google Scholar 

  19. Margulies M, Egholm M, Altman WE, Attiya S, Bader JS, Bemben LA: Genome sequencing in microfabricated high-density picolitre reactors. Nature. 2005, 437: 376-380.

    PubMed Central  CAS  PubMed  Google Scholar 

  20. Bentley DR, Balasubramanian S, Swerdlow HP, Smith GP, Milton J, Brown CG: Accurate whole human genome sequencing using reversible terminator chemistry. Nature. 2008, 456: 53-59. 10.1038/nature07517.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  21. McKernan KJ, Peckham HE, Costa GL, McLaughlin SF, Fu Y, Tsung EF: Sequence and structural variation in a human genome uncovered by short-read, massively parallel ligation sequencing using two-base encoding. Genome Res. 2009, 19: 1527-1541. 10.1101/gr.091868.109.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  22. Harris TD, Buzby PR, Babcock H, Beer E, Bowers J, Braslavsky I: Single-molecule DNA sequencing of a viral genome. Science. 2008, 320: 106-109. 10.1126/science.1150427.

    Article  CAS  PubMed  Google Scholar 

  23. Drmanac R, Sparks AB, Callow MJ, Halpern AL, Burns NL, Kermani BG: Human genome sequencing using unchained base reads on self-assembling DNA nanoarrays. Science. 2010, 327: 78-81. 10.1126/science.1181498.

    Article  CAS  PubMed  Google Scholar 

  24. Eck SH, Benet-Pagès A, Flisikowski K, Meitinger T, Fries R, Strom TM: Whole genome sequencing of a single Bos taurus animal for single nucleotide polymorphism discovery. Genome Biol. 2009, 10: R82-10.1186/gb-2009-10-8-r82.

    Article  PubMed Central  PubMed  Google Scholar 

  25. Kawahara-Miki R, Tsuda K, Shiwa Y, Arai-Kichise Y, Matsumoto T, Kanesaki Y: Whole-genome resequencing shows numerous genes with nonsynonymous SNPs in the Japanese native cattle Kuchinoshima-Ushi. BMC Genomics. 2011, 12: 103-10.1186/1471-2164-12-103.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  26. Zhan B, Fadista J, Thomsen B, Hedegaard J, Panitz F, Bendixen C: Global assessment of genomic variation in cattle by genome resequencing and high-throughput genotyping. BMC Genomics. 2001, 12: 557-

    Article  Google Scholar 

  27. Stothard P, Choi JW, Basu U, Sumner-Thomson JM, Meng Y, Liao X: Whole genome resequencing of black Angus and Holstein cattle for SNP and CNV discovery. BMC Genomics. 2011, 12: 559-10.1186/1471-2164-12-559.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  28. Canavez FC, Luche DD, Stothard P, Leite KR, Sousa-Canavez JM, Plastow G: Genome sequence and assembly of Bos indicus. J Hered. 2012, 103: 342-348. 10.1093/jhered/esr153.

    Article  CAS  PubMed  Google Scholar 

  29. Bickhart DM, Hou Y, Schroeder SG, Alkan C, Cardone MF, Matukumalli LK: Copy number variation of individual cattle genomes using next-generation sequencing. Genome Res. 2012, 22: 778-790. 10.1101/gr.133967.111.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  30. Larkin DM, Daetwyler HD, Hernandez AG, Wright CL, Hetrick LA, Boucek L: Whole-genome resequencing of two elite sires for the detection of haplotypes under selection in dairy cattle. Proc Natl Acad Sci U S A. 2012, 109: 7693-7698. 10.1073/pnas.1114546109.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  31. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008, 5: 621-628. 10.1038/nmeth.1226.

    Article  CAS  PubMed  Google Scholar 

  32. Cloonan N, Forrest AR, Kolle G, Gardiner BB, Faulkner GJ, Brown MK: Stem cell transcriptome profiling via massive-scale mRNA sequencing. Nat Methods. 2008, 5: 613-619. 10.1038/nmeth.1223.

    Article  CAS  PubMed  Google Scholar 

  33. Canovas A, Rincon G, Islas-Trejo A, Wickramasinghe S, Medrano JF: SNP discovery in the bovine milk transcriptome using RNA-Seq technology. Mamm Genome. 2010, 21: 592-598. 10.1007/s00335-010-9297-z.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  34. Huang W, Nadeem A, Zhang B, Babar M, Soller M, Khatib H: Characterization and comparison of the leukocyte transcriptomes of three cattle breeds. PLoS One. 2012, 7: e30244-10.1371/journal.pone.0030244.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  35. Wickramasinghe S, Rincon G, Islas-Trejo A, Medrano JF: Transcriptional profiling of bovine milk using RNA sequencing. BMC Genomics. 2012, 13: 45-10.1186/1471-2164-13-45.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  36. Baldwin RL, Wu S, Li W, Li C, Bequette BJ, Li RW: Quantification of transcriptome responses of the rumen epithelium to butyrate infusion using RNA-seq technology. Gene Regul Syst Bio. 2012, 6: 67-80.

    PubMed Central  CAS  PubMed  Google Scholar 

  37. Li RW, Rinaldi M, Capuco AV: Characterization of the abomasal transcriptome for mechanisms of resistance to gastrointestinal nematodes in cattle. Vet Res. 2011, 42: 114-10.1186/1297-9716-42-114.

    Article  PubMed Central  PubMed  Google Scholar 

  38. Driver AM, Peñagaricano F, Huang W, Ahmad KR, Hackbart KS, Wiltbank MC: RNA-Seq analysis uncovers transcriptomic variations between morphologically similar in vivo- and in vitro-derived bovine blastocysts. BMC Genomics. 2012, 13: 118-10.1186/1471-2164-13-118.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  39. Hu ZL, Fritz ER, Reecy JM: AnimalQTLdb: a livestock QTL database tool set for positional QTL information mining and beyond. Nucleic Acids Res. 2007, 35: D604-D609. 10.1093/nar/gkl946.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  40. Peng Z, Cheng Y, Tan BC, Kang L, Tian Z, Zhu Y: Comprehensive analysis of RNA-Seq data reveals extensive RNA editing in a human transcriptome. Nat Biotechnol. 2012, 30: 253-260. 10.1038/nbt.2122.

    Article  CAS  PubMed  Google Scholar 

  41. Pastinen T: Genome-wide allele-specific analysis: insights into regulatory variation. Nat Rev Genet. 2010, 11: 533-538.

    Article  CAS  PubMed  Google Scholar 

  42. Gautier M, Laloe D, Moazami-Goudarzi K: Insights into the genetic history of French cattle from dense SNP data on 47 worlwide breeds. PLoS One. 2010, 5: e13038-10.1371/journal.pone.0013038.

    Article  PubMed Central  PubMed  Google Scholar 

  43. Blott SC, Williams JL, Haley CS: Genetic relationships among European cattle breeds. Anim Genet. 1998, 29: 273-282. 10.1046/j.1365-2052.1998.00327.x.

    Article  CAS  PubMed  Google Scholar 

  44. Amigues Y, Boitard S, Bertrand C, Sancristobal M, Rocha D: Genetic characterization of the Blonde d’Aquitaine cattle breed using microsatellite markers and relationship with three other French cattle populations. J Anim Breed Genet. 2011, 128: 201-208. 10.1111/j.1439-0388.2010.00890.x.

    Article  PubMed  Google Scholar 

  45. Grobet L, Martin LJ, Poncelet D, Pirottin D, Brouwers B, Riquet J: A deletion in the bovine myostatin gene causes the double-muscled phenotype in cattle. Nat Genet. 1997, 17: 71-74. 10.1038/ng0997-71.

    Article  CAS  PubMed  Google Scholar 

  46. Kambadur R, Sharma M, Smith TP, Bass JJ: Mutations in myostatin (GDF8) in double-muscled Belgian Blue and Piedmontese cattle. Genome Res. 1997, 7: 910-916.

    CAS  PubMed  Google Scholar 

  47. McPherron AC, Lee SJ: Double muscling in cattle due to mutations in the myostatin gene. Proc Natl Acad Sci U S A. 1997, 94: 12457-12461. 10.1073/pnas.94.23.12457.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  48. Grobet L, Poncelet D, Royo LJ, Brouwers B, Pirottin D, Michaux C: Molecular definition of an allelic series of mutations disrupting the myostatin function and causing double-muscling in cattle. Mamm Genome. 1989, 9: 210-213.

    Article  Google Scholar 

  49. Smith JA, Lewis AM, Wiener P, Williams JL: Genetic variation in the bovine myostatin gene in UK beef cattle: allele frequencies and haplotype analysis in the South Devon. Anim Genet. 2000, 31: 306-309. 10.1046/j.1365-2052.2000.00521.x.

    Article  CAS  PubMed  Google Scholar 

  50. Dunner S, Miranda ME, Amigues Y, Cañón J, Georges M, Hanset R: Haplotype diversity of the myostatin gene among beef cattle breeds. Genet Sel Evol. 2003, 35: 103-118. 10.1186/1297-9686-35-1-103.

    Article  PubMed Central  PubMed  Google Scholar 

  51. Marchitelli C, Savarese MC, Crisà A, Nardone A, Marsan PA, Valentini A: Double muscling in Marchigiana beef breed is caused by a stop codon in the third exon of myostatin gene. Mamm Genome. 2003, 14: 392-395. 10.1007/s00335-002-2176-5.

    Article  CAS  PubMed  Google Scholar 

  52. Jiang MS, Liang LF, Wang S, Ratovitski T, Holmstrom J, Barker C: Characterization and identification of the inhibitory domain of GDF-8 propeptide. Biochem Biophys Res Commun. 2004, 315: 525-531. 10.1016/j.bbrc.2004.01.085.

    Article  CAS  PubMed  Google Scholar 

  53. Vankan DM, Waine DR, Fortes MR: Real-time PCR genotyping and frequency of the myostatin F94L mutation in beef cattle breeds. Animal. 2010, 4: 530-534. 10.1017/S175173110999139X.

    Article  CAS  PubMed  Google Scholar 

  54. Sellick GS, Pitchford WS, Morris CA, Cullen NG, Crawford AM, Raadsma HW: Effect of myostatin F94L on carcass yield in cattle. Anim Genet. 2007, 38: 440-446. 10.1111/j.1365-2052.2007.01623.x.

    Article  CAS  PubMed  Google Scholar 

  55. Esmailizadeh AK, Bottema CD, Sellick GS, Verbyla AP, Morris CA, Cullen NG: Effects of the myostatin F94L substitution on beef traits. J Anim Sci. 2008, 86: 1038-1046.

    Article  CAS  PubMed  Google Scholar 

  56. Alexander LJ, Kuehn LA, Smith TP, Matukumalli LK, Mote B, Koltes JE: A Limousin specific myostatin allele affects longissimus muscle area and fatty acid profiles in a Wagyu-Limousin F2 population. J Anim Sci. 2009, 87: 1576-1581. 10.2527/jas.2008-1531.

    Article  CAS  PubMed  Google Scholar 

  57. Lines DS, Pitchford WS, Kruk ZA, Bottema CD: Limousin myostatin F94L variant affects semitendinosus tenderness. Meat Sci. 2009, 81: 126-131. 10.1016/j.meatsci.2008.07.007.

    Article  CAS  PubMed  Google Scholar 

  58. Ciobanu D, Bastiaansen J, Malek M, Helm J, Woollard J, Plastow G: Evidence for new alleles in the protein kinase adenosine monophosphate-activated gamma(3)-subunit gene associated with low glycogen content in pig skeletal muscle and improved meat quality. Genetics. 2001, 159: 1151-1162.

    PubMed Central  CAS  PubMed  Google Scholar 

  59. McKay SD, White SN, Kata SR, Loan R, Womack JE: The bovine 5′ AMPK gene family: mapping and single nucleotide polymorphism detection. Mamm Genome. 2003, 14: 853-858. 10.1007/s00335-003-2276-x.

    Article  CAS  PubMed  Google Scholar 

  60. Yu SL, Kim JE, Chung HJ, Jung KC, Lee YJ, Yoon DH: Molecular cloning and characterization of bovine PRKAG3 gene: structure, expression and single nucleotide polymorphism detection. J Anim Breed Genet. 2005, 122: 294-301. 10.1111/j.1439-0388.2005.00545.x.

    Article  CAS  PubMed  Google Scholar 

  61. Roux M, Nizou A, Forestier L, Ouali A, Levéziel H, Amarger V: Characterization of the bovine PRKAG3 gene: structure, polymorphism, and alternative transcripts. Mamm Genome. 2006, 17: 83-92. 10.1007/s00335-005-0093-0.

    Article  CAS  PubMed  Google Scholar 

  62. Reardon W, Mullen AM, Sweeney T, Hamill RM: Association of polymorphisms in candidate genes with colour, water-holding capacity, and composition traits in bovine M. longissimus and M. semimembranosus. Meat Sci. 2010, 86: 270-275. 10.1016/j.meatsci.2010.04.013.

    Article  CAS  PubMed  Google Scholar 

  63. Yang F, Wang QP, He K, Wang MH, Pan YC: Association between gene polymorphisms of propanoate metabolism pathway and meat quality as well as carcass traits in pigs. Yi Chuan. 2012, 34: 872-878. 10.3724/SP.J.1005.2012.00872.

    Article  CAS  PubMed  Google Scholar 

  64. McClure MC, Morsci NS, Schnabel RD, Kim JW, Yao P, Rolf MM: A genome scan for quantitative trait loci influencing carcass, post-natal growth and reproductive traits in commercial Angus cattle. Anim Genet. 2010, 41: 597-607. 10.1111/j.1365-2052.2010.02063.x.

    Article  CAS  PubMed  Google Scholar 

  65. Choi I, Bates RO, Raney NE, Steibel JP, Ernst CW: Evaluation of QTL for carcass merit and meat quality traits in a US commercial Duroc population. Meat Sci. 2012, 92: 132-138. 10.1016/j.meatsci.2012.04.023.

    Article  PubMed  Google Scholar 

  66. MacNeil MD, Grosz MD: Genome-wide scans for QTL affecting carcass traits in Hereford x composite double backcross populations. J Anim Sci. 2002, 80: 2316-2324.

    CAS  PubMed  Google Scholar 

  67. Allais S, Levéziel H, Payet-Duprat N, Hocquette JF, Lepetit J, Rousset S: The two mutations Q204X and nt821 of the myostatin gene affect carcass and meat quality in heterozygous young bulls of French beef breeds. J Anim Sci. 2010, 88: 446-454. 10.2527/jas.2009-2385.

    Article  CAS  PubMed  Google Scholar 

  68. Li H, Durbin R: Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics. 2009, 25: 175417-60.

    Google Scholar 

  69. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N: The Sequence alignment/map (SAM) format and SAMtools. Bioinformatics. 2009, 25: 2078-2079. 10.1093/bioinformatics/btp352.

    Article  PubMed Central  PubMed  Google Scholar 

  70. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ: Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotech. 2010, 28: 511-515. 10.1038/nbt.1621.

    Article  CAS  Google Scholar 

  71. Zimin AV, Delcher AL, Florea L, Kelley DR, Schatz MC, Puiu D: A whole-genome assembly of the domestic cow. Bos taurus. Genome Biol. 2009, 10: R42-10.1186/gb-2009-10-4-r42.

    Article  PubMed  Google Scholar 

  72. Li H: A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011, 27: 2987-2993. 10.1093/bioinformatics/btr509.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  73. McLaren W, Pritchard B, Rios D, Chen Y, Flicek P, Cunningham F: Deriving the consequences of genomic variants with the Ensembl API and SNP Effect Predictor. Bioinformatics. 2010, 26: 2069-2070. 10.1093/bioinformatics/btq330.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  74. Ng PC, Henikoff S: Predicting deleterious amino acid substitutions. Genome Res. 2001, 11: 863-874. 10.1101/gr.176601.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  75. Adzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, Bork P: A method and server for predicting damaging missense mutations. Nat Methods. 2010, 7: 248-249. 10.1038/nmeth0410-248.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  76. Morgulis A, Coulouris G, Raytselis Y, Madden TL, Agarwala R, Schäffer AA: Database indexing for production MegaBLAST searches. Bioinformatics. 2008, 24: 1757-1764. 10.1093/bioinformatics/btn322.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  77. Al-Shahrour F, Minguez P, Tárraga J, Montaner D, Alloza E, Vaquerizas JMM: BABELOMICS: a systems biology perspective in the functional annotation of genome-scale experiments. Nucleic Acids Res. 2006, 34: W472-W476. 10.1093/nar/gkl172.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  78. GENETIX v. 4.05. http://www.univ-montp2.fr/~genetix/genetix/genetix.htm,

  79. Raymond M, Rousset F: GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism. J Hered. 1995, 86: 248-249.

    Google Scholar 

  80. Guo SW, Thompson EA: Performing the exact test of Hardy-Weinberg proportion for multiple alleles. Biometrics. 1992, 48: 361-372. 10.2307/2532296.

    Article  CAS  PubMed  Google Scholar 

  81. Cockerham CC, Weir BS: Covariances of relatives stemming from a population undergoing mixed self and random mating. Biometrics. 1984, 40: 157-164. 10.2307/2530754.

    Article  CAS  PubMed  Google Scholar 

  82. Reynolds J, Weir BS, Cockerham CC: Estimation of the coancestry coefficient: basis for a short-term genetic distance. Genetics. 1983, 105: 767-779.

    PubMed Central  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

The authors would like to thank Hubert Levéziel for his help, the different cattle breeding societies that provided semen and blood samples for the animals analysed in this study and Yves Amigues and colleagues at LABOGENA for their help in DNA preparation. The RNA-Seq work was funded by the INRA Animal Genetics Department (BovRNA-Seq project). The sampling of the Limousin Longissimus thoraci biopsies was part of the Qualvigène project, funded by Agence Nationale de la Recherche (contracts ANR-05-GANI-005 and ANR-05-GANI-017-01) and APIS GENE (contract 01-2005-QualviGenA-02). The authors wish to thank the anonymous reviewers for their valuable comments and suggestions, which were helpful in improving our manuscript.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Dominique Rocha.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

AD carried out the bioinformatic analysis, under the supervison of CK. DE performed the RNA-Seq experiment. BW and MB contributed to the data analyses. DE and FM performed the SNP genotyping. CM prepared the RNA samples. DR conceived the study, analysed the data and drafted the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material

Additional file 1: Table S1: Pearson correlation coefficient between individuals. (DOCX 11 KB)

Additional file 2: Table S2: List of damaging SNPs predicted by SIFT and PolypPhen-2. (XLSX 11 KB)

Additional file 3: Table S3: Enrichment of SNP-containing contigs in GO terms. (XLSX 18 KB)

Additional file 4: Table S4: List of putative SNPs located within known QTL regions. (XLSX 15 MB)

Additional file 5: Table S5: Chi-squared test details. (XLSX 10 MB)

Additional file 6: Table S6: List of the high-confidence SNPs and annotation. (XLSX 568 KB)

Additional file 7: Table S7: Details on the observed and expected heterozygosities. (DOCX 26 KB)

12864_2012_5274_MOESM8_ESM.docx

Additional file 8: Table S8: Genetic differentiation (FST) between pairs of cattle populations (above the diagonal) and Reynold’s genetic distance (DR) between pairs of cattle populations (below diagonal) as observed in this study. (DOCX 12 KB)

12864_2012_5274_MOESM9_ESM.docx

Additional file 9: Figure S1: Principal Component Analysis. Per cent value in each axis indicates contribution to the total genetic variation. (DOCX 57 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Djari, A., Esquerré, D., Weiss, B. et al. Gene-based single nucleotide polymorphism discovery in bovine muscle using next-generation transcriptomic sequencing. BMC Genomics 14, 307 (2013). https://doi.org/10.1186/1471-2164-14-307

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1471-2164-14-307

Keywords