BACKGROUND
the pig has been providing with large scale products for human consumption.
it also services as an important animal model for human diseases including obesity, diabetes, cardiovascular disease and endocrinology because of the similarity in physiology, organ development and disease progression  <cit> .
however, until now, only several causative mutations have been isolated for porcine complex traits, e.g.
igf <dig> for muscle growth and fat deposit  <cit> , prkag <dig> for glycogen content in skeletal muscle  <cit>  and nr6a <dig> for number of vertebrae  <cit> .
elucidating the complexity of the pig transcriptome is not only essential for interpreting the functional elements of the genome, but also benefits the understanding of human related complex traits such as fat deposition, metabolism and growth.
for example, utr lengths are correlated with gene function, localization and requirements for regulation  <cit> ; alternative splicing plays a major part in biological complexity in humans  <cit> ; and non-coding rnas  are crucial for multiple biological processes  <cit> .
the global transcriptome profilings in abdominal fat, induced fat cells, muscle, liver and pituitary gland have been investigated in pigs by long sage analysis, full-length enriched cdna library and microarray  <cit> .
however, the complexity of the porcine transcriptome is not yet fully elucidated.
the massively parallel sequencing using next generation sequencer  provides a huge potential to revolutionize the field of pig transcriptome, owning to its abilities to discover extensive alternative splicing and identify large-scale novel transcripts at single-nucleotide resolution  <cit> .
moreover, the paired-end sequencing strategy of rna-seq further improves sequencing efficiency and extends short read lengths for better understanding transcriptome  <cit> .
rna-seq has opened a new horizon for our understanding of global gene expression and has demonstrated the complexity of mammalian transcriptome vastly underestimated before.
mirnas are a class of small rnas that regulate gene expression by decreasing the target mrna levels or inhibiting translation of protein encoding transcripts.
global mirna abundance has been measured in skeletal muscle by microarray to evaluate the roles of mirnas in pig development and meat production  <cit> .
mirbase  <dig>  database has collected  <dig> porcine mature mirnas.
rna-seq also provides valuable insights into mirna transcriptome, especially into those mirnas insufficiently detected by microarray.
until now, porcine mirna transcriptome has been investigated by next generation sequencer in intestine  <cit> , pre- and postnatal piglet  <cit> , developing brain  <cit>  and skeletal muscle  <cit> .
we herein performed a global transcriptome analysis on three types of tissues related to metabolism, meat production and fat deposition: liver, longissimus dorsi muscle  and abdominal fat  harvested from a full-sib f <dig> female pair with extreme phenotypes by rna-seq.
the results allowed us to investigate large-scale alternative splicing events, identify extensive number of novel transcript units, determine gene boundaries at single nucleotide resolution and comprehensively survey porcine micrornas in the tested tissues.
to our knowledge, this study presented the first systematical investigation on the complexity of porcine transcriptome with nucleotide resolution.
moreover, this study identified many important candidate genes related to growth, meat quality, blood physiological and biochemical parameters by the integrated analysis of genome-wide association study and differential gene expression.
RESULTS
sequencing and mapping of the porcine transcriptome
we sequenced cdna libraries from  <dig> tissues of the full-sib female pair using high-seq  <dig> at bgi-shenzhen, china.
the data set was analyzed according to the bgi bioinformatics protocols for rna-seq.
the sequence reads have been submitted to the ncbi gene expression omnibus under accession no.
gse <dig>  in total, we acquired  <dig> ,956- <dig> , <dig> paired-end reads of  <dig> bp.
the total read length was  <dig>  gigabases , representing about 8-fold of the porcine genome size.
we technically replicated the rna-seq experiments in the  <dig> samples with  <dig>  <r <dig> <  <dig>  .
furthermore, the expression patterns of  <dig> randomly selected transcripts between two individuals were validated by qrt-pcr with a relative coefficient of r <dig> =  <dig>  .
the technical replicates and qrt-pcr confirmed the high reproducibility of rna-seq in this study.
we aligned all short reads onto the whole reference genome .
tolerances were set to allow at most two mismatches for  <dig> bp reads in each alignment.
about  <dig> - <dig> % of reads were mapped to the pig reference genome, of which  <dig> - <dig> % fell in annotated exons;  <dig> - <dig> % located in introns;  <dig> - <dig> % overlapped with exons, and the remaining  <dig> - <dig> % were assigned to intergenic regions .
total  <dig> - <dig> % of reads had a uniquely genomic location, and  <dig> - <dig> % of reads corresponded to reference genes with  <dig> - <dig> % of uniquely matched reads.
unmapped or multi-position matched reads  were excluded from further analyses .
identification of an extensive number of novel transcript units
the uniquely mapped reads  covered  <dig> %  of the annotated transcripts in ucsc pig genome database  <cit>  by at least one sequence read.
a total of  <dig>  transcripts were expressed in all three tissues, and  <dig>   <dig> and  <dig>  transcripts were discovered exclusively in af, ld and liver, respectively .
we quantified the gene expression level by counting the number of reads per kilobase per million mapped reads .
about  <dig> - <dig> % of the annotated transcripts showed expression with >  <dig>  rpkm.
the obtained transcription fragments of more than  <dig> % of the detected transcripts were >  <dig> bp in length .
the percentage of a gene covered by reads was defined as gene coverage.
extensive read mapping revealed about  <dig> - <dig> % of the detected transcripts with more than  <dig> % in gene coverage.
we detected an extensive number of novel transcript units by the procedures described in methods and zhang et al.
<cit> .
in total, we obtained  <dig> ,  <dig> ,  <dig> ,  <dig> ,  <dig>  and  <dig>  novel transcript units, respectively, in the six tested samples with a mean length of  <dig> bp and a size range from  <dig> to  <dig>  bp .
many novel transcript units  have >  <dig> exons and the largest one is  <dig>  bp in length containing  <dig> exons .
the novel transcript units were often identified in clusters, indicating that closely spaced transcript units are likely to merge into longer transcripts with increasing sequencing depth.
by scanning the genomic location of each novel transcribed unit, we clustered adjacent fragments  into one transcriptionally active region.
in this way, we identified  <dig> - <dig>  novel transcriptionally active regions in the six tested samples .
a proportion  of the regions was comprised of single exon.
to investigate whether the identified novel transcript units were non-coding rnas, we aligned the sequences of the novel transcript units to non-coding rna precursor sequences in rfam database  <cit> .
we found that  <dig> - <dig> % of novel transcript units were non-coding rna precursors.
of these,  <dig> % were in average trna precursors,  <dig> % were mirna precursors and  <dig> % were snorna precursors.
alternative splicing events in pig transcriptome
to investigate alternative splicing, we identified the sequence reads that were mapped to the regions of computationally determined theoretical splicing junctions.
four known types of alternative splicing models including alternative 3' splicing site, alternative 5' splicing site, exon skipping and intron retention were considered in this study.
the distribution of alternative splicing events is shown in additional file  <dig>  the alternative splicing events from  <dig> tissues of one individual were pooled and the redundancy was removed to get a final set of alternative splicing events.
we found that up to  <dig>  genes accounting for  <dig> % of the known genes had undergone alternative splicing in these tissues, displaying  <dig>  alternative events in individual  <dig>  similarly, the number of alternative splicing genes and events in individual  <dig> were  <dig>   and  <dig> , respectively .
figure 2a and 2b show an example of intron retention and alternative 5' splicing site.
we found that about  <dig> % of the alternative spliced genes underwent multiple alternative splicing events , illustrating the complexity of porcine transcriptome.
alternative 3' splicing site is the most common type of alternative splicing events accounting for  <dig> % of all alternative splicing events in pigs, while intron retention and exon-skipping only constituted  <dig> % and  <dig> % in individual  <dig>  or  <dig> % and  <dig> % in individual  <dig>  the average size of retained introns is  <dig> bp, ranging from  <dig> bp to  <dig> bp.
extension of annotated gene boundaries
the extensions of 5' and 3' boundaries were determined by comparison of the gene models obtained by rna-seq with the existing gene annotations.
in the six tested samples, a total of 587- <dig> genes were extended at the 5' end, of which more than  <dig> % had an extension of at least  <dig> bp in length.
in comparison,  <dig> - <dig>  genes were extended at their 3' end, of which more than  <dig> % were extended by at least  <dig> bp.
furthermore, 411- <dig> genes were extended at both ends .
in individual  <dig>  total  <dig>  annotated genes were extended at the 5' end, of which  <dig> %  were observed in all three tissues; and  <dig>  annotated genes were extended at the 3' end, of which  <dig> %  had the extended 3' boundary in all three tissues.
similarly, in individual  <dig>  the percentage of the extended 5' or 3' genes shared in all three tissues was  <dig> %  and  <dig> % , respectively.
comprehensive survey of porcine microrna by deep sequencing
to get a comprehensive view of mirna transcriptome in pigs, we carried out deep sequencing of small rna  using the tested samples mentioned above.
the obtained mirna sequence reads have been submitted to the ncbi gene expression omnibus under accession no.
gse <dig>  we obtained a total of  <dig> ,703- <dig> , <dig> raw reads.
after removing low quality reads and corrupted adapter sequences,  <dig> ,608- <dig> , <dig> clean reads were retained for further analyses.
the majority of small rna was 20- <dig> nt for all libraries , with  <dig> nt small rna being the most abundant , which is in agreement with the common nucleotide length of mirnas.
we mapped  <dig> - <dig> % of clean reads to the reference genome .
chromosomes  <dig>   <dig>   <dig>   <dig>   <dig> and  <dig> harbored clean reads of >  <dig> , <dig> tags .
after further removal of trna, rrna, scrna, snrna, snorna, exon rna, intron rna and repeat regions, a total of  <dig> ,788- <dig> , <dig> mirna sequences were obtained .
aligning mirna sequences to the porcine mature mirnas in mirbase  <dig>  database  <cit>  revealed  <dig>  - <dig>  % of mature mirnas in each sample.
we pooled all identified mature mirnas and found that only five mature mirnas in mirbase  <dig>  database were not detected in this study .
the expression levels of the identified mirnas displayed a very large range, as reflected by the number of sequence reads, which varied from single counts for rare mirnas to several hundred thousand reads for the most abundant mirnas .
as many as  <dig> - <dig> % of the clean reads belong to the mirnas ranked top  <dig> in expression levels in each sample .
to identify potential novel mirna, we further analyzed the small rna tags that could not match known mirnas and were mapped to intergenic or intronic regions of the reference genome.
the characteristic of hairpin structure of mirna precursor was used to predict novel mirna by exploring the secondary structure.
in total,  <dig> to  <dig> potential novel mirnas each supported by at least five sequence reads, were identified in the  <dig> tested samples .
the potential novel mirnas of  <dig> samples were pooled and the redundancy was removed to get a final set of  <dig> unique porcine putative novel mirnas.
the vast majorities of these mirnas were expressed at low levels.
nevertheless,  <dig> mirnas showed relatively high expression levels representing by more than  <dig>  sequence reads, and  <dig> by more than  <dig> sequence reads .
to determine the genomic locations of these potential novel mirnas, their precursor sequences were blasted against the porcine reference genome sequence .
as shown in figure 3c,  <dig> % of the novel mirnas were assigned to the reference genome.
most of the precursors were located on chromosome  <dig>  followed by chromosomes × and  <dig> , and none of the novel mirnas were mapped to chromosome  <dig>
cross-tissue comparison of differential transcription events
deep rna sequencing in three types of tissues allows us to investigate the tissue-differentially transcriptional events.
of the  <dig>  identified transcripts,  <dig>   <dig> and  <dig>  transcripts were discovered exclusively in af, ld and liver, respectively .
interestingly, as reported in humans  <cit> , the majority of alternative splicing events showed clear tissue specificity, demonstrating the importance of alternative splicing in tissue specific programs of gene expression and its major roles in functional complexity.
functional annotations of alternative splicing genes showed that they play important roles in their corresponding tissues.
the porcine ryr <dig> gene  which had the largest number of ld-specific alternative splicing events has been associated with malignant hyperthermia and has significant effects on pig meat quality and carcass leanness  <cit> .
the alb  gene having the most number of liver-specific alternative splicing events plays important roles in transportation of fatty acids  <cit> .
sla- <dig>  implicating in immune and type i diabetes  <cit>  had the most number of alternative splicing events in af that has been known as an important immune organ.
other tissue-specifically alternative splicing genes also show major roles in their corresponding tissues.
for instance, myhc-2a , ppap2c  and apoh  had tissue-specific alternative splicing in ld, af and liver, respectively, corresponding to their roles in normal muscle development and function  <cit> , converting phosphatidic acid to diacylglycerol  <cit> , transportation of fatty acids  <cit> , and lipoprotein metabolism  <cit>  .
the expression abundance of many mirnas also showed apparently tissue-differential patterns.
the mature mirnas with high abundance in each tissue are conserved in mammals and likely related to important biological functions.
mir- <dig> is the most abundant mirna in human liver and also had the most abundance of expressed read counts in porcine liver .
mir- <dig> plays a crucial role in cholesterol, fatty acid and lipid metabolism  <cit> .
the mir- <dig> and mir- <dig> are key mediators in proper skeletal and cardiac muscle development and function, myogenesis during embryonic development and muscle cell differentiation  <cit> .
the two mirnas, respectively, had the most and third abundance of expression level in ld in this study.
for those potential novel mirnas identified in this study,  <dig> of  <dig> putative novel mirnas were specifically identified in af,  <dig> and  <dig> in liver and ld, respectively .
their biological functions need further investigation.
differentially expressed transcripts and mirnas between two individuals
overall, there were clear linear relationships in the gene expression levels  between two individuals in all three tissues.
the number of unique reads mapped to different genes ranged from  <dig> to  <dig> .
the differentially expressed genes were selected based on the expression profiles and the following criteria:  if the fold change in gene expression levels between  <dig> and  <dig> was more than or equal to two fold  and  if the false discovery rate value was less than  <dig> .
with this, we identified  <dig> ,  <dig>  and  <dig> differentially expressed genes in af, liver and ld, respectively.
of these, the expressions of  <dig> ,  <dig> and  <dig> genes were up-regulated in  <dig> with respect to individual  <dig> .
to further investigate the biological relationships of differentially expressed genes with phenotypes, we performed the gene ontology  analysis by querying each differentially expressed gene into the records of the go database  <cit> .
the results of go functional annotations are presented in figure 5d-5f.
the main functional groups of differentially expressed genes in af are related to lipid metabolic process, alcohol metabolic process, lipid binding and protein homooligomerization.
the functions of differentially expressed genes in liver are enriched in metabolic process, catalytic activity and oxidoreductase activity.
and the differentially expressed genes in ld are mainly associated with metabolic process, cellular metabolic process, catalytic activity and oxidoreductase activity.
differentially expressed mirnas between two individuals were identified by comparing the normalized expression data of the mature mirnas.
in total,  <dig> differentially expressed mirnas  ≥  <dig> or fold-change  ≤ -1; p-value <  <dig> ) were identified in liver,  <dig> and  <dig> in ld and af, respectively.
most of the differentially expressed mirnas had relatively low expression levels .
interestingly, some differentially expressed mirnas are involved in the pathway relevant to development and diabetes.
for instance, mir- <dig> differentially expressed in ld enables precisely specific the muscle cell types by sharpening cellular responses to hedgehog in zebrafish  <cit> ; a liver-differentially expressed mirna of mir-10b is predicted to regulate genes in pathways relevant to type  <dig> diabetes  <cit> .
investigation of candidate genes for related phenotypes by integrated analysis of genome-wide association study  and differential gene expression
we selected total  <dig> most differentially expressed genes including  <dig> from liver and af, respectively,  <dig> from ld for further functional annotations .
only  <dig> genes could be mapped to the pig reference genome  and had functional annotations in mammals.
of the  <dig> genes,  <dig> have the description of phenotypes in knocked-out mice  <cit> .
a genome-wide association study using the pig 60k snp chip has been performed in the current white duroc × erhualian f <dig> resource population .
seven,  <dig> and  <dig> differentially expressed genes in liver, af and ld are located within  <dig>  mb around the snps that were most significantly associated with phenotypes in the gwas and the associated phenotypes were also observed in the corresponding gene-deficient mice .
here, we present the identification results of interesting candidate genes igf <dig>  cyp1a <dig>  ckm and ces <dig> for the related phenotypes by integrating the analysis of gwas and differential gene expression in figure  <dig>
confirmation of igf <dig> as a causative gene for heart weight
igf <dig> is the causative gene underlying the qtl for muscle growth, fat deposition and heart weight on ssc2p  <cit> .
a nucleotide substitution in intron  <dig> of igf <dig> abrogates interaction with zbed <dig>  resulting in threefold increase of igf <dig> messenger rna expression  <cit> .
in the current f <dig> population, a significant qtl for heart weight at day  <dig> was located at chr2:  <dig>  mb with a corrected p-value of  <dig>  × 10- <dig> in gwas .
this region overlapped with the igf <dig> locus.
in the rna-seq analysis, individual  <dig> with the heavier heart weight  showed a  <dig>  fold increase in igf <dig> expression level compared with individual  <dig>  the result was perfectly consistent with the previous result  <cit> , and is another example demonstrating the causality of this gene in the qtl.
cyp1a <dig> as a strong candidate gene for hemoglobin
two founder breeds of the f <dig> resource population are white duroc and chinese erhualian, which are divergent in hemoglobin concentration  <cit> .
blood parameters of hemoglobin and hematocrit at day  <dig> were measured in the f <dig> population.
we identified a significant qtl for both hemoglobin and hematocrit content at day  <dig>  which is located at chr7:  <dig>  mb in gwas .
a differentially expressed gene of cyp1a <dig> in liver by rna-seq is located within this qtl region.
the cyp1a <dig> knocked-out mice had the increased hemoglobin content  <cit> .
a quantitative trait transcript  analysis in this f <dig> population also revealed that the cyp1a <dig> expression level was associated with hemoglobin at day  <dig> .
these results strongly support cyp1a <dig> as a candidate for the qtl effect on hemoglobin content.
ckm as a candidate gene for pork ph value
creatine kinase, muscle  is a cytoplasmic enzyme involved in energy homeostasis.
ckm-deficient mice show abnormal muscle physiology and an increased skeletal muscle glycogen level  <cit> .
excess glycogen content in pig skeletal muscle is associated with low meat ph value, leading to bad meat quality  <cit> .
a significant qtl for meat ph value at  <dig> min,  <dig> h and  <dig> h after slaughter was detected around the region of chr2:  <dig>  mb with a p-value of  <dig>  × 10- <dig> in the f <dig> cross.
ckm gene locates at the chr2:  <dig>  mb and was differentially expressed in this full-sib pair with different meat ph values.
this indicates that ckm is an important candidate gene for meat ph value in pigs.
ces <dig> as an important candidate gene for serum cholesterol and triglyceride level
the liver is an important organ for lipid metabolism.
gene expression in liver influences the circulating cholesterol level.
ces <dig> that participates in fatty acyl and cholesterol ester metabolism  <cit>  was one of the differentially expressed genes in liver.
our previous qtt study in liver showed that the expression level of ces <dig> is significantly associated with serum total cholesterol  and triglyceride level  in the f <dig> population .
a significant qtl for both low density lipoprotein cholesterol and total cholesterol was identified at chr6:  <dig>  mb containing the ces <dig> gene in gwas.
the concordant results of qtt, gwas and differential expression suggest that ces <dig> is strong candidate gene for qtl affecting serum total cholesterol on chromosome  <dig>  the supporting evidence was also from the ces <dig> knocked-out mice that showed decreased circulating cholesterol level  <cit> .
discussion
in this study, we presented the systematical transcriptome profiling of pigs on three tissues related to metabolism, meat production and fat deposition using high throughput rna-seq technology.
this efficient deep sequencing not only allows us to analyze novel transcribed regions and mirnas, but also improves gene annotations at single nucleotide resolution.
integrated analysis of genome-wide association study and differential gene expression between two individuals revealed important candidate genes for related phenotypes.
a total of  <dig> ,956- <dig> , <dig> reads were obtained from rna-seq in six samples with  <dig> - <dig> % of reads uniquely mapped to the reference genome.
this data set provided comprehensive starting resources for improving the gene annotations across the porcine genome.
high r <dig> values between technical replicates showed high reproducibility of rna-seq in pigs.
it should be mentioned that there are still  <dig> - <dig> % of reads that cannot be matched to the reference genome.
this could be caused by low sequence coverage of the reference genome, reference errors, sequencing errors and defined mapping criterions.
a proportion  of mapped reads with multiple positions was discarded, attributing to known duplicated genes or chromosome segments.
a majority of the annotated transcripts  at ucsc database  <cit>  was covered by sequence reads, showing the sensitivity of rna-seq in transcript discovery even for lowly expressed genes  <cit> .
furthermore, we identified large number of novel transcript units which improved the gene annotations of the porcine genome and transcriptome.
alternative splicing is an important model of gene expression regulation and has not been generally accessible for microarray or sage methods in pigs.
some genes showed all four types of alternative splicing models , revealing the complexity of alternative splicing in pigs.
we found that alternative 3' splicing is the most common type of alternative splicing events in pigs.
this is in contrast to the report in human and yeast where exon-skipping is the most prevalent mechanism  <cit> , and it is also different from rice in which intron retention is the primary alternative splicing type  <cit> .
the percentage of alternative 5' or 3' splicing in total alternative splicing events in this study were higher than that reported in lim et al.
<cit>  where  <dig> % of alternative 5' splicing and 3% of alternative 3' splicing were observed.
we discovered that more than  <dig> % of the detected genes were alternatively spliced.
this number is much lower than the reported  <dig> % in human  <cit>  and  <dig> % in rice  <cit> .
three types of alternative splicing models including alternative first exon, alternative last exon and mutually exclusive exon were excluded from analyses because of currently unperfected algorithms.
more alternative spliced genes would be discovered if these alternative splicing models were considered.
as an important regulator of gene expression, mirna regulates the gene expression through decreasing the target mrna levels or repressing the translation  <cit> .
in this study, the high identification rate illustrates that the six small rna libraries from the tested tissues almost encompass the entire repertoire of known mirnas.
consistent with the results in li et al.
<cit> , porcine mature mirnas had a broad range of expression levels.
the highly expressed mirnas are known to have important regulatory functions in corresponding tissues.
the mir-1/ <dig> showed high abundance of expression levels in muscle .
interestingly, muscular hypertrophy in texel sheep has been shown to be caused by a mutation that creates an illegitimate binding site for mir-1/ <dig> in the 3' utr of the myostatin gene, leading to efficient translational inhibition of the myostatin gene and an increase in muscularity  <cit> .
the mir- <dig> had the most abundant expression level in liver.
it has the diversity of its roles in liver, e.g.
metabolism, hepatocarcinogenesis  <cit> .
just as in humans  <cit> , most of the potential novel mirnas discovered in this study are expressed at low levels , explaining why they were not discovered in previous efforts and showing the advantage of rna-seq in transcriptome analysis.
furthermore, we found that most of the novel mirnas were tissue-specifically identified .
we cautioned that some could be artifact due to the low expression level, the insufficient depth of sequencing and limited tested tissues.
mapping genetic factors that underlie quantitative traits in farm animals has been a challenging task  <cit> .
the recent wave of genome-wide association studies in human showed that a majority of snps associated with disease traits locate in regulated regions  <cit> .
integration of gene expression with genotype and phenotype data to elucidate the network of molecular interactions that underlie complex traits can facilitate the identification of variants that contribute to phenotypes  <cit> .
in this study, high-throughput sequencing allowed us to digitally discover an extensive number of differentially expressed genes in the full-sib pair with different phenotypes.
gene ontology analysis indicated that the differentially expressed genes had enrichment on functions related to metabolic process, catalytic activity and lipid binding.
the results suggest that these differentially expressed genes are likely related to the phenotypes of growth, metabolism or fat deposition.
by integrating gwas, differentially expressed genes and altered phenotypes in knocked-out mice, we found that many differentially expressed genes are the important candidate genes related to the phenotypes of serum cholesterol, growth traits, hemoglobin at day  <dig>  fatty acid level and muscle ph value .
for example, activating transcription factor  <dig>  gene differentially expressed in muscle and located in the qtl region for body weight at day  <dig> is related to decreased body weight in atf <dig> knocked-out mice   <cit> ; fatty acid binding protein  <dig>  located in the qtl region for fat deposition on ssc <dig> and related to increased white adipose tissue amount in fabp <dig> knocked-out mice was differentially expressed in abdominal fat in this study   <cit> .
the findings provide important clues for further dissecting of the responsible genes and variants.
CONCLUSIONS
this study provides a global view of the complexity of the pig mrna and mirna transcriptome, gives an extensive new knowledge about alternative splicing, novel transcript units, gene boundaries and novel mirnas in pigs.
the cross-tissue comparison identified lots of tissue-differential transcription events.
integrated analysis of gwas and differential gene expression allowed us to detect important candidate genes related to growth, meat quality, serum lipids and fatness.
the findings significantly enhance the current genome annotation of pigs and improve our understanding of complex traits.

 the pig has been providing with large scale products for human consumption. it also services as an important animal model for human diseases including obesity, diabetes, cardiovascular disease and endocrinology because of the similarity in physiology, organ development and disease progression. however, until now, only several causative mutations have been isolated for porcine complex traits, e.g. igf for muscle growth and fat deposit, prkag for glycogen content in skeletal muscle and nr6a for number of vertebrae. elucidating the complexity of the pig transcriptome is not only essential for interpreting the functional elements of the genome, but also benefits the understanding of human related complex traits such as fat deposition, metabolism and growth. for example, utr lengths are correlated with gene function, localization and requirements for regulation; alternative splicing plays a major part in biological complexity in humans; and non-coding rnas are crucial for multiple biological processes. the global transcriptome profilings in abdominal fat, induced fat cells, muscle, liver and pituitary gland have been investigated in pigs by long sage analysis, full-length enriched cdna library and microarray. however, the complexity of the porcine transcriptome is not yet fully elucidated. the massively parallel sequencing using next generation sequencer provides a huge potential to revolutionize the field of pig transcriptome, owning to its abilities to discover extensive alternative splicing and identify large-scale novel transcripts at single-nucleotide resolution. moreover, the paired-end sequencing strategy of rna-seq further improves sequencing efficiency and extends short read lengths for better understanding transcriptome. rna-seq has opened a new horizon for our understanding of global gene expression and has demonstrated the complexity of mammalian transcriptome vastly underestimated before. mirnas are a class of small rnas that regulate gene expression by decreasing the target mrna levels or inhibiting translation of protein encoding transcripts. global mirna abundance has been measured in skeletal muscle by microarray to evaluate the roles of mirnas in pig development and meat production. mirbase database has collected porcine mature mirnas. rna-seq also provides valuable insights into mirna transcriptome, especially into those mirnas insufficiently detected by microarray. until now, porcine mirna transcriptome has been investigated by next generation sequencer in intestine, pre- and postnatal piglet, developing brain and skeletal muscle. we herein performed a global transcriptome analysis on three types of tissues related to metabolism, meat production and fat deposition: liver, longissimus dorsi muscle and abdominal fat harvested from a full-sib f female pair with extreme phenotypes by rna-seq. the results allowed us to investigate large-scale alternative splicing events, identify extensive number of novel transcript units, determine gene boundaries at single nucleotide resolution and comprehensively survey porcine micrornas in the tested tissues. to our knowledge, this study presented the first systematical investigation on the complexity of porcine transcriptome with nucleotide resolution. moreover, this study identified many important candidate genes related to growth, meat quality, blood physiological and biochemical parameters by the integrated analysis of genome-wide association study and differential gene expression.  sequencing and mapping of the porcine transcriptome we sequenced cdna libraries from tissues of the full-sib female pair using high-seq at bgi-shenzhen, china. the data set was analyzed according to the bgi bioinformatics protocols for rna-seq. the sequence reads have been submitted to the ncbi gene expression omnibus under accession no. gse in total, we acquired,956-, paired-end reads of bp. the total read length was gigabases, representing about 8-fold of the porcine genome size. we technically replicated the rna-seq experiments in the samples with <r <. furthermore, the expression patterns of randomly selected transcripts between two individuals were validated by qrt-pcr with a relative coefficient of r =. the technical replicates and qrt-pcr confirmed the high reproducibility of rna-seq in this study. we aligned all short reads onto the whole reference genome. tolerances were set to allow at most two mismatches for bp reads in each alignment. about - % of reads were mapped to the pig reference genome, of which - % fell in annotated exons; - % located in introns; - % overlapped with exons, and the remaining - % were assigned to intergenic regions. total - % of reads had a uniquely genomic location, and - % of reads corresponded to reference genes with - % of uniquely matched reads. unmapped or multi-position matched reads were excluded from further analyses. identification of an extensive number of novel transcript units the uniquely mapped reads covered % of the annotated transcripts in ucsc pig genome database by at least one sequence read. a total of transcripts were expressed in all three tissues, and and transcripts were discovered exclusively in af, ld and liver, respectively. we quantified the gene expression level by counting the number of reads per kilobase per million mapped reads. about - % of the annotated transcripts showed expression with > rpkm. the obtained transcription fragments of more than % of the detected transcripts were > bp in length. the percentage of a gene covered by reads was defined as gene coverage. extensive read mapping revealed about - % of the detected transcripts with more than % in gene coverage. we detected an extensive number of novel transcript units by the procedures described in methods and zhang et al.. in total, we obtained,,,, and novel transcript units, respectively, in the six tested samples with a mean length of bp and a size range from to bp. many novel transcript units have > exons and the largest one is bp in length containing exons. the novel transcript units were often identified in clusters, indicating that closely spaced transcript units are likely to merge into longer transcripts with increasing sequencing depth. by scanning the genomic location of each novel transcribed unit, we clustered adjacent fragments into one transcriptionally active region. in this way, we identified - novel transcriptionally active regions in the six tested samples. a proportion of the regions was comprised of single exon. to investigate whether the identified novel transcript units were non-coding rnas, we aligned the sequences of the novel transcript units to non-coding rna precursor sequences in rfam database. we found that - % of novel transcript units were non-coding rna precursors. of these, % were in average trna precursors, % were mirna precursors and % were snorna precursors. alternative splicing events in pig transcriptome to investigate alternative splicing, we identified the sequence reads that were mapped to the regions of computationally determined theoretical splicing junctions. four known types of alternative splicing models including alternative 3' splicing site, alternative 5' splicing site, exon skipping and intron retention were considered in this study. the distribution of alternative splicing events is shown in additional file the alternative splicing events from tissues of one individual were pooled and the redundancy was removed to get a final set of alternative splicing events. we found that up to genes accounting for % of the known genes had undergone alternative splicing in these tissues, displaying alternative events in individual similarly, the number of alternative splicing genes and events in individual were and, respectively. figure 2a and 2b show an example of intron retention and alternative 5' splicing site. we found that about % of the alternative spliced genes underwent multiple alternative splicing events, illustrating the complexity of porcine transcriptome. alternative 3' splicing site is the most common type of alternative splicing events accounting for % of all alternative splicing events in pigs, while intron retention and exon-skipping only constituted % and % in individual or % and % in individual the average size of retained introns is bp, ranging from bp to bp. extension of annotated gene boundaries the extensions of 5' and 3' boundaries were determined by comparison of the gene models obtained by rna-seq with the existing gene annotations. in the six tested samples, a total of 587- genes were extended at the 5' end, of which more than % had an extension of at least bp in length. in comparison, - genes were extended at their 3' end, of which more than % were extended by at least bp. furthermore, 411- genes were extended at both ends. in individual total annotated genes were extended at the 5' end, of which % were observed in all three tissues; and annotated genes were extended at the 3' end, of which % had the extended 3' boundary in all three tissues. similarly, in individual the percentage of the extended 5' or 3' genes shared in all three tissues was % and %, respectively. comprehensive survey of porcine microrna by deep sequencing to get a comprehensive view of mirna transcriptome in pigs, we carried out deep sequencing of small rna using the tested samples mentioned above. the obtained mirna sequence reads have been submitted to the ncbi gene expression omnibus under accession no. gse we obtained a total of,703-, raw reads. after removing low quality reads and corrupted adapter sequences,,608-, clean reads were retained for further analyses. the majority of small rna was 20- nt for all libraries, with nt small rna being the most abundant, which is in agreement with the common nucleotide length of mirnas. we mapped - % of clean reads to the reference genome. chromosomes and harbored clean reads of >, tags. after further removal of trna, rrna, scrna, snrna, snorna, exon rna, intron rna and repeat regions, a total of,788-, mirna sequences were obtained. aligning mirna sequences to the porcine mature mirnas in mirbase database revealed - % of mature mirnas in each sample. we pooled all identified mature mirnas and found that only five mature mirnas in mirbase database were not detected in this study. the expression levels of the identified mirnas displayed a very large range, as reflected by the number of sequence reads, which varied from single counts for rare mirnas to several hundred thousand reads for the most abundant mirnas. as many as - % of the clean reads belong to the mirnas ranked top in expression levels in each sample. to identify potential novel mirna, we further analyzed the small rna tags that could not match known mirnas and were mapped to intergenic or intronic regions of the reference genome. the characteristic of hairpin structure of mirna precursor was used to predict novel mirna by exploring the secondary structure. in total, to potential novel mirnas each supported by at least five sequence reads, were identified in the tested samples. the potential novel mirnas of samples were pooled and the redundancy was removed to get a final set of unique porcine putative novel mirnas. the vast majorities of these mirnas were expressed at low levels. nevertheless, mirnas showed relatively high expression levels representing by more than sequence reads, and by more than sequence reads. to determine the genomic locations of these potential novel mirnas, their precursor sequences were blasted against the porcine reference genome sequence. as shown in figure 3c, % of the novel mirnas were assigned to the reference genome. most of the precursors were located on chromosome followed by chromosomes × and, and none of the novel mirnas were mapped to chromosome cross-tissue comparison of differential transcription events deep rna sequencing in three types of tissues allows us to investigate the tissue-differentially transcriptional events. of the identified transcripts, and transcripts were discovered exclusively in af, ld and liver, respectively. interestingly, as reported in humans, the majority of alternative splicing events showed clear tissue specificity, demonstrating the importance of alternative splicing in tissue specific programs of gene expression and its major roles in functional complexity. functional annotations of alternative splicing genes showed that they play important roles in their corresponding tissues. the porcine ryr gene which had the largest number of ld-specific alternative splicing events has been associated with malignant hyperthermia and has significant effects on pig meat quality and carcass leanness. the alb gene having the most number of liver-specific alternative splicing events plays important roles in transportation of fatty acids. sla- implicating in immune and type i diabetes had the most number of alternative splicing events in af that has been known as an important immune organ. other tissue-specifically alternative splicing genes also show major roles in their corresponding tissues. for instance, myhc-2a, ppap2c and apoh had tissue-specific alternative splicing in ld, af and liver, respectively, corresponding to their roles in normal muscle development and function, converting phosphatidic acid to diacylglycerol, transportation of fatty acids, and lipoprotein metabolism. the expression abundance of many mirnas also showed apparently tissue-differential patterns. the mature mirnas with high abundance in each tissue are conserved in mammals and likely related to important biological functions. mir- is the most abundant mirna in human liver and also had the most abundance of expressed read counts in porcine liver. mir- plays a crucial role in cholesterol, fatty acid and lipid metabolism. the mir- and mir- are key mediators in proper skeletal and cardiac muscle development and function, myogenesis during embryonic development and muscle cell differentiation. the two mirnas, respectively, had the most and third abundance of expression level in ld in this study. for those potential novel mirnas identified in this study, of putative novel mirnas were specifically identified in af, and in liver and ld, respectively. their biological functions need further investigation. differentially expressed transcripts and mirnas between two individuals overall, there were clear linear relationships in the gene expression levels between two individuals in all three tissues. the number of unique reads mapped to different genes ranged from to. the differentially expressed genes were selected based on the expression profiles and the following criteria: if the fold change in gene expression levels between and was more than or equal to two fold and if the false discovery rate value was less than. with this, we identified, and differentially expressed genes in af, liver and ld, respectively. of these, the expressions of, and genes were up-regulated in with respect to individual. to further investigate the biological relationships of differentially expressed genes with phenotypes, we performed the gene ontology analysis by querying each differentially expressed gene into the records of the go database. the results of go functional annotations are presented in figure 5d-5f. the main functional groups of differentially expressed genes in af are related to lipid metabolic process, alcohol metabolic process, lipid binding and protein homooligomerization. the functions of differentially expressed genes in liver are enriched in metabolic process, catalytic activity and oxidoreductase activity. and the differentially expressed genes in ld are mainly associated with metabolic process, cellular metabolic process, catalytic activity and oxidoreductase activity. differentially expressed mirnas between two individuals were identified by comparing the normalized expression data of the mature mirnas. in total, differentially expressed mirnas ≥ or fold-change ≤ -1; p-value < ) were identified in liver, and in ld and af, respectively. most of the differentially expressed mirnas had relatively low expression levels. interestingly, some differentially expressed mirnas are involved in the pathway relevant to development and diabetes. for instance, mir- differentially expressed in ld enables precisely specific the muscle cell types by sharpening cellular responses to hedgehog in zebrafish; a liver-differentially expressed mirna of mir-10b is predicted to regulate genes in pathways relevant to type diabetes. investigation of candidate genes for related phenotypes by integrated analysis of genome-wide association study and differential gene expression we selected total most differentially expressed genes including from liver and af, respectively, from ld for further functional annotations. only genes could be mapped to the pig reference genome and had functional annotations in mammals. of the genes, have the description of phenotypes in knocked-out mice. a genome-wide association study using the pig 60k snp chip has been performed in the current white duroc × erhualian f resource population. seven, and differentially expressed genes in liver, af and ld are located within mb around the snps that were most significantly associated with phenotypes in the gwas and the associated phenotypes were also observed in the corresponding gene-deficient mice. here, we present the identification results of interesting candidate genes igf cyp1a ckm and ces for the related phenotypes by integrating the analysis of gwas and differential gene expression in figure confirmation of igf as a causative gene for heart weight igf is the causative gene underlying the qtl for muscle growth, fat deposition and heart weight on ssc2p. a nucleotide substitution in intron of igf abrogates interaction with zbed resulting in threefold increase of igf messenger rna expression. in the current f population, a significant qtl for heart weight at day was located at chr2: mb with a corrected p-value of × 10- in gwas. this region overlapped with the igf locus. in the rna-seq analysis, individual with the heavier heart weight showed a fold increase in igf expression level compared with individual the result was perfectly consistent with the previous result, and is another example demonstrating the causality of this gene in the qtl. cyp1a as a strong candidate gene for hemoglobin two founder breeds of the f resource population are white duroc and chinese erhualian, which are divergent in hemoglobin concentration. blood parameters of hemoglobin and hematocrit at day were measured in the f population. we identified a significant qtl for both hemoglobin and hematocrit content at day which is located at chr7: mb in gwas. a differentially expressed gene of cyp1a in liver by rna-seq is located within this qtl region. the cyp1a knocked-out mice had the increased hemoglobin content. a quantitative trait transcript analysis in this f population also revealed that the cyp1a expression level was associated with hemoglobin at day. these results strongly support cyp1a as a candidate for the qtl effect on hemoglobin content. ckm as a candidate gene for pork ph value creatine kinase, muscle is a cytoplasmic enzyme involved in energy homeostasis. ckm-deficient mice show abnormal muscle physiology and an increased skeletal muscle glycogen level. excess glycogen content in pig skeletal muscle is associated with low meat ph value, leading to bad meat quality. a significant qtl for meat ph value at min, h and h after slaughter was detected around the region of chr2: mb with a p-value of × 10- in the f cross. ckm gene locates at the chr2: mb and was differentially expressed in this full-sib pair with different meat ph values. this indicates that ckm is an important candidate gene for meat ph value in pigs. ces as an important candidate gene for serum cholesterol and triglyceride level the liver is an important organ for lipid metabolism. gene expression in liver influences the circulating cholesterol level. ces that participates in fatty acyl and cholesterol ester metabolism was one of the differentially expressed genes in liver. our previous qtt study in liver showed that the expression level of ces is significantly associated with serum total cholesterol and triglyceride level in the f population. a significant qtl for both low density lipoprotein cholesterol and total cholesterol was identified at chr6: mb containing the ces gene in gwas. the concordant results of qtt, gwas and differential expression suggest that ces is strong candidate gene for qtl affecting serum total cholesterol on chromosome the supporting evidence was also from the ces knocked-out mice that showed decreased circulating cholesterol level. discussion in this study, we presented the systematical transcriptome profiling of pigs on three tissues related to metabolism, meat production and fat deposition using high throughput rna-seq technology. this efficient deep sequencing not only allows us to analyze novel transcribed regions and mirnas, but also improves gene annotations at single nucleotide resolution. integrated analysis of genome-wide association study and differential gene expression between two individuals revealed important candidate genes for related phenotypes. a total of,956-, reads were obtained from rna-seq in six samples with - % of reads uniquely mapped to the reference genome. this data set provided comprehensive starting resources for improving the gene annotations across the porcine genome. high r values between technical replicates showed high reproducibility of rna-seq in pigs. it should be mentioned that there are still - % of reads that cannot be matched to the reference genome. this could be caused by low sequence coverage of the reference genome, reference errors, sequencing errors and defined mapping criterions. a proportion of mapped reads with multiple positions was discarded, attributing to known duplicated genes or chromosome segments. a majority of the annotated transcripts at ucsc database was covered by sequence reads, showing the sensitivity of rna-seq in transcript discovery even for lowly expressed genes. furthermore, we identified large number of novel transcript units which improved the gene annotations of the porcine genome and transcriptome. alternative splicing is an important model of gene expression regulation and has not been generally accessible for microarray or sage methods in pigs. some genes showed all four types of alternative splicing models, revealing the complexity of alternative splicing in pigs. we found that alternative 3' splicing is the most common type of alternative splicing events in pigs. this is in contrast to the report in human and yeast where exon-skipping is the most prevalent mechanism, and it is also different from rice in which intron retention is the primary alternative splicing type. the percentage of alternative 5' or 3' splicing in total alternative splicing events in this study were higher than that reported in lim et al. where % of alternative 5' splicing and 3% of alternative 3' splicing were observed. we discovered that more than % of the detected genes were alternatively spliced. this number is much lower than the reported % in human and % in rice. three types of alternative splicing models including alternative first exon, alternative last exon and mutually exclusive exon were excluded from analyses because of currently unperfected algorithms. more alternative spliced genes would be discovered if these alternative splicing models were considered. as an important regulator of gene expression, mirna regulates the gene expression through decreasing the target mrna levels or repressing the translation. in this study, the high identification rate illustrates that the six small rna libraries from the tested tissues almost encompass the entire repertoire of known mirnas. consistent with the results in li et al., porcine mature mirnas had a broad range of expression levels. the highly expressed mirnas are known to have important regulatory functions in corresponding tissues. the mir-1/ showed high abundance of expression levels in muscle. interestingly, muscular hypertrophy in texel sheep has been shown to be caused by a mutation that creates an illegitimate binding site for mir-1/ in the 3' utr of the myostatin gene, leading to efficient translational inhibition of the myostatin gene and an increase in muscularity. the mir- had the most abundant expression level in liver. it has the diversity of its roles in liver, e.g. metabolism, hepatocarcinogenesis. just as in humans, most of the potential novel mirnas discovered in this study are expressed at low levels, explaining why they were not discovered in previous efforts and showing the advantage of rna-seq in transcriptome analysis. furthermore, we found that most of the novel mirnas were tissue-specifically identified. we cautioned that some could be artifact due to the low expression level, the insufficient depth of sequencing and limited tested tissues. mapping genetic factors that underlie quantitative traits in farm animals has been a challenging task. the recent wave of genome-wide association studies in human showed that a majority of snps associated with disease traits locate in regulated regions. integration of gene expression with genotype and phenotype data to elucidate the network of molecular interactions that underlie complex traits can facilitate the identification of variants that contribute to phenotypes. in this study, high-throughput sequencing allowed us to digitally discover an extensive number of differentially expressed genes in the full-sib pair with different phenotypes. gene ontology analysis indicated that the differentially expressed genes had enrichment on functions related to metabolic process, catalytic activity and lipid binding. the results suggest that these differentially expressed genes are likely related to the phenotypes of growth, metabolism or fat deposition. by integrating gwas, differentially expressed genes and altered phenotypes in knocked-out mice, we found that many differentially expressed genes are the important candidate genes related to the phenotypes of serum cholesterol, growth traits, hemoglobin at day fatty acid level and muscle ph value. for example, activating transcription factor gene differentially expressed in muscle and located in the qtl region for body weight at day is related to decreased body weight in atf knocked-out mice; fatty acid binding protein located in the qtl region for fat deposition on ssc and related to increased white adipose tissue amount in fabp knocked-out mice was differentially expressed in abdominal fat in this study. the findings provide important clues for further dissecting of the responsible genes and variants.  this study provides a global view of the complexity of the pig mrna and mirna transcriptome, gives an extensive new knowledge about alternative splicing, novel transcript units, gene boundaries and novel mirnas in pigs. the cross-tissue comparison identified lots of tissue-differential transcription events. integrated analysis of gwas and differential gene expression allowed us to detect important candidate genes related to growth, meat quality, serum lipids and fatness. the findings significantly enhance the current genome annotation of pigs and improve our understanding of complex traits. 
here we used massive parallel high-throughput rna sequencing to generate a high-resolution map of the porcine mrna and mirna transcriptome in liver, longissimus dorsi and abdominal fat from two full-sib f hybrid pigs with segregated phenotypes on growth, blood physiological and biochemical parameters, and fat deposition. elucidation of the pig transcriptome is essential for interpreting functional elements of the genome and understanding the genetic architecture of complex traits such as fat deposition, metabolism and growth. we obtained,418-, uniquely mapped reads that covered % of the current annotated transcripts and identified - novel transcript fragments, which constituted - novel transcriptional active regions in six tested samples. cross-tissue comparison revealed that many transcriptional events are tissue-differential and related to important biological functions in their corresponding tissues. this study provides a global view of the complexity of the pig transcriptome, and gives an extensive new knowledge about alternative splicing, gene boundaries and mirnas in pigs. integrated analysis of genome-wide association study and differential gene expression revealed interesting candidate genes for complex traits, such as igf cyp1a ckm and ces for heart weight, hemoglobin, pork ph value and serum cholesterol, respectively. integrated analysis of genome wide association study and differential gene expression allows us to find important candidate genes for porcine complex traits. we found that about % of the annotated genes showed alternative splicing patterns, and alternative 3' splicing is the most common type of alternative splicing events in pigs. we also detected a total of potential novel mirnas, most of which were tissue-specifically identified.
 elucidation of the pig transcriptome is essential for interpreting functional elements of the genome and understanding the genetic architecture of complex traits such as fat deposition, metabolism and growth.  here we used massive parallel high-throughput rna sequencing to generate a high-resolution map of the porcine mrna and mirna transcriptome in liver, longissimus dorsi and abdominal fat from two full-sib f hybrid pigs with segregated phenotypes on growth, blood physiological and biochemical parameters, and fat deposition. we obtained,418-, uniquely mapped reads that covered % of the current annotated transcripts and identified - novel transcript fragments, which constituted - novel transcriptional active regions in six tested samples. we found that about % of the annotated genes showed alternative splicing patterns, and alternative 3' splicing is the most common type of alternative splicing events in pigs. cross-tissue comparison revealed that many transcriptional events are tissue-differential and related to important biological functions in their corresponding tissues. we also detected a total of potential novel mirnas, most of which were tissue-specifically identified. integrated analysis of genome-wide association study and differential gene expression revealed interesting candidate genes for complex traits, such as igf cyp1a ckm and ces for heart weight, hemoglobin, pork ph value and serum cholesterol, respectively.  this study provides a global view of the complexity of the pig transcriptome, and gives an extensive new knowledge about alternative splicing, gene boundaries and mirnas in pigs. integrated analysis of genome wide association study and differential gene expression allows us to find important candidate genes for porcine complex traits. novel transcript unitalternative splicinggene boundarymirnadifferentially expressed genecomplex traitpig 