Systematic microRNAome profiling reveals the roles of microRNAs in milk protein metabolism and quality: insights on low-quality forage utilization

In this study, we investigated the molecular regulatory mechanisms of milk protein production in dairy cows by studying the miRNAomes of five key metabolic tissues involved in protein synthesis and metabolism from dairy cows fed high- and low-quality diets. In total, 340, 338, 337, 330, and 328 miRNAs were expressed in the rumen, duodenum, jejunum, liver, and mammary gland tissues, respectively. Some miRNAs were highly correlated with feed and nitrogen efficiency, with target genes involved in transportation and phosphorylation of amino acid (AA). Additionally, low-quality forage diets (corn stover and rice straw) influenced the expression of feed and nitrogen efficiency-associated miRNAs such as miR-99b in rumen, miR-2336 in duodenum, miR-652 in jejunum, miR-1 in liver, and miR-181a in mammary gland. Ruminal miR-21-3p and liver miR-2285f were predicted to regulate AA transportation by targeting ATP1A2 and SLC7A8, respectively. Furthermore, bovine-specific miRNAs regulated the proliferation and morphology of rumen epithelium, as well as the metabolism of liver lipids and branched-chain AAs, revealing bovine-specific mechanisms. Our results suggest that miRNAs expressed in these five tissues play roles in regulating transportation of AA for downstream milk production, which is an important mechanism that may be associated with low milk protein under low-quality forage feed.

enrichment in the urine and fecal matter 12 . Therefore, improvement of feed and N efficiency is vital for the dairy industry to meet the challenges in food security and food quality that China is facing.
The dairy industry in China has heavily relied on imported alfalfa hay (AL) to maintain a high production yield in recent years. Alfalfa hay, as high-quality forage, can enhance the feed efficiency and N efficiency in lactating dairy cows 13,14 , however, China does not produce sufficient AL locally. Therefore, the enhancement of the utilization of local forage sources such as corn stover (CS) or rice straw (RS) has been one of the main tasks for China's dairy industry. More than 2.0×10 11 kg of CS are produced annually in China, most of which is wasted, and very little is used as feed sources for livestock 15 . Some have speculated that cows undergo physiological changes when fed low-quality CS, due to the reduced milk productivity and quality, compared to those fed high-quality AL 13,14 . This suggests that a fundamental understanding of the molecular mechanisms of milk production and milk protein synthesis processes under low-quality forages is essential to develop strategies to enhance the potentials of CS as feed for dairy cows.
Therefore, we aimed to investigate the systematic molecular regulatory mechanisms of milk protein production in dairy cows by studying the microRNAomes of five key metabolic tissues, the rumen, duodenum, jejunum, liver and mammary gland, and their roles in dairy feed efficiency and N efficiency because miRNAs are non-coding RNAs that regulate gene expression 16 . We hypothesized that different miRNAs are involved in regulating feed efficiency and N efficiency in these tissues, which can also be influenced by the different forage sources.

Materials and Methods
The experimental protocol was approved by the Animal Care Committee, Zhejiang University, Hangzhou, P. R. China, and all procedures were conducted in accordance with the approved protocol.
Animal study and sample collection. A total of 18 cows were fed three diets with different forage sources: AL (n = 6), RS (n = 6), or CS (n = 6), for 90 days. The details of ingredients and nutritional composition of diets have been described in our previous study 14 . Briefly, diets were isonitrogenous, with a forage-to-concentrate ratio of 45:55 [dry matter (DM) basis] and contained identical concentrate mixtures and 15% corn silage, with different forage sources (on a DM basis): 23% alfalfa hay and 7% Chinese wild rye hay (AL), 30% corn stover (CS), and 30% rice straw (RS). All animals were slaughtered as described in a previous study 17 . Samples of five tissues were collected immediately after slaughter following the methods described by Ash & Baird 18 for rumen epithelium, Burrin et al. 19 for duodenum epithelium and jejunum epithelium, and Laffitte et al. 20 for liver and mammary gland. Tissue samples were immediately frozen in liquid nitrogen and then moved to − 80 °C until they were used for RNA extraction.
Total RNA extraction and miRNA-Seq library construction. Total RNA was isolated from ~100 mg tissue using the mirVana miRNA Isolation Kit (Ambion ® , Life Technologies, Carlsbad, CA, USA) following manufacturer's protocol. The RNA concentration was determined using a Qubit 2.0 Fluorometer (Life Technologies). The quality of RNA was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) with the RNA 6000 NanoLabchip Kit. RNA samples with the RNA integrity number (RIN) more than 7.0 were used for small RNA library construction and quantitative real-time PCR validation.
Equal amounts (1000 ng) of total RNA from each sample were used to construct small RNA libraries using the TruSeq Small RNA Sample Preparation kit (Illumina, San Diego, CA, USA). PCR amplification was performed for 11 cycles. Overall, 90 small RNA libraries were constructed and pooled in equal amounts for gel purification and sequencing, as described in a previous study 21 . All libraries were sequenced at Génome Québec (Montréal, Canada) using the HiSeq 2000 system (Illumina) to generate 50 bp single reads.
Sequencing data analysis. Small RNA sequencing data were analyzed using the methods as outlined by Liang et al. 22 . Briefly, the low-quality reads were discarded from the raw data with CASAVA 1.8 based on chastity, while the reads with acceptable quality were subjected to 3′ adaptor sequence trimming. The sequences with size from 18 to 30 nt were then mapped to the ncRNA sequences (Rfam) to remove non-miRNA sequences (including tRNA, snoRNA, rRNA, and other non-coding RNAs), and the known and novel miRNAs were determined using miRDeep2 23 . Briefly, the filtered sequences were then aligned against the corresponding known miRNA precursor sequences (miRBase release version 21) with the module of quantifier.pl in miRDeep2 (default parameters) to quantify known miRNAs.
The conservation of known miRNAs was described based on the TargetScan definitions for "highly conserved", "conserved", "poorly conserved", and "bovine-specific" 24 . A highly conserved miRNA is conserved across most vertebrates; a conserved miRNA is conserved across most mammals but usually not beyond placental mammals; a poorly conserved miRNA not present in the above three groups; and a bovine-specific miRNA is only identified in bovine species.
Identification of differentially expressed miRNAs. The expression of known miRNAs in each library was normalized to counts per million reads (CPM) as follows: CPM = (miRNA reads number/total mapped reads number per library) × 1,000,000. The miRNAs with CPM > 1 in at least 50% of the samples were defined to be expressed miRNAs in each tissue. Tissue-specific miRNAs were defined as miRNAs uniquely expressed in one tissue. A miRNA expressed in all 5 tissues was defined as a commonly expressed miRNA, and miRNAs that were expressed significantly higher or lower in one tissue than any others (P < 0.05, log 2 (fold change (FC)) > 1 or < − 1, ANOVA) were defined as tissue-dependent differentially expressed (DE) miRNAs. Two-way ANOVA (considering diet and tissue as main effects) and the PROC GLM procedure in SAS (SAS Institute, 2000) were applied to identify tissue-dependent DE miRNAs 25 . The diet effect on miRNAs expression was determined by Scientific RepoRts | 6:21194 | DOI: 10.1038/srep21194 a pair-wise comparison (AL vs CS, AL vs RS, and CS vs RS) using Student's t-test with significance indicated at P < 0.05 and log 2 (FC) > 1 or < − 1.
Correlation analysis. Under each diet, the possible relationships between miRNA expression (CPM) in each tissue and phenotypic traits [feed efficiency (milk yield/DM intake, kg/kg), N efficiency (milk protein yield/crude protein intake (DM basis), kg/kg) and MUN content (mg/dL, an important indicator for N efficiency)] were identified by Pearson's correlation by using R software. The significant correlations were declared at P < 0.05, R > 0.81 (positive correlation) or R < − 0.81 (negative correlation). The miRNAs positively associated with feed and N efficiency (R > 0.81, P < 0.05) and negatively associated with MUN (R < − 0.81, P < 0.05) were considered as positive for dairy performance. These miRNAs and those negatively associated with feed and N efficiency (R < − 0.81, P < 0.05) and positively associated with MUN (R > 0.81, P < 0.05) were subjected to further functional analysis. Correlation between the expression of miRNAs and the expression of their predicted targeted were conducted with Pearson's correlation with R < − 0.81 and P < 0.05 as criteria of significantly targeting. miRNA targets prediction and functional analysis. Target genes for selected miRNAs were predicted by using TargetScan Release 6.0 (http://www.targetscan.org/) and miRanda (http://www.microrna.org/microrna). The target genes that were predicted by both TargetScan (default parameters) and miRanda (Total score ≥ 145, Total energy ≤ − 10) for each miRNA were further analyzed through IPA (Ingenuity Systems, www.ingenuity. com). The multiple testing corrected P value calculated by Benjamini-Hochberg method (FDR) was used to determine the significance of the predicted function in IPA. A threshold of FDR < 0.05 and enriched gene number ≥ 2 were applied to enrich significant biological functions of each miRNA.
Validation of selected miRNA expression results using stem-loop RT-qPCR. The TAQMAN miRNA assay (Applied Biosystems, Carlsbad, CA) was used to validate miRNA expression detected by RNA-Seq following the methods described in a previous study 22 . Briefly, cDNAs were reverse-transcribed from 10 ng total RNA using specific miRNART primers (Applied Biosystems), and DNA was amplified using a 20× TAQMANmiRNA assay (Applied Biosystems). A StepOnePlus ™ Real-Time PCR System (Applied Biosystems) was used to detect the fluorescence signal. The calculation of relative expression of selected miRNAs were conducted as following equation: Δ Ct target miRNA = Ct target miRNA -Ct U6 , with U6 snRNA as internal control 22 . One-way ANOVA was used to compare the differences among groups, and P < 0.05 was considered significant.
Results miRNAome profiling of five bovine metabolic tissues. The miRNA profiles were generated by sequencing 90 small RNA libraries prepared from rumen epithelium, duodenum epithelium, jejunum epithelium, liver, and mammary gland tissues collected from eighteen mid-lactating cows (164 ± 25 days in milk) fed AL (n = 6), RS (n = 6) or CS diet (n = 6) for 90 days. In total, 546 million out of 606 million reads (average 6.73 M per tissue) were mapped to the known miRNAs database (miRBase release version 21), resulting in 340, 338, 337, 330 and 328 known miRNAs and 32, 37, 37, 26, 36 novel miRNAs identified in rumen, duodenum, jejunum, liver, and mammary gland, respectively (Table 1). In total, 388 known miRNAs were expressed in these five tissues. The principal component analysis (PCA) showed clear separations in the miRNA profiles between any two tissues except between duodenum and jejunum. In addition, no clear separations were observed among the three diets (Fig. 1A).
Tissue commonly and specifically expressed miRNAs. In total, 289 miRNAs were commonly expressed in all five tissues (Fig. 1B). Among the top 20 most highly expressed miRNAs in the five tissues, miR-143 was highly expressed (~16%) in all tissues (Table 1), while the expression of some miRNAs varied between tissues, such as miR-122 and miR-205 ( Fig. 2A). The predicted functions of the top 20 commonly expressed miR-NAs included cell death and survival, cell morphology, cellular growth and proliferation, cellular movement and organismal survival ( Fig. 2B and Table S1).
For each tissue, the top 10 most highly expressed bovine-specific miRNAs covered more than 80% of the bovine-specific reads (data not shown). Based on the functional analysis, the three most predominant function categories targeted by the top 10 most highly expressed bovine-specific miRNAs in the rumen, duodenum and jejunum were cellular growth and proliferation, organismal survival, and gene expression ( Fig. 5B and Table S6), with specific predicted functions in regulating abnormal morphology of epithelial tissue (FDR = 1.98E-05) and proliferation of epithelial cells (FDR = 1.46E-04) in rumen (Table S6). For duodenum and jejunum, the main predicted functions were organism survival, cellular growth and proliferation and gene expression (Fig. 5B), with specific functions in regulating epithelial cell death (FDR = 9.74E-04 (duodenum) and 1.05E-03 (jejunum)) and   Table S6). For liver and mammary gland, the main predicted functions were organismal survival, cellular growth and proliferation, and cell death and survival (Fig. 5B), with specific functions in regulating fatty acid metabolism (FDR = 0.01) and accumulation of lipids (FDR = 4.38E-04), and gluconeogenesis (FDR = 1.56E-03) in the liver. The functions targeted by bovine-specific miRNAs included synthesis of lipid, synthesis of carbohydrate, and catabolism of AAs (FDR = 0.02) in mammary gland (Table S6).

Predicted target genes of identified tissue specific and DE miRNAs.
To further predict the function of identified miRNAs that were expressed uniquely among tissues (DE or specific), RNA-Seq data (data not shown) were used to determine the association between the miRNA levels and their predicted targets ( Table 2). Among these genes, the expression of the ATPase Na + /K + transportation alpha 2 polypeptide (ATP1A2) gene was predicted to be negatively correlated with miR-21-3p (R = − 0.54, P = 0.02). Additionally, the solute carrier family 7, member 8 (SLC7A8) transcript was predicted to be targeted by miR-2285f (R = − 0.53, P = 0.02). These two genes were predicted to be involved in AA transportation by Panther.
Quantitative real-time PCR validation of expression of dairy efficiency-associated miRNAs. MiR-21-3p (rumen up-regulated) and miR-2285f (bovine-specific and liver up-regulated) were selected for qPCR validation due to their targeting of genes involved in amino acids transport. The qPCR analysis confirmed that the expression of miR-21-3p was higher in the rumen and miR-2285f was up-regulated in the liver, compared with the other four tissues (P < 0.05, Fig. 6). Additionally, miR-21-3p was down-regulated by the AL diet in the rumen and mammary gland compared with the RS and CS diets (P < 0.05, Fig. 6).
Effects of diet on dairy efficiency and expression of dairy efficiency-associated miRNAs. The expression of some dairy efficiency-associated miRNAs was also affected by the diet. When cows were fed AL diet, the expression of some of the dairy efficiency positively associated miRNAs was different comparing to those under RS and/or CS diets. For example, the expression of rumen specific miRNA miR-2299-3p (one of the N efficiency positively associated miRNAs) tended to be lower than that that under RS (P = 0.094) or CS (P = 0.056) diet (Table S7), while the expression of miR-103 in the rumen under RS (P = 0.004) or CS (P = 0.034) diet was higher than that under AL diet (Table S7).

Discussion
The functions of miRNA in mammalian animals have been widely studied 26 , including their potential roles in regulation of nutrient metabolism in different tissues in humans 27 , pigs 28 and cows 29,30 . MiRNAs may play important roles in lactation initiation and milk component (protein and fat) synthesis because certain miRNAs are highly expressed in the mammary gland of dairy cows 31,32 and expression of mammary gland miRNAs can be affected by the diet 33 . Thus, miRNA-mediated regulation in the mammary gland is considered important for the lactation process 31,32 . However, most of these studies only focused on mammary gland tissue and lactation, which precludes a comprehensive understanding of the molecular mechanisms of milk production, especially milk protein synthesis.
In the current study, we focus on studying the miRNAs and their roles in protein metabolism from various metabolic tissues in dairy cows. To our knowledge, this is the first study to compare miRNA expression profiles and their functions in different tissues of the same animal and their changes under different diets using RNA-Seq. Among all the tissues of the cows, rumen, duodenum, jejunum, liver, and mammary gland are considered to be the most important organs for protein metabolism [33][34][35][36] . Rumen is the specific organ in ruminant animals where microbial crude protein synthesis is performed and N sources (AAs, peptides and ammonia) are absorbed 36 . Duodenum and jejunum are sites for protein digestion and absorption 35 , while liver is the key site for AA re-synthesis and urea cycle 34 . Finally, milk protein is synthesized in mammary gland 33 .
In this study, the 289 commonly expressed miRNAs identified across these five tissues suggest that the functions involved in the same biological or physiological processes in these tissues may be partially regulated by these miRNAs. For both tissue-DE and tissue-specific miRNAs, the most commonly affected enriched functions were basic cellular functions, such as cellular morphology and cellular movement, which were similar to each other, suggesting basic cellular functions could be mediated by different miRNAs. At the same time, the identification of  tissue-DE and tissue-specific miRNAs suggests that these miRNAs may play critical roles in the specific functions of these cell types 37 , contributing to the different physiological functions of the five tissues [33][34][35][36] . More fundamental studies should be performed to investigate how tissue-DE and tissue-specific miRNAs play their roles in each tissue. Regardless, abundant miRNAs from these five tissues suggest their unique roles in the digestive tract. From detected miRNA in mammary gland, only a few of them overlapped between our current study and previous reports 31,32 . This could possibly due to different physiological stages (dry and early lactating 31,32 vs. mid-lactating in this study) and different methods for sampling tissues (biopsy 31,32 vs. slaughter in this study). However, some of the previously reported miRNAs that regulate lactation initiation in mammary gland 31,32 were detected in all five metabolic tissues and were likely to be involved in feed and N efficiency regulation of the dairy cows. For example, miR-125b, miR-141 and miR-181a have been associated with milk lipid synthesis in mammary gland of lactating cows 31 , and the expression of miR-221 and miR-205 is correlated with lactating initiation 32 . In our study, the expression of miR-125b, miR-141, miR-181a, miR-221 and miR-15b were detected in all five tissues (Table S2) and they were associated with feed or N efficiency under different diets ( Table 2 and Table S7). Cell proliferation is an important function in lactation initiation, milk precursor absorption and milk component synthesis in mammary gland as well as in rumen and adipose tissues 38,39 . There is evidence that miR-125b, miR-141, miR-181a, miR-221 and miR-205 are involved in cellular proliferation [40][41][42][43][44] , suggesting that these miRNAs play roles in lactation initiation, feed efficiency and N efficiency in dairy cows at different physiological stages by regulating cell proliferation. Moreover, the expression of miR-125b in jejunum and miR-141 in mammary gland were impacted by the diet (Table S7), suggesting that nutritional changes could influence the expression of those miRNAs which could regulate milk initiation as reported as well as feed efficiency, and N efficiency as observed in this study.
Our results on identification of miRNAs associated with feed efficiency, N efficiency or MUN content in all five tissues suggested that miRNAs may be involved in nutrient utilization (indicated by feed efficiency) and N source utilization (indicated by N efficiency or MUN content) through the gut, in addition to mammary gland. For example, miR-130b, a dairy efficiency negatively associated miRNA with lower expression under AL comparing to both RS and CS in the rumen (Table S7), has been reported to regulate cellular proliferation in muscle. Another dairy efficiency negatively associated miRNA, miR-378, which had highest expression in the rumen under CS (Table S7), is associated with cellular proliferation and differentiation 46 . In addition, hepatic miRNAs such as miR-493 (under RS diet) and mammary let-7b (higher in AL diet, compared with CS diet) have been reported to be involved in cellular proliferation 47,48 . Thus, let-7b may regulate AAs metabolism by impacting cell proliferation in mammary gland. These further suggest that enhanced cellular proliferation and differentiation in various tissues by AL could be an import mechanism that could contribute to the variation in diary efficiency. Some miRNAs were associated with more than one of the dairy efficiency traits in different tissues. For example, under RS diet, the expression of miR-125b was positively associated with feed efficiency in duodenum and jejunum (Table S7). MiR-125b were previously reported to regulate cellular proliferation 49 , so the expression change of this miRNA may impact the proliferation rates of the duodenum and jejunum epithelium and further regulate nutrients utilization. Similarly, miR-15b may impact cellular proliferation 50 in mammary gland by regulating feed efficiency under CS diet. Moreover, AA metabolism in dairy cow impacted N efficiency 51 . When cows were fed AL diet, lower expression of ruminal miR-16b, -181d, indicating their less inhibitory effect on the genes that may regulate the AA metabolism. Overall, the data show that miRNAs may regulate protein metabolism in all tissues, suggesting the miRNAs have systematic regulatory functions in feed efficiency and N efficiency ( Table 2). Our results suggest that milk production is a process that could be regulated by miRNAs in all five of these tissues. Moreover, regulating AAs transport function was a method to regulate N efficiency 52 , and phosphorylation of AAs was one mechanism for regulation of AAs transporter 53 . Our study revealed that systematic phosphorylation (duodenal miR-1296, jejunal miR-30b-3p, hepatic miR-199a-3p and mammary miR-196 (data not shown)) and transport of AAs (ruminal miR-130a, duodenal miR-1296, jejunal miR-30b-3p and mammary miR-196 (data not shown)) may be conducted by mediated-miRNAs, suggesting that miRNA-mediated mechanisms regulated N efficiency by regulating AAs transport and phosphorylation throughout all these five tissues. More studies should be conducted to reveal the mechanism of the systematic association between miRNAs and efficiency (feed and N) in different tissues.
It is known that transport of AAs plays an important role in dairy feed and N efficiency 55 . MiRNAs associated with N or feed efficiency in the five tissues targeted genes involved in AA transportation and metabolism, some of which are tissue-specific or tissue-DE. Different tissues may have unique miRNA-mediated regulation mechanisms for the same function, possibly induced by tissue-DE or tissue-specific miRNAs 21 . For example, ATP1A2 may regulate Na + /K + pump and hydrolysis of ATP 56 , which is important in AA absorption of gut epithelial cells 36 . The negative association (R = −0.87, P < 0.05) between miR-21-3p and ATP1A2 identified in the current study suggests the importance of miR-21-3p in AA transportation through the rumen epithelium tissue, which could impact milk production. In addition, the relatively higher expression of miR-21-3p in rumen under RS diet (Fig. 6A) indicate the diet may affect this function, which may lead to the altered AA metabolism contributing lower milk protein quality. Furthermore, functions including transport, metabolism and oxidation of AAs regulated by dietary DE (between high-and low-quality forage diets) miRNAs were associated with N efficiency or MUN content. For example, miR-21-3p was down-regulated in both rumen (Fig. 6C) and mammary gland (Fig. 6E) in the AL compared with the RS and CS groups, suggesting that miRNAs may respond systematically to dietary changes and target similar functions. In summary, miRNA-mediated regulation may play a role in AA metabolism in the five tissues, but the mechanisms could be varied among tissues.
Due to anatomical differences, many physiological metabolic processes are different between ruminants and monogastric animals, and these processes can be regulated by host-specific miRNAs. The current study indicates that bovine-specific miRNAs regulate abnormal morphology and proliferation of epithelial tissue, which may contribute to the nutrient absorption by the rumen 57 . Therefore, although bovine-specific miRNAs had low expression (accounting for 0.3% of the total reads), we speculate that they are involved in ruminant-specific functions. Branched chain amino acid (BCAA) (such as leucine) catabolism is lower in ruminants than monogastric animals due to a lack of BCAA aminotransferase 58 . Our results showed that a liver-up-regulated bovine-specific miRNA, miR-2285f (Fig. 6B), targets SLC7A8, a gene involved in high-affinity transport of small and large neutral AAs, including leucine. This suggests that miR-2285f may inhibit leucine membrane transportation by down-regulating SLC7A8 to reduce the potential accumulation of leucine in the liver 58 . In addition to AA metabolism, bovine-specific miRNAs may regulate lipid and fatty acids metabolism in liver, which is also different from monogastric animals, especially in terms of the relatively lower-density lipoprotein in dairy cows 59 . Some of the bovine-specific miRNAs, such as miR-2285t (Higher expressed under RS, compared with CS) and miR-2419-5p (lower expressed under AL diet, compared with CS ) in jejunum, were also involved in AA transport. These results suggest the potential functions of bovine-specific miRNAs in epithelium proliferation/morphology, AA metabolism, and other specific functions.

Conclusions
By using the high throughput RNA-Seq, our results showed how miRNAs are involved in biological process or physiological functions in five key metabolic tissues in the dairy cow, considering the animal as a whole biological system. By profiling miRNAs in key tissues of dairy protein metabolism (rumen, duodenum, jejunum, liver and mammary gland), our study reveals that 1) in addition to mammary gland, miRNAs from rumen, duodenum, jejunum and liver could contribute to feed and N efficiency systematically in mid-lactating dairy cows; 2) each tissue has its functional uniqueness due to their physiological and functional differences, which could be regulated by miRNAs, and miRNA-mediated mechanisms could be involved in AA metabolism in the five tissues; 3) bovine-specific miRNAs may be involved in unique functions, such as BCAA metabolism; 4) diet can significantly change the expression of dairy efficiency associated miRNAs which suggests that nutritional alterations could be a method to regulate the miRNA profile in these protein metabolism-related tissues, and may help with the feeding management of the dairy cows; 5) Low-quality forage diets can impact AA metabolism in the gut by altering miRNA expression, and post-transcriptional regulation patterns in regulation AAs metabolism were differed among different diets, suggesting potential ways to improve milk protein quality by manipulation of miRNAs under low-quality forage. By identifying miRNAs that regulate feed and N efficiency and measuring their changes under diets with different forages, our results provide fundamental understanding on the molecular regulatory mechanism through miRNA regulation which could explain the reduced feed and N efficiency under low-quality forage diets. Whether these miRNAs regulate these traits are undefined and future studies using in vitro system such as cell-lines (mammary gland and/or epithelial cells) or in vivo feeding trials (feeding with nucleotides based on identified miRNAs) are needed to verify the findings from current study.