Skip to main content

Deep RNA-Seq reveals miRNome differences in mammary tissue of lactating Holstein and Montbéliarde cows



Genetic polymorphisms are known to influence milk production and composition. However, the genomic mechanisms involved in the genetic regulation of milk component synthesis are not completely understood. MicroRNAs (miRNAs) regulate gene expression. Previous research suggests that the high developmental potential of the mammary gland may depend in part on a specific miRNA expression pattern. The objective of the present study was to compare the mammary gland miRNomes of two dairy cow breeds, Holstein and Montbéliarde, which have different mammogenic potentials that are related to differences in dairy performance.


Milk, fat, protein, and lactose yields were lower in Montbéliarde cows than in Holstein cows. We detected 754 distinct miRNAs in the mammary glands of Holstein (n = 5) and Montbéliarde (n = 6) midlactating cows using RNA-Seq technology, among which 738 were known and 16 were predicted miRNAs. The 25 most abundant miRNAs accounted for 90.6% of the total reads. The comparison of their abundances in the mammary glands of Holstein versus Montbéliarde cows identified 22 differentially expressed miRNAs (Padj ≤ 0.05). Among them, 11 presented a fold change ≥2, and 2 (miR-100 and miR-146b) were highly expressed. Among the most abundant miRNAs, miR-186 is known to inhibit cell proliferation and epithelial-to-mesenchymal transition.

Data mining showed that 17 differentially expressed miRNAs with more than 20 reads were involved in the regulation of mammary gland plasticity. Several of them may potentially target mRNAs involved in signaling pathways (such as mTOR) and lipid metabolism, thereby indicating that they could influence milk composition.


We found differences in the mammary gland miRNomes of two dairy cattle breeds. These differences suggest a potential role for miRNAs in mammary gland plasticity and milk component synthesis, both of which are related to milk production and composition. Further research is warranted on the genetic regulation of miRNAs and their role in milk synthesis.


The mammary gland (MG) is a complex secretory organ composed of different cell types that interact to ensure proper mammary function and milk synthesis. The rate of MG development and dairy potential differ among breeds. For instance, milk, fat, and protein yield were greater for Holstein (22.7, 0.8, and 0.7 kg/day, respectively) than for Montbéliarde cows (17.6, 0.6, and 0.5 kg/day, respectively) in mountain grazing conditions [1]. Similarly, a comparison between Holstein and Montbéliarde cows fed a maize silage-based diet showed higher milk, fat, and protein yield in Holstein (35.7, 1.4 and 1.1 kg/day, respectively) than in Montbéliarde cows (28.1, 1.1, and 0.9 kg/day, respectively) [2]. The fat yield differences were accompanied by differences in milk fatty acid composition [2]. Genetic and, more recently, genomic selection of dairy cows has led to increased milk production resulting from the increased secretory capacity of the MG [3]. The mammogenic potential in dairy heifers in comparison with beef heifers is facilitated by an increased number of mammary stem cells and the differential expression of genes involved in the development of the mammary stem cell niche. These genes influence the proliferation, migration, differentiation, and remodeling of mammary tissue and the regulation of adipocyte transdifferentiation [4]. As a result, genetic polymorphisms that differentiate species and breeds also influence milk production and composition. Despite the increase in knowledge about the genetics of dairy cows in recent decades, the genomic mechanisms influencing milk secretion and composition are not fully understood.

MicroRNAs (miRNAs) are small noncoding RNAs with 18–25 nucleotides. They regulate gene expression by base-pairing with mRNA to induce their degradation or to inhibit their translation [5]. MiRNAs are thought to regulate at least 60% of genes. Therefore, they are involved in many different cellular processes [6,7,8,9]. Nevertheless, the genetic regulation of miRNAs is poorly understood. A recent study identified 125 differentially expressed miRNAs in the kidney among 3 distinct porcine breeds. This suggests that miRNAs could be used for the study of the genetic variability underlying complex traits [10]. Similarly, the identification of 50 miRNAs that were differentially abundant in serum from Warmblood horses (Arabian, Anglo-Arabian, Selle Français, Cob Normand, French trotter and Trait du Nord breeds) and ponies (Shetland and Welsh breeds) suggests the potential of miRNA to serve as biomarkers of different equine breeds [11]. In cows, the miRNomes of ruminant MGs were investigated using RNA-Seq, which is a powerful tool to characterize a large miRNA repertoire [12,13,14,15,16]. The roles of miRNAs during mammary development have also been reported [17,18,19]. While the physiological [15] and nutritional [20,21,22] regulation of MG miRNomes in bovine and caprine species have been reported, the genetic influence on miRNA gene expression is still limited. A comparison of MGs from different bovine breeds with contrasted mammogenic potential found 54 differentially expressed miRNAs (DEMs) in mammary tissues in dairy (Holstein-Friesian) and beef (Limousin) postpubertal heifers. These results suggest that the high developmental potential of the MG in dairy cattle, leading to high milk production, may depend on the specific miRNA expression pattern [23]. Thus, the objective of the present study was to compare MG miRNomes in Holstein and Montbéliarde cows, which are two dairy breeds with different dairy performances.

Results & discussion

Milk production and composition

Milk production and composition were measured before mammary biopsies to characterize the performance of cows in the present study. The milk yield was higher in Holstein than in Montbéliarde cows (29.1 vs 23.9 kg/day; Fig. 1a). Fat, protein and lactose yields were higher in Holstein (972, 869, and 1460 g/day, respectively) than in Montbéliarde cows (846, 756, and 1179 g/day, respectively; Fig. 1b). These differences are in agreement with those found in previous studies comparing lactating Holstein and Montbéliarde cows fed diets based on maize silage [2] and grazed grass [1].

Fig. 1

Milk production (a) and composition (b). A comparison between Holstein (n = 5) lactating cows and Montbéliarde (n = 6) lactating cows. *** P < 0.001. ** P < 0.01. * P < 0.05. Bars are SEM

Global description of mammary miRNomes

Five and six libraries were constructed using RNA extracted from the MG of lactating Holstein and Montbéliarde cows, respectively. The high-throughput sequencing performed allowed us to obtain more than 12 million raw reads on average (Table 1). After cleaning (poly-A stretches and adaptors were removed), 11.0 and 8.5 million cleaned reads were mapped onto the Bos taurus genome from the Holstein and Montbéliarde libraries, respectively. The percentage of total mapped reads ranged from 96.1 to 98.3% and was comparable between the libraries (Additional file 1: Figure S1). The correlation between the libraries that was calculated using the log2 of the normalized counts of expressed miRNAs was R = 0.96, indicating a strong correlation between the Holstein and Montbéliarde libraries (Additional file 1: Figure S2). A total of 754 miRNAs were identified using miRDeep2 software, 738 of which were known and 16 were predicted miRNAs. The latter may be considered to be potentially novel miRNAs. These results are in accordance with data previously obtained using NGS technology in bovine MG that reported 487 known bovine miRNAs, 167 of which were miRNAs that were already known in other species [14]. However, the number of predicted novel miRNAs was higher in the study by Le Guillou et al. [14] than in our study (679 vs 16, respectively). In the present study, we used the latest miRBase version, which includes a large number of known miRNAs, and this could explain the lower number of predicted miRNAs. The 25 most abundant miRNAs in the MG represented 90.6% of the total reads. Among them, four (miR-143, miR-30a-5p, miR-148a, and miR-10b) represented more than 50% of the total reads (Fig. 2).

Table 1 Sequencing data from Holstein (n = 5) and Montbéliarde (n = 6) lactating cows. Clean reads were after size (15 to 40 nt), adapter and soft-clipping cleaning
Fig. 2

Top 25 most abundant miRNAs in the mammary gland of lactating cows. Representation of the percentage of reads of each miRNA relative to the total reads of the Holstein (n = 5) libraries and Montbéliarde (n = 6) libraries

Differentially expressed miRNAs in the mammary gland of Holstein and Montbéliarde cows

A comparison of the MG miRNomes of Holstein and Montbéliarde cows using the DESeq2 package allowed the identification of 22 miRNAs that were differentially expressed (Padj ≤ 0.05; Table 2) between the breeds. Among them, 11 were up- and 11 were down-regulated in Holstein cows compared to Montbéliarde cows, and 17 were found to correspond to more than 20 reads on average. Six (miR-16a, miR-186, miR-25, miR-100, mir-30e-5p and miR-146b) corresponded to more than 5,000 reads on average. Among these 22 DEMs, 11 presented a fold change ≥2, and miR-100 and miR-146b were highly expressed miRNAs. Breed-specific patterns in miRNomes have also been observed in porcine kidney [10], equine serum [11], and porcine and bovine MG [23, 24]. In porcine MG, a miRNome comparison between the Jinhua and Yorkshire breeds identified 391 DEMs [24]. This higher number of DEMs might be due to the large differences in the numbers of predicted miRNAs. Indeed, the stringency of our analysis allowed us to detect only 16 predicted miRNAs, whereas the use of a BLAST strategy identified 2823 and 2286 potential miRNAs in Jinhua and Yorkshire breeds, respectively. In addition, the breed comparison in the porcine study showed large differences in terms of genetic selection and therefore in lactation performance. The Jinhua is a traditional breed that produces less milk than the selected Yorkshire breed. The large divergence in phenotype and the difference in the predicted number of miRNAs could in part explain the detection of a higher number of DEMs in the porcine study. Nevertheless, five DEMs (miR-7, miR-19a, miR-19b, miR-30e-5p, and miR-1271) were identified in both the present and the porcine study, including miR-30e-5p, which was highly expressed in MG in cows (Fig. 3). In bovine MG, Wicik et al. [23] identified 54 DEMs in beef and dairy heifer miRNome comparisons. In the present study, Holstein and Montbéliarde cows have closer milk performance than Holstein and beef Limousin cattle compared by Wicik et al. [23]. Therefore, it is not surprising that the MG miRNomes were more similar in Holstein and Montbéliarde cows than when comparing beef and dairy breeds. Among the 22 DEMs identified in the present study, six DEMs (miR-16a, miR-186, miR-25, miR-409a, mir-199c and miR-146b) were among the 54 DEMs detected by Wicik et al. [23] (Fig. 3), and four (miR-16a, miR-186, miR-25, and miR-146b) were highly expressed in our study. The comparison between the DEMs found in the studies on the effects of breed on the MG miRNomes in husbandry animals (from present study, Peng et al. [24] and Wicik et al. [23]) revealed two miRNAs, miR-25 and miR-186, that are shared. Therefore, these two miRNAs are of interest as their expression could be related to the selection of genes involved in milk production.

Table 2 Differentially expressed miRNAs. MiRNAs whose expression was different in lactating mammary glands of Holstein (H) and Montbéliarde (M) cows. afor n = 5 and bfor n = 6. MiRNAs with higher abundance in Holstein cows are in italic and miRNAs with more than 20 reads on average are in bold. Means correspond to the mean of normalized read counts
Fig. 3

Venn diagram of DEMs in different mammary gland studies. The list of the comparison between Jinhua and Yorshine swine and between Holstein and Limousin heifers were from Peng et al. [23] and Wicik et al. [24], respectively. The name of the common miRNAs with the present study was indicated. In bold there were those highly expressed

Putative functions of the 17 differentially expressed miRNAs

Bioinformatic analyses were performed to identify the biological processes potentially regulated by the 17 DEMs with more than 20 reads (on average) and one predicted miRNA (chr21_22372_star), which was not used for the functional analysis. Among the first 20 identified process networks, 12 were linked to development, the regulation of the cytoskeleton, the establishment and regulation of tissues, the cell cycle, and apoptosis (Table 3). These results suggest a role played by the 17 DEMs in the morphology and functioning of mammary tissue. In particular, one process network was linked to epithelial-to-mesenchymal transition (EMT), which has already been reported to be important for MG development and remodeling [25]. Indeed, the polarization of mammary tissue architecture is crucial for the maintenance of transcription factor activation, chromatin organization, and tissue-specific gene expression, which influence milk synthesis [26].

Table 3 Biological processes potentially regulated by the 17 DEMs. The 17 miRNAs were selected as having more than a mean of 20 reads. Process Networks analysis using Metacore™ software. * indicate processes linked with cytoskeleton and cell life and § indicate signal transduction processes. FDR: False Discovery Rate

The first modified process was NOTCH signaling transduction. This signaling pathway is involved in cell differentiation and influences mammary development by promoting mammary epithelial stem cell activity [27]. The second identified signal transduction process was Wnt signaling. Wnt proteins influence the self-renewal of stem cells in the MG [28], and recent data indicate that Wnt signaling is also involved in mammary stem cell maintenance by generating MG plasticity in mammary lineage cells [29]. The interactions between epithelial cells and adipocytes coordinate MG development and influence milk production [30]. During lactation, cell turnover was estimated to be ~ 50% due to cell proliferation and apoptosis [3], thereby showing the importance of cell life during lactation. More recently, two mammary epithelial cell populations were linked to milk production, showing the importance of these types of cells [31]. These previously reported data reinforce the importance of MG plasticity for efficient MG development and lactation. Our results based on algorithmic prediction suggest that miRNAs may influence the cell life cycle, cell differentiation and transdifferentiation (e.g., EMT). The occurrence of these mechanisms in MG are still poorly documented and need further research.

Role of the 6 most differentially expressed miRNAs

The hierarchical classification of DEMs using PermutMatrix software allowed the classification of cows according to their breed and the clustering of miRNAs according to their level of expression (Fig. 4). Three groups were identified. The first group corresponded to miR-186 and miR-30e-5p, both of which showed a high expression level. The second group comprised miR-25, miR-16a, miR-100, and miR-146b, which showed an intermediate expression level. The third group comprised DEMs with low expression. Only the first and second groups (containing 6 miRNAs: miR-100, miR-146b, miR-186, miR-30e-5p, miR-25, and miR-16a) were considered for further investigation and discussion. RT-qPCR analyses were performed for miR-25, miR-146b, miR-100, and miR-186, and miR-25, miR-146b, and miR-186 showed the same expression pattern as that obtained using RNA-Seq (Additional file 1: Table S1).

Fig. 4

Two-way hierarchical clustering dendrogram of differentially expressed miRNAs. Clustering performed with data of the 11 mammary samples of the 19 DEMS (Padj < 0.05), and using PermutMatrix software (Caraux and Pinloche, 2005) using Euclidean distance calculation and Ward’s minimum variance method. Columns represent each mammary sample for Holstein (H, samples 1 to 5) and Montbéliarde (M, samples 1 to 6) cows. Rows represent each miRNA grouped according to abundance (from Group 1 with higher to Group 3 with lesser abundance)

Putative roles of miR-100 and miR-146b, the two most highly expressed miRNAs with fold-changes > 2 in the mammary gland

MiR-100 and miR-146b exhibited fold-changes (FC) > 2 and mean counts > 8,000. MiR-146b is involved in the control of the immune system [32] by suppressing the expression of inflammatory cytokines and inducing the inhibition of autophagy by targeting the PTEN/Akt/mTOR signaling pathway [33]. MiR-100 is identified in numerous cancers [34,35,36] and acts via mTOR signaling in human mesenchymal stem cells [37]. Elsewhere, mTOR signaling regulates cell growth, proliferation, and cell life, as well as gene transcription and protein synthesis in the MG and could be related to lactation [38, 39]. In addition, mTOR was reported to affect lipogenic gene networks in bovine mammary epithelial cells, indicating a potentially important role in the regulation of milk fat synthesis [40]. The potential regulation of mTOR signaling by miR-100 and miR-146b and their differential expression in the present study could be related to the differences in milk fat and protein yield between Holstein and Montbéliarde cows. However, this indirect relationship has to be studied further. In addition, miR-100 is reported to induce EMT in a model of mammary epithelial cells, which is in line with the identified biological processes (Table 3) and mammosphere formation [41].

Putative roles of the four DEMs with high expression: miR-186, miR-30e-5p, miR-25 and miR-16a

MiR-186, miR-30e-5p, and miR-16a were more abundant and miR-25 was less abundant in Holstein cows than in Montbéliarde cows. In particular, miR-186 was one of the 25 most abundant miRNAs in MG tissue (Table 2). MiR-186 inhibits cell proliferation in prostate [42], gastric [43], melanoma [44], and colorectal [45] cancers. In contrast, it has been reported that miR-186 promotes cell proliferation in human melanoma [46]. The mechanism of its action is still unclear, but miR-186 seems to influence cell proliferation. In addition, the downregulation of miR-186 expression was associated with the EMT phenotype in cisplatin-resistant ovarian cancer [47] and colorectal cancer [45]. MiR-186 was also reported to regulate glucose uptake and lactate production by targeting the 3′-UTR of glucose transporter 1 (Glut1) mRNA in cancer-associated fibroblasts [48]. Glucose plays a key role in milk synthesis as a precursor for the synthesis of lactose, which is a major driver of milk volume [49]. However, the upregulation of miR-186 in the MG in Holstein cows is not consistent with the higher lactose yield observed in Holstein cows. Thus, the potential role of miR-186 in the regulation of lactose synthesis and milk yield needs further investigation.

MiR-30e-5p was the second most abundant miRNAs among the DEMs and showed a higher abundance in the MG in Holstein cows. The expression of miR-30e-5p in goat mammary epithelial cells was associated with the upregulation of the expression of genes involved in lipid metabolism, such as PPARγ, LPL, DGAT1, and CD36, to promote milk fat synthesis [50].

MiR-25 was the third most highly expressed miRNAs among the DEMs and had greater expression in Montbéliarde cows than in Holstein cows. This miRNA represses triacylglycerol synthesis and lipid accumulation in goat mammary epithelial cells [51]. Indeed, miR-25 directly targets peroxisome proliferative activated receptor gamma coactivator 1 beta (PGC-1beta) and reduces the mRNA levels of SREBP1, FASN, GPAM and PPARG [51]. The expression of miR-16a was lower in Montbéliarde cows than in Holstein cows. MiR-16a was reported to target genes involved in lipid metabolism [52].

Links between milk production and composition and DEMs

Bibliographic analyses indicated that the EMT process, which could be related to the structure of mammary epithelium-synthetizing milk components, could be modulated by highly expressed DEMs. Elsewhere, genes involved in the mTOR signaling pathway, reported to be related to protein and lipid synthesis [38,39,40], were predicted to be targeted by highly expressed DEMs. Lipid metabolism was previously shown to be regulated by miR-30e-5p, miR-25 and miR-16a. Therefore, we can suggest a potential role for the highly expressed DEMs in milk component synthesis.

Furthermore, Spearman correlation analysis was performed showing significant correlations between miR-186, miR-30e-5p, and miR-16a expression and milk, protein, fat and lactose secretion yields (Additional file 1: Table S2; P ≤ 0.05). These results supported the hypothesis regarding the potential role of several DEMs in milk production and composition. However, the direct role of these miRNAs in milk synthesis needs to be confirmed by functional analyses, which would provide compelling evidence of their role.


We found 22 DEMs by comparing the MG miRNomes between two dairy cattle breeds with different dairy performance (Holstein cows and Montbéliarde cows). The bioinformatic analyses of their predicted target genes identified genes involved in signaling networks and tissue structure. We identified 6 miRNAs (miR-100, miR-146b, miR-186, miR-30e-5p, miR-25, and miR-16a) which were highly expressed in lactating MG in both breeds. These microRNAs are predicted to target genes involved in the synthesis of milk constituents. In addition, several DEMs are predicted to target the mTOR signaling pathway and genes involved in lipid metabolism, which could be related to the lower milk fat yield in Montbéliarde cows compared to that in Holstein cows (Fig. 5). Among the 6 highly expressed DEMs, we found a correlation between 3 DEMs and milk production and composition supports this hypothesis. Nonetheless, further research is warranted on the genetic regulation of miRNAs and their precise functional roles.

Fig. 5

Biological processes and genes potentially modulated. Biological processes potentially regulated (Padj ≤ 0.05) by the 19 DEMS and genes potentially targeted by the 6 most expressed and/or differentially expressed miRNAs. The modified process networks are represented in green and the selected miRNA targeting genes are represented in blue


Animals and sampling

Animal procedures were performed in compliance with Regional Animal Care Committee guidelines CEMEAA: Auvergne, French Ministry of Agriculture and European Union guidelines for animal research C2EA-02. All procedures were approved by the ethics committee on animal experimentation (3737-2015043014541577 V2).

Twelve cows (6 Holstein and 6 Montbéliarde), with an average of 3.3 ± 1.5 as lactation number were studied at mid-lactation (165 ± 21 days in milk (DIM)). However, one Holstein cow was excluded from the experiment due to mastitis. The diet was the same for all animals and was composed of forage and concentrate (74:26 of dry matter). Milk yield and composition were recorded on 9 consecutive days before biopsies. The milk was analyzed via mid-infrared spectroscopy to determine the fat and protein content (Galilait, Theix, France).

Mammary biopsies were performed as previously described [53]. The biopsy site was selected at a midpoint on a rear quarter. The collected tissue, approximately 600–650 mg, was rinsed in sterile 0.9% saline solution, inspected to verify tissue homogeneity, and snap-frozen in liquid nitrogen and stored at − 80 °C until RNA extraction.

RNA preparation and RNA-Seq analysis

Total RNA was extracted from on average 50 mg MG biopsies (n = 5 Holstein and n = 6 Montbéliarde) using miRVana kit (Thermo Fisher Sciences, USA). The concentration and purity of the RNA was estimated by spectrophotometry NanodropTH (ND-1000, NanoDrop Technologies LLC, Wilmington, DE, USA) and using the Bioanalyzer 2100 (Agilent Technologies Inc., Santa Clara, CA, USA), respectively. Means of RIN were 8.0 and 7.8 in Holstein and Montbéliarde, respectively. The preparation of the libraries were performed by the IGBMC Microarray and Sequencing Platform (Strasbourg, France). Small RNA-Seq librairies were generated from 2000 ng of total RNA using TruSeq Small RNA Library Prep Kit (Illumina, San Diego, CA), according to manufacturer’s instructions. The protocol takes advantage of the natural structure common to most microRNA molecules that have a 3′ hydroxyl group resulting from enzymatic cleavage by Dicer or other RNA processing enzymes. Briefly, during the first step, RNA adapters were sequentially ligated to each end of the RNA; firstly the 3′ RNA adapter (5′ TGGAATTCTCGGGTGCCAAGG 3′) which is specifically designed to target microRNAs and other small RNAs, then the 5′ RNA adapter (5′ GTTCAGAGTTCTACAGTCCGACGATC 3′). Small RNA ligated with 3′ and 5′ RNA adapters were reverse transcribed and PCR amplified (30 s at 98 °C; [10 s at 98 °C, 30 s at 60 °C, 15 s at 72 °C] × 13 cycles; 10 min at 72 °C) to obtain cDNA. Acrylamide gel purifications of 140–160 nt amplified cDNA (corresponding to cDNA obtained from small RNA + 120 nt from the adapters) were performed. The final cDNA libraries were checked for quality and quantified using capillary electrophoresis. Libraries were loaded in the flowcell at 2.8 nM and clusters were generated using Cbot and sequenced on HiSeq 4000 (Illumina) as single-end 50 base reads, according to the manufacturer’s instructions. The quantity and quality of reads for each library are shown in Additional file 1: Table S1. After trimming of adaptor sequences and removal of reads containing ambiguous base calls (FASTX-Toolkit,, reads were filtered according to their size (15–40 nt). Reads quantification and annotation were performed using the ncPRO-seq pipeline [54]. The sequence reads were aligned against the Bos taurus btau5.0.1 genome as miRBase_v22.1 using the miRDeep2 package [55]. Precursors and mature miRNAs were identified using the miRDeep2 core module, Potential miRNAs were annotated accordingly against ortholog miRNAs in goat, sheep and humans (miRBase release 22.1). We used a miRDeep2 score ≥ 0 as a cut-off threshold. The accession number of the RNA-Seq data is GSE131057. The X- (Holstein) and Y- (Montbéliarde) axes show the Log2 of the mean of the normalized counts of expressed miRNAs (Additional file 1: Figure S2).

RT-qPCR analyses

RT-qPCR were performed to confirm the expression pattern, using TaqManTM Advanced miRNA cDNA synthesis and TaqManTM Advanced Master Mix (Life Technologies, Foster City, CA, USA), according to the manufacturer’s instructions. Reverse-transcriptions were achieved on 10 ng of total RNA. PCR were performed using StepOnePlus (Applied Biosystems, Foster City, CA, USA) at 95 °C for 20 s, pursued by 40 cycles of 95 °C for 1 s and 60 °C for 20 s. Used amplification systems (miR-100, -25, -146b, -186, -26a, -92a, and -191) are presented in Additional file 1: Table S1. All miRNA levels were normalized to the mean of the values of 3 internal control miRNAs (miR-26a [56], -92a [57], and -191 [58]). These miRNAs were identified as suitable internal controls using geNorm (Additional file 1: Figure S3). The results were expressed as log2 of fold changes of threshold cycle (Ct) values relative to the control [59].

Statistical and bioinformatics analyses

Statistical analyses to compare breed effects on the milk production means per cow were performed using mixed models of SAS (version 9.4; SAS Institute INC, Cary, NC). The normalization and differential expression analysis of the known and predicted miRNAs were conducted with the DESeq2 R package v1.18.1. Significance was considered at Padj ≤ 0.05. Relationships among production variables and DEMs were explored by Spearman correlations and significance were considered at P ≤ 0.05. Target genes of studied miRNAs and corresponding putative pathways were investigated using the miRWalk database (version 3.0) and Metacore™ software (release 6). Clustering was performed using the reads count for each animal of the DEMs using Permutmatrix software using Euclidean distance calculation and Ward’s minimum variance method [60]. Venn diagram was obtained using Venny 2.1 (

Availability of data and materials

The data supporting the conclusions of this article are within the paper and its additional files. The sequence data generated is publicly available in the NCBI GEO database with the accession GSE131057:




Differentially expressed miRNAs


Epithelial-to-mesenchymal transition


Fold change


Mammary gland






  1. 1.

    Pomies D, Martin B, Chilliard Y, Pradel P, Remond B. Once-a-day milking of Holstein and Montbeliarde cows for 7 weeks in mid-lactation. Animal. 2007;1(10):1497–505.

    CAS  Article  PubMed  Google Scholar 

  2. 2.

    Ferlay A, Martin B, Lerch S, Gobert M, Pradel P, Chilliard Y. Effects of supplementation of maize silage diets with extruded linseed, vitamin E and plant extracts rich in polyphenols, and morning v. evening milking on milk fatty acid profiles in Holstein and Montbeliarde cows. Animal. 2010;4(4):627–40.

    CAS  Article  Google Scholar 

  3. 3.

    Capuco AV, Wood DL, Baldwin R, McLeod K, Paape MJ. Mammary cell number, proliferation, and apoptosis during a bovine lactation: relation to milk production and effect of bST. J Dairy Sci. 2001;84(10):2177–87.

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    Osinska E, Wicik Z, Godlewski MM, Pawlowski K, Majewska A, Mucha J, Gajewska M, Motyl T. Comparison of stem/progenitor cell number and transcriptomic profile in the mammary tissue of dairy and beef breed heifers. J Appl Genet. 2014;55(3):383–95.

    CAS  Article  PubMed  Google Scholar 

  5. 5.

    Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2):281–97.

    CAS  Article  Google Scholar 

  6. 6.

    Ambros V. The functions of animal microRNAs. Nature. 2004;431(7006):350–5.

    CAS  Article  PubMed  Google Scholar 

  7. 7.

    Ameres SL, Zamore PD. Diversifying microRNA sequence and function. Nat Rev Mol Cell Biol. 2013;14(8):475–88.

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Bushati N, Cohen SM, et al. Annu Rev Cell Dev Biol. 2007;23:175–205.

    CAS  Article  PubMed  Google Scholar 

  9. 9.

    Friedman RC, Farh KK, Burge CB, Bartel DP. Most mammalian mRNAs are conserved targets of microRNAs. Genome Res. 2009;19(1):92–105.

    CAS  Article  PubMed  Google Scholar 

  10. 10.

    Timoneda O, Balcells I, Nunez JI, Egea R, Vera G, Castello A, Tomas A, Sanchez A. miRNA expression profile analysis in kidney of different porcine breeds. PLoS One. 2013;8(1):e55402.

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Pacholewska A, Mach N, Mata X, Vaiman A, Schibler L, Barrey E, Gerber V. Novel equine tissue miRNAs and breed-related miRNA expressed in serum. BMC Genomics. 2016;17(1):831.

    Article  PubMed  Google Scholar 

  12. 12.

    Ji Z, Wang G, Xie Z, Wang J, Zhang C, Dong F, Chen C. Identification of novel and differentially expressed MicroRNAs of dairy goat mammary gland tissues using solexa sequencing and bioinformatics. PLoS One. 2012;7(11):e49463.

    CAS  Article  PubMed  Google Scholar 

  13. 13.

    Ji Z, Wang G, Xie Z, Zhang C, Wang J. Identification and characterization of microRNA in the dairy goat (Capra hircus) mammary gland by Solexa deep-sequencing technology. Mol Biol Rep. 2012;39(10):9361–71.

    CAS  Article  Google Scholar 

  14. 14.

    Le Guillou S, Marthey S, Laloe D, Laubier J, Mobuchon L, Leroux C, Le Provost F. Characterisation and comparison of lactating mouse and bovine mammary gland miRNomes. PLoS One. 2014;9(3):e91938.

    Article  PubMed  Google Scholar 

  15. 15.

    Li Z, Liu H, Jin X, Lo L, Liu J. Expression profiles of microRNAs from lactating and non-lactating bovine mammary glands and identification of miRNA related to lactation. BMC Genomics. 2012;13:731.

    CAS  Article  PubMed  Google Scholar 

  16. 16.

    Mobuchon L, Marthey S, Boussaha M, Le Guillou S, Leroux C, Le Provost F. Annotation of the goat genome using next generation sequencing of microRNA expressed by the lactating mammary gland: comparison of three approaches. BMC Genomics. 2015;16:285.

    Article  PubMed  Google Scholar 

  17. 17.

    Avril-Sassen S, Goldstein LD, Stingl J, Blenkiron C, Le Quesne J, Spiteri I, Karagavriilidou K, Watson CJ, Tavare S, Miska EA, et al. Characterisation of microRNA expression in post-natal mouse mammary gland development. BMC Genomics. 2009;10:548.

    Article  PubMed  Google Scholar 

  18. 18.

    Le Guillou S, Sdassi N, Laubier J, Passet B, Vilotte M, Castille J, Laloe D, Polyte J, Bouet S, Jaffrezic F, et al. Overexpression of miR-30b in the developing mouse mammary gland causes a lactation defect and delays involution. PLoS One. 2012;7(9):e45727.

    Article  PubMed  Google Scholar 

  19. 19.

    Tanaka T, Haneda S, Imakawa K, Sakai S, Nagaoka K. A microRNA, miR-101a, controls mammary gland development by regulating cyclooxygenase-2 expression. Differentiation. 2009;77(2):181–7.

    CAS  Article  PubMed  Google Scholar 

  20. 20.

    Li X, Lian F, Liu C, Hu KQ, Wang XD. Isocaloric pair-fed high-carbohydrate diet induced more hepatic steatosis and inflammation than high-fat diet mediated by miR-34a/SIRT1 Axis in mice. Sci Rep. 2015;5:16774.

    CAS  Article  PubMed  Google Scholar 

  21. 21.

    Mobuchon L, Le Guillou S, Marthey S, Laubier J, Laloe D, Bes S, Le Provost F, Leroux C. Sunflower oil supplementation affects the expression of miR-20a-5p and miR-142-5p in the lactating bovine mammary gland. PLoS One. 2017;12(12):e0185511.

    Article  PubMed  Google Scholar 

  22. 22.

    Mobuchon L, Marthey S, Le Guillou S, Laloe D, Le Provost F, Leroux C. Food deprivation affects the miRNome in the lactating goat mammary gland. PLoS One. 2015;10(10):e0140111.

    Article  PubMed  Google Scholar 

  23. 23.

    Wicik Z, Gajewska M, Majewska A, Walkiewicz D, Osinska E, Motyl T. Characterization of microRNA profile in mammary tissue of dairy and beef breed heifers. J Anim Breed Genet. 2016;133(1):31–42.

    CAS  Article  Google Scholar 

  24. 24.

    Peng J, Zhao JS, Shen YF, Mao HG, Xu NY. MicroRNA expression profiling of lactating mammary gland in divergent phenotype swine breeds. Int J Mol Sci. 2015;16(1):1448–65.

    CAS  Article  PubMed  Google Scholar 

  25. 25.

    Ye X, Tam WL, Shibue T, Kaygusuz Y, Reinhardt F, Ng Eaton E, Weinberg RA. Distinct EMT programs control normal mammary stem cells and tumour-initiating cells. Nature. 2015;525(7568):256–60.

    CAS  Article  PubMed  Google Scholar 

  26. 26.

    Xu R, Nelson CM, Muschler JL, Veiseh M, Vonderhaar BK, Bissell MJ. Sustained activation of STAT5 is essential for chromatin remodeling and maintenance of mammary-specific function. J Cell Biol. 2009;184(1):57–66.

    CAS  Article  PubMed  Google Scholar 

  27. 27.

    Bouras T, Pal B, Vaillant F, Harburg G, Asselin-Labat ML, Oakes SR, Lindeman GJ, Visvader JE. Notch signaling regulates mammary stem cell function and luminal cell-fate commitment. Cell Stem Cell. 2008;3(4):429–41.

    CAS  Article  Google Scholar 

  28. 28.

    Brennan KR, Brown AM. Wnt proteins in mammary development and cancer. J Mammary Gland Biol Neoplasia. 2004;9(2):119–31.

    Article  Google Scholar 

  29. 29.

    Yu QC, Verheyen EM, Zeng YA. Mammary Development and Breast Cancer: A Wnt Perspective. Cancers (Basel). 2016;8(7):65.

    Article  PubMed  Google Scholar 

  30. 30.

    Prokesch A, Smorlesi A, Perugini J, Manieri M, Ciarmela P, Mondini E, Trajanoski Z, Kristiansen K, Giordano A, Bogner-Strauss JG, et al. Molecular aspects of adipoepithelial transdifferentiation in mouse mammary gland. Stem Cells. 2014;32(10):2756–66.

    CAS  Article  PubMed  Google Scholar 

  31. 31.

    Perruchot MH, Arevalo-Turrubiarte M, Dufreneix F, Finot L, Lollivier V, Chanat E, Mayeur F, Dessauge F. Mammary epithelial cell hierarchy in the dairy cow throughout lactation. Stem Cells Dev. 2016;25(19):1407–18.

    CAS  Article  Google Scholar 

  32. 32.

    Testa U, Pelosi E, Castelli G, Labbaye C. miR-146 and miR-155: two key modulators of immune response and tumor development. Noncoding RNA. 2017;3(3):22.

    Article  PubMed  Google Scholar 

  33. 33.

    Gao S, Zhao Z, Wu R, Wu L, Tian X, Zhang Z. MiR-146b inhibits autophagy in prostate cancer by targeting the PTEN/Akt/mTOR signaling pathway. Aging. 2018;10(8):2113–21.

    Article  PubMed  Google Scholar 

  34. 34.

    Chen JF, Wu P, Xia R, Yang J, Huo XY, Gu DY, Tang CJ, De W, Yang F. STAT3-induced lncRNA HAGLROS overexpression contributes to the malignant progression of gastric cancer cells via mTOR signal-mediated inhibition of autophagy. Mol Cancer. 2018;17:6.

    Article  PubMed  Google Scholar 

  35. 35.

    Rapa I, Votta A, Gatti G, Izzo S, Buono NL, Giorgio E, Vatrano S, Napoli F, Scarpa A, Scagliotti G, et al. High miR-100 expression is associated with aggressive features and modulates TORC1 complex activation in lung carcinoids. Oncotarget. 2018;9(44):27535–46.

    Article  PubMed  Google Scholar 

  36. 36.

    Sun X, Liu X, Wang Y, Yang S, Chen Y, Yuan T. miR-100 inhibits the migration and invasion of nasopharyngeal carcinoma by targeting IGF1R. Oncol Lett. 2018;15(6):8333–8.

    PubMed Central  PubMed  Google Scholar 

  37. 37.

    Frith JE, Kusuma GD, Carthew J, Li F, Cloonan N, Gomez GA, Cooper-White JJ. Mechanically-sensitive miRNAs bias human mesenchymal stem cell fate via mTOR signalling. Nat Commun. 2018;9(1):257.

    Article  PubMed  Google Scholar 

  38. 38.

    Zarogoulidis P, Lampaki S, Turner JF, Huang H, Kakolyris S, Syrigos K, Zarogoulidis K. mTOR pathway: a current, up-to-date mini-review (review). Oncol Lett. 2014;8(6):2367–70.

    Article  PubMed  Google Scholar 

  39. 39.

    Li R, Dudemaine PL, Zhao X, Lei C, Ibeagha-Awemu EM. Comparative analysis of the miRNome of bovine Milk fat, Whey and Cells. PloS One. 2016;11(4):e0154129.

    Article  PubMed  Google Scholar 

  40. 40.

    Bionaz M, Loor JJ. Gene networks driving bovine milk fat synthesis during the lactation cycle. BMC Genomics. 2008;9:366.

    Article  PubMed  Google Scholar 

  41. 41.

    Chen D, Sun Y, Yuan Y, Han Z, Zhang P, Zhang J, You MJ, Teruya-Feldstein J, Wang M, Gupta S, Hung MC, Liang H, Ma L. miR-100 induces epithelial-mesenchymal transition but suppresses tumorigenesis, migration and invasion. PLoS Genet. 2014;10(2):e1004177.

    Article  PubMed  Google Scholar 

  42. 42.

    Hua X, Xiao Y, Pan W, Li M, Huang X, Liao Z, Xian Q, Yu L. miR-186 inhibits cell proliferation of prostate cancer by targeting GOLPH3. Am J Cancer Res. 2016;6(8):1650–60.

    CAS  PubMed Central  PubMed  Google Scholar 

  43. 43.

    Cao C, Sun D, Zhang L, Song L. miR-186 affects the proliferation, invasion and migration of human gastric cancer by inhibition of Twist1. Oncotarget. 2016;7(48):79956–63.

    Article  PubMed  Google Scholar 

  44. 44.

    Su BB, Zhou SW, Gan CB, Zhang XN. MiR-186 inhibits cell proliferation and invasion in human cutaneous malignant melanoma. J Cancer Res Ther. 2018;14(Supplement):S60–4.

    Article  PubMed  Google Scholar 

  45. 45.

    Li J, Xia L, Zhou Z, Zuo Z, Xu C, Song H, Cai J. MiR-186-5p upregulation inhibits proliferation, metastasis and epithelial-to-mesenchymal transition of colorectal cancer cell by targeting ZEB1. Arch Biochem Biophys. 2018;640:53–60.

    CAS  Article  PubMed  Google Scholar 

  46. 46.

    Qiu H, Yuan S, Lu X. miR-186 suppressed CYLD expression and promoted cell proliferation in human melanoma. Oncol Lett. 2016;12(4):2301–6.

    CAS  Article  PubMed  Google Scholar 

  47. 47.

    Zhu X, Shen H, Yin X, Long L, Xie C, Liu Y, Hui L, Lin X, Fang Y, Cao Y, et al. miR-186 regulation of Twist1 and ovarian cancer sensitivity to cisplatin. Oncogene. 2016;35(3):323–32.

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    Sun P, Hu JW, Xiong WJ, Mi J. miR-186 regulates glycolysis through Glut1 during the formation of cancer-associated fibroblasts. Asian Pac J Cancer Prev. 2014;15(10):4245–50.

    Article  PubMed  Google Scholar 

  49. 49.

    Zhao FQ, Keating AF. Expression and regulation of glucose transporters in the bovine mammary gland. J Dairy Sci. 2007;90(Suppl 1):E76–86.

    Article  PubMed  Google Scholar 

  50. 50.

    Chen Z, Qiu H, Ma L, Luo J, Sun S, Kang K, Gou D, Loor JJ. miR-30e-5p and miR-15a Synergistically Regulate Fatty Acid Metabolism in Goat Mammary Epithelial Cells via LRP6 and YAP1. Int J Mol Sci. 2016;17(11):1909.

    Article  PubMed  Google Scholar 

  51. 51.

    Ma L, Qiu H, Chen Z, Li L, Zeng Y, Luo J, Gou D. miR-25 modulates triacylglycerol and lipid accumulation in goat mammary epithelial cells by repressing PGC-1beta. J Anim Sci Biotechnol. 2018;9:48.

    Article  PubMed  Google Scholar 

  52. 52.

    Naeem A, Zhong K, Moisa SJ, Drackley JK, Moyes KM, Loor JJ. Bioinformatics analysis of microRNA and putative target genes in bovine mammary tissue infected with streptococcus uberis. J Dairy Sci. 2012;95(11):6397–408.

    CAS  Article  Google Scholar 

  53. 53.

    Bernard L, Richard C, Gelin V, Leroux C, Heyman Y. Milk fatty acid composition and mammary lipogenic genes expression in bovine cloned and control cattle. Livest Sci. 2015;176:188–95.

    Article  Google Scholar 

  54. 54.

    Chen CJ, Servant N, Toedling J, Sarazin A, Marchais A, Duvernois-Berthet E, Cognat V, Colot V, Voinnet O, Heard E, et al. ncPRO-seq: a tool for annotation and profiling of ncRNAs in sRNA-seq data. Bioinformatics. 2012;28(23):3147–9.

    CAS  Article  Google Scholar 

  55. 55.

    Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37–52.

    Article  Google Scholar 

  56. 56.

    Bae IS, Seo KS & Kim SH. Identification of endogenous microRNA references in porcine serum for quantitative real-time PCR normalization. Mol Biol Rep. 2018;45:943–9.

  57. 57.

    Lai YC, Fujikawa T, Ando T, Kitahara G, Koiwa M, Kubota C, Miura N. Rapid communication: MiR-92a as a housekeeping gene for analysis of bovine mastitis-related microRNA in milk. J Anim Sci. 2017;95(6):2732–5.

    CAS  Article  Google Scholar 

  58. 58.

    Li D, Liu H, Li Y, Yang M, Qu C, Zhang Y, Liu Y, Zhang X. Identification of suitable endogenous control genes for quantitative RT-PCR analysis of miRNA in bovine solid tissues. Mol Biol Rep. 2014;41(10):6475–80.

    CAS  Article  Google Scholar 

  59. 59.

    Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C(T)) method. Methods. 2001;25:402–8.

    CAS  Article  Google Scholar 

  60. 60.

    Caraux G, Pinloche S. PermutMatrix: a graphical environment to arrange gene expression profiles in optimal linear order. Bioinformatics. 2005;21(7):1280–1.

    CAS  Article  Google Scholar 

Download references


The authors thank the staff at Herbipôle INRA (UE1414, Theix, France), Simon Collange, Julien Mongiat and Denys Durand for performing mammary gland biopsies, Sébastien Bes, Karol Pawlowski, Arnaud Delavaud, and David Chadeyron for technical assistance and for sample collection; Sébastien Bes for RT-qPCR analyses. Sequencing was performed by the IGBMC Microarray and Sequencing platform, a member of the ‘France Génomique’ consortium (ANR-10-INBS-0009).


Billa PA grant was funded by Region & FEDER under the ARTH-INNOV project and this research was funded in part by the GISA program of INRA (LongHealth project, N° 00000290).

Author information




CL, YF, and JP conceived, designed and initiated the project. PAB, TY, YF, and MC performed analyses. CL and PAB wrote the original draft. JP, YXF, and FLP helped to revise the manuscript. CL supervised the project. All authors read and approved the final version.

Corresponding author

Correspondence to C. Leroux.

Ethics declarations

Ethics approval and consent to participate

Animal procedures were performed in compliance with the Regional Animal Care Committee guidelines CEMEAA: Auvergne, French Ministry of Agriculture and European Union guidelines for animal research C2EA-02. All procedures were approved by the ethics committee on animal experimentation (under n° 3737-2015043014541577 V2).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional file

Additional file 1:

Table S1. Individual analysis of the expression of four miRNA in mammary gland. RT-qPCR analyses were performed from 11 cows (5 Holstein (H) and 6 Montbéliarde (M)). All miRNA levels were normalized to the values of 3 internal control miRNAs (miR-26a, −92a, − 191) and the results expressed as fold changes of threshold cycle (Ct) values relative to the control using the 2-ΔΔCt method. Results are presented as log2 ratio between M/H. TaqMan advance references are given for ech system. Table S2. Relationships among milk, fat, protein and lactose yields and 6 discussed differentially expressed miRNAs. The relationships were explored by Spearman correlations using 5 Holstein and 6 Montbéliarde cows. *** P < 0.001. ** P < 0.01. * P < 0.05. Figure S1. Comparison of percentage of match miRNA reads between libraries. Libraries were constructed from total RNA from mammary gland biopsies of 5 Holstein (H) and 6 Montbéliarde (M) cows. Blue correspond to percentage of unique-mapped, orange to multi-mapped and grey to unmapped reads. Figure S2. Correlation between Holstein and Montbéliarde miRNAs libraries. The correlation was calculated using log2 of the normalized counts of expressed miRNAs. Figure S3. Choice of 3 internal controls. Expression stability values (M) and rankings of the reference genes are determined by geNorm software. (DOCX 208 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Billa, P.A., Faulconnier, Y., Ye, T. et al. Deep RNA-Seq reveals miRNome differences in mammary tissue of lactating Holstein and Montbéliarde cows. BMC Genomics 20, 621 (2019).

Download citation


  • MiRNome
  • Mammary gland
  • Dairy cows
  • Breeds
  • RNA-Seq
  • Differentially expressed