Bovine hepatic miRNAome profiling and differential miRNA expression analyses between beef steers with divergent feed efficiency phenotypes

MicroRNAs (miRNAs) are small RNA molecules involved in regulation of multiple biological processes through modulating expression of their target genes. Here we employed RNAseq to profile liver tissue miRNAome of 60 steers from Angus, Charolais, and Kinsella Composite (KC) populations. Of these animals, 36 animals (n = 12 for each breed) were utilized to identify differentially expressed (DE) miRNAs between animals with high (n = 6) or low (n = 6) phenotypic values of residual feed intake (RFI), a common measurement of feed efficiency. At a threshold of fold-change > 1.5 and P-value < 0.05, we detected 12 (7 up- and 5 downregulated in low-RFI animals), 18 (12 up- and 6 downregulated), and 13 (8 up- and 5 downregulated) DE miRNAs for Angus, Charolais, and KC steers, respectively. Most of the DE miRNAs were breed specific, with bta-miR-449a and bta-miR-AB-2 being differentially expressed in all three breeds. The predicted target genes of the identified DE miRNA are mainly involved in cell cycle, cell death and survival, cell signaling, cellular growth and proliferation, protein trafficking, cell morphology, cell-to-cell signaling and interaction, cellular development, molecular transport, post-translational modification, as well as nutrient metabolism (lipids, carbohydrates, protein and amino acid). Our results provide insights into the bovine hepatic miRNAome and their potential roles in molecular regulation of RFI in beef cattle.

Genetic selection and breeding for more feed efficient beef animals is of great interest to beef producers, since increased feed efficiency will reduce the cost of beef production as feed and feeding related expenditures contribute up to 75% of the total variable production costs 1 . Additionally, studies have shown that breeding for more feed efficient animals can potentially reduce methane emissions from beef cattle 2,3 , which could consequently lower the environmental carbon footprint of beef animals. Understanding the genetic control of complex traits such as feed efficiency can help enhance the rate of genetic improvement of the traits via the designing of more effective genetic or genomic selection tools 4 . In this regard, multiple genome wide DNA marker association studies [5][6][7][8][9] , and global transcriptomic profiling studies [10][11][12][13][14][15][16][17][18] have endeavored to identify markers, genes and biological functions associated with feed efficiency in beef cattle.
MicroRNAs are a group of small RNAs with an average length of approximately 22 nucleotides, mainly biosynthesized through the enzymatic cleavage of longer RNA molecules by DROSHA (nucleic) and DICER (cytoplasmic) RNase endonucleases [19][20][21] . The synthesized miRNA then combines with Argonaute proteins to form a miRNA-induced silencing complex (miRISC) 22 , which under the guidance of the miRNA predominantly binds to the seed regions in the 3′ untranslated regions (UTRs) of the target mRNA, and hence leading to downregulation of the target gene 22,23 . In mammalian cells, miRNAs target and regulate expression of up to 60% of the transcribed genes 24 . Therefore, they are involved in multiple biological functions including cell proliferation, cell cycle, cell development, apoptosis, metabolism of amino acids, metabolism of lipids, metabolism of carbohydrates, and metabolism of minerals and vitamins 25 . In the liver, miRNAs have been implicated in regulating Results miRNA sequence data and alignment quality. On average the Illumina next generation sequencing yielded over 9 M (million) high quality raw reads per sample for Angus and Charolais, and over 11 M for KC samples (Table S1 in Supplementary File 1). After 3′ adaptor clipping, an average of 36.78% of the reads were removed as long reads (> 28 bp), and an average of 8.19% of the reads were removed as short reads (< 15 bp) (Table S1 in Supplementary File 1). Additionally, on average 0.24%, 0.14%, 0.03%, and 0.03% of the reads were removed since they aligned to bovine rRNAs, tRNAs, snRNAs, and snoRNAs, respectively (Table S1 in Supplementary File 1). An average of 5.5 M reads were retained for miRNA profiling analysis by mirDeep2 (Table S1 in Supplementary File 1). Most of the reads ranged between 20 and 24 bp in length as shown in Fig. 1a, with an average length of 21 bp. The retained reads were of high quality as depicted by high average Phred scores in Fig. 1b. Of the retained reads, 74.78% (SD = 2.31%) mapped to the UMD3.1 bovine reference genome on average, ranging from 72.47% (SD = 1.2%) for Charolais to 77.09% (SD = 1.34%) for KC (Table S1 in Supplementary File 1).
Known miRNA expression and novel miRNA profiles. We identified 541, 551, and 575 expressed known miRNAs in Angus, Charolais, and KC samples, respectively. Of all these, 528 unique known miRNAs (approximately 90%) were common to all the three breeds as shown in Fig. 2a. Among the highly expressed miR-NAs, bta-miR-192 was the most abundant miRNA in all the three breeds with an average of 867,342, 1,060,828, and 1,272,798 aligned reads per sample from Angus, Charolais, and KC populations, respectively. Ten miRNAs showed predominantly high expression including bta-miR-192, bta-miR-143, bta-miR-148a, bta-miR-26a, bta-miR-30a-5p, bta-miR-22-3p, bta-miR-27b, bta-let-7f, bta-miR-27a-3p, and bta-miR-101, and they accounted for an average of 78.4%, 78.3%, and 77.9% of the total aligned sequence reads in Angus, Charolais, and KC animals, respectively. The top 20 expressed miRNAs across studied samples from each of the breeds are presented in Table 1, while all the expressed known miRNAs identified, and their average aligned read counts in each breed are presented in Supplementary File 2.
At a mirDeep2 score ≥ 4, an estimated probability that the predicted miRNA candidate is a true positive is greater than 70%, and with a significant Randfold p-value suggesting that the miRNA's precursor sequence could be folded into a thermodynamically stable hairpin, we identified 126 (from 129 precursors), 101 (from 103 precursors), and 119 (from 125 precursors) novel miRNAs in Angus, Charolais, and KC samples, respectively. The identified novel miRNAs were largely expressed in one breed, with only 31 of them being commonly expressed across all the three breeds (Fig. 2b). Of the 241 unique novel miRNAs, bta-miR-AB-10 and bta-miR-AB-9 were the most expressed ones across the three breeds with an average of 52,817 and 46,932 reads aligned to these miRNAs, respectively ( Table 2). The hairpin structure of the precursor and the read alignment distribution (i.e. alignment to the mature, star and loop sequences) of bta-miR-AB-10 across the three breeds are presented in Fig. S1 in Supplementary File 1, indicating that the majority of the reads were aligned to mature micro RNA portion of the precursor. The top 20 expressed novel miRNAs and their miRDeep2 prediction scores are presented in Table 2 for each breed population, while all the identified novel miRNAs and their miRDeep2 prediction scores are provided in Supplementary File 3. miRNA differential expression. Differential miRNA expression analysis was performed between the six (n = 6) low-and six (n = 6) high-RFI steer groups from the profiled animals for each breed population as selec-  www.nature.com/scientificreports/ tion of these 6 low-and 6 high-RFI resulted in an increased divergence between the low and high RFI animal groups. The low-and high-RFI steer groups of the three breed populations were all significantly different in their average RFI phenotypic values at P-value < 0.0042 after the Bonferroni correction for multiple comparisons ( Table S2 in Supplementary File 1). For other traits, low-RFI animals consumed significantly less feed per day than their high-RFI counterparts in Charolais and KC. In Angus, low-RFI animals on average also consumed less feed as compared to the high-RFI animals although the difference did not reach the level of statistical significance (i.e. P-value > 0.0042). The average phenotypic values of other traits including animal slaughter age were not significantly different between the high and low RFI groups (Table S1 in Supplementary File 1). At a fold change (FC) > 1.5 and a P-value < 0.05, we identified 12 DE miRNAs in the liver tissue of Angus including 10 known miRNAs and two novel miRNAs (Table 3). Of these DE miRNAs, five were downregulated and seven were upregulated in low-RFI animals. In Charolais, we identified 18 DE miRNAs including 16 known miRNAs and two novel miRNAs, of which six were downregulated and 12 were upregulated in low-RFI Charolais steers (Table 3). In KC, 13 DE miRNAs including 10 known miRNAs and three novel miRNAs were identified, with five downregulated and eight upregulated in low-RFI steers (Table 3). Of all the identified DE miRNAs, bta-miR-449a and bta-miRNA-2 were common to all three breeds (Fig. 3). Of the two common DE miRNAs, bta-miR-449a was upregulated in the low-RFI steers of all three breed populations whereas bta-miRNA-AB-2 was upregulated in Charolais, and downregulated in both Angus and KC low-RFI steers (Fig. 3).
Target gene prediction and functional enrichment analyses for the most abundant known and novel miRNAs. We performed target gene prediction for 18 miRNAs (16 known and 2 novel miRNAs bta-miR-AB-10 and bta-miR-AB-9) that were identified as the most expressed across the three breed populations for functional enrichment analyses (Tables 1, 2). At a threshold cut off of not less than the 99th context++ score percentile, we identified 1022 target genes. A list of all these identified targets, their respective TargetScan scores for the fourteen parameters, and the mRNA-miRNA correlation coefficients in the samples of the three populations are provided in Supplementary File 4. Of these target genes, 1008 mapped to the IPA ingenuity database, and they are mainly involved in cell morphology, cellular assembly and organization, cellular growth and proliferation, free radical scavenging, and cell death and survival biological functions. All metabolic and cellular functions significantly enriched by the identified targets are presented in Fig. S2 of Supplementary File 1. For cell morphology, the target genes are mainly involved in maintaining morphology of cellular organelles such as nucleus and cellular the transmembrane potential. For cellular growth and proliferation, the target genes are involved in proliferation of liver cells, and proliferation and development of different immune cells such as macrophages, lymphocytes, and natural killer cells, whereas for cell death and survival biological function, the target genes are primarily involved in necrosis and apoptosis. Also, IPA identified sirtuin signaling, 1L-12 signaling and production in macrophages, tryptophan degradation III, receptor recognition of bacteria and viruses and senescence pathway, as presented in Fig. S3 in Supplementary File 1.

DE miRNAs target gene prediction.
For the differentially expressed known and novel miRNAs within each breed (Table 3), we identified 767, 1667, and 787 target genes for Angus, Charolais, and KC, respectively, at a minimum context++ score percentile of 99. The targets, their targeting DE miRNAs, the target-miRNA TargetScan estimated scores and the mRNA-miRNA correlation coefficients in the studied animals from each breed are presented in Supplementary File 5. Functional enrichment analysis for the target genes identified for DE miRNAs in Angus are mainly involved in regulating the gene expression of cell cycle, gene expression modulation, cellular assembly and organization, DNA replication, recombination and repair and RNA post-transcription modification (Fig. S4 in Supplementary File 1). Additionally, the major pathways enriched by these target genes included HIPPO signaling, GADD45 signaling, ATM signaling, endoplasmic reticulum stress pathway, and Granulocyte adhesion and diapedesis pathway as shown in Fig. S5 in Supplementary File 1. For Charolais, the identified targets are majorly involved in protein synthesis, protein trafficking, RNA post-translation modification, lipid metabolism, and molecular transport (Fig. S6 in Supplementary File 1). Estrogen receptor signaling, granzyme A signaling, sirtuin signaling pathway, protein ubiquitination pathway, and coagulation system  Table 3. Differentially expressed known and novel micro RNAs (miRNAs) between high and low-RFI animals within each breed population [Angus, Charolais, and Kinsella Composite (KC)], with a differential expression threshold of P-value < 0.05 and Fold-change (FC) > 1. 5. The sign of logFC shows the direction of miRNA expression in low-RFI steers relative to high-RFI animals. The miRNAs that were significantly differentially expressed at a more stringent threshold of Bonferroni false discovery rate (FDR) correction < 0.05 and Foldchange > 1.5 are highlighted in bold and a "**" was added to their P-values. logCPM = log 2 (counts per million), LR = likelihood ratio test static value, and logFC = log 2 (Fold-change). ). The miRNA bta-miR-449a, which was identified as a common and upregulated DE in low-RFI steers of all the three breeds, was predicted to target 16, 11, and 35 DE genes in Angus, Charolais, and KC steers, respectively, of which only three target genes including (SERPINA3, TP53INP1 and LPIN1) were common to all the three breeds. Of the 16 target genes of bta-miR-449a identified in Angus, 12 and 4 were downregulated and upregulated respectively in the liver tissue of low-RFI animals (Fig. 4). Of the 11 targets identified for bta-miR-449a in Charolais, 4 and 7 were downregulated and upregulated respectively in the liver tissue of low-RFI animals (Fig. 5). In KC, of the 35 targets by bta-miR-449a 23 and 12 were downregulated and upregulated, respectively, in the liver tissue of low-RFI animals (Fig. 6). In Angus, RAB30 and FKBP5 genes were found to be the major target genes for the identified DE miRNAs, with each predicted to be targeted by six miRNAs (Fig. 4). The other six major target genes including COL1A1, ELOVL5, SCD, TLE1, TP53INP1, and TTC39C, were predicted to be regulated by five DE miRNAs each, as shown in Fig. 4. Of these targets, FKBP5, COL1A1, SCD, TLE1, and TP53INP1 were downregulated in low-RFI animals, whereas RAB30, ELOVL5, and TTC39C were upregulated in the same animals. Additionally, the miRNAs and DE target gene interaction networks showed bta-miR-2285n, bta-miR-2285u, and bta-miR-449a (all upregulated) as the major miRNAs, targeting 23, 19, and 16 DE genes, respectively (Fig. 4). Other major miRNAs including bta-miR-AB-2, bta-miR-AB-47, and bta-miR-424-3p, are predicted to regulate 14, 13, and 12 DE genes, respectively. We also identified that 39 DE genes were targeted by one or multiple upregulated (in low-RFI animals) DE miRNAs in Angus, of which 23 were downregulated in low-RFI animals whereas 16 were upregulated in the animals of low-RFI.
For Charolais, the major targeted genes included SIRPA (predicted to be targeted by 10 DE miRNAs), ABCC4 (predicted to be targeted by nine DE miRNAs), DLK1, and TP53INP1 (each predicted to be targeted by eight DE miRNAs), SCD, SLC7A5, and THEM4 (each predicted to be targeted by seven DE miRNAs), and AK4 (predicted to be regulated by six DE miRNAs) (Fig. 5). SLC7A5, TP53INP1, SCD, and THEM4 were downregulated, www.nature.com/scientificreports/ whereas AK4, SIRPA, DLK1, and ABCC4 were upregulated in low-RFI animals. Furthermore, bta-miR-2285ai-5p (14 targets), bta-miR-7859 (12 targets), bta-miR-2284ac (11 targets), bta-miR-AB-2 (11 targets), bta-miR-449a (10 targets) were predicted to be major regulators among the identified DE miRNAs and were all upregulated in liver tissue of low-RFI steers, as shown in Fig. 5. SIRPA (upregulated in low-RFI steers) was predicted to be co-regulated by all the five major DE miRNAs.  (Fig. 6). Of the 106 DE genes predicted to be targeted by a single or multiple upregulated (in low-RFI steers) DE miRNAs, 62 were upregulated in the liver tissue of the same animals. However, the other 44 DE genes were upregulated despite being targeted by either a single or multiple upregulated miRNAs (Fig. 6).
Functional enrichment analysis of the DE miRNA targeted DE genes identified multiple biological functions that greatly overlapped with those identified for all the DE genes reported by Mukiibi et al. 18 . For Angus, the DE miRNA target DE genes are mainly involved in lipid metabolism, molecular transport, small molecule biochemistry, energy production, and carbohydrate metabolism, as shown in Fig. S10 in Supplementary File 1. For Charolais, lipid metabolism, molecular transport, small molecule biochemistry, cellular movement, and cell-to-cell signaling and interaction were identified as the main biological functions enriched by the identified DE targets, as shown in Fig. S11 in Supplementary File 1. For KC, the DE genes targeted by the DE miRNA are primarily involved in amino acid metabolism, small molecule biochemistry, lipid metabolism, vitamin and mineral metabolism, and cell death and survival, as presented in Fig. S12 in Supplementary File 1. All the target DE genes involved in each of the top five significantly enriched cellular and metabolic functions for the three breeds are presented in Table S3 in Supplementary File 1.

Discussion
RNA sequencing offers a greater resolution to profile miRNAs, even at very low levels of expression in the cells 31 , and allows for the parallel profiling of abundance of known miRNAs as well as the identification of novel miRNAs 32 . Additionally, the availability of miRNA sequences enables the prediction of potential target genes for both known and novel miRNAs 31 . In the current study we employed deep sequencing of small RNAs to profile miRNA expression in the liver tissues of 60 beef steers from three distinct beef breed populations. We obtained high quality sequence data as revealed by our sequencing quality results, with an average Phred quality score of www.nature.com/scientificreports/ 37.3 across the samples of the three breeds. Additionally, sequence data quality control processing (including removal of other small miRNAs) resulted in read sequences with an average length of 21 bp, and the majority of the reads ranged between 20 and 24 bp across the samples from the three breeds, as shown in Fig. 1a, providing high quality and reliable reads for downstream alignment and abundancy profiling of miRNAs 33 . Furthermore, to validate the RNAseq miRNA profiling results, we conducted TaqMan qPCR differential expression analysis on six selected miRNAs (Table S4 in Supplementary File 1). It showed that all expression levels in low-RFI animals were consistent between the RNAseq and TaqMan qPCR ( Fig. S13 and Table S4 in Supplementary File 1), with the P-values of significance tests of qPCR between high and low-RFI animals groups ranging from a significant P-value (< 0.05) for bta-miR-133a (P-value = 0.003) to suggestive P-values (P-value = 0.054 and P-value = 0.086) for bta-miR-223 and bta-miR-424-5p, respectively, to relatively higher P-values ranging from 0.144 to 0.314 for the other three miRNAs, as shown in Fig. S13 and Table S4 in Supplementary File 1. Although the differential expression P-values for qPCR results were not as low as those for RNAseq, which might be related to sensitivity differences between the two methods regarding profiling gene expression 31 , we observed that the Pearson correlation coefficients were high (0.73 -0.99) between the qPCR and RNAseq profiles and significant (6.79E−03 to 6.00E−10) ( Table S5 in the Supplementary File 1), which indicates that the RNAseq results were highly reproducible.
With a high average alignment rate close to 75%, we found that approximately 90% of the identified miRNAs were expressed in all the three breed populations, indicating a high level of similarity among the breeds in terms of hepatic miRNA expression. Similarly, high similarity rates in expression patterns were also observed for the protein coding genes in the same populations, where over 96% of the expressed genes were commonly expressed among the liver tissues of the same beef cattle populations 18 . Of the expressed miRNAs, 10 were highly expressed including bta-miR-192, bta-miR-143, bta-miR-148a, bta-miR-26a, bta-miR-30a-5p, bta-miR-22-3p, bta-miR-27b, bta-let-7f, bta-miR-27a-3p, and bta-miR-101, and these accounted for an average of 78.2% of the aligned reads for the profiled liver tissues of the studied animals across the three breeds ( Table 1). The miRNAs bta-miR-101, bta-miR-143, bta-miR-30a-5p, bta-let-7f, bta-miR-192, and bta-miR-148a were previously reported among the 10  35 . In sheep, miR-192 has been reported to be involved in regulating the growth and development of skeletal muscle 36 .
To explore the biological importance of the highly expressed miRNAs across the liver tissues of the studied animals, we performed target gene prediction of 18 miRNAs that were highly expressed in steers from the three breeds. Functional enrichment analysis of the target genes revealed that the candidate target genes were involved in key biological processes including maintaining cellular homeostasis, immune functions, proliferation of liver cells, and apoptosis of cells. Consistent with our results, some of the most abundant miRNAs have been identified as important modulators of liver cellular metabolic homeostasis, liver cell proliferation and development, liver cell death and regeneration in different species 26 . For example, miR-143 which was the second most abundant miRNA in our studied samples was reported to be involved in glucose and insulin metabolism in mice 37 . miR-148 and miR-26a are involved in the regulation of mice hepatocyte proliferation 38,39 , a key process in liver tissue regeneration. In human and mouse liver tissue, miR-148a 40 and miR-27 41 , respectively, have been identified as regulators of liver detoxification. Based on our results and the conserved nature of the miRNA-mRNA interactions across www.nature.com/scientificreports/ mammalian species 24 , we speculate that these highly expressed miRNAs in the bovine liver might play similar functions as those highlighted in other species. Furthermore, the involvement of these miRNAs in proliferation of different cells as well as cell death and survival is plausible, since the liver is in a constant state of self-regeneration or regrowth to recover hepatic tissue lost due to assault by pathogens, toxins, and exogenous antigens 42,43 . Liver regeneration is a complex and highly regulated process that includes the initiation phase, the cell proliferation phase, and the regeneration termination phase, all of which are likely modulated by miRNAs 44,45 . However, further molecular studies beyond Insilco predictions are needed to precisely validate the target genes of these miRNAs and their respective functions in cattle given the physiological differences between the species. The liver is a central metabolic organ serving major biological functions in the mammalian body including nutrient (lipids, carbohydrates, proteins/amino acid, and vitamins and minerals) metabolism, xenobiotics and toxin metabolism, in addition to pathogen processing and growth regulation 46,47 . MicroRNAs are known to modulate all these functions 26 . Therefore, differential hepatic miRNA expression between efficient and inefficient animals can potentially contribute to the molecular variability for feed efficiency in beef cattle. In the current study we identified 39 differentially expressed (DE) known and novel miRNAs between the efficient and inefficient steer groups in all the three studied populations. However, most of the identified DE miRNAs (94.9%) were breed specific. These findings were in line with the results in our mRNA differential expression studies 18,27 that involved the same animals, where most of the DE genes (from 63.4% for Angus to 84.6% for KC) associated with RFI 18 and its component traits 27 were breed specific as well, which we speculate may be due to genetic uniqueness of the studied breeds. However, since we used a lower threshold or raw P-value < 0.05 to identify the differentially expressed miRNAs, further studies involving larger beef cattle populations and more divergent phenotype animal groups would be required to validate the current results.
Within each breed, most (i.e. 58.3% for Angus, 66.7% for Charolais, and 61.5% for KC) of the differentially expressed miRNAs were upregulated in low-RFI animals, and hence suggests a general expectation of reduced expression of their target genes. Compared with previous studies that have investigated the association of miRNA expression with feed efficiency 28,29 , bta-miR-424-5p was the only miRNA that was commonly DE in our study and in the Austrian Angus population 28 . However, in our study, bta-miR-424-5p was upregulated in the liver tissue of low-RFI animals (KC), whereas it was identified as downregulated in the liver tissue of low-RFI Angus bulls 28 . This difference could be due to genetic, physiological or environmental differences between the animals studied in the two studies. It is also worth noting that Al-Hussein et al. sequenced two cDNA libraries of pooled RNA from efficient and inefficient bulls 28 , whereas in our study we sequenced each individual cDNA library for each of the 60 studied steers, and this could also contribute to the differences between the two studies.
To investigate the potential biological importance of the DE miRNAs related to RFI within each breed, we predicted their target genes at a TargetScan context++ score percentile threshold of 99 for each breed. The key biological functions associated with the target genes included cell cycle, cellular growth and proliferation, cellular assembly and organization, lipid metabolism, protein synthesis, protein trafficking, cell-to-cell signaling and interaction, molecular transport, post-translational modification, and RNA post-transcription and translation modification. Of these functions, cell signaling, cellular growth and proliferation, lipid metabolism, protein synthesis, cellular development, cell death and survival molecular transport and protein synthesis have been previously reported to be associated with feed efficiency in beef cattle through liver transcriptomic studies 10,11,15,18 . To further investigate the interactions between the identified DE miRNAs and the 72, 41, and 175 differentially expressed (DE) genes that were previously identified in the same liver tissues of the same three animal populations 18 , we identified only 1 (GNAZ), 2 (THEM4 and CES1), and 5 (PALMD, C12orf45, IRAK3, CITED4, and IL20RA) out of the 72, 41, and 175 DE genes as targeted DE genes for the DE miRNAs, respectively, in the Angus, Charolais, and KC breed populations using a TargetScan context++ score percentile threshold of 99 for each breed populaton. The low number of targeted DE genes by the DE miRNAs for each breed population would not allow for meaningful gene network and functional analysis. Thereore, we further used a lower threshold of context++ score percentile > 50 to predict targeted DE genes out of the 72, 41, and 175 DE genes of each breed population since a context++ score percentile > 50 is also considered as a high context score to shows that a specific site is more favorable than most other sites of the miRNA for target gene identifiction 48 . With a context++ score percentile > 50, it was revealed that 61.1% (44 DE genes), 75.6% (31 DE genes) and 73.7% of the previously identified DE genes were predicted as potential targets of the DE miRNAs in Angus, Charolais, and KC steers, respectively. These target DE genes are mainly involved in lipid metabolism, molecular transport, small molecule biochemistry, energy production, carbohydrate metabolism, cellular movement, cell-to-cell signaling and interaction, cell death and survival, amino acid metabolism and vitamin and mineral metabolism, implying that these identified DE miRNAs influence feed efficiency through differential modulation of these biological functions in the liver. From the DE miRNA-DE mRNA interaction networks in Figs. 4, 5 and 6, it can be observed that some miRNAs were predicted to target multiple DE genes, and some single genes were predicted to be targets for multiple miRNAs. These complex forms of miRNA-mRNA interactions emanate from the fact that a single miRNA using its seed region can bind to multiple sites in the 3′-UTRs of different genes (mRNAs), and also one target can have multiple binding sites for several miRNAs 23,49 , hence allowing miRNAs to modulate multiple biological processes even though they are small in number compared with the mRNAs that they regulate.
To a large extent, we observed contrasting expression patterns in DE miRNAs compared with DE genes in the liver tissues as expected since the main mode of action of miRNAs is through promoting deadenylation of their target transcripts which conquently results into accelerated degradation of those target mRNAs 50 . As an example, of the identified DE target genes for the upregulated miRNAs bta-miR-2285n, bta-miR-2285u, bta-miR-449a and bta-miR-47 in Angus low-RFI animals, 52%, 63%, 75% and 62% were downregulated, respectively. For Charolais, of the predicted DE targets of upregulated miRNAs bta-miR-2285ai-5p, bta-miR-7859, bta-miR-2284ac and bta-miR-449a in low-RFI animals, nearly 50%, 50%, 64% and 36%, respectively, were downregulated. In KC, of the predicted DE targets of upregulated DE miRNAs bta-miR-424-5p, bta-miR-223, bta-miR- www.nature.com/scientificreports/ and bta-miR-363 in low RFI animals, approximately 62%, 67%, 66% and 54%, respectively, were downregulated. However, we also observed a significant number of upregulated predicted target genes despite being targeted by multiple upregulated miRNAs. These observations could be attributed to the different mechanisms of miRNA gene regulation in mammalian cells including the augmentation of mRNA degradation through deadenylation and translation (proteins) repression when the miRNAs bind to the 3′ UTRs of their targets 50 . Therefore, the observation that some DE genes predicted to be targeted by upregulated DE miRNAs were of downregulated in the liver tissue of low-RFI animals could due to the augmented mRNA degradation mode of regulation. On the other hand, protein translation repression could be the main miRNA gene expression modulation mode for those genes that retained high expression profiles despite being targeted by multiple upregulated DE miRNAs.
To better understand the relationship of miRNA and mRNA, the correlation coefficients between miRNA and mRNA of the target genes were also investigated. In general the negative correlation coefficients between the DE miRNAs and their DE target genes largely confirmed the DE miNRA-target mRNA regulatory interaction patterns in Figs. 4, 5, and 6, while the magnitudes of correlation coefficients shed some light on strength of their DE miNRA-target mRNA regulatory interactions. For Angus, miRNA bta-miR-2285n has the strongest negative correlation with the target gene SCD (r = − 0.49, Supplementary File 6), which indicates that increased expression of the miRNA would reduce SCD mRNA to a greater extent. Similar strong negative correlations were also found between bta-miR-2346 and target gene CYP2C19 for Charolais, and between bta-miR-2483-5p and target gene ZBTB16 for KC. However, strong positive correlations were also detected between certain miRNA and target DE genes in the three breeds (Supplementary File 6), which suggests that these targeted genes are under different mechanisms of miRNA gene regulation to a greater extent. Although we used a context++ percentile score from greater than 50 to 100 as suggested by the Targetscan 7 documentation 48 to predict targeted genes of miRNAs in this study, the interactions of the DE miNRA and target mRNAs should be interpreted with caution due to the nature of Insilco prediction that is subject to prediction accuracy. Hence, further molecular experiments are warranted to validate these interactions.

Conclusions
In the current study we employed RNAseq to perform hepatic miRNAome profiling of beef steers from Angus, Charolais, and KC populations. We identified a total of 588 expressed known bovine miRNAs, of which 89.8% were commonly expressed in the liver tissue of the animals from the three populations. Of these miRNAs, bta-miR-192, bta-miR-143, bta-miR-148a, bta-miR-26a, bta-miR-30a-5p, bta-miR-22-3p, bta-miR-27b, bta-let-7f., bta-miR-27a-3p and bta-miR-101 had the greatest expression levels in all three breeds. We also detected 241 novel bovine miRNAs expressed in the liver tissue, with 69.3% identified as expressed in only one of the three breeds whereas 12.9% were identified as expressed in all of the three populations. Differential miRNA expression analyses found 39 miRNAs that were associated with RFI, including five novel miRNAs (bta-miR-AB-2, bta-miR-AB-47, bta-miR-AB-15, bta-miR-AB-63, and bta-miR-AB-225). Most of the DE miRNAs were breed specific with only two miRNAs (bta-miR-449a and bta-miR-AB-2) being differentially expressed in all three breed populations. The predicted target genes of the identified DE miRNAs are involved in multiple cellular and molecular functions including cell cycle, cell death and survival, cell signaling, cellular growth and proliferation, lipid metabolism, protein synthesis, protein trafficking, cell morphology, cell-to-cell signaling and interaction, cellular development, molecular transport, post-translational modification, amino acid and carbohydrate metabolism. Additionally, the identified DE miRNAs were found to target approximately 70% of the previously identified DE genes in the liver tissue of feed efficient and inefficient animals from the same populations. These target DE genes are mainly involved in lipid metabolism, molecular transport, small molecule biochemistry, energy production, carbohydrate metabolism, cellular movement, cell-to-cell signaling and interaction, cell death and survival, amino acid metabolism, and vitamin and mineral metabolism. Our results provide further insight into hepatic miRNAome expression profiles of beef cattle and their potential molecular regulatory mechanisms relating to feed efficiency in beef cattle.

Materials and methods
Animal populations and management. The population description and management practices of the experimental animals used in the current study have been presented in our recent studies 18,27 . Briefly, the animals were raised and managed following the Canadian Council of Animal Care (CCAC) guidelines on the care and use of farm animals in research teaching and testing 51 , and all the experimental procedures applied to the animals were approved by the University of Alberta Animal Care and Use Committee (AUP00000777). The animal populations included a total of 256 beef steers from two purebred populations (i.e. Angus and Charolais), and one Kinsella Composite (KC) beef steer population. The steers were born, raised and managed similarly at the Roy Berg Kinsella Ranch, University of Alberta, Canada. All the purebred Angus and Charolais cows were serviced through estrous synchronization and artificial insemination, followed by natural service by purebred Angus or Charolais bulls whose pedigree records were maintained by the Canadian Angus or Charolais Association, respectively. The KC animals were produced through crossing Angus, Charolais, or Alberta Hybrid bulls with the University of Alberta's hybrid dam line. The crossbreeding design used to generate the University of Alberta's hybrid dam line was previously described by Goonewardene et al. 52 . Additionally, commercial crossbred bulls have been added to the KC herd since 2012 for natural service. The animals used in the current study were born between the months of April and May of 2014 and were castrated immediately after birth. The steer calves were maintained with their dams on pasture and weaned at an average age of about six months. The weaned animals were transitioned to a backgrounding diet composed of 80% barley silage, 17% barley grain, and 3% rumensin pellet supplement, and thereafter were fed set-up diets with gradually decreasing barley silage and increasing barley grain proportions for 3 weeks. Subsequently, the animals were introduced to the finishing Scientific Reports | (2020) 10:19309 | https://doi.org/10.1038/s41598-020-73885-5 www.nature.com/scientificreports/ diet of 75% barley grain, 20% barley silage, and 5% rumensin pellet supplement (as fed basis) after the transition period of 3 weeks. The nutrient composition of the finishing diet is provided in the

GrowSafe feedlot test, phenotype measurement and RFI calculations. Between the months of
April and August in 2015, 50 Angus, 48 Charolais, and 158 KC steers were measured for individual feed intake using the GrowSafe System (GrowSafe Systems Ltd., Airdrie, Alberta, Canada). A detailed description on measuring each individual animal's daily feed intake using the GrowSafe automated system has been described by Mao et al. 53 and also provided in our previous reports 18,27 . Briefly, the animals were tested for feed intake for a period ranging from 70 to 73 days, during which animals were fed on the finishing diet described above. Daily dry matter intake (DMI) of each animal was calculated as the average of the daily feed intake records over the test period and was standardized to 12 MJ ME per kg dry matter based on the energy content of the diet. Initial body weight (BW) and average daily gain (ADG) for each animal were obtained from a linear regression between serial body weight measurements and time (days) that were recorded on two consecutive days at the start of feedlot test, at approximately 14-day intervals during the feedlot test, and on two consecutive days prior to the end of the test. Metabolic mid-weight (MWT) was calculated as midpoint BW 0.75 , where midpoint BW was computed as the sum of the initial BW of the animal and the product of its ADG multiplied by half the number of days under the feedlot test. RFI was calculated as the difference between the actual standardized daily DMI and the expected DMI that was predicted using the regression intercept and regression coefficients of ADG and MWT on actual standardized daily DMI as proposed by Koch et al. 54 .
Liver tissue collection. Tissue collection and processing procedures were previously described in our recent studies 18,27 . Briefly, all animals used in the current study were slaughtered at the Agriculture and Agri-Food Canada (AAFC) Lacombe Research Centre (Lacombe, AB) between July and September of 2015. Animals were considered ready for slaughter at an average back-fat thickness of 8 mm between the 12th and 13th ribs, which was measured using an Aloka 500 diagnostic real time ultrasound machine with a 17 cm 3.5 MHz linear array transducer. For slaughter, animals were first stunned by captive bolt and then exsanguinated. The animals were on average slaughtered at the age of 494 ± 3, 518 ± 4, and 457 ± 4 days for Angus, Charolais, and KC, respectively. The liver of each animal was collected immediately after slaughter and dissected at relatively the same location on the right lobe, and the fibrous capsule was removed from the sliced liver tissue samples. The sample tissues were further dissected into smaller portions that were bagged separately in plastic re-closable bags, labelled, flash frozen in liquid nitrogen, and transported on dry ice to the laboratory within 6 h, where they were stored at − 80 °C until total RNA extraction.
Total RNA extraction. Total RNA extraction was performed on 20 of the liver tissue samples (10 with the highest and 10 with the lowest RFI phenotype values) from each breed as described in our recent studies 18,27 . Each of the selected liver tissue samples were pulverized into a fine powder using liquid nitrogen and a prechilled mortar and pestle on dry ice. www.nature.com/scientificreports/ Bioinformatic sequence data processing and miRNA expression profiling. Raw sequence reads were initially assessed for sequencing quality using FASTQC version 0.11.7 55 as described in our recent studies 18,27 . The reads were evaluated for quality based on multiple parameters such as average read length, adaptor content, per sequence GC content, and per base sequence quality scores. Thereafter, the Illumina 3′ adaptor sequence (TGG AAT TCT CGG GTG CCA AGG) was clipped from all the raw read sequences using cutadapt version 1.16 56 . In this study, reads of lengths shorter than 15 bp, or longer than 28 bp were removed as short or long reads, respectively. The retained reads were filtered for other bovine short RNA species including ribosomal RNAs (rRNAs), transfer RNAs (tRNAs), small nuclear RNAs (snRNAs), and small nucleolar RNAs (snoRNAs), downloaded from the Rfam database via RNAcentral a non-coding RNA database (https ://rnace ntral .org/searc h?q=exper t_db:%22Rfa m%22, accessed July 2018). The final processed sequence reads were re-evaluated for quality using FASTQC version 0.11.7 55 . To profile both novel and known miRNA expression in the samples from the cleaned sequence data, the miR-Deep2 package version 2.0.0.8 modules 33 were used together with the UMD3.1 bovine genome from Ensembl version 93, the known bovine mature miRNA sequences and their precursor sequences from the miRBase database (release 22) 57 . The mapper module (mapper.pl) with default parameters was used to collapse reads of the sequences into clusters, and then bowtie-1.1.1 short sequence aligner 58 was employed to align the collapsed reads to the indexed UMD3.1 reference genome. Using default parameters and input files including the reference genome, collapsed reads versus reference genome alignment, known bovine (and human) mature miRNAs and their precursor sequences (including the hairpin structures), and Bos taurus (bta) as the species of interest, all known bovine miRNAs were quantified by the miRDeep2 module (miRDeep2.pl), hence producing read counts for each individual sample.
Subsequently, miRDeep2 was used to predict possible novel miRNAs and their respective precursors based on their read alignment to the bovine reference genome 33 . Genomic regions stacked with aligned reads were excised as potential precursors and evaluated by the RNAfold tool within ViennaRNA-1.8.4 59 for their potential to form stable secondary structures (hairpins), their ability to be partitioned into mature, loop and star strand, and their base pairing in the mature miRNA region. Subsequently, for each predicted miRNA, a mature miRNA consensus sequence, a precursor sequence, RNAfold P-value, the miRDeep2 score, and the probability that the predicted miRNA was a true positive were estimated and produced as outputs. Additionally, for each predicted novel miRNA, mirDeep2 outputs their hairpin structure of the precursor sequence and the read counts of each sample.
Differential miRNA expression analysis. Initially, counts for each mature miRNA from more than one precursor were averaged. Thereafter, all miRNAs that had less than 10 total read counts across the studied samples within each breed were filtered out. Then miRNA expression variation patterns between 20 samples in each breed were visualized through a principle component analysis of the read counts from the mirDeep2 module using the DESeq2 Bioconductor package 60 and the ggplot2 R 61 package. We performed differential miRNA expression for all expressed known mature miRNAs and top 25 expressed novel miRNAs in each breed. Twelve samples including six samples with extreme high and six samples with extreme low-RFI phenotypes that showed consistent miRNA expression (as compared to all other samples in the same breed) were considered for differential miRNA expression using the egdeR Bioconductor package in R 62 . To reduce false positive rates of the analyses, miRNAs within samples from each of the breeds that had less than one count per million (CPM) in at least six samples (half of the analyzed samples) were filtered out from the analyses, as proposed by Anders et al. 63 . For the retained miRNAs, their counts were normalized using the TMM method 64 . To test for differential miRNA expression between high and low-RFI steer groups from each breed using the egdeR package in R 62 , the normalized counts were modeled using a generalized linear model under a negative binomial distribution with the high-RFI group as a reference. MicroRNAs were deemed differentially expressed (DE) at a P-value less than 0.05, and fold change (FC) greater than 1.5.

Validation of differentially expressed miRNAs. Six differentially expressed miRNAs with relatively
high expression levels per sample were selected for validation of the small RNAseq results. These included bta-miR-2415-3p, bta-miR-133a, and bta-miR-2419-5p for Charolais, and bta-miR-424-5p, bta-miR-223, and bta-miR-155 for KC. bta-miR-192 and bta-miR-93 were selected as endogenous controls for Charolais whereas bta-miR-2284 × and bta-let-7b were selected as reference miRNA genes for KC based on their expression abundance and stability across samples (average M values of 0.18 for KC and 0.22 for Charolais), which were determined using geNorm in the GenEx Software v.5.2.7.44 (MultiD Analyses AB, Gothenburg, Sweden). Relative expressions of the selected miRNAs were obtained through stem-loop RT-TaqMan qPCR 65 from the same total RNA as was used for the small RNA sequencing. Reverse-transcription (RT) stem-loop primers and TaqMan qPCR assays (containing the probe and forward and reverse primers) were purchased from Thermo Fisher Scientific (https ://www.therm ofish er.com). RT primer IDs and TaqMan qPCR assay IDs for each validated internal control are provided in Table S7 of Supplementary File 1. Serial dilutions of pooled cDNA samples were used to determine amplification efficiencies using the equation E = − 1 + 10 (−1/slope) , and the slope was calculated by plotting the linear curve of cycle threshold (C T ) values against the log dilutions 66 . Primers had PCR efficiencies between 89 and 110%.
The reverse transcription reactions for each sample including no RNA template controls were performed using the TaqMan  Threshold cycle (C T ) values from the Real-Time PCR thermocycler were then imported into GenEx Software v.5.2.7.44 (MultiD Analyses AB, Göteborg, Sweden). The C T values were adjusted to account for inter-plate variation using the inter-plate calibrator sample included on the plates and to account for amplification efficiencies. The adjusted C T values of the replicates were averaged, and then normalized to the endogenous controls of miRNA's (reference genes) expression. Lastly Log 2 relative quantities were calculated to the average (for all test samples) C T value. The relative quantities were then analyzed for differential miRNA expression between the high and low-RFI steers using a two tailed t-test.
miRNA target genes prediction and functional enrichment analyses. Target genes prediction was performed using three TargetScan version 7.0 Perl scripts 48 downloaded from https ://www.targe tscan .org/ cgi-bin/targe tscan /data_downl oad.vert7 2.cgi (Accessed July 2018) for the top expressed known miRNAs and top most expressed novel miRNAs that were commonly expressed across the three breeds. TargetScan predicts miRNA target genes based on a quantitative model that scores candidate target genes based on 14 features including 3′-UTR target-site abundance, predicted seed-pairing stability, identity of the nucleotide at position 1 of the small RNA (sRNA), identity of the nucleotide at the 8th position of the sRNA, identity of the nucleotide at the 8th position of the target site, local AU content near the target site, supplementary pairing at the miRNA 3′ end, predicted structural accessibility, minimum distance of the site from the stop codon or polyadenylation site, probability of target site conservation, open reading frame (ORF) length, 3′-UTR length, number of offset-6mer sites in the 3′ UTR, and the number of 8mer sites in the ORF 48 . Firstly, we predicted the conserved and non-conserved target sites using the targetscan_70.pl script by providing all the known gene transcripts UTR sequence alignments and the miRNA family information as inputs. The miRNA family information file included the miRNA family IDs, the seed sequence (from 2nd position nucleotide to 8th position nucleotide 5′ of miRNA sequences), and the NCBI IDs of 8 species (cow, sheep, domestic goat, horse, human, mouse, rat and pig). Thereafter, we used the targetscan_70_BL_bins.pl script and the targetscan_70_BL_PCT.pl script to calculate the branch length and the conservation probability of the conserved target sites. Finally, in combination with RNAplfold-2.4.11 59 , the targetscan_70_context_scores.pl script was used to calculate context++ scores for the miRNA target genes based mainly on the 14 attributes mentioned above. Subsequently, potential target genes for the most abundant miRNAs, and DE miRNAs of the 99th context++ score or > 50 context++ score percentile rank were considered for functional enrichment analysis using Ingenuity Pathway Analysis (IPA) web-based software (Redwood City, CA; https ://www.qiage nbioi nform atics .com/produ cts/ingen uity-pathw ayanaly sis/, IPA Spring 2019 release). In addition, Pearson correlation coefficients of most abundant miRNAs or DE miRNAs expression with mRNA expression of their predicted target genes that were expressed in the same animals of our previous study 18 were calculated.

Interaction network analyses between DE miRNA and DE genes.
To explore potential interactions between DE miRNAs identified in this study and the DE genes between high and low-RFI animals that were reported in the same populations reported by Mukiibi et al. 18 , a context++ score greater than the 50 percentile rank for the DE miRNA target genes was used. For each breed, the DE genes identified as targets of the DE miRNAs at this threshold were retained for DE mRNA-DE miRNA (Target-miRNA) interaction network visualization using Cytoscape version 3.7.1 67 , and Pearson correlation coefficients of the DE miRNA expression with mRNA expression of the targeted DE genes were calculated. Furthermore, IPA was used to perform functional enrichment analyses for the targeted DE genes to identify major biological functions that are potentially differentially modulated by the identified DE miRNAs.