# Food deprivation affects the miRNome in the lactating goat mammary gland

Lenha Mobuchon, Sylvain S. Marthey, Sandrine Le Guillou, Denis Laloë, Fabienne Le Provost, Christine Leroux

► To cite this version:

Lenha Mobuchon, Sylvain S. Marthey, Sandrine Le Guillou, Denis Laloë, Fabienne Le Provost, et al. Food deprivation affects the miRNome in the lactating goat mammary gland. PLoS ONE, 2015, 10 (10), 20 p. 10.1371/journal.pone.0140111. [DOI: 10.1371/journal.pone.0140111](https://doi.org/10.1371/journal.pone.0140111)

HAL Id: hal-02630964

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

Submitted on 27 May 2020

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>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Food Deprivation Affects the miRNome in the Lactating Goat Mammary Gland

Lenha Mobuchon¹,²,³, Sylvain Marthey¹, Sandrine Le Guillou¹, Denis Laloë¹, Fabienne Le Provost¹, Christine Leroux²,³*

1 INRA, UMR1313 Génétique Animale et Biologie Intégrative, F-78350, Jouy-en-Josas, France, 2 INRA, UMR1213 Herbivores, F-63122, Saint Genès Champanelle, France, 3 Clermont Université, VetAgro Sup, UMR1213 Herbivores, BP 10448, F-63000, Clermont-Ferrand, France

* Christine.Leroux@clermont.inra.fr

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

CrossMark
click for updates

## Abstract

### Background

Nutrition affects milk composition thus influencing its nutritional properties. Nutrition also modifies the expression of mammary genes, whose regulation is not fully understood. MicroRNAs (miRNA) are small non coding RNA which are important post-transcriptional regulators of gene expression by targeting messenger RNAs. Our goal was to characterize miRNA whose expression is regulated by nutrition in the lactating goat mammary gland, which may provide clues to deciphering regulations of the biosynthesis and secretion of milk components.

### Methodology/principal findings

Using high-throughput sequencing technology, miRNomes of the lactating mammary gland were established from lactating goats fed ad libitum or deprived of food for 48h affecting milk production and composition. High throughput miRNA sequencing revealed 30 miRNA with an expression potentially modulated by food deprivation; 16 were down-regulated and 14 were up-regulated. Diana-microT predictive tools suggested a potential role for several nutriregulated miRNA in lipid metabolism. Among the putative targets, 19 were previously identified as differently expressed genes (DEG). The functions of these 19 DEG revealed, notably, their involvement in tissue remodelling.

### Conclusion/significance

In conclusion, this study offers the first evidence of nutriregulated miRNA in the ruminant mammary gland. Characterization of these 30 miRNA could contribute to a clearer understanding of gene regulation in the mammary gland in response to nutrition.

## Introduction

The mammary gland synthesizes and secretes large quantities of milk components during lactation; these include proteins, lactose and lipids which have a considerable effect on the

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

OPEN ACCESS

**Citation:** Mobuchon L, Marthey S, Le Guillou S, Laloë D, Le Provost F, Leroux C (2015) Food Deprivation Affects the miRNome in the Lactating Goat Mammary Gland. PLoS ONE 10(10): e0140111. doi:10.1371/journal.pone.0140111

**Editor:** Christian Schönbach, Nazarbayev University, KAZAKHSTAN

**Received:** March 17, 2015

**Accepted:** September 22, 2015

**Published:** October 16, 2015

**Copyright:** © 2015 Mobuchon 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:** RNA sequencing data were deposited in the Gene Expression Omnibus (GEO): GSE61025.

**Funding:** This study was supported by the ApisGene (France) and INRA. Lenha Mobuchon's grant was financially supported by INRA and ApisGene. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing Interests:** The authors have declared that no competing interests exist.
nutritional, technological and sensory properties of milk. Ruminant milk production and composition are linked to both intrinsic and extrinsic factors, such as nutrition [1]. The synthesis and secretion of milk components by mammary epithelial cells involve numerous genes. Interestingly, nutrition affects the expression of genes encoding important factors for milk production in ruminants [2–4]. More precisely, 48 hours of food deprivation applied to lactating goats was shown to cause a drop in milk production and component secretion, associated with the altered expression of 161 genes, including those coding for lipogenic enzymes and major milk proteins [2]. However, the mechanisms underlying the regulation by nutrition of these genes remain poorly documented.

MicroRNAs (miRNA) are small non-coding RNA of ~22 nucleotides in length that principally regulate gene expression by altering mRNA stability or translation. MiRNA are thought to regulate at least ~60% of genes and are involved in all cellular processes [5–8]. Increasing numbers of miRNA are now being identified thanks to the latest advances in deep sequencing which enables access to an overview of the miRNA expressed in a tissue (referred to a miRNome).

However, few studies have reported the role of miRNA in development of the mammary gland [9–13] and in the lactation function [14, 15]. In particular, miRNA have been characterized in the ruminant mammary gland [16–19]. Recent *in vitro* studies have notably revealed potential roles for *miR-103* and *miR-27a* in the regulation of milk fat synthesis in goat mammary epithelial cells [20, 21].

Moreover, miRNA expression and function have been reported to be modulated by diets involving a deficiency or augmented intake [22, 23]. For example, caloric restriction up-regulates the expression of *miR-140-3p* in the epiphyseal growth plate of rats, determined using microarray and RT-qPCR analyses [24]. Ørom *et al*. [25] demonstrated that caloric restriction for six months had a considerable effect on the miRNA profile of the mouse mammary gland, determined by microarray. In addition, the expression of about fifteen known and predicted miRNA was changed in the muscle of monkeys subjected to caloric restriction when compared to those receiving an *ad libitum* diet [26]. Nevertheless, the nutritional regulation of ruminant miRNome is still poorly documented. In particular, in the ruminant mammary gland, no studies offering evidence of the effects of nutrition on miRNA expression have yet been performed.

Therefore, in order to better understand the mechanisms underlying the regulation of mammary gene expression in response to dietary factors, the aim of this study was to obtain an overview of mammary miRNA of goat which received two different diets and to identify miRNA differently expressed in these two conditions. Moreover, we studied the relationship between differently expressed miRNA and genes in the same samples. To our knowledge, this is the first study to have addressed the impacts of nutrition on the expression of miRNA in the mammary gland of ruminants and on lactation.

# Materials and Methods

## Ethics Statement

This study was performed at the experimental unit of INRA Research Center of Theix in 2003. At that time we did not require to submit each animal trail to ethical committee, the institute had its recommendations of the Animal Care as well as the Ethics Committee for Animal Experimentation of our region (Auvergne: CEMEAA number 02) that have been strictly followed. To minimize suffering, all goats were euthanized using captive bolt gun followed by exsanguination at the abattoir of the research center under accreditation number 63 345 001. All samples have been collected once the animals have been slaughtered and not while they were alive.
## Animals and RNA preparation

Twelve Alpine goats at peak lactation (48±2 days post-partum at the beginning of the experiment) from the Lusignan experimental station (France) were chosen based on their homogeneity of their milk yield, number of lactations and their genotype at the α-S1-casein locus which has an impact on milk protein content and indirectly on milk fat content. Goats were fed twice daily an orchard grass hay-based diet for two weeks during the pre-experimental period. Forty-eight hours before slaughtering, six goats were fed *ad libitum* (Control), whereas the six others were food deprived (FD). Milk was collected 15 min before euthanization (S1 Fig). Tissues samples were collected within 30 min after euthanization. All animal manipulations and samples collection were the same as described by Ollier *et al.* [2]. Mammary gland tissues were isolated just after slaughtering, snap frozen immediately and stored at -80°C. RNA were extracted using TRIZol® Reagent (Life Technologies) and further purified using the SV Total RNA isolation system (Promega) to eliminate contaminating genomic DNA. The concentration and quality of RNA preparations were assessed using spectrophotometry (Nanodrop™, ND-1000) and 2100 Bioanalyzer Instrument (Agilent). Samples with an RNA Integrity Number (RIN) higher than eight were used. For high throughput sequencing, RNA from four Control goats and six FD goats were used.

## Sequencing library preparation and data processing

The preparation and sequencing of libraries were performed by the IGBMC Microarray and Sequencing Platform (Strasbourg, France). Briefly, 7 µg of total RNA were used to generate small RNA libraries using the TruSeq™ Small RNA kit protocol (Illumina). In the first step, RNA adapters were ligated sequentially to each end of the RNA: first the 3' RNA adapter (5' TGGAATTCTCGGGTGCCAAGG 3') that was specifically modified to target small RNA, then the 5' RNA adapter (5' GTTCAGAGTTCTACAGTCCGACGATC 3'). Small RNA ligated with 3' and 5' adapters were reverse transcribed and PCR amplified (30 sec at 98°C; [10 sec at 98°C, 30 sec at 60°C, 15 sec at 72°C] x 13 cycles; 10 min at 72°C) to obtain cDNA constructs. Finally, acrylamide gel purification of the 140-150 nt amplified cDNA constructs (corresponding to cDNA inserts from small RNA + 120 nt from the adapters) was performed. The libraries were checked for quality and then quantified using a 2100 Bioanalyzer Instrument (Agilent). The libraries were loaded in the flowcell at a concentration of 8 pM and clusters were generated using the Cbot before being sequenced on HiSeq 2500 (Illumina) as single-end 50 base reads, according to Illumina's instructions.

For analysis, sequences of adapters were removed using Cutadapt [27] and reads were filtered according to their size (17–28 nt). Data analyses were processed using mostly miRDeep2 software [28] as described by Le Guillou *et al.* [19]. The cleaned sequences were clustered into unique reads and then mapped to the reference Goat genome (CHIR_1.0, http://goat.kiz.ac.cn/GGD/, [29]). Precursors and miRNA were identified using the miRDeep2 core module, miR-Deep2.pl. Potential miRNA datasets were created by adding known miRNA in goat (miRBase v21) to miRNA associated with predicted precursors with a miRDeep2 score ≥0. This number scored for the likelihood for a miRNA to be a real precursor according to the combined of energetic stability, positions and frequencies of read with Dicer processing. This score is comprised between 10 and -10, so a score ≥0 has been chosen representing because it allowed the discovery of a maximum of miRNA minimizing the rate of false positives. Then, the same operation was performed to create a data set of potential precursors. The quantifier.pl miRDeep2 module was then used to map unique reads, the set of potential miRNA and all known miRNA (miR-Base v21) on the set of potential precursors enabling the annotation of miRNA. The quantification results produced by the quantifier.pl module were then filtered with a custom perl script
parse_mIRDeep2 outputs.pl (https://mulcyber.toulouse.inra.fr/projects/bioinfutils/) to eliminate any redundancy between known and predicted miRNA. RNA sequencing data were deposited in the Gene Expression Omnibus (GEO): GSE61025.

## Quantitative RT-PCR

Seven miRNA were chosen for RT-qPCR validation among known miRNA which expression was affected by the food deprivation and based on their level of expression; miR-99a-5p (TaqMan® ID 006254 Mat, Applied Biosystems), miR-126-3p (TaqMan® ID 008451 Mat), miR-140-3p (TaqMan® ID 471823 Mat), miR-222-3p (TaqMan® ID 000525), miR-223-3p (TaqMan® ID 002295), miR-204-5p (TaqMan® ID 000508), miR-409-3p (TaqMan® ID 002332). Reverse transcription was achieved on 10 ng of total RNA using the TaqMan® MicroRNA Reverse Transcription kit (Applied Biosystems, Foster City, CA, USA), according to the manufacturer's instructions. In the thermal cycler (StepOne+, Applied Biosystems, Foster City, CA, USA), each 15 µL RT reaction followed 30 min at 16°C, 30 min at 42°C and 5 min at 85°C. Then, 1.3 µL of miRNA-specific cDNA from the reaction were amplified using TaqMan® Small RNA Assays (Applied Biosystems, Foster City, CA, USA) according to the manufacturer's instructions. Amplification was performed at 95°C for 10 min, followed by 40 cycles of 95°C for 15 sec and 60°C for 1 min. All miRNA levels were normalized to the values of U6 snoRNA [30, 31] and the results expressed as fold changes of threshold cycle (Ct) values relative to the control using the $2^{-\Delta\Delta Ct}$ method [32].

## Statistical analysis

For sequencing data, a comprehensive description of miRNA expression patterns was first of all made using a specific principal component analysis (PCA). Each goat is represented by a point, each colour corresponds to a diet and goats which have received the same diet are connected to their barycentre and an inertia ellipse.

Statistical analysis to determine differential expression was performed with R version 3.0.2 (R Development Core Team, 2013, http://www.R-project.org/) using the Bioconductor package DESeq2 version 1.0.17 [33]. DESeq2 utilizes a negative binomial distribution to model read counts per miRNA and then implements a method to normalise the counts. This normalisation procedure uses the library median of the ratios between the read count and the geometric mean of each miRNA as a scaling factor for each library. Fold changes were estimated using an empirical Bayes shrinkage procedure.

This procedure helps to moderate the broad spread in fold changes for genes with low counts, while it has negligible effect on genes with high counts. Since hypothesis tests are performed for miRNA-by-miRNA differential analyses, the p-values obtained need to be adjusted to correct for multiple testing. However, procedures to adjust p-values in order to control the number of false positives found often lead to a loss of power to detect truly differentially expressed miRNA because of the large number of hypothesis tests performed. To reduce the impact of such procedures, the filtering method described by Rau *et al.* [34] was used to remove genes that appeared to generate an uninformative signal. This method identifies a filtering threshold that maximizes so-called filtering similarity among replicates. Tests for differential expression were only applied to miRNA whose maximum count across all ten samples was higher than its threshold. This method was implemented under the Bioconductor HTSFilter package, version 1.0.0. The threshold value was found to be equal to 28. The p-values were adjusted for multiple testing using the Benjamini and Hochberg method [35], and those with an adjusted p-value <0.1 were considered to be significant.
For RT-qPCR analyses, data are represented as log base 2 of the fold change between Control and FD goats. The data were analysed using the Mann-Whitney test performed using R version 3.0.2 (R Development Core Team, 2013, http://www.R-project.org/).

## Analysis of targeted pathways

Putative targets for differently expressed miRNA were predicted with a high degree of accuracy based on DIANA-microT-CDS version 5.0 [36] (http://diana.imis.athena-innovation.gr/). Targets of predicted miRNA were searched using TargetScan Custom release version 5.2 [37] (http://www.targetscan.org/). Putative targeted pathways were investigated through the use of QIAGEN's Ingenuity Pathway Analysis (IPA, QIAGEN Redwood City, www.qiagen.com/ingenuity).

## Results and Discussion

### Global description of mammary miRNomes in lactating goats

Ten libraries were constructed from the mammary gland of four Control and six FD goats using the Illumina/Solexa technology. A mean of about 14.2 and 16.4 million reads were obtained from the Control and FD libraries, respectively (Table 1). After removing sequencing adapters and filtering reads by their size, 6.8 and 8.5 million reads of 17–28 nt were obtained for the Control and FD libraries, respectively. The reads were mapped to the goat genome (CHIR 1.0), and clustered in 67,971 and 75,474 unique sequences, on average, in the Control and FD libraries, respectively. Using the goat precursors reported in miRBase version 21 (v21), 1,070 precursors were first of all predicted by miRDeep2 for Control and FD goats. Then, using all matures from miRBase v21, it was possible to identify 1,804 miRNA. This step allowed the detection of 539 known miRNA, corresponding to 422 miRNA known in goat and 117 miRNA not yet described in the goat but known in other species. In addition, 1,265 predicted miRNA were characterized, corresponding to miRNA never described in any species.

By performing expression pattern analysis using PCA, the ten libraries appeared to be clearly discriminated according to the treatments applied to the goats (Fig 1).

Using HTSFilter [34], a statistical threshold was applied and reduced the 1,804 miRNA detected to 470, 343 were known and 127 predicted (S1 Table). At this step, lists of miRNA were conserved between the different nutritional conditions.

<table><caption>**Table 1. Summary of sequencing data.**</caption><thead><tr><th></th><th>Control<sup>a</sup></th><th>FD<sup>b</sup></th></tr></thead><tbody><tr><td>Raw reads</td><td>14,214,146</td><td>16,415,279</td></tr><tr><td>Cleaned reads<sup>c</sup></td><td>13,787,000</td><td>15,681,069</td></tr><tr><td>Sized reads<sup>d</sup></td><td>6,838,797</td><td>8,464,633</td></tr><tr><td>Sized and unique sequences processed</td><td>215,663</td><td>204,005</td></tr><tr><td>Reads mapped<sup>e</sup></td><td>5,784,949</td><td>7,278,339</td></tr><tr><td>Unique sequences corresponding to mapped reads</td><td>67,971</td><td>75,474</td></tr></tbody></table>

doi:10.1371/journal.pone.0140111.t001

<footnote><sup>a</sup>means of data for the 4 Control goat libraries</footnote>

<footnote><sup>b</sup>means of data for the 6 FD goat libraries</footnote>

<footnote><sup>c</sup>sequencing adapters removed</footnote>

<footnote><sup>d</sup>17-28 nt size filter, used by the miRDeep2 software</footnote>

<footnote><sup>e</sup>reads with at least one and at most five reported alignments, used by the miRDeep2 software.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig 1.** Principal component analyses of miRNA in Control and FD goats. L0X: represent 1 library, FD: Food deprived. Analyses were performed on miRNA for libraries according to the diet received by the goats. The two different diets are labelled within their 95% inertia ellipse (Red: Control; Green: FD).</figcaption></figure>

doi:10.1371/journal.pone.0140111.g001

Identification of miRNA differentially expressed in the mammary gland of goats after food deprivation

To identify miRNA significantly affected by food deprivation, a differential analysis was performed. Using the DESeq2 package, of the 470 miRNA remaining after the HTSFilter analysis, 30 miRNA were found to be differently expressed between Control and FD goats, of which 14 were up-regulated and 16 were down-regulated (Table 2). Among the potentially nutriregulated miRNA, 19 were known miRNA and 11 were predicted (Table 2 and S2 Table; S2 Fig). Among these 30 potentially nutriregulated miRNA pointed out by high-throughput sequencing, the expression of seven miRNA have been studied by RT-qPCR analysis on 12 goats (six Control and six FD) (S3 Fig) with the aim to reinforce the sequencing data. Probably due to the difference of sensibility and normalization, only the expression of *miR-204-5p* was found significantly affected when analysed using RT-qPCR. Nevertheless, the orientation of regulation of all miRNA tested was in accordance with data obtained by sequencing data. It should be noted that nutritional intervention studies aim to measure low changes in gene expression. In these conditions, it has been suggested that PCR should not be considered as validation tool but as complementary technology [38]. The authors concluded that PCR approach allowed picking up the expression of genes that were not measurable on high throughput technology but also vice versa, and that both techniques have their own (dis)advantages and specificities. For less pronounced changes, both technologies may be useful as complementation rather than validation.
<table><caption>Table 2. MiRNA whose expression was affected by food deprivation in lactating goat mammary gland.</caption><thead><tr><th>miRNA name</th><th>Family</th><th>Sequence</th><th>Control*</th><th>FD*</th><th>Fold change (FD/control)</th><th>P-value adjusted</th></tr></thead><tbody><tr><td colspan="7">Known miRNA</td></tr><tr><td>miR-6119-5p</td><td></td><td>AGAGGTAAAAAATTGATTTGACT</td><td>19,082</td><td>14,396</td><td>0.77</td><td>0.07</td></tr><tr><td>miR-126-3p</td><td>mir-126</td><td>TCGTACCGTGAGTAATAATGC</td><td>17,218</td><td>22,489</td><td>1.28</td><td>0.06</td></tr><tr><td>miR-660-5p</td><td>mir-188</td><td>TACCCATTGCATATCGGAGCTGT</td><td>7,698</td><td>6,144</td><td>0.81</td><td>0.09</td></tr><tr><td>miR-99a-5p</td><td>mir-10</td><td>AACCCGTAGATCCGATCTTGT</td><td>5,533</td><td>7,709</td><td>1.35</td><td>0.06</td></tr><tr><td>miR-451-5p</td><td>mir-451</td><td>AAACCGTTACCATTACTGA</td><td>4,957</td><td>6,855</td><td>1.33</td><td>0.1</td></tr><tr><td>let-7c-5p</td><td>let-7</td><td>TGAGGTAGTAGGTTGTATGGTT</td><td>4,934</td><td>6,427</td><td>1.28</td><td>0.06</td></tr><tr><td>miR-140-3p</td><td>mir-140</td><td>ACCACAGGGTAGAACCACGGAC</td><td>3,798</td><td>4,888</td><td>1.27</td><td>0.09</td></tr><tr><td>miR-409-3p</td><td>mir-154</td><td>GAATGTTGCTCGGTGAACCCC</td><td>1,954</td><td>1,159</td><td>0.66</td><td>0.06</td></tr><tr><td>miR-125b-3p</td><td>mir-10</td><td>ACAAGTCAGGCTCTTGGGACC</td><td>1,918</td><td>2,496</td><td>1.27</td><td>0.06</td></tr><tr><td>miR-222-3p</td><td>mir-221</td><td>AGCTACATCTGGCTACTGGGTCTC</td><td>812</td><td>1,408</td><td>1.47</td><td>0.09</td></tr><tr><td>miR-204-5p</td><td>mir-204</td><td>TTCCCTTTGTCATCCTATGCCT</td><td>541</td><td>1,113</td><td>1.8</td><td>0.00</td></tr><tr><td>miR-99a-3p</td><td>mir-10</td><td>CAAGCTCGCTTCTATGGGTCTGT</td><td>472</td><td>615</td><td>1.27</td><td>0.08</td></tr><tr><td>miR-196a-5p</td><td>mir-196</td><td>TAGGTAGTTTCATGTTGTTGG</td><td>413</td><td>663</td><td>1.47</td><td>0.06</td></tr><tr><td>miR-494-3p</td><td>mir-154</td><td>TGAAACATACACGGGAAACCTCT</td><td>311</td><td>197</td><td>0.67</td><td>0.05</td></tr><tr><td>miR-223-3p</td><td>mir-223</td><td>TGTCAGTTTGTCAAATACCCCA</td><td>245</td><td>552</td><td>1.65</td><td>0.06</td></tr><tr><td>miR-188-5p</td><td>mir-188</td><td>CATCCCTTGCATGGTGGAGGG</td><td>117</td><td>81</td><td>0.73</td><td>0.09</td></tr><tr><td>miR-671-5p</td><td>mir-671</td><td>AGGAAGCCCTGGAGGGGCTGGAGG</td><td>73</td><td>50</td><td>0.73</td><td>0.09</td></tr><tr><td>miR-541-5p</td><td>mir-541</td><td>AAAGGATCTGCTGTCGGTCCCACT</td><td>29</td><td>15</td><td>0.64</td><td>0.08</td></tr><tr><td>miR-223-5p</td><td>mir-223</td><td>TGTGTATTTGACAAGCTGAGTTG</td><td>17</td><td>34</td><td>1.6</td><td>0.06</td></tr><tr><td colspan="7">Predicted miRNA</td></tr><tr><td>chr12_17655-3p</td><td></td><td>CCCGGGTTTCGGCACCA</td><td>611</td><td>337</td><td>0.66</td><td>0.09</td></tr><tr><td>chr19_26739-5p</td><td></td><td>GAGGGTTTGGGTTTGGTCGTGGGA</td><td>508</td><td>289</td><td>0.67</td><td>0.09</td></tr><tr><td>chr3_4386-5p</td><td></td><td>ATAGTTCATTCAGGTTTTTCTG</td><td>420</td><td>227</td><td>0.66</td><td>0.09</td></tr><tr><td>chr27_34538-5p</td><td></td><td>AAAGTTCATTCAGGTTTTTCTG</td><td>412</td><td>224</td><td>0.66</td><td>0.09</td></tr><tr><td>chr3_4386-3p</td><td></td><td>GAAGAATCTGAATGAACATTT</td><td>253</td><td>152</td><td>0.66</td><td>0.05</td></tr><tr><td>chr23_30758-5p</td><td></td><td>GGGGATGTAGCTCAGTGGTAGA</td><td>133</td><td>68</td><td>0.66</td><td>0.09</td></tr><tr><td>chr9_13534-5p</td><td></td><td>GTACATGATGACAACTGGCTC</td><td>83</td><td>113</td><td>1.33</td><td>0.06</td></tr><tr><td>chr12_18027-5p</td><td></td><td>GAAAGGTTCATTTGGGTTTTT</td><td>56</td><td>97</td><td>1.51</td><td>0.08</td></tr><tr><td>chr22_29775-3p</td><td></td><td>ATCAGAGTCACGGCACCA</td><td>37</td><td>16</td><td>0.59</td><td>0.05</td></tr><tr><td>chr3_3319-5p</td><td></td><td>GAAAGTTTGTTTGGGTTTTTC</td><td>39</td><td>13</td><td>0.63</td><td>0.08</td></tr><tr><td>chr4_6064-5p</td><td></td><td>CCCGATCTTGTCTGATC</td><td>29</td><td>10</td><td>0.63</td><td>0.08</td></tr></tbody></table>

*Mean of normalized read counts.

The seed region of each miRNA is underlined. A positive log2 fold change means the corresponding miRNA is more highly expressed in FD goats than in Controls and vice versa. P-values were adjusted using the Benjamini-Hochberg [35] correction at 10%.

doi:10.1371/journal.pone.0140111.t002

Looking at the abundance of the 30 nutriregulated miRNA, two known miRNA (miR-126-3p and miR-6119-5p) were very highly expressed with normalized read counts greater than 10,000, while seven (let-7c-5p, miR-99a-5p, miR-125b-3p, miR-140-3p, miR-409-3p, miR-451-5p and miR-660-5p) were highly expressed and represented by normalized read counts of between 1,000 and 10,000 (Table 2). In addition, seven known (miR-99a-3p, miR-188-5p, miR-196a-5p, miR-204-5p, miR-222-3p, miR-223-3p and miR-494-3p) and six predicted miRNA (chr3_4386-3p, chr3_4386-5p, chr12_17655-3p, chr19_26739-5p, chr23_30758-5p and chr27_34538-5p) were expressed at moderate levels with normalized read counts of between
100 and 1,000, and three known (*miR-223-5p*, *miR-541-5p* and *miR-671-5p*) and five predicted (*chr3_3319-5p*, *chr4_6064-5p*, *chr9_13534-5p*, *chr12_18027-5p* and *chr23_30758-5p*) were weakly expressed with normalized read counts below 100. This last category may have a lesser impact in cells because of their weak expression.

Interestingly, the expression of both matures of the same precursor (*mir-99*, *mir-223* and *chr3_4384*) were affected by nutrition. This case has only ever been reported once for *miR-9/9*\*, whose expression was dramatically decreased after maternal supplementation with retinoic acid in the rat embryo spina bifida model [39]. The regulation may occur on the promoter of miRNA gene leading to deregulate both -5p and -3p strands. In all other cases in the literature, only one strand of the precursor has been found to be affected by nutrition. These observations suggest that two different regulatory mechanisms may exist. In this case, the regulation may happen on the stability of one or the other strand leading to the deregulation of only one strand by nutrition which seems to be the main mechanism in our study. These two different mechanisms of regulation require further investigations.

Otherwise, two clusters, defined as miRNA genes with an inter-distance ≤10 kb [40], were altered by food deprivation. Indeed, *miR-494-3p*, *miR-409-3p* and *miR-541-5p* are members of the greatest and conserved miRNA cluster located on goat chromosome 21. These three miRNA were down-regulated by food deprivation and their expressions were well correlated with a high mean Pearson correlation coefficient of 0.8. Similarly, *let-7c-5p* and *mir-99*, which clustered on the goat's X chromosome, were up-regulated in FD goats and displayed a strong mean Pearson correlation coefficient of 0.8. These results are consistent with the previously reported co-regulation of clustered miRNA, where they were transcribed as a long polycistron [41].

Among the nutriregulated miRNA, five (*miR-126-3p*, *miR-140-3p*, *miR-223-3p*, *miR-409-3p*, *miR-451-5p*) have previously been described as being affected by caloric restriction in others species and tissues. The expressions of *miR-223-3p* and *miR-409-3p* were increased, whereas that of *miR-451-5p* was decreased in skeletal muscle of monkeys receiving a restricted diet [26]. Conversely, *miR-409-3p* was down-regulated in the blood of mice receiving about 60% of their normal diet, whereas the expression of *miR-451-5p* was increased [42]. Most notable was the fact that the expressions of *miR-409-3p* (down-regulated in FD) and *miR-451-5p* (up-regulated in FD) appeared to be sensitive to caloric restriction. This apparent discrepancy between findings may have been due to the type of restriction, and probably to the tissues and species studied. In addition, the expressions of *miR-126-3p* and *miR-140-3p*, which were up-regulated in our study, were down-regulated in the epiphyseal growth plate of Sprague-Dawley rats by 60% food restriction for 10 days [24]. Although this nutritional model appears to be the closest to ours, deprivation in our case was much more drastic because the animals received no nutrients at all. Finally, among the predicted miRNA, the *chr3_4386-5p* and *chr27_34538-5p* sequences were found to be almost similar except at position 2, suggesting that these miRNA may belong to the same miRNA family. In addition some nutriregulated miRNA belonged to the same family. Indeed, *miR-409-3p* and *miR-494-3p* belonged to the *mir-154* family and *miR-188-5p* and *miR-660-5p* to the *mir-188* family (Table 2). Interestingly, it is known that members of a miRNA family are designed to execute similar biological functions [43]. Therefore, nutriregulated miRNA sharing the same family may have linked roles in the lactating mammary gland.

## Putative functions of nutriregulated miRNA

Among the 30 miRNA that were regulated in our model, only a few have previously been described in the mammary gland. In the mouse mammary gland, *miR-126-3p* may inhibit cell proliferation; it may also regulate the expression of CSN2 coding for β-casein, one of the major milk proteins, and PGR (ProGesteron Receptor) protein [44]. In addition, *miR-99a-5p* has
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Molecular and cellular functions of potentially targeted genes by the 30 nutriregulated miRNA. Data were analysed through the use of QIAGEN's Ingenuity Pathway Analysis (IPA®), QIAGEN Redwood City, www.qiagen.com/ingenuity).</figcaption></figure>

doi:10.1371/journal.pone.0140111.g002

been described as a modulator of the TGF β pathway affecting epithelial-to-mesenchyme transition in normal mouse mammary gland [45]. Therefore, in order to investigate the functional role of the 30 nutriregulated miRNA in the mammary gland, target gene predictions were performed based on miRNA/mRNA interactions using Diana-microT v5.0 for known miRNA and TargetScan Custom v5.2 for predicted miRNA. Because no predictive tools have so far been designed for ruminants, conservation between goat miRNA and humans was checked, particularly concerning the seed region that determines miRNA interactions with mRNA targets. Thus, 7,129 different putative targets for the 30 nutriregulated miRNA were identified. Investigating the molecular and cellular functions of these targets using Ingenuity Pathway Analysis revealed that “gene expression”, “cellular development” and “cellular growth and proliferation” were the three pathways most significantly targeted (Fig 2). These findings suggest that food deprivation might lead to a change in gene expression through the actions of miRNA linked to cellular growth and proliferation, as well as the remodelling of mammary cells, as has been observed in this organ in the event of nutrient restriction [46].

In addition, some putative targets of nutriregulated miRNA in this study were associated with involution in the mammary gland; these included factors promoting involution such as *CASP3*, *IL6* and *STAT3*, or factors inhibiting involution such as *BCL2* [47, 48] (S3 Table). Markers of autophagosome formation, *MAP1LC3B2* and *ATG7*, were also found to be potentially targeted by *miR-188-5p* and *miR-223-3p*, respectively. These observations were in line with the drop in milk production and component synthesis previously observed in FD goats [2].

## Potential involvement of nutriregulated miRNA in the mammary biosynthesis and secretion of milk components

**Potential involvement in mammary lipid metabolism.** Of all the 7,129 putative targets of nutriregulated miRNA, particular efforts were made to identify those involved in lipid metabolism. Very interestingly, nutriregulated miRNA may act in fatty acid availability (*VLDLR*), synthesis (ACACA), transport (FABP), activation (ACS) and desaturation (SCD), triglyceride synthesis (AGPAT) and in genes coding for proteins associated with the milk fat globule membrane (BTN, XDH) (Fig 3) [49]. Several nutriregulated miRNA also target transcription factors
such as the SREBF1 and PPARγ cofactors, *PPARGC1A* and *PPARGC1B* that govern the expres-
sion of essential lipogenic enzymes in the mammary gland (Fig 3). Moreover, other genes
related to milk lipid metabolism, such as *INSIG2* (INSulin_Induced Genes 2) which codes for
an endoplasmic reticulum protein that blocks the processing of SREBP, is predicted to be tar-
geted by *chr12_18027-5p* [50]. Elsewhere, *miR-204-5p* may also target *LPIN2* (LiPIN 2),
recently characterized in the bovine mammary gland [49] and involved in triacylglycerol accu-
mulation [51]. The gene coding for FADS1 (Fatty Acid DeSaturase 1), which plays a role in
triacylglycerol synthesis [52], is potentially targeted by *miR-671-5p*. Furthermore, *miR-125b-
3p*, *miR-494-3p* and *chr3_3319-5p* may conjointly target ABCA1 (ATP-Binding Cassette sub-
family A member 1) which has been hypothesized to play a role in the transfer, storage and
removal of cholesterol in the mammary gland [53]. Finally, five nutriregulated miRNA: *miR-
222-3p*, *miR-188-5p*, *miR-541-5p*, *miR-494-3p* and *chr3_3319-5p*, may target *PTEN* (Phospha-
tase and TENsin homologue) (Fig 4). In particular, *miR-494-3p* binding has been confirmed
experimentally in bronchial epithelial cells [54]. The role of *PTEN* has been characterized in
dairy cow mammary epithelial cells (DCMECs), revealing that it plays a crucial role in the via-
bility and proliferation capacity of cells, as well as in the secretion of β casein, triglycerides and
lactose. In particular, *PTEN* may regulate the expression of key lactation-related pathways such
as *PPARγ*, *SREBF1*, *mTOR*, *PRLR* and *GLUT1*, suggesting a pleiotropic role for *PTEN* in the
lactation process [55].

Among the predicted miRNA, *chr3_3319-5p* appeared to be of great interest. Indeed, *chr3_3319-5p* potentially targets numerous genes involved in mammary metabolism. Its seed region is similar to that of *hsa-miR-561-3p*. This miRNA has not yet been reported in either the goat or cattle, but its identical seed region with *hsa-miR-561* suggests that this miRNA is part of the *mir-561* family, even though the function of this family is not yet known in the mammary gland (Fig 4, S4 and S5 Tables).

Taken together, these observations highlight a potentially synergistic action of miRNA in regulating mammary metabolic pathways, and offer perspectives for the study of new types of regulation in the lactating mammary gland in response to dietary treatments. Those miRNA which expression was affected by nutrition may modulate that of the genes responsible for fatty acid metabolism, lipid droplet formation and/or milk fat globule secretion at the same time, and thus consequently have a crucial role to play in the synthesis and secretion of milk components. In particular, *miR-409-3p* and *miR-494-3p* may target numerous genes involved in lipid metabolism and its regulation. It was notable that they belong to the same cluster and the same family, *mir-154*, which could be of considerable interest in future studies.

**Potential involvement in the metabolism of mammary protein and lactose.** Proteins are also a principal constituent of milk and are synthesized by mammary epithelial cells during lactation. Amino acids are supplied to the mammary epithelial cells by a system of amino acid membrane transporters [56]. Among the genes putatively targeted by nutriregulated miRNA, several amino acid transporters were identified (Table 3). The availability of amino acids in the mammary gland is important for the regulation of translation, and also appears to be one of the major limitations affecting milk protein synthesis [57], and nutriregulated miRNA may regulate this availability. Caseins (αs1, αs2, β andκ), the major milk proteins, are phosphorylated in mammary epithelial cells. This phosphorylation is involved in aggregation of the four caseins in micelles before their secretion into milk. Despite this key role, the phosphorylation step in mammary epithelial cells is still poorly documented. Interestingly, *miR-222-3p*, *miR-409-3p*, *miR-541-5p*, *miR-6119-5p*, *miR-660-5p*, *chr23_30758-5p* and *chr3_4386-5p* potentially target several genes coding for casein kinases 1 and 2 alpha and gamma (S4 Table). In addition, *miR-204-5p* potentially targets *ELF5* (E74-like factor 5) which belongs to the JAK-STAT signaling pathway and is essential for the expression of milk-related genes, particularly those coding
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Potential involvement of nutriregulated miRNA in milk fat synthesis. In red, miRNA which expression was up-regulated in FD goat mammary glands compared with Controls, while those whose expression was down-regulated are in blue.</figcaption></figure>

doi:10.1371/journal.pone.0140111.g003
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Potential contribution of nutriredulated miRNA to regulating PTEN and associated networks. GLUT1: GLUcose Transporter 1, mTOR: mechanistic Target Of Rapamycin, PTEN: Phosphatase and TENsin homolog, SREBF1: Sterol REsponse Binding Factor 1.</figcaption></figure>

doi:10.1371/journal.pone.0140111.g004

<table><caption>**Table 3.** Amino acid transporters potentially targeted by nutriredulated miRNA.</caption><thead><tr><th>miRNA</th><th>Gene</th><th>Protein</th><th>Associated transport system</th><th>Description</th><th>References</th></tr></thead><tbody><tr><td>miR-223-5p, miR-541-5p, miR-671-5p</td><td>SLC1A1</td><td>GLT-1;EAAT1</td><td>X<sup>+</sup><sub>AG</sub></td><td>Na+-dependent system for anionic amino acids</td><td>[61, 62]</td></tr><tr><td>miR-223-5p</td><td>SLC1A3</td><td>GLAST; EAAT1</td><td>X<sup>+</sup><sub>AG</sub></td><td>Na+-dependent system for anionic amino acids</td><td>[61]</td></tr><tr><td>miR-6119-5p</td><td>SLC1A4</td><td>ASCT1-SATT</td><td>ASC</td><td>Na+-dependent transporter, particularly linear dipolar amino acids (L-alanine, L-serine, L-cysteine)</td><td>[56]</td></tr><tr><td>miR-409-3p, miR-541-5p, chr3_3319-5p</td><td>SLC6A6</td><td>Taut</td><td>System Gly</td><td>Na+-Cl<sup>-</sup>-dependent system specific to β-amino acids</td><td>[61]</td></tr><tr><td>miR-671-5p</td><td>SLC7A1</td><td>CAT-1</td><td>y+</td><td>Na+-independent system specific to cationic amino acids</td><td>[63]</td></tr><tr><td>miR-126-3p</td><td>SLC7A5</td><td>LAT1</td><td>L</td><td>Na+-dependent electroneutral transport mechanism for neutral amino acids</td><td>[59]</td></tr><tr><td>miR-196a-5p, miR-223-3p</td><td>SLC7A8</td><td>LAT2</td><td>L</td><td>Na+-dependent electroneutral transport mechanism for neutral amino acids</td><td>[62]</td></tr><tr><td>miR-409-3p, miR-494-3p</td><td>SLC15A2</td><td>PEPT2</td><td>H+-peptide cotransporter 2</td><td>Proton-peptide electrogenic transporter</td><td>[64]</td></tr><tr><td>miR-671-5p</td><td>SLC36A1</td><td>PAT1</td><td></td><td>Proton/amino acid symporter</td><td>[59]</td></tr><tr><td>let-7c-5p, miR-409-3p, miR-494-3p</td><td>SLC38A2</td><td>SNAT2</td><td>System A</td><td>Na+-dependent system for amino neutral acids</td><td>[62]</td></tr></tbody></table>

Putative targets were predicted from DIANA microT v5.0 [36] for known miRNA and Custom Target Scan v5.2 [37] for predicted miRNA.

doi:10.1371/journal.pone.0140111.t003
for major milk proteins [58, 59]. Above all, *miR-99a-5p* was found to potentially target *mTOR* (mechanistic Target Of Rapamycin) and its binding site has been confirmed experimentally using the Luciferase Reporter Assay in breast cancer cell lines [60]. The regulation of protein synthesis in all mammalian tissues is under the control of *mTOR*, and especially the translation step. In addition, *mTOR* may also have a function in milk protein synthesis [59] (S4 Table), suggesting a potential role of *miR-99a-5p* in controlling this process.

Regarding differentially expressed miRNA targeting the SLC (Solute Carrier) family (S4 Table), *miR-223-5p* was found to target *SLC2A1* (GLUT1). Glucose is the primary precursor for lactose synthesis, the main milk carbohydrate. In the mammary gland of lactating cows, *GLUT1* mRNA and protein were shown to be strongly expressed whereas they were barely detectable in dry cows [65]. In addition, GLUT1 protein represents about half of glucose transporters in the Golgi membranes of lactating rat epithelial cells [66]. Consequently, *miR-223-5p* may have a role in regulating glucose uptake and lactose synthesis, a hypothesis which requires further investigation.

## Putative targets differently expressed in FD goats

Of all the putative targets predicted for nutriregulated miRNA, 47 were among the genes previously identified as being significantly affected (named DEG) by food deprivation using a microarray technique [2] (S5 Table). One to six miRNA were predicted to targets these 47 DEG. In addition, when the way of variations in DEG and nutriregulated miRNA was taken account, 19 DEG were found to be potentially targeted by nutriregulated miRNA, including 18 down-regulated DEG which might be targeted by one to three up-regulated miRNA (Table 4) one up-regulated DEG (CD24) might be targeted by *chr3_3319-5p* which was down-regulated in FD goats. By investigating the cellular and molecular functions of these 19 DEG, “cellular death and survival”, “cellular morphology” and “cell-to-cell interaction” were found to be among the five most targeted functions (data not shown). This suggests that the deregulation of miRNA may explain the phenotype of mammary tissue remodelling previously hypothesized in FD goats [2]. Among the DEG that may be targeted by three nutriregulated miRNA, *ESR1* is a transcription factor that regulates genetic programming of cell cycle progression and growth in the mammary gland. The knockout of *ESR1* in mice suggests its crucial role in alveologenesis during lactation [67]. Interestingly, *ESR1* was down-regulated by food deprivation and potentially targeted by *miR-125b-3p*, *miR-222-3p*, and *chr19_26739-5p*, and up-regulated in the mammary gland of FD goats. In addition *miR-222-3p* has been confirmed as directly targeting *ESR1* in breast cancer cell lines [68]. Its binding site (Fig 5) is conserved among different species, including cattle and goats. These observations suggest that the up-regulation of *miR-222-3p* by food deprivation may be responsible for the down-regulation of *ESR1*.

Furthermore, 11 nutriregulated miRNA (Table 4) were found to potentially target DEG in which three predicted miRNA *chr3_3319-5p*, *chr9_13534-5p* and *chr12_18027-5p*. Of the 11 nutriregulated miRNA that could target DEG, *let-7c-5p*, *miR-223-3p* and *miR-204-5p* may target seven, five and four DEG, respectively (Table 4). Notably, *miR-204-5p* and *miR-223-3p* are nutriregulated miRNA with the highest fold change in high throughput sequencing. In addition, *miR-223-3p* may target five DEG, three of which (*ANKFY1*, *PLAGL2* and *ZCCHC14*) encode for zinc finger proteins that bind DNA, suggesting that this miRNA may indirectly regulate gene transcription. Consequently, the deregulation of *miR-204-5p* and *miR-223-3p* by food deprivation may have a significant impact in the mammary gland through the regulation of DEG.

## Co-localization of nutriregulated miRNA with QTL associated with traits

Quantitative Trait Loci (QTL) have been associated with economically important traits such as health through the detection of clinical mastitis (somatic cell score) as well as milk production
<table><caption>Table 4. Differently expressed genes potentially targeted by nutriregulated miRNA.</caption><thead><tr><td>miRNA</td><td>Differently Expressed Genes</td></tr></thead><tbody><tr><td>let-7c-5p</td><td>ABCC5, ANKFY1, DCBLD1, DCUN1D3, DPAGT1, LEPREL2, PLAGL2</td></tr><tr><td>miR-99a-3p</td><td>ZCCHC14</td></tr><tr><td>miR-125b-3p</td><td>ESR1</td></tr><tr><td>miR-140-3p</td><td>PLAGL2</td></tr><tr><td>miR-204-5p</td><td>IRF2BP2, KCNH8, PHLDB1, TSPAN31</td></tr><tr><td>miR-222-3p</td><td>ESR1, INSIG1</td></tr><tr><td>miR-223-3p</td><td>ANKFY1, ARFIP1, MYO1B, PLAGL2, ZCCHC14</td></tr><tr><td>miR-223-5p</td><td>EFHC1</td></tr><tr><td>chr3_3319-5p</td><td>CD24</td></tr><tr><td>chr9_13534-5p</td><td>ESR1, ZCCHC14</td></tr><tr><td>chr12_18027-5p</td><td>RAD9A, TMEM2</td></tr></tbody></table>

In bold, miRNA or DEG whose expression are upregulated by the food deprivation. ABCC5: ATP-binding cassette, sub-family C (CFTR/MRP), member 5, ANKFY1: Ankyrin repeat and FYVE domain-containing protein 1, ARFIP1: Arfaptin-1, CD24: Signal transducer CD24 precursor, DCBLD1: Discoidin, CUB and LCCL domain-containing protein 1 precursor, DCUN1D3: DCN1, defective in cullin neddylation 1, domain containing 3, DPAGT1: UDP-N-acetylglucosamine-dolichyl-phosphate N-acetylglucosaminephosphotransferase, EFHC1: EF-hand domain-containing protein 1, ESR1: Oestrogen receptor alpha, INSIG1: Insulin-induced gene 1 protein, IRF2BP2: Interferon regulatory factor 2 binding protein 2, KCNH8: Potassium voltage-gated channel subfamily H member 8, LEPREL2: Prolyl 3-hydroxylase 3 precursor, MYO1B: Myosin Ib, PHLDB1: Pleckstrin homology-like domain family B member 1, PLAGL2: Zinc finger protein PLAGL2, RAD9A: Cell cycle checkpoint control protein Rad9A, TMEM2: Transmembrane protein 2, TSPAN31: Tetraspanin-31, ZCCHC14: Zinc finger CCHC domain-containing protein 14.

doi:10.1371/journal.pone.0140111.t004

and composition (lactose, protein and fat) in livestock, particularly cattle [69]. However no QTL are yet available in the goat. Using conservation between species, cow miRNA expressed during lactation and conserved in the goat have been located in cattle QTL associated with milk production and composition (personal communication). Among the 30 miRNA nutriregulated in the present study, *mir-671* was located in QTL linked to somatic cells on bovine chromosome 4 (Table 5). It is also located in a QTL region associated with both milk protein yield and milk fat percentage and content. In addition, *mir-140* is located in QTL associated with the somatic cell score and milk yield on bovine chromosome 18. Furthermore, *mir-409*, *mir-494* and *mir-541* belonging to a large and conserved cluster in mammals have also been located in QTL associated with milk yield and milk fat percentage and content on bovine chromosome 21. As well as having many putative targets involved in lipid metabolism (Fig 3, S4 Table), *mirR-494-3p* and *mirR-409-3p* are also positioned in QTL regions associated with milk fatty acids, thus rendering these two miRNA good candidates to study their involvement in the regulation of milk lipid synthesis. The co-localization of nutriregulated miRNA with QTL associated with milk composition may provide clues to their role in the regulation of milk traits.

## Conclusion

The findings reported here highlight for the first time using high throughput miRNA sequencing, 30 known or predicted miRNA in the ruminant mammary gland which expression seems to be nutriregulated. Among these 30 potentially interesting miRNA, members of clusters and similar families were found. Five of them have already been reported as being affected by
<figure><img src="image_8.png" /><figcaption>**Fig 5.** Nutriregulated miRNA binding sites on *ESR1* mRNA. Analyses were performed using DIANA microT v5.0 [36] for miR-125b-3p and miR-222-3p and Custom Target Scan v5.2 [37] for chr19_26739-5p. The interaction between 3'UTR of ESR1 and the seed region of miR-222-3p of one binding site is described in the box.</figcaption></figure>

doi:10.1371/journal.pone.0140111.g005

nutrition in other species and tissues. By investigating the genes they might target, some nutr下一 regulated miRNA appear to be good candidates for further studies in mammary gland biology. For example, *miR-494-3p*, located in QTL associated with milk fat in cows, may be preponderant for milk lipid synthesis, notably through the regulation of *PTEN*, and in targeting *mTOR*, *miR-99a-5p* may regulate milk protein synthesis. The links between miRNA and DEG previously identified by transcriptomic assays show that 19 DEG may be targeted by nutriregulated miRNA. These observations highlighted two miRNA, *miR-204-5p* and *miR-223-3p*, whose expression was the most markedly affected by food deprivation and which may target several DEG, suggesting important roles for both miRNA in the nutritional regulation of gene expression in the mammary gland. Nevertheless, further studies are necessary to support our results and to decipher the function of these miRNA in mammary epithelial cells.

<table><caption>**Table 5.** Nutriregulated miRNA which bovine equivalent is located in QTL associated with milk production and composition.</caption><thead><tr><th rowspan="2">QTL</th><th colspan="2">Associated with milk components</th><th rowspan="2">Associated with mammary tissue health</th></tr><tr><th>Fat percentage and content</th><th>Protein yield</th></tr></thead><tbody><tr><td>Bovine chromosome</td><td></td><td></td><td></td></tr><tr><td>BTA 4</td><td>mir-671</td><td>mir-671</td><td>mir-671</td></tr><tr><td>BTA 18</td><td></td><td></td><td>mir-140</td></tr><tr><td>BTA 21</td><td>mir-494, mir-541, mir-409a</td><td></td><td></td></tr></tbody></table>

Conservation between nutriregulated and bovine miRNA was checked and their localization in QTL (CattleQTLdb, http://www.animalgenome.org/cgi-bin/QTLdb/BT/) associated with milk was investigated.

doi:10.1371/journal.pone.0140111.t005
MiRNA offer new insights to understand the impact of nutrition on the mammary gland and may help to clarify mammary gene regulation in order to modulate milk quality via nutrition. The nutrition model of food deprivation used in this study was a large challenge, so further studies should be performed with a diet that could be applied by livestock farmers in order to identify the nutriregulation of miRNA under breeding conditions.

**Supporting Information**

**S1 Fig. Milk production, fat, protein and lactose in Control and FD goats.** Data were extracted from [2] and were analysed with a Mann-Whitney test (n = 6). **: p<0.01.d: day. (TIF)

**S2 Fig. Structure of predicted precursors whose miRNA are differently expressed in FD versus Control goats.** Putative targets were predicted from DIANA microT v5.0 [36] for known miRNA and Custom Target Scan v5.2 [37] for predicted miRNA. (TIF)

**S3 Fig. Quantitative RT-PCR validation of NGS data.** The expression of miRNA was normalized with the expression of U6 in Control and FD goats. The relative abundance of the miARN is represented with boxplots (n = 6). (TIF)

**S1 Table. miRNA expressed in lactating goat mammary glands.** (XLSX)

**S2 Table. Characteristics of nutriregulated predicted miRNA.** (XLSX)

**S3 Table. Example of putative targets involved in mammary involution.** (XLSX)

**S4 Table. Selected putative targets involved in mammary metabolism.** (XLSX)

**S5 Table. Differently expressed genes potentially targeted by nutriregulated miRNA.** (XLSX)

**Acknowledgments**

Authors are grateful to L. Bernard for tissues isolation at slaughterhouse. They also specially thank J. Domagalski and J. Laubier for technical assistance for RNA extraction and RT-qPCR advices, respectively, as well as A. De La Foye for her help in statistics and Y. Chilliard and J. Rouel for their advices in nutrition.

**Author Contributions**

Conceived and designed the experiments: CL FLP. Performed the experiments: LM. Analyzed the data: LM SM SLG DL. Contributed reagents/materials/analysis tools: LM FLP CL DL SLG SM. Wrote the paper: LM FLP CL.

**References**

1. Chilliard Y, Ferlay A. Dietary lipids and forages interactions on cow and goat milk fatty acid composition and sensory properties. Reprod Nutr Dev. 2004; 44(5):467-92. PMID: 15636165
2. Ollier S, Robert-Granie 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(3):560-7. PMID: 17311940

3. Bernard L, Leroux C, Chilliard Y. Expression and nutritional regulation of lipogenic genes in the rumi- nant lactating mammary gland. Adv Exp Med Biol. 2008;(606:):67-108.

4. Bionaz M, Loor JJ. Ruminant metabolic systems biology: reconstruction and integration of transcrip- tume dynamics underlying functional responses of tissues to nutrition and physiological state. Gene Regul Syst Bio. 2012; 6:109-25. doi: 10.4137/GRSB.S9852 PMID: 22807626

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

6. Ambros V. The functions of animal microRNAs. Nature. 2004; 431(7006):350-5. PMID: 15372042

7. Bushati N, Cohen SM. microRNA functions. Annu Rev Cell Dev Biol. 2007; 23:175-205. PMID: 17506695

8. Friedman RC, Farh KK, Burge CB, Bartel DP. Most mammalian mRNAs are conserved targets of micro-RNAs. Genome Res. 2009; 19(1):92-105. doi: 10.1101/gr.082701.108 PMID: 18955434

9. Ibarra I, Erlich Y, Muthuswamy SK, Sachidanandam R, Hannon GJ. A role for microRNAs in mainte- nance of mouse mammary epithelial progenitor cells. Genes Dev. 2007; 21(24):3238-43. PMID: 18079172

10. Greene SB, Gunaratne PH, Hammond SM, Rosen JM. A putative role for microRNA-205 in mammary epithelial cell progenitors. J Cell Sci. 2010; 123(Pt 4):606-18. doi: 10.1242/jcs.056812 PMID: 20103531

11. 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. doi: 10.1371/journal.pone.0045727 PMID: 23029204

12. Llobet-Navas D, Rodriguez-Barrueco R, Castro V, Ugalde AP, Sumazin P, Jacob-Sendler D, et al. The miR-424(322)/503 cluster orchestrates remodeling of the epithelium in the involuting mammary gland. Genes Dev. 2014; 28(7):765-82. doi: 10.1101/gad.237404.114 PMID: 24636986

13. Yoo KH, Kang K, Feuermann Y, Jang SJ, Robinson GW, Hennighausen L. The STAT5-regulated miR- 193b locus restrains mammary stem and progenitor cell activity and alveolar differentiation. Dev Biol. 2014; 395(2):245-54. doi: 10.1016/j.ydbio.2014.09.012 PMID: 25236432

14. Li HM, Wang CM, Li QZ, Gao XJ. Mir-15a decreases bovine mammary epithelial cell viability and lacta- tion and regulates growth receptor expression. Molecules. 2012; 17(10):12037-48. doi: 10.3390/ molecules171012037 PMID: 23085654

15. Wang J, Bian Y, Wang Z, Li D, Wang C, Li Q, et al. MicroRNA-152 Regulates DNA Methyltransferase 1 and Is Involved in the Development and Lactation of Mammary Glands in Dairy Cows. PLoS One. 2014; 9(7):e101358. doi: 10.1371/journal.pone.0101358 PMID: 24987964

16. Gu Z, Eleswarapu S, Jiang H. Identification and characterization of microRNAs from the bovine adipose tissue and mammary gland. FEBS letters. 2007; 581(5):981-8. PMID: 17306260

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

18. 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. PLoS One. 2012; 7(12): e52388. doi: 10.1371/journal.pone.0052388 PMID: 23300659

19. Le Guillou S, Marthey S, Laloe D, Laubier J, Mobuchon L, Leroux C, et al. Characterisation and com- parison of lactating mouse and bovine mammary gland miRNomes. PLoS One. 2014; 9(3):e91938. doi: 10.1371/journal.pone.0091938 PMID: 24658750

20. 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. doi: 10.1371/journal.pone.0079258 PMID: 24244462

21. Lin XZ, Luo J, Zhang LP, Wang W, Shi HB, Zhu JJ. mir-27a suppresses triglycerides accumulation and affects gene mRNA expression associated with fat metabolism in dairy goat mammary gland epithelial cells. Gene 2013; 521(1):15-23. doi: 10.1016/j.gene.2013.03.050 PMID: 23537996

22. Shah MS, Davidson LA, Chapkin RS. Mechanistic insights into the role of microRNAs in cancer: influ- ence of nutrient crosstalk. Frontiers in genetics. 2012; 3:1-14.

23. Garcia-Segura L, Perez-Andrade M, Miranda-Rios J. The Emerging Role of MicroRNAs in the Regula- tion of Gene Expression by Nutrients. J Nutrigenet Nutrigenomics. 2013; 6(1):16-31. doi: 10.1159/ 000345826 PMID: 23445777
24. Pando R, Even-Zohar N, Shtaif B, Edry L, Shomron N, Phillip M, et al. MicroRNAs in the growth plate are responsive to nutritional cues: association between miR-140 and SIRT1. J Nutr Biochem. 2012; 23 (11):1474-81. doi: 10.1016/j.jnutbio.2011.09.010 PMID: 22402365

25. Orom UA, Lim MK, Savage JE, Jin L, S A.D., Lisanti MP, et al. MicroRNA-203 regulates caveolin-1 in breast tissue during caloric restriction. Cell Cycle. 2012; 11(7):1291-5. doi: 10.4161/cc.19704 PMID: 22421148

26. Mercken EM, Majounie E, Ding J, Guo R, Kim J, Bernier M, et al. Age-associated miRNA alterations in skeletal muscle from rhesus monkeys reversed by caloric restriction. Aging. 2013; 5(9):692-703. PMID: 24036467

27. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnetjour- nal. 2011; 17(1):10-2.

28. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012; 40(1):37-52. doi: 10.1093/nar/gkr688 PMID: 21911355

29. Dong Y, Xie M, Jiang Y, Xiao N, Du X, Zhang W, et al. Sequencing and automated whole-genome opti- cal mapping of the genome of a domestic goat (Capra hircus). Nat Biotechnol. 2013; 31(2):135-41. doi: 10.1038/nbt.2478 PMID: 23263233

30. Masotti A, Caputo V, Da Sacco L, Pizzuti A, Dallapiccola B, Bottazzo GF. Quantification of small non-coding RNAs allows an accurate comparison of miRNA expression profiles. J Biomed Biotechnol. 2009; 659028(10):1.

31. Li D, Liu H, Li Y, Yang M, Qu C, Zhang Y, et al. Identification of suitable endogenous control genes for quantitative RT-PCR analysis of miRNA in bovine solid tissues. Mol Biol Rep. 2014; 41(10):6475-80. doi: 10.1007/s11033-014-3530-x PMID: 24973888

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

33. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010; 11 (10):2010-11.

34. Rau A, Gallopin M, Celeux G, Jaffrezic F. Data-based filtering for replicated high-throughput transcrip- tome sequencing experiments. Bioinformatics. 2013; 29(17):2146-52. doi: 10.1093/bioinformatics/ btt350 PMID: 23821648

35. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a pratical and powerful approach to mul- tiple testing. J R Statist Soc B. 1995; 57(1):289-300.

36. Paraskevopoulou MD, Georgakilas G, Kostoulas N, Vlachos IS, Vergoulis T, Reczko M, et al. DIANA- microT web server v5.0: service integration into miRNA functional analysis workflows. Nucleic Acids Res. 2013; 41(Web Server issue):W169-73. doi: 10.1093/nar/gkt393 PMID: 23680784

37. Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets: Cell. 2005; 120(1):15-20. PMID: 15652477

38. Gaj S, Eijssen L, Mensink RP, Evelo CT. Validating nutrient-related gene expression changes from microarrays using RT(2) PCR-arrays. Genes Nutr. 2008; 3(3-4):153-7. doi: 10.1007/s12263-008- 0094-1 PMID: 19034552

39. Zhao JJ, Sun DG, Wang J, Liu SR, Zhang CY, Zhu MX, et al. Retinoic acid downregulates microRNAs to induce abnormal development of spinal cord in spina bifida rat model. Childs Nerv Syst. 2008; 24 (4):485-92. PMID: 17962954

40. Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequenc- ing data. Nucleic Acids Res. 2014; 42 part 1:1-10. doi: 10.1093/nar/gkt1181 PMID: 24275495

41. Baskerville S, Bartel DP. Microarray profiling of microRNAs reveals frequent coexpression with neigh- boring miRNAs and host genes. RNA. 2005; 11(3):241-7. PMID: 15701730

42. Dhahbi JM, Spindler SR, Atamna H, Yamakawa A, Guerrero N, Boffelli D, et al. Deep sequencing iden- tifies circulating mouse miRNAs that are functionally implicated in manifestations of aging and respon- sive to calorie restriction. Aging. 2013; 5(2):130-41. PMID: 23470454

43. Ding J, Zhou S, Guan J. miRFam: an effective automatic miRNA classification method based on n- grams and a multiclass SVM. BMC Bioinformatics. 2011; 12(216):1471-2105.

44. Cui W, Li Q, Feng L, Ding W. MiR-126-3p regulates progesterone receptors and involves development and lactation of mouse mammary gland. Mol Cell Biochem. 2011; 355(1-2):17-25. doi: 10.1007/ s11010-011-0834-1 PMID: 21526342

45. Turcatel G, Rubin N, El-Hashash A, Warburton D. MIR-99a and MIR-99b modulate TGF-beta induced epithelial to mesenchymal plasticity in normal murine mammary gland cells. PLoS One. 2012; 7(1): e31032. doi: 10.1371/journal.pone.0031032 PMID: 22299047
46. 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. J Dairy Sci. 2011; 94(9):4623-35. doi: 10.3168/jds.2010-4012 PMID: 21854935

47. Zarzynska J, Motyl T. Apoptosis and autophagy in involuting bovine mammary gland. J Physiol Pharmacol. 2008; 59 Suppl 9:275-88. PMID: 19261986

48. Watson CJ, Kreuzaler PA. Remodeling mechanisms of the mammary gland during involution. Int J Dev Biol. 2011; 55(7-9):757-62. doi: 10.1387/ijdb.113414cw PMID: 22161832

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

50. Baylin A, Deka R, Tuitele J, Viali S, Weeks DE, McGarvey ST. INSIG2 variants, dietary patterns and metabolic risk in Samoa. Eur J Clin Nutr. 2013; 67(1):101-7. doi: 10.1038/ejcn.2012.124 PMID: 22968099

51. Phan J, Peterfy M, Reue K. Lipin expression preceding peroxisome proliferator-activated receptor-gamma is critical for adipogenesis in vivo and in vitro. J Biol Chem. 2004; 279(28):29558-64. PMID: 15123608

52. Han LQ, Li HJ, Wang YY, Zhu HS, Wang LF, Guo YJ, et al. mRNA abundance and expression of SLC27A, ACC, SCD, FADS, LPIN, INSIG, and PPARGC1 gene isoforms in mouse mammary glands during the lactation cycle. Genet Mol Res. 2010; 9(2):1250-7. doi: 10.4238/vol9-gmr814 PMID: 20603810

53. Mani O, Korné M, Sørensen MT, Sejrsen K, Wotzkow C, Ontsouka CE, et al. Expression, localization, and functional model of cholesterol transporters in lactating and nonlactating mammary tissues of murine, bovine, and human origin. Am J Physiol Regul Integr Comp Physiol. 2010; 299(2):R642-54. doi: 10.1152/ajpregu.00723.2009 PMID: 20445153

54. Liu L, Jiang Y, Zhang H, Greenlee AR, Han Z. Overexpressed miR-494 down-regulates PTEN gene expression in cells transformed by anti-benzo(a)pyrene-trans-7,8-dihydrodiol-9,10-epoxide. Life Sci. 2010; 86(5-6):192-8. doi: 10.1016/j.lfs.2009.12.002 PMID: 20006626

55. Wang Z, Hou X, Qu B, Wang J, Gao X, Li Q. Pten regulates development and lactation in the mammary glands of dairy cows. PLoS One. 2014; 9(7):e102118. doi: 10.1371/journal.pone.0102118 PMID: 25009983

56. Shennan DB, Boyd CA. The functional and molecular entities underlying amino acid and peptide transport by the mammary gland under different physiological and pathological conditions. J Mammary Gland Biol Neoplasia. 2014; 19(1):19-33. doi: 10.1007/s10911-013-9305-5 PMID: 24158403

57. Shennan DB, Peaker M. Transport of milk constituents by the mammary gland. Physiol Rev. 2000; 80 (3):925-51. PMID: 10893427

58. Harris J, Stanford PM, Sutherland K, Oakes SR, Naylor MJ, Robertson FG, et al. Socs2 and elf5 mediate prolactin-induced mammary gland development. Mol Endocrinol. 2006; 20(5):1177-87. PMID: 16469767

59. Bionaz M, Loor JJ. Gene networks driving bovine mammary protein synthesis during the lactation cycle. Bioinform Biol Insights. 2011; 5:83-98. doi: 10.4137/BBI.S7003 PMID: 21698073

60. Hu Y, Zhu Q, Tang L. MiR-99a antitumor activity in human breast cancer cells through targeting of mTOR expression. PLoS One. 2014; 9(3):e92099. doi: 10.1371/journal.pone.0092099 PMID: 24637915

61. Kansal VK, Sharma R, Rehan G. Characterization of anionic amino acid transport systems in mouse mammary gland. Indian J Exp Biol. 2000; 38(11):1097-103. PMID: 11395952

62. Baik M, Etchebarne BE, Bong J, VandeHaar MJ. Gene expression profiling of liver and mammary tissues of lactating dairy cows. Asian Australas J Anim Sci. 2009; 22(6):871-84.

63. Baumrucker CR. Cationic amino acid transport by bovine mammary tissue. J Dairy Sci. 1984; 67 (11):2500-6. PMID: 6440910

64. Zhou MM, Wu YM, Liu HY, Zhao K, Liu JX. Effects of tripeptides and lactogenic hormones on oligopeptide transporter 2 in bovine mammary gland. J Anim Physiol Anim Nutr. 2011; 95(6):781-9.

65. Komatsu T, Itoh F, Kushibiki S, Hodate K. Changes in gene expression of glucose transporters in lactating and nonlactating cows. J Anim Sci. 2005; 83(3):557-64. PMID: 15705752

66. Madon RJ, Martin S, Davies A, Fawcett HA, Flint DJ, Baldwin SA. Identification and characterization of glucose transport proteins in plasma membrane- and Golgi vesicle-enriched fractions prepared from lactating rat mammary gland. Biochem J. 1990; 272(1):99-105. PMID: 2264840

67. Feng Y, Manka D, Wagner KU, Khan SA. Estrogen receptor-alpha expression in the mammary epithelium is required for ductal and alveolar morphogenesis in mice. Proc Natl Acad Sci U S A. 2007; 104 (37):14718-23. PMID: 17785410
68. Zhao JJ, Lin J, Yang H, Kong W, He L, Ma X, et al. MicroRNA-221/222 negatively regulates estrogen receptor alpha and is associated with tamoxifen resistance in breast cancer. J Biol Chem. 2008; 283 (45):31079-86. doi: 10.1074/jbc.M806041200 PMID: 18790736

69. Ashwell MS, Heyen DW, Sonstegard TS, Van Tassell CP, Da Y, VanRaden PM, et al. Detection of quantitative trait loci affecting milk production, health, and reproductive traits in Holstein cattle. J Dairy Sci. 2004; 87(2):468-75. PMID: 14762090