Transcriptome analysis of miRNA and mRNA in the livers of pigs with highly diverged backfat thickness

Fat deposition is very important in pig production, and its mechanism is not clearly understood. MicroRNAs (miRNAs) play critical roles in fat deposition and energy metabolism. In the current study, we investigated the mRNA and miRNA transcriptome in the livers of Landrace pigs with extreme backfat thickness to explore miRNA-mRNA regulatory networks related to lipid deposition and metabolism. A comparative analysis of liver mRNA and miRNA transcriptomes from pigs (four pigs per group) with extreme backfat thickness was performed. We identified differentially expressed genes from RNA-seq data using a Cufflinks pipeline. Seventy-one differentially expressed genes (DEGs), including twenty-eight well annotated on the porcine reference genome genes, were found. The upregulation genes in pigs with higher backfat thickness were mainly involved in fatty acid synthesis, and included fatty acid synthase (FASN), glucokinase (GCK), phosphoglycerate dehydrogenase (PHGDH), and apolipoprotein A4 (APOA4). Cytochrome P450, family 2, subfamily J, polypeptide 34 (CYP2J34) was lower expressed in pigs with high backfat thickness, and is involved in the oxidation of arachidonic acid. Moreover, 13 differentially expressed miRNAs were identified. Seven miRNAs were associated with fatty acid synthesis, lipid metabolism, and adipogenic differentiation. Based on comprehensive analysis of the transcriptome of both mRNAs and miRNAs, an important regulatory network, in which six DEGs could be regulated by differentially expressed miRNAs, was established for fat deposition. The negative correlate in the regulatory network including, miR-545-5p and GRAMD3, miR-338 and FASN, and miR-127, miR-146b, miR-34c, miR-144 and THBS1 indicate that direct suppressive regulation may be involved in lipid deposition and energy metabolism. Based on liver mRNA and miRNA transcriptomes from pigs with extreme backfat thickness, we identified 28 differentially expressed genes and 13 differentially expressed miRNAs, and established an important miRNA-mRNA regulatory network. This study provides new insights into the molecular mechanisms that determine fat deposition in pigs.

regulator of adipocyte differentiation via interaction with the 3ʹ-UTR of peroxisome proliferator activated receptor gamma (PPARγ) mRNA. MiR-27a was also able to accelerate the hydrolysis of triglyceride (TG), and miR-143 was shown to promote TG synthesis during lipid metabolism in porcine adipocytes 9 . In the adipose tissue, miR-NAs are differentially expressed in pigs with different genotypes (lean and obese), with miR-9, miR-124a, miR-1a, miR-133a, miR-122, miR-204, and miR-183 all being upregulated in obese pigs. Conversely, the expression of miR-215, miR-135, miR-224, and miR-146b was higher in lean pigs than in obese pigs 10,11 .
The pig (Sus scrofa) is an important agricultural animal, and a cost-effective source of meat for human consumption, particularly in China, where pork is the most common meat consumed. Levels of fatness and leanness are very important in pig production because of their ability to affect productive performance, meat quality, and reproductive performance 12,13 . Due to their anatomical, physiological, and genetic similarities to humans, and their high level of fat deposition, pigs represent a suitable model for use in the study of human obesity and energy metabolism 14 . Landrace pigs are a typical lean-type western breed, which have been intensively selected over the past few decades to increase lean meat production and reduce fat deposition 15 . Therefore, fat deposition and metabolism in pigs are always an immediate area of research focus.
The liver is the largest digestive gland and most important metabolic organ. It plays key roles in the regulation of appetite, body weight, and several metabolic processes 16 . In pigs, the main site for de novo fatty acid synthesis is the adipose tissue 17 . However, polyunsaturated fatty acid (PUFA) synthesis, de novo cholesterol synthesis, and fatty acid oxidation mainly occur in the liver 18,19 . The liver is also critical for glucose metabolism. When blood glucose levels are low, the liver releases glycogen 20 . Liver is suitable tissue to study the molecular mechanisms of fat deposition and metabolism in pigs.
To understand the genes and miRNAs involved in the molecular mechanism of fat deposition in pigs, in the present study we used mRNA and miRNA sequencing to identify differences in liver tissue transcriptomes between pigs with opposing backfat thickness phenotypes. A number of significantly differentially expressed genes and miRNAs were found in the liver tissue from pigs with high and low backfat thickness. Subsequently, the expression patterns and co-expression of differentially expressed mRNA were analysed. Furthermore, we assessed mRNA-miRNA interactions using computational prediction and expression relationship analysis. Finally, integrated transcriptome analysis of mRNA and miRNA in livers was performed in order to identify important networks for differential backfat thickness in pigs.

Results
Animal phenotypes. In this study, 132 female landrace pigs, which were on average 186-days-old with an average live body weight of 93.38 kg, were used. Backfat thickness measured 5.76 ± 1.75 mm on average (Table S1). From these data, two groups were selected with low or high backfat thickness respectively. The genealogy information of selected animals is shown in Supplementary Table S2. Backfat thickness of the BFH group (8.88 ± 0.80 mm) was significantly higher than that of the BFL group (3.58 ± 0.39 mm). Additional traits assessed related to fat deposition, such as carcass backfat thickness, kidney fat weight, and intramuscular fat, also exhibited significant differences between the two groups. Individual differences in genetic is the important factor for leading the difference in live fat thickness so significant. There was no difference in body weight between the two groups, indicating that changes in fat deposition were not driven by significant changes in body weight (Table 1). mRNA expression profiles. In this study, eight cDNA libraries, which included four BFH pigs and four BFL pigs, were sequenced from the liver tissue using Illumina HiSeq. 2000. The deep sequencing data of total RNA have been submitted to National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under Accession no. SRP117778, Bioproject: PRJNA407236. After removing the adaptors and low-quality reads, 23,102,658−27,426,569 clean paired reads, were obtained, with 82.10−86.20% of clean reads being mapped to the porcine reference genome (Table S3). After normalization by FPKM, 19,279 genes were found to be expressed in the liver tissue of all eight pigs. Seventy-one genes were found to be differently expressed between two groups using Cuffdiff (Table S4), with the criteria of at least a 2-fold difference and a q-value of less than 0.05 (Fig. 1). Twenty-eight of the differentially expressed genes were well annotated on the porcine reference genome, including 15 upregulated and 13 downregulated genes for the BFH group (Table 2). Those DEGs not annotated transcript were analyzed through homology with reference genome (GRCh38/hg38). Twenty one of not annotated transcripts were mapped in human reference genome, and eleven genes were found in those locations (Table S5). We also identified 127 DEGs between pigs in half-sib pairs, and 114, 68, and 63 DEGs between pigs in each  Figure S1). However, 29 common genes, such as APOA4, CXCL2, FCN, SAA and so on, were found in more than one pairs. Because of the sample limitations of the DEGs for each pair, we only present DEGs and their functional analyses obtained by treating animals as four biological replicates.
Functional analysis of differentially expressed genes. Following conversion to human ortholog ensemble gene IDs, 25 DEGs were annotated in DAVID for pathway and GO analysis (Table S7). Focal adhesion, vascular smooth muscle contraction, and tight junction were identified by KEGG pathway analysis. As little of DEGs, none pathway was found significantly (P < 0.05). For GO functional enrichment analysis, a total of 10 GO terms were enriched (Table 3) in categories related to lipid metabolism (e.g., regulation of lipid transport, regulation of lipoprotein metabolic process, and high-density lipoprotein particle remodelling), and general metabolism (e.g., carbohydrate binding, sugar binding, and regulation of steroid metabolic processes). Gene interaction networks were constructed using STRING online software (Fig. 2). All DEGs were input to the STRING, and those discrete genes in the network was removed.
miRNA expression profiles. The deep sequencing data of total RNA have been submitted to NCBI SRA under Accession no. PRJNA407236, Bioproject: PRJNA407256. An average of 19.49 million raw reads were obtained from each individual. After sequencing analysis, 334 miRNAs from all liver samples were identified. In our sequencing data, the miRNAs that were common to the top most abundant miRNAs of all liver samples were: ssc-miR-192, ssc-miR-26a, ssc-miR-148a-3p, ssc-miR-22-3p, ssc-miR-143-3p, ssc-miR-122, ssc-miR-101, ssc-miR-30a-5p, ssc-miR-10a-5p, ssc-miR-27b-3p, ssc-miR-126-5p, ssc-miR-191, ssc-let-7a, ssc-miR-181a, ssc-miR-27a, ssc-let-7f, ssc-miR-92a, ssc-miR-21, ssc-miR-30e-5p, and ssc-miR-194a (Fig. 3). Compared with the BFL group, 13 miRNAs were differentially expressed, including four that were upregulated and nine that were downregulated in the BFH group (Table 4). To further elucidate the functions of the differentially expressed miRNAs, we searched their potential target genes. A total of 3997 putative target genes of 13 DEM were predicted using TARGETSCAN (Table S8). Some putative target genes were repeated. Overall, 2827 target genes were found, and were remained for further analysis. After GO and KEGG pathway analysis, 44 KEGG pathways were identified as being significantly related to genes targeted by differentially expressed miRNAs. Almost all of these pathways are associated with activities that are fundamental to life. However, multiple pathways and GO terms were associated with energy metabolism and fat deposition, such as adipocytokine signalling pathways, regulation of glucose metabolic processes, and lipid binding (Fig. 4). www.nature.com/scientificreports www.nature.com/scientificreports/ Validation of DEGs and DE miRNAs by qPCR. We selected seven DEGs (FASN, THBS1, ACTA2, LIPG, GCK, CYP2J34, and APOA4) and six DE miRNAs (miR-338, miR-127, miR-451, miR-34c, miR-452, and miR-136) randomly to validate the accuracy of the RNA-seq and miRNA-seq using qPCR, including the up-and downregulated genes and miRNAs. Furthermore, the fold changes of the seven genes and the six miRNAs in the qPCR and in the RNA-seq or miRNA-seq showed the same trends (Fig. 5). Most of them were significant differently expressed between two groups (p < 0.05). The correlation between RNA-seq and qPCR is 0.84 (p < 0.05), and the correlation between microRNA-seq and qPCR is 0.91 (p < 0.05). Those results indicated that the DEGs and DE miRNAs identified with NGS were reliable and efficient.
Construction of the miRNA/mRNA network. To further analyse the relationship between the node mRNAs and non-coding miRNAs, a potential network of DEGs and differentially expressed miRNAs interacting in pig liver that might affect fat deposition was constructed (Fig. 6). The node mRNAs and non-coding miRNAs  www.nature.com/scientificreports www.nature.com/scientificreports/ indicate that those mRNAs and non-coding miRNAs have potential regulatory relationships. In the miRNA/ mRNA network, six DEGs, including GRAM domain containing 3 (GRAMD3), FASN, desmin (DES), thrombospondin 1 (THBS1) and lipase G, endothelial type (LIPG), could be regulated by differentially expressed miRNAs. MiR-338 targeted FASN, and their expression level shown negative relation (correlation coefficient = −0.52). The correlation coefficient between miR-127, miR-146b, miR-144, miR-34c and their potential targeted gene THBS1 Figure 2. STRING analysis shows that differentially expressed genes are involved in known and predicted proteinprotein interactions. STRING is used to analyse differentially expressed genes in the livers of pigs from the BFH and BFL groups. Those up-regulated RNAs were red, and those down-regulated RNAs were green to BFL.

Discussion
In this study, three full-sib pairs of pigs and one half-sib pair of individuals with extreme backfat thickness phenotypes were used. Landrace pigs is an lean-type pigs, and have an low level of fat deposition as an whole. But phenotypic variation of fatty traits such as backfat thickness still was found in this population. There were no differences in the body weight of pigs from different groups; however, significant differences in traits related to fat deposition such as carcass backfat thickness, kidney fat weight, and intramuscular fat, were observed. This suggested that the significant changes in fat deposition were not mediated by changes in body weight. To our knowledge, this is the first study to perform parallel comparisons between mRNA and miRNA transcriptomes in the liver in order to investigate the function of mRNAs, the roles of miRNAs, and their targets in lipid deposition and metabolism. In contrast with previous studies [21][22][23] , the use of more full-sib pairs may reduce the noise associated with differences in the digenetic background and the generation of false-positive results. Comparing the result of each pair, there were more DEGs in half-sib pair than that in any full-sib pair. Thus, selected full-sib pigs represent a better sampling strategy to assess differences in gene expression.  Table 4. Thirteen differentially expressed miRNAs in Landrace pigs with low and high fat deposition. www.nature.com/scientificreports www.nature.com/scientificreports/ After removal of low-quality reads, between 82.1% and 86.2% of the obtained sequence reads, each sample was mapped to the porcine reference genome. These values were higher than those reported in previous studies investigating the porcine liver transcriptome: 72.63 24 , 64.46 25 , 71.42-77.75 26 and 78.29-79.60% 27 . This represents a good result for RNA-seq data analysis. The total number of genes expressed in the liver was similar between groups (BFH = 18,708, BFL = 18,640), and 18,051 common genes were expressed in both groups. This suggests that most of the liver transcriptome is conserved between animals of the two groups. Only 28 DEGs were identified in the livers of pigs with significantly different backfat thickness. The number of DEGs was lower than that reported in several studies that performed porcine liver transcriptomic analyses using RNA-seq to identify genes related to fat deposition: 169 23   The transcript expression fold changes measured by qPCR and NGS are indicated by dark grey and light grey columns, respectively. Asterisks on the qPCR values indicate significant differences between BFH and BFL at p < 0.05. Figure 6. Key network of genes and miRNAs differentially expressed in the livers between pigs with extreme backfat thickness. The network diagram was generated using Cytoscape. Those up-regulated RNAs or miRNAs were red, and those down-regulated RNAs or miRNAs were green to BFL. general metabolism. Due to the low number of DEGs, nominal p-values (P < 0.05) rather than corrected p value (P < 0.05) is the threshold value as significance level in GO enrichment analysis. Several GO terms, such as regulation of lipid transport, regulation of lipoprotein metabolic process, and carbohydrate binding, were identified as significant categories. APOA4 and LIPG are involved in processes related to several GO terms, such as those related to lipid metabolism and general metabolism [29][30][31] , and were up-regulated in the BFH group in this study. This suggests that fat content in the liver is subjected to constant regulation by negative feedback. GCK, which was more highly expressed in the BFH groups, is involved in the first step of glucose metabolism, and controls gluconeogenesis and glycogen synthesis in the liver 32 . This shows that the livers of pigs with high backfat thickness possess a high level of fatty acid synthesis. Hence, differences in the rates of de novo fatty acid synthesis, and lipid metabolism in the liver between leaner pigs, compared with fatter pigs, are attributable to genetic factors resulting in differential fat deposition.
In this study, liver microRNAomes were compared between BFH and BFL, and 13 differentially expressed miRNAs were identified. Several differentially expressed miRNAs were considered to play regulatory roles in fatty acid synthesis, lipid metabolism, and adipogenic differentiation, and similar results have been reported previously. For example, miR-338 was specifically detected in the adipose tissue of leaner pigs 33 , and was also found to be significantly down-regulated in the adipose tissue of obese pregnant women 34 . MiRNAs differentially expressed in Large White pigs (lean type pig) and Meishan pigs (Chinese indigenous fatty pig) have been identified through a high-throughput Solexa sequencing approach. The results showed that the expression of miR-136 in Meishan pigs was significantly higher than that in Large White pigs 10 . Expression of miR-144 has been reported to cause diseases related to lipid metabolism 35,36 . Abundance expression of miR-4332 was shown to be higher in lean pigs (Landrace and Yorkshire) than in obese pigs ((Diannan Small-ear pig and Tibetan pig) 37 . These data are consistent with the results of the present study.
The miRNA-mRNA regulatory networks identified in this study provide a comprehensive profile that can help to delineate the mechanisms of fat deposition. MiRNA acts a negative regulator of gene expression in fat formation 38 . GRAMD3 has been reported as a candidate gene associated with ectopic fat 39 , and is a potential target gene of miR-545-3p. GRAMD3 and miR-545-3p display opposing patterns of expression in the two groups. Highly negative correlation between GRAMD3 and miR-545-3p (−0.83) means miR-545-3p may regulate fat deposition through targeted GRAMD3. FASN is a key lipogenic enzyme and the rate-limiting step in de novo fatty acid synthesis 40 . Consistent with previous reports 41 , FASN expression was higher in the fatter pigs in our study. MiR-338 was upregulated in the liver of animals in the BFL group, and is predicted to target FASN. This implies the expression of FASN could be negatively regulated by miR-338. THBS1 represents the hub molecules in the network, plays a causal role in the pathogenesis of insulin resistance and adipose tissue inflammation, and serves as a biological marker of obesity and the metabolic syndrome 42 . MiR-127, miR-146b, miR-34c and miR-144 have been predicted to regulate the level of THBS1 expression. In our study, miR-127, miR-146b, miR-34c and miR-144 were upregulation in the liver of animals in the BFH group, whereas THBS1 was downregulated in animals of this group. Those potential negative regulation between miR-127, miR-34c, miR-144 and THBS1 may play important roles in fat deposition. As it represents the target off our identified differentially expressed miRNAs, THBS1 represents a new candidate gene influencing backfat thickness in pigs. However, the expression level of proteins and the specific functions of these crucial genes and miRNAs have not been studied in depth, so further studies are needed to understand fat deposition in pigs.

Conclusion
In the present study, we identified 28 genes and 13 miRNAs that were differentially expressed in the livers of Landrace pigs with extreme and divergent backfat thickness in this population. Genes functionally related to fatty acid synthesis, glucose metabolism, and fatty acid degradation were found to be differentially expressed. Additionally, a few miRNAs that target genes involved in fatty acid synthesis, lipid metabolism, and adipogenic differentiation were also found to be differentially expressed. Based on the comprehensive transcriptome analysis of both mRNAs and miRNAs, we conclude that the levels of fatty acid synthesis and degradation in liver are due to the degree of fat deposition. Adipogenic differentiation in the liver regulated by miRNAs could influence body fat of pigs. Overall, we modelled the miRNA-mRNA regulatory network related to fatty acid synthesis and lipid transport using data on differentially expressed miRNAs and mRNAs from Landrace pigs with extreme and divergent backfat thickness. The negative correlate including, miR-545-5p and GRAMD3, miR-338 and FASN, and miR-127, miR-146b, miR-34c, miR-144 and THBS1 indicate a possible direct suppressive regulation between miRNAs and their target mRNAs. The overall results of this study facilitate our understanding of the molecular mechanisms regulating fat deposition in pigs and provide an insight into a regulatory role for miRNAs in energy metabolism.

Materials and Methods
Animals and phenotypes. A landrace female pig resource population was housed in the Tianjin Ninghe Primary Pig Breeding Farm (Ninghe, China) under consistent and standard environmental conditions. Animals were fed three times a day and had free access to water. A total of 132 individuals were measured at the last 3rd and 4th ribs using real-time B-mode ultrasonography (HS1500 convex scanner, Honda Electronics, Toyohashi, Japan) to determine live backfat thickness. Age, body weight, and pedigree information were also available for the resource population. Selection methods and standards included two criteria: live backfat thickness in pigs from the high backfat thickness group (BFH) was at least twice that of pigs from the low backfat thickness group (BFL), and the pairs of pigs with divergent backfat thickness, which were full-sibs, were the first choice for selection. Four pairs of pigs with divergent backfat thickness phenotypes, three of which were full-sibs, were selected, slaughtered, and used for the collection of livers. Each pairs were traced three generations in pedigree, and had no relationship between ancestors. (2019) 9:16740 | https://doi.org/10.1038/s41598-019-53377-x www.nature.com/scientificreports www.nature.com/scientificreports/ After being stunned with a captive bolt and exsanguinated, four pairs of selected pigs were slaughtered in a commercial abattoir (Beijing Huadu Sunshine Food Co., Ltd., Beijing, China). All chosen pigs were stunned with a captive bolt, exsanguinated and slaughtered in commercial abattoir called Beijing Huadu Sunshine Food co., LTD. The chosen pigs were slaughtered according to guidelines of operating procedures of pig-slaughtering (GB/T 17236-2008), which was promulgated by General Administration of Quality Supervision, Inspection and Quarantine of the People's Republic of China (AQSIQ) and Standardization Administration of the People's Republic of China (SAC). All efforts to minimize animal suffering were made during the study. The whole procedure for collection of the tissue samples of all animals was by our researchers. This study was approved specifically by the Animal Welfare Committee of China Agricultural University (Permit number: DK996). All protocols and procedures involving animals were performed in accordance with the Regulations for the Administration of Affairs Concerning Experimental Animals. The liver tissue was collected aseptically and stored in liquid nitrogen immediately after slaughter until required for RNA isolation. The backfat thickness of carcasses at the last 3-4 th ribs and leaf fat weight were measured on slaughtered pigs. The percentage of intramuscular fat in the longissimus muscle was quantified by a near infrared technique Foodscan (Foss CO. Ltd., Sweden) and the value was expressed as the weight percentage of wet muscle tissue. The significance of the difference in fatness traits between the two groups was assessed using t-test by R 43 .
RNA isolation, library preparation, and sequencing. Total RNA was isolated from eight livers using Hiseq. 2000 sequencing system to generate 2 × 90 PE reads. High-quality reads of PE were obtained for each sample. Adaptor sequences, reads with more than 10% unknown sequences, "N", and low-quality sequences (the percentage of low-quality bases with a threshold quality score <20) were removed. The obtained sequence reads were quality-checked by FastQC 44 . Clean reads were aligned to the porcine reference genome sequence (Sscrofa11.1) using TopHat v2.0.1 software 45 . The porcine genome (Sscrofa11.1) was obtained from Ensembl database. Reads aligned to the reference genome were assembled by Cufflinks software 46 . Cuffdiff, a part of the Cufflinks package, was used to identify differentially expressed genes (DEGs) and transcripts from two groups of pigs with opposing backfat thickness performance 47 . The threshold value for selection of differentially expressed genes is q-value ≤ 0.05 and fold change (FC) ≥ 2 or ≤ 0.5, as |Log 2 FC| ≥ 1. Moreover, DEGs between pigs in each pair were also analyzed using the same analysis pipeline. miRNA sequencing and statistical analysis. Eight small RNA libraries were prepared using the TruSeq ® Small RNA Sample Prep Kit (Illumina ® ). All of the procedures and standards were performed as described in the instructions supplied with the kit. After quality control, eight libraries were sequenced on a HiSeq. 2000 platform (Illumina) and 50-bp single-end reads were obtained. Before reads were analysed, the adapters were removed and low-quality reads were discarded using the Fastx-toolkit for quality control 48 . The CAP-miRSeq pipeline was used for miRNA analysis 49 . MiRNA data were converted from fatsq to fasta format using awk in the Linux system. MirDeep2 was used to analyse miRNA sequencing data 50 . The miRNA reference was downloaded from the miRBase database 4 . The expression of microRNAs in different libraries was normalised by Trimmed Mean of M-values (TMM) 51 , and differentially expressed miRNAs (P-value ≤ 0.05, and FC ≥ 2 or ≤ 0.5, as |Log2 FC| ≥ 1) were identified using the edgeR package 52 .
Quantitative real-time PCR (qPCR) of mRNA and miRNA. Reverse transcription quantitative real-time PCR were used to confirm the sequencing data by Light Cycler ® 480 Real-Time PCR System (Roche, USA). Seven differentially expressed genes and six differentially microRNAs were selected randomly. Total RNA was isolated from four pairs pigs and converted into cDNA using the Revert Aid ™ First Strand cDNA Synthesis Kit (Thermo Fisher Scientific Inc, USA). The reaction volume of 20 μL was used in the qPCR reactions according to the manufacturer's protocol. The primers used are shown in Table S9. Triplicate qPCRs were performed for each cDNA and the average Ct was used for further analysis. The 2 −ΔΔCt method was used to determine the relative mRNA abundance 53 . To confirm changes in miRNA levels between BJ and HL groups, parallel qPCRs were also used to quantify relevant miRNAs. Six differential expressed microRNAs were selected for qPCR analysis. Chicken U6 was chosen as an endogenous control to correct for analytical variations. This experiment was performed by Beijing SinoGene Scientific Co., Ltd.
Functional analysis. Gene IDs were converted to homologous human Ensembl Gene Ids using bioDBnet 54 (https://biodbnet-abcc.ncifcrf.gov/) and subjected to functional analysis. The Database for Annotation, Visualization and Integrated Discovery (DAVID) 55 bioinformatics resources (http://david.abcc.ncifcrf.gov/) were used for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) enrichment analysis of differentially expressed genes using nominal p-values (P < 0.05) including molecular functions (MFs), biological processes (BPs), and cellular components (CCs). An interaction network of DEGs was structured