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

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.


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 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][9][10][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][14][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.

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. 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).  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.

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).

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 6 .

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).

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.

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.

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].

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.

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].

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.  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).

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.  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.

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).

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).

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.

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.  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).

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.