BACKGROUND
over the past decade, genome-wide transcriptional studies discovered that a large fraction of the eukaryotic genomes is transcribed in a heterogeneous population of noncoding rnas .
these are transcripts that are not translated into a protein but act as functional rnas.
ncrnas shorter than  <dig> nucleotides are usually identified as small/short ncrna and include piwi-interacting rnas , endogeneous small interfering rnas  and micrornas  but also classical ncrnas, such as ribosomal rnas , transfer rnas  and small nucleolar rnas ; whereas those longer than  <dig> nucleotides are classified as long ncrnas .
lncrnas can be classified as lincrnas  that are transcribed adjacent to protein-coding genes, ernas , intronic lncrnas  and antisense lncrnas   <cit> .
in the past few years, an increasing number of lncrnas have been discovered in eukaryotic organisms, ranging from nematodes to humans .
for example, the most recent report of the encode  project published in september  <dig>  described  <dig>  lncrna loci in comparison to  <dig>  protein-coding genes in  <dig> human cell lines .
lncrnas can be polyadenylated or non-polyadenylated  <cit> , spliced or mono-exonic unspliced and the expression level of individual lncrnas is generally lower than the level of expression of the typical protein-coding mrnas , and some lncrnas have high tissue specificity .
despite the fact that only few lncrnas have been characterized experimentally in detail to date, it is already known that they can act via diverse mechanisms  <cit>  and can play regulatory and structural roles in almost every important biological process, such as x-chromosome inactivation and genomic imprinting, nuclear compartmentalization and architecture, cell fate specification, rna splicing, translational control, and chromatin modification  <cit> .
because of the key role of lncrnas in regulation of gene expression and therefore possible impact on phenotypes, it is important to identify most lncrnas.
catalogues of lncrna have been established for many species, including cattle .
for example, huang et al.
have identified a total of  <dig> putative lncrnas located in  <dig> intergenic regions using public bovine-specific expressed sequence tags sequences  <cit> .
more recently, weikard et al.
predicted more than  <dig>  potential lncrnas in bovine skin using rna-seq data  <cit> .
the current number of bovine lncrnas identified is rather low compared to more than  <dig>  lncrnas found in human, suggesting that more efforts are needed to discover all bovine lncrnas.
in the present study, we identified lncrnas 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 nine different limousin bulls were analysed.
we have identified more than  <dig> different lincrnas and  <dig> out of the  <dig> selected lincrnas were validated experimentally.
the rna-seq data and the collection of newly discovered lincrnas improve the genomic resources available for cattle, especially for beef breeds.
this collection of lincrnas may prove useful to study their link with genetic variability of meat quality traits.
RESULTS
rna sequencing and assembly of a muscle transcriptome
to identify lincrnas expressed in the bovine longissimus thoraci, we used paired-end rna sequencing  from nine limousin bull calves.
we used already published data from three limousin animals  <cit>  and poly-enriched mrna from six new limousin bull calves were retrotranscribed and subjected to high-throughput sequencing.
the six rna-seq libraries were barcode-tagged and sequenced on two lanes  of an illumina hiseq <dig> sequencer.
the reads were then de-multiplexed to assign reads to each sequenced sample according to its barcode index.
sequencing of all nine cdna libraries generated a total of  <dig> , <dig> raw paired-end reads with a length of  <dig> bases, resulting in a total of  <dig> gigabases.
approximately  <dig> to  <dig> million paired-end reads were obtained for each library.
the reads were then aligned using tophat  <cit>  onto the bovine umd <dig>  reference genome sequence.
65% to 75% of the reads were aligned onto the bovine genome, and 82% to 92% of the mapped reads were aligned properly paired .
transcripts were reconstructed using cufflinks  <cit> , resulting in the assembly of a total number of  <dig>  transcripts  with at least one paired-end read.table  <dig>
summary of reads mapping to the bovine transcriptomes
similar rna-seq read mapping rates were obtained in other rna-seq bovine studies .
for example, baldwin and collaborators found by sequencing the rumen epithelium that ~71% of the reads uniquely mapped to specific regions of the bovine genome  <cit> .
interestingly a comparable number of genes has been detected in bovine skin in a rna-seq project using a similar sequencing coverage and bioinformatics pipeline  <cit> .
raw gene expression levels were estimated by measuring the normalised count number for each transcript .
the five most frequent transcripts are shown in table  <dig>  these five genes  represent nearly 20% of all sequencing reads mapped to the bovine genome and are all associated with muscle structure.
these results were consistent with the physiological role of genes expected in the surveyed tissue.table  <dig>
top five transcripts with most assigned reads
not all genes were expressed among the nine selected samples.
transcripts corresponding to  <dig>  different genes were detected in all nine samples, while approximately 32% of the genes were expressed in only some of the samples, including  <dig>  different genes  only expressed in one sample.
transcript models predicted with cufflinks for each sample were sorted into different categories using the bovine genome annotation .
approximately between 57-66%  of the transcripts correspond to already known transcripts.
surprisingly, between ~18-22%  of the transcript models generated correspond to novel isoforms of known genes presumably from alternative splicing events .
this indicates that the bovine genome remains poorly annotated and that a large number of new transcript isoforms are still to be described.
interestingly, more than 2-4%  of the transcript models are predicted as unknown intergenic transcripts .
it is notable that some transcripts are classified into different categories.
for example, a transcript model could be part of an already known transcript and at the same time it could be part of the novel isoform category, as the rna-seq data indicates a novel exon.table  <dig>
number of mapped sequencing reads for each different class of assembled transcripts
code
e
=
x
s
j
c
p
u
.
o
i
e, single exon overlapping a reference exon and at least 10 bp of a reference intron indicating a possible pre-mrna fragment; =, complete match of intron chain; x, exonic overlap with reference on the opposite strand; s, an intron of the transfrag overlaps a reference intron on the opposite strand ; j, potentially novel isoform; c, contained in reference; p, possible polymerase run-on fragment ; u, unknown intergenic transcript; ., tracking file only, indicates multiple classifications; o, generic exonic overlap with a reference transcript; i, a single exon transcript falling entirely within a reference intron.
identification of putative lincrnas
to identify lncrnas, we developed a stringent filtering pipeline to discard transcripts with evidence for protein-coding potential.
we identify putative lncrnas by considering their open reading frame, their phylogenetic conservation across species and homology with known proteins and protein domains.
the reads were generated from non-directional rna-seq libraries, we therefore focus our effort on unknown intergenic transcripts and therefore could only identify putative lincrnas.
first, a minimal transcript size criterion was applied.
transcripts with multiple exons and larger than 200 nt were used.
this analysis resulted in the identification of  <dig>  putative multi-exonic intergenic transcripts .
second, we used phylocsf to score the coding potential of unknown multi-exonic transcripts using multi-species alignments.
phylocsf scores were calculated for the  <dig>  putative multi-exonic intergenic transcripts and two control sets, one of  <dig>  known protein-coding genes found in our rna-seq libraries and one with  <dig> already known bovine lncrnas  <cit> .
we set the phylocsf threshold empirically to a value retaining 73% of the known bovine ncrnas while removing 7% of protein-coding transcripts.
this filter retained  <dig>  putative non-coding transcripts .
third, we used cpat on the same  <dig>  putative multi-exonic intergenic transcripts, in order to assess their coding potential with a second prediction method.
to determine the optimum cut-off value, cpat was trained using a set of  <dig>  bovine known protein-encoding transcripts and a set of  <dig>  bovine non-coding sequences and a 10-fold cross-validation analysis was performed to estimate the prediction accuracy.
a cut-off value of  <dig>  was selected, maximising specificity and sensitivity  .
this procedure identified  <dig>  transcripts  as potential non-coding rnas.
the intersection of phylocsf and cpat predicted  <dig>  transcripts  as potential non-coding rnas.
finally, we removed any remaining transcripts of uncertain coding potential that had similarity to known proteins or protein domains recorded in the pfam database.
the resulting set contained  <dig>  transcript models corresponding to  <dig> different putative non-coding genes .
there is a possibility that some real lincrnas would have not been detected because of our stringent selection criteria.
for example, some real lincrnas could be lost owing to the chosen minimum transcript size.
we might also have missed some lowly expressed lincrna genes due to our moderate sequencing depth.
comparison of the genomic position of the  <dig> different genes encoding putative lincrnas found in bovine longissimus thoraci with mapping positions of previously identified bovine lncrnas publicly available in the noncode database   <cit>  show that  <dig>  of our lincrnas overlap with previously described bovine non-coding rna genes .
characterisation of identified lincrnas
the chromosomal location of the genes encoding these  <dig> putative lincrnas is presented in the additional file 4: figure s <dig>  the chromosomal distribution usually reflects the gene content of the chromosomes: larger chromosomes have more lincrna loci than shorter chromosomes.
for example, we found  <dig> genes encoding putative lincrnas on bta <dig>  whereas only  <dig> on btax.
previous studies have shown that genes encoding lncrnas are shorter in length, have shorter transcripts and have fewer exons than protein-coding genes  <cit> .
to determine whether the bovine muscle lincrnas we detected have the same features, we compared the size of  <dig> lincrna genes to  <dig>  protein-encoding genes detected in our rna-seq data.
the lincrnas represent much shorter gene length on average than protein-encoding genes .
we compared also the size of the  <dig>  transcripts corresponding to the  <dig> lincrna genes to the size of  <dig>  transcripts of the selected protein-encoding genes.
the lincrnas have shorter transcript on average than protein-encoding genes .
in addition, lincrna genes show also fewer exons than protein-encoding genes .
to determine whether the bovine muscle lincrnas we detected have the same expression feature, we compared the normalised quantified expression levels of the  <dig> lincrnas to that of the  <dig>  known selected protein-encoding genes.
our comparison indicates that the identified bovine lincrnas do show significant lower expression than the protein-encoding genes .
previous studies also showed that lincrnas are expressed at significantly lower levels than are protein-coding transcripts  <cit> .
lincrna gene expression
not all  <dig> lincrna genes were expressed among the nine selected samples.
<dig> genes were detected in all samples, while  <dig> genes were detected in less than five samples .
the highest number of putative lincrnas  was found in sample lim <dig>
the sequencing reads derived from the  <dig> different putative lincrna made up a bit more than 1% of all the paired-end reads mapped onto the bovine genome.
three genes encoding putative lincrnas had each more than  <dig> % of the total mapped paired-end sequencing reads.
the gene encoding a putative lincrna with the most sequencing reads  was represented by  <dig> % of the total sequencing reads.
these lincrnas with a relatively high expression level in bovine longissimus thoraci might play an important role in muscle function.
as our study identified lincrnas from different samples but from the same tissue, it is possible to analyse the inter-individual variability of the gene expression level of the identified lincrnas.
we calculated the coefficient of variation for each lincrna, from the expression level measured with the nine different longissimus thoraci samples and compared it with the gene expression variation of the selected protein-coding genes.
interestingly, there was a significant difference when we compared the expression variation of the  <dig> lincrnas with the  <dig>  known selected protein-encoding genes found in all nine muscle samples .
the higher averaged expression variation of the identified lincrnas suggests a loosened gene expression regulation of these genes compared to protein-encoding genes.
to our knowledge, it is the first time that this is described.
the observed higher gene expression variation of lincrna genes could be due to a higher number of polymorphisms within the regulatory regions.
to support this hypothesis we compared the snp density  of the regulatory region of the lincrna genes and protein-coding genes.
we mapped snps from ensembl  to the upstream  and downstream  regions of the lincrna genes and protein-encoding genes.
there was no significant difference when we compared the averaged snp density in the downstream regions of the  <dig> lincrnas with the one of the  <dig>  known selected protein-encoding genes .
in addition, there was no significant difference when we compared the averaged snp density in the upstream regions of the group of genes .
the lack of difference in snp densities do not rule out that the observed higher gene expression variation of lincrna genes could be due to the effect of polymorphisms within the regulatory regions.
however, the snp densities should be determined after sequencing the whole-genome or the regulatory region of the lincrna genes and protein-coding genes of the nine animals for which the expression data was generated.
the expression variation of the identified genes might be affected by the moderate sequencing depth obtained for some samples.
more work is therefore required to confirm the higher averaged expression variation seen with lincrna genes compared to protein-coding genes.
co-expression analysis
we have in our study the expression levels of the identified lincrnas and of known protein-encoding genes in nine samples from the same tissue.
it is therefore possible to analyse the co-expression of lincrnas with protein-encoding genes.
using the normalised expression levels, we calculated the spearman's rank correlation coefficient for each lincrna with each protein-encoding gene.
we found after correction for multiple testing  <dig>  pairs of lincrna/protein-encoding genes showing a highly correlated expression  .
<dig> different lincrna genes and  <dig> different protein-encoding genes had their expression highly correlated.
<dig> lincrnas showed correlations with more than  <dig> protein-encoding genes.
four lincrna genes  had the most correlations, with  <dig> different protein-encoding genes.
all lincrna/protein-encoding gene pairs showing highly correlated expression had the same correlation and p-values .
no anti-correlated lincrna/protein-coding gene pairs were found.
there were neither no pairs of lincrna/protein-encoding genes located at less than 2 mb apart, suggesting a lack of cis-regulation among the paired genes we detected in bovine muscle.
however, we found  <dig>  co-expression correlations without the bonferroni correction for multiple testing  between  <dig> different lincrnas and  <dig>  different protein-encoding genes, including  <dig>  cis and  <dig>  anti- correlations.
validation of novel lincrnas
to confirm that the identified bovine lincrnas are transcribed in vivo,  <dig> lincrnas detected in all  <dig> samples were randomly selected for rt-pcr validation.
<dig> out of the  <dig> selected lincrnas could be amplified using total rna from longissimus thoraci, as shown in figure  <dig>  all amplification products have the expected size; however for lincrna xloc_ <dig> we obtain an extra band of ~ <dig>  bp.
the high percentage of validation suggests that most putative lincrnas might be truly expressed in vivo.figure  <dig>
validation of selected lincrnas using rt-pcr.
selected lincrnas are numbered from  <dig> to  <dig> as detailed in additional file 6: table s <dig>  the first and second lanes are with pcr products using cdna or without cdna , respectively.
*indicates faint but specific amplification products.
functional lincrna candidates
the positions of the  <dig> genes encoding putative linrnas were compared to the position on the umd <dig>  bovine genome assembly of know quantitative trait loci  deposited in the public database animalqtldb  <cit> .
<dig> loci were located in  <dig>  different qtl regions, including  <dig> lincrnas within  <dig> qtls for meat quality/muscle-related traits .
for example,  <dig> different putative lincrnas are found in  <dig> qtl regions for meat tenderness; whereas  <dig> putative lincrnas are within  <dig> qtls for marbling score.
qtls were sorted into two groups  and the number of lincrnas found in these two groups were counted.
we then performed a chi-squared test and found a significant difference  in the number of lincrnas between the two groups , suggesting an enrichment of snps in meat/muscle related qtls.
the high number of putative lincrnas located within known qtl regions, particularly in chromosomal regions harbouring qtls for meat quality-related traits, indicates that the collection of lincrnas found in the longissimus thoraci transcriptome may prove useful to elucidate their effect on mechanisms underlying the genetic variability of meat quality traits.
CONCLUSIONS
the present study represents the first analysis of large intergenic non-coding genes discovered in bovine muscle.
using a computational pipeline that we developed to analyse rna-seq data, we identified  <dig> different novel putative lincrnas.
we could validate by rt-pcr  <dig> out of fourteen selected putative lincrnas, suggesting that most putative lincrnas might be truly expressed in vivo.
the identified putative bovine lincrna genes share most features with mammalian counterparts.
as our study identified lincrnas in longissimus thoraci from nine different limousin animals, we could analyse the variability of the gene expression level of the identified lincrnas and compared it with the gene expression variation of the known selected protein-coding genes.
we found a higher averaged expression variation for the identified lincrnas suggesting a loosened gene expression regulation.
we also analysed the co-expression of lincrnas with protein-encoding genes and found  <dig>  pairs of lincrna/protein-encoding genes showing a highly correlated expression.
some of the lincrnas expressed in muscle are located within quantitative trait loci for meat quality traits.
future experimental studies are required to unravel the function of these molecules and to elucidate their effect on mechanisms underlying the genetic variability of meat quality traits.
the lincrnas identified here will complement the catalog of lincrnas already discovered in cattle and therefore will help to better annotate the bovine genome.

 over the past decade, genome-wide transcriptional studies discovered that a large fraction of the eukaryotic genomes is transcribed in a heterogeneous population of noncoding rnas. these are transcripts that are not translated into a protein but act as functional rnas. ncrnas shorter than nucleotides are usually identified as small/short ncrna and include piwi-interacting rnas, endogeneous small interfering rnas and micrornas but also classical ncrnas, such as ribosomal rnas, transfer rnas and small nucleolar rnas; whereas those longer than nucleotides are classified as long ncrnas. lncrnas can be classified as lincrnas that are transcribed adjacent to protein-coding genes, ernas, intronic lncrnas and antisense lncrnas. in the past few years, an increasing number of lncrnas have been discovered in eukaryotic organisms, ranging from nematodes to humans. for example, the most recent report of the encode project published in september described lncrna loci in comparison to protein-coding genes in human cell lines. lncrnas can be polyadenylated or non-polyadenylated, spliced or mono-exonic unspliced and the expression level of individual lncrnas is generally lower than the level of expression of the typical protein-coding mrnas, and some lncrnas have high tissue specificity. despite the fact that only few lncrnas have been characterized experimentally in detail to date, it is already known that they can act via diverse mechanisms and can play regulatory and structural roles in almost every important biological process, such as x-chromosome inactivation and genomic imprinting, nuclear compartmentalization and architecture, cell fate specification, rna splicing, translational control, and chromatin modification. because of the key role of lncrnas in regulation of gene expression and therefore possible impact on phenotypes, it is important to identify most lncrnas. catalogues of lncrna have been established for many species, including cattle. for example, huang et al. have identified a total of putative lncrnas located in intergenic regions using public bovine-specific expressed sequence tags sequences. more recently, weikard et al. predicted more than potential lncrnas in bovine skin using rna-seq data. the current number of bovine lncrnas identified is rather low compared to more than lncrnas found in human, suggesting that more efforts are needed to discover all bovine lncrnas. in the present study, we identified lncrnas 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 nine different limousin bulls were analysed. we have identified more than different lincrnas and out of the selected lincrnas were validated experimentally. the rna-seq data and the collection of newly discovered lincrnas improve the genomic resources available for cattle, especially for beef breeds. this collection of lincrnas may prove useful to study their link with genetic variability of meat quality traits.  rna sequencing and assembly of a muscle transcriptome to identify lincrnas expressed in the bovine longissimus thoraci, we used paired-end rna sequencing from nine limousin bull calves. we used already published data from three limousin animals and poly-enriched mrna from six new limousin bull calves were retrotranscribed and subjected to high-throughput sequencing. the six rna-seq libraries were barcode-tagged and sequenced on two lanes of an illumina hiseq sequencer. the reads were then de-multiplexed to assign reads to each sequenced sample according to its barcode index. sequencing of all nine cdna libraries generated a total of, raw paired-end reads with a length of bases, resulting in a total of gigabases. approximately to million paired-end reads were obtained for each library. the reads were then aligned using tophat onto the bovine umd reference genome sequence. 65% to 75% of the reads were aligned onto the bovine genome, and 82% to 92% of the mapped reads were aligned properly paired. transcripts were reconstructed using cufflinks, resulting in the assembly of a total number of transcripts with at least one paired-end read.table summary of reads mapping to the bovine transcriptomes similar rna-seq read mapping rates were obtained in other rna-seq bovine studies. for example, baldwin and collaborators found by sequencing the rumen epithelium that ~71% of the reads uniquely mapped to specific regions of the bovine genome. interestingly a comparable number of genes has been detected in bovine skin in a rna-seq project using a similar sequencing coverage and bioinformatics pipeline. raw gene expression levels were estimated by measuring the normalised count number for each transcript. the five most frequent transcripts are shown in table  these five genes represent nearly 20% of all sequencing reads mapped to the bovine genome and are all associated with muscle structure. these results were consistent with the physiological role of genes expected in the surveyed tissue.table top five transcripts with most assigned reads not all genes were expressed among the nine selected samples. transcripts corresponding to different genes were detected in all nine samples, while approximately 32% of the genes were expressed in only some of the samples, including different genes only expressed in one sample. transcript models predicted with cufflinks for each sample were sorted into different categories using the bovine genome annotation. approximately between 57-66% of the transcripts correspond to already known transcripts. surprisingly, between ~18-22% of the transcript models generated correspond to novel isoforms of known genes presumably from alternative splicing events. this indicates that the bovine genome remains poorly annotated and that a large number of new transcript isoforms are still to be described. interestingly, more than 2-4% of the transcript models are predicted as unknown intergenic transcripts. it is notable that some transcripts are classified into different categories. for example, a transcript model could be part of an already known transcript and at the same time it could be part of the novel isoform category, as the rna-seq data indicates a novel exon.table number of mapped sequencing reads for each different class of assembled transcripts code e = x s j c p u. o i e, single exon overlapping a reference exon and at least 10 bp of a reference intron indicating a possible pre-mrna fragment; =, complete match of intron chain; x, exonic overlap with reference on the opposite strand; s, an intron of the transfrag overlaps a reference intron on the opposite strand; j, potentially novel isoform; c, contained in reference; p, possible polymerase run-on fragment; u, unknown intergenic transcript;., tracking file only, indicates multiple classifications; o, generic exonic overlap with a reference transcript; i, a single exon transcript falling entirely within a reference intron. identification of putative lincrnas to identify lncrnas, we developed a stringent filtering pipeline to discard transcripts with evidence for protein-coding potential. we identify putative lncrnas by considering their open reading frame, their phylogenetic conservation across species and homology with known proteins and protein domains. the reads were generated from non-directional rna-seq libraries, we therefore focus our effort on unknown intergenic transcripts and therefore could only identify putative lincrnas. first, a minimal transcript size criterion was applied. transcripts with multiple exons and larger than 200 nt were used. this analysis resulted in the identification of putative multi-exonic intergenic transcripts. second, we used phylocsf to score the coding potential of unknown multi-exonic transcripts using multi-species alignments. phylocsf scores were calculated for the putative multi-exonic intergenic transcripts and two control sets, one of known protein-coding genes found in our rna-seq libraries and one with already known bovine lncrnas. we set the phylocsf threshold empirically to a value retaining 73% of the known bovine ncrnas while removing 7% of protein-coding transcripts. this filter retained putative non-coding transcripts. third, we used cpat on the same putative multi-exonic intergenic transcripts, in order to assess their coding potential with a second prediction method. to determine the optimum cut-off value, cpat was trained using a set of bovine known protein-encoding transcripts and a set of bovine non-coding sequences and a 10-fold cross-validation analysis was performed to estimate the prediction accuracy. a cut-off value of was selected, maximising specificity and sensitivity. this procedure identified transcripts as potential non-coding rnas. the intersection of phylocsf and cpat predicted transcripts as potential non-coding rnas. finally, we removed any remaining transcripts of uncertain coding potential that had similarity to known proteins or protein domains recorded in the pfam database. the resulting set contained transcript models corresponding to different putative non-coding genes. there is a possibility that some real lincrnas would have not been detected because of our stringent selection criteria. for example, some real lincrnas could be lost owing to the chosen minimum transcript size. we might also have missed some lowly expressed lincrna genes due to our moderate sequencing depth. comparison of the genomic position of the different genes encoding putative lincrnas found in bovine longissimus thoraci with mapping positions of previously identified bovine lncrnas publicly available in the noncode database show that of our lincrnas overlap with previously described bovine non-coding rna genes. characterisation of identified lincrnas the chromosomal location of the genes encoding these putative lincrnas is presented in the additional file 4: figure s the chromosomal distribution usually reflects the gene content of the chromosomes: larger chromosomes have more lincrna loci than shorter chromosomes. for example, we found genes encoding putative lincrnas on bta whereas only on btax. previous studies have shown that genes encoding lncrnas are shorter in length, have shorter transcripts and have fewer exons than protein-coding genes. to determine whether the bovine muscle lincrnas we detected have the same features, we compared the size of lincrna genes to protein-encoding genes detected in our rna-seq data. the lincrnas represent much shorter gene length on average than protein-encoding genes. we compared also the size of the transcripts corresponding to the lincrna genes to the size of transcripts of the selected protein-encoding genes. the lincrnas have shorter transcript on average than protein-encoding genes. in addition, lincrna genes show also fewer exons than protein-encoding genes. to determine whether the bovine muscle lincrnas we detected have the same expression feature, we compared the normalised quantified expression levels of the lincrnas to that of the known selected protein-encoding genes. our comparison indicates that the identified bovine lincrnas do show significant lower expression than the protein-encoding genes. previous studies also showed that lincrnas are expressed at significantly lower levels than are protein-coding transcripts. lincrna gene expression not all lincrna genes were expressed among the nine selected samples. genes were detected in all samples, while genes were detected in less than five samples. the highest number of putative lincrnas was found in sample lim the sequencing reads derived from the different putative lincrna made up a bit more than 1% of all the paired-end reads mapped onto the bovine genome. three genes encoding putative lincrnas had each more than % of the total mapped paired-end sequencing reads. the gene encoding a putative lincrna with the most sequencing reads was represented by % of the total sequencing reads. these lincrnas with a relatively high expression level in bovine longissimus thoraci might play an important role in muscle function. as our study identified lincrnas from different samples but from the same tissue, it is possible to analyse the inter-individual variability of the gene expression level of the identified lincrnas. we calculated the coefficient of variation for each lincrna, from the expression level measured with the nine different longissimus thoraci samples and compared it with the gene expression variation of the selected protein-coding genes. interestingly, there was a significant difference when we compared the expression variation of the lincrnas with the known selected protein-encoding genes found in all nine muscle samples. the higher averaged expression variation of the identified lincrnas suggests a loosened gene expression regulation of these genes compared to protein-encoding genes. to our knowledge, it is the first time that this is described. the observed higher gene expression variation of lincrna genes could be due to a higher number of polymorphisms within the regulatory regions. to support this hypothesis we compared the snp density of the regulatory region of the lincrna genes and protein-coding genes. we mapped snps from ensembl to the upstream and downstream regions of the lincrna genes and protein-encoding genes. there was no significant difference when we compared the averaged snp density in the downstream regions of the lincrnas with the one of the known selected protein-encoding genes. in addition, there was no significant difference when we compared the averaged snp density in the upstream regions of the group of genes. the lack of difference in snp densities do not rule out that the observed higher gene expression variation of lincrna genes could be due to the effect of polymorphisms within the regulatory regions. however, the snp densities should be determined after sequencing the whole-genome or the regulatory region of the lincrna genes and protein-coding genes of the nine animals for which the expression data was generated. the expression variation of the identified genes might be affected by the moderate sequencing depth obtained for some samples. more work is therefore required to confirm the higher averaged expression variation seen with lincrna genes compared to protein-coding genes. co-expression analysis we have in our study the expression levels of the identified lincrnas and of known protein-encoding genes in nine samples from the same tissue. it is therefore possible to analyse the co-expression of lincrnas with protein-encoding genes. using the normalised expression levels, we calculated the spearman's rank correlation coefficient for each lincrna with each protein-encoding gene. we found after correction for multiple testing pairs of lincrna/protein-encoding genes showing a highly correlated expression. different lincrna genes and different protein-encoding genes had their expression highly correlated. lincrnas showed correlations with more than protein-encoding genes. four lincrna genes had the most correlations, with different protein-encoding genes. all lincrna/protein-encoding gene pairs showing highly correlated expression had the same correlation and p-values. no anti-correlated lincrna/protein-coding gene pairs were found. there were neither no pairs of lincrna/protein-encoding genes located at less than 2 mb apart, suggesting a lack of cis-regulation among the paired genes we detected in bovine muscle. however, we found co-expression correlations without the bonferroni correction for multiple testing between different lincrnas and different protein-encoding genes, including cis and anti- correlations. validation of novel lincrnas to confirm that the identified bovine lincrnas are transcribed in vivo, lincrnas detected in all samples were randomly selected for rt-pcr validation. out of the selected lincrnas could be amplified using total rna from longissimus thoraci, as shown in figure  all amplification products have the expected size; however for lincrna xloc_ we obtain an extra band of ~  bp. the high percentage of validation suggests that most putative lincrnas might be truly expressed in vivo.figure validation of selected lincrnas using rt-pcr. selected lincrnas are numbered from to as detailed in additional file 6: table s the first and second lanes are with pcr products using cdna or without cdna, respectively. *indicates faint but specific amplification products. functional lincrna candidates the positions of the genes encoding putative linrnas were compared to the position on the umd bovine genome assembly of know quantitative trait loci deposited in the public database animalqtldb. loci were located in different qtl regions, including lincrnas within qtls for meat quality/muscle-related traits. for example, different putative lincrnas are found in qtl regions for meat tenderness; whereas putative lincrnas are within qtls for marbling score. qtls were sorted into two groups and the number of lincrnas found in these two groups were counted. we then performed a chi-squared test and found a significant difference in the number of lincrnas between the two groups, suggesting an enrichment of snps in meat/muscle related qtls. the high number of putative lincrnas located within known qtl regions, particularly in chromosomal regions harbouring qtls for meat quality-related traits, indicates that the collection of lincrnas found in the longissimus thoraci transcriptome may prove useful to elucidate their effect on mechanisms underlying the genetic variability of meat quality traits.  the present study represents the first analysis of large intergenic non-coding genes discovered in bovine muscle. using a computational pipeline that we developed to analyse rna-seq data, we identified different novel putative lincrnas. we could validate by rt-pcr out of fourteen selected putative lincrnas, suggesting that most putative lincrnas might be truly expressed in vivo. the identified putative bovine lincrna genes share most features with mammalian counterparts. as our study identified lincrnas in longissimus thoraci from nine different limousin animals, we could analyse the variability of the gene expression level of the identified lincrnas and compared it with the gene expression variation of the known selected protein-coding genes. we found a higher averaged expression variation for the identified lincrnas suggesting a loosened gene expression regulation. we also analysed the co-expression of lincrnas with protein-encoding genes and found pairs of lincrna/protein-encoding genes showing a highly correlated expression. some of the lincrnas expressed in muscle are located within quantitative trait loci for meat quality traits. future experimental studies are required to unravel the function of these molecules and to elucidate their effect on mechanisms underlying the genetic variability of meat quality traits. the lincrnas identified here will complement the catalog of lincrnas already discovered in cattle and therefore will help to better annotate the bovine genome. 
interestingly, there was a significant difference when we compared the expression variation of the lincrnas with the known selected protein-encoding genes found in all muscle samples. using a computational pipeline, we defined a stringent set of different lincrnas with lincrnas found in all nine muscle samples. as for the first time, our study identified lincrnas from nine different samples from the same tissue, it is possible to analyse the inter-individual variability of the gene expression level of the identified lincrnas. to identify lincrnas, we used paired-end rna sequencing to explore the transcriptomes of longissimus thoraci from nine limousin bull calves. some of the lincrnas expressed in muscle are located within quantitative trait loci for meat quality traits. it may prove useful to elucidate their effect on mechanisms underlying the genetic variability of meat quality traits. in the present study, we performed a survey in bovine longissimus thoraci of lincrnas. bovine lincrnas share characteristics seen in their mammalian counterparts: relatively short transcript and gene lengths, low exon number and significantly lower expression, compared to protein-encoding genes. to our knowledge, this represents the first such study in bovine muscle. our study provides a glimpse into the lincrna content of bovine muscle and will facilitate future experimental studies to unravel the function of these molecules. this catalog will complement the list of lincrnas already discovered in cattle and therefore will help to better annotate the bovine genome.
 the advent of large-scale gene expression technologies has helped to reveal in eukaryotic cells, the existence of thousands of non-coding transcripts, whose function and significance remain mostly poorly understood. among these non-coding transcripts, long non-coding rnas are the least well-studied but are emerging as key regulators of diverse cellular processes. in the present study, we performed a survey in bovine longissimus thoraci of lincrnas. to our knowledge, this represents the first such study in bovine muscle.  to identify lincrnas, we used paired-end rna sequencing to explore the transcriptomes of longissimus thoraci from nine limousin bull calves. approximately 14– million paired-end reads were obtained per library. a total of different transcripts were identified. using a computational pipeline, we defined a stringent set of different lincrnas with lincrnas found in all nine muscle samples. bovine lincrnas share characteristics seen in their mammalian counterparts: relatively short transcript and gene lengths, low exon number and significantly lower expression, compared to protein-encoding genes. as for the first time, our study identified lincrnas from nine different samples from the same tissue, it is possible to analyse the inter-individual variability of the gene expression level of the identified lincrnas. interestingly, there was a significant difference when we compared the expression variation of the lincrnas with the known selected protein-encoding genes found in all muscle samples. in addition, we found pairs of lincrna/protein-encoding genes showing a highly significant correlated expression. fourteen lincrnas were selected and were validated by rt-pcr. some of the lincrnas expressed in muscle are located within quantitative trait loci for meat quality traits.  our study provides a glimpse into the lincrna content of bovine muscle and will facilitate future experimental studies to unravel the function of these molecules. it may prove useful to elucidate their effect on mechanisms underlying the genetic variability of meat quality traits. this catalog will complement the list of lincrnas already discovered in cattle and therefore will help to better annotate the bovine genome. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords cattlemusclerna-seqbeeflong non-coding rnaissue-copyright-statement© the author 2014 