# Nutrigenomic analyses reveal miRNAs and mRNAs affected by feed restriction in the mammary gland of midlactation dairy cows

Pierre-Alexis Billa, Yannick Faulconnier, Tao Ye, Céline Bourdon, José Pires, Christine Leroux

► To cite this version:

Pierre-Alexis Billa, Yannick Faulconnier, Tao Ye, Céline Bourdon, José Pires, et al.. Nutrigenomic analyses reveal miRNAs and mRNAs affected by feed restriction in the mammary gland of midlactation dairy cows. PLoS ONE, 2021, 16 (4), 10.1371/journal.pone.0248680 . hal-03200342

HAL Id: hal-03200342

https://hal.inrae.fr/hal-03200342

Submitted on 16 Apr 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
RESEARCH ARTICLE

# Nutrigenomic analyses reveal miRNAs and mRNAs affected by feed restriction in the mammary gland of midlactation dairy cows

Pierre-Alexis Billa¹, Yannick Faulconnier¹, Tao Ye², Céline Bourdon³, José A. A. Pires¹,
Christine Leroux¹*

1 INRAE, Université Clermont Auvergne, VetAgro Sup, UMR Herbivores, Saint-Genès-Chamanelle, France, 2 Institut de Génétique et de Biologie Moléculaire et Cellulaire (IGBMC), Centre National de la Recherche Scientifique, UMR7104, Institut National de la Santé et de la Recherche Médicale, U964, Université de Strasbourg, Illkirch, France, 3 INRAE, AgroParisTech, Université Paris-Saclay, UMR Génétique Animale et Biologie Intégrative, Jouy-en-Josas, France

* christine.leroux@inrae.fr

<figure><img src="image_2.png" /></figure>

## Abstract

The objective of this study was to investigate the effects of feed restriction on mammary miRNAs and coding gene expression in midlactation cows. Five Holstein cows and 6 Montbéliarde cows underwent 6 days of feed restriction, during which feed allowance was reduced to meet 50% of their net energy for lactation requirements. Mammary biopsies were performed before and at the end of the restriction period. Mammary miRNA and mRNA analyses were performed using high-throughput sequencing and microarray analyses, respectively. Feed restriction induced a negative energy balance and decreased milk production and fat and protein yields in both breeds. Feed restriction modified the expression of 27 miRNAs and 374 mRNAs in mammary glands from Holstein cows, whereas no significant miRNA change was observed in Montbéliarde cows. Among the 27 differentially expressed miRNAs, 8 miRNAs were associated with dairy QTL. Analysis of target genes indicate that the 8 most abundantly expressed miRNAs control transcripts related to lipid metabolism, mammary remodeling and stress response. A comparison between the mRNAs targeted by the 8 most strongly expressed miRNAs and 374 differentially expressed mRNAs identified 59 mRNAs in common. The bioinformatic analyses of these 59 mRNAs revealed their implication in lipid metabolism and endothelial cell proliferation. These effects of feed restriction on mammary miRNAs and mRNAs observed in Holstein cows suggest a potential role of miRNAs in mammary structure and lipid biosynthesis that could explain changes in milk production and composition.

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Billa P-A, Faulconnier Y, Ye T, Bourdon C, Pires JAA, Leroux C (2021) Nutrigenomic analyses reveal miRNAs and mRNAs affected by feed restriction in the mammary gland of midlactation dairy cows. PLoS ONE 16(4): e0248680. https://doi.org/10.1371/journal.pone.0248680

**Editor:** Angel Abuelo, Michigan State University, UNITED STATES

**Received:** June 5, 2020

**Accepted:** March 3, 2021

**Published:** April 15, 2021

**Copyright:** © 2021 Billa et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All RNAseq files are available from the GEO database (accession number(s) GSE146762). Microarray data were deposited in the Gene Expression Omnibus (GEO): GSE159302.

**Funding:** We thank the Région Auvergne S3 and FEDER for the funding of P-A Billa's thesis grants (project nr. 23000794/1130 entitled "Innovative biomarkers for the quality products of tomorrow's cattle farming ") awarded to PAB, and ApisGene for the funding of Céline Bourdon's thesis grants. This

## Introduction

The mammary gland (MG) is a complex secretory organ with the ability to synthesize and secrete large quantities of milk components, including proteins, lactose and lipids. Genetics, nutritional and physiological studies of MG functioning have shown that many genes regulate
research was funded by GISA meta-program of INRAE (Ruminflame and Longhealth projects).

**Competing interests:** NO authors have competing interests.

the synthesis and secretion of milk components [1]. Forty-eight hours of feed deprivation induced a severe negative energy balance (NEB) and altered the mammary transcriptome, milk production and composition in lactating goats [2]. However, the molecular mechanisms underlying the nutritional regulation of gene expression have not been fully determined.

MicroRNAs (miRNAs), small noncoding RNAs of 18–25 nucleotides, regulate at least 60% of gene expression by base-pairing with mRNAs, inducing their degradation or the inhibition of their translation [3]. Therefore, the RNAs are involved in many different cellular processes [4–6]. MiRNAs were studied in the MG of different species, including ruminants [7–12], and their potential role in mammary development [8] and milk synthesis [13] was demonstrated. Breed effects on miRNA expression were also studied in different tissues, including MG [14–16]. In a previous study, the comparison of mammary miRNomes from Holstein (HOLS) and Montbéliarde (MONT) cows revealed 22 differentially expressed miRNAs (DEMs), suggesting a potential role of miRNAs in MG plasticity and milk component synthesis [15].

In ruminants, few studies have investigated the effect of nutrition on miRNA expression in the MG. Feed deprivation leading to severe NEB modified the expression of 30 miRNAs in MG of lactating goats [10].

We hypothesized that bovine mammary miRNome would be modified by feed restriction. Thus, the objective was to determine the effects of NEB induced by feed restriction on MG miRNome, assess potential differences between midlactation HOLS and MONT cows, and determine the role of miRNAs in the regulation of milk production and milk component synthesis.

# Materials and methods

## Ethics statement

The experiment was conducted at the INRAE Herbipole experimental farm of Marcenat (https://doi.org/10.15454/1.5572318050509348E12). Animal procedures were performed in compliance with Regional Animal Care Committee guidelines CEMEAA: Auvergne, French Ministry of Agriculture and European Union guidelines for animal research C2EA-02. All mammary biopsies were performed with relevant national legislation and were performed by an accredited person (N° of certification: 63-20). All procedures involving animals were approved by the Ethics Committee of the French Ministry of Higher Education, Research and Innovation (APAFIS # 3737-2015043014541577v2).

## Animals and sampling

Twelve multiparous midlactation HOLS (n = 6; 2.5 ± 0.8 lactations and 170 ± 11 days in milk (DIM)) and MONT (n = 6; 3.3±2.0 lactations and 156 ± 26 DIM) cows were used to study the effects of 6 days of feed restriction on milk production and mammary miRNomes. One HOLS cow was excluded from the experiment due to mastitis. During the control period (CONT) before restriction, all cows were allowed *ad libitum* intake of a total mixed ration. During the restriction period (d1 to 6, REST), feed allowance was reduced to meet 50% of NE<sub>L</sub> requirements calculated from body weight, feed intake and milk production and composition recorded before restriction [17]. Cows had free access to water and were housed in a freestall barn equipped with automatic feed bunks that control individual access and weight feed intake (CRFI, Biocontrol, Rakkestad, Norway). Cows had free access to water and were milked twice daily at approximately 6:30 and 16:00. The nutrient composition of the experimental diets was previously described [18].

MG samples were collected by biopsies for miRNome analyses on day 0 and day 6 relative to initiation of feed restriction. Mammary biopsies were performed as previously described
[19]. Briefly, the biopsy site was selected at a midpoint on alternate rear quarters of mammary gland. Approximately 600–650 mg of tissue was collected, rinsed in sterile 0.9% saline solution, snap-frozen in liquid nitrogen, and stored at –80°C until RNA extraction.

Milk samples were collected on day -1 and on day 5 relative to initiation of feed restriction. Milk composition was determined by mid-infrared spectroscopy (LIAL, Aurillac, France), and their means were computed according to PM/AM production and composition. Milk production and energy balance data were analyzed using mixed models of SAS (version 9.4; SAS Institute Inc., Cary, NC) that included day, breed, and their interaction as fixed effects, cow as random effect, and Kenward-Roger adjustment for calculation of degrees of freedom.

## RNA preparation

Total RNA was extracted from approximately 50 mg of MG biopsies using the miRVana kit (Thermo Fisher Scientific, USA). The concentration and quality of the RNAs were estimated by spectrophotometry Nanodrop TH (ND-1000, NanoDrop Technologies LLC, Wilmington, DE, USA) and using the Bioanalyzer 2100 (Agilent Technologies Inc., Santa Clara, CA, USA), respectively. The means of the obtained RIN were 8.3 ± 0.6 and 8.3 ± 0.9 in HOLS and MONT, respectively.

## RNA-seq analysis

RNA-seq analyses were performed on RNAs from MG of five HOLS and six MONT cows before and at the end of the restriction period, by the IGBMC Microarray and Sequencing Platform (Strasbourg, France), as previously described [15]. Briefly, small RNA-seq libraries were generated from 2 µg of total RNA using the TruSeq Small RNA Library Prep Kit (Illumina, San Diego, CA) according to the manufacturer's instructions. Specific miRNA and small RNA adapters were added to each end of the miRNA followed by reverse transcription and PCR amplification to obtain cDNA. Acrylamide gel purification of small RNAs was performed. The final cDNA libraries were sequenced on HiSeq 4000 (Illumina) as single-end 50 base reads according to the manufacturer's instructions. Adaptor sequences and reads with undetected bases were removed (FASTX-Toolkit, http://hannonlab.cshl.edu/fastx_toolkit/index.html), and all remaining reads were filtered according to their size. Read quantification and annotation were performed using the ncPRO-seq pipeline [20]. The sequence reads were aligned against the *Bos taurus* btau5.0.1 genome as miRBase_v22.1 using the miRDeep2 package [21]. Precursors and mature miRNAs were identified using the miRDeep2 core module, miRDeep2.pl. Potential miRNAs were annotated accordingly against orthologous miRNAs in goat, sheep and humans (miRBase release 22.1). We used a miRDeep2 score ≥ 0 as a cut-off threshold. The normalization and differential expression analyses of miRNAs were conducted using the DESeq2 R package v1.18.1. Significance was considered at $P_{adj} \le 0.05$. The accession number of the RNA-seq data is GSE146762.

## Microarray analyses

Due to the miRNA results, mRNA analyses were performed only in HOLS. Fifty nanograms of total RNA from the MG samples (n = 12) already used for miRNA analyses were analyzed using a three bovine 4 × 44K microarray (Agilent Technologies, Inc., Santa Clara, CA, USA) leading to 12 arrays. All procedures were performed according to the manufacturer's instructions. Briefly, total RNA was linearly amplified and labeled with Cy3 using the one-color Low Input Quick Amp Labeling Kit (Agilent Technologies, Inc. Santa Clara, CA, USA). Cy3-labeled cRNA (1650 ng) was hybridized using the Gene Expression Hyb Kit (Agilent Technologies, Inc., Santa Clara, CA, USA). Hybridization was performed for 17 h at 65°C in a rotating
hybridization oven at 10 rpm. Following hybridization, all microarrays were washed and scanned using the Agilent Microarray Scanner G2565A (Agilent Technologies, Inc. Santa Clara, CA, USA). The resulting TIFF images (Tagged Image File Format) were processed using Feature Extraction software Version 11. Microarrays were analyzed to determine differential gene expression (DEG) at the mRNA level using GeneSpring software, and data were normalized with 75th percentile shift followed by *t*-tests moderated with the Benjamini-Hochberg (FDR) correction [22], and significance was considered at $P_{adj} < 0.10$. Microarray data were deposited in the Gene Expression Omnibus (GEO): GSE159302.

## Bioinformatic analyses

The genomic localization of DEMs was compared with the genomic location of milk-associated QTL. All bovine QTL downloaded from CattleQTLdb (AnimalQTLdb release 40, December 29, 2019, https://www.animalgenome.org) were obtained from the "All data by bp (on UMD3.1)" link in bed format. These localizations were used to identify the mammary miRNAs hosted in these regions. The genomic positions of the miRNAs in the UMD3.1 reference genome permitted their screening in the QTL regions using a homemade bioinformatics script in Python language (S1 Fig). The QTL related to milk production and composition were then extracted.

Target genes of DEMs and corresponding putative biological processes ($P \le 0.05$) were investigated using the miRWalk database (version 3.0) and Enrichr platform [23, 24], integrating several tools for pathway visualizations, including KEGG, BioCarta, for more efficient analyses of genes in networks. Modulated biological processes were grouped using the REVIGO platform (http://revigo.irb.hr/) [25]. Metacore™ software (release 6), which uses text mining algorithms to identify interactions and build pathways/networks, was also used to corroborate bioinformatics investigations. However, these algorithms do not use the exact representativeness corresponding to the background of the used whole genome microarray. A Venn diagram was obtained using Venny 2.1 (http://bioinfogp.cnb.csic.es/tools/venny/) to find common genes between miRNA target genes and DEGs identified by microarray analyses. Biological processes ($P \le 0.05$) of these common genes were investigated using the Enrichr platform for both up and down-regulated genes together.

## Results and discussion

### Animal responses and milk analyses

Energy balance, milk production and composition were analyzed to characterize cow response to the feed restriction model (Table 1). Feed restriction induced NEB in both HOLS and MONT cows (-35 and -29 ± 4 MJ/day, respectively) due to the experimental design. Milk production, fat and protein yields decreased by 40, 33 and 43% in HOLS and 43, 34 and 49% in MONT, respectively, during REST (Table 1). Production and metabolic responses are detailed elsewhere [18]. The effects of restriction on milk production are in agreement with previous studies in mid- and late-lactation HOLS cows [26, 27].

A significant breed effect was observed for energy balance and production variables. Montbéliarde cows had lower milk, fat and protein yields than HOLS cows, which is in agreement with previous studies [28, 29]. A significant period by breed interaction was observed for energy balance and because energy balance was lower in MONT during the CONT period compared to HOLS cows. This could be due to lower intake in MONT than in HOLS cows during this period [18].
<table><caption>Table 1. Effects of feed restriction on energy balance, DMI, milk, milk fat and protein yields in midlactation Holstein (HOLS) and Montbéliarde (MONT) cows during control (CONT) and restriction (REST) periods.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Breed</th><th colspan="2">Period</th><th rowspan="2">SEM</th><th colspan="3">P-value</th></tr><tr><th>CONT</th><th>REST</th><th>Breed</th><th>Period</th><th>Breed x Period</th></tr></thead><tbody><tr><td rowspan="2">Energy balance (MJ/d)</td><td>HOLS</td><td>56.5<sup>a</sup></td><td>-35.1</td><td>4.2</td><td rowspan="2">0.03</td><td rowspan="2">&lt;0.001</td><td rowspan="2">&lt;0.001</td></tr><tr><td>MONT</td><td>26.0<sup>b</sup></td><td>-28.9</td><td>3.9</td></tr><tr><td rowspan="2">Milk yield (kg/d)</td><td>HOLS</td><td>30.8</td><td>18.5</td><td>1.0</td><td rowspan="2">&lt;0.001</td><td rowspan="2">&lt;0.001</td><td rowspan="2">0.38</td></tr><tr><td>MONT</td><td>25.1</td><td>14.2</td><td>0.9</td></tr><tr><td rowspan="2">Fat yield (g/d)</td><td>HOLS</td><td>1 030.1</td><td>694.4</td><td>32.1</td><td rowspan="2">0.005</td><td rowspan="2">&lt;0.001</td><td rowspan="2">0.62</td></tr><tr><td>MONT</td><td>896.2</td><td>590.3</td><td>29.3</td></tr><tr><td rowspan="2">Protein yield (g/d)</td><td>HOLS</td><td>929.3</td><td>529.1</td><td>33.9</td><td rowspan="2">0.008</td><td rowspan="2">&lt;0.001</td><td rowspan="2">0.82</td></tr><tr><td>MONT</td><td>795.6</td><td>407.4</td><td>30.9</td></tr><tr><td rowspan="2">Lactose yield (g/d)</td><td>HOLS</td><td>1 565.3</td><td>970.7</td><td>50.8</td><td rowspan="2">&lt;0.001</td><td rowspan="2">&lt;0.001</td><td rowspan="2">0.96</td></tr><tr><td>MONT</td><td>1 263.6</td><td>665.8</td><td>44.3</td></tr></tbody></table>

SEM corresponds to the standard error of the mean. a, b correspond to breed LSMEANS not sharing a common superscript differ between CONT and REST periods ($P \le 0.05$) and are presented when Breed x Period effect was significant.

https://doi.org/10.1371/journal.pone.0248680.t001

## Mammary gland miRNome responses

**Summary of miRNA sequencing data.** Ten and twelve libraries were constructed using RNA extracted from the MG of HOLS and MONT cows at CONT and at REST periods. More than 12 million raw reads and more than 8.5 million clean reads were obtained after removing adapters, no-clipped reads and reads that were too short. The percentage of clean reads mapped into the *Bos taurus* genome was 97.4 on average for all libraries. In HOLS and MONT cows, the log₂ of the normalized counts of the CONT and REST libraries had a strong correlation of R ≥ 0.96. MiRdeep2 software analysis identified 754 miRNAs (737 known and 17 predicted) for HOLS libraries and 759 miRNAs (743 known and 16 predicted) for MONT libraries. These results are in accordance with other studies on bovine MG using RNA-seq technology, reporting 654 [9] and 324 [11] known miRNAs. However, these studies presented a higher number of predicted miRNAs (679 in Le Guillou et al. [9] and 176 in Li et al. [11]) than in our study. The difference could be due to the version of miRBase used (miRBase V20, V21 and V22 in Le Guillou et al. [9], Li et al. [11] and our study, respectively). In the present study, the output ratio of CONT compared to the REST period was closer in HOLST than in MONT libraries (Δ1.6% and Δ6.4%, respectively; Table 2).

**Differentially expressed miRNAs between the control and restriction periods.** The expression of 27 miRNAs was significantly modified by feed restriction in HOLS cows ($P_{adj} \le 0.05$; Table 3), but no significant miRNA change was observed in MONT cows. The absence of DEMs in MONT cows may be due to the lower quality of their libraries compared to those of

<table><caption>Table 2. Effect of feed restriction on miRNA RNA-seq data in mammary glands from midlactation Holstein (n = 5) and Montbéliarde (n = 6) cows in control (CONT) and restriction (REST) periods.</caption><thead><tr><th>Breed</th><th>Condition</th><th>Raw read</th><th>Too short</th><th>Adapter</th><th>Clean read</th><th>Output ratio (%)</th></tr></thead><tbody><tr><td rowspan="2">Holstein</td><td>CONT</td><td>13,393,626</td><td>724,538</td><td>56,192</td><td>11,027,924</td><td>82.2</td></tr><tr><td>REST</td><td>17,101,298</td><td>794,008</td><td>61,674</td><td>14,352,279</td><td>83.8</td></tr><tr><td rowspan="2">Montbéliarde</td><td>CONT</td><td>11,953,902</td><td>1,405,670</td><td>824,394</td><td>8,500,865</td><td>71.0</td></tr><tr><td>REST</td><td>14,841,154</td><td>1,354,192</td><td>622,388</td><td>11,539,527</td><td>77.4</td></tr></tbody></table>

Clean reads were obtained after size (18 to 25 nt), adapter and soft-clipping cleaning.

https://doi.org/10.1371/journal.pone.0248680.t002
<table><caption>Table 3. Differentially expressed miRNAs in the mammary gland of midlactation Holstein cows during restriction (REST) compared to the control (CONT) period.</caption><thead><tr><th rowspan="2">Name</th><th colspan="2">Mean read</th><th rowspan="2">Expression ratio</th><th rowspan="2">FDR</th><th rowspan="2">Expression change</th></tr><tr><th>CONT</th><th>REST</th></tr></thead><tbody><tr><td>bta-miR-143</td><td>2,292,707.1</td><td>3,119,812.5</td><td>1.37</td><td>1.4E-02</td><td>Up</td></tr><tr><td>bta-miR-181a</td><td>197,294.7</td><td>270,887.9</td><td>1.37</td><td>2.9E-02</td><td>Up</td></tr><tr><td>bta-miR-26b</td><td>97,181.9</td><td>138,114.2</td><td>1.43</td><td>1.4E-02</td><td>Up</td></tr><tr><td>bta-miR-200c</td><td>65,416.1</td><td>48,881.0</td><td>0.75</td><td>3.6E-02</td><td>Down</td></tr><tr><td>bta-miR-25</td><td>19,640.0</td><td>24,430.3</td><td>1.25</td><td>1.7E-02</td><td>Up</td></tr><tr><td>bta-miR-200b</td><td>17,718.1</td><td>12,055.1</td><td>0.68</td><td>1.2E-02</td><td>Down</td></tr><tr><td>bta-miR-181b</td><td>8,203.8</td><td>11,208.3</td><td>1.37</td><td>2.6E-02</td><td>Up</td></tr><tr><td>bta-miR-155</td><td>1,178.0</td><td>1,724.1</td><td>1.47</td><td>3.9E-02</td><td>Up</td></tr><tr><td>bta-miR-193a-3p</td><td>889.5</td><td>597.8</td><td>0.67</td><td>2.2E-06</td><td>Down</td></tr><tr><td>bta-miR-669</td><td>695.1</td><td>537.3</td><td>0.78</td><td>2.6E-02</td><td>Down</td></tr><tr><td>bta-miR-374b</td><td>694.3</td><td>533.5</td><td>0.77</td><td>9.3E-03</td><td>Down</td></tr><tr><td>bta-miR-379</td><td>377.7</td><td>526.7</td><td>1.39</td><td>5.3E-02</td><td>Up</td></tr><tr><td>bta-miR-222</td><td>266.7</td><td>370.5</td><td>1.39</td><td>2.2E-03</td><td>Up</td></tr><tr><td>bta-let-7a-3p</td><td>348.6</td><td>207.4</td><td>0.60</td><td>3.3E-03</td><td>Down</td></tr><tr><td>bta-miR-326</td><td>198.5</td><td>138.9</td><td>0.70</td><td>5.8E-03</td><td>Down</td></tr><tr><td>bta-miR-2898</td><td>199.2</td><td>117.8</td><td>0.59</td><td>5.8E-03</td><td>Down</td></tr><tr><td>bta-miR-500</td><td>174.1</td><td>104.1</td><td>0.60</td><td>9.3E-03</td><td>Down</td></tr><tr><td>bta-miR-29d-5p</td><td>98.8</td><td>67.1</td><td>0.68</td><td>3.5E-02</td><td>Down</td></tr><tr><td>bta-miR-495</td><td>70.5</td><td>26.2</td><td>0.37</td><td>1.3E-08</td><td>Down</td></tr><tr><td>bta-miR-299</td><td>74.9</td><td>19.4</td><td>0.26</td><td>1.3E-08</td><td>Down</td></tr><tr><td>bta-miR-299-2</td><td>71.9</td><td>17.7</td><td>0.25</td><td>1.3E-08</td><td>Down</td></tr><tr><td>bta-miR-671</td><td>50.8</td><td>28.4</td><td>0.56</td><td>5.8E-03</td><td>Down</td></tr><tr><td>chr7_8584_mature</td><td>55.5</td><td>17.5</td><td>0.32</td><td>6.9E-04</td><td>Down</td></tr><tr><td>bta-miR-296-5p</td><td>47.8</td><td>25.0</td><td>0.52</td><td>9.5E-03</td><td>Down</td></tr><tr><td>chr10_11366_mature</td><td>31.3</td><td>12.6</td><td>0.40</td><td>3.5E-03</td><td>Down</td></tr><tr><td>bta-miR-18b</td><td>29.3</td><td>12.5</td><td>0.44</td><td>1.4E-02</td><td>Down</td></tr><tr><td>bta-miR-2483-5p</td><td>17.2</td><td>6.6</td><td>0.40</td><td>2.4E-02</td><td>Down</td></tr></tbody></table>

Means correspond to the mean of normalized read counts. The bold and italic line corresponds to all miRNAs with more than 1,000 reads and known miRNAs, respectively. FDR: False discovery rate.

https://doi.org/10.1371/journal.pone.0248680.t003

HOLS cows (with a greater number of short reads and adapters and thus a smaller output ratio in MONT; Table 3). Another hypothesis is that MONT cows presented a lower energy value during the CONT period than HOLS cows. Montbéliarde cows presented a lower variation of energy balance between the CONT and REST periods (26 MJ in CONT and -28.9 MJ in REST: Δ54.6 MJ) compared to HOLS cows (56.5 MJ in CONT and -35.1 MJ in REST: Δ91.6 MJ). In addition, MONT cows had a lower decrease in the milk, fat, and protein yields than HOLS. Moreover, decades of genetic selection in HOLS have reduced genetic variation in this breed, which may have allowed detecting significant differences according to the diet. These hypotheses or a combination of these may partly explain the absence of significant modification of miRNA gene expression. Thus, the following results and discussion are specific to Holstein breed.

Among the 27 DEMs affected by the restriction in HOLS, 25 are known and 2 predicted with 19 down- and 8 upregulated by feed restriction. Eight DEMs were highly expressed with more than 1,000 reads (Table 3). Thirty DEMs were modulated by 48 h of feed deprivation in
the MG of lactating goats [10]. The number of DEMs in the current study is similar to those reported in the study of Mobuchon et al. [10] using similar analyses expected that the cut-off was 0.1 in the study performed by Mobuchon et al. and 0.05 in the present study. Two miRNAs are common to our study: *miR-222*, which was upregulated with a similar fold change (FC) between studies (FC: 1.47 vs 1.39 in caprine and bovine, respectively), and *miR-671*, which was downregulated (FC: 0.73 vs 0.56 in caprine and bovine, respectively). However, whereas *miR-671* presented a similar level of expression in bovine and caprine MG, *miR-222* showed more reads in caprine (812 during control and 1,408 in feed deprived goats) than in bovine (267 in CONT and 371 in REST), suggesting a higher level of *miR-222* expression in caprine MG compared to bovine.

Early lactation is known to be a natural period of NEB for dairy ruminants. A comparison of DEMs obtained during an induced NEB in cows (present study) and during early lactation in goats [30] highlighted similarities between these two species. For instance, Ji et al. [30] observed 378 DEMs in goat MG between early and late lactation periods using all mammalian miRNA sequences for mapping. Among them, 4 DEMs were common to the present study: *miR-181a* and *-299*, *-2483-5P* (mapped on ruminant genome) and *miR-200c* (mapped on Tasmanian devil [30]). Except for *miR-200c*, all other miRNAs presented a similar manner of regulation as our study. Among the common miRNAs, *miR-181a* and *miR-200c* are highly expressed miRNAs in our study. Wang et al. [7] showed that the expression of 13 miRNAs in bovine MG was greater during NEB (i.e., early lactation) than during the positive energy balance period (i.e., dry period). Despite the dissimilarity of physiological states between Wang et al. [7] (early vs non-lactation period) and our study, two differentially expressed miRNAs (*miR-155* and *-181a*) are common and have a similar sense of regulation by energy balance. These miRNAs belong to the highly expressed DEMs in our study.

**Differentially expressed miRNAs located in milk QTL.** In cattle, at least 50,178 quantitative trait loci (QTL) out of a total of 130,471 have been associated with milk production traits (https://www.animalgenome.org). Analyses of genome positions of the 25 known DEMs in the current study revealed that 8 DEMs are associated with a milk QTL (Table 4). One (*miR-155*) was within a QTL associated with milk α-lactalbumin content, and 7 miRNAs were associated with milk fat composition (Table 4). Among these 7 miRNAs, 5 (*miR-200c*, *-2898*, *-181a*, *-181b*, and *-296-5p*) were located in QTL associated with saturated fatty acid content, and 4 (*miR-26b*, *-200c*, *-296-5p*, and *-669*) were located in QTL linked to monounsaturated fatty acid content. The difference in the number of miRNAs observed in QTL of saturated vs. unsaturated fatty acids may be due to the difference in the abundance of QTL associated with two

<table><caption>Table 4. Location of nutriregulated miRNAs in QTL associated with milk production and composition.</caption><thead><tr><th>miRNA name</th><th>Chromosome</th><th>Start</th><th>End</th><th>Character</th></tr></thead><tbody><tr><td>bta-miR-155</td><td>1</td><td>10227315</td><td>10227337</td><td>Milk alpha-lactalbumin percentage</td></tr><tr><td>bta-miR-26b</td><td>2</td><td>107133408</td><td>107133429</td><td>Milk oleic acid content</td></tr><tr><td>bta-miR-200c</td><td>5</td><td>103859438</td><td>103859460</td><td>Milk myristic acid content; Milk palmitoleic acid content</td></tr><tr><td>bta-miR-2898</td><td>8</td><td>74354027</td><td>74354044</td><td>Milk capric acid content; Milk myristic acid content</td></tr><tr><td>bta-miR-181a</td><td>11</td><td>95709449</td><td>95709472</td><td>Milk butyric acid content; Milk caproic acid content</td></tr><tr><td>bta-miR-181b</td><td>11</td><td>95710641</td><td>95710664</td><td>Milk butyric acid content; Milk caproic acid content</td></tr><tr><td>bta-miR-296-5p</td><td>13</td><td>58091948</td><td>58091966</td><td>Milk capric acid content; Milk caproic acid content; Milk caprylic acid content; Milk myristoleic acid content; Milk palmitoleic acid content</td></tr><tr><td>bta-miR-669</td><td>16</td><td>63239604</td><td>63239625</td><td>Milk palmitoleic acid content</td></tr></tbody></table>

Characters are from https://www.animalgenome.org.

https://doi.org/10.1371/journal.pone.0248680.t004
classes of fatty acids [31]. In bovines, the close genome positions of *miR-181a* and *miR-181b* (< 10 kb; referred to miRBase criteria) on chromosome 11 suggest that they are in a cluster. This finding is in line with the results in mice reporting that they belong to the same family and are in a cluster [32]. Their cluster is located within the same QTL associated with the regulation of butyric and caproic acid content in milk (Table 4). Moreover, the *miR-181a* and *miR-181b* upregulation and their high level of expression during the REST period may be related, at least partly, to the decrease in short chain fatty acid content. Indeed, in our study, ΣC4 to C9 was 5.41 vs 3.43 g/100 g total fatty acids during CONT and REST (S1 Table), respectively (P < 0.001).

**Analyses of the 8 highly expressed and differential miRNAs.** We focused on highly expressed miRNAs to decipher the mechanisms underlying the response of MG to restriction. Among the 25 known DEMs in Holstein cows, 8 (*miR-143*, *miR-181a*, *miR-26b*, *miR-200c*, *miR-25*, *miR-200b*, *miR-181b*, and *miR-155*) presented more than 1,000 reads, representing 99% of the total reads of the DEMs (Table 3). *MiR-143* was the most abundant miRNA in this study, possessing more than 2 million reads, which is in agreement with previous studies in bovine MG [33] and bovine liver [34, 35]. *MiR-143*, known to regulate lipid metabolism, was downregulated in bovine liver during severe NEB [34]. In the present study, *miR-143* was upregulated in MG during the NEB period. This difference in regulation may be due to the studied tissue (liver vs MG). Indeed, the NEB period is characterized by strong adipose tissue mobilization, increased NEFA esterification, triacylglycerol accumulation, and ketone synthesis in the liver, as well as decreased lipid metabolism in the MG [36].

*MiR-26b* and *-181a* are, respectively, the second and third most expressed miRNAs in our study, possessing more than 100,000 reads. *MiR-26b* was identified as playing a role in inflammation [37]. The upregulation of *miR-26b* in our study might be linked to inflammation due to restriction. Indeed, lipid mobilization during NEB may lead to inflammation [38]. *MiR-181a* was reported to influence lipid and glucose metabolism in humans, mice and bovine liver [39]. *MiR-181a* targets *ACSL1* mRNA, which is involved in lipid synthesis in bovine mammary epithelial cells [13]. The upregulation of *miR-181a* in our study may be related to the decrease in milk fat yield observed during the NEB period (Table 1) and the diminution of short- (data not shown) and medium-chain fatty acids [18]. Ma et al. [40] showed that *miR-25* was involved in triacylglycerol synthesis and lipid accumulation in goat mammary epithelial cells. This miRNA directly targets *PPARG* and regulates the expression of *SREBP1*, *FASN*, and *GPAM* involved in lipid synthesis. The observed upregulation of *miR-181a* and *miR-25* may be related to the decrease in fat yield during the feed restriction period. Among the highly expressed DEMs, *miR-200b* and *miR-200c* belong to the same family, have the same seed sequence, and regulate epithelial-to-mesenchymal transition in breast cancer [41, 42]. The downregulation of *miR-200b* and *miR-200c* in our study may suggest a link with MG structural changes due to feed restriction [43].

*MiR-181b* and *miR-155* were the least expressed of the 8 miRNAs and were upregulated by feed restriction. The modulation of their expression may be linked to the stress induced by feed restriction. Indeed, *miR-181b* was upregulated in the serum of heat-stressed Frieswal bulls compared with summer and winter periods [44] and after chirurgical stress in Angus cattle muscle [45]. In addition, *miR-181b* and *miR-155* were also reported to play a role in the immune response [46, 47]. LPS stimulation of mouse macrophages *in vitro* resulted in the upregulation of *miR-155*, and similar changes also occurred after intraperitoneal LPS injection in C57BL/6 mice [47]. In addition, *miR-155* promoted the production of the proinflammatory cytokine TNF-α in human B cells [47]. Thus, the upregulation of *miR-181b* and *miR-155* observed in the present study may be explained by inflammation that can be affected by NEB, as demonstrated in other studies [38].
Differentially expressed mRNAs in Holstein cows

To elucidate the mechanisms of nutritional regulation of MG function, mRNA analysis was performed. Considering the results obtained by miRnome studies, mRNA analysis was achieved only in HOLS cows. Microarray analyses identified 374 differentially expressed genes (DEGs, P<0.1) between the CONT and REST periods with 120 upregulated and 254 down-regulated genes during the REST period. Ten DEGs presented an FC ≥ 3. In caprine MG, Ollier et al. [2] identified 161 DEGs whose expression was altered by 48 h of feed deprivation using the same bovine microarray. The lower number of DEGs reported by Ollier et al. [2] might be due to the heterogeneous system used (caprine mRNA with a bovine array). Bioinformatic analysis of the 374 DEGs showed 256 altered biological processes (BPs; P<0.05). The analyses of the first 25 BPs revealed that 12 BPs were linked to lipid metabolism and involved 28 DEGs (23 downregulated and 5 upregulated), and 9 BPs were related to respiration metabolism and involved 33 DEGs (30 downregulated and 3 upregulated; Table 5). In addition, metabolic process (S2 Table) and SCAP/SREBP, cholesterol biosynthesis and regulation of lipid metabolism (S3 Table) were also highlighted using Metacore™ analysis. Despite the fact that the results were corroborated between different bioinformatic analyses, it should be noted that

<table><caption>Table 5. Top 25 biological processes potentially regulated by the 374 DEGs in Holstein cows.</caption><thead><tr><th></th><th>Term</th><th>Overlap</th><th>P-value</th><th>FDR</th></tr></thead><tbody><tr><td>%</td><td>mitochondrial ATP synthesis coupled electron transport</td><td>27/86</td><td>4.7E-26</td><td>8.9E-23</td></tr><tr><td>%</td><td>respiratory electron transport chain</td><td>27/95</td><td>9.6E-25</td><td>9.1E-22</td></tr><tr><td>%</td><td>mitochondrial electron transport, NADH to ubiquinone</td><td>21/47</td><td>2.2E-24</td><td>1.4E-21</td></tr><tr><td>%</td><td>mitochondrial respiratory chain complex I biogenesis</td><td>20/65</td><td>2.1E-19</td><td>6.5E-17</td></tr><tr><td>%</td><td>NADH dehydrogenase complex assembly</td><td>20/65</td><td>2.1E-19</td><td>6.5E-17</td></tr><tr><td>%</td><td>mitochondrial respiratory chain complex I assembly</td><td>20/65</td><td>2.1E-19</td><td>6.5E-17</td></tr><tr><td>%</td><td>mitochondrial respiratory chain complex assembly</td><td>22/98</td><td>5.9E-18</td><td>1.6E-15</td></tr><tr><td>*</td><td>cholesterol biosynthetic process</td><td>11/36</td><td>3.2E-11</td><td>7.7E-09</td></tr><tr><td>*</td><td>sterol biosynthetic process</td><td>11/41</td><td>1.6E-10</td><td>3.3E-08</td></tr><tr><td>*</td><td>secondary alcohol biosynthetic process</td><td>10/37</td><td>1.0E-09</td><td>1.9E-07</td></tr><tr><td>*</td><td>regulation of cholesterol biosynthetic process</td><td>9/41</td><td>5.1E-08</td><td>8.8E-06</td></tr><tr><td>*</td><td>regulation of cholesterol metabolic process</td><td>9/42</td><td>6.4E-08</td><td>1.0E-05</td></tr><tr><td>*</td><td>regulation of alcohol biosynthetic process</td><td>8/35</td><td>2.1E-07</td><td>3.0E-05</td></tr><tr><td>*</td><td>long-chain fatty acid transport</td><td>7/26</td><td>3.6E-07</td><td>4.9E-05</td></tr><tr><td>*</td><td>cholesterol metabolic process</td><td>10/69</td><td>5.8E-07</td><td>7.3E-05</td></tr><tr><td>*</td><td>regulation of steroid biosynthetic process</td><td>8/45</td><td>1.6E-06</td><td>1.9E-04</td></tr><tr><td>*</td><td>sterol metabolic process</td><td>8/60</td><td>1.5E-05</td><td>1.7E-03</td></tr><tr><td>%</td><td>mitochondrial transmembrane transport</td><td>6/33</td><td>2.9E-05</td><td>3.1E-03</td></tr><tr><td>*</td><td>fatty acid metabolic process</td><td>10/107</td><td>3.2E-05</td><td>3.2E-03</td></tr><tr><td>*</td><td>carnitine shuttle</td><td>4/12</td><td>5.2E-05</td><td>4.7E-03</td></tr><tr><td>*</td><td>fatty acid transmembrane transport</td><td>4/12</td><td>5.2E-05</td><td>4.7E-03</td></tr><tr><td>*</td><td>lipid biosynthetic process</td><td>8/73</td><td>6.3E-05</td><td>5.5E-03</td></tr><tr><td>%</td><td>cellular respiration</td><td>7/58</td><td>9.8E-05</td><td>8.1E-03</td></tr><tr><td>*</td><td>intracellular lipid transport</td><td>4/16</td><td>1.8E-04</td><td>1.4E-02</td></tr><tr><td>*</td><td>fatty-acyl-CoA biosynthetic process</td><td>5/31</td><td>2.5E-04</td><td>1.9E-02</td></tr></tbody></table>

Process network analysis was performed using Enrichr software.

% indicates processes associated with respiration, and * indicates lipid metabolism processes. The overlap corresponds to the ratio of the DEG number on the number of genes known in the BP. FDR: false discovery rate.

https://doi.org/10.1371/journal.pone.0248680.t005
a limitation of these algorithms is that they did not consider the background corresponding to
the exact gene representation of the whole genome microarray used.

The effects of feed restriction on genes involved in lipid metabolism in the MG of lactating Holstein cows (Fig 1) were in agreement with previous studies. The downregulation of *ACACA*, *FABP3* and *LPL* in our study was in accordance with the reduction of their expression after 4 days of restriction (cows fed at 60% of their *ad libitum* DMI; [48]). The downregulation of *ACSL1*, *FABP3*, and *LPIN1* in mammary tissue of dairy cows is in line with the decrease in their expression during the natural NEB period of early lactation (day 15 of lactation) compared to midlactation (day 60 of lactation), which is associated with a lower fat yield [49]. Taken together, these results indicate that the downregulation of the DEGs (*ACACA*, *FABP3*, *LPL*, *LPIN1* and *ACSL1*) involved in lipid metabolism observed in our study (Fig 1) could be related to the decrease in milk fat yield observed during REST.

In the present study, feed restriction downregulated 21 DEGs (e.g., *NDUFA3*, *NDUFB1*, and *NDUFS2*) which are related to the NADH pathway and respiration (Table 5, S2 Table). A downregulation of several NADH dehydrogenase mRNAs was previously observed in bovine liver during heat stress [50], suggesting that feed restriction may affect energy production through mitochondrial respiration, which was already identified among the most inhibited functions in bovine liver [51].

Analysis of mRNAs common between the DEG list and mRNAs targeted by
the 8 major DEMs

Using miRWalk analyses, we identified 2934 mRNAs potentially targeted by the 8 most expressed DEMs in Holstein cows. A comparison between those potential targets and the 374 DEGs identified by microarray analysis showed 59 common mRNAs (Fig 2). Among these mRNAs, 37 were downregulated, and 22 were upregulated. Two DEGs (*LPIN1* and *LPL*) had an absolute FC ≥ 3, and both were downregulated. *LPIN1* and *LPL* are involved in lipid metabolism [48, 49]. Analysis of the 59 common mRNAs by the Enrichr platform identified 266 BPs (P ≤ 0.05). The elimination of redundancy using the REVIGO platform highlighted 146 BPs, which were regrouped into 13 categories denoting a semantic similarity between BPs (Figs 2 and 3), and 6 of these categories represented 89% of BPs. Two categories (positive regulation of endothelial cell proliferation and lipid biosynthesis) contained more than 20% BPs each. Positive regulation of endothelial cell proliferation was the first category, with 38 of the 146 BPs being regulated. Among those 38 BP, 19 are predicted to be regulated by PROX1 (Prospero Homeobox 1; S4 Table), which is a downregulated DEG that is targeted by *miR-181a*, a highly expressed miRNA that was determined to be upregulated by REST in HOLS cows. PROX1 is a transcriptional factor influencing cell fate, gene transcriptional regulation and progenitor cell regulation, playing a role in lymphatic system development. Wilge and Oliver [53] showed that knockouts in *PROX1* mouse embryos lacked lymphatic system development. Elsewhere, a decrease in lymphatic vasculature (identified using *PROX1* and *CD31*) in MG was observed across the lactation of mice [54]. In our study, the upregulation of *miR-181a* and the downregulation of *PROX1* by feed restriction suggest a modification of lymphatic vasculature. In addition, *PROX1* was hosted in a QTL associated with milk palmitoleic acid content and cheese-making properties [55]. The modulation of *miR-181a* and *PROX1* gene expression is consistent with the milk fatty acid changes observed during feed restriction [18].

Lipid biosynthesis had 29 of the 146 BPs and was the second category identified by REVIGO analysis. Furthermore, three BPs (lipid biosynthetic process, triglyceride biosynthetic process and fatty acid metabolic process) were also identified among the top 15 BPs identified by Enrichr analyses. In these three BPs, four DEGs (*ACACA*, *CYP51A1*, *LPL*, and *LPIN1*) were
<figure><img src="image_4.png" /><figcaption>Fig 1. MiRNAs and coding genes involved in the milk fat synthesis pathway. Genes in the red circle correspond to DEGs identified in our study, and miRNAs in the purple box correspond to DEMs targeting the DEGs (adapted from Bionaz and Loor [52]).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0248680.g001
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Comparison between mRNA targeted by the 8 most expressed DEMs and DEGs in Holstein cows. Venn diagram between the potential targets of DEMs (blue circle) and the 374 DEGs (yellow circle). Red arrows indicate the workflow used to analyze biological processes (BPs) common mRNA with the Enrichr and REVIGO platforms.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0248680.g002

<figure><img src="image_6.png" /><figcaption>**Fig 3.** REVIGO categories obtained from the Enrichr biological process list using the 59 common mRNAs (P ≤ 0.05). In brackets: Numbers of BPs in each category.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0248680.g003
downregulated by restriction. A downregulation of *CYP51A1*, a gene involved in cholesterol synthesis, was already reported during feed restriction in chicken liver [56] and in adipose tissue of beef cows [57]. ACACA, *CYP51A1*, *LPL*, and *LPIN1* are predicted to be targeted by 4 DEMs (*miR-155*, *-181a*, *-181b*, and *-200b*). *MiR-181a* and *miR-155* target *LPL*, *miR-181b* target ACACA, and *miR-200b* target *LPIN1* and *CYP51A1*. The upregulation of *miR-181a*, *-181b*, and *-155* was in line with the downregulation of their targeted genes. A study demonstrated that the absence of *miR-155* led to an increase in *LPL* expression in the liver of *miR-155*⁻/⁻ mice [58].

## Conclusions

Feed restriction modified MG miRNome and transcriptome in Holstein but not in Montbéliarde cows. The differentially expressed miRNAs highlight a potential role of miRNAs in MG remodeling and lipid synthesis. Moreover, these roles were supported by the analyses of common differential mRNAs and predicted targets of miRNAs sharing the same function. In addition, the location of differentially expressed miRNAs in QTL associated with milk fatty acid composition was consistent with the observed reduction in short- and medium-chain milk fatty acid content in milk. The decrease in *de novo* milk fatty acid synthesis observed during restriction mediated in part by the regulation of the miRNAs affecting mRNAs involved in lipid metabolism and suggest a role of DEMs in this biological process.

## Supporting information

S1 Fig. Homemade bioinformatics script in Python language for screening QTL.
(TIF)

S1 Table. Effects of feed restriction on milk fatty acid concentrations (g/100 g of FA) in mid-lactation cows.
(DOCX)

S2 Table. TOP 20 of the biological processes affected by restriction in mammary gland of mid-lactation cows. Enrichment by GO Processes using Metacoret™ software using the 374 differentially expressed genes after restriction. FDR: false Discovery Rate.
(DOCX)

S3 Table. TOP 20 of the pathways affected by restriction in mammary gland of mid-lactation cows. Enrichment by Pathway maps using Metacoret™ software using the 374 differentially expressed genes after restriction. FDR: false Discovery Rate.
(DOCX)

S4 Table. Biological processes predicted to be regulated by Prox1 identified by bioinformatics analyses of the 59 genes differentially expressed and predicted to be targeted by the 8 major differentially expressed miRNAs highlighted the regulation of endothelial cell proliferation containing 38 biological processes among which 19 were predicted to be regulated by Prox1.
(DOCX)

## Acknowledgments

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

**Author Contributions**

**Conceptualization:** Yannick Faulconnier, José A. A. Pires, Christine Leroux.

**Formal analysis:** Pierre-Alexis Billa, Tao Ye, Céline Bourdon.

**Funding acquisition:** Yannick Faulconnier, José A. A. Pires, Christine Leroux.

**Investigation:** Pierre-Alexis Billa, Yannick Faulconnier, José A. A. Pires, Christine Leroux.

**Project administration:** Christine Leroux.

**Resources:** Yannick Faulconnier, José A. A. Pires, Christine Leroux.

**Supervision:** Christine Leroux.

**Visualization:** Pierre-Alexis Billa.

**Writing – original draft:** Pierre-Alexis Billa, Yannick Faulconnier, Christine Leroux.

**Writing – review & editing:** Pierre-Alexis Billa, Yannick Faulconnier, José A. A. Pires, Christine Leroux.

**References**

1. Anderson SM, Rudolph MC, McManaman JL, Neville MC. Key stages in mammary gland development. Secretory activation in the mammary gland: it's not just about milk protein synthesis! Breast Cancer Research. 2007; 9: 204. https://doi.org/10.1186/bcr1653 PMID: 17338830

2. Ollier S, Robert-Granié C, Bernard L, Chilliard Y, Leroux C. Mammary Transcriptome Analysis of Food-Deprived Lactating Goats Highlights Genes Involved in Milk Secretion and Programmed Cell Death. J Nutr. 2007; 137: 560–567. https://doi.org/10.1093/jn/137.3.560 PMID: 17311940

3. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. cell. 2004; 116: 281–297. https://doi.org/10.1016/s0092-8674(04)00045-5 PMID: 14744438

4. Ambros V. The functions of animal microRNAs. Nature. 2004; 431: 350–355. https://doi.org/10.1038/ nature02871 PMID: 15372042

5. Friedman RC, Farh KK, Burge CB, Bartel DP. Most mammalian mRNAs are conserved targets of microRNAs. Genome Res. 2009;19. https://doi.org/10.1101/gr.100073.109 PMID: 19887576

6. Ameres SL, Zamore PD. Diversifying microRNA sequence and function. Nature Reviews Molecular Cell Biology. 2013; 14: 475–488. https://doi.org/10.1038/nrm3611 PMID: 23800994

7. 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. Journal of Dairy Science. 2012; 95: 6529–6535. https://doi.org/10.3168/jds.2012-5748 PMID: 22959945

8. Li Z, Lan X, Guo W, Sun J, Huang Y, Wang J, et al. Comparative Transcriptome Profiling of Dairy Goat MicroRNAs from Dry Period and Peak Lactation Mammary Gland Tissues. Qu L-H, editor. PLoS ONE. 2012; 7: e52388. https://doi.org/10.1371/journal.pone.0052388 PMID: 23300659

9. 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. Pfeffer S, editor. PLoS ONE. 2014; 9: e91938. https://doi.org/10.1371/journal.pone.0091938 PMID: 24658750

10. Mobuchon L, Marthey S, Le Guillou S, Laloë D, Le Provost F, Leroux C. Food Deprivation Affects the miRNome in the Lactating Goat Mammary Gland. Schönbach C, editor. PLOS ONE. 2015; 10: e0140111. https://doi.org/10.1371/journal.pone.0140111 PMID: 26473604

11. Li R, Beaudoin F, Ammah AA, Bissonnette N, Benchaar C, Zhao X, et al. Deep sequencing shows microRNA involvement in bovine mammary gland adaptation to diets supplemented with linseed oil or
safflower oil. BMC Genomics. 2015; 16: 884. https://doi.org/10.1186/s12864-015-1965-7 PMID: 26519053

12. Mobuchon L, Guillou SL, Marthey S, Laubier J, Laloë D, Bes S, et al. Sunflower oil supplementation affects the expression of miR-20a-5p and miR-142-5p in the lactating bovine mammary gland. PLOS ONE. 2017; 12: e0185511. https://doi.org/10.1371/journal.pone.0185511 PMID: 29281677

13. Lian S, Guo JR, Nan XM, Ma L, Loor JJ, Bu DP. MicroRNA Bta-miR-181a regulates the biosynthesis of bovine milk fat by targeting ACSL1. Journal of Dairy Science. 2016; 99: 3916–3924. https://doi.org/10.3168/jds.2015-10484 PMID: 26971144

14. Wicik Z, Gajewska M, Majewska A, Walkiewicz D, Osińska E, Motyl T. Characterization of microRNA profile in mammary tissue of dairy and beef breed heifers. Journal of Animal Breeding and Genetics. 2016; 133: 31–42. https://doi.org/10.1111/jbg.12172 PMID: 26060050

15. Billa PA, Faulconnier Y, Ye T, Chervet M, Le Provost F, Pires JAA, et al. Deep RNA-Seq reveals miR-Nome differences in mammary tissue of lactating Holstein and Montbéliarde cows. BMC Genomics. 2019; 20: 621. https://doi.org/10.1186/s12864-019-5987-4 PMID: 31362707

16. Le Guillou S, Leduc A, Laubier J, Barbey S, Rossignol M-N, Lefebvre R, et al. Characterization of Holstein and Normande whole milk miRNomes highlights breed specificities. Sci Rep. 2019; 9: 1–13. https://doi.org/10.1038/s41598-018-37186-2 PMID: 30626917

17. INRA. Alimentation des bovins, ovins et caprins: besoins des animaux, valeurs des aliments. Editions Quae Versailles, France; 2007.

18. Billa PA, Faulconnier Y, Larsen T, Leroux C, Pires J a. A. Milk metabolites as noninvasive indicators of nutritional status of mid-lactation Holstein and Montbéliarde cows. Journal of Dairy Science. 2020;0. https://doi.org/10.3168/jds.2019-17466

19. Bernard L, Richard C, Gelin V, Leroux C, Heyman Y. Milk fatty acid composition and mammary lipo-genic genes expression in bovine cloned and control cattle. Livestock Science. 2015; 176: 188–195. https://doi.org/10.1016/j.livsci.2015.03.016

20. Chen C-J, Servant N, Toedling J, Sarazin A, Marchais A, Duvernois-Berthet E, et al. ncPRO-seq: a tool for annotation and profiling of ncRNAs in sRNA-seq data. Bioinformatics. 2012; 28: 3147–3149. https://doi.org/10.1093/bioinformatics/bts587 PMID: 23044543

21. 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: 37–52. https://doi.org/10.1093/nar/gkr688 PMID: 21911355

22. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society: Series B (Methodological). 1995; 57: 289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x

23. Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 2013; 14: 128. https://doi.org/10.1186/1471-2105-14-128 PMID: 23586463

24. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016; 44: W90–97. https://doi.org/10.1093/nar/gkw377 PMID: 27141961

25. Supek F, Bošnjak M, Škunca N, Šmuc T. REVIGO Summarizes and Visualizes Long Lists of Gene Ontology Terms. PLOS ONE. 2011; 6: e21800. https://doi.org/10.1371/journal.pone.0021800 PMID: 21789182

26. Gross J, van Dorland HA, Bruckmaier RM, Schwarz FJ. Performance and metabolic profile of dairy cows during a lactational and deliberately induced negative energy balance with subsequent realignment. J Dairy Sci. 2011; 94: 1820–1830. https://doi.org/10.3168/jds.2010-3707 PMID: 21426971

27. Pires JAA, Stumpf LF, Soutullo ID, Pescara JB, Stocks SE, Grummer RR. Effects of abomasal infusion of nicotinic acid on responses to glucose and β-agonist challenges in underfed lactating cows. J Dairy Sci. 2016; 99: 2297–2307. https://doi.org/10.3168/jds.2015-10308 PMID: 26774717

28. Pomiès D, Martin B, Chilliard Y, Pradel P, Rémond B. Once-a-day milking of Holstein and Montbéliarde cows for 7 weeks in mid-lactation. Animal. 2007;1. https://doi.org/10.1017/S1751731107658467 PMID: 22444203

29. Ferlay A, Martin B, Lerch S, Gobert M, Pradel P, Chilliard Y. Effects of supplementation of maize silage diets with extruded linseed, vitamin E and plant extracts rich in polyphenols, and morning v. evening milking on milk fatty acid profiles in Holstein and Montbéliarde cows. Animal. 2010; 4: 627–640. https://doi.org/10.1017/S1751731109991224 PMID: 22444050

30. Ji Z, Liu Z, Chao T, Hou L, Fan R, He R, et al. Screening of miRNA profiles and construction of regulation networks in early and late lactation of dairy goat mammary glands. Sci Rep. 2017; 7: 1–13. https://doi.org/10.1038/s41598-016-0028-x PMID: 28127051
31. Govignon-Gion A, Fritz S, Larroque H, Brochard M, Chantry C, Lahalle F, et al. QTL detection for milk fatty acids in French dairy cattle. Proceedings, 10th World Congress of Genetics Applied to Livestock Production 2014; 10 World Congress of Genetics Applied to Livestock Production, Vancouver, CAN, 2014-08-17-2014-08-22, 1-3. 2014 [cited 18 Jan 2020]. Available: http://agris.fao.org/agris-search/search.do?recordID=LV2016006646

32. Williams A, Henao-Mejia J, Harman CCD, Flavell RA. miR-181 and Metabolic Regulation in the Immune System. Cold Spring Harb Symp Quant Biol. 2013; 78: 223-230. https://doi.org/10.1101/sqb.2013.78. 020024 PMID: 24163395

33. Gu Z, Eleswarapu S, Jiang H. Identification and characterization of microRNAs from the bovine adipose tissue and mammary gland. FEBS Letters. 2007; 581: 981-988. https://doi.org/10.1016/j.febslet.2007.01.081 PMID: 17306260

34. Fatima A, Lynn DJ, O'Boyle P, Seoighe C, Morris D. The miRNAome of the postpartum dairy cow liver in negative energy balance. BMC Genomics. 2014; 15: 279. https://doi.org/10.1186/1471-2164-15-279 PMID: 24725334

35. Al-Husseini W, Chen Y, Gondro C, Herd RM, Gibson JP, Arthur PF. Characterization and Profiling of Liver microRNAs by RNA-sequencing in Cattle Divergently Selected for Residual Feed Intake. Asian-Australas J Anim Sci. 2016; 29: 1371-1382. https://doi.org/10.5713/ajas.15.0605 PMID: 26954124

36. Bell AW. Regulation of organic nutrient metabolism during transition from late pregnancy to early lactation. Journal of Animal Science. 1995; 73: 2804. https://doi.org/10.2527/1995.7392804x PMID: 8582872

37. Zhang L, Huang C, Guo Y, Gou X, Hinsdale M, Lloyd P, et al. MicroRNA-26b Modulates the NF-κB Pathway in Alveolar Macrophages by Regulating PTEN. The Journal of Immunology. 2015; 195: 5404-5414. https://doi.org/10.4049/jimmunol.1402933 PMID: 26503952

38. Sordillo LM, Mavangira V. The nexus between nutrient metabolism, oxidative stress and inflammation in transition cows. Anim Prod Sci. 2014; 54: 1204. https://doi.org/10.1071/AN14503

39. Du X, Yang Y, Xu C, Peng Z, Zhang M, Lei L, et al. Upregulation of miR-181a impairs hepatic glucose and lipid homeostasis. Oncotarget. 2017; 8: 91362-91378. https://doi.org/10.18632/oncotarget.20523 PMID: 29207650

40. Ma L, Qiu H, Chen Z, Li L, Zeng Y, Luo J, et al. miR-25 modulates triacylglycerol and lipid accumulation in goat mammary epithelial cells by repressing PGC-1beta. Journal of Animal Science and Biotechnology. 2018; 9: 48. https://doi.org/10.1186/s40104-018-0262-0 PMID: 29946461

41. Radisky DC. miR-200c at the nexus of epithelial-mesenchymal transition, resistance to apoptosis, and the breast cancer stem cell phenotype. Breast Cancer Research. 2011; 13: 110. https://doi.org/10.1186/bcr2885 PMID: 21682933

42. Hilmarsdottir B, Briem E, Bergthorsson JT, Magnusson MK, Gudjonsson T. Functional Role of the microRNA-200 Family in Breast Morphogenesis and Neoplasia. Genes. 2014; 5: 804-820. https://doi.org/10.3390/genes5030804 PMID: 25216122

43. Dessauge F, Lollivier V, Ponchon B, Bruckmaier R, Finot L, Wiart S, et al. Effects of nutrient restriction on mammary cell turnover and mammary gland remodeling in lactating dairy cows. Journal of Dairy Science. 2011; 94: 4623-4635. https://doi.org/10.3168/jds.2010-4012 PMID: 21854935

44. Sengar GS, Deb R, Singh U, Raja TV, Kant R, Sajjanar B, et al. Differential expression of microRNAs associated with thermal stress in Frieswal (Bos taurus x Bos indicus) crossbred dairy cattle. Cell Stress Chaperones. 2018; 23: 155-170. https://doi.org/10.1007/s12192-017-0833-6 PMID: 28776223

45. Zhao C, Tian F, Yu Y, Liu G, Zan L, Updike M, et al. miRNA-dysregulation associated with tenderness variation induced by acute stress in Angus cattle. J Anim Sci Biotechnol. 2012; 3: 12. https://doi.org/10.1186/2049-1891-3-12 PMID: 22958451

46. Zheng Y, Chen K, Zheng X, Li H, Wang G. Identification and bioinformatics analysis of microRNAs associated with stress and immune response in serum of heat-stressed and normal Holstein cows. Cell Stress and Chaperones. 2014; 19: 973-981. https://doi.org/10.1007/s12192-014-0521-8 PMID: 24917036

47. Tili E, Michaille J-J, Cimino A, Costinean S, Dumitru CD, Adair B, et al. Modulation of miR-155 and miR-125b Levels following Lipopolysaccharide/TNF-α Stimulation and Their Possible Roles in Regulating the Response to Endotoxin Shock. J Immunol. 2007; 179: 5082-5089. https://doi.org/10.4049/jimmunol.179.8.5082 PMID: 17911593

48. Abdelatty AM, Iwaniuk ME, Garcia M, Moyes KM, Teter BB, Delmonte P, et al. Effect of short-term feed restriction on temporal changes in milk components and mammary lipogenic gene expression in mid-lactation Holstein dairy cows. J Dairy Sci. 2017; 100: 4000-4013. https://doi.org/10.3168/jds.2016-11130 PMID: 28237595
49. Bionaz M, Loor JJ. ACSL1, AGPAT6, FABP3, LPIN1, and SLC27A6 Are the Most Abundant Isoforms in Bovine Mammary Tissue and Their Expression Is Affected by Stage of Lactation. J Nutr. 2008; 138: 1019-1024. https://doi.org/10.1093/jn/138.6.1019 PMID: 18492828

50. Skibiel AL, Zachut M, do Amaral BC, Levin Y, Dahl GE. Liver proteomic analysis of postpartum Holstein cows exposed to heat stress or cooling conditions during the dry period. Journal of Dairy Science. 2018; 101: 705-716. https://doi.org/10.3168/jds.2017-13258 PMID: 29128212

51. Akbar H, Bionaz M, Carlson DB, Rodriguez-Zas SL, Everts RE, Lewin HA, et al. Feed restriction, but not l-carnitine infusion, alters the liver transcriptome by inhibiting sterol synthesis and mitochondrial oxidative phosphorylation and increasing gluconeogenesis in mid-lactation dairy cows. Journal of dairy science. 2013; 96(4):2201-13. https://doi.org/10.3168/jds.2012-6036 PMID: 23462174

52. Bionaz M. Loor JJ. Gene networks driving bovine milk fat synthesis during the lactation cycle. BMC Genomics. 2008; 9: 366. https://doi.org/10.1186/1471-2164-9-366 PMID: 18671863

53. Wigle JT, Oliver G. Prox1 Function Is Required for the Development of the Murine Lymphatic System. Cell. 1999; 98: 769-778. https://doi.org/10.1016/s0092-8674(00)81511-1 PMID: 10499794

54. Betterman KL, Paquet-Fifield S, Asselin-Labat M-L, Visvader JE, Butler LM, Stacker SA, et al. Remodeling of the Lymphatic Vasculature during Mouse Mammary Gland Morphogenesis Is Mediated via Epi-theial-Derived Lymphangiogenic Stimuli. The American Journal of Pathology. 2012; 181: 2225-2238. https://doi.org/10.1016/j.ajpath.2012.08.035 PMID: 23063660

55. Sanchez M-P, Ramayo-Caldas Y, Wolf V, Laithier C, El Jabri M, Michenet A, et al. Sequence-based GWAS, network and pathway analyses reveal genes co-associated with milk cheese-making properties and milk composition in Montbéliarde cows. Genetics Selection Evolution. 2019; 51: 34. https://doi.org/10.1186/s12711-019-0473-7 PMID: 31262251

56. Désert C, Duclos MJ, Blavy P, Lecerf F, Moreews F, Klopp C, et al. Transcriptome profiling of the feeding-to-fasting transition in chicken liver. BMC Genomics. 2008; 9: 611. https://doi.org/10.1186/1471-2164-9-611 PMID: 19091074

57. Cunningham HC, Cammack KM, Hales KE, Freetly HC, Lindholm-Perry AK. Differential transcript abundance in adipose tissue of mature beef cows during feed restriction and realimentation. PLOS ONE. 2018; 13: e0194104. https://doi.org/10.1371/journal.pone.0194104 PMID: 29534082

58. Miller AM, Gilchrist DS, Nijjar J, Araldi E, Ramirez CM, Lavery CA, et al. MiR-155 Has a Protective Role in the Development of Non-Alcoholic Hepatosteatosis in Mice. PLOS ONE. 2013; 8: e72324. https://doi.org/10.1371/journal.pone.0072324 PMID: 23991091