Integrated mRNA and miRNA profile expression in livers of Jinhua and Landrace pigs

Objective To explore the molecular mechanisms of fat metabolism and deposition in pigs, an experiment was conducted to identify hepatic mRNAs and miRNAs expression and determine the potential interaction of them in two phenotypically extreme pig breeds. Methods mRNA and miRNA profiling of liver from 70-day Jinhua (JH) and Landrace (LD) pigs were performed using RNA sequencing. Blood samples were taken to detect results of serum biochemistry. Bioinformatics analysis were applied to construct differentially expressed miRNA-mRNA network. Results Serum total triiodothyronine and total thyroxine were significantly lower in Jinhua pigs, but the content of serum total cholesterol (TCH) and low-density lipoprotein cholesterol were strikingly higher. A total of 467 differentially expressed genes (DEGs) and 35 differentially expressed miRNAs (DE miRNAs) were identified between JH and LD groups. Gene ontology analysis suggested that DEGs were involved in oxidation-reduction, lipid biosynthetic and lipid metabolism process. Interaction network of DEGs and DE miRNAs were constructed, according to target prediction results. Conclusion We generated transcriptome and miRNAome profiles of liver from JH and LD pig breeds which represent distinguishing phenotypes of growth and metabolism. The potential miRNA-mRNA interaction networks may provide a comprehensive understanding in the mechanism of lipid metabolism. These results serve as a basis for further investigation on biological functions of miRNAs in the porcine liver.


INTRODUCTION
The pig is an important livestock animal and a major meat producer in China. At the same time, pigs are one of the ideal animal models for human medical research, because of high similarity in structure of physiological, anatomical, physiological and biochemical index, food structure and drug metabolism with human [1]. The Jinhua (JH) pig is a famous Chinese indigenous pig breed, being known for its slow growth and high fat deposition [2] and is widely used for ham production. On the contrary, the Landrace (LD) pig represents the fast growing lean type selected for high carcass yield [3]. Therefore, divergent performance makes the two pig breeds ideal comparisons to investigate genetic diversity and molecular mechanisms in the metabolic processes in the liver.
Liver is a critical organ in the metabolism of carbohydrates, fatty acids, amino acids and so on. It undertakes most of the synthesis, decomposition, transformation, excretion and other metabolic processes. MicroRNAs (miRNAs) are important regulators in gene expression which play key roles in biological process such as cell proliferation and differentiation [4], and pathogenesis and disease prevention [5]. The regulation of miR-122 has been as-sociated with lipid metabolism and liver diseases [6,7]. Mir-224-5p directly targets to early growth response 2 (EGR2) and acyl-CoA synthetase long-chain family member (ACSL4), which is a novel negative regulator of adipocyte differentiation through post-transcription regulation of EGR2 during early adipogenesis and regulates fatty acid metabolism at terminal differentiation through ACSL4 [8]. MiR-130 and miR-27 suppress adipogene sis by inhibiting peroxisome proliferator-activated receptor γ (PPARγ) expression [9,10]. Although the regulatory mechanisms of metabolism in other breeds have been studied extensively, the differences between JH and LD pigs remain unreported.
The intention of this study was to characterize and compared the transcriptome profiles and of the liver in two porcine breeds with distinct phenotypes using RNA-seq technology. This study aimed to construct mRNA-miRNA regulatory network related in breed specific growth and fat deposition.

Ethics statement
The experiment was approved by the Animal Care and Use Committee of Zhejiang University, Zhejiang, China (ZJU 20160346). All procedures complied with the guidelines and the regulations for the Administration of Affairs Concerning Experimental Animals.

Sample collection and RNA extraction
The liver tissue samples were collected from three uncastrated male JH pigs and three uncastrated male LD pigs respectively at 70-day. The liver samples were immediately frozen in liquid nitrogen and stored at -80°C until RNA extraction. Total RNA was isolated from liver tissue samples using TRIzolTM reagent (Invitrogen, Carlsbad, CA, USA), according to the manufacturer's protocol. RNA quantity and purity were analyzed with the ND-1000 Nanodrop (Thermo Fisher, Wilmington, DE, USA). The criteria used to select the RNA for following analysis were A260/A280≥1.8, A260/A230≥2.0 and RNA integrity number >7.0.

Serum biochemical parameters assays
The blood samples were collected from the precaval vein of each animal, kept at 37°C for 2 h, and then centrifuged at 4°C for 10 min (3,000×g). Serum glucose, triglyceride (TG), total triiodothyronine (TT3), total thyroxine (TT4), serum total cholesterol (TCH), high density lipoprotein cholesterol, low density lipoprotein cholesterol (LDLC) were determined by Affiliated Hospital of Hangzhou Normal University.

Small RNA sequencing and data analysis
Approximately 1 μg RNA of total RNA were used to prepare a small RNA library according to instructions of NEBNext Multiplex Small RNA Library Prep Set for Illumina (Illumina, San Diego, CA, USA). The purified library products were evaluated using the Agilent 2200 TapeStation and diluted to 10 pM for cluster generation in situ on the HiSeq2500 singleend flow cell followed by sequencing (1×50 bp) on HiSeq 2500. After removing adapter sequences, reads containing poly-N and low quality reads, all clean reads were mapped to Rfam12.1 (rfam.xfam.org), pirnabank (pirnabank.ibab.ac.in), miRBase 21 (www.miRBase.org) and Mireap to annotate rRNA, tRNA, snRNA, snoRNA, piRNA, mature miRNAs and porcine novel miRNA. The miRNA expression levels were calculated by reads per million (RPM) values (RPM = [number of reads mapping to miRNA/number of reads in Clean data]×10 6 ). Differential expression analysis between two breeds was determined by edger v3.10.0 and |log2 (fold change)| >1 and p<0.05 were set as the threshold for significance. Small RNA-seq data have been deposited in the gene expression omnibus (GEO) database and are available through the series accession numbers GSE124484.

RNA sequencing and data analysis
Approximately 1 μg RNA per sample were used to construct the complementary (cDNA) library with NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA) according to the manufacturer's instructions. After cluster generation using TruSeq Rapid SR Cluster Kit V2 (Illumina, USA), six libraries were sequenced on Illumina Hiseq 2500 platform and 50 bp single reads were generated. After filtering out adaptor sequences and removing low quality reads from raw data, the clean reads were aligned to the reference genome (Sus scrofa 10.2) using Tophat v2.0.13. Gene expression level was calculated by reads per kilobase per million reads (RPKM) and the numbers of reads mapped to each gene were counted by gfold v1.1.2. DEGSeq R package (1.18.0) was applied to determine differentially expressed genes (DEGs) and |log2 (fold change)| >1 and q<0.05 were set as the threshold for significance. RNA-seq data have been deposited in the GEO database and are available through the series accession numbers GSE124484.

Real-time quantitative polymerase chain reaction validation of miRNA and mRNA
Real-time polymerase chain reaction (PCR) was performed on an ABI Step One Plus system (Applied Biosystem, Carlsbad, CA, USA) using SYBR Premix Ex Taq kit (TaKaRa, Dalian, China) with specific primers (Supplementary Table  S1). Glyceraldehyde 3 phosphate dehydrogenase and met-tRNA were chosen as a control of mRNA and miRNA, respectively. Three biological replicates were used for each of the miRNAs and mRNAs. The method of 2 -ΔΔCt was used to calculate fold changes of miRNA and mRNA expression.

Metabolic characteristics of the two porcine breeds
As shown in Table 1, body weight and liver weight of LD pigs were significantly higher than those in JH pigs (p<0.05), however the liver index remained unchanged (p>0.05). The levels of TT4 and TT3 were significantly higher in LD pigs (p<0.05). Serum TCH and LDLC in Jinhua pigs were significantly higher than those in LD pigs (p<0.05), while high density lipoprotein showed no difference between two breeds. Serum glucose and serum TG did not differ between the two pig breeds, while serum insulin levels were significantly low in Jinhua pigs.
Hepatic mRNA expression patterns of two breeds in porcine liver tissue A total of six libraries JH1, JH2, JH3, LD1, LD2, and LD3 were constructed to identify the mRNA expression profiles in JH and LD pigs. As shown in Supplementary  Table S4). In total, 16,051 genes were found to be expressed in the liver of two pig breeds which calculated by RPKM and counted by gfold v1.1.2.
Differentially-expressed genes in the liver tissues of LD (control group) and JH pig breeds were acquired by DEGseq. There were 172 genes significantly upregulated by comparing the JH samples with the LD samples, whereas 295 genes were significantly down regulated (|log2(fold change)| >1 and q<0.05). The volcano map displayed DEGs (Supplementary Figure S2) significantly differentially expressed in the two breeds.  Figure 1.

Functional enrichment of differentially expressed genes
The GO enrichment analysis was applied to the DEGs to explore the biological function between two pig breeds ( Figure  2A). The enriched GO targets were mainly associated with oxidation-reduction, lipid biosynthetic and lipid metabolic process. The DEGs were classified according to KEGG function annotations to identify the pathways ( Figure 2B). The involved pathways included: metabolic pathways, fatty acid metabolism, PPAR signaling pathway, fat digestion and absorption and protein digestion and absorption.

DISCUSSION
The liver is important to vertebrates and some other animals as it plays a major role in carbohydrate, protein, amino acid and lipid metabolism [11]. In this study, significant differences in physiological and biochemical traits were observed between growing LD and JH pigs. The body weight and liver weight indicated discrepancy between cell proliferation and differentiation in two breeds. The levels of serum TT3 and TT4 were significantly higher in LD pigs, which was in agreement with study of Jinhua pigs and Yorkshire pigs [12]. TT3 and TT4 play a role in accelerating the decomposition of fat in adipose tissue and increase the content of free fatty acids in blood [13] which also enhance the function of growth hormone to activate the adenylate cyclase in order to increase steatolysis through cAMP-protein kinase system [14]. The results showed that the content of serum TCH and LDLC were strikingly higher in Jinhua pigs which indicated higher rate of cholesterol metabolism. What's more, leptin produced by adipose tissue participates in the metabolism of lipid in the liver which influences the synthesis of phosphoenolpyruvate carboxylase and gluconeogenesis in order to restrict the synthesis of TG and to enhance the sensitivity of insulin to the liver [15]. We identified 467 DEGs between JH and LD pig breeds. Various of them were associated with lipid metabolic and biosynthetic processes (Figure 2A). There were 6 DEGs involved in PPAR signal pathway through KEGG pathway analysis. Among them, acyl-CoA dehydrogenase long chain and solute carrier family 27 member 4 had higher expression in JH liver which were essential for oxidizing unsaturated fatty acids [16] and transporting fatty acids across the membrane [17]. The higher expression of apolipo protein A4 indicated antioxidant properties and reducing LDLC peroxidation [18], which may respond to serum LDLC in-creased significantly in JH pigs. Differences in biochemical performances and gene expression profiles showed the comprehensive scope of the liver metabolism processes varied between the two breeds.
MiRNAs play important roles in the post-transcription regulation which is involved in various developmental and physiological processes [19]. In our study we found out that several DE miRNAs were related to fat metabolism. Among them, miR-374b directly targets CCAAT enhancer binding protein beta which is an important liver transcription factor involved in the effect of maternal dietary protein on offspring lipid metabolism in Meishan pigs [20]. MiR-370 can control expression of miR-122 and carnitine palmitoyltransferase 1A and then affects lipid metabolism [21]. Inhibition of miR-145 which was abundant expressed in LD pigs regulates ATP binding cassette subfamily A member 1 expression and promotes high-density lipoprotein biogenesis in the liver [22]. MiR-34a acts an inhibitor of beige and brown fat formation by suppressing adipocyte fibroblast growth factor 21 signaling and sirtuin 1 function [23] and our observations indicated that miR-34a was up-regulated in the LD group.
MiRNA-mRNA co-regulatory networks have provided a comprehensive result for understanding the mechanism of lipid metabolism. In the present study, despite the fact that the differences in growth and fat deposition are significant between the two pig breeds, there are still no effective networks research on it. After filtering the target genes predicted by bioinformatics software which differentially expressed in the two pig breeds, miRNA-mRNA potential networks were constructed for more understanding. As shown in the regulatory network (Figure 1), MiR-204 which was down-regulated in JH pigs was reported to be up-regulated during adipogenesis of human bone marrow stem cells and favors adipogenesis when overexpressed [24]. ELL associated factor 1 (EAF1), the predicted target of miR-204 which is expressed higher in JH pigs, negatively regulates canonical Wnt/β-catenin signaling  [25] involved in porcine adipose tissue development [26]. The GLi1 is a transcription factor of the Hedgehog (Hh) signaling pathway which inhibits fat accumulation in drosophila and mammalian models [27]. Acetyl-CoA carboxylase beta plays an important role in fatty acids oxidation [28]. What's more, stearoyl-CoA desaturase (SCD) is related to intramuscular fat deposition through synthesis and desaturation of fatty acids in pigs [29]. MiR-429 which interacted with SCD in our prediction is upregulated during adipogenic differentiation in mouse ST2 mesenchymal stem cells [30]. In our study, miR-370 expressed higher in LD group which influence on oxidation of fatty acids and synthesis of TG [21]. Our findings indicated that the differential expression of genes and miRNAs might contribute to lipid metabolism and fat deposition in pigs.

CONFLICT OF INTEREST
We certify that there is no conflict of interest with any financial organization regarding the material discussed in the manuscript.