Transcriptome analysis of mRNA and miRNA in skeletal muscle indicates an important network for differential Residual Feed Intake in pigs

Feed efficiency (FE) can be measured by feed conversion ratio (FCR) or residual feed intake (RFI). In this study, we measured the FE related phenotypes of 236 castrated purebred Yorkshire boars, and selected 10 extreme individuals with high and low RFI for transcriptome analysis. We used RNA-seq analyses to determine the differential expression of genes and miRNAs in skeletal muscle. There were 99 differentially expressed genes identified (q ≤ 0.05). The down-regulated genes were mainly involved in mitochondrial energy metabolism, including FABP3, RCAN, PPARGC1 (PGC-1A), HK2 and PRKAG2. The up-regulated genes were mainly involved in skeletal muscle differentiation and proliferation, including IGF2, PDE7A, CEBPD, PIK3R1 and MYH6. Moreover, 15 differentially expressed miRNAs (|log2FC| ≥ 1, total reads count ≥ 20, p ≤ 0.05) were identified. Among them, miR-136, miR-30e-5p, miR-1, miR-208b, miR-199a, miR-101 and miR-29c were up-regulated, while miR-215, miR-365-5p, miR-486, miR-1271, miR-145, miR-99b, miR-191 and miR-10b were down-regulated in low RFI pigs. We conclude that decreasing mitochondrial energy metabolism, possibly through AMPK - PGC-1A pathways, and increasing muscle growth, through IGF-1/2 and TGF-β signaling pathways, are potential strategies for the improvement of FE in pigs (and possibly other livestock). This study provides new insights into the molecular mechanisms that determine RFI and FE in pigs.

Feed accounts for more than 60% of the costs for pig production, therefore improving feed efficiency (FE) is one of the major ways to reduce costs in the pig industry. FE can be measured as feed conversion ratio (FCR) or residual feed intake (RFI) 1 . FCR is the feed intake divided by the weight gained during a specified period. RFI is defined as the difference between the actual and the predicted dry matter (DM) intakes of each animal, based on its metabolic body weight and average weight gain during a specified period 2 . Thus, animals with higher RFI/FCR are less efficient at converting feed into body mass, whereas those with lower RFI/FCR are more efficient. Previous studies indicated that the heritability of RFI is 0.14-0.40 and FCR is 0.13-0.31 [3][4][5] , and a strong correlation exists between them (R equals 0.76-0.99) 3 .
With microsatellite typing based QTL mapping, Zhang and his colleagues identified three genomic regions on SSC2, SSC7 and SSC9 in a White Duroc × Chinese Erhualian F2 segregated population 6 , associated with the feed consumption and feeding behavior traits, average daily feed intake (ADFI), feed conversion ratio (FCR), number of visits to the feeder per day (NVD) and average feeding rate (AFR).

Results
Animal performance. An original 238 castrated purebred Yorkshire boars were grown from 30 to 90 kg (average body weight) on a commercial pig feed (see Materials and Methods). Daily feed intake (DFI) and average daily gain (ADG) were measured and the FCR and RFI determined. From this data, two sub-groups were selected with Low or High RFI (RFI_L and RFI_H respectively). Importantly, the total weight gained and the time taken to do so was not significantly different between the two groups. The High RFI group (RFI_H) had RFI and FCR values of 0.43 ± 0.19 (Kg/day) and 3.07 ± 0.15 respectively, compared with − 0.28 ± 0.07 (Kg/day) and 2.19 ± 0.06 for the more efficient, Low RFI group (RFI_L, p < 0.05 for both RFI and FCR, Table 1). The lower FCR seen in the RFI_L group was achieved by a reduced DFI (p = 0.02) and, interestingly, a trend for an increase in ADG (p = 0.08) compared with the RFI_H group (Table 1). Although the average metabolic body weight gain (AMBW) was higher in the RFI_L group, it was not statistically significant (p = 0.65). Two assessments of body fat, average back fat (ABF) thickness and intramuscular fat (IMF) content, showed no significant differences, suggesting that changes in RFI or FCR were not mediated by significant changes in body fat.
RNA sequencing data mapping and annotation. A total of 6 cDNA libraries were sequenced from the longissimus dorsi muscle of RFI_L and RFI_H groups (n = 3 from each). The RNA sequence reads have been submitted to the NCBI Gene Expression Omnibus under accession E-MTAB-2946. After removing the adaptors and filtering, RNA-seq yielded from 7.8 M to 9.6 M single-end reads for all 6 samples, more than 99.99% reads were qualified (Supplementary Table S1).
After mapping clean reads to the porcine genome, 89.03%-91.11% reads were successfully aligned, with 51.5-58.4% of reads mapped to CDS regions, 0.53-0.59% of reads mapped to 5'UTR regions, 12.36-14.37% of reads mapped to 3'UTR regions and 10.14-12.66% of reads mapped to introns or intergenic regions, while 8.83-11.00% of reads had multiple alignments (Supplementary Figure S1a). The reads distribution in each chromosome were analyzed according to the porcine genome, and about 10.19-21.29% of the total reads mapped to the mitochondrial genome, which was higher than the proportion of reads mapped to other chromosomes (Supplementary Figure S1b). The mitochondrial genome distribution of the RFI_H group was 19.73 ± 4.98%, compared with 14.33 ± 3.06% in the RFI_L group (Supplementary Figure S1c). This indicates that genes encoded for by the mitochondrial DNA account for a large proportion of the genes differentially expressed in skeletal muscle.
Scientific RepoRts | 5:11953 | DOi: 10.1038/srep11953 miRNA sequencing data mapping and annotation. A total of 2 cDNA libraries were sequenced from the longissimus dorsi muscle of RFI_L and RFI_H groups (1 pool of n = 5 for each group), each cDNA library included the 3 samples used for RNA sequencing). After sequencing, a total of 13,307,888 reads were obtained from the RFI_H group and 11,744,483 reads from the RFI_L group. The miRNA sequence reads have been submitted to the NCBI Gene Expression Omnibus under accession E-MTAB-2949. After removing reads with non-canonical letters or with low quality, the 3' adapter was trimmed and the sequences shorter than 18 nt discarded. In total, 12,503,358 (RFI_H), and 10,408,712 (RFI_L) clean reads were obtained, which corresponded to 93.95% and 88.63% respectively of the raw reads from each small RNA library. The length distribution of clean reads showed that most of the reads were between 21-23 nt in length, and read counts with 22 nt were highest (Supplementary Figure S2a).
A total of 402 mature miRNAs were identified. Among them, 152 were annotated porcine miRNAs already present in miRbase v20, 153 were miRNAs homologous to human or mouse, while 97 were novel miRNAs not homologous to any other species (Supplementary Figure S2b). The distribution of miRNAs on each chromosome depended on the number of annotated miRNAs on each chromosome (Supplementary Figure S2c).
Differentially expressed genes between low and high RFI pigs. In the present RNA-seq study, 30,484 genes were detected in the skeletal muscle of all 6 individuals. A total of 645 genes were differentially expressed, with a criteria of at least a 2 fold difference and a p-value less than 0.05 (|log2FC| ≥ 1, p < 0.05), of which 99 genes had a q-value ≤ 0.05. Of the 99 differentially expressed genes, 60 genes were well annotated on the Sus Scrofa genome (version Sscrofa10.2.72), with 45 genes being up-regulated and 54 genes down-regulated in the low-RFI group (Supplementary Table S2). Table 2 shows the top 20 differentially expressed (DE) genes, the top 10 genes with either higher or lower expression in low RFI compared to high RFI pigs.
To validate the differential expression of genes, eight genes were selected for qRT-PCR analysis. Compared with the RFI_H group, expression of FABP3, RCAN, PPARGC1, and PRKAG2 mRNA were all lower in RFI_L muscles, whereas expression of IGF2 mRNA was higher in RFI_L muscles (Fig. 1a). Expression of 4 of the selected genes (IGF2, RCAN, PPARGC1 and PRKAG2) showed significant differences between the RFI_H and RFI_L groups. Hence the qRT-PCR analyses largely confirmed the RNA-seq data, with the correlation coefficient of the fold-change (FC) values from the two methods being 0.99 and the R 2 for the linear regression also being 0.99, indicating the reliability of the RNA-seq analysis.
In addition, we analyzed the expression levels of various mitochondrial genes. There were 17 mitochondrial coding genes detected (log2FC ≥ 0.5) and 16 genes (94.1%) were down regulated in the skeletal muscle from low RFI pigs (Fig. 1b Table 1. Animal performance of Yorkshire pigs used in RNA and miRNA sequencing. DFI -daily feed intake. ADG -average daily gain over the assessed feeding period. BW -body weight. TBG -total body weight gain (Kg) during the assessed feeding period. ABF -average of back fat thicknesses (mm) measured at three points between 6 th and 7 th ribs (6 th -7 th BF) and at the10th rib (10 th BF). LMA -loin muscle area (cm 2 ) measured between the 10 th and 11 th . IMF -intramuscular fat, percentage IMF determined by the petroleum ether extraction method. p-value as calculated by t-test. reads count no less than 20 (|logFC| ≥ 1, total ≥ 20, p ≤ 0.05). Of these 25 miRNAs, 14 were up-regulated and 11 down-regulated in RFI_L pigs ( Table 3). Six of them were not found in the Sus scrofa miRNA database (miRbase v20), but are homologous with human. Six novel miRNAs were identified as being differentially expressed; the secondary structures of which were predicted by mirDeep2 (Supplementary Figure S3) and the predicted scores from mirDeep2 listed in Supplementary Table S4. To validate the differential expression identified by the miRNA sequencing, miR-1, miR-30e, miR-10b and miR-145 were selected for qRT-PCR analysis. Compared to RFI_H muscles, expression of miR-1 and miR-30e miRNAs was higher, whereas expression of miR-10b, and miR-145 was lower in RFI_L muscles (Fig. 1c). All the selected miRNAs showed significant differences between RFI_H and RFI_L groups, confirming the results of miRNA-seq, and the correlation coefficient for the fold-change values (FC) from the two methods was 0.99 and the R 2 for linear regression was also 0.99, indicating the reliability of the miRNA-seq analysis. The differentially expressed miRNAs were functionally involved in energy metabolism processes and pathways, which were all down regulated in low RFI pigs. Pathway analysis of differentially expressed genes. For gene ontology (GO) biological processes analysis, 54 well annotated genes among the 99 significant DE genes were submitted and 45 of these genes were available in DAVID v6.7.
A total of 26 GO terms were enriched (Fisher Exact Probability Test, p < 0.05), which can be divided into 5 major groups: 1) Genes relating to general metabolism, 2) Genes relating to energy metabolism, 3) Genes relating to lipid metabolism, 4) Genes relating to cell differentiation, and 5) Genes relating to biosynthetic processes (Table 4).
Most genes related to energy metabolism, cell differentiation and biosynthesis were down regulated in pigs with low RFI. For example PPARGC1, PRKAG2, HK2 are related to energy metabolism; NOR-1, RCAN, ESRRB, ESRRG and JUNB are related to cell differentiation; and LIPG and GM2A are related to lipid catabolism. The significant down regulation (p ≤ 0.05) of all those genes indicates that low RFI pigs may be more efficient due to (i) a decrease in energy expenditure because of an inhibition of metabolic processes and/or (ii) an increase in skeletal muscle growth. Pathway analysis of differentially expressed miRNAs. To better understand the biological functions of the 25 differentially expressed miRNAs identified, we predicted the potential target genes of these miRNAs. There were 12,687 terms of target genes returned, including 6,744 unique genes (Supplementary  Table S5). Through the miRNA-targeted pathway union analysis, we found 55 KEGG pathways significantly (Fisher Exact Probability Test, p < 0.05) related with genes targeted by up or down-regulated miRNAs (Supplementary Table S6). A lot of pathways were involved in energy metabolism and skeletal muscle growth, including the TGF-beta signaling pathway, PI3K-Akt signaling pathway, mTOR signaling pathway, GnRH signaling pathway and Hypertrophic cardiomyopathy (HCM).
To further classify and predict the function of differentially expressed miRNAs, we also performed hierarchical clustering of differentially expressed miRNAs and their target pathways (Fig. 2). Some miR-NAs with the same regulation pattern or similar function were clustered together, for example miR-130a-3p and miR-301b-3p were clustered together. After checking the mature sequences of these two miRNAs, we found they have exactly the same seed sequences and both of them are in the same miRNA family.
A key network of mRNA and miRNA in pig muscle potentially regulates RFI in pigs. Cytoscape v3.0.1 was used to integrate a potential network of differentially expressed genes and miRNAs interacting in pig skeletal muscle that might lead to differences in RFI (Fig. 3). When looking into the differentially expressed genes, we found that most of those genes were involved in mitochondrial activity, glycolysis or myogenesis pathways and were actually connected directly or indirectly through just one or two genes. In this network, the mitochondrial activity was separated into 3 parts: the uncoupling reaction, the mitochondria respiratory control and the mitochondria transcriptional control.
In mitochondrial activity, the expression levels of PGC-1, ESRRB and TFAM were down-regulated in the low RFI pigs more than other genes, which suggests that the down-regulation of mitochondrial activity might possibly be through mitochondria transcriptional control.
In the glycolysis pathway, the upstream genes of CREB were also involved in mitochondria transcriptional control. Most of these genes were down-regulated in low RFI pigs, such as PGC1, PPAR and PRKAG2. HK2, the gene downstream of CREB, was also down-regulated in low RFI pigs. HK2 is a hexokinase that phosphorylates glucose to produce glucose 6-phosphate 18 . Hexokinase regulates glycolysis as the first rate-limiting enzyme 19 . Moreover, the miRNAs, which were reported or predicted to target and therefore inhibit these genes, were all up-regulated in low RFI pigs, including mir-30e, mir-301b, mir-130a, mir-335 and mir-199b. This suggests an overall decrease in glycolysis and mitochondrial activity in skeletal muscle of low RFI pigs.
In contrast, most genes involved in the myogenesis pathway were up-regulated in low RFI pigs, such as IGF-2, ITGA5, PIK3R1 and MYHCs. Also, TGF-β, which is an inhibitor of myogenesis, was down-regulated in low RFI pigs. Interestingly, mir-141, which targets IGF-2, was down-regulated and   mir-30e, which targets TGF-β , was upregulated in low RFI pigs. These results indicate that the increased efficiency of low RFI pigs is associated with enhanced skeletal muscle growth.

Discussion
Two groups of pigs were identified which had either high or low residual feed intake. Those pigs in the RFI_L group were more efficient, having a lower FCR than the RFI_H group, due to lower DFI and slightly higher ADG. A common observation in these types of study is that selection for greater feed efficiency often targets reduced feed intake, which has been reported as undesirable as it can reduce growth and reproductive performance 20,21 . However in this study there was trend for a positive effect on ADG in the RFI_L pigs, despite a reduction in feed intake. From the estimates of fat deposition, there were  no changes in body fat, although the values for back fat and IMF tended to be lower in the RFI_L pigs.
In addition the AMBW of the RFI_L was higher than RFI_H, but this was not statistically significant. When combined with the assessment of body fat, this suggests an increase in lean tissue growth and a corresponding decrease in fat in the more efficient RFI_L pigs. We identified 99 differentially expressed genes and 25 differentially expressed miRNAs in LD muscle from pigs with significantly different RFI. These genes and miRNAs corresponded to two key pathways/ functions, one related to mitochondria and energy metabolism and the other related to skeletal muscle growth. The energy metabolism and growth of skeletal muscle may be the two key factors responsible for low RFI and therefore increased efficiency of pigs. Mitochondria produce more than 90% of cellular energy through oxidative phosphorylation (OXPHOS), but also waste energy and generate heat through uncoupling reactions 22 . In our study, both the uncoupling reactions and OXPHOS were down-regulated in low RFI pigs. Uncoupling proteins (UCP) play a critical role in energy-dissipating uncoupling reactions through [H+ ] leakage in mitochondria. It has been reported that overexpression of UCP3 in skeletal muscle increases energy expenditure and decreases feed efficiency in mice 23,24 . In our study, we found UCP2, the ubiquitously expressed isoform of UCP, was down regulated in muscle from low RFI pigs. Also miR-30e, which directly targets UCP2, was up regulated in RFI_L pigs. We also found that genes associated with OXPHOS and ATP synthesis and mitochondria transcriptional control were all down-regulated in low RFI pigs, including COX-I, COX-III, COX-IV, COX-V, PGC-1, PRKAG2, ESRRGB, and HK2. MiR-338 has been confirmed to inhibit COXIV at mRNA and protein levels 25 and miR-338 was up-regulated in low RFI pigs.
PGC-1, also called PPARGC1 or PPARGC1A, plays an important role in mitochondrial biogenesis, by activating cAMP response element binding protein (CREB) and nuclear respiratory factors (NRFs). These nuclear factors then increase the transcription of mitochondrial transcription factors (TFAM, TFB1M, TFB2M) which promote mitochondrial biogenesis 26 . Overexpression of PGC-1 has been reported to increase energy expenditure and mitochondrial number 27 . In our study, PGC-1 and all these down-stream genes related to mitochondrial biogenesis were down-regulated in the low RFI pigs. Also, we found that expression of almost all the genes located in the mitochondrial DNA were down-regulated in low RFI pigs, suggesting either a reduction in the number of mitochondria or a general down-regulation of mitochondrial function. In addition, CREB has been reported to be targeted by miR-335, which was miRNA/Pathway  up-regulated in low RFI pigs 28 . Therefore, these results indicate that mitochondrial biogenesis/function and energy expenditure were reduced in low RFI pigs. AMPK is a key regulator of cellular and whole-body energy balance, which is activated by an increase in the AMP/ATP ratio 29,30 . AMPK can also increase mitochondrial proteins of oxidative metabolism, as well as promote the expression of Hexokinase II (HK2) through CREB in skeletal muscle 31 . In our study, we found that HK2 and Protein kinase, AMP-activated, gamma 2 non-catalytic subunit (PRKAG2), which is a member of AMPK gamma subunit family, were both down regulated in low RFI pigs. Also, miR-144, which has been reported to inhibit the phosphorylation of AMPK alpha (AMPKα), was increased in low RFI pigs 32 . These results indicate that the level of energy metabolism in skeletal muscle was probably reduced in low RFI pigs compared to high RFI pigs.
Further KEGG pathway analyses indicated that PGC-1 and PRKAG2 are involved in the adipocytokine signaling pathway (KEGG hsa04920). In this pathway, mitochondrial fatty acid β -oxidation is increased by the expression of PGC-1, while PRKAG2 increases appetite by up-regulating AGRP and NPY. Recently, Lindholm-Perry found the expression level PRKAG2 in the rumen was associated with ADFI in cattle 33 , although the direction of the association depended upon season. Thus, PGC-1 and PRKAG2 are potentially key genes in the molecular regulation of feed efficiency. Although this pathway is termed the "adipocytokine signaling pathway", the genes are expressed in lots of tissues/ cell types where they have similar functions.
As a PPARα coactivator, PGC-1 plays a key role in the transcriptional control of genes encoding mitochondrial fatty acid β -oxidation (FAO) 34 . It is reported that PPARα/PGC-1 signaling pathway can be activated by chronic overnutrition and obesity, resulting in the up regulation of fatty acid β -oxidation related genes, such as FATP1, FACS1, UCP2 and UCP3 35 . PRKAG2 also plays a critical role in regulating cellular fatty acid metabolic pathways 36 . Thus, the down regulation of these genes indicate that mitochondrial fatty acid β -oxidation was reduced in low RFI pigs compared to high RFI pigs. Hence we propose that the pathway is more to do with regulating whole body energy balance via effects on energy expenditure or appetite, rather than relating to adipose tissue metabolism per se.
A number of miRNAs related to skeletal muscle growth and development were differentially expressed between low and high RFI pigs, including miR-208b, miR-499, miR-29c, miR-1 and miR-99b. The transforming growth factor-beta (TGF-β) signalling pathway, which includes Myostatin (MSTN), is considered the most potent negative regulator of skeletal muscle growth and development 37 . It has been reported that, miR-29 and miR-30b are both inhibitors of TGF-beta [38][39][40] , while miR-208b and miR-499 inhibit the MSTN gene 41 . Also, overexpression of miR-29 promotes myogenic differentiation in C2C12 cells, due to the reduction in TGF-beta, which inhibits differentiation 42 . Furthermore, miR-99b has been reported to be directly stimulated by the TGF-β signaling pathway 43 and miR-99b is reported to inhibit IGF-1/mTOR signaling by targeting AKT1, IGF1R and mTOR 44,45 . In the present study, we found miR-29c, miR-30b, miR-208b and miR-499 were all up-regulated in low RFI pigs, while miR-99b was down-regulated. All these results suggest that muscle growth and development was increased in low RFI pigs via the inhibition of the TGF-β signaling pathway and stimulation of the IGF-1/mTOR signaling pathway. We also found IGF-2 and MYHCs were up-regulated in low RFI pigs. It appears from the pig growth performance data that these changes were associated with a slightly increased ADG, even though DFI was reduced, suggesting that the RFI_L pigs were able to sustain their growth whilst reducing their feed intake.
Thus, based on the entire expression profiles of both mRNAs and miRNAs, we conclude that the improvements in feed efficiency in low RFI pigs are due to inhibition of skeletal muscle mitochondrial activity through PGC-1/TFAM and PRKAG2 combined with the stimulation of muscle growth through TGF-β and IGF-1/2 signaling pathways.

Conclusions
Overall, we identified 99 mRNAs and 25 miRNAs that were differentially expressed in skeletal muscle from pigs with different RFI. These genes were functionally related to metabolism, particularly energy and lipid metabolism, as well as biosynthetic processes and muscle cell growth and differentiation. A number of genes involved in energy metabolism were down-regulated, whereas quite a few miRNAs that target energy metabolism genes were up-regulated in muscle from low RFI pigs. Similarly, a number of genes and miRNAs which stimulate skeletal muscle differentiation and proliferation were up-regulated. We propose that feed efficiency in pigs can be improved by reducing energy metabolism in muscle, particularly mitochondrial metabolism, and/or by enhancing skeletal muscle growth and we identify a number of miRNAs and genes that might be targets for manipulation. This study enhances our understanding of molecular mechanisms regulating feed efficiency in pigs.

Materials and Methods
Animals and tissues. In this study, 238 castrated purebred Yorkshire boars were grown from 30 to 90 Kg, the average period of study was 67.35 days. Pigs were slaughtered at 90 kg according to a standard procedure approved by Biological Studies Animal Care and Use Committee 46 , Hubei Province, P. R. China. They were individually fed ad libitum a complete mixed commercial feedlot ration (see Supplementary Table S8), using ACEMA 64 automated individual feeding systems in the Agricultural Ministry Breeding Swine Quality Supervision Inspecting and Testing Center (Wuhan). All the methods in this study were carried out in accordance with the approved guidelines from Regulation of the Standing Committee of Hubei People's Congress . All experimental protocols were approved by the Ethics Committee of Huazhong Agricultural University.
For RNA sequencing, the 3 pigs with the highest RFI (named RFI_H group) and the 3 pigs with lowest RFI (named RFI_L group) were selected from the 238 pigs, each RFI group having no difference in starting body weight (Table 1). For miRNA sequencing, the 5 pigs with the highest RFI (named RFI_H group) and the 5 pigs with the lowest RFI (named RFI_L group) were selected.
Within 30 minutes after slaughter, a piece of longissimus dorsi muscle of each animal was sampled at the thoracolumbar junction. All tissue samples were immediately frozen in liquid nitrogen and stored at − 80 °C for RNA isolation.
Phenotypes. RFI was calculated by a linear regression model according to the records of daily feed intake (DFI), average daily gain (ADG) and mid-test metabolic body weight (MBW) of all the pigs.
The base model used was Y j = β 0 + β 1 MBW j + β 2 ADG j + e j , where Y j is the DMI of the j th animal, β 0 is the regression intercept, β 1 is the regression coefficient on MBW, β 2 is the regression coefficient on ADG, and e j is the uncontrolled error of the j th animal. Analyses of RNA-Seq data. Having transferred the RNA-Seq results from Illumina fastq format to standard Sanger fastq format with fq_all2std.pl, data were processed with the Tophat-Cufflinks pipeline 47 . The porcine reference genome and gtf annotation file were downloaded from Ensembl (Sscrofa10.2.72) and build index with bowtie version 2.1.0. TOPHAT (version 2.0.9) was used to align reads to the genome with the option --library-type fr-firststrand. Cufflinks (version 2.1.1) was used for transcriptome assembling, and Cuffdiff script from Cufflinks was used for gene expression analysis with the option -classic-fpkm. The expression level of each gene was represented by the FPKM value, which means fragments per kilobase of exon per million fragments mapped, and was calculated by the following formula 48  Analyses of miRNA-Seq data. Firstly, miRNA-Seq data were transferred from Illumina fastq format to fasta format and then the miRNA-seq datasets analyzed with mirDeep (v2.0.0.5) 49,50 . The porcine genome (Sscrofa10.2.72) was downloaded from Ensembl, and the miRNA reference was obtained from the miRBase database (version 20) [51][52][53][54][55] .
The miRNA expression level of each library was normalized by the following formula 56 : Reads count total clean reads count Normalized reads count 1 000 000 2 The P-value between the two libraries was calculated using the following formulas 57 : N1 and N2 represent the total numbers of clean reads in the two small RNA libraries. The q-value was calculated by fdrtool in R 58 . Finally, |log2FC| ≥ 1, total reads count ≥ 20, p ≤ 0.05 was set as the threshold for selection of differentially expressed (DE) miRNA.

Q-RT-PCR validation of differentially expressed genes and miRNAs.
Relative expression levels of the differentially expressed genes and miRNAs in muscle were quantified by real-time PCR. The RPL4 gene and 18S RNA were selected as the internal controls for qRT-PCR validation because of their stable expression in skeletal muscle tissues. The poly(A)-tailed RT-PCR method 59 was performed for miRNA reverse transcription. Primer sequences and PCR conditions for analyzed genes and miRNAs are listed in Supplementary Table S9. The reactions were performed on a Roche Lightcycler 480 Sequence Detection System using SYBR Green PCR Master Mix (TOYOBO, QPK201), following the instruction manual. The 50 μ L reaction mixture consisted of 5 μ L cDNA, 25 μ L 2 × SYBR Green PCR Master Mixture, 0.5 μ M each primers and water. The qPCR profiles began with initial denaturation at 95 °C for 10 min, and then followed by 40 cycles of 95 °C denaturation for 15s, annealed at 60 °C for 15s, and 72 °C for 15s extension. A dissociation curve was generated at the end of the last cycle by collecting the fluorescence data from 58 °C to 95 °C. The 2-∆ ∆ Ct method was employed for relative gene expression level analysis. For each gene, the average ∆ Ct value of the RFI_L group was used as reference to calculate the ∆ ∆ Ct value, and Student's t-test was used to analyze the expression difference between the 2 groups 59 .
Potential target gene prediction of miRNAs. To explore the potential function of miRNAs with significant differential expression in the two groups, potential target genes and pathways of miRNAs were predicted by DIANA miRPath (v.2.0) 60 . As porcine genes were not included in the current version of DIANA miRPath, prediction was performed using human miRNAs. The P-value threshold was 0.05 and MicroT threshold was 0.8 60 (http://www.microrna.gr/miRPathv2).
Gene ontology and pathway analyses. The human homologous Ensembl Gene IDs of the identified DE genes and miRNA target genes were utilized for the following bioinformatics analysis. Gene enrichment in gene ontology (GO) biological processes and pathways were performed with the DAVID Bioinformatics Resources v6.7 (http://david.abcc.ncifcrf.gov/) 61,62 . Cutoff criteria were EASE scores less than 0.01 (GO enrichment) or 0.05 (Gene enrichment in pathways). EASE score were given by DAVID, which is a modified Fisher's exact test. Cytoscape (v3.0.1) was used to create the potential important network(s) of DE genes and miRNAs 63 .