Comprehensive MicroRNA Expression Profile of the Mammary Gland in Lactating Dairy Cows With Extremely Different Milk Protein and Fat Percentages

Abstract

A total of 31 differentially expressed genes in the mammary glands were identified in our previous study using RNA sequencing (RNA-Seq), for lactating cows with extremely high and low milk protein and fat percentages. To determine the regulation of milk composition traits, we herein investigated the expression profiles of microRNA (miRNA) using small RNA sequencing based on the same samples as in the previous RNA-Seq experiment. A total of 497 known miRNAs (miRBase, release 22.1) and 49 novel miRNAs among the reads were identified. Among these miRNAs, 71 were found differentially expressed between the high and low groups (p < 0.05, q < 0.05). Furthermore, 21 of the differentially expressed genes reported in our previous RNA-Seq study were predicted as target genes for some of the 71 miRNAs. Gene ontology and KEGG pathway analyses showed that these targets were enriched for functions such as metabolism of protein and fat, and development of mammary gland, which indicating the critical role of these miRNAs in regulating the formation of milk protein and fat. With dual luciferase report assay, we further validated the regulatory role of 7 differentially expressed miRNAs through interaction with the specific sequences in 3′UTR of the targets. In conclusion, the current study investigated the complexity of the mammary gland transcriptome in dairy cattle using small RNA-seq. Comprehensive analysis of differential miRNAs expression and the data from previous study RNA-seq provided the opportunity to identify the key candidate genes for milk composition traits.

Background

MicroRNAs (miRNAs), which are a class of non-coding small RNA (sRNA) molecules with the length of 18-24 nucleotides, are important regulators of gene expression. They can play important roles in a wide range of biological processes, including animal and plant development, cell differentiation, proliferation, apoptosis, and metabolism (Martello et al., 2010; Chen et al., 2012; Rottiers and Näär, 2012; Almughlliq et al., 2019; Barbu et al., 2020). In animal cells, miRNAs interact with a specific sequence in mRNA of the target gene and post-transcriptionally negatively regulate the expression of target genes by inhibiting their translation or inducing degradation of the target mRNAs (Huntzinger and Izaurralde, 2011; Barbu et al., 2020). MiRNAs have emerged as new potential biomarkers for miRNA-gene interactions and gene networks responsible for human diseases and economically important traits in livestock. Several diseases and conditions have been reported to be linked with the abnormal expression in miRNAs relating with differentiation, apoptosis and development (Lewis et al., 2005; Berezikov et al., 2006b; Lee et al., 2007). Many experimental techniques and computational methods have been developed to identify miRNAs (Aravin and Tuschl, 2005; Berezikov et al., 2006a; Landgraf et al., 2007), and large number of miRNAs have been identified in primates, rodents, birds, fish, and plants (Lagos-Quintana et al., 2003; Chen et al., 2005; Finucane et al., 2008; Glazov et al., 2008).

The bovine mammary gland is a complex organ which grows and develops after calving and is able to produce more than 30,000 kg of milk in a complete lactation cycle (Hennighausen and Robinson, 2005; Muroya et al., 2019). Because of its important functions, the mammary gland, especially mammary epithelial cells, has been used as the target tissue for gene expression profiling in order to identify key genes underlying milk production traits in dairy cattle (Silveri et al., 2006; Bionaz and Loor, 2008, 2011; Bionaz et al., 2012; Zhang et al., 2016; Pu et al., 2017; Cai et al., 2018; Ju et al., 2018; Yang et al., 2018; Billa et al., 2019; Li et al., 2020). However, only a few studies have been reported related to the miRNAs in the bovine mammary gland. A total of 798 mature bovine miRNAs have been deposited in miRBase (Luoreng et al., 2018), Release 22.1 (October 2018) and 55 of them were detected in the mammary gland. Li et al. (2012b) reported 283 known miRNAs and 74 novel miRNAs in the mammary gland of Holstein cows, among which 56 miRNAs were differentially expressed between lactating and non-lactating cows and might be involved in regulating lactation. Shen et al. (2016) identified 292 known miRNAs and 116 novel miRNAs in the bovine mammary epithelial cells, and three of them (bta-miR-33a, bta-miR-152 and bta-miR-224) might be involved in milk fat metabolism. Li et al. (2015) detected 370 known and 341 novel miRNAs in the bovine mammary gland infected with Staphylococcus aureus, and 358 known and 232 novel miRNAs in control group, 77 of which were differentially expressed between infected and healthy Holstein cows. In addition, Le Guillou et al. (2012) found that the overexpression of miR-30b caused a defect in lactation and delayed involution in mouse mammary gland.

In a previous study from our lab (Cui et al., 2014), 31 differentially expressed genes was identified by using RNA sequencing (RNA-Seq) to investigate the mammary gland epithelial tissues of four lactating Holstein cows with extremely high and low milk protein (PP) and fat percentages (FP). The objectives of the present study were to investigate the miRNA expression profiles in the same mammary gland samples that were used in the previous RNA-Seq study to identify known and novel miRNAs, and to perform an analysis of the differentially expressed miRNAs and previously identified genes. Some candidate miRNAs and their target genes that may be involved in milk protein and fat metabolism were identified.

Materials and Methods

Animals and Mammary Gland Tissue Samples

In the current study, the mammary gland epithelium samples of four lactating Chinese Holstein cows (high group vs. low group) same as our previous RNA-Seq experiment (Cui et al., 2014) were used. These four cows were selected from 30,000 Holstein cows in Beijing Sanyuanlvhe Dairy Farming Center, and the average PP and FP were 3.1% (2.7–3.8%) and 3.6% (3.1–4.5%) in this population. In order to keep the environmental factors identical, these four cows in almost the same period of lactation (353, 341, 377, and 325 days) were collected from the same farm possessing a total of 800 Holstein cows. Selected cows were divided into two groups according to the phenotypic values for PP and FP: two cows (high group) had high PP (3.6% and 3.8%) and FP (3.9% and 4.5%); the other two cows (low group) showed low PP (3.0%, 2.9%) and FP (3.2%, 3.1%).

The cows were killed by electroshock, and then they were bled, skinned, and dismembered in the same slaughterhouse. The rear mammary gland from each individual was harvested within 30 min after slaughtered. White mammary ducts and pink epithelium tissue were clearly observed when the right rear quarter of the mammary gland was cut in half lengthways from the teat and some milk were flowed out. Five pieces of epithelium tissue samples per cow were carefully collected and placed into a clean RNAse-free Eppendorf tube, and then stored in liquid nitrogen for subsequent RNA isolation. All procedures of collecting samples were carried out in strict accordance with the protocol approved by the Animal Welfare Committee of China Agricultural University (Permit Number: DK996). Total RNA was extracted from one piece of mammary gland epithelium samples from each cow and quality was controlled according to the protocols described by Cui et al. (2014). The value of RNA integrity number (RIN) from each sample was above 8.0.

Small RNA Sample Preparation and Sequencing

The preparation of small RNA library, including quality control and sequencing, was performed by Novogene (Beijing, China). The preparation of library was performed on 3 μg total RNA per sample using an IlluminaTruSeq™Small RNA Sample Preparation Kit (Illumina, San Diego, CA, United States). The samples were indexed using four codes in order to facilitate sequencing of these samples on one flow cell channel. Quality control in library preparation showed that adapter-adapter contamination was <5% and 85% of the sequences were miRNAs. The samples were subsequently sequenced on the Illumina Hiseq2000 platform and 50-bp single-end reads were obtained.

Sequencing Data Analysis

The sequencing data were obtained in the format of Illumina FASTQ (Illumina). The procedure of data filtering included removing low quality reads, reads containing poly-N stretches, reads with 5′primer contaminants, reads with 3′primers or the insert tag, and reads with poly-A, T, G, or C stretches. Thereafter, the sRNA tags within a certain range (18-30 nt) were retained for the successive steps. The Q20, Q30, and GC-content of the cleaned reads were calculated to evaluate the quality of data. Then, the sRNA tags were mapped to the bovine genome assembly (UMD3.1.66) using Bowtie (Langmead et al., 2009), no mismatches were allowed and the “seed” region size was set at 8 (Gupta et al., 2012; Giurato et al., 2013; Aggarwal et al., 2014; Kuksa et al., 2018). The mapped sRNA tags were aligned to the 798 bovine miRNA precursor sequences in miRBase (Release 22.1) to identify the known miRNA in the sRNA libraries allowing one mismatch. The sRNA tags that matched known miRNAs from species other than bovine may be novel bovine miRNAs, and were predicted the secondary structure, the Dicer cleavage site, and the minimum free energy of the mapped sRNA sequences using the miREvo (Wen et al., 2012) and miRDeep2 (Friedländer et al., 2012) software packages.

The expression of miRNA was measured as counts per million (CPM) using the following formula: normalized expression = mapped read count/total reads × 1000000 (Zhou et al., 2010), and DESeq2 R package (1.8.3) (Anders and Huber, 2010; Trapnell et al., 2013) was used to identify significantly differentially expressed miRNAs between high and low groups of cows. The threshold for differential expression was—log2 (FC)— > 1 and FDR p < 0.05 when using DESeq2 R package for differential expression miRNA analysis so that miRNAs with—log2 (FC)— > 1 and adjusted FDR p < 0.05 were designated as differentially expressed.

Furthermore, two cows in the same group were used to eliminate the background noise of individual-specific transcription by applying a pairwise approach, which enabled acquisition of more relevant data from the two groups.

Target Prediction, Pathway, and Annotation Analysis

TargetScan 6.2 and MiRanda (Enright et al., 2003) were used to predict putative target genes with the established miRNA seed database and the bovine genome sequence (UMD3.1.66). TargetScan 6.2 predicts targets by searching for the presence of conserved 8mer, 7mer, and 6mer sites that match the seed region of each miRNA. MiRanda predicts targets based on a development of the miRanda algorithm which incorporates current biological knowledge on target rules and on the use of an up-to-date compendium of mammalian miRNAs.

Gene ontology (GO) functional enrichment analysis was used for the candidate target genes of the miRNAs. GOseq with the Wallenius non-central hyper-geometric distribution (Young et al., 2010), which can adjust for the bias in gene length, was implemented for the GO enrichment analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008) pathways analysis was performed using KOBAS 2.0 (Mao et al., 2005) software to test the statistical enrichment of the candidate target genes in the KEGG pathways.

Quantitative Real Time PCR

Expression levels of selected miRNAs were confirmed by quantitative real-time PCR (qRT-PCR) using the DyNAmo SYBR Green PCR kit (Applied Biosystems, Foster City, CA, United States) on a LightCycler480 (Roche Applied Science, Penzberg, Germany). qRT-PCR of target mRNAs was performed using specific miRNA stem-loop primers (Supplementary Table 8) and all reactions were run in triplicate. Relative quantification of miRNA was quantified using the 2^–ΔΔCT method and normalized against the U6 gene (ssD0904071006: Guangzhou RiboBio, Guangzhou, China) for each sample.

Plasmid Construction and Site-Directed Mutagenesis of 3′UTR in Predicted Target Genes

The 3′un-translated region (UTR) of four predicted target genes for the identified miRNAs, TRIB3, M-SAA3.2, PTHLH, and VEGFA, were PCR amplified using DNA collected from the bovine mammary gland samples applied for sequencing in this study as a template, and connected into pmirGLO Dual-Luciferase miRNA Target Expression Vector (pmirGLO, Promega) (Figure 1), respectively. The primers were listed in the Table 1. Afterward, the connected products were transfected into Escherichia coli, and then verified the correct sequence and orientation by sequencing. The QuikChange site-directed mutagenesis kit (Stratagene, La Jolla, CA, United States) was used to generate the 3′UTR variants of TRIB3, M-SAA3.2, PTHLH, and VEGFA where seed sequences recognized by microRNAs were deleted (Figures 2, 3). After the point mutation, same way was applied in order to find the correct mutant sequences for such four genes.

FIGURE 1

FIGURE 2

FIGURE 3

Luciferase Reporter Assays

To further explore the repressing mechanism of miRNAs on the expression of 4 target genes (TRIB3, M-SAA3.2, PTHLH, and VEGFA) expression, the full-length TRIB3, M-SAA3.2, PTHLH and VEGFA 3′UTRs and the corresponding mutant version (the seed sequences were deleted) were transfected into human embryonic kidney HEK293 cells (GM-070001H: Shanghai, China), respectively. These cells were cultured at 37°C with 5% CO₂ in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 4.5 g/liter glucose, 5% fetal bovine serum (Invitrogen), 2 mmol/liter glutamine, and antibiotics. Before transfection, HEK293 cells were plated into 24-well plates at 1.0 × 10⁵ cells/well 24 h. 30 ng empty pmirGLO vector, pmirGLO-TRIB3/M-SAA3.2/PTHLH/VEGFA-3′UTR with 50 μl opti-MEM (Invitrogen) and 30 nM (final concentration) mimic miRNA, inhibitor miRNA, control miRNA (GenePharma) were co-transfected into each well with 1 μl Lipofectamine 2000 (Invitrogen). 30 ng mutants of the TRIB3/M-SAA3.2/PTHLH/VEGFA 3′UTR with 50 μl opti-MEM (Invitrogen) and 30 nM (final concentration) mimic miRNA, control miRNA (GenePharma) were co-transfected into each well with 1 μl Lipofectamine 2000 (Invitrogen). Relative firefly luciferase activities (normalized to Renilla luciferase activities) were measured 24 h after transfection with the Dual-Luciferase Reporter Assay Kit (Promega) on TECAN Infinite 200 multifunctional microplate reader (TECAN). All experiments were performed in triplicate so that data averaged from three independent experiments.

Results

Sequencing and Mapping of the sRNA Tags

Four new miRNA libraries were constructed using sRNA isolated from bovine mammary glands and sequenced using Illumina next-generation sequencing. A total of 10,538,878 (high milk PP and FP), 12,745,512 (high milk PP and FP), 9,744,027 (low milk PP and FP), and 9,682,136 (low milk PP and FP) high-quality cleaned reads were obtained from the four sRNA libraries (Supplementary Table 1; NCBI SRA accession numbers: SRR3631014, SRR3631016, SRR3631053, and SRR3631054). Distribution of the length for reads showed that most of the generated reads had 21 (>24%), 22 (>30%), and 23 (>13%) nucleotides (Supplementary Figure 1), which is the size of most known mature miRNAs. When aligning the sequenced reads against the bovine genome assembly (UMD3.1.66), it was found that 77.57%, 76.93%, 80.88%, and 78.15% of them uniquely aligned from the four libraries, respectively (Supplementary Table 1); 55-57% of them were aligned in the same direction as the reference genome sequence, and 20-25% were aligned in the opposite direction (Supplementary Table 2). The correlation coefficient (R²) between the two individuals within the high and low groups for milk PP and FP was calculated based on the CPM mapped fragment of each cow and was shown to be 0.988 and 0.980, respectively. This indicated that the similarity of the two biological replicates within each group was sufficiently high (Supplementary Figure 2).

MicroRNAs Identification and Target Prediction

Among the uniquely aligned reads across the four samples and six downloaded miRNA libraries (Cai et al., 2018), 24,320,809 (54.4%) matched known miRNAs in miRBase (Release21.0), which resulted in 497 known bovine miRNAs and 49 novel bovine miRNAs were identified (Supplementary Tables 3, 4). Subsequently, two well-established target prediction tools, TargetScan and miRanda, were used to predict target mRNAs of the miRNAs, and a total of 12,202 target genes were commonly predicted for the known and novel miRNAs (Supplementary Table 5). It is noteworthy that some well-known genes associated with milk composition traits were included such as β-casein (CSN2), κ-casein (CSN3), α-lactalbumin (LALBA), diacylglycerol O-acyltransferase 2 (DGAT2), growth hormone receptor (GHR), signal transducer and activator of transcription 5B (STAT5B), and stearoyl-coenzyme A desaturase (SCD) etc. This finding implied that the identified mammary miRNAs in this study were involved in metabolism of milk protein and lipid through the regulation of key genes affecting these traits.

Differentially Expressed miRNAs Between the High and Low Groups for Milk PP and FP and Target Prediction

The miRNAs that differed between the high and low PP and FP groups were determined in this study. A total of 71 top half miRNAs displayed significantly differential expression between the high and low groups using the DEseq2 algorithm (p < 0.05, FDR q < 0.05), with 35 were up-regulated and 36 were down-regulated in the high milk PP and FP group compared with the low group (Table 2). Subsequently, a total of 5,634 target genes were commonly obtained for these differentially expressed miRNAs by TargetScan and miRanda (Supplementary Table 7).

Afterward, the results of the sequencing were validated with an independent method of real-time PCR assay. By using the same four mammary gland samples as used for sequencing, eight known miRNAs and seven novel miRNAs identified in the present study were randomly chosen for validation. It was found that the expression levels of miR-125a, miR-2904, miR-345-5p, miR-378c and Novel-18 were significantly higher in the high milk PP and FP group than in the low group (p < 0.05), and the expression levels of miR-21-3p, miR-29c, miR-106b and miR-190a were lower in the high group than in the low group (p < 0.05). Whereas, Novel-13, Novel-2, Novel-22, Novel-32, Novel-4 and Novel-42 did not display significant differences on miRNA levels between the two groups (p > 0.05) (Figure 4). Such expression patterns were exactly consistent with those shown by small-RNA sequencing data.

FIGURE 4

Gene Ontology Enrichment and Pathway Analysis

To further investigate the functional associations of the target genes, gene ontology (GO) annotation analysis was performed. It was found that these targets have a wide range of diverse functions, among which most were involved in protein and lipid metabolism, mammary gland development and differentiation, and immune functions (p < 0.01, FDR q < 0.01). Under the GO biological process category, the enriched terms related to lipid and protein metabolisms and cell growth were included such as protein binding, protein localization, protein transport, protein complex, regulation of protein metabolic process, lipid biosynthetic process, programmed cell death, protein targeting, lipid metabolic process, amino acid transport, regulation of protein kinase activity, and cellular response to mechanical stimulus (Supplementary Table 6).

A KEGG metabolic pathway analysis was also performed to identify functions that associate with the predicted target genes using KOBAS. Targets were enriched for functions such as mitogen-activated protein kinases (MAPK), adipocytokine, mammalian target of rapamycin (mTOR), glycosphingolipid biosynthesis, glycerophospholipid metabolism, hypoxia inducible factor-1 (HIF-1), and phosphatidylinositol 3 kinase-protein kinase B (PI3K-Akt) signaling pathways (Table 3).

For the 71 top half differentially expressed miRNAs, 5,634 target genes were obtained and the targets were highly enriched in biological process consisting of synthesis and metabolism of protein and energy metabolism, as well as pathways mainly related to synthesis and metabolism of lipid and protein including glutathione metabolism, NF-kappa B signaling pathway, mTOR signaling pathway, fatty acid degradation, fatty acid metabolism and protein processing in endoplasmic reticulum (Table 4).

Comparison of the Target Genes of the Differentially Expressed miRNAs and the Differentially Expressed Genes Reported Previously

In our previous study (Cui et al., 2014), 21 of the target genes, which are listed in Table 5, were found to be differentially expressed between the high and low groups using the same four mammary gland samples in the current study. Among the 21 differentially expressed target genes, the expressions of only six down-regulated genes and one up-regulated gene matched the expression profiles of the differentially expressed miRNAs that targeted them. While 5 down-regulated genes were targeted by at least one up-regulated miRNA each, and 10 genes were targeted by both up-regulated and down-regulated miRNAs. Especially, 7 of the 21 differentially expressed target genes were the most promising candidate genes affecting milk protein and fat percentage identified by integrated analysis of differential gene expression, previously reported quantitative trait loci (QTLs) and genome-wide association studies (GWAS) (Cui et al., 2014), including tribbles homolog 3 (TRIB3), serum amyloid A1 (SAA1), serum amyloid A3 (SAA3), mammary serum amyloid A3 (M-SAA3.2), vascular endothelial growth factor A (VEGFA), parathyroid hormone-like hormone (PTHLH) and ribosomal protein L23A (RPL23A). In addition, KEGG pathway analysis using KOBAS, showed that two of the 21 target genes, DNA-damage-inducible transcript 3 (DDIT3) and nuclear receptor subfamily 4, group A, member 1 (NR4A1), were involved in the MAPK signaling pathway that plays critical role in protein synthesis and metabolism and fatty acid metabolism pathway (p < 0.05; Tables 3, 4), and 2 other genes, vascular endothelial growth factor A (VEGFA) and cyclin-dependent kinase inhibitor 1A (CDKN1A), were involved in the mTOR, HIF-1, PI3K-Akt, p53 and duct acid secretion signaling pathways, which are mostly related to synthesis and metabolism of protein and fat (p < 0.05; Tables 3, 4).

MicroRNAs Repress the Expression of Target Genes Through the Binding of a Specific Target Sequence in Their mRNA 3′UTR

To study the regulatory functions of the identified miRNAs, four differentially expressed genes were chosen including TRIB3, M-SAA3.2, PTHLH and VEGFA, which having the expression pattern negatively correlated with their targeting miRNAs. Using the dual luciferase reporter assays, whether miR-2904, miR-339b/miR-146b/miR-339a, miR-29c/miR-106b/miR-190a, and miR-2904/miR-106b/miR-21-3p regulated the expression of the TRIB3, M-SAA3.2, PTHLH and VEGFA, respectively, were detected. Consequently, it was found that the luciferase level in HEK293 cells with mimics of miR-2904 decreased 40% relative to those with the empty vector, respectively (p < 0.05), while the inhibitor of miR-2904 yielded the same luciferase level as negative control (p > 0.05) (Figure 5). However, when the predicted binding sites of such miRNA seed sequences were mutated, luciferase activity was efficiently restored to the control levels (p > 0.05; Figure 5). Such results clearly indicated the notable regulatory role of the miR-2904 on the expression of TRIB3 by directly targeting its 3′UTR. Similarly, with regard to M-SAA3.2, it was also found that the overexpression of miR-146b, miR-339a and miR-339b decreased the luciferase levels in HEK293 cells by 80%, 72% and 74% after transfecting these mimics compared with the negative controls, respectively (p < 0.05), and the depressed expression of such miRNAs did not change the luciferase level in HEK293 cells transfected with their inhibitors (p > 0.05), respectively (Figure 6). When the mutant 3′UTR of M-SAA3.2 and mimics of the 3 miRNAs were co-transfected, the luciferase activity was same as the control level (p > 0.05; Figure 6). For PTHLH, the luciferase level in HEK293 cells transfected with the mimics of miR-29c, miR-106b and miR-190a was decreased by 37%, 49%, and 50% relative to the negative control, respectively (p < 0.05; Figure 7), however, the same level was kept by transfecting the inhibitors of such miRNAs (p > 0.05), respectively (Figure 7). When the mutant version of the PTHLH 3′UTR and mimics of miR-29c, miR-106b and miR-190a were co-transfected, respectively, the luciferase activities were same as the control levels (p > 0.05; Figure 7). Whereas, the expression of VEGFA was not affected by miR-2904, miR-106b and miR-21-3p (p > 0.05, Figure 8).

FIGURE 5

FIGURE 6

FIGURE 7

FIGURE 8

Discussion

The current study is the first comparative profiles of the mRNA and miRNA transcriptome in the mammary gland epithelium of dairy cows to the best of our knowledge. In this study, we generated an extensive miRNA expression profile of the mammary glands from lactating cows with extremely high and low milk PP and FP, and identified a total of 497 known bovine miRNAs and 49 novel bovine miRNAs. In previous studies, bovine miRNAs were identified using computational and direct cloning approaches (Coutinho et al., 2007; Gu et al., 2007; Jin et al., 2009, 2010; Long and Chen, 2009; Li et al., 2012b, 2015; Shen et al., 2016). Li et al. (2012b) identified 298 known miRNAs in lactating and non-lactating mammary gland of Holstein cows using miRNA-seq; 204 of them were among the 497 known miRNAs identified in the current study. Furthermore, 9 of the 71 differentially expressed miRNAs (miR-100, miR-10a, miR-133a, miR-1, miR-146b, miR-148a, miR-221, miR-30f, and miR-339b) identified in the current study were also reported by Li et al. (2012b) as differentially expressed between lactating and non-lactating bovine mammary glands. Gu et al. (2007) identified 31 distinct miRNAs in the mammary glands of Holstein cows, and all of these miRNAs was detected in the present study except miR-142b. Shen et al. (2016) identified 292 known miRNAs in the bovine primary mammary cells, among which 217 miRNAs and 38 differentially expressed miRNAs were also identified in the current study. For the 30 differentially expressed miRNAs in the lactating goat mammary gland fed ad libitum or deprived of food affecting milk composition reported by Mobuchon et al. (2015), only 6 miRNAs, including miR-660-5p, miR-451-5p, miR-125b, miR-196a, miR-223-3p, and miR-223-5p were detected as well in the current study. miR-30b related to lactation in mouse (Le Guillou et al., 2012) was also detected in this study, but did not show differential expression between high and low groups. The reason could be due to the mammary gland tissues were collected from different time points of lactation between the previous (Le Guillou et al., 2012) and the current studies.

miR-15a has been reported to be critical in cell development (Bonci et al., 2008), cell cycle (Bandi et al., 2009), and death (Cimmino et al., 2005; Aqeilan et al., 2010). Li et al. (2012a) found that miR-15a can inhibit the viability of mammary epithelial cells as well as the mRNA and protein expression of GHR, which is a major gene for milk composition traits (Bonci et al., 2008). In the current study, we also detected miR-15a and predicted that it may target GHR as well as candidate genes for milk PP and FP identified in our previous study, namely activating transcription factor 3 (ATF3), VEGFA, parathyroid hormone-like hormone (PTHLH), cation transport regulator homolog 1 (CHAC1), and NR4A1. Therefore, miR-15a was considered may affect milk composition by regulating the expression of these genes, although miR-15a was not one of the differentially expressed miRNA identified in this study. It was reported that miR-23b inhibited the expression of the transforming growth factor-beta (TGF-β) signaling (Finnerty et al., 2010). In the current study, miR-23b and 5 other miRNAs (miR-2454-3p, miR-496, miR-503-3p, miR-6520, and novel-6) were predicted to regulate STAT5B, which is known to be involved in TGF-β signaling (Passerini et al., 2008; Hosui et al., 2009). In addition, genes that are known to affect milk traits (CSN3, CSN2, LALBA, DGAT2, STAT5B, and SCD) were predicted to be targets of some of the identified miRNAs, which implied that they may play critical regulatory roles in mammary gland development and milk composition.

It was found that 21 of 31 differentially expressed genes detected in our previous study (Cui et al., 2014) were the predicted targets for some of the 71 differentially expressed miRNAs detected in the present study. Serum amyloid A1 (SAA1), serum amyloid A1 (SAA3), and mammary serum amyloid A3.2 (M-SAA3.2) were predicted to be regulated by miR-146b (SAA1 was also regulated by miR-125a and miR-125b); VEGFA was regulated by miR-125a, miR-125b, miR-106b, and miR-2904; and ribosomal protein L23a(RPL23A), tribbles homolog 3 (TRIB3), and PTHLH were regulated by miR-378c, miR-2904, and miR-106b, respectively. Moreover, Cai et al., performed RNA sequencing with mammary gland tissue samples from six Chinese Holstein cows with three extremely high and three low milk protein percentage phenotypes and miR-2904, miR-339b, miR-146b, miR-339a, miR-29c, miR-106b, miR-190a, miR-21-3p, miR-15a, miR-486, miR-135, miR-101a, miR-152 and miR-139 were found differentially expressed, which were also identified in our study and targeted on four differentially expressed genes (TRIB3, PTHLH, VEGFA, and M-SAA3.2). These seven genes represent the most promising candidates may affect milk PP and FP in dairy cattle (Cui et al., 2014). Specifically, miR-146b was reported to be involved mainly in leukemia, epidermal growth factor receptor (EGFR), MAPK, and nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) signaling pathways (Mathews et al., 2004; Taganov et al., 2006; Xiang et al., 2014). The EGFR and MAPK signaling pathways have been demonstrated to be related to adipocyte differentiation (Devaraj et al., 2009; Gao and Bing, 2011) and the NF-κB pathway controls the DNA transcription protein complexes. In human study, miR-146b was shown to regulate the NF-κB signaling pathway in which breast cancer metastasis suppressor 1 (BRMS1) has already been implicated, and inhibited both migration and invasion related to metastasis (Taganov et al., 2006; Xiang et al., 2014). Members of the miR-125 family were reported to be implicated in a variety of carcinomas and other diseases as either repressors or promoters. Sun et al. (2013) found that up-regulated miR-125 significantly inhibited the expression of VEGFA both in vitro and in vivo (Jiang et al., 1997). The miR-125 family was found to be a NF-κB-dependent gene in the study by Kim et al. (2012). miR-378c was shown to be involved in the regulation of RPL23A, which plays a critical role in translation and participates in apoptosis, cell division, and differentiation (Wool, 1996; Fang et al., 2012; Knezevic et al., 2012). This is consistent with previous reported study where miR-378c was found associated with apoptosis (Lee et al., 2007; Fang et al., 2012; Knezevic et al., 2012; Wang et al., 2014).

The GO and KEGG pathway analyses indicated that VEGFA, NR4A1, DDIT3, and CDKN1A were involved in the MAPK, mTOR, HIF-1, and PI3K-Akt signaling pathways, respectively. These four genes were predicted as target genes for miR-106b, miR-2904, miR-125a(b), miR-21-3p, miR-224, miR-31, miR-345-5p, and miR-3431. mTOR signaling is known as playing a fundamental role in adipogenesis (Laplante and Sabatini, 2009), which is the process that leads to the formation of adipose tissue and the most important energy storage site in mammals. It has been demonstrated that mTORC1 positively regulates the activity of sterol regulatory element binding protein 1 (SREBP1) and peroxisome proliferator-activated receptor gamma (PPARG) (Benmoussa et al., 2020), which are two transcription factors that control the expression of genes encoding proteins involved in lipid and cholesterol homeostasis (Kim and Chen, 2004; Porstmann et al., 2008; Kim et al., 2012). HIF-1 is a heterodimeric transcription factor that increases the phosphorylation of signal transducer and activator of transcription 5A (STAT5A) in mammary epithelial cells, and the phosphorylation of STAT5 is known to play important roles in the regulation of milk protein gene expression and mammary development (Shao and Zhao, 2014; Benmoussa et al., 2020). Several studies have shown that hypoxia causes mammary epithelial disorganization and induces a cancer cell-like phenotype in human mammary epithelial cells (MECs) (Whelan et al., 2010; Whelan and Reginato, 2011; Vaapil et al., 2012). The PI3K-Akt pathway has important functions in mammary gland development and function (Wickenden and Watson, 2010). One of the most important functions of Akt is the regulation of glucose homeostasis and metabolism, particularly in muscle and fat tissues (Enright et al., 2003). Therefore, these miRNAs could play critical roles in regulating formation of milk composition trait.

Considering that microRNAs regulate gene expression by targeting specific sequences in the 3′UTR of their cognate genes (Lewis et al., 2005; Friedman et al., 2009), the regulatory roles of some miRNAs on their predicted targets were verified using dual luciferase report assay transfected with mimics, inhibitors and mutants of seed sequences. The results demonstrated that miR-2904, miR-29c/miR-146b/miR-339a, miR-339b/miR-106b/miR-190a indeed down-regulated the expression of the TRIB3, M-SAA3.2 and PTHLH, respectively. The molecular mechanisms of how these miRNAs regulate their targets will be further validated through RNAi and over-expression in bovine mammary epithelial cell lines. In addition, it is generally recognized that miRNAs regulate the expression of target genes by inhibiting their translation or inducing degradation of the target mRNAs in animal cells. However, several predicted target genes were regulated in the same direction of expression as those of the corresponding miRNAs between high and low groups. The reason could be due to either target prediction error of the current commonly used prediction softwares (TargetScan 6.2 and MiRanda) or some unknown biological mechanisms. Actually, target prediction was only the first step for studies on interaction between miRNA and their targets. The miRNAs and targets with reverse expression patterns will be considered as the key components for further validation.

In this study, only two biological replicates, which were the same as in our previous RNA-seq investigation (Cui et al., 2014), were used for each condition due to the availability of mammary gland sample from lactating cows, especially high production ones. In order to minimize false-positive errors and ensure substantial detection power and accuracy, two strategies were applied to detect the differentially expressed miRNAs between milking Holstein cows with high PP and FP and cows with low PP and FP, by controlling the critical influencing factors. Small RNA transcripts were deeply sequenced (9-10G data per transcriptome), and only those differentially expressed miRNAs ranked in the top half of the expressed miRNAs were considered, as suggested by Rapaport et al. (2013), Trapnell et al. (2013). Rapaport et al. (2013) investigated the impact of different sequencing depths and number of replicates on the identification of differentially expressed genes, where the authors demonstrated that with most methods, over 90% of differently expressed genes at the top expression levels could be detected with using two replicates and 5% of the reads (Rapaport et al., 2013; Trapnell et al., 2013). The differentially expressed miRNAs expressed in the bottom half level were eliminated to ensure the power in detection. Although mRNA sequencing data was used in this study, detection of differentially expressed miRNAs is based on same statistical theory and software (Anders and Huber, 2010). However, more biological replications are still preferred and recommended in order to provide broader application (Rapaport et al., 2013; Trapnell et al., 2013). The more replicates are performed, the more the detection power is improved. The potential regulatory roles on target genes from such differentially expressed miRNAs will be validated further by performing more in-depth investigation.

Conclusion

Using sRNA sequencing, 497 known bovine miRNAs and 49 novel bovine miRNAs were identified in the mammary glands of lactating dairy cows. Among all these miRNAs, 71 were differentially expressed between cows with the high and low milk PP and FP. Combined with our previous RNA-Seq data, 21 differentially expressed genes were predicted as the targets for some of the 71 differentially expressed miRNAs. Biological processes related to protein metabolism, fat metabolism, and mammary gland development were enriched for some of the identified miRNAs, which indicated that they may play critical roles in regulating of milk protein and fat traits in dairy cattle. Expression of the TRIB3, M-SAA3.2 and PTHLH were significantly down-regulated by miR-2904, miR-29c/miR-146b/miR-339a and miR-339b/miR-106b/miR-190a through binding the specific target sequences in 3′UTR of these genes, respectively.

References

• 1

  AggarwalP.TurnerA.MatterA.KattmanS. J.StoddardA.LorierR.et al (2014). Rna expression profiling of human ipsc-derived cardiomyocytes in a cardiac hypertrophy model.PLoS One9:e108051. 10.1371/journal.pone.0108051

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

  AlmughlliqF. B.KohY. Q.PeirisH. N.VaswaniK.HollandO.MeierS.et al (2019). Circulating exosomes may identify biomarkers for cows at risk for metabolic dysfunction.Sci. Rep.9:13879.

  • Google Scholar

• 3

  AndersS.HuberW. (2010). Differential expression analysis for sequence count data.Genome Biol.11:R106.

  • Google Scholar

• 4

  AqeilanR. I.CalinG. A.CroceC. M. (2010). miR-15a and miR-16-1 in cancer: discovery, function and future perspectives.Cell Death Differ.17215–220. 10.1038/cdd.2009.69

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  AravinA.TuschlT. (2005). Identification and characterization of small RNAs involved in RNA silencing.FEBS Lett.5795830–5840. 10.1016/j.febslet.2005.08.009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  BandiN.ZbindenS.GuggerM.ArnoldM.KocherV.HasanL. (2009). miR-15a and miR-16 are implicated in cell cycle regulation in a Rb-dependent manner and are frequently deleted or down-regulated in non-small cell lung cancer.Cancer Res.695553–5559. 10.1158/0008-5472.can-08-4277

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  BarbuM. G.CondratC. E.ThompsonD. C.BugnarO. L.CretoiuD.ToaderO. D.et al (2020). MicroRNA Involvement in Signaling Pathways During Viral Infection.Front. Cell Dev. Biol.8:143. 10.3389/fcell.2020.00143

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  BenmoussaA.LaugierJ.BeauparlantC. J.LambertM.DroitA.ProvostP. (2020). Complexity of the microRNA transcriptome of cow milk and milk-derived extracellular vesicles isolated via differential ultracentrifugation.J. Dairy Sci.10316–29. 10.3168/jds.2019-16880

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  BerezikovE.CuppenE.PlasterkR. H. (2006a). Approaches to microRNA discovery.Nat. Genet.38(Suppl.)S2–S7.

  • Google Scholar

• 10

  BerezikovE.Van TeteringG.VerheulM.Van De BeltJ.Van LaakeL.VosJ.et al (2006b). Many novel mammalian microRNA candidates identified by extensive cloning and RAKE analysis.Genome Res.161289–1298. 10.1101/gr.5159906

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  BillaP. A.FaulconnierY.YeT.ChervetM.Le ProvostF.PiresJ. A. A.et al (2019). Deep rna-seq reveals mirnome differences in mammary tissue of lactating holstein and montbéliarde cows.BMC Genomics20:621. 10.1186/s12864-019-5987-4

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  BionazM.LoorJ. J. (2008). Gene networks driving bovine milk fat synthesis during the lactation cycle.BMC Genomics9:366. 10.1186/1471-2164-9-366

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  BionazM.LoorJ. J. (2011). Gene networks driving bovine mammary protein synthesis during the lactation cycle.Bioinform. Biol. Insights583–98.

  • Google Scholar

• 14

  BionazM.PeriasamyK.Rodriguez-ZasS. L.EvertsR. E.LewinH. A.HurleyW. L.et al (2012). Old and new stories: revelations from functional analysis of the bovine mammary transcriptome during the lactation cycle.PLoS One7:e33268. 10.1371/journal.pone.0033268

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  BonciD.CoppolaV.MusumeciM.AddarioA.GiuffridaR.MemeoL. (2008). The miR-15a-miR-16-1 cluster controls prostate cancer by targeting multiple oncogenic activities.Nat. Med.141271–1277. 10.1038/nm.1880

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  CaiW.LiC.LiuS.ZhouC.YinH.SongJ. (2018). Genome Wide Identification of Novel Long Non-coding RNAs and Their Potential Associations With Milk Proteins in Chinese Holstein Cows.Front. Genet.9:281. 10.3389/fgene.2018.00281

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

  ChenB.LiH.ZengX.YangP.LiuX.ZhaoX. (2012). Roles of microRNA on cancer cell metabolism.J. Transl. Med.10:228. 10.1186/1479-5876-10-228

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  ChenP. Y.ManningaH.SlanchevK.ChienM.RussoJ. J.JuJ.et al (2005). The developmental miRNA profiles of zebrafish as determined by small RNA cloning.Genes Dev.191288–1293. 10.1101/gad.1310605

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  CimminoA.CalinG. A.FabbriM.IorioM. V.FerracinM.ShimizuM.et al (2005). miR-15 and miR-16 induce apoptosis by targeting BCL2.Proc. Natl. Acad. Sci. U S A.10213944–13949. 10.1073/pnas.0506654102

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  CoutinhoL. L.MatukumalliL. K.SonstegardT. S.Van TassellC. P.GasbarreL. C.CapucoA. V.et al (2007). Discovery and profiling of bovine microRNAs from immune-related and embryonic tissues.Physiol. Genomics2935–43. 10.1152/physiolgenomics.00081.2006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  CuiX.HouY.YangS.XieY.ZhangS.ZhangY.et al (2014). Transcriptional profiling of mammary gland in Holstein cows with extremely different milk protein and fat percentage using RNA sequencing.BMC Genomics15:226. 10.1186/1471-2164-15-226

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22

  DevarajS.YunJ. M.Duncan-StaleyC.JialalI. (2009). C-reactive protein induces M-CSF release and macrophage proliferation.J. Leukoc Biol.85262–267. 10.1189/jlb.0808458

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23

  EnrightA. J.JohnB.GaulU.TuschlT.SanderC.MarksD. S. (2003). MicroRNA targets in Drosophila.Genome Biol.5:R1.

  • Google Scholar

• 24

  FangJ.SongX. W.TianJ.ChenH. Y.LiD. F.WangJ. F.et al (2012). Overexpression of microRNA-378 attenuates ischemia-induced apoptosis by inhibiting caspase-3 expression in cardiac myocytes.Apoptosis17410–423. 10.1007/s10495-011-0683-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25

  FinnertyJ. R.WangW. X.HebertS. S.WilfredB. R.MaoG.NelsonP. T. (2010). The miR-15/107 group of microRNA genes: evolutionary biology, cellular functions, and roles in human diseases.J. Mol. Biol.402491–509. 10.1016/j.jmb.2010.07.051

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  FinucaneK. A.McfaddenT. B.BondJ. P.KennellyJ. J.ZhaoF. Q. (2008). Onset of lactation in the bovine mammary gland: gene expression profiling indicates a strong inhibition of gene expression in cell proliferation.Funct. Integr. Genomics8251–264. 10.1007/s10142-008-0074-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  FriedländerM. R.MackowiakS. D.LiN.ChenW.RajewskyN. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades.Nucl. Acids Res.4037–52. 10.1093/nar/gkr688

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  FriedmanR. C.FarhK. K.BurgeC. B.BartelD. P. (2009). Most mammalian mRNAs are conserved targets of microRNAs.Genome Res.1992–105. 10.1101/gr.082701.108

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  GaoD.BingC. (2011). Macrophage-induced expression and release of matrix metalloproteinase 1 and 3 by human preadipocytes is mediated by IL-1β via activation of MAPK signaling.J. Cell Physiol.2262869–2880. 10.1002/jcp.22630

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  GiuratoG.De FilippoM. R.RinaldiA.HashimA.NassaG.RavoM.et al (2013). Imir: An integrated pipeline for high-throughput analysis of small non-coding rna data obtained by smallrna-seq.BMC Bioinformatics14:362. 10.1186/1471-2105-14-362

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  GlazovE. A.CotteeP. A.BarrisW. C.MooreR. J.DalrympleB. P.TizardM. L. (2008). A microRNA catalog of the developing chicken embryo identified by a deep sequencing approach.Genome Res.18957–964. 10.1101/gr.074740.107

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  GuZ.EleswarapuS.JiangH. (2007). Identification and characterization of microRNAs from the bovine adipose tissue and mammary gland.FEBS Lett.581981–988. 10.1016/j.febslet.2007.01.081

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  GuptaV.MarkmannK.PedersenC. N.StougaardJ.AndersenS. U. (2012). Shortran: A pipeline for small rna-seq data analysis.Bioinformatics282698–2700. 10.1093/bioinformatics/bts496

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  HennighausenL.RobinsonG. W. (2005). Information networks in the mammary gland.Nat. Rev. Mol. Cell Biol.6715–725.

  • Google Scholar

• 35

  HosuiA.KimuraA.YamajiD.ZhuB. M.NaR.HennighausenL. (2009). Loss of STAT5 causes liver fibrosis and cancer development through increased TGF-{beta} and STAT3 activation.J. Exp. Med.206819–831. 10.1084/jem.20080003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  HuntzingerE.IzaurraldeE. (2011). Gene silencing by microRNAs: contributions of translational repression and mRNA decay.Nat. Rev. Genet.1299–110. 10.1038/nrg2936

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  JiangH.LinJ. J.TaoJ.FisherP. B. (1997). Suppression of human ribosomal protein L23A expression during cell growth inhibition by interferon-beta.Oncogene14473–480. 10.1038/sj.onc.1200858

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  JinW.DodsonM. V.MooreS. S.BasarabJ. A.GuanL. L. (2010). Characterization of microRNA expression in bovine adipose tissues: a potential regulatory mechanism of subcutaneous adipose tissue development.BMC Mol. Biol.11:29. 10.1186/1471-2199-11-29

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  JinW.GrantJ. R.StothardP.MooreS. S.GuanL. L. (2009). Characterization of bovine miRNAs by sequencing and bioinformatics analysis.BMC Mol. Biol.10:90. 10.1186/1471-2199-10-90

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  JuZ.JiangQ.LiuG.WangX.LuoG.ZhangY.et al (2018). Solexa sequencing and custom microrna chip reveal repertoire of micrornas in mammary gland of bovine suffering from natural infectious mastitis.Anim. Genet.493–18. 10.1111/age.12628

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  KanehisaM.ArakiM.GotoS.HattoriM.HirakawaM.ItohM.et al (2008). KEGG for linking genomes to life and the environment.Nucl. Acids Res.36D480–D484.

  • Google Scholar

• 42

  KimJ. E.ChenJ. (2004). regulation of peroxisome proliferator-activated receptor-gamma activity by mammalian target of rapamycin and amino acids in adipogenesis.Diabetes532748–2756. 10.2337/diabetes.53.11.2748

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  KimS. W.RamasamyK.BouamarH.LinA. P.JiangD.AguiarR. C. (2012). MicroRNAs miR-125a and miR-125b constitutively activate the NF-κB pathway by targeting the tumor necrosis factor alpha-induced protein 3 (TNFAIP3. A20).Proc. Natl. Acad. Sci. U S A.1097865–7870. 10.1073/pnas.1200081109

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  KnezevicI.PatelA.SundaresanN. R.GuptaM. P.SolaroR. J.NagalingamR. S.et al (2012). A novel cardiomyocyte-enriched microRNA, miR-378, targets insulin-like growth factor 1 receptor: implications in postnatal cardiac remodeling and cell survival.J. Biol. Chem.28712913–12926. 10.1074/jbc.m111.331751

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  KuksaP. P.Amlie-WolfA.KatanicŽValladaresO.WangL. S.LeungY. Y. (2018). Spar: Small rna-seq portal for analysis of sequencing experiments.Nucl. Acids Res.46W36–W42.

  • Google Scholar

• 46

  Lagos-QuintanaM.RauhutR.MeyerJ.BorkhardtA.TuschlT. (2003). New microRNAs from mouse and human.Rna9175–179. 10.1261/rna.2146903

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  LandgrafP.RusuM.SheridanR.SewerA.IovinoN.AravinA.et al (2007). A mammalian microRNA expression atlas based on small RNA library sequencing.Cell1291401–1414.

  • Google Scholar

• 48

  LangmeadB.TrapnellC.PopM.SalzbergS. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome.Genome Biol.10:R25.

  • Google Scholar

• 49

  LaplanteM.SabatiniD. M. (2009). An emerging role of mTOR in lipid biosynthesis.Curr. Biol.19R1046–R1052.

  • Google Scholar

• 50

  Le GuillouS.SdassiN.LaubierJ.PassetB.VilotteM.CastilleJ.et al (2012). Overexpression of miR-30b in the developing mouse mammary gland causes a lactation defect and delays involution.PLoS One7:e45727. 10.1371/journal.pone.0045727

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 51

  LeeD. Y.DengZ.WangC. H.YangB. B. (2007). MicroRNA-378 promotes cell survival, tumor growth, and angiogenesis by targeting SuFu and Fus-1 expression.Proc. Natl. Acad. Sci. U S A.10420350–20355. 10.1073/pnas.0706901104

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 52

  LewisB. P.BurgeC. B.BartelD. P. (2005). Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets.Cell12015–20. 10.1016/j.cell.2004.12.035

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 53

  LiH. M.WangC. M.LiQ. Z.GaoX. J. (2012a). MiR-15a decreases bovine mammary epithelial cell viability and lactation and regulates growth hormone receptor expression.Molecules1712037–12048. 10.3390/molecules171012037

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 54

  LiQ.QiaoJ.ZhangZ.ShangX.ChuZ.FuY.et al (2020). Identification and analysis of differentially expressed long non-coding rnas of chinese holstein cattle responses to heat stress.Anim. Biotechnol.319–16. 10.1080/10495398.2018.1521337

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  LiR.ZhangC. L.LiaoX. X.ChenD.WangW. Q.ZhuY. H. (2015). Transcriptome microRNA profiling of bovine mammary glands infected with Staphylococcus aureus.Int. J. Mol. Sci.164997–5013. 10.3390/ijms16034997

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  LiZ.LiuH.JinX.LoL.LiuJ. (2012b). Expression profiles of microRNAs from lactating and non-lactating bovine mammary glands and identification of miRNA related to lactation.BMC Genomics13731. 10.1186/1471-2164-13-731

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 57

  LongJ. E.ChenH. X. (2009). Identification and characteristics of cattle microRNAs by homology searching and small RNA cloning.Biochem. Genet.47329–343. 10.1007/s10528-009-9234-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 58

  LuorengZ. M.WangX. P.MeiC. G.ZanL. S. (2018). Comparison of microrna profiles between bovine mammary glands infected with staphylococcus aureus and escherichia coli.Int. J. Biol. Sci.1487–99. 10.7150/ijbs.22498

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 59

  MaoX.CaiT.OlyarchukJ. G.WeiL. (2005). Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary.Bioinformatics213787–3793. 10.1093/bioinformatics/bti430

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 60

  MartelloG.RosatoA.FerrariF.ManfrinA.CordenonsiM.DupontS.et al (2010). A MicroRNA targeting dicer for metastasis control.Cell1411195–1207. 10.1016/j.cell.2010.05.017

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 61

  MathewsD. H.DisneyM. D.ChildsJ. L.SchroederS. J.ZukerM.TurnerD. H. (2004). Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure.Proc. Natl. Acad. Sci. U S A.1017287–7292. 10.1073/pnas.0401799101

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 62

  MobuchonL.MartheyS.Le GuillouS.LaloeD.Le ProvostF.LerouxC. (2015). Food Deprivation Affects the miRNome in the Lactating Goat Mammary Gland.PLoS One10:e0140111. 10.1371/journal.pone.0140111

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 63

  MuroyaS.OgasawaraH.NoharaK.OeM.OjimaK.HojitoM. (2019). Coordinate alteration of mRNA-microRNA transcriptomes associated with exosomes and fatty acid metabolism in adipose tissue and skeletal muscle in grazing cattle.Asian Australas J. Anim. Sci.331824–1836. 10.5713/ajas.19.0682

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 64

  PasseriniL.AllanS. E.BattagliaM.Di NunzioS.AlstadA. N.LevingsM. K.et al (2008). STAT5-signaling cytokines regulate the expression of FOXP3 in CD4+CD25+ regulatory T cells and CD4+CD25- effector T cells.Int. Immunol.20421–431. 10.1093/intimm/dxn002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 65

  PorstmannT.SantosC. R.GriffithsB.CullyM.WuM.LeeversS. (2008). SREBP activity is regulated by mTORC1 and contributes to Akt-dependent cell growth.Cell Metab.8224–236. 10.1016/j.cmet.2008.07.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 66

  PuJ.LiR.ZhangC.ChenD.LiaoX.ZhuY. (2017). Expression profiles of mirnas from bovine mammary glands in response to streptococcus agalactiae-induced mastitis.J. Dairy Res.84300–308. 10.1017/s0022029917000437

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 67

  RapaportF.KhaninR.LiangY.PirunM.KrekA.ZumboP.et al (2013). Comprehensive evaluation of differential gene expression analysis methods for RNA-seq data.Genome Biol.14:R95.

  • Google Scholar

• 68

  RottiersV.NäärA. M. (2012). MicroRNAs in metabolism and metabolic disorders.Nat. Rev. Mol. Cell Biol.13239–250. 10.1038/nrm3313

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 69

  ShaoY.ZhaoF. Q. (2014). Emerging evidence of the physiological role of hypoxia in mammary development and lactation.J. Anim. Sci. Biotechnol.5:9.

  • Google Scholar

• 70

  ShenB.ZhangL.LianC.LuC.ZhangY.PanQ.et al (2016). Deep Sequencing and Screening of Differentially Expressed MicroRNAs Related to Milk Fat Metabolism in Bovine Primary Mammary Epithelial Cells.Int. J. Mol. Sci.17:200. 10.3390/ijms17020200

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 71

  SilveriL.TillyG.VilotteJ. L.Le ProvostF. (2006). MicroRNA involvement in mammary gland development and breast cancer.Reprod. Nutr. Dev.46549–556. 10.1051/rnd:2006026

  • CrossRef
  • Google Scholar

• 72

  SunY. M.LinK. Y.ChenY. Q. (2013). Diverse functions of miR-125 family in different cell contexts.J. Hematol. Oncol.6:6. 10.1186/1756-8722-6-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 73

  TaganovK. D.BoldinM. P.ChangK. J.BaltimoreD. (2006). NF-kappaB-dependent induction of microRNA miR-146, an inhibitor targeted to signaling proteins of innate immune responses.Proc. Natl. Acad. Sci. U S A.10312481–12486. 10.1073/pnas.0605298103

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 74

  TrapnellC.HendricksonD. G.SauvageauM.GoffL.RinnJ. L.PachterL. (2013). Differential analysis of gene regulation at transcript resolution with rna-seq.Nat. Biotechnol.3146–53. 10.1038/nbt.2450

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 75

  VaapilM.HelczynskaK.VilladsenR.PetersenO. W.JohanssonE.BeckmanS.et al (2012). Hypoxic conditions induce a cancer-like phenotype in human breast epithelial cells.PLoS One7:e46543. 10.1371/journal.pone.0046543

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 76

  WangK. Y.MaJ.ZhangF. X.YuM. J.XueJ. S.ZhaoJ. S. (2014). MicroRNA-378 inhibits cell growth and enhances L-OHP-induced apoptosis in human colorectal cancer.IUBMB Life66645–654. 10.1002/iub.1317

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 77

  WenM.ShenY.ShiS.TangT. (2012). miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments.BMC Bioinformatics13:140. 10.1186/1471-2105-13-140

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 78

  WhelanK. A.ReginatoM. J. (2011). Surviving without oxygen: hypoxia regulation of mammary morphogenesis and anoikis.Cell Cycle102287–2294. 10.4161/cc.10.14.16532

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 79

  WhelanK. A.CaldwellS. A.ShahriariK. S.JacksonS. R.FranchettiL. D.JohannesG. J.et al (2010). Hypoxia suppression of Bim and Bmf blocks anoikis and luminal clearing during mammary morphogenesis.Mol. Biol. Cell213829–3837. 10.1091/mbc.e10-04-0353

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 80

  WickendenJ. A.WatsonC. J. (2010). Key signalling nodes in mammary gland development and cancer. Signalling downstream of PI3 kinase in mammary epithelium: a play in 3 Akts.Breast Cancer Res.12:202.

  • Google Scholar

• 81

  WoolI. G. (1996). Extraribosomal functions of ribosomal proteins.Trends Biochem. Sci.21164–165. 10.1016/s0968-0004(96)20011-8

  • CrossRef
  • Google Scholar

• 82

  XiangM.BirkbakN. J.VafaizadehV.WalkerS. R.YehJ. E.LiuS.et al (2014). STAT3 induction of miR-146b forms a feedback loop to inhibit the NF-κB to IL-6 signaling axis and STAT3-driven cancer phenotypes.Sci. Signal7:ra11. 10.1126/scisignal.2004497

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 83

  YangB.JiaoB.GeW.ZhangX.WangS.ZhaoH.et al (2018). Transcriptome sequencing to detect the potential role of long non-coding rnas in bovine mammary gland during the dry and lactation period.BMC Genomics19:605. 10.1186/s12864-018-4974-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 84

  YoungM. D.WakefieldM. J.SmythG. K.OshlackA. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias.Genome Biol.11:R14.

  • Google Scholar

• 85

  ZhangC.WuH.WangY.ZhuS.LiuJ.FangX.et al (2016). Circular rna of cattle casein genes are highly expressed in bovine mammary gland.J. Dairy Sci.994750–4760. 10.3168/jds.2015-10381

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 86

  ZhouL.ChenJ.LiZ.LiX.HuX.HuangY.et al (2010). Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma.PLoS One5:e15224. 10.1371/journal.pone.0015224

  • Pubmed Abstract
  • CrossRef
  • Google Scholar