# Differences in the microRNAs Levels of Raw Milk from Dairy Cattle Raised under Extensive or Intensive Production Systems

Loubna Abou El Qassim, Jaime Alonso, Ke Zhao, Sandrine Le Guillou, Jorge

Diez, Fernando Vicente, Manuel Fernández-Sanjurjo, Eduardo

Iglesias-Gutiérrez, Leluo Guan, Luis Royo

► To cite this version:

Loubna Abou El Qassim, Jaime Alonso, Ke Zhao, Sandrine Le Guillou, Jorge Diez, et al.. Differences in the microRNAs Levels of Raw Milk from Dairy Cattle Raised under Extensive or Intensive Production Systems. Veterinary Sciences, 2022, 9 (12), pp.661. 10.3390/vetsci9120661. hal-03934186

HAL Id: hal-03934186

https://hal.science/hal-03934186

Submitted on 11 Jan 2023

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
Article

# Differences in the microRNAs Levels of Raw Milk from Dairy Cattle Raised under Extensive or Intensive Production Systems

Loubna Abou el qassim <fnref n="1" />,Jaime Alonso <fnref n="2" />, Ke Zhao <fnref n="3" />, Sandrine Le Guillou <fnref n="4" />, Jorge Diez <fnref n="2" />, Fernando Vicente <fnref n="1" />, Manuel Fernández-Sanjurjo <fnref n="5" />,<fnref n="6" />, Eduardo Iglesias-Gutiérrez <fnref n="5" />,<fnref n="6" />, Leluo Guan <fnref n="7" /> and Luis J. Royo <fnref n="1" />,*

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

check for
updates

**Citation:** Abou el qassim, L.; Alonso, J.; Zhao, K.; Le Guillou, S.; Diez, J.; Vicente, F.; Fernández-Sanjurjo, M.; Iglesias-Gutiérrez, E.; Guan, L.; Royo, L.J. Differences in the microRNAs Levels of Raw Milk from Dairy Cattle Raised under Extensive or Intensive Production Systems. *Vet. Sci.* **2022**, *9*, 661. https://doi.org/10.3390/vetsci9120661

Academic Editor: Xinwei Li

Received: 16 September 2022

Accepted: 21 November 2022

Published: 27 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2022 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

**Simple Summary:** Extensive animal production systems are generally considered more sustainable and beneficial for the environment and for maintaining rural populations. However, there is no defined concept of what extensive milk production is. It is assumed to be a kind of production based on pasture and forage, with animals spending part of their daytime free and with a low stocking density. In order to increase consumer confidence in this type of product, we have studied markers based on molecules naturally present in the milk that allow us to differentiate the production system in which the cows that have produced the milk have come from. In addition, we are attempting to determine whether the milk production system can have benefits for consumer health.

**Abstract:** Studying microRNA (miRNAs) in certain agri-food products is attractive because (1) they have potential as biomarkers that may allow traceability and authentication of such products; and (2) they may reveal insights into the products’ functional potential. The present study evaluated differences in miRNAs levels in fat and cellular fractions of tank milk collected from commercial farms which employ extensive or intensive dairy production systems. We first sequenced miRNAs in three milk samples from each production system, and then validated miRNAs whose levels in the cellular and fat fraction differed significantly between the two production systems. To accomplish this, we used quantitative PCR with both fractions of tank milk samples from another 20 commercial farms. Differences in miRNAs were identified in fat fractions: overall levels of miRNAs, and, specifically, the levels of *bta-mir-215*, were higher in intensive systems than in extensive systems. Bovine mRNA targets for *bta-miR-215* and their pathway analysis were performed. While the causes of these miRNAs differences remain to be elucidated, our results suggest that the type of production system could affect miRNAs levels and potential functionality of agri-food products of animal origin.

**Keywords:** milk; miRNA; dairy production systems; biomarker

1. Introduction

Consumers' growing concern about food characteristics has contributed to the creation of a new concept of quality, particularly for animal products. This concept includes traditional attributes related to nutritional value, flavor, aroma, and color, together with new indicators related to ethical aspects, such as animal welfare and environmental impact of the production system [1]. Consumers assume that products from cows raised under

<footnote n="1">¹ Servicio Regional de Investigación y Desarrollo Agroalimentario (SERIDA), Carretera AS-267, PK. 19, 33300 Villaviciosa, Spain</footnote>

<footnote n="2">² Centro de Inteligencia Artificial, Universidad de Oviedo, Campus de Viesques, 33271 Gijón, Spain</footnote>

<footnote n="3">³ Institute of Food Science, Zhejiang Academy of Agricultural Sciences, Hangzhou 310004, China</footnote>

<footnote n="4">⁴ INRAE, AgroParisTech, GABI, Université Paris-Saclay, 78350 Jouy-en-Josas, France</footnote>

<footnote n="5">⁵ Department of Functional Biology, Physiology, University of Oviedo, 33006 Oviedo, Spain</footnote>

<footnote n="6">⁶ Instituto de Investigación Sanitaria del Principado de Asturias (ISPA), 33011 Oviedo, Spain</footnote>

<footnote n="7">⁷ Department of Agricultural, Food and Nutritional Science, University of Alberta, Edmonton, AB T6G 2R3, Canada</footnote>

<footnote>* Correspondence: ljroyo@serida.org</footnote>
pasture and/or grazing are more natural and better for meeting animal welfare demands than those from the cows raised under the cereal-rich diets typical of intensive production systems. Intensive systems are considered less sustainable than production based on pastures, particularly because they have a larger ecological footprint and because they divert cereals from human consumption. These considerations highlight the need to compare whether agri-food products obtained from one production system or the other may have advantages, so that consumers can make fully informed purchasing decisions [2].

Milk is a complex secretory product and a source of nutrients for children and adults. Milk also contains numerous biomolecules, including microRNAs (miRNAs), a group of non-coding RNAs that bind to specific regions within messenger RNA to regulate gene expression post-transcriptionally [3]. In vitro models [4] have suggested that exogenous miRNAs, such as bovine milk miRNAs, may influence the health of humans due to their resistance to stomach digestion, and that they may affect cellular function. Nevertheless, the potential bioactivity of dietary miRNAs is still under study, and several aspects remain poorly understood, including the dose needed to produce biological effects [5]. The miRNA level of milk depends on genetic factors, such as the animal's breed and physiology [6,7], as well as environmental factors, such as the animal's lifestyle, pathological state, and diet [8–11]. Furthermore, the miRNA level depends on the milk fractions, milk fat, whey, and cells [12]. The miRNAs found in milk are highly resistant to acidic pH and to enzymes such as RNases [13], which implies that they are very stable and resistant to industrial treatments [13–15]. For these reasons, miRNAs have been proposed as biomarkers of quality control in dairy products, and they have been used to control for fraud in the labeling of milk powder [14].

Changes in diet have been shown to affect the expression of genes in the mammary gland [16,17]. Considering that miRNAs are essential regulators of gene expression, miRNA profiles in milk may also change with diet. However, few studies have investigated the effect of diet on miRNA profiles in bovine milk [18], such as through comparisons of the quite different diets in intensive or extensive production systems. Therefore, our objective was to evaluate differences in miRNA levels of bovine milk produced under extensive or intensive systems. Characterizing production systems as intensive or extensive is not always straightforward, given the lack of regulatory definitions and the complexity of the factors involved [19]. In the present study, we defined intensive systems as those without grazing and with high amounts of maize silage and concentrate in the diet. We defined extensive systems as those with grazing and a diet with a low amount of concentrate and no maize silage. Our work was carried out in Asturias, north of Spain, where the edapho-climatic conditions permit the coexistence of different milk production systems, varying from intensive systems, with animals fed indoors, to extensive pastured-based systems. This results in a fairly wide range of production systems covering all management possibilities [20]. Furthermore, we also explored the use of miRNA markers for the characterization and traceability of agri-food products. We used bioinformatics to predict what target genes may be regulated by miRNAs whose levels differ between the two production systems, in an effort to examine how certain miRNAs may reveal the functional properties of agri-food products.

**2. Materials and Methods**

2.1. *Sample Collection and Preparation*

In order to maximize potential differences in milk composition between intensive and extensive production systems, for sequencing, we sampled six farms at the end of one spring season (May 2016). On three farms (extensive production), animals grazed for at least for 12 h each day and ate a diet based on fresh grass with a small amount of concentrates. On three farms (intensive production), animals did not graze, instead eating a diet based on conserved feed and concentrates. Farm characteristics and details of the diets are described in Table 1. All animals on all farms were Holstein cows.
<table><caption>**Table 1.** Characteristics and diet of cows on the dairy farms where milk was sampled for RNA sequencing.</caption><thead><tr><th>Production System</th><th>Number of Cows</th><th>Milk Production (L/Day/Cow)</th><th>Grazing (h/Day)</th><th>Grass Silage (kg F */Day/Cow)</th><th>Maize Silage (kg F */Day/cow)</th><th>Hay (kg F */Day/Cow)</th><th>Concentrate (kg F */Day/Cow)</th></tr></thead><tbody><tr><td rowspan="3">Intensive (3 farms)</td><td>51</td><td>30</td><td>0</td><td>17</td><td>10.0</td><td>2</td><td>10</td></tr><tr><td>65</td><td>28</td><td>0</td><td>10</td><td>15.0</td><td>10</td><td>11</td></tr><tr><td>90</td><td>29</td><td>0</td><td>16</td><td>20.0</td><td>2</td><td>12</td></tr><tr><td rowspan="3">Extensive (3 farms)</td><td>24</td><td>21</td><td>20</td><td>10</td><td>0.0</td><td>2</td><td>7</td></tr><tr><td>14</td><td>31</td><td>&gt;12</td><td>14</td><td>0.0</td><td>6</td><td>6</td></tr><tr><td>15</td><td>29</td><td>18</td><td>15</td><td>0.0</td><td>3</td><td>6</td></tr></tbody></table>

For validation by RT-qPCR, twenty farms were sampled during autumn 2017 and spring 2018. Ten dairy farms applied an extensive system in which animals grazed and consumed a diet based on fresh grass, without maize silage and with a low amount of concentrates (4–8 kg/animal/day), while the other ten applied an intensive system in which animals were housed and ate a diet consisting of a high dry matter intake of maize silage (16–30 kg/animal/day) and a high amount of concentrate (>10 kg/animal/day) (Table 2).

<table><caption>**Table 2.** Characteristics and diet from twenty dairy farms featuring extensive or intensive production systems, where milk was sampled for quantitative real-time PCR validation.</caption><thead><tr><th>Production System</th><th>Number of Cows</th><th>Milk Production (L/Day/Cow)</th><th>Grazing (h/Day)</th><th>Grass Silage (kg F */Day/Cow)</th><th>Maize Silage (kg F */Day/Cow)</th><th>Hay (kg F */Day/Cow)</th><th>Concentrate (kg F */Day/Cow)</th></tr></thead><tbody><tr><td rowspan="10">Intensive (10 farms)</td><td>124</td><td>37.4</td><td>0</td><td>8.0</td><td>30.0</td><td>0.8</td><td>11.5</td></tr><tr><td>116</td><td>37.0</td><td>0</td><td>10.0</td><td>30.0</td><td>0.0</td><td>10.5</td></tr><tr><td>90</td><td>29.0</td><td>0</td><td>16.0</td><td>20.0</td><td>3.0</td><td>12.0</td></tr><tr><td>240</td><td>36.0</td><td>0</td><td>10.0</td><td>16.0</td><td>2.5</td><td>12.0</td></tr><tr><td>250</td><td>38.0</td><td>0</td><td>12.0</td><td>30.0</td><td>0.9</td><td>12.3</td></tr><tr><td>37</td><td>27.0</td><td>0</td><td>14.0</td><td>28.0</td><td>0.0</td><td>10.5</td></tr><tr><td>110</td><td>30.0</td><td>0</td><td>16.0</td><td>20.0</td><td>0.0</td><td>11.0</td></tr><tr><td>72</td><td>28.0</td><td>0</td><td>15.0</td><td>20.0</td><td>2.5</td><td>12.0</td></tr><tr><td>118</td><td>36.0</td><td>0</td><td>16.0</td><td>22.0</td><td>0.0</td><td>10.5</td></tr><tr><td>124</td><td>37.0</td><td>0</td><td>11.0</td><td>20.0</td><td>4.5</td><td>12.0</td></tr><tr><td rowspan="8">Extensive (10 farms)</td><td>20</td><td>21.0</td><td>6</td><td>10.0</td><td>0.0</td><td>6.8</td><td>6.5</td></tr><tr><td>24</td><td>21.0</td><td>20</td><td>12.0</td><td>0.0</td><td>0.0</td><td>7.0</td></tr><tr><td>12</td><td>26.2</td><td>22</td><td>0.0</td><td>0.0</td><td>6.0</td><td>4.0</td></tr><tr><td>8</td><td>18.8</td><td>21</td><td>0.0</td><td>0.0</td><td>4.6</td><td>4.1</td></tr><tr><td>35</td><td>19.5</td><td>20</td><td>14.0</td><td>0.0</td><td>0.0</td><td>4.7</td></tr><tr><td>15</td><td>29.0</td><td>18</td><td>33.0</td><td>0.0</td><td>0.7</td><td>6.2</td></tr><tr><td>30</td><td>27.0</td><td>18</td><td>0.0</td><td>0.0</td><td>4.0</td><td>8.0</td></tr><tr><td>7</td><td>20.0</td><td>22</td><td>6.0</td><td>0.0</td><td>0.0</td><td>5.0</td></tr><tr><td>16</td><td>23.0</td><td>21</td><td>0.0</td><td>0.0</td><td>0.0</td><td>5.0</td></tr><tr><td>22</td><td>24.0</td><td>20</td><td>0.0</td><td>0.0</td><td>0.0</td><td>6.0</td></tr></tbody></table>

<footnote>* kg F, kg of fresh matter.</footnote>

<footnote>* kg F, kg of fresh matter.</footnote>
From each type of farm, bulk tank milk was sampled after an even number of milkings (to avoid differences due to afternoon and morning milk composition). Samples were maintained at 4 °C and immediately transported to the laboratory for processing. Milk was well mixed, and 50 mL of each sample was centrifuged at 1900 × g for 20 min. The fat in the upper phase was transferred to a new 50-mL RNase-free tube. Then, 7.5 mL of QIAzol lysis reagent (Qiagen, Crawley, UK) was added, and the emulsion was vigorously mixed until the fat was well dispersed. The pellet (cellular fraction) from the initial centrifugation was washed twice with phosphate-buffered saline, then homogenized with 1 mL of QIAzol lysis reagent. All samples were stored at -80 °C until RNA extraction.

## 2.2. RNA Isolation

For each sample, total RNA was extracted from 2 mL of milk fat in QIAzol lysis reagent, and from 1 mL of cell pellet resuspended in QIAzol lysis reagent. RNA was extracted using the miRVana miRNA isolation kit (Applied Biosystems, Foster City, CA, USA) following the manufacturer's instructions, then stored at -80 °C. The concentration and integrity of RNA (RIN: RNA Integrity Number) for sequencing was further determined on an Agilent 2100 Bioanalyzer using an RNA 6000 Pico kit (both from Agilent Technologies, Santa Clara, CA, USA).

## 2.3. *Search for miRNAs Candidates from Sequencing*

### 2.3.1. RNA Sequencing

In order to identify miRNAs differing between the two production systems, 12 libraries corresponding to the fat or cellular fractions of milk from six samples of bulk tank milk (three per production system) were prepared and sequenced using the Illumina platform (Illumina, San Diego, CA, USA) as 50 bp single reads. The raw sequence data were processed for quality control, and low quality reads were removed from raw data using CASAVA 1.8 based on chastity. Then, adaptors were trimmed, and sequences with read lengths between 15 to 40 nt were mapped to the bovine genome (bostau 7) and the miRNA database (miRBase, release version 21) in order to identify the known miRNAs. Expression levels of each miRNA were estimated based on the frequency of reads, and results were normalized to the number of reads per million (RPM) using the following formula: RPM = (specific miRNA reads number/total mapped miRNA reads per library) × 10⁶.

### 2.3.2. Identification of Reference miRNAs for qPCR Normalization

In order to select the miRNAs to be used as candidates for normalizers in RT-qPCR, we chose those miRNAs with more stable expression among samples for each milk fraction, that is, miRNAs with the smallest coefficient of variation (CV = standard deviation/mean).

### 2.3.3. Identification of miRNAs Whose Levels Differed between Production Systems

In order to identify those miRNAs whose levels differ between production systems, the results from miRNA sequencing were analyzed using three statistical tests. One test was the ratio of the difference between the means to the sum of the standard deviations between the two production systems: value = |(mean1-mean2)| / (standard deviation1 + standard deviation2) (strict Cohen's d). This ratio indicates how many sums of deviations fit between the means; the higher the ratio is, the greater the difference between the means [21]. Another test was Student's t test, for which lower t values indicated greater differences between the means for the two production systems [22]. The third test was the absolute value of the correlation coefficient; the higher this value was, the greater the difference between the two systems [23]. Afterward, to choose the miRNAs whose levels differed between production systems, they were ranked first according to each statistical test, and then according to the average of the classification of the three tests.
## 2.4. Validation of Candidate miRNAs Using RT-qPCR

### 2.4.1. RT-qPCR Analysis

A subset of the miRNAs which were determined to differ significantly between the two production systems identified in Section 3.2 were validated using quantitative RT-PCR and milk samples from twenty dairy farms featuring extensive or intensive production systems (Table 2).

Total RNA was used for cDNA synthesis using the TaqMan Advanced miRNA cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA), and the resulting cDNA was stored at -20 °C until use. Levels of miRNAs were determined using quantitative RT-PCR (TaqMan Advanced miRNA Assays; ThermoFisher Scientific, Waltham, MA, USA) in a StepOne thermocycler (Applied Biosystems, Foster City, CA, USA). The final reaction solution contained 10 µL of 2 × TaqMan Fast Advanced Master mix (ThermoFisher Scientific, Waltham, MA, USA), 1 µL of 20 × TaqMan Advanced miRNA Assay (ThermoFisher Scientific, Waltham, MA, USA), 4 µL of RNase free water, and 5 µL of cDNA (diluted 1:10). The thermocycling program was set at 95 °C for 20 s, followed by 40 cycles at 95 °C for 1 s and 60 °C for 20 s. All PCR reactions were performed in duplicate, and a maximum of 0.5 threshold cycles were permitted between duplicates.

### 2.4.2. Selection of Stable Reference miRNAs. GeNorm Analysis

Normalization is an essential component of a reliable qPCR assay. geNorm [24] is one of the most popular algorithms to find stable reference genes from a set of tested candidate reference genes in a given experimental condition. We used geNorm to find the optimal number and choice of reference genes for normalization, using miRNAs identified in Section 3.2 in 22 tank milk samples, representing the experimental variation of the dairy production systems existing in the area of study [25].

### 2.4.3. miRNAs Levels Normalization and Estimation

Levels of miRNA were normalized based on the geometric mean of the selected reference miRNAs selected by geNorm, estimated using QBase+ 3.1 software, Biogazelle (Gent, Belgium) [26] and expressed in base log10. Unless otherwise noted, results were reported as mean ± standard deviation. Mean miRNAs levels between extensive and intensive production systems were compared using Student's t test in R-Commander 2.7-1. Significance was defined as *p* < 0.05.

## 2.5. Prediction and Functional Analysis of Genes Targeted by miRNAs

The Target Scan 7.2 bioinformatics tool [27] was used to predict bovine mRNA targets of candidate miRNAs. A pathway analysis of targeted genes was performed using Panther bioinformatics tool version 16.0 (http://www.pantherdb.org/ accessed on 8 June 2021) based on Gene Ontology classification [28].

# 3. Results

## 3.1. miRNAs Levels in Fat and Cellular Fractions of Milk

The total mean RNA concentration was 192.4 ± 37.4 ng/µL in milk fat and 30.6 ± 15.8 ng/µL (mean ± SD) in milk cells. The RIN value was 2.6 ± 0.15 for RNA from milk fat, and 6.5 ± 0.3 for RNA from milk cells.

The six libraries from the cellular fraction of milk yielded a mean of 24,017,290.17 reads, significantly more than the 6,964,122.33 reads from the six libraries from the fat fraction (p = 0.004, Table 3). Almost half of the reads came from small RNAs, of which the most abundant were transfer RNAs (tRNAs) in the cellular fraction and non-coding RNAs (ncRNAs) in the fat fraction. Significant differences were found between the fat and cellular fractions for a percentage of all small RNAs, except for small nuclear RNAs (snRNAs) and miRNAs. Ribosomal RNAs (rRNAs), small nucleolar RNAs (snoRNAs), and ncRNAs were more abundant in the fat fraction than in the cellular fraction, while the converse was true for tRNAs. The levels of miRNAs in the fat fraction differed significantly between extensive
and intensive systems (p = 0.040, Table 4). Intensive production was associated with higher miRNAs levels.

<table><caption>**Table 3.** Read mapping statistics.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th rowspan="2">Total Reads</th><th rowspan="2">Small RNA Reads</th><th colspan="7">Small RNAs (%)</th></tr><tr><th>Total</th><th>rRNA</th><th>snRNA</th><th>snoRNA</th><th>tRNA</th><th>miRNA</th><th>ncRNA</th></tr></thead><tbody><tr><td rowspan="6">Cellular fraction</td><td rowspan="3">Intensive Farms</td><td>32,560,190</td><td>16,693,358</td><td>0.5</td><td>12.7</td><td>3.3</td><td>5.7</td><td>45.0</td><td>2.2</td><td>31.1</td></tr><tr><td>26,678,803</td><td>13,626,331</td><td>0.5</td><td>8.3</td><td>2.0</td><td>4.7</td><td>64.3</td><td>0.9</td><td>19.8</td></tr><tr><td>23,616,528</td><td>11,871,084</td><td>0.5</td><td>16.5</td><td>3.7</td><td>5.5</td><td>37.6</td><td>2.0</td><td>34.8</td></tr><tr><td rowspan="3">Extensive Farms</td><td>19,944,125</td><td>9,753,050</td><td>0.5</td><td>13.6</td><td>3.1</td><td>5.6</td><td>51.5</td><td>1.4</td><td>24.7</td></tr><tr><td>21,166,161</td><td>9,650,109</td><td>0.5</td><td>14.9</td><td>3.9</td><td>6.6</td><td>45.4</td><td>1.4</td><td>27.8</td></tr><tr><td>20,137,934</td><td>10,053,989</td><td>0.5</td><td>15.9</td><td>3.9</td><td>6.1</td><td>36.8</td><td>2.3</td><td>34.9</td></tr><tr><td rowspan="6">Fat fraction</td><td rowspan="3">Intensive Farms</td><td>7,620,977</td><td>3,867,636</td><td>0.5</td><td>20.6</td><td>3.4</td><td>6.1</td><td>29.8</td><td>1.7</td><td>38.4</td></tr><tr><td>7,209,138</td><td>3,139,312</td><td>0.4</td><td>21.9</td><td>3.7</td><td>6.6</td><td>25.3</td><td>1.9</td><td>40.6</td></tr><tr><td>6,180,232</td><td>3,115,110</td><td>0.5</td><td>26.0</td><td>3.8</td><td>6.0</td><td>15.8</td><td>1.6</td><td>46.9</td></tr><tr><td rowspan="3">Extensive Farms</td><td>7,056,743</td><td>2,599,438</td><td>0.4</td><td>19.4</td><td>3.7</td><td>8.4</td><td>24.8</td><td>2.2</td><td>41.5</td></tr><tr><td>6,682,264</td><td>3,279,814</td><td>0.5</td><td>25.0</td><td>4.1</td><td>7.1</td><td>16.9</td><td>2.4</td><td>44.5</td></tr><tr><td>7,035,380</td><td>3,854,280</td><td>0.6</td><td>25.9</td><td>4.3</td><td>7.0</td><td>19.2</td><td>2.3</td><td>41.4</td></tr></tbody></table>

rRNA: ribosomal ribonucleic acid; snRNA: small nuclear RNA; snoRNA: small nucleolar RNA; tRNA: transfer RNA; miRNA: microRNA; ncRNA: non-coding RNA.

<table><caption>**Table 4.** Differences in abundance of small RNA classes depending on milk fraction and production system.</caption><thead><tr><th rowspan="2">Small RNA Class</th><th colspan="3">p-Value Based on Student's t Test</th></tr><tr><th>Cellular vs. Fat Fraction *</th><th>Cellular Fraction: Extensive vs. Intensive</th><th>Fat Fraction: Extensive vs. Intensive</th></tr></thead><tbody><tr><td>rRNA</td><td>0.004</td><td>0.513</td><td>0.827</td></tr><tr><td>snRNA</td><td>0.196</td><td>0.268</td><td>0.184</td></tr><tr><td>snoRNA</td><td>0.024</td><td>0.127</td><td>0.050</td></tr><tr><td>tRNA</td><td>0.004</td><td>0.827</td><td>0.513</td></tr><tr><td>miRNA</td><td>0.260</td><td>0.825</td><td>0.040</td></tr><tr><td>Non-coding RNA</td><td>0.004</td><td>0.827</td><td>0.513</td></tr></tbody></table>

We identified 518 known miRNAs in the cellular fraction of milk and 477 in the fat fraction. Most of these miRNAs (454) were present in both fractions (Table S1).

## 3.2. Validation of miRNAs Whose Levels Differed between Intensive and Extensive Production

The first 10 miRNAs in the fat fraction whose levels differed between the two production systems are ranked in Table 5, and those in the cellular fraction are ranked in Table 6. Levels of the first five miRNAs from each fraction were validated using quantitative RT-PCR and milk samples from another 20 farms. The following miRNAs in the fat fraction were subjected to validation: bta-miR-215, bta-miR-369-3p, bta-miR-6520, bta-miR-7863, and bta-miR-133a. The following miRNAs in the cellular fraction were subjected to validation: bta-miR-574, bta-miR-3432a, bta-miR-2285e, bta-miR-197, and bta-miR-2284y.

<footnote>* Both intensive and extensive production systems together. rRNA: ribosomal ribonucleic acid; snRNA: small nuclear RNA; snoRNA: small nucleolar RNA; tRNA: transfer RNA; miRNA: microRNA; ncRNA: non-coding RNA. In bold are significant differences.</footnote>
<table><caption>**Table 5. Ranking average of the first ten differentially expressed miRNAs in the fat fraction, the three applied tests' values, and their rankings.**</caption><thead><tr><th rowspan="2">miRNA</th><th colspan="3">Result for</th><th colspan="3">Ranking According to</th><th rowspan="2">Average Ranking</th></tr><tr><th>Test 1 <sup>a</sup></th><th>Test 2 <sup>b</sup></th><th>Test 3 <sup>c</sup></th><th>Test 1</th><th>Test 2</th><th>Test 3</th></tr></thead><tbody><tr><td>bta-miR-215 *</td><td>3.180</td><td>0.010</td><td>0.960</td><td>1</td><td>1</td><td>1</td><td>1.0</td></tr><tr><td>bta-miR-369-3p *</td><td>1.760</td><td>0.020</td><td>0.900</td><td>3</td><td>2</td><td>2</td><td>2.3</td></tr><tr><td>bta-miR-6520 *</td><td>1.360</td><td>0.030</td><td>0.850</td><td>4</td><td>3</td><td>4</td><td>3.7</td></tr><tr><td>bta-miR-7863 *</td><td>1.970</td><td>0.080</td><td>0.860</td><td>2</td><td>7</td><td>3</td><td>4.0</td></tr><tr><td>bta-miR-133a *</td><td>1.300</td><td>0.040</td><td>0.840</td><td>5</td><td>5</td><td>5</td><td>5.0</td></tr><tr><td>bta-miR-532</td><td>1.260</td><td>0.040</td><td>0.840</td><td>6</td><td>4</td><td>6</td><td>5.3</td></tr><tr><td>bta-miR-148a</td><td>1.210</td><td>0.120</td><td>0.780</td><td>7</td><td>13</td><td>7</td><td>9.0</td></tr><tr><td>bta-miR-138</td><td>1.000</td><td>0.070</td><td>0.770</td><td>22</td><td>6</td><td>8</td><td>12.0</td></tr><tr><td>bta-miR-450a</td><td>1.190</td><td>0.140</td><td>0.760</td><td>8</td><td>18</td><td>10</td><td>12.0</td></tr><tr><td>bta-miR-6527</td><td>1.010</td><td>0.090</td><td>0.770</td><td>21</td><td>8</td><td>9</td><td>12.7</td></tr></tbody></table>

<table><caption>**Table 6. Ranking average of the first ten differentially expressed miRNAs in the cellular fraction, the three applied tests' values, and their rankings.**</caption><thead><tr><th rowspan="2">miRNA</th><th colspan="3">Result for</th><th colspan="3">Ranking According to</th><th rowspan="2">Average Ranking</th></tr><tr><th>Test 1 <sup>a</sup></th><th>Test 2 <sup>b</sup></th><th>Test 3 <sup>c</sup></th><th>Test 1</th><th>Test 2</th><th>Test 3</th></tr></thead><tbody><tr><td>bta-miR-574 *</td><td>5.770</td><td>0.000</td><td>0.990</td><td>1</td><td>1</td><td>1</td><td>1.0</td></tr><tr><td>bta-miR-3432a *</td><td>5.520</td><td>0.010</td><td>0.980</td><td>2</td><td>3</td><td>2</td><td>2.3</td></tr><tr><td>bta-miR-2285e *</td><td>2.540</td><td>0.010</td><td>0.950</td><td>5</td><td>2</td><td>3</td><td>3.3</td></tr><tr><td>bta-miR-197 *</td><td>1.970</td><td>0.010</td><td>0.920</td><td>6</td><td>4</td><td>5</td><td>5.0</td></tr><tr><td>bta-miR-2284y *</td><td>2.750</td><td>0.020</td><td>0.940</td><td>3</td><td>8</td><td>4</td><td>5.0</td></tr><tr><td>bta-miR-219</td><td>1.740</td><td>0.010</td><td>0.910</td><td>9</td><td>5</td><td>7</td><td>7.0</td></tr><tr><td>bta-miR-2397-3p</td><td>1.770</td><td>0.020</td><td>0.900</td><td>8</td><td>7</td><td>8</td><td>7.7</td></tr><tr><td>bta-miR-2308</td><td>2.560</td><td>0.050</td><td>0.910</td><td>4</td><td>14</td><td>6</td><td>8.0</td></tr><tr><td>bta-miR-2419-5p</td><td>1.620</td><td>0.020</td><td>0.890</td><td>11</td><td>6</td><td>9</td><td>8.7</td></tr><tr><td>bta-miR-2409</td><td>1.790</td><td>0.040</td><td>0.890</td><td>7</td><td>12</td><td>10</td><td>9.7</td></tr></tbody></table>

Six miRNAs in each milk fraction were chosen as candidates for normalization (Table 7). Following GeNorm, normalization of levels in the fat fraction was optimal using the geometric mean of bta-miR-151-3p and bta-miR-30a-5p. In the case of cellular fraction, GeNorm recommended the use of the geometric mean of the three most stable miRNAs (bta-miR-103, bta-miR-107, bta-miR-28).

Based on analysis of the normalized miRNAs levels (Figures 1 and 2), the only miRNA in the fat or cellular fraction that differed significantly between intensive and extensive production systems was bta-miR-215, in the fat fraction. This miRNA was significantly more abundant in the fat fraction of milk from intensive production (p = 0.030). The miRNAs bta-miR-2284y and bta-miR-2285e in the cellular fraction showed a trend towards lower levels in milk from extensive production, but the differences did not reach statistical significance.

<footnote>* These miRNAs were validated using quantitative RT-PCR. <sup>a</sup> Strict Cohen's d test. <sup>b</sup> Student's t test. <sup>c</sup> Correlation coefficient test.</footnote>

<footnote>* These miRNAs were validated using quantitative RT-PCR. <sup>a</sup> Strict Cohen's d test. <sup>b</sup> Student's t test. <sup>c</sup> Correlation coefficient test.</footnote>
<table><caption>**Table 7. Ranking of the most stable miRNAs according to the coefficient of variation.**</caption><thead><tr><th>Milk Fraction</th><th>miRNA</th><th>Coefficient of Variation</th></tr></thead><tbody><tr><td rowspan="6">Fat</td><td>bta-miR-532</td><td>0.060</td></tr><tr><td>bta-miR-151-3p</td><td>0.070</td></tr><tr><td>bta-miR-27b</td><td>0.090</td></tr><tr><td>bta-miR-103</td><td>0.090</td></tr><tr><td>bta-miR-30a-5p</td><td>0.090</td></tr><tr><td>bta-miR-99a-3p</td><td>0.090</td></tr><tr><td rowspan="7">Cellular</td><td>bta-miR-103</td><td>0.080</td></tr><tr><td>bta-miR-107</td><td>0.090</td></tr><tr><td>bta-miR-181a</td><td>0.090</td></tr><tr><td>bta-miR-28</td><td>0.100</td></tr><tr><td>bta-miR-345-3p</td><td>0.100</td></tr><tr><td>bta-miR-28342</td><td>0.100</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Average levels of miRNAs estimated by qRT-PCR in fat fractions of milk from extensive and intensive farms. On the X axis, production systems are represented; on the Y axis, average levels of each miRNA are represented. Levels of bta-miR-369-3p showed a skewed distribution, so they are shown using a box-and-whisker plot. The errors bars in the other plots indicate standard deviation.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Average levels of miRNAs estimated by qRT-PCR in cellular fractions of milk from extensive and intensive farms. On the X axis, production systems are represented; on the Y axis, average levels of each miRNA are represented. The errors bars indicate standard deviation.</figcaption></figure>

### 3.3. *Putative Target Gene and Pathway Analyses*

Using Target Scan, 143 potential target genes were identified for bta-miR-215, which was abundant in intensive compared to extensive dairy production (Table S2). Among these targets, one gene was particularly involved in lipid metabolism and energy metabolism: the fatty acid binding protein 3 (Fabp3).

The pathway analysis of bta-miR-215 target genes allowed for the identification of 41 associated biological pathways. The gonadotropin-releasing hormone receptor pathway and the Transforming Growth Factor β (TGF-β) signaling pathway (Figure 3) were highlighted, as they include the most bta-miR-215 target genes, i.e., Activin Receptor Type-2B (Acvr2b), Activin A Receptor, Type 2A (Acvr2a), Mitogen-Activated Protein Kinase 1 (Mapk1), and bone morphogenetic protein receptor type II (Bmpr2).
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Functional Gene Ontology classification of bta-miR-215 pathway analysis of targeted genes.</figcaption></figure>
## 4. Discussion

The objective of this study was to identify a set of miRNAs in milk whose levels differed with the production system, in order to (1) determine whether miRNAs profiling can be used to authenticate milk from a given production system, and (2) investigate whether the production system can influence the functional properties of agri-food products such as bovine milk. To address these questions, we sampled milk from dairy farms in Asturias in northern Spain that applied intensive or extensive production practices, which differ significantly in numerous characteristics that can alter milk quality, including feeding management, animal density, and access to exercise through grazing [29].

The origin of miRNAs in milk is controversial, and most miRNAs in milk are not found in blood [30]. In addition, the miRNA profile in milk differs across the fractions of fat, whey, and cells [12]. In the present work, we decided not to study the whey fraction because its miRNA content is lower, its miRNA profile is highly similar to that of milk fat [12], and it is not exploited in certain dairy industry practices, such as cheese production.

We successfully isolated total RNA from fat and cellular fractions of bovine milk, especially considering that the samples came from milk tanks on commercial dairy farms. The RIN was low for RNA from milk fat, which likely reflects its abundant content of low-molecular-weight RNA, which is different from ribosomal RNA [12].

RNA sequencing analysis confirmed different RNA profiles in the fat and cellular fractions of bovine milk (Table 3). Interestingly, we did not identify significant differences in miRNAs levels between the two fractions, although we did obtain what appears to be the first evidence that miRNAs in the fat fraction differ between milk from intensive or extensive production systems. Given that milk has been proposed as a major epigenetic modulator of the gene expression of the milk recipient [31], and its modulation can be dependent on the amount of miRNAs [32], our results imply that the animal production system can influence the functional properties of agri-food products of animal origin.

Sequencing results did not allow for the identification of miRNAs that were specific to a given production system in either the fat or cellular fraction. Therefore, we focused on miRNAs whose levels differed between the two systems, and we validated a subset of the promising miRNAs using quantitative RT-PCR. The only miRNA that we validated to differ significantly between the two production systems was bta-miR-215 in the fat fraction (Table 5). This miRNA was upregulated in milk from intensive production. At the moment, we can only speculate as to why intensive production might upregulate this miRNA. One possibility is that its upregulation somehow compensates for poor efficiency of feed conversion into milk: dairy cows with medium potential show lower conversion efficiency during indoor feeding than during grazing in pastures [33], and Angus cows, less efficient at feed conversion, show upregulation of bta-miR-215 [34]. Consistent with this possibility is that heat and oxidative stress, which can easily occur in intensive production systems, upregulate bta-miR-215 in the serum of Holstein cows [35].

Among the 143 genes identified by bioinformatic tools, Fabp3 stands out for its involvement in fatty acid transport and activation in the bovine mammary gland [36]. Vargas-Bello-Perez et al. (2020) [37] demonstrated that when the diet of Holstein dairy cows is supplemented with hydrogenated vegetable oil, Fabp3 is downregulated in milk somatic cells. Likewise, when fed a high-concentrate diet, the fatty acid transporter Fabp3 is inhibited in mammary glands [38]. In our study, bta-miR-215 increased in milk in response to intensive farming conditions, which are generally associated with high-energy diets, which is consistent with its likely involvement in the downregulation of Fabp3. A direct relationship between bta-miR-215 and Fabp3 has, indeed, been validated in bone marrow mesenchymal stem cells [39].

Apart from Fabp3, another two genes are of particular interest: Mapk1 and Bmpr2. Mapk1 is involved in the regulation of milk protein synthesis [40], and Bmpr2 is associated with the glucose metabolism and insulin response. In fact, when Bmpr2 is altered, it likely blunts glucose response and lipid uptake [41]. Altogether, they likely have a negative impact on protein and lactose synthesis in animals bred under intensive dairy systems.
We also identified the miRNAs bta-mir-2284y and bta-mir-2285e in the cellular fraction of milk whose levels tended to differ between intensive and extensive production systems. Further analysis showed that these differences became significant when we compared animals on diets containing more or fewer than 10 kg of concentrates per day (data not shown). These results suggest that our failure to detect significant differences in levels of bta-mir-2284y and bta-mi-2285e in the overall analysis may reflect milk sampling from commercial farms, as well as the difficulty of defining intensive and extensive production. Further studies should examine these two miRNAs as additional potential biomarkers for authentication and functional analysis of agri-foods.

Finally, our study presents several limitations. The study was carried out in commercial farms, where managing conditions and diets are quite different even within the same group. Validation in controlled conditions at experimental farms, which reduce internal variance of the experimental groups, should help to identify miRNAs as putative biomarkers for dairy production systems. In vitro experiments will help to elucidate the functional features of milk produced under different dairy systems.

**5. Conclusions**

We investigated differences in miRNA profiles of raw cow tank milk from commercial farms which applied either intensive or extensive production systems. We identified bta-miRNA-215 in the fat fraction of milk as a possible biomarker of milk from intensive production systems. Our results imply that the type of production system can influence miRNA levels, and, therefore, functional properties of bovine milk, as well as, potentially, other agri-food products of animal origin.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/vetsci9120661/s1, Table S1: Sequencing RPM; Table S2: Selection of Predicted bta-miR-215 Targets from Target Scan database.

**Author Contributions:** Conceptualization, L.J.R., L.A.e.q., and L.G.; methodology, J.A., J.D., L.J.R., L.A.e.q., S.L.G., and K.Z.; formal analysis, L.A.e.q., S.L.G., M.F.-S., E.I.-G., and L.J.R.; investigation, L.A.e.q., L.J.R.; resources, L.A.e.q. and F.V.; data curation, L.J.R., K.Z., and L.A.e.q.; writing—original draft preparation, L.A.e.q. and L.J.R.; writing—review and editing, L.G., K.Z., S.L.G., M.F.-S., E.I.-G., and F.V.; visualization, L.A.e.q. and L.J.R.; supervision, L.J.R. and L.G.; project administration, L.J.R.; funding acquisition, L.J.R., L.G., and F.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Principado de Asturias Regional Government (project IDI/2021/000102), co-financed by the European Union through the ERDF (European Regional Development Fund) and Ministry of Science and Innovation (AEI) project PID2021-126010OR-I00. Loubna Abou el qassim was founded by FICYT Severo Ochoa Grant (BP17-49).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We want to acknowledge all dairy farmers that contribute to this study.

**Conflicts of Interest:** The authors declare no conflict of interest.

**References**

1. Luykx, D.M.A.M.; van Ruth, S.M. An Overview of Analytical Methods for Determining the Geographical Origin of Food Products. *Food Chem.* **2008**, *107*, 897–911. [CrossRef]

2. Valenti, B.; Martin, B.; Andueza, D.; Leroux, C.; Labonne, C.; Lahalle, F.; Larroque, H.; Brunschwig, P.; Lecomte, C.; Brochard, M.; et al. Infrared Spectroscopic Methods for the Discrimination of Cows' Milk According to the Feeding System, Cow Breed and Altitude of the Dairy Farm. *Int. Dairy J.* **2013**, *32*, 26–32. [CrossRef]

3. He, L.; Hannon, G.J. MicroRNAs: Small RNAs with a Big Role in Gene Regulation. *Nat. Rev. Genet.* **2004**, *5*, 522–531. [CrossRef] [PubMed]
4. Thorning, T.K.; Raben, A.; Tholstrup, T.; Soedamah-Muthu, S.S.; Givens, I.; Astrup, A. Milk and Dairy Products: Good or Bad for Human Health? An Assessment of the Totality of Scientific Evidence. *Food Nutr. Res.* **2016**, *60*, 32527. [CrossRef]

5. López de Las Hazas, M.C.; Del Pozo-Acebo, L.; Hansen, M.S.; Gil-Zamorano, J.; Mantilla-Escalante, D.C.; Gómez-Coronado, D.; Marín, F.; Garcia-Ruiz, A.; Rasmussen, J.T.; Dávalos, A. Dietary bovine milk miRNAs transported in extracellular vesicles are partially stable during GI digestion, are bioavailable and reach target tissues but need a minimum dose to impact on gene expression. *Eur. J. Nutr.* **2022**, *61*, 1043–1056. [CrossRef]

6. Billa, P.A.; Faulconnier, Y.; Ye, T.; Chervet, M.; Le Provost, F.; Pires, J.A.A.; Leroux, C. Deep RNA-Seq Reveals MiRNome Differences in Mammary Tissue of Lactating Holstein and Montbéliarde Cows. *BMC Genom.* **2019**, *20*, 621. [CrossRef]

7. Le Guillou, S.; Leduc, A.; Laubier, J.; Barbey, S.; Rossignol, M.N.; Lefebvre, R.; Marthey, S.; Laloë, D.; Le Provost, F. Characterization of Holstein and Normande Whole Milk MiRNomes Highlights Breed Specificities. *Sci. Rep.* **2019**, *9*, 20345. [CrossRef]

8. Miretti, S.; Lecchi, C.; Ceciliani, F.; Baratta, M. MicroRNAs as Biomarkers for Animal Health and Welfare in Livestock. *Front. Vet. Sci.* **2020**, *7*, 578193. [CrossRef]

9. Tzelos, T.; Ho, W.; Charmana, V.I.; Lee, S.; Donadeu, F.X. MiRNAs in Milk Can Be Used towards Early Prediction of Mammary Gland Inflammation in Cattle. *Sci. Rep.* **2022**, *12*, 131. [CrossRef]

10. Muroya, S.; Ogasawara, H.; Hojito, M. Grazing Affects Exosomal Circulating MicroRNAs in Cattle. *PLoS ONE* **2015**, *10*, e0136475. [CrossRef]

11. Muroya, S.; Shibata, M.; Hayashi, M.; Oe, M.; Ojima, K. Differences in Circulating MicroRNAs between Grazing and Grain-Fed Wagyu Cattle Are Associated with Altered Expression of Intramuscular MicroRNA, the Potential Target PTEN, and Lipogenic Genes. *PLoS ONE* **2016**, *11*, e0162496. [CrossRef] [PubMed]

12. Li, R.; Dudemaine, P.L.; Zhao, X.; Lei, C.; Ibeagha-Awemu, E.M. Comparative Analysis of the MiRNome of Bovine Milk Fat, Whey and Cells. *PLoS ONE* **2016**, *11*, e0154129. [CrossRef] [PubMed]

13. Zhang, L.; Hou, D.; Chen, X.; Li, D.; Zhu, L.; Zhang, Y.; Li, J.; Bian, Z.; Liang, X.; Cai, X.; et al. Exogenous Plant MIR168a Specifically Targets Mammalian LDLRAP1: Evidence of Cross-Kingdom Regulation by MicroRNA. *Cell Res.* **2012**, *22*, 107–126. [CrossRef]

14. Chen, X.; Gao, C.; Li, H.; Huang, L.; Sun, Q.; Dong, Y.; Tian, C.; Gao, S.; Dong, H.; Guan, D.; et al. Identification and Characterization of MicroRNAs in Raw Milk during Different Periods of Lactation, Commercial Fluid, and Powdered Milk Products. *Cell Res.* **2010**, *20*, 1128–1137. [CrossRef] [PubMed]

15. Izumi, H.; Kosaka, N.; Shimizu, T.; Sekine, K.; Ochiya, T.; Takase, M. Bovine Milk Contains MicroRNA and Messenger RNA That Are Stable under Degradative Conditions. *J. Dairy Sci.* **2012**, *95*, 4831–4841. [CrossRef]

16. Bauman, D.E.; Harvatine, K.J.; Lock, A.L. Nutrigenomics, Rumen-Derived Bioactive Fatty Acids, and the Regulation of Milk Fat Synthesis. *Annu. Rev. Nutr.* **2011**, *31*, 299–319. [CrossRef] [PubMed]

17. Ollier, S.; Leroux, C.; de la Foye, A.; Bernard, L.; Rouel, J.; Chilliard, Y. Whole Intact Rapeseeds or Sunflower Oil in High-Forage or High-Concentrate Diets Affects Milk Yield, Milk Composition, and Mammary Gene Expression Profile in Goats. *J. Dairy Sci.* **2009**, *92*, 5544–5560. [CrossRef] [PubMed]

18. Quan, S.; Nan, X.; Wang, K.; Zhao, Y.; Jiang, L.; Yao, J.; Xiong, B. Replacement of Forage Fiber with Non-Forage Fiber Sources in Dairy Cow Diets Changes Milk Extracellular Vesicle-MiRNA Expression. *Food Funct.* **2020**, *11*, 2154–2162. [CrossRef] [PubMed]

19. Marin, E.C. Sistema de Explotación Ganadera: Notas En Torno a Su Concepto. *Lurralde* **1996**, *19*, 89–104.

20. Ruiz, J.; Herrera, M.P.; Barba, R.; Busqué, J. *Definición y Caracterización de la Extensividad en las Explotaciones Ganaderas en España. Situación de la Ganadería Extensiva en España* (I); Ministerio de Agricultura y Pesca, Alimentación y Medio Ambiente: Madrid, Spain, 2017.

21. Cohen, J. *Statistical Power Analysis for the Behavioral Sciences*, 2nd ed.; Erlbaum, L., Ed.; Routledge: London, UK, 1988; ISBN 0805802835.

22. Student. The Probable Error of a Mean. *Biometrika* **1908**, *6*, 1–25. [CrossRef]

23. Pearson, K., VII. Note on Regression and Inheritance in the Case of Two Parents. *Proc. R. Soc. Lond.* **1895**, *58*, 240–242.

24. Vandesompele, J.; De Preter, K.; Pattyn, F.; Poppe, B.; Van Roy, N.; De Paepe, A.; Speleman, F. Accurate Normalization of Real-Time Quantitative RT-PCR Data by Geometric Averaging of Multiple Internal Control Genes. *Genome Biol.* **2002**, *3*, research0034.1. [CrossRef] [PubMed]

25. Abou el qassim, L. Variaciones de Los Perfiles de MicroARN en la Leche Cruda de Vaca Según el Sistema de Alimentación. Master's Thesis, University of Zaragoza, Zaragoza, Spain, 2017.

26. Hellemans, J.; Mortier, G.; De Paepe, A.; Speleman, F.; Vandesompele, J. QBase Relative Quantification Framework and Software for Management and Automated Analysis of Real-Time Quantitative PCR Data. *Genome Biol.* **2008**, *8*, R19. [CrossRef] [PubMed]

27. Agarwal, V.; Bell, G.W.; Nam, J.-W.; Barte, D.P. Predicting Effective MicroRNA Target Sites in Mammalian MRNAs. *eLife* **2015**, *4*, e05005. [CrossRef] [PubMed]

28. Mi, H.; Ebert, D.; Muruganujan, A.; Mills, C.; Albou, L.P.; Mushayamaha, T.; Thomas, P.D. PANTHER Version 16: A Revised Family Classification, Tree-Based Classification Tool, Enhancer Regions and Extensive API. *Nucleic Acids Res.* **2021**, *49*, D394–D403. [CrossRef] [PubMed]

29. Vicente, F.; Santiago, C.; Jiménez-Calderón, J.D.; Martínez-Fernández, A. Capacity of Milk Composition to Identify the Feeding System Used to Feed Dairy Cows. *J. Dairy Res.* **2017**, *84*, 254–263. [CrossRef] [PubMed]
30. Alsweed, M.; Lai, C.T.; Hartmann, P.E.; Geddes, D.T.; Kakulas, F. Human Milk MiRNAs Primarily Originate from the Mammary Gland Resulting in Unique MiRNA Profiles of Fractionated Milk. *Sci. Rep.* **2016**, *6*, 20680. [CrossRef] [PubMed]

31. Melnik, B.C.; Kakulas, F. Milk Exosomes and MicroRNAs: Potential Epigenetic Regulators. In *Handbook of Nutrition, Diet, and Epigenetics*; Springer International Publishing: Berlin/Heidelberg, Germany, 2019; Volume 3, pp. 1467–1494. ISBN 9783319555300.

32. Del Pozo-Acebo, L.; de las Hazas, M.C.L.; Tomé-Carneiro, J.; Gil-Cabrerizo, P.; San-Cristobal, R.; Busto, R.; García-Ruiz, A.; Dávalos, A. Bovine Milk-Derived Exosomes as a Drug Delivery Vehicle for Mirna-Based Therapy. *Int. J. Mol. Sci.* **2021**, *22*, 1105. [CrossRef]

33. Knaus, W. Perspectives on Pasture versus Indoor Feeding of Dairy Cows. *J. Sci. Food Agric.* **2016**, *96*, 9–17. [CrossRef]

34. Al-Husseini, W.; Chen, Y.; Gondro, C.; Herd, R.M.; Gibson, J.P.; Arthur, P.F. 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. [CrossRef]

35. Zheng, Y.; Chen, K.L.; Zheng, X.M.; Li, H.X.; Wang, G.L. Identification and Bioinformatics Analysis of MicroRNAs Associated with Stress and Immune Response in Serum of Heat-Stressed and Normal Holstein Cows. *Cell Stress Chaperones* **2014**, *19*, 973–981. [CrossRef] [PubMed]

36. Zhang, H.; Shen, Z.; Yang, Z.; Jiang, H.; Chu, S.; Mao, Y.; Li, M.; Chen, Z.; Aboragah, A.; Loor, J.J.; et al. Abundance of Solute Carrier Family 27 Member 6 (SLC27A6) in the Bovine Mammary Gland Alters Fatty Acid Metabolism. *Food Funct.* **2021**, *12*, 4909–4920. [CrossRef] [PubMed]

37. Vargas-Bello-pérez, E.; Geldsetzer-Mendoza, C.; Cancino-Padilla, N.; Morales, M.S.; Leskinen, H.; Garnsworthy, P.C.; Loor, J.J.; Romero, J. Effects of Dietary Vegetable Oils on Mammary Lipid-Related Genes in Holstein Dairy Cows. *Animals* **2020**, *10*, 57. [CrossRef]

38. Ma, N.; Abaker, J.A.; Wei, G.; Chen, H.; Shen, X.; Chang, G. A High-Concentrate Diet Induces an Inflammatory Response and Oxidative Stress and Depresses Milk Fat Synthesis in the Mammary Gland of Dairy Cows. *J. Dairy Sci.* **2022**, *105*, 5493–5505. [CrossRef]

39. Zhou, Q.; Yang, H.; Pan, H.; Pan, H.; Zhou, J. Exosomes Isolated from the MiR-215-Modified Bone Marrow Mesenchymal Stem Cells Protect H2O2-Induced Rat Myoblasts via the MiR-215/FABP3 Pathway. *Exp. Mol. Pathol.* **2021**, *119*, 104608. [CrossRef]

40. Lu, L.M.; Li, Q.Z.; Huang, J.G.; Gao, X.J. Proteomic and Functional Analyses Reveal MAPK1 Regulates Milk Protein Synthesis. *Molecules* **2013**, *18*, 263–275. [CrossRef]

41. Hemnes, A.; Fessel, J.; Chen, X.; Zhu, S.; West, J.; Talati, M. BMPR2 Dysfunction Impairs Insulin Signaling and Glucose Homeostasis in Cardiomyocytes. *Am. J. Physiol. Lung Cell. Mol. Physiol.* **2019**, *318*, L429–L441. [CrossRef] [PubMed]