Transcriptome Analysis Reveals that Vitamin A Metabolism in the Liver Affects Feed Efficiency in Pigs

Feed efficiency (FE) is essential for pig production. In this study, 300 significantly differentially expressed (DE) transcripts, including 232 annotated genes, 28 cis-natural antisense transcripts (cis-NATs), and 40 long noncoding RNAs (lncRNAs), were identified between the liver of Yorkshire pigs with extremely high and low FE. Among these transcripts, 25 DE lncRNAs were significantly correlated with 125 DE annotated genes at a transcriptional level. These DE genes were enriched primarily in vitamin A (VA), fatty acid, and steroid hormone metabolism. VA metabolism is regulated by energy status, and active derivatives of VA metabolism can regulate fatty acid and steroid hormones metabolism. The key genes of VA metabolism (CYP1A1, ALDH1A2, and RDH16), fatty acid biosynthesis (FASN, SCD, CYP2J2, and ANKRD23), and steroid hormone metabolism (CYP1A1, HSD17B2, and UGT2B4) were significantly upregulated in the liver of high-FE pigs. Previous study with the same samples indicated that the mitochondrial function and energy expenditure were reduced in the muscle tissue of high-FE pigs. In conclusion, VA metabolism in liver tissues plays important roles in the regulation of FE in pigs by affecting energy metabolism, which may mediate fatty acid biosynthesis and steroid hormone metabolism. Furthermore, our results identified novel transcripts, such as cis-NATs and lncRNAs, which are also involved in the regulation of FE in pigs.

. In young ferrets, chronic supplementation with b-carotene, which is a natural material for VA synthesis, increase their body weight and subcutaneous fat mass (Murano et al. 2005). The short-term atRA treatment in adulthood is associated with a decrease in adiposity in ferrets (Sánchez et al. 2009). The effect of VA and its metabolites on energy metabolism may be dependent on developmental stage.
Molecular mechanism studies have indicated that VA, especially its metabolites, can regulate several genes related to energy metabolism. atRA and 9-cis RA can bind to retinoic acid receptors (RARs) with high affinity in vitro. The 9-cis RA specifically binds to retinoid X receptors (RXRs) in vitro (Blomhoff and Blomhoff 2006). RAR:RXR heterodimers regulate typical RA target genes by binding to RA response elements, which include phosphoenolpyruvate carboxykinase (PEPCK) (Cadoudal et al. 2008) and stearoyl-CoA desaturase 1 (SCD1) (Samuel et al. 2001;Zolfaghari and Ross 2003). Moreover, RXR can form heterodimers with liver X receptor (LXR), which can positively regulate key lipogenic genes, including SREBP-1C, SCD, and FASN, in mouse liver cells and human hepatoma cells (Mukherjee et al. 1997;Roder et al. 2007). Therefore, VA and its metabolites play important roles in energy metabolism.
Long noncoding RNAs (lncRNAs) and natural antisense transcripts (NATs) have emerged as essential components of regulatory factors in several physiological processes. lncRNAs are involved in adipogenesis, hepatic lipid metabolism, energy balance, and so on (Zhao and Lin 2015). The knockdown of the liver-enriched lncRNA, lncLSTR, can reduce plasma triglyceride levels in a hyperlipidemic mouse model . sno-lncRNAs from the Prader-Willi syndrome locus modulate energy balance in mice (Yin et al. 2012;Powell et al. 2013). In humans, the lncRNA HULC modulates lipid metabolism in hepatocellular carcinoma by activating the acyl-CoA synthetase subunit ACSL1 (Cui et al. 2015). In previous studies, NATs, which include trans-and cis-NATs, have been identified in humans, mice, and pigs, etc. (Katayama et al. 2005;Zhang et al. 2006;Zhao et al. 2016). In eukaryotic, NATs-mediated gene expression regulation mechanisms include chromatin remodeling, transcriptional interference, RNA masking, and doublestranded RNA-dependent mechanisms (Lapidot and Pilpel 2006;Faghihi and Wahlestedt 2009). The antisense transcript of apolipoprotein A1 (APOA1), referred to as ApoA1-AS, is an lncRNA that negatively regulates APOA1 expression (Halley et al. 2014). However, the roles of lncRNAs and NATs in the FE of pigs remain largely unknown.
On the basis of high-throughput RNA sequencing, we systematically analyzed differentially expressed (DE) annotated genes, lncRNAs, and cis-NATs in liver tissues between low-and high-FCR pigs. Gene ontology (GO) and pathway analysis revealed that the VA metabolism pathway in liver tissues was related to the variation of FE in pigs.

MATERIALS AND METHODS
Animals, tissues, and RNA extraction In this study, the feed intake of 236 purebred castrated boars from a Yorkshire pig population was detected using an ACEMA 64 automated individual feeding system at the Agricultural Ministry Breeding Swine Quality Supervision Inspecting and Testing Center (Wuhan, China) (Jing et al. 2015). The FCR of each individual was analyzed in the R environment using the following formula:

FCR ¼
Feed intake Average daily gain : Furthermore, the FCR value between the three lowest pigs and the three highest pigs was significantly different (P , 0.05) (Jing et al. 2015). Those six castrated individuals were selected for RNA sequencing; the selected animals were not sibships. For tissue samples, pigs were slaughtered at 90 kg according to a standard procedure approved by guidelines from Regulation of the Standing Committee of Hubei People's Congress (Hubei Province, P. R. n  Figure 1 Annotation of the uniquely mapped reads of RNA-seq in liver tissue of pigs. (A) Distribution of uniquely mapped reads in the pig genome. The percentages in this pie chart represent the mean of all six RNA-seq data. On average, over 50% uniquely mapped reads were located in the CDS region. The untranscriptional regions were defined as genome region of annotated genes but without coverage of annotated transcripts. (B) Reads distributed on the sense and antisense strands of annotated genes. Approximately 2% reads were mapped on antisense strands of annotated genes.
China). Liver tissues were sampled and snap-frozen in liquid nitrogen within half an hour after slaughter before storage at -80°. Briefly, total RNA was isolated from frozen liver tissues with the TRIzol reagent (Invitrogen). All experimental protocols were approved by the Ethics Committee of Huazhong Agricultural University (HZAUMU2013-0005).

Library construction and sequencing
For each liver sample, equal quantities of RNA were sent to Genergy Biotechnology (Shanghai, China) for library construction. The RNA-seq library of each sample was prepared by with the TruSeq Stranded Total RNA Sample Preparation kit (Illumina). The Second Strand Marking Mix was also used during second strand cDNA synthesis to replace dTTP with dUTP to ensure strand specificity. After quality control, sequencing was performed with Illumina HiSeq2000.

RNA sequencing analysis
TopHat (version 2.0.9) (Trapnell et al. 2009) software was used to align reads to the pig reference genome (NCBI Sscrofa10.2). Up to two mismatches were allowed in reads mapping. Multiple-mapped reads were discarded; uniquely mapped reads were compared with the gff3 file of the Sscrofa10.2 genome with in-house Perl scripts. The reference genome sequence and its gff3 file were downloaded from the NCBI genome database (ftp://ftp.ncbi.nlm.nih.gov/genomes/Sus_scrofa/).
Furthermore, mapped reads in the intergenic region were also compared to annotated lncRNAs (Zhou et al. 2014).

Genome-wide identification of cis-NATs
Mapped reads were pooled into one SAM file, which was used to identify novel transcripts by Cufflinks (Trapnell et al. 2012) with the option -library-type fr-firststrand. The transcripts identified by Cufflinks were compared with the gff3 file of the Sscrofa10.2 genome. Cis-NATs were identified on the basis of the following criteria: (i) located on the antisense strand of annotated genes, and (ii) novel transcripts but not belonging to the annotated genes. Cis-NATs were extracted from novel transcripts by in-house Perl scripts. The identified cis-NATs were used for further analysis.

Differential expression analysis
The count of reads located in the exon regions of each annotated gene, cis-NAT, and lncRNA was calculated using in-house Perl scripts. Genes with CPM (at least one count per million) . 1 in at least four samples were kept for further analysis. Both normalization of expression profiling for all expressed transcripts, and identification of DE transcripts were performed by edgeR (Robinson et al. 2010) in the R environment. Transcripts were determined as significantly DE with FDR , 0.05, and with an upregulated or downregulated fold-change (FC) of $ 2.5 between low-FCR and high-FCR pigs.

Correlation analysis
Pearson correlation analysis was performed to identify the correlatively expressed SA (sense-antisense) pairs and DE annotated gene-lncRNA pairs in the R environment. The cis-NATs, annotated genes, and lncRNAs that were detected at least in four samples were chosen for correlation analysis. The criteria for significantly correlated SA pairs was P , 0.05 and |R| . 0.8, whereas those for lncRNA-annotated gene pairs was P , 0.01 and |R| . 0.95.

q-PCR validation of DE genes
Total RNA was reverse-transcribed into cDNA using a PrimeScript RT reagent kit (Takara Bio Inc.). Oligonucleotide primers for six DE genes were designed with oligo7 software; the primer sequences are available in Supplemental Material, Table S1. Relative expression levels of these genes in liver were quantified by qPCR. The reactions were conducted on a BIO-RAD CFX384 Real-Time System with SYBR Green PCR Master Mix (Bio-Rad) as described in the manufacturer's instruction manual. The 10 ml reaction mixture consisted of 1 ml cDNA, with 5 ml 2· SYBR Green PCR Master Mixture, 0.2 ml each of the forward and reverse primers, and 3.6 ml of RNase-free water. Samples were preincubated at 95°for 3 min, followed by 40 PCR amplification cycles (denaturation: 95°for 20 sec, annealing: 60°for 20 sec, elongation: 72°for 15 sec). A dissociation curve was generated at the end of the last cycle by collecting the fluorescence data from 60°to 95°. Relative gene expression levels were normalized to the RPL32 gene, which has stable expression in liver tissues (Ostrowska et al. 2014), by the 2 2DDCt method. Student's t-test was used to analyze the expression difference between the low-FCR and high-FCR groups.
Gene ontology enrichment and pathway analysis GO enrichment analysis was performed with DAVID Bioinformatics Resources 6.7 (Da Wei Huang and Lempicki 2008). DE genes were sorted by their involvement in significantly enriched biological process GO terms. The pathway that involved the significant GO-terms-enriched genes was structured by referring to the KEGG pathway database and published articles. cis-NATs and DE lncRNAs that were correlated with the expression of genes involved in structured pathways, were also investigated and displayed. Visualization of pathway analysis was implemented in Cytoscape (Saito et al. 2012). The log 2 FC of annotated genes, cis-NATs, and lncRNAs was calculated as log 2 FC ¼ log 2 ðLow=HighÞ; where Low and High represented expression profile of the low-FCR and high-FCR groups, respectively.

Data availability
The raw data of RNA-seq were submitted to NCBI Sequence Read Archive database (SRA) under series SRP076030 which will be release in June 2017. DE transcripts list can be found in File S2.

RESULTS
RNA sequencing data mapping and annotation RNA of liver tissues of three high-FCR (H) and three low-FCR (L) pigs was isolated for RNA-seq. The RNA-seq raw data have been submitted to the NCBI Sequence Read Archive (SRA) under series SRP076030. After removing adaptors and filtering low quality reads, 7.1-10.0 million clean reads were yielded (Table 1). In these clean reads, $87% (6.1-9.6 million) were mapped on the Sscrofa10.2 genome, and over 90% of them were unique mappings ( Table 1). The majority of uniquely mapped reads (90.72%) were located in the annotated gene region, whereas . 50% were located in the CDS region ( Figure 1A). However, nearly 20% uniquely mapped reads were located in the intron region (9.28%), or other nontranscriptional regions (12.24%; Figure 1A). We also identified 1.8-2.5% uniquely mapped reads that were distributed on the antisense strand of annotated genes ( Figure 1B).

Identification of cis-NATs
Based on the strand-specific RNA-seq reads, cis-NATs could be identified directly. Reads that were mapped on the antisense strand of the annotated gene were identified and defined as antisense reads. A total of 0.78 million antisense reads were identified from the uniquely mapped reads. Finally, 1769 cis-NATs were identified (File S1).
To understand the structural characteristics of the cis-NATs, we further investigated their length distribution and exon number. Statistical results showed that the majority (78.80%) of cis-NATs were shorter than 2 kb (Figure 2A). The statistical results of reported lncRNAs showed that $50% of the lncRNAs were shorter than 5 kb ( Figure 2B). Both cis-NATs and lncRNAs differed with the detected annotated genes. For the annotated genes, . 80% were . 5 kb ( Figure 2C). Moreover, the overwhelming majority of cis-NATs (91.22%) and lncRNAs (93.84%) contained fewer than three exons ( Figure 2D). Besides, 78.83% of cis-NATs contained only one exon ( Figure 2D). Most of the detected annotated genes contained more than three exons ( Figure 2D). Therefore, the structural characteristics of cis-NATs were similar to those of lncRNAs, and both of them were different from the detected annotated genes.

DE liver transcripts between high and low FCR pigs
To compare the transcriptome alteration between high-and low-FCR pigs, the DE annotated genes, cis-NATs, and lncRNAs were determined using the edgeR package. A total of 300 significantly DE transcripts were identified, which included 232 annotated genes, 28 cis-NATs, and 40 lncRNAs, respectively (up or downregulated FC $ 2.5 and FDR , 0.05; File S2). Among these 232 annotated genes, 126 genes were upregulated, and 106 genes were downregulated in the low-FCR pigs (Figure 3 and File S2). Moreover, the majority of the DE cis-NATs (35 of 40) and lncRNAs (20 of 28) were downregulated in the low-FCR pigs (Figure 3 and File S2). The top 10 DE annotated genes cis-NATs and lncRNAs are listed in Table 2. qPCR was performed to validate the DE transcripts that were identified by RNA-seq. A total of 10 low-FCR individuals and 10 high-FCR individuals, which contained the individuals for RNAseq, were chosen for qPCR analysis. The difference between these lowand high-FCR pigs was significant. Here, six DE transcripts were chosen for qPCR analysis. Three DE genes (HSD17B2, CYP1A1, and FASN) were selected from the GO enrichment analysis, and another three DE transcripts (SLC27A6, linc-ssct5281, and linc-ssct3489) were selected randomly from DE analysis results (File S2). The qPCR results shows that all six selected genes were validated as significantly DE genes in n low-FCR vs. high-FCR. Moreover, the correlation coefficient of the log 2 FC values between RNA-seq and qPCR was 0.93 (P , 0.05).

Identification of correlated expression pairs
To understand function of cis-NATs, correlation analysis was performed between the expression of cis-NATs and their corresponding sense annotated genes. A total of 91 cis-NATs were significantly (P , 0.05 and |R| . 0.8) correlated with the expression of their sense annotated genes ( Figure 4A and File S3). Among these pairs, 61 (67.4%) were positively correlated ( Figure 4A and File S3). We also analyzed the correlation between the expression of DE lncRNAs and their annotated genes. In total, 337 significant correlation pairs were identified (P , 0.01 and |R| . 0.95), which included 26 DE lncRNAs and 126 DE annotated genes ( Figure 4B and File S3). Furthermore, the majority (19 in 25) of the lncRNAs were correlated with at least two annotated genes ( Figure 4B and File S3). Moreover, 333 of 337 correlated pairs (98.81%) were positively correlated.

Function annotation
Function enrichment analysis was performed to detect significant biological processes GO terms of the DE annotated genes and cis-NATs according to the DAVID Bioinformatics Resources 6.7. A total of 39 significantly (EASE Score , 0.1) enriched GO terms were identified ( Figure 5A and Table S2). The GO enrichment results showed the GO terms could be clustered based on the same DE annotated genes or similar biological functions ( Figure 5, A and B). Based on that, the resultant 39 GO terms were clustered into four categories: CYP1A1or ALDH1A2-related, fatty acid metabolic processes, immunerelated, and other GO terms ( Figure 5, A and C, left pie plot). Furthermore, 11 of 16 CYP1A1-or ALDH1A2-related GO terms contained both genes. Distribution investigation of annotated genes, among significantly GO terms, showed that CYP1A1 and ALDH1A2 were the top two most enriched genes ( Figure 5B). Each was involved at least 12 GO terms ( Figure 5B). Moreover, five of the 16 CYP1A1-or ALDH1A2-related GO terms were related to vitamin metabolism ( Figure 5C).

Pathway analysis of DE genes and lncRNAs
Based on functional enrichment analysis, CYP1A1 and ALDH1A2 were the two most important genes, mainly participating in VA-and hormone-related processes. Furthermore, we drew the VA and steroid hormone metabolism pathways according to the KEGG pathway database and previous studies. The results showed that CYP1A1, ALDH1A2, and RDH16 were both involved in VA anabolism and catabolism (Figure 6, right). Moreover, fatty acid and steroid hormone metabolic processes were located downstream of the VA metabolism pathway. These processes were regulated by RA, which is a metabolite of VA. The varied expression of the key genes and lncRNAs involved in these pathways were labeled by different colors. In Figure 6, genes in red and pink were upregulated in low-FCR pigs, whereas genes in green and blue were downregulated. The CYP1A1, ALDH1A2, RDH16, HSD17B2, UGT2B4, FASN, SCD, CYP2J2, and ANKRD23 genes were upregulated, in low-FCR pigs; the DE lncRNAs linc-ssct5500, linc-ssct5436, and linc-ssct0074 were also positively correlated with SCD, CYP2J2, and HSD17B2, respectively ( Figure 6). Other VA metabolisminvolved genes (ADH5, CYP26A1, AOX1, and RXRB) were also slightly upregulated in low-FCR pigs.

DISCUSSION
The improvement of FE is one of the most efficient ways to improve the benefits of pig production. Thus, investigation of the mechanisms of FE is very important for pig breeding. In this study, we systematically  analyzed the transcription profiles of liver tissue in high-and low-FCR pigs. We found that the VA metabolism pathway in liver tissues was important for FE in pigs. The upregulation of genes related to VA metabolism in the liver tissue was positively correlated with FE in pigs.
The liver is one of the main tissues for VA storage in mammals. In this study, the key genes of VA metabolism, namely, CYP1A1, ALDH1A2, and RDH16, were all significantly upregulated in the liver of high-FE (low-FCR) pigs. This result indicated that the pathways of VA metabolism differed between high-and low-FE pigs. A recent study indicated that RA biosynthesis is regulated by energy status at the rate-limiting step catalyzed by retinol dehydrogenases (RDH) (Obrochta et al. 2015). The retinol dehydrogenase family members RDH1, RDH10, and RDH16 had significantly decreased transcriptional levels after refeeding, oral gavage with glucose, or injection with insulin (Obrochta et al. 2015). These results implied a negative relationship between RA biosynthesis and energy levels. Moreover, previous studies indicated that mitochondrial reactive oxygen species, and mitochondrial uncoupling reactions involving genes of high-FE pigs, were lower than those of low-FE pigs (Grubbs et al. 2013a;Jing et al. 2015). However, the antioxidant defenses of high-FE pigs were higher than those of low-FE pigs (Grubbs et al. 2013b). The oxidation and energy loss of high-FE pigs was lower than that of low-FE pigs. Thus, the upregulation of genes involving in VA metabolism in high-FE pigs was possibly caused by the relatively lower oxidation and energy loss in this study. Besides, the key genes of fatty acid biosynthesis SCD (Samuel et al. 2001;Zolfaghari and Ross 2003;Zhang et al. 2012;Weiss et al. 2014) and FASN (Zhou et al. 2010;Zhang et al. 2012) were induced by RA at the transcriptional level and significantly upregulated in high-FE pigs. SCD encodes stearoyl-CoA desaturase, which converts saturated fatty acids into monounsaturated fatty acids, and influences fatty acid partitioning in the liver by promoting fatty acid synthesis but decreasing oxidation (Flowers and Ntambi 2009). Furthermore, other genes related to fatty acid biosynthesis, namely, CYP2J2 and ANKRD23 were also upregulated. These results indicated that fatty acid biosynthesis in the liver was increased in high-FE pigs. Therefore, the VA pathway could affect FE by involving the energy metabolism of pigs and promoting fatty acid biosynthesis. Moreover, the moderate activity of the VA metabolic pathway in liver tissue may benefit FE in pigs.
Previous studies also indicated that RA could induce and maintain testosterone production (Tucci et al. 2008). Testosterone is a steroid hormone in vertebrates and belongs to the androgen group. This hormone plays a key role in male reproductive tissue development, as well as in increasing muscle mass, muscle strength, bone mass, and bone mineral density (van den Beld et al. 2000;Sinha-Hikim et al. 2002). In this study, genes involved in testosterone metabolism, namely, CYP1A1, HSD17B2, and UGT2B4, were significantly upregulated in low-FCR pigs. Moreover, HSD17B2 is positively regulated by RA (Ito et al. 2001). HSD3B1 and AKR1C2 are involved in testosterone metabolism. and were also slightly upregulated in high-FE pigs. These results suggested that the metabolism of testosterone was increased in high-FE pigs; VA may also affect FE by affecting steroid hormone metabolism.
In addition, 40 lncRNAs were differentially expressed between highand low-FCR pigs. lncRNAs are involved in several bioprocesses, including energy metabolism (Yin et al. 2012;Powell et al. 2013;Cui et al. 2015;Li et al. 2015) and growth (Dey et al. 2014;Mueller et al. 2015). The DE lncRNAs linc-ssct5500, linc-ssct5436, and linc-ssct0074 were upregulated in the livers of low-FCR pigs, and positively correlated with SCD, CYP2J2, and HSD17B2, respectively. Previous studies indicated that SCD and CYP2J2 (Wang et al. 2004) are involved in fatty acid metabolism, whereas HSD17B2 is related to testosterone metabolism Figure 6 Potential pathway of the annotated DE genes and lncRNAs in the liver tissues of high and low FCR pigs. Red and pink indicates upregulation in low FCR pigs (red, log2FC . 0.5; pink, 0 , log2FC # 0.5); green and blue indicates downregulation in low FCR pigs (green, log 2 FC , 20.5; blue, 20.5 # log 2 FC , 0). (Labrie et al. 1997). Therefore, these lncRNAs may also affect the FE of pigs by altering the metabolism of fatty acids and steroid hormones.
We also identified 1769 cis-NATs in pig liver tissue. Among them, 28 were significantly different between low and high FE pigs. Correlation analysis showed 61 (67.4%) SA pairs from 91 significantly correlated SA pairs were positively correlated. The coexpression (positively) between SA pairs was reported in humans (Chen et al. 2005), Drosophila melanogaster (Okamura et al. 2008), and pigs (Zhao et al. 2016). These results indicate that cis-NATs are involved in the regulation of FE of pigs. However, the mechanisms of the regulatory roles of cis-NATs in FE of pigs remained largely unknown. Furthermore, some immunerelated signaling pathways were also identified to relate to FE. Therefore, the FE of pigs may also be regulated by ncRNAs or other pathways.
In summary, annotated DE genes, cis-NATs, and lncRNAs in the liver tissues of high-and low-FE pigs were analyzed comparatively. Our results revealed that different expression was enriched mainly in VA, fatty acid, and steroid hormone metabolism. VA metabolism in liver tissues involved regulation of FE in pigs by mediating fatty acid biosynthesis and steroid hormone metabolism. CYP1A1, ALDH1A2, and RDH16 are the key upstream genes of VA metabolism, and possible candidate genes for FE in pigs. Some lncRNAs and cis-NATs were also related to the FE of pigs. Author contributions: X.L. and S.Z. conceived and supervised the project. Y.Z., F.L., Y.L., A.L., and L.J. undertook the analysis. Y.Z. and YH wrote the paper. Y.H., C.Z., and Y.M. performed the tissue collection and PCR experiments. All authors read and approved the final manuscript. The authors declare that they have no competing interests.