Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Deep sequencing shows microRNA involvement in bovine mammary gland adaptation to diets supplemented with linseed oil or safflower oil

BMC Genomics volume 16, Article number: 884 (2015) Cite this article

Abstract

Background

Bovine milk fat composition is responsive to dietary manipulation providing an avenue to modify the content of fatty acids and especially some specific unsaturated fatty acid (USFA) isomers of benefit to human health. MicroRNAs (miRNAs) regulate gene expression but their specific roles in bovine mammary gland lipogenesis are unclear. The objective of this study was to determine the expression pattern of miRNAs following mammary gland adaptation to dietary supplementation with 5 % linseed or safflower oil using next generation RNA-sequencing.

Methods

Twenty-four Canadian Holstein dairy cows (twelve per treatment) in mid lactation were fed a control diet (total mixed ration of corn:grass silages) for 28 days followed by a treatment period (control diet supplemented with 5 % linseed or safflower oil) of 28 days. Milk samples were collected weekly for fat and individual fatty acid determination. RNA from mammary gland biopsies harvested on day-14 (control period) and on days +7 and +28 (treatment period) from six randomly selected cows per treatment was subjected to small RNA sequencing.

Results

Milk fat percentage decreased significantly (P < 0.001) during treatment with the two diets as compared to the control period. The individual saturated fatty acids C4:0, C6:0, C8:0, C14:0 and C16:0 decreased significantly (P < 0.05) while five USFAs (C14:1, C18:1n11t, C20:3n3, C20:5n3 and CLA:t10c12) increased remarkably (P < 0.05) in response to both treatments. Analysis of 361 million sequence reads generated 321 known bovine miRNAs and 176 novel miRNAs. The expression of fourteen and twenty-two miRNAs was affected (P < 0.05) by linseed and safflower oil treatments, respectively. Seven miRNAs including six up-regulated (bta-miR-199c, miR-199a-3p, miR-98, miR-378, miR-148b and miR-21-5p) and one down-regulated (bta-miR-200a) were found to be regulated (P < 0.05) by both treatments, and thus considered core differentially expressed (DE) miRNAs. The gene targets of core DE miRNAs have functions related to gene expression and general cellular metabolism (P < 0.05) and are enriched in four pathways of lipid metabolism (3-phosphoinositide biosynthesis, 3-phosphoinositide degradation, D-myo-inisitol-5-phosphate metabolism and the superpathway of inositol phosphate compounds).

Conclusion

Our results suggest that DE miRNAs in this study might be important regulators of bovine mammary lipogenesis and metabolism. The novel miRNAs identified in this study will further enrich the bovine miRNome repertoire and contribute to understanding mammary gland biology.

Background

Milk fat determines the physical properties and quality of milk and milk products and is also the main energy source in whole milk [1]. Bovine milk fat typically consists of 70 % saturated fatty acids (SFAs) which could exert negative effects when consumed in excess, as well as 25 % monounsaturated fatty acids (MUFAs) and 5 % polyunsaturated fatty acids (PUFAs) both of which have potential positive effects on human health [2, 3].

In ruminants, milk fat has been found to be especially responsive to dietary manipulation [4], providing an avenue to modify the fatty acid (FA) profile of milk fat in favor of health promoting isomers, and also a great model for scientists to study the mechanisms of mammary lipogenesis. Dietary supplementation with plant oils that are rich in unsaturated fatty acids (USFAs) can lead to significant reductions in milk fat yield of up to 50 % [5] as well as increased concentrations of conjugated linoleic acids (CLA) in bovine milk [6]. Therefore, supplementation of cow diets with USFAs from plant oils has been considered as a feasible dietary strategy to improve beneficial FA contents of milk and milk products [7]. Linseed oil (rich in α-linolenic acid, about 57 % of total fat, C18:3n3) is one of the typical USFA-rich plant oils that are frequently used as a dietary supplement to manipulate milk FA composition. Numerous studies have shown that dietary supplementation with linseed oil can reduce milk fat yield through repression of de novo FA synthesis, resulting in increased milk CLA, omega-3 FA levels as well as other USFAs [8–11]. Safflower oil/seed (rich in linoleic acid, about 76 % of total fat, C18:2n6) supplementation can also affect FA composition of ruminant milk, resulting in increased concentrations of CLA in bovine milk [6, 12] and in lean tissue of sheep [13]. Recently, it was found that addition of 3 % safflower oil to the diet can induce a decrease in milk fat yield, short-chain FA yield and C16:0 yield and an increase in milk C18:1n10t and CLA:10t12c yield [14].

Milk fat is derived from de novo synthesized FAs by mammary epithelial cells and preformed FAs taken up from blood circulation [15]. Milk fat synthesis requires the coordinated regulation of enzymatic activities in pathways of metabolite transport, de novo lipogenesis, FA transport, desaturation, and esterification, and milk fat globule formation and secretion [16]. Milk fat synthesis is driven by a coordinated participation of key lipogenic enzymes. For example, stearoyl-CoA desaturase 1 (SCD1) is responsible for MUFA synthesis, while FA desaturases (FADS1, FADS2 and FADS3) are involved in synthesis of long chain FAs [17]. Long term dietary supplementation with a mixture of linseed oil and algae resulted in down-regulation of SCD1, FA synthase (FASN) and the regulatory element binding transcription factor (SREBF1) [18]. Abomasal infusion with linseed oil also affected expression of lipogenic genes in the bovine mammary gland, including decreased abundance of SCD [19]. However, the molecular mechanism by which dietary USFAs alter transcription activities of lipogenic genes and thus affect the pathways of lipogenesis as well as individual FA profiles is not well established.

MicroRNAs (miRNAs) are small regulatory RNA molecules that have been shown to be involved in a wide range of biological pathways by modulating expression of specific mRNAs [20]. The role of miRNAs in lipid metabolism has been recently demonstrated [21] and lipid homeostasis was found to be governed in part by an intricate web of miRNA activity [22]. In bovine mammary gland, the expression of numerous miRNAs was increased during the postpartum or the early lactation periods as compared to the dry period [23]. Comparative transcriptome profiling using high throughput sequencing has also revealed considerable differentially expressed miRNAs between lactating and non-lactating bovine mammary glands [24] thus suggesting a critical role of miRNAs in mammary gland development and lactation. Even the overexpression of one miRNA (miR-30d) could cause dysregulation of lactation and delay mammary gland involution in mouse [25]. Another study found that ten miRNAs in early lactation could be essential in mammary lipid biosynthesis of rats by putatively targeting down-regulated genes expressed in the mammary gland [19]. In dairy goats, miR-27a and miR-103 have been found to function as regulators of milk fat metabolism and lactation cycle [26, 27]. Furthermore, studies have found that miRNA expression profile can be manipulated by dietary fat composition in adipose tissues of lambs [28] and cattle [29]. Thus, we suppose that miRNAs could play critical roles in milk fat synthesis in bovine mammary glands and contribute to the process of milk fat depression as well as increased milk beneficial FAs in response to diets rich in USFAs.

In order to explore the potential role of miRNAs in mammary lipogenesis, miRNA expression profiles of bovine mammary glands from cows fed diets supplemented with 5 % (on dry matter bases) linseed oil (6 cows) or safflower oil (6 cows) were investigated using next generation small RNA sequencing technology. The high throughput sequencing data enabled us to examine the bovine mammary gland miRNome profile and uncovered differentially expressed miRNAs in response to the dietary treatments, which will help further understanding of miRNA involvement in bovine mammary gland lipogenesis.

Methods

Ethics statement

All the experimental procedures were approved by the Animal Care and Ethics Committee of Agriculture and Agri-Food Canada.

Animals, diets and samples collection

Twenty-four high producing Canadian Holstein cows in mid-lactation were blocked according to parity and days in milk and randomly assigned to one of two treatments (12 cows/treatment), namely linseed oil treatment (rich in α-linolenic acid, C18:3n3, 57 % of total fat) and safflower oil treatment (rich in linoleic acid, C18:2n6, 76 % of total fat) (Additional file 1). The experiment started with a control period of 28 days during which time the diet consisted of total mixed ration of corn:grass silages (50:50) and concentrates (control diet). This was followed by a period of supplementation of 28 days (treatment period) during which cows were fed the control diet supplemented with 5 % linseed oil (DM) or 5 % safflower oil. Mammary gland biopsies were collected from 12 cows (6/treatment, randomly selected) at three different time points during the experiment: day-14 (middle of control period), day + 7 (seven days after onset of treatment), and day + 28 (28 days after onset of treatment) following an established protocol [30]. Biopsy was performed on the same quarter and approximately the same site each time. About 400 mg of tissue (70 x 4 mm in diameter) was collected from each cow per biopsy, snap-frozen in liquid nitrogen and stored at −80 °C until isolation of total RNA.

Milk samples were collected from cows weekly (day-14, day-7, day-1, day + 7, day + 14, day + 21 and day + 28). At each sampling, a volume of 40 ml of milk from each cow was collected at morning and evening milking and thoroughly mixed to get a composite sample. One portion of the composite milk was sent to Valacta (St-Anne-de-Bellevue, Quebec, www.valacta.com) for determination of test day milk fat percentage. The other portion was used for determination of individual FA composition. Fatty acid methyl esters (FAME) were prepared from milk fat according to O’Fallon et al. [31]. FA composition of milk fat was determined using Agilent 6890 N gas chromatograph (Agilent Technologies, USA) equipped with flame ionisation detector.

Statistical analysis

The response of test day milk fat percentage and individual FAs to treatments was analyzed with SAS software (version 9.3, SAS Institute Inc., USA). A completely randomized design with repeated measures and mixed effects ANOVA model (F-test) with Tukey adjustment was used to determine the effects of treatments on test day milk fat percentage and FA composition by comparing early treatment period (day + 7) vs. control (day-14), late treatment period (day + 28) vs. control and late treatment period vs. early treatment period. For some FAs, the data was asymmetrically distributed and thus log transformed to fit a normal distribution before applying the ANOVA model. The mean values of the log transformed data were back transformed with a 95 % confidence interval for result interpretation.

RNA extraction

Total RNA from mammary gland biopsy (50 mg/sample) at day-14 (control period), day + 7 and day + 28 (treatment period) were purified using miRNeasy Kit (Qiagen Sciences, USA) following manufacturer’s instructions. The quantity and quality of purified RNA were respectively determined using NanoDrop 1000 (NanoDrop Technologies, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, USA). RNA samples for sequencing had RIN (RNA integrity number) values from 7.5 to 9.0.

miRNA library preparation and sequencing

A total of thirty-six libraries were prepared and barcoded for sequencing according to a previous protocol [32] with minor modifications. Briefly, total RNA was first ligated to a primer at the 3’end (3’adaptor) by T4 RNA Ligase 22tr K227Q (New England Biolabs Inc., Canada) which was then annealed to a reverse transcription primer to prevent undesirable dimerization of 3’ and 5’ adaptors in the following step. Before reverse transcription, the 5’adaptor was ligated to the 5’end of the RNA by T4 RNA Ligase 1 (Enzymatics Inc., USA). This RNA:DNA hybrid was reverse transcribed into cDNA using Superscript III Kit (Life Technologies, USA). The barcoding of the different libraries was finally achieved by PCR using primers which consisted of specific barcodes for each library. Size separation of the miRNA libraries was performed by polyacrylamide gel electrophoresis. The libraries were then eluted from the gel using an elution buffer (10 mM Tris–HCl pH 7.5; 50 mM NaCl, 1 mM EDTA) and concentrated using DNA clean and concentrator-5 from Zymo Research (Zymo Research, USA). Finally, the concentration of the purified libraries was evaluated by Picogreen assay (Life Technologies) on a Nanodrop 3300 fluorescent spectrophotometer.

Multiplexed libraries (12 libraries per lane) were subjected to 50 bp single end sequencing on an Illumina HiSeq 2000 system (Illumina Inc. USA) by McGill University and Genome Quebec Innovation Centre (Montreal, QC, Canada) using TruSeq v3 (Illumina Inc.) reagents. Barcodes (indexes) and adaptor sequences for multiplexed samples are shown in additional file 2.

Small RNA sequence data analysis

The raw data of 36 fastq files were checked for sequencing quality with FastQC program version 0.10.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Cutadapt v1.2.2 (http://code.google.com/p/cutadapt/) was used to trim 3’ adaptor sequences and filter 5’ adaptor contaminants and repeats. Reads shorter than 18 nucleotides after trimming or having a low Phred score of less than 20 for at least 50 % of the bases were discarded. Low quality bases with Phred score <20 corresponding to a probability of error higher than 0.01 [33] as well as reads containing unknown bases were removed using FASTQ Quality Filter tool of FASTX-toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). For novel miRNA discovery, clean reads from the 35 files that passed all filtering criteria were parsed into one file and mapped to the bovine genome (Bta_4.6.1) using bowtie 1.0.0 [34]. Reads that mapped to more than five positions of the genome were discarded. Furthermore, reads that mapped to bovine rRNA, tRNA, snRNA and snoRNA in the Rfam RNA family database (http://rfam.sanger.ac.uk/) were annotated using blastn of blast + (v2.2.28) package [35].

Identification of known miRNA and novel miRNA discovery

The identification of known miRNA and discovery of novel miRNA were performed using miRDeep2 v2.0.0.5 [36] which uses a probabilistic algorithm based on the miRNA biogenesis model and designed to detect miRNAs from deep sequencing reads. The core module of miRDeep2 was used to identify known miRNA and novel miRNA in the pooled dataset of all the libraries. The Quantifier module of miRDeep2 was used to profile the detected miRNAs in each library. The miRDeep2 software outputs a scored list of known and novel miRNAs with log-odds score to help in the determination of true known and novel miRNAs. In this study, only miRNAs with a greater than one count per million (CPM) in at least ten libraries were considered as true known miRNAs. In order to predict novel miRNAs with high confidence, only those with a miRDeep2 score higher than five and present in at least three libraries with at least 10 total counts were retained as true novel miRNAs. The other potential novel miRNAs that didn’t meet all the above criteria but with a miRDeep2 score higher than five or at least 100 reads were classified as novel miRNA candidates. Possible modifications within miRNA sequences were analyzed with Chimira v1 program [37] (http://wwwdev.ebi.ac.uk/enright-srv/chimira/about.php).

Differential expression analysis

The R (v3.0.1) package edgeR (v2.4.6) [38] which uses a negative binomial model was used to identify significantly differentially expressed miRNAs (including known and novel miRNAs) in response to the treatments. The trimmed mean of Mvalues (TMM) normalization method, which assumes that the majority of the genes are not differentially expressed [39], was used for normalization of the miRNA expression count data. Significantly differentially expressed miRNAs were defined as having a Benjamini and Hochberg [40] corrected P value ≤ 0.05.

Function analysis of significantly differentially expressed miRNAs

The target genes of DE miRNAs were predicted using Ingenuity Pathway Analysis (IPA) software (Ingenuity Systems Inc., USA). MiRNA target prediction information in IPA database is very comprehensive as it includes not only bioinformatics predictions using TargetScan (www.targetscan.org), but also from experimentally validated information on gene-miRNA interactions from TarBase database (http://www.microrna.gr/tarbase) and miRecords (http://c1.accurascience.com/miRecords/). The miRNA target filter function of IPA enabled us to focus on the target genes that were experimentally observed or predicted with high confidence. For those miRNAs (bta-miR-199c, miR-3431, miR-2299-5p, miR-152 and miR-1388-5p) without target information in IPA database, we used the perl scripts from the TargetScan website (http://targetscan.org) to predict (targetscan_60.pl) and calculate the context scores (targetscan_61_context_scores.pl) of their gene targets. Predicted targets with context + scores above 95th percentile were further used. Predicted gene targets of DE miRNAs were then subjected to function and pathway analysis using IPA core analysis function. Targets with potential associations with lipid metabolism were used to produce gene networks.

Quantitative PCR and data analysis

Quantitative RT-PCR was used to validate the expression of select differentially expressed miRNAs (bta-miR-199a-3p, miR-378, miR-34a and miR-98). Sequences of the most dominant isomiR of these miRNAs in our study were submitted to Exiqon (Exiqon Inc., USA) for use in designing the miRCURY LNA™ miRNA qPCR assays (product numbers 205968 [bta-miR-199a-3p], 204486 [bta-miR-34a], 205946 [bta-miR-378] and 204640 [bta-miR-98] used. Twelve samples (six/treatment) of day-14 and day + 28 were used for validation, as most of the miRNAs were significantly regulated at day + 28 when compared with day-14. Total RNA was reverse transcribed using Universal cDNA Synthesis Kit II (Exiqon). The cDNA was then quantified with the ExiLENT SYBR® Green Master Mix Kit using the miRCURY LNA™ assay according to manufacturer’s instructions. Bta-miR-103 (product No. 204063) was used as internal control as it was one of the most stable miRNA in our study (C.V = 11.4 %) and was also one of the reference assays suggested by Exiqon.

The expression of miRNA in a sample was normalized to bta-miR-103 and calculated relative to day-14 as follows: ΔCttarget miRNA = Cttarget miRNA–Ctbta-miR-103; relative expression of target miRNA = 2-ΔΔCt, where ΔΔCt = ΔCttarget miRNA, day+28 - ΔCttarget miRNA, day-14. Significant differential expression was declared at P < 0.05. A repeated measures ANOVA model was used:

$$ {\mathrm{Y}}_{\mathrm{ijk}} = \upmu + {\mathrm{cow}}_{\left(\mathrm{i}\right)} + {\mathrm{time}}_{\mathrm{k}} + {\mathrm{e}}_{\mathrm{ik}} $$

where: μ = general mean; cow(i) = random effect of cow j (experimental unit on which three measurements were made); timek = time effect (k: −14, 28); eik = residual error.

Results

Milk fat percentage and individual fatty acid proportions during the period of linseed/safflower oil supplementation

Dietary supplementation of cow diets with 5 % linseed oil reduced (P < 0.001) milk fat percentage from 3.711 % ± 0.183 in the control period (day-14) to 2.411 % ± 0.167 in the late treatment period (day + 28) while diets with 5 % safflower oil caused a reduction (P < 0.001) from 3.617 % ± 0.134 in the control period (day-14) to 2.536 % ± 0.121 in the late treatment period (day + 28) (Fig. 1). Five (C4:0, C6:0, C8:0, C14:0 and C16:0) out of seven saturated FAs measured and one USFA (C18:2n6tt) were significantly decreased by the two treatments (P < 0.05 to 0.0001) (Table 1). The concentration of C17:0 was also decreased by both treatments but significantly by linseed treatment only (P = 0.015) while C18:0 increased remarkable seven days after onset of treatment (significant for safflower treatment, P = 006) and dropped thereafter. The levels of five (C14:1, C18:1n11t, CLA:10t12c, C20:3n3 and C20:5n3) out of 13 USFAs measured increased significantly (P < 0.05 to ) in response to both treatments (Table 1). Significant (P < 0.05) increase in CLA:9c11t was recorded for linseed oil treatment only. The concentrations of C18:1n9c, C18:2n6cc and C18:3n3 were increased by both treatments but not significantly. The proportions of C14:1 t, C18:1n9t, C22:5n3 and C22:5n6 were unchanged by treatments.

Fig. 1
figure 1

Response of test day milk fat percentage to 5 % linseed oil or 5 % safflower oil treatments (***P < 0.001; **P < 0.01, ns P > 0.05)

Full size image
Table 1 Change in milk fatty acid profiles (g/100 g) in response to 5 % linseed oil treatment or 5 % safflower oil treatment
Full size table

High-throughput sequencing of small RNA libraries prepared from bovine mammary gland tissues

Thirty five libraries out of 36 were successfully sequenced, yielding more than 387 million raw reads with an average of 11 million reads per sample (Fig. 2a and Additional file 3). One library (cow 62 at day + 28) with only 15,149 raw reads was discarded and not further processed. Following adaptor removal, quality filtering and removal of short sequences (<18 nt), more than 361 million reads (361,597,513 reads) were retained for further analysis. Out of this number, more than 260 million reads (about 67.2 % of the total reads) were successfully mapped to the bovine genome (Btau_4.6.1) (Fig. 2a). A small proportion of the mapped reads (5.6 %) were further classified into different small RNA categories (rRNA, tRNA, snRNA and snoRNA) (Fig. 2b). The remaining reads (245,563,876) were used to profile known miRNAs and to discover novel miRNAs using miRDeep2 [36]. The length distribution of the mapped reads showed a sharp peak at 22 nt (Fig. 2c).

Fig. 2
figure 2

Raw data filtering and mapping statistics. a The statistics of raw reads from the preprocessing and mapping steps. A total of 387.4 million raw reads were processed. A total of 11.9 million reads were filtered during the preprocessing step (aqua bar). Reads that could not map to the bovine genome were 115.3 million (purple bar). This number (115.3 million) includes reads which mapped to more than 5 genomic locations. Finally, reads that mapped successfully to the bovine genome (Bta_4.6.1) were 260.2 million (red bar). b The relative abundance of different classes of small RNAs in the total reads that were successfully mapped to the bovine genome; (c) Length distribution of the mapped reads across all the libraries

Full size image

Known miRNA identification and genomic expression pattern

Based on miRBase Release 21, we identified a total of 321 miRNAs with an abundance of more than one CPM in at least ten libraries (Additional file 4). MiRNAs were the predominant small RNA species, accounting for 82.2 % of mapped reads (Fig. 2b). Ten miRNAs were highly expressed, accounting for about 70.5 % of all the known miRNAs (Fig. 3) with similar expression levels across all the libraries (Additional file 5), while the remaining miRNAs were expressed at much lower levels. The genomic distribution of the detected miRNAs was diverse among all the chromosomes (Chr). A large number of miRNAs were located on Chr 21 (44 miRNAs), X (41 miRNAs), and 19 (30 miRNAs), while only 11 miRNAs in total were mapped on Chr 6, 17, 27 and 28 (Fig. 4a). Most miRNAs (195) were located within intergenic regions, followed by intronic regions (118) while only a small number (7) were found within exons (Fig. 4b).

Fig. 3
figure 3

Top10 highly expressed miRNAs in bovine mammary gland following supplemental feeding with 5 % linseed oil and 5 % safflower oil

Full size image
Fig. 4
figure 4

Distribution of known miRNAs (blue bar) and novel miRNAs (red bar) across bovine chromosomes (a) and bovine genomic regions (b)

Full size image

Furthermore, we compared the dominant isomiRs of each miRNA in our study with those in miRBase (Release 21). Our results indicate that, the dominant isomiRs of more than half the miRNAs in our study (194 out of 321) differed from the consensus sequences deposited in miRBase (Release 21) (Additional file 6), while only 127 (out of 321) miRNAs share the same consensus sequences as in miRBase (Release 21). To understand the reason for this discrepancy, we conducted global modification analysis [37] of the dominant isomiR sequences in our study and found out that single nucleotide 3’extensions were the most prevalent type of modifications while 5’end modifications were much lower (Additional file 7a). The prevalent 3’extensions were U and A nucleotides. Internal modification analysis showed that deamination of adenosine to inosine bases at positions 4 to 20 by adenosine deaminases was the only type of internal modification present (Additional file 7b). We also examined the arm selection of the precursors of miRNAs in our study using the formula ω = 5p/(5p + 3p) as a measure for arm selection [41]. Out of the 293 precursors encoding the 321 miRNAs identified in this study, 114 were predominantly expressed at the 5’ arm (0.9 ≤ ω ≤ 1) and 107 at the 3’ arm (0 ≤ ω ≤ 0.1) as well as 72 precursors co-expressed at both arms (0.1 < ω < 0.9) (Additional file 8).

Novel miRNA discovery

In this study, miRDeep2 score of 5 which produced a signal-to-noise ratio of 15.7 was used as the cut-off value for novel miRNA prediction [42] (Additional file 9a). Using this cut-off value and miRBase (Release 21), 261 novel miRNAs were initially detected by miRDeep2. In order to identify novel miRNAs with high confidence, only those novel miRNAs with more than 10 total counts and present in at least three libraries were considered true novel miRNAs. Based on our criteria, 187 novel miRNA precursors coding for 176 putatively novel mature miRNAs (Additional file 9b) were identified in this study. Novel miRNAs that did not fit our criteria were considered miRNA candidates (Additional file 9c). The highest number of novel miRNAs were found on Chr X (12 miRNAs), followed by Chr 13 (11 miRNAs) with the least number on Chr 8 (6 miRNAs) (Fig. 4a). Chr 27 was the only chromosome devoid of novel miRNAs. Furthermore, we also identified five novel and two already reported [43] miRNA precursors that were reversely complementary to known bovine miRNA precursors (Additional file 9d). Additionally, two novel mature miRNAs, miR-nv163 and miR-nv185, were reversely complementary to known miRNAs (bta-miR-194 and miR-338, respectively).

Differentially expressed miRNAs in response to linseed oil supplementation

The correlations between libraries (day-14, day + 7 and day+28) of each cow were calculated using the normalized counts of expressed miRNAs. The lowest correlation between different time points of each cow in linseed oil treatment was R2 = 0.956, indicating a very high correlation between all the libraries (Additional file 10a). The expression of 14 miRNAs were significantly altered (P < 0.05) at day + 28 (treatment period) as compared to the control period (day-14) (Table 2). Out of this number, 11 were up-regulated (bta-miR-4286, miR-885, miR-199c, miR-199a-3p, miR-3431, miR-98, miR-196a, miR-378, miR-23b-3p, miR-148b and miR-21-5p) while only 3 were down-regulated (miR-200a, miR-335 and miR-2299-5p) (Table 2). All the significantly DE miRNAs demonstrated a subtle change ranging from 1.28 to 1.73 fold. No significantly DE miRNAs were identified at day+7 (one week after onset of treatment) as compared to Day-14.

Table 2 Differentially expressed miRNAs in response to 5 % linseed oil or 5 % safflower oil treatments
Full size table

Differentially expressed miRNAs in response to safflower oil supplementation

For cows supplemented with safflower oil, there was also a high correlation (R2 ≥ 0.964) between libraries (Additional file 10b). When compared with the control period (day-14), we identified a total of 22 DE miRNAs at day+28 including 10 up-regulated (bta-miR-199c, miR-199a-3p, miR-98, miR-378, miR-21-5p, miR-148b, miR-34a, miR-152, miR-16a, and miR-28) and 12 down-regulated (bta-miR-200a, miR-145, miR-99a-5p, miR-125b, miR-99b, miR-125a, miR-96, miR-484, miR-1388-5p, miR-342, miR-486 and miR-1271) (Table 2). Among the DE miRNAs, bta-miR-1271 was the most affected showing a 2.26 fold down-regulation followed by bta-miR-486 with a 1.75 fold decrease. Except for bta-miR-1271 and miR-486, the majority of the DE miRNAs showed subtle changes in fold gene expression. Additionally, only the expression of bta-miR-486 was significantly affected at day+7 (7 days following onset of treatment) as compared to the control period (day-14), showing a 2.84 fold down-regulation. Seven of the DE miRNAs by safflower oil treatment (6 up-regulated: bta-miR-199c, miR-199a-3p, miR-98, miR-378, miR-148b, miR-21-5p; one down-regulated: bta-miR-200a) were also significantly affected by linseed oil supplementation. However, the effect of safflower oil treatment on miRNA expression in this study was more pronounced (22 DE miRNAs) as compared to linseed oil treatment (14 DE miRNAs) (Table 2).

Real-time quantitative PCR validation of identified miRNA expression profiles

Real-time quantitative PCR (qPCR) was used to validate the expression levels of selected DE miRNAs (bta-miR-199a-3p, miR-34a, miR-378 and miR-98) identified in this study (Fig. 5). Expression trends of selected miRNAs were generally similar to the results from miRNA-sequencing. For example, bta-miR-199a and bta-miR-378 were differentially expressed at day+28 compared with day-14 in both treatments (P < 0.05) after RNA-sequencing and confirmed by qPCR. Bta-miR-34a in linseed treatment was not significantly differentially expressed using qPCR and miRNA-sequencing. The expression level of bta-miR-34a was higher with qPCR detection but only approached significance (P = 0.056) as compared to RNA-sequencing results. The only exception was the qPCR results of bta-miR-98, which showed a reversed trend from RNA-sequencing results in linseed treatment.

Fig. 5
figure 5

Results of qPCR validation of the expression of four miRNAs in linseed (a) and safflower (b) treatments, when comparing day+28 with day-14

Full size image

Target gene prediction of differentially expressed miRNAs and functional annotation

In order to explore possible roles of the DE miRNAs in mammary lipogenesis, we used Ingenuity Pathway Analysis (IPA) software to predict DE miRNA gene targets and to investigate the enriched functions which might be impacted by the alterations of miRNA expression by diets rich in USFAs. Since there is a high false positive rate for miRNA target prediction, we only analyzed the targets with high confidence (highly predicted or experimentally observed).

Our analysis indicates that, about 3594 genes could be targeted by the eleven up-regulated miRNAs (bta-199a-3p, miR-98, miR-378, miR-21-5p, miR-148b, miR-4286, miR-885, miR-196a, miR-23b-3p, bta-miR-199c and miR-3431) whereas 1163 genes could be targeted by the three down-regulated miRNAs (bta-miR-335, miR-200a and bta-miR-2299-5p) in linseed oil-treated cows. Predicted gene targets were further subjected to core analysis to explore their functions. The targets of up-regulated and down-regulated miRNAs were associated (P < 0.05) with 18 and 16 molecular/cellular functions respectively (Additional file 11). The most significantly enriched functions for all the gene targets for both up-regulated and down-regulated DE miRNAs were related to gene expression, followed by functions associated with general cellular metabolism such as cellular growth and proliferation, cellular development, and cellular morphology.

In the safflower oil-treated group, the ten up-regulated (bta-miR-199a-3p, miR-98, miR-378, miR-21-5p, miR-148b, miR-34a, miR-16a, miR-28, miR-199c and miR-152) and nine down-regulated (bta-miR-200a, miR-145, miR-99b, miR-125a, miR-484, miR-342, miR-486, miR-1271 and miR-1388-5p) miRNAs could target respectively 1610 and 1300 genes with high confidence while three miRNAs (bta-miR-99a-5p, miR-125b and miR-96) only have targets with low confidence in IPA knowledge base and were not considered for further analysis. Gene targets for both up- and down-regulated miRNAs were most enriched in functions related to gene expression (P < 0.05), followed by general cellular metabolism functions similar with those of linseed treatment (Additional file 11).

Functional annotation of core differentially expressed miRNAs shared by the two treatments

The expression of seven miRNAs including six up-regulated (bta-miR-199c, miR-199a-3p, miR-98, miR-378, miR-148b and miR-21-5p) and one down-regulated (bta-miR-200a) were significantly affected by both treatments. Since both treatments significantly induced milk fat depression, the common DE miRNAs identified in this study were more likely to have critical roles in the regulation of mammary lipogenesis and thus considered core DE miRNAs. The seven core DE miRNAs could target 2922 genes with high confidence. The highest number of gene targets (1180) was recorded for bta-miR-98 while miR-199a-3p only targeted 81 genes (with high confidence). The most significantly enriched function by target genes of core DE miRNAs was cellular growth and proliferation, followed by gene expression (Fig. 6a). We further explored the pathways associated with lipid metabolism. Among all the affected canonical pathways, four including 3-phosphoinositide biosynthesis, 3-phosphoinositide degradation, D-myo-inisitol-5-phosphate metabolism and the superpathway of inositol phosphate compounds, were related with lipid metabolism (P < 0.05) (Fig. 6b). The IPA filter function enabled us to further identify the target genes that are associated with lipid metabolism. The largest number of target genes (39) were associated with two functions (synthesis of lipid and concentration of lipid) related with lipogenesis (Fig. 7).

Fig. 6
figure 6

Functional annotation and enriched canonical pathways of core differentially expressed miRNA targets. a Functional annotation of all the targets by core differentially expressed miRNAs (b) Enriched canonical pathways related with lipid metabolism

Full size image
Fig. 7
figure 7

Two functions of core miRNAs with the highest number of targets related with lipid metabolism. Red color indicates up-regulated miRNAs; Green color indicates down-regulated miRNAs; Yellow color shows genes targeted by more than one miRNAs

Full size image

Discussion

In this study, we examined the effect of linseed/safflower oil treatments on the miRNome expression in bovine mammary glands using high throughput miRNA sequencing method. Feeding ingredients with high concentrations of USFAs (linseed oil, soybean oil, corn oil, safflower oil, etc.) have been considered as an efficient strategy to alter the milk FA profiles of bovine milk in favor of health promoting components, mainly the increase in USFAs including CLA, α-linolenic acid and linoleic acid [5, 8–10, 44]. These reports are supported by our findings. The increased beneficial USFA contents recorded in this study could be associated with the altered expression of mammary lipogenic enzymes [45, 46]. The concentration of trans-FA and CLA isomers in milk are negatively correlated with the abundance of key lipogenic genes and transcription factors such as peroxisome proliferator-activated receptor gamma (PPARγ) and sterol regulatory element binding factor 1 (SREBP1) [47]. MiRNAs are universal regulators of gene expression at post-translational level and thus are potentially involved in the mammary adaptation process to increased availability of dietary USFAs and modulation of the pathways of lipogenesis. Deep sequencing of six samples per treatment produced an average of more than 10 million reads per sample, enabling us to profile the miRNA expression with high depth and identify miRNAs significantly affected by the diets with high confidence.

We identified 321 miRNAs with an average count greater than one CPM in at least ten libraries, accounting for 40.35 % of all known bovine miRNAs deposited in miRBase (Release 21). Six (bta-miR-148a, miR-26a, miR-21-5p, miR-27b, le-7f and let-7a-5p), four (bta-miR-30a-5p, miR-26a, miR-21-5p and let-7a-5p) and five (bta-miR-148a, miR-26a, let-7a-5p, miR-143 and miR-21-5p) of the highly expressed miRNAs in our study are also among the top 10 highly expressed miRNAs detected respectively in bovine mammary epithelial cells (MAC-T) [43] and lactating glands [24, 48]. Our results and these reports suggest that these highly abundant miRNAs may be implicated in the lactation process, milk synthesis and in mammary gland biology. Previous studies have already unveiled the roles of some of these highly expressed miRNAs in mammary glands. For example, miR-148a and miR-143 are highly expressed in both bovine and goat mammary glands during lactation [49]; miR-148a and miR-26a have been shown to demonstrate consistent expression patterns in bovine milk throughout the lactation period [50] and miR-21-5p increased in expression remarkably at fresh period compared with dry period [23]. Accumulating evidence support the notion that miR-148a can promote cell proliferation and differentiation [51, 52]. Another highly expressed miRNA (miR-21-5p) displayed a marked increase in expression during early lactation indicating that it might be involved in promoting mammary cell proliferation during early lactation [23]. Furthermore, miR-148a can repress WNT (Wingless/INT-1) signaling and thus promote adipogenesis [53] while miR-27b can repress human adipocyte differentiation by directly targeting PPARγ [54]. Repression of miR-143 by its antisense sequence inhibited differentiation of preadipocytes by down regulating ERK5 (also known as MAPK7-mitogen-activated protein kinase 7), indicating a role for this miRNA in adipocyte differentiation [55].

A total of 187 novel miRNA precursors encoding 176 novel miRNAs were detected with high confidence in this study as well as 46 novel miRNAs candidates (had miRDeep2 score higher than 5 but present in less than 3 libraries). Novel miRNAs discovered by high throughput sequencing data have become the main contributors to the miRNome repertoire since 2008 [56]. Considering the number of bovine miRNAs (793 miRNAs) in the miRBase (Release 21), our results will greatly increase the repertoire of bovine miRNAs. Besides, we also found that eight novel precursors were reversely complementary to known bovine miRNA precursors and two novel mature miRNAs were reversely complementary to the known mature miRNAs, which confirmed previous results [43]. This kind of complementation cannot be clearly represented by the current miRNA naming scheme [56–58] which needs to be improved. In this study, we found out that U and A modifications at the 3’end of sequences are the dominant variations, and could explain the differences between the sequences of the dominant isomiRs of 60 % (194 out of 321 miRNAs) of the miRNAs in our study as compared to their miRBase consensus sequences. Previous studies have reported that post-transcriptional alterations modify the 3’ends of miRNAs via mono/poly-uridylation, polydenylation, adenylation as well as other forms of modifications [59–62]. It is also thought that such modifications may modulate miRNA targeting effectiveness [59]. Taken together, the modifications observed in our data may have functional significance which however needs further investigation.

We have investigated the response of miRNA expression of bovine mammary glands to diets supplemented with linseed oil and safflower oil. Linseed oil is rich in α-linolenic acid (57 % of total fat) while safflower oil is rich in linoleic acid (76 % of total fat). Nevertheless, dietary supplementation with the two oil types introduced similar levels of milk fat reductions of ~30 % which is in concordance with previous studies that saw reductions of up to 50 % in milk fat yield [6, 9, 63]. When compared with the control period, only one miRNA was significantly affected at day+7 with safflower supplementation while no miRNAs were significantly affected at day+7 with linseed supplementation. This might be the result of subtle changes at day+7 that could not be detected. Results on milk fat production showed no significant reductions in test day milk fat percentage seven days after onset of treatment with linseed and safflower oils. Subsequently, significant reductions in test day fat percentage were evident 28 days after onset of treatment which coincided with the time when the most miRNAs were significantly differentially regulated by treatments. The seven core DE miRNAs were common to both treatments, making them more likely to be involved in the regulation of lipogenesis in response to USFAs. Functional annotation of their targets showed that these miRNAs were highly enriched in functions associated with cellular metabolism as well as gene expression. Furthermore, four pathways relating to lipid metabolism were also significantly enriched (P < 0.05), implying that these miRNAs could be critical players in these pathways. It should be noted that some of the core DE miRNAs in this study (miR-98, miR-21-5p and miR-200a) have been previously suggested to play important roles in the adipogenesis network [64, 65]. miR-21-5p can negatively regulate PPARα, which is a transcription factor and a major regulator of genes encoding lipid metabolizing enzymes [66]. Elevated expression of miR-21-5p can repress the lipid metabolism pathway regulated by PPARα [67]. Furthermore, miR-200a/miR-8 cluster has been shown to negatively regulate WNT signaling which blocks adipogenic differentiation in multipotent mesenchymal stem cells [68]. Moreover, some important lipogenic enzymes were predicted to be potential targets of the core DE miRNAs in this study. Target analysis showed that stearoyl-CoA desaturases, SCD1 and SCD5, which are involved in FA biosynthesis, are targeted by bta-miR-200a and miR-199a-3p respectively. FADS2, which cause desaturation of FAs is a direct target of bta-miR-98a (IPA knowledge base) while miR-378 can regulate adipocyte differentiation by directly targeting PPARγ and C/EBPα (CCAAT/enhancer-binding protein α), which promote lipogenesis stimulation and increase lipid droplet size in developing adipocytes when overexpressed [69]. Interestingly, six of the seven core DE miRNAs were up-regulated, from which we could expect the repression of a large number of mRNAs related with FA synthesis. This expectation agrees with a previous study which found that feeding conjugated linoleic acid (CLA) repressed milk fat synthesis accompanied by inhibition of the expression of many genes involved in milk lipid synthesis [70]. SCD1 gene, presumably targeted by bta-miR-199a-3p is a key gene with role in the synthesis of USFAs, was reported to be down-regulated in response to linseed oil supplementation [18]. Although we constrained the targets of core DE miRNAs to a small number in this study and demonstrated that they might regulate FA synthesis by targeting lipogenic genes, the involvement of most of the miRNAs in bovine mammary lipogenesis still need to be further validated.

In this study, safflower oil treatment was accompanied by more DE miRNAs (22 DE miRNAs) than linseed oil treatment (14 DE miRNAs) suggesting that there might be treatment-specific effects on the expression of miRNAs, which can further influence expression of lipogenic genes and pathways. Dietary lipids entering the rumen are first transformed to free FAs by microbial lipases, which produce the intermediates for biohydrogenation [71]. Due to the distinctions in chemical structure, the transformation of α-linolenic acid (3 double bonds in their structure) and linoleic acid (2 double bonds) during the biohydrogenation process may require different metabolic pathways [72]. Despite the fact that their biohydrogenation pathways are not fully clarified, it is well acknowledged that they produce diverse intermediates of biohydrogenation [73], which may affect the mammary milk fat synthesis pathway differently. These treatment-specific effects could explain the differential number of significantly regulated miRNAs by the two treatments (Table 2). The functions of most of the DE miRNAs specific to each treatment remain largely unknown. miR-335 is up-regulated in response to lipid loading and is highly abundant in liver and adipose tissues of obese mice, suggesting a role in lipid metabolism [74]. miR-196 has a higher expression in visceral than subcutaneous fat in beef cattle, showing tissue specificity, which could be associated with adipogenesis [29]. miR-342 can inhibit the expression of SREBP, resulting in down-regulation of FASN and 3-hydroxy-3-methylglutaryl CoA reductase (HMGCR) and inhibition of FA biosynthesis in prostate cancer cells [75]. It should also be noted that two of the DE miRNAs in linseed group (bta-miR-2299-5p and miR-3431) are ruminant specific [76, 77]. Bta-miR-2299-5p is one of bovine specific miRNAs associated with insulin signaling pathway [76] whereas involvement of bta-miR-3431 in postnatal muscle development has been demonstrated [77]. These bovine specific miRNAs shape ruminant specific phenotypes [76], indicating that they might also play roles in ruminant specific lipid synthesis pathways, which however need further verification.

Conclusion

In this study, dietary supplementation with 5 % linseed oil or 5 % safflower oil induced significant milk fat reductions, decreased significantly the proportions of some individual SFAs while the proportions of some individual USFAs were significantly increased. Twenty-two miRNAs were significantly differentially expressed in response to the treatments, indicating that miRNAs could be important regulators of bovine mammary lipogenesis. Core differentially regulated miRNAs by both treatments suggest roles in regulating lipogenic pathways or mammary gland functions. The abundant novel miRNAs we identified will further enrich the bovine miRNome repertoire and contribute in the understanding of mammary gland biology.

References

  1. German JB, Dillard CJ. Composition, structure and absorption of milk lipids: a source of energy, fat-soluble nutrients and bioactive molecules. Crit Rev Food Sci Nutr. 2006;46(1):57–92.

    Article  CAS  PubMed  Google Scholar 

  2. MÃ¥nsson HL. Fatty acids in bovine milk fat. Food Nutr Res. 2008;52.

  3. Salter A. Dietary fatty acids and cardiovascular disease. Animal. 2013;7(s1):163–71.

    Article  CAS  PubMed  Google Scholar 

  4. Shingfield KJ, Bonnet M, Scollan ND. Recent developments in altering the fatty acid composition of ruminant-derived foods. Animal. 2013;7(Supplements1):132–62.

    Article  CAS  PubMed  Google Scholar 

  5. Bauman DE, Griinari JM. Regulation and nutritional manipulation of milk fat: low-fat milk syndrome. Livest Prod Sci. 2001;70(1–2):15–29.

    Article  Google Scholar 

  6. Bell JA, Griinari JM, Kennelly JJ. Effect of Safflower Oil, Flaxseed Oil, Monensin, and Vitamin E on Concentration of Conjugated Linoleic Acid in Bovine Milk Fat. J Dairy Sci. 2006;89(2):733–48.

    Article  CAS  PubMed  Google Scholar 

  7. Kouba M, Mourot J. A review of nutritional effects on fat composition of animal products with special emphasis on n-3 polyunsaturated fatty acids. Biochimie. 2011;93(1):13–7.

    Article  CAS  PubMed  Google Scholar 

  8. Bauman DE, Griinari JM. Nutritional regulation of milk fat synthesis. Annu Rev Nutr. 2003;23(1):203–27.

    Article  CAS  PubMed  Google Scholar 

  9. Flowers G, Ibrahim SA, AbuGhazaleh AA. Milk Fatty Acid Composition of Grazing Dairy Cows When Supplemented with Linseed Oil. J Dairy Sci. 2008;91(2):722–30.

    Article  CAS  PubMed  Google Scholar 

  10. Benchaar C, Romero-Pérez G, Chouinard P, Hassanat F, Eugene M, Petit H, et al. Supplementation of increasing amounts of linseed oil to dairy cows fed total mixed rations: Effects on digestion, ruminal fermentation characteristics, protozoal populations, and milk fatty acid composition. J Dairy Sci. 2012;95(8):4578–90.

    Article  CAS  PubMed  Google Scholar 

  11. Lerch S, Ferlay A, Shingfield KJ, Martin B, Pomiès D, Chilliard Y. Rapeseed or linseed supplements in grass-based diets: Effects on milk fatty acid composition of Holstein cows over two consecutive lactations. J Dairy Sci. 2012;95(9):5221–41.

    Article  CAS  PubMed  Google Scholar 

  12. Alizadeh AR, Alikhani M, Ghorbani GR, Rahmani HR, Rashidi L, Loor JJ. Effects of feeding roasted safflower seeds (variety IL-111) and fish oil on dry matter intake, performance and milk fatty acid profiles in dairy cattle. J Anim Physiol Anim Nutr. 2012;96(3):466–73.

    Article  CAS  Google Scholar 

  13. Boles JA, Kott RW, Hatfield PG, Bergman JW, Flynn CR. Supplemental safflower oil affects the fatty acid profile, including conjugated linoleic acid, of lamb. J Anim Sci. 2005;83(9):2175–81.

    CAS  PubMed  Google Scholar 

  14. Stoffel CM, Crump PM, Armentano LE. Effect of dietary fatty acid supplements, varying in fatty acid composition, on milk fat secretion in dairy cattle fed diets supplemented to less than 3 % total fatty acids. J Dairy Sci. 2015;98(1):431–42.

    Article  CAS  PubMed  Google Scholar 

  15. Barber MC, Clegg RA, Travers MT, Vernon RG. Lipid metabolism in the lactating mammary gland. Biochim Biophys Acta Lipids Lipid Metabol. 1997;1347(2):101–26.

    Article  CAS  Google Scholar 

  16. Clegg R, Barber M, Pooley L, Ernens I, Larondelle Y, Travers M. Milk fat synthesis and secretion: molecular and cellular aspects. Livest Prod Sci. 2001;70(1):3–14.

    Article  Google Scholar 

  17. Bionaz M, Loor J. Gene networks driving bovine milk fat synthesis during the lactation cycle. BMC Genomics. 2008;9(1):366.

    Article  PubMed Central  PubMed  Google Scholar 

  18. Angulo J, Mahecha L, Nuernberg K, Nuernberg G, Dannenberger D, Olivera M, et al. Effects of polyunsaturated fatty acids from plant oils and algae on milk fat yield and composition are associated with mammary lipogenic and SREBF1 gene expression. Animal. 2012;6(12):1961–72.

    Article  CAS  PubMed  Google Scholar 

  19. Palin M-F, Côrtes C, Benchaar C, Lacasse P, Petit HV. mRNA Expression of lipogenic enzymes in mammary tissue and fatty acid profile in milk of dairy cows fed flax hulls and infused with flax oil in the abomasum. Br J Nutr. 2014;111(06):1011–20.

    Article  CAS  PubMed  Google Scholar 

  20. Bartel DP. MicroRNAs: Target Recognition and Regulatory Functions. Cell. 2009;136(2):215–33.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  21. Fernández-Hernando C, Suárez Y, Rayner KJ, Moore KJ. MicroRNAs in lipid metabolism. Curr Opin Lipidol. 2011;22(2):86.

    Article  PubMed Central  PubMed  Google Scholar 

  22. Sethupathy P. Needles in the genetic haystack of lipid disorders: single nucleotide polymorphisms in the microRNA regulome. J Lipid Res. 2013;54(5):1168–73.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  23. Wang M, Moisá S, Khan MJ, Wang J, Bu D, Loor JJ. MicroRNA expression patterns in the bovine mammary gland are affected by stage of lactation. J Dairy Sci. 2012;95(11):6529–35.

    Article  CAS  PubMed  Google Scholar 

  24. 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(1):731.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  25. Le Guillou S, Sdassi N, Laubier J, Passet B, Vilotte M, Castille J, 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 Central  PubMed  Google Scholar 

  26. Lin X, Luo J, Zhang L, Wang W, Gou D. MiR-103 Controls Milk Fat Accumulation in Goat (Capra hircus) Mammary Gland during Lactation. PLoS One. 2013;8(11), e79258.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  27. Lin X-z, Luo J, Zhang L-P, Wang W, Shi H-B, Zhu J-J. miR-27a suppresses triglyceride accumulation and affects gene mRNA expression associated with fat metabolism in dairy goat mammary gland epithelial cells. Gene. 2013;521(1):15–23.

    Article  CAS  PubMed  Google Scholar 

  28. Meale SJ, Romao JM, He ML, Chaves AV, McAllister TA, Guan LL. Effect of diet on microRNA expression in ovine subcutaneous and visceral adipose tissues. J Anim Sci. 2014;92(8):3328–37.

    Article  CAS  PubMed  Google Scholar 

  29. Romao JM, Jin W, He M, McAllister T, Guan LL. Altered MicroRNA Expression in Bovine Subcutaneous and Visceral Adipose Tissues from Cattle under Different Diet. PLoS One. 2012;7(7), e40605.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  30. Farr V, Stelwagen K, Cate L, Molenaar A, McFadden T, Davis S. An improved method for the routine biopsy of bovine mammary tissue. J Dairy Sci. 1996;79(4):543–9.

    Article  CAS  PubMed  Google Scholar 

  31. O’fallon J, Busboom J, Nelson M, Gaskins C. A direct method for fatty acid methyl ester synthesis: application to wet meat tissues, oils, and feedstuffs. J Anim Sci. 2007;85(6):1511–21.

    Article  PubMed  Google Scholar 

  32. Vigneault F, Ter-Ovanesyan D, Alon S, Eminaga SC, Christodoulou D, Seidman JG, et al. High-Throughput Multiplex Sequencing of miRNA. In: Current Protocols in Human Genetics. vol. 73 (11.12). New Jersey, USA: John Wiley & Sons, Inc.; 2012. p. 1-11

  33. Ewing B, Hillier L, Wendl MC, Green P. Base-Calling of Automated Sequencer Traces UsingPhred. I. Accuracy Assessment. Genome Res. 1998;8(3):175–85.

    Article  CAS  PubMed  Google Scholar 

  34. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25.

    Article  PubMed Central  PubMed  Google Scholar 

  35. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10(1):421.

    Article  PubMed Central  PubMed  Google Scholar 

  36. Friedländer MR, Chen W, Adamidi C, Maaskola J, Einspanier R, Knespel S, et al. Discovering microRNAs from deep sequencing data using miRDeep. Nat Biotechnol. 2008;26(4):407–15.

    Article  PubMed  Google Scholar 

  37. Vitsios DM, Enright AJ. Chimira: analysis of small RNA sequencing data and microRNA modifications. Bioinformatics 2015. doi:10.1093/bioinformatics/btv380

  38. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  39. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3):R25.

    Article  PubMed Central  PubMed  Google Scholar 

  40. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc B. 1995;289–300. doi:10.2307/2346101

  41. Hu W, Wang T, Yue E, Zheng S, Xu J-H. Flexible microRNA arm selection in rice. Biochem Biophys Res Commun. 2014;447(3):526–30.

    Article  CAS  PubMed  Google Scholar 

  42. Friedländer 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  PubMed Central  PubMed  Google Scholar 

  43. Jin W, Ibeagha-Awemu EM, Liang G, Beaudoin F, Zhao X, Guan L. Transcriptome microRNA profiling of bovine mammary epithelial cells challenged with Escherichia coli or Staphylococcus aureus bacteria reveals pathogen directed microRNA expression profiles. BMC Genomics. 2014;15(1):181.

    Article  PubMed Central  PubMed  Google Scholar 

  44. Luna P, Bach A, Juárez M, de la Fuente MA. Effect of a Diet Enriched in Whole Linseed and Sunflower Oil on Goat Milk Fatty Acid Composition and Conjugated Linoleic Acid Isomer Profile. J Dairy Sci. 2008;91(1):20–8.

    Article  CAS  PubMed  Google Scholar 

  45. Mach N, Jacobs AAA, Kruijt L, van Baal J, Smits MA. Alteration of gene expression in mammary gland tissue of dairy cows in response to dietary unsaturated fatty acids. Animal. 2011;5(08):1217–30.

    Article  CAS  PubMed  Google Scholar 

  46. Li XZ, Yan CG, Lee HG, Choi CW, Song MK. Influence of dietary plant oils on mammary lipogenic enzymes and the conjugated linoleic acid content of plasma and milk fat of lactating goats. Anim Feed Sci Technol. 2012;174(1–2):26–35.

    Article  CAS  Google Scholar 

  47. Mach N, Goselink RMA, van Baal J, Kruijt L, van Vuuren AM, Smits MA. Relationship between milk fatty acid composition and the expression of lipogenic genes in the mammary gland of dairy cows. Livest Sci. 2013;151(1):92–6.

    Article  Google Scholar 

  48. Le Guillou S, Marthey S, Laloë D, Laubier J, Mobuchon L, Leroux C, et al. Characterisation and Comparison of Lactating Mouse and Bovine Mammary Gland miRNomes. PLoS One. 2014;9(3), e91938.

    Article  PubMed Central  PubMed  Google Scholar 

  49. Ji Z, Wang G, Xie Z, Wang J, Zhang C, Dong F, et al. Identification of Novel and Differentially Expressed MicroRNAs of Dairy Goat Mammary Gland Tissues Using Solexa Sequencing and Bioinformatics. PLoS One. 2012;7(11), e49463.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  50. Chen X, Gao C, Li H, Huang L, Sun Q, Dong Y, et al. Identification and characterization of microRNAs in raw milk during different periods of lactation, commercial fluid, and powdered milk products. Cell Res. 2010;20(10):1128–37.

    Article  CAS  PubMed  Google Scholar 

  51. Guo S-L, Peng Z, Yang X, Fan K-J, Ye H, Li Z-H, et al. miR-148a promoted cell proliferation by targeting p27 in gastric cancer cells. Int J Biol Sci. 2011;7(5):567–74.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  52. Zhang J, Ying Z-Z, Tang Z-l, Long L-Q, Li K. MicroRNA-148a promotes myogenic differentiation by targeting the ROCK1 gene. J Biol Chem. 2012;287(25):21093–101.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  53. Qin L, Chen Y, Niu Y, Chen W, Wang Q, Xiao S, et al. A deep investigation into the adipogenesis mechanism: Profile of microRNAs regulating adipogenesis by modulating the canonical Wnt/beta-catenin signaling pathway. BMC Genomics. 2010;11(1):320.

    Article  PubMed Central  PubMed  Google Scholar 

  54. Karbiener M, Fischer C, Nowitsch S, Opriessnig P, Papak C, Ailhaud G, et al. microRNA miR-27b impairs human adipocyte differentiation and targets PPARγ. Biochem Biophys Res Commun. 2009;390(2):247–51.

    Article  CAS  PubMed  Google Scholar 

  55. Esau C, Kang X, Peralta E, Hanson E, Marcusson EG, Ravichandran LV, et al. MicroRNA-143 Regulates Adipocyte Differentiation. J Biol Chem. 2004;279(50):52361–5.

    Article  CAS  PubMed  Google Scholar 

  56. Kozomara A, Griffiths-Jones S. miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res. 2011;39 suppl 1:D152–7.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  57. Griffiths-Jones S, Grocock RJ, van Dongen S, Bateman A, Enright AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34 suppl 1:D140–4.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  58. Ambros V, Bartel B, Bartel DP, Burge CB, Carrington JC, Chen X, et al. A uniform system for microRNA annotation. RNA. 2003;9(3):277–9.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  59. Burroughs AM, Ando Y, de Hoon MJL, Tomaru Y, Nishibu T, Ukekawa R, et al. A comprehensive survey of 3′ animal miRNA modification events and a possible role for 3′ adenylation in modulating miRNA targeting effectiveness. Genome Res. 2010;20(10):1398–410.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  60. Heo I, Joo C, Kim Y-K, Ha M, Yoon M-J, Cho J, et al. TUT4 in Concert with Lin28 Suppresses MicroRNA Biogenesis through Pre-MicroRNA Uridylation. Cell. 2009;138(4):696–708.

    Article  CAS  PubMed  Google Scholar 

  61. Heo I, Ha M, Lim J, Yoon M-J, Park J-E, Kwon SC, et al. Mono-Uridylation of Pre-MicroRNA as a Key Step in the Biogenesis of Group II let-7 MicroRNAs. Cell. 2012;151(3):521–32.

    Article  CAS  PubMed  Google Scholar 

  62. Li S-C, Tsai K-W, Pan H-W, Jeng Y-M, Ho M-R, Li W-H. MicroRNA 3' end nucleotide modification patterns and arm selection preference in liver tissues. BMC Syst Biol. 2012;6 Suppl 2:S14.

    Article  Google Scholar 

  63. Loor JJ, Ferlay A, Ollier A, Doreau M, Chilliard Y. Relationship Among Trans and Conjugated Fatty Acids and Bovine Milk Fat Yield Due to Dietary Concentrate and Linseed Oil. J Dairy Sci. 2005;88(2):726–40.

    Article  CAS  PubMed  Google Scholar 

  64. Chartoumpekis DV, Zaravinos A, Ziros PG, Iskrenova RP, Psyrogiannis AI, Kyriazopoulou VE, et al. Differential expression of microRNAs in adipose tissue after long-term high-fat diet-induced obesity in mice. PLoS One. 2012;7(4), e34872.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  65. Romao J, Jin W, He M, McAllister T, Guan L. MicroRNAs in bovine adipogenesis: genomic context, expression and function. BMC Genomics. 2014;15(1):137.

    Article  PubMed Central  PubMed  Google Scholar 

  66. Kida K, Nakajima M, Mohri T, Oda Y, Takagi S, Fukami T, et al. PPARα is regulated by miR-21 and miR-27b in human liver. Pharm Res. 2011;28(10):2467–76.

    Article  CAS  PubMed  Google Scholar 

  67. Kim YJ, Hwang SJ, Bae YC, Jung JS. MiR-21 Regulates Adipogenic Differentiation through the Modulation of TGF-β Signaling in Mesenchymal Stem Cells Derived from Human Adipose Tissue. Stem Cell. 2009;27(12):3093–102.

    CAS  Google Scholar 

  68. Kennell JA, Gerin I, MacDougald OA, Cadigan KM. The microRNA miR-8 is a conserved negative regulator of Wnt signaling. Proc Natl Acad Sci. 2008;105(40):15417–22.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  69. Gerin I, Bommer GT, McCoin CS, Sousa KM, Krishnan V, MacDougald OA. Roles for miRNA-378/378* in adipocyte gene expression and lipogenesis. Am J Physiol Endocrinol Metabol. 2010;299(2):E198–206.

    CAS  Google Scholar 

  70. Baumgard LH, Matitashvili E, Corl BA, Dwyer DA, Bauman DE. trans-10, cis-12 Conjugated Linoleic Acid Decreases Lipogenic Rates and Expression of Genes Involved in Milk Lipid Synthesis in Dairy Cows1. J Dairy Sci. 2002;85(9):2155–63.

    Article  CAS  PubMed  Google Scholar 

  71. Bauman DE, Harvatine KJ, Lock AL. Nutrigenomics, Rumen-Derived Bioactive Fatty Acids, and the Regulation of Milk Fat Synthesis. Annu Rev Nutr. 2011;31(1):299–319.

    Article  CAS  PubMed  Google Scholar 

  72. Kellens M, Goderis H, Tobback P. Biohydrogenation of unsaturated fatty acids by a mixed culture of rumen microorganisms. Biotechnol Bioeng. 1986;28(8):1268–76.

    Article  CAS  PubMed  Google Scholar 

  73. Jenkins T, Wallace R, Moate P, Mosley E. Board-invited review: Recent advances in biohydrogenation of unsaturated fatty acids within the rumen microbial ecosystem. J Anim Sci. 2008;86(2):397–412.

    Article  CAS  PubMed  Google Scholar 

  74. Nakanishi N, Nakagawa Y, Tokushige N, Aoki N, Matsuzaka T, Ishii K, et al. The up-regulation of microRNA-335 is associated with lipid metabolism in liver and white adipose tissue of genetically obese mice. Biochem Biophys Res Commun. 2009;385(4):492–6.

    Article  CAS  PubMed  Google Scholar 

  75. Li X, Chen Y-T, Josson S, Mukhopadhyay NK, Kim J, Freeman MR, et al. MicroRNA-185 and 342 Inhibit Tumorigenicity and Induce Apoptosis through Blockade of the SREBP Metabolic Pathway in Prostate Cancer Cells. PLoS One. 2013;8(8), e70987.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  76. Bao H, Kommadath A, Sun X, Meng Y, Arantes AS, Plastow GS, et al. Expansion of ruminant-specific microRNAs shapes target gene expression divergence between ruminant and non-ruminant species. BMC Genomics. 2013;14:609.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  77. Sun J, Sonstegard TS, Li C, Huang Y, Li Z, Lan X, et al. Altered microRNA expression in bovine skeletal muscle with age. Anim Genet. 2015;46(3):227–38.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgement

Authors thank Steve Méthot for statistical analysis of data and farm staff of Agriculture and Agri-Food Canada’s Dairy and Swine Research and Development Center for caring for animals during the experimental period. Ran Li was supported by a doctoral fellowship from the China Scholarship Council (201206300128). This work was supported by funding from Agriculture and Agri-Food Canada.

Author information

Authors and Affiliations

  1. Agriculture and Agri-Food Canada, Dairy and Swine Research and Development Centre, Sherbrooke, Quebec, J1M 0C8, Canada

    Ran Li, Frédéric Beaudoin, Adolf A. Ammah, Nathalie Bissonnette, Chaouki Benchaar & Eveline M. Ibeagha-Awemu

  2. College of Animal Science and Technology, Northwest A&F University, Xi’an, Shaanxi, 712100, China

    Ran Li & Chuzhao Lei

  3. Department of Animal Science, McGill University, Ste-Anne-de-Bellevue, Quebec, H9X 3V9, Canada

    Xin Zhao

Authors
  1. Ran Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Frédéric Beaudoin
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Adolf A. Ammah
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Nathalie Bissonnette
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Chaouki Benchaar
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Xin Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Chuzhao Lei
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Eveline M. Ibeagha-Awemu
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Eveline M. Ibeagha-Awemu.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

EIA designed the study with inputs from CB, NB and XZ. FB and AAA participated in the animal phase of the experiment. FB constructed miRNA sequencing libraries. AAA analyzed FA profiles by gas chromatography. RL performed bioinformatics analysis of data, qPCR verifications and drafted the manuscript. All authors contributed in writing manuscript and also approved the final manuscript.

Additional files

Additional file 1:

Composition of the experimental diets. (DOCX 16 kb)

Additional file 2:

Library barcodes of miRNA-seq library preparation. (DOCX 16 kb)

Additional file 3:

Mapping statistics of reads by library. (XLSX 16 kb)

Additional file 4:

Known miRNAs and number of reads per library per sampling point. (XLSX 78 kb)

Additional file 5:

Top expressed miRNAs at different time points according to treatments. (DOCX 17 kb)

Additional file 6:

Summary of most abundant isomiRs of detected miRNAs. (XLSX 36 kb)

Additional file 7:

Summary of modification analysis of dominant isomiRs that differed from consensus miRBase sequences. (A) 5’end (5p) and 3’ end (3p) modifications of isomiRs showed that single nucleotide 3’extensions are dominated by U and A modifications while 5’end modifications are much lower. (B) Internal modifications are due to deamination of adenosine to inosine by adenosine deaminases (ADAR). The start of a miRNA on the x-axis is at index ‘0’ (5' end) while its end is at index ‘22’ (3' end). Note: U, A, C and G refer to additions; G_adar refers to deamination of an adenosine to inosine by adenosine deaminases ; U_snp, A_snp, C_Snp and G_snp refer to snps which are different from the genome. (PNG 212 kb)

Additional file 8:

Summary of dominant arm selection of detected miRNAs. (XLSX 34 kb)

Additional file 9:

mirDeep2 results of known miRNA identification and novel miRNA discovery; ( A ) miRDeep2 performance for score cut-offs −10 to 10, ( B ) novel miRNA information, ( C ) candidate novel miRNAs and ( D ) novel Precursors which were reversely complementary to known precursors. (XLSX 76 kb)

Additional file 10:

Pearson correlation of the libraries in ( A ) safflower and ( B ) linseed oil treatments. (XLSX 58 kb)

Additional file 11:

Functional annotation of the targets of ( A ) up-regulated and ( B ) down regulated miRNAs in linseed oil treatment, and ( C ) up-regulated and ( D ) down-regulated miRNAs in safflower oil treatment. (PDF 290 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), 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 (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Li, R., Beaudoin, F., Ammah, A.A. et al. Deep sequencing shows microRNA involvement in bovine mammary gland adaptation to diets supplemented with linseed oil or safflower oil. BMC Genomics 16, 884 (2015). https://doi.org/10.1186/s12864-015-1965-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-015-1965-7

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us