Grazing affects metabolic pattern of individual cow milk

Effective traceability tools able to characterize milk from pasture are important to safeguard low-input farming systems, niche dairy products, and local tradi-tions. The aims of the present study were to investigate the ability of proton nuclear magnetic resonance ( 1 H NMR) spectroscopy to discriminate between milk produced from cows before and after the beginning of the grazing season, and to assess the effects of grazing on milk metabolites. The research trial involved a single alpine holding with 72 lactating cows. Individual milks were repeatedly sampled from the same animals before (i.e., d −3 and −1) and after (i.e., d 2, 3, 7, 10, and 14) the onset of the grazing period. One-dimensional 1 H NMR spectra of milk extracts were collected through a Bruker spectrometer. Random forest discriminant analysis was applied to 1 H NMR spectra to predict the period of collection for each sample. Data concerning the relative abundance of milk metabolites were analyzed through a linear mixed model, which included the fixed effects of period of sampling, cow breed, stage of lactation, and parity, and the random effect of cow nested within breed. The random forest model exhibited great accuracy (93.1%) in discriminating between samples collected on d −3, −1, 2, and 3 and those collected on d 7, 10, and 14. Univariate analysis performed on the 40 detected metabolites highlighted that milk samples from pasture had lower levels of 14 compounds (with fumarate being the most depressed metabolite) and greater levels of 15 compounds (with methanol and hippurate being the most elevated metabolites). Results indicate that milk 1 H NMR spectra are promising to identify milk produced in different conditions. Also, our study highlights that grazing is associated with significant changes of milk metabolic profile, suggesting the potential use of several metabolites as indicators of farm management.


INTRODUCTION
The use of pasture as a source of fresh forages for cattle is commonly adopted in the context of dairy farming in mountainous areas. Fresh herbages are administered as animal feed through different strategies, depending on farm location, management, and facilities. Usually, in intensive farming systems, summer forages are mowed, transported to the barn, and mixed with other feed ingredients as fresh herbages, hay, or silage. By contrast, in low-input farming systems that are widely diffused in marginal or mountainous areas, dairy farmers allow the herd to have direct access to pasture. Depending on the geographical position of the barn, grazing can be practiced according to different strategies: the first is adopted when pastures are adjacent to the barn, allowing animals to autonomously move from barn to pastures after morning and evening milking (Niero et al., 2021a); the second, known as transhumance, is adopted when pastures are far from the barn, and includes active moving or passive transport of animals from the lowland barn to highland pastures (Niero et al., 2018). Regardless of the scenario, it is acknowledged that the inclusion of fresh forages in animal diet is threatened by high labor demands and relatively low financial returns. Still, this practice remains of strategic significance due to its direct and indirect implications, ranging from environmental sustainability and animal welfare (Burow et al., 2013;Byrnes et al., 2018) to the quality of animal-derived products (Niero et al., 2021b), the touristic appeal of the territory (Bele et al., 2017), consumers' appreciation of the final product (Jackson et al., 2020), and the maintenance of local traditions, biodiversity, and landscapes (Niero et al., 2018).
Given the notable importance of this topic, it is worth developing and implementing analytical methods based on inherent food properties that are able to distinguish between dairy foods originating from conventional and pasture-based conditions. For instance, it is well documented that both gas chromatography and mid-infrared spectroscopy techniques can be used to effectively monitor the modification of detailed milk fatty acid composition, where an increased concentration of polyunsaturated fatty acids can be used as a possible biomarker for pasture-based farming systems (Leiber et al., 2005;De Marchi et al., 2018). Mid-infrared spectra can be also used as a fingerprinting technique to discriminate between milk produced under intensive and extensive farming conditions (Capuano et al., 2014). Furthermore, the profiling of stable isotopes and isotope ratios has been proposed as a traceability tool able to discriminate between organically and conventionally produced milk (Molkentin and Giesemann, 2010) and between milks obtained through different feeding regimens (Camin et al., 2008).
In the realm of food science, proton nuclear magnetic resonance ( 1 H NMR) spectroscopy is one of the most potent analytical methods that may be used to analyze liquid or solid materials. Proton nuclear magnetic resonance has enormous potential for analyzing complex matrices and offers opportunity to advance the field of food science due to its nondestructive nature, high accuracy, and reproducibility, which are frequently achieved without using any separation or purification steps. The area of the 1 H NMR signal is exactly proportional to the number of nuclei that create the signal in 1 H NMR spectroscopy, which is a quantitative method that may be used to analyze a variety of nuclei based on the type of food being analyzed and the information being sought. Moreover, even though conventional quality, safety, and authenticity control in food is based on targeted strategies, high-resolution 1 H-NMR, performed under well-defined instrumental specifications, offers several advantages, including producing an extremely reproducible food fingerprint and fully quantitative data by means of a single experiment (Vignoli et al., 2019). In this scenario, 1 H NMR spectroscopy has been effectively used in the dairy sector for different purposes and applications, primarily traceability of milk and cheese produced according to different farming systems (Tenori et al., 2018;Segato et al., 2019), detection of bovine milk adulteration in caprine milk (Rysova et al., 2021), and characterization of metabolites in milk from cows with subclinical and clinical mastitis (Luangwilai et al., 2021). Still, a paucity of information exists regarding the possibility to use 1 H NMR spectroscopy as an analytical traceability tool for the distinction between intensive and extensive farming conditions. Similarly, few studies investigating the effect of grazing pasture on milk metabolome exist so far.
These research questions have been targeted through the trials of the present study, which involved dairy cows farmed on a single alpine herd adopting summer grazing. Specific aims were to (1) assess the ability of 1 H NMR spectroscopy to discriminate between milk produced from cows before and after the beginning of the grazing season and (2) investigate how, and to what extent, grazing can affect milk metabolic profile.

Herd
The trials were performed during routine milking procedures and were not invasive; therefore, animal welfare committee authorization was not required. The commercial dairy farm involved in this study was located in Veneto Region (Pian del Cansiglio, Belluno, Italy) at 1,100 m above sea level. The farm was enrolled in the official milk recording system of the Breeders Association of Veneto Region (Vicenza, Italy). The herd comprised Simmental (SI), Holstein (HO), and HO × SI crossbred (CR) lactating cows (n = 72). Between October and May, animals were housed indoors in a freestall barn and were fed with TMR based on grass silage, alfalfa hay, high-moisture corn, cereal meal, and a protein-mineral-vitamin mix (9.5 kg as fed, 3 kg, 7 kg, 3 kg, and 3.5 kg, respectively). Milking took place twice per day, in the morning (0500 h) and in the afternoon (1700 h). Between late spring and the beginning of autumn, animals had access to the pasture located near the farm. Cows grazed fresh herbage following the Voisin rotational system (Voisin, 1959), with animals being moved to a different paddock every day. Cows received 2 kg of high-moisture corn and cereal meals as an energy supplement during each milking event. Further details on chemical composition of the fresh herbage grazed by the animals can be retrieved from Niero et al. (2021a).

Sample Collection and Gross Chemical Composition
Individual milk samples from the entire herd were repeatedly collected during the afternoon milking from late May to June 2020. Milk sampling started 3 d and 1 d before the beginning of grazing, to characterize milk obtained under indoor farming conditions. Milk samples were also collected 2, 3, 7, 10, and 14 d after the beginning of grazing, to characterize the milk obtained at pasture (Figure 1). After sample collection, 200 µL of Bronopol (2-bromo-2-nitropropan-1,3-diol; ANA. LI.TIK. Austria) preservative were added to 40 mL of milk. Samples were transferred at 4°C to the laboratory of the Breeders Association of Veneto Region (Padova, Italy) and analyzed within 12 h for gross chemical composition (fat, protein, casein, and lactose percentages) and urea content (mg/dL), using a MilkoScan FT6000 (Foss Analytical A/S). A Fossomatic 7 DC (Foss Analytical A/S) was used to determine SCC (cells/ µL) and differential SCC (%). Values of SCC were log-transformed to SCS using the formula SCS = 3 + log 2 (SCC/100) to achieve the normality and homogeneity of variances (Ali and Shook, 1980). A total of 451 milk samples from 72 cows were collected for 1 H NMR analysis (291 samples from 49 SI cows; 21 samples from 3 HO cows; 139 samples from 20 CR cows).
One-dimensional (1D) 1 H NMR spectra of milk aqueous extracts were recorded on a Bruker spectrometer operating at 600.13-MHz proton Larmor frequency and equipped with a 5-mm PATXI 1 H-13C-15N probe including a z-axis gradient coil, automatic tuningmatching, and an automatically refrigerated samplechanger (SampleJet, Bruker BioSpin). A BTO 2000 thermocouple provided thermal stabilization of the samples, ensuring a maximum fluctuation of approximately 0.1 K. For temperature equilibration, samples were maintained inside the NMR probe for 5 min prior to measurement (310 K). Nuclear Overhauser effect spectroscopy (NOESY) pulse sequence (noesygppr1d, Bruker BioSpin) was used to acquire 1D 1 H NMR spectra with 64 scans, 98 k data points, a spectral width of 18,028 Hz, an acquisition time of 2.7 s, a relaxation time of 4 s, and a mixing time of 10 ms. The NOESY pulse sequence suppresses solvent signal and yields a spectrum in which both the signals of metabolites and high molecular weight molecules are visible (McKay, 2011). Before applying Fourier transformation, free induction decays were multiplied by an exponential function equivalent to a 0.3-Hz line-broadening factor. TopSpin software (Bruker) was used to automatically adjust for phase and baseline aberrations in transformed spectra. To calibrate spectra, the α-lactose doublet (5.24 δ 1 H ppm) was used. Each 1D spectra in the range of 0.02 to 10.00 δ 1 H ppm was segmented into 0.05 δ 1 H ppm chemical shift buckets for multivariate analysis. Buckets are spectral areas calculated under spectral segments, reported in arbitrary units. Prior to performing the statistical analyses, the portions of NO-ESY spectra corresponding to the water (4.61-4.77 δ 1 H ppm) and dichloromethane (5.30-5.33 and 5.42-5.65 δ 1 H ppm) resonating regions were removed and the buckets were normalized using probabilistic quotient normalization (Dieterle et al., 2006), resulting in a data matrix of 451 rows (samples) and 228 columns (buckets). A total of 40 signals were characterized in the 1 H NMR spectra (Supplemental Table S1, https: / / figshare .com/ articles/ figure/ Supplemental _Table _S1 _pdf/ 20712427; Visentin, 2022a), resulting in a data matrix of 451 rows (samples) and 40 columns (metabolites). Signal identification was performed using a library of NMR spectra of pure organic compounds (Assure NMR 2.2 software, Bruker BioSpin), public databases storing references (FooDB, https: / / foodb .ca/ ; Milk Composition Database, http: / / www .mcdb .ca/ ), and literature data (Tenori et al., 2018;Meoni et al., 2020).

Data Editing and Statistical Analyses
Among the 451 acquired records, only samples that satisfied specific inclusion criteria were used for statistical analyses. In particular, (1) DIM and parity were restricted to range from 5 to 500 d and between 1 and 8, respectively; (2) records with milk yield (MY) <2 kg/milking were discarded; and (3) only samples collected from animals that were present in at least 5 of the 7 test-days (days before and after the beginning of grazing) were kept. Additionally, for MY and milk quality traits, samples exceeding 3 standard deviations from the respective mean were treated as missing values. Following these criteria, 421 milk samples from 72 animals (256 samples from 49 SI cows; 21 samples from 3 HO cows; 135 samples from 20 CR cows) remained for statistical analysis. A total of 118 individual milks were collected 3 d and 1 d before the beginning of grazing, and 303 samples were obtained 2, 3, 7, 10, and 14 d after the beginning of grazing. Descriptive statistics of production-related traits, MY, and quality traits are presented in Supplemental Table S2 (https: / / figshare .com/ articles/ figure/ Untitled _Item/ 20712499; Visentin, 2022b).
Data analyses were performed using R software (version 3.5.3; R Core Team, 2020). Principal component analysis was used as the first exploratory unsupervised analysis on bucketed spectra, to evaluate the presence of a net discrimination between the days of sampling and potential presence of outliers. Random forest (RF) discriminant analysis ("Random Forest" R package) was applied to bucketed spectra. Here, each tree of the forest tree was used to predict the day of collection for each sample. For this study, RF was used with double cross-validation, and samples from various collection times (days) of the same animal were deleted from the training at each test cycle of the double cross-validation to obtain an unbiased predictive accuracy. For each computation, a forest of 5,000 trees was employed.
The univariate analysis was carried out on the 40 identified metabolites integrated from 1D 1 H-NMR spectra. Data were processed according to the following linear mixed model, using PROC MIXED of SAS software v. 9.4 (SAS Institute Inc.): where y ijklm is a milk metabolite; period i is the fixed effect of the ith period of sampling [2 classes: d −3 to d 3 (representing TMR and early-stage grazing diets) and d 7 to d 14 (representing prolonged grazing diet)]; breed j is the fixed effect of the jth cow breed (3 classes: SI, HO, and CR); stage k is the fixed effect of the kth stage of lactation (8 classes: 5-45, 46-90, 91-135, 136-180, 181-225, 226-270, 271-315, and >315 d); parity l is the fixed effect of the lth parity (6 classes: 1, 2, 3, 4, 5, and >5); cow m is the random effect of the mth cow nested within the jth breed; and e ijklm is the random residual term. Variation of milk metabolites in the 2 sampling periods were reported as log 2 (fold change), and differences were tested with the Wilcoxon test with Bonferroni correction for P-value determination. Furthermore, for each metabolite, the Cliff's delta effect size was calculated by means of the R package "effsize" (Torchiano, 2020).

Classification of Samples Based on 1 H NMR Spectra
The capacity of the RF model to discriminate between the different days of sampling is reported in Table 1. The overall accuracy of the model was 58.7%, and the proportion of samples correctly classified by the model (i.e., the sensitivity) ranged from 37.1% (d 3) to 80.3% (d 14). In general, the proportion of samples wrongly attributed was greater between adjacent days of sampling. For example, 35.6% of samples collected in d −3 (i.e., 3 d before access to pasture) were attributed to d −1 (i.e., the day before access to pasture), whereas almost none of such samples were attributed to samples collected from d 7 onward (Table 1). It is worth noting that samples collected from d −3 to d 3 tended to be confounded among each other, as well as samples collected from d 7 to 14. This suggests that milk NMR spectra, and therefore the metabolic profile of dairy cows, is to some extent influenced by the diet of the dairy cows, which become evident only after some days of continued access to pasture and fresh herbage (i.e., from d 7 onward). This is also in agreement with previous studies, which reported that several days of metabolic and behavioral adaptation are needed when animals move from TMR to pasture (Schären et al., 2016). Accordingly, studies involving Latin square designs commonly include a transition period before starting new sampling with different nutritional regimen, to avoid the influence of the previous treatment.
We therefore applied the RF model to only 2 periods of days of sampling, created by grouping together samples collected in d −3, −1, 2, and 3 (period 1) and samples collected in d 7, 10, and 14 (period 2; Table 2). In this case, the overall accuracy of the model increased to 93.1%, with 12.9% of samples in period 2 attributed to period 1, and only 2.4% of samples in period 1 attributed to period 2 (Table 2). It is clear, in particular from Table 2, that exploitable information exists in milk NMR spectra to discriminate between milk samples collected from cows before moving to pasture and cows at pasture. The effect of feeding regimen on the metabolic profile of dairy cows has been reported in previous studies (Tenori et al., 2018;Billa et al., 2020). For example, a restriction in concentrate administration is associated with increased milk glucose-6-phosphate and isocitrate, and with decreased milk BHB, glucose, glutamate, uric acid, and free amino acids concentrations (Billa et al., 2020). Tenori et al. (2018) applied a canonical analysis to partial least squares latent variables to discriminate milk samples based on their herd of origin. This is extremely useful, for example, for food traceability purposes or in a production context in which the administration of some feedstuffs is prohibited by specific voluntary regulations, as in the case of the Parmigiano Reggiano cheese area. In this production area, located in some provinces of the Emilia-Romagna and Lombardy regions (Northern Italy), silage is not permitted as feed for lactating animals, to minimize the risk of anaerobic bacteria presence in milk, which could potentially lead to undesired fermentation during cheese ripening (Buonaiuto et al., 2021). Exploitable information to discriminate between diets offered to dairy cattle is also available in spectra generated from milk mid-infrared spectroscopy analyses. Frizzarin et al. (2021) applied different machine-learning algorithms to more than 4,000 milk spectra collected from grazing or TMR-fed dairy cows and demonstrated that linear discriminant analysis and partial least squares discrimi-  The diagonal of the confusion matrix reports the sensitivity (%) for the classification of each animal. Overall predictive accuracy = 93.1%.
2 Period 1 refers to 3 and 1 d before the beginning of grazing and 2 and 3 d after the beginning of grazing; period 2 refers to 7, 10, and 14 d after the beginning of grazing. nant analysis were the most effective algorithms, with an overall accuracy of 0.968 in both cases.

Phenotypic Correlations
Pearson correlation coefficients (r) between milk metabolites, milk chemical composition, and SCS are presented in Figure 2. Only 2 out of the 40 detected metabolites (namely citrate and galactose-1-phosphate) were positively associated with MY. The remaining metabolites were negatively correlated with MY, meaning that higher MY was associated with a decreased relative abundance of these compounds, probably due to a dilution effect.
N-Acetyl carbohydrates were found to positively correlate with both protein (r = 0.61; P < 0.001) and CN content (r = 0.62; P < 0.001). The specific role of N-acetyl carbohydrates in milk is still largely unknown (Sischo et al., 2017;Giorgio et al., 2018;Rysova et al., 2021), even if the physiological and biochemical patterns regarding these compounds are gradually being elucidated. Consistent with the negative association between N-acetyl carbohydrates and MY observed in this study, Tomassini et al. (2019) and Zhu at al. (2020) reported an increasing trend of N-acetyl carbohydrates across the lactation period, probably due to qualitative and quantitative changes in the synthesis of oligosaccharides in different lactation stages. It has been reported that N-acetyl carbohydrates can originate from catabolism of glycoproteins. In particular, during processes involving proteolysis, glycans become detached and are released from proteins, thereby becoming detectable in 1 H NMR spectroscopy (Rysova et al., 2021). This argument can explain the positive association between protein content and relative abundance of N-acetyl carbohydrates observed in the present study.
Citrate exhibited weak negative correlations with both protein and CN content (r = −0.22; P < 0.001). Milk citrate acts as a calcium chelating agent, creating soluble calcium-citrate complexes that are in competition with colloidal calcium-phosphate complexes. Therefore, increased concentrations of soluble citrate are likely to decrease the level of colloidal calcium phosphate, weakening the stability of casein micelles and, ultimately, deteriorating milk coagulation properties (Sundekilde et al., 2011).
Lecithin was positively correlated with fat (r = 0.32; P < 0.001), protein (r = 0.52; P < 0.001), and CN content (r = 0.51; P < 0.001). The positive association between lecithin and milk fat was expected, due to the fat-soluble properties of this compound, which is found in liposome structures (Holmes et al., 2000). As in the present study, Holmes et al. (2000) reported a direct association between milk lecithin and milk protein content, even if the physiological reasons at the basis of this phenomenon need further detailed study. Moreover, milk lecithin has been discussed by Bobbo et al. (2022), who observed that milk samples with greater SCC had lower levels of lecithin.
As expected, a strong positive correlation (r = 0.76; P < 0.001) was observed between the relative abundance of lactose determined through H 1 NMR and lactose content predicted through mid-infrared spectroscopy, documenting a good level of agreement between the 2 analytical techniques.
Similar to acetate and lecithin, galactose exhibited positive correlations with both protein (r = 0.50; P < 0.001) and CN contents (r = 0.53; P < 0.001). Li et al. (2017) suggested that different polymorphisms of the UDP-galactose-4-epimerase bovine gene, which codes for an enzyme regulating galactose metabolism, can be used as a molecular marker implicated in milk protein concentration in dairy cattle. This may corroborate the association between milk galactose and milk protein contents observed in the present study.
Milk samples with greater SCC were characterized by lower levels of sarcosine and dimethylamine (r = −0.22; P < 0.001), lactose (r = −0.27; P < 0.001), and mannose (r = −0.21; P < 0.001). By contrast, valine, lactate, and O-acetilcarnitine had weak positive correlations with SCS (r = 0.27, 0.24, and 0.23, respectively; P < 0.01). The decreased content of lactose in milk with high SCC is well documented in literature and can be ascribed to (1) an increased transfer of lactose from milk to blood, to keep osmotic pressure constant, and (2) an impaired biosynthesis of lactose at the mammary gland level (Costa et al., 2019). Furthermore, in accordance with the results of the present study, Luangwilai et al. (2021) indicated that increased levels of valine and lactate (together with acetate and phenylalanine) are associated with greater milk SCC and mastitis events.

Effect of Grazing on Milk Metabolites
Results from the ANOVA of milk metabolites are summarized in Table 3. The model adjusted for the fixed effects of breed, parity, DIM (results not shown), and period of sampling. Out of 40 detected metabolites, the period of sampling significantly affected 29 compounds (P < 0.05). Among the latter, and in terms of effect size, 8 metabolites showed negligible variation, 12 small variation, 3 medium variation, and 6 large variation. Figure 3 shows the log2(fold change) of the area under the peaks of the identified metabolites. Out of the 40 detected compounds, 14 metabolites were significantly lower in stable grazing conditions (i.e., 1 wk after the beginning of the grazing period), including fumarate, glucose, galactose-1-phosphate, orotate, glutamate, cis-aconitate, valine, phosphocreatine+creatine, 2-oxoglutarate, N-acethylcarbohydrates, 2 unknown species of aldehyde, and 1 unknown metabolite (P < 0.05). In contrast, grazing resulted in a significant increase of 15 milk metabolites, comprising methanol, hippurate, acetate, 3-hydroxybuthirate, sarcosine+dimethylamine, N-dimethylglycine, dimethylsulfone, butyrate, lactate, carnitine, galactose, acetone, lecithin, succinate, and 1 unknown species of phenyl (P < 0.05).
In accordance with results of the current study, previous authors have observed that milk from pasture-based systems has significantly greater amount of hippuric acid, which has been associated with the presence of caffeoylquinic compounds in fresh forages (Besle et al., 2010;Carpio et al., 2010). This variation was confirmed for both cow (Besle et al., 2010;O'Callaghan et al., 2018) and goat milk (Carpio et al., 2010). Therefore, based on findings of the present study and in accordance with previous literature, it seems appropriate to propose hippurate as a robust marker for milk produced under pasture conditions. Increased levels of 3-hydroxybutyrate in milk obtained under pasture conditions may be due to greater fat mobilization in grazing animals as a result of lower energy intake and higher motility (Benedet et al., 2019). Ashokan et al. (2021) reported that grazing significantly affected 35 milk metabolites. Among them, tyrosyl-threonine, histidinyl-cysteine, 1-methyladenine, cysteine, and selenocysteine showed a sharp increase in milk of grazing cows. Tenori et al. (2018) reported significant variations on bulk milk metabolites depending on different feeding regimens. In particular, the inclusion of silages in animal feed rations caused lower levels of choline, methionine, and hippurate, and greater levels of creatinine, lactate, and an unknown compound. Still, a direct comparison with the results of the present study is not always feasible, due to differences in test hypothesis, experimental design (i.e., bulk milk samples and individual milk samples), and analytical conditions (i.e., protocols for sample preparation and metabolite extraction).
Evidence from the present research indicates that 1 H NMR can be exploited to guarantee food origin, but also that some milk metabolites can potentially be used as markers to differentiate among diets. Although costly and therefore challenging to be implemented on a routine basis, this laboratory technique could be proposed as a food authentication technique, in particular for dairy foods sold at a greater marker price because of processing from grazing cows' milk. This is the case, for example, for Asiago d' Allevo Protected Designation of Origin cheese (Segato et al., 2019), which, in some cases, is produced only from grazing animals and marketed with a premium market value. Indeed, it has been demonstrated that consumers were willing to pay more for dairy products with certified practices of improved animal welfare, including pasture access (Olynk and Ortega, 2013). This would encourage the development of a geographical traceability system for authentication of milk origin, whose cost could be, at least partially, covered by the additional market value of dairy foods from pasture-based cows' milk. Such a strategy would increase the breadth, depth, and precision of a geographical traceability system, all of which are desirable characteristics in terms of establishing products' origins and histories (Dalvit et al., 2007), and consumers would benefit in terms of increased confidence when purchasing and paying a higher price for a certified labeled food.

CONCLUSIONS
Results of the present study highlighted that the RF models applied to 1 H NMR spectra were able to distinguish between samples obtained under different farming conditions with an overall accuracy from 58.7 to 93.1%. Moreover, our study suggests that grazing conditions are associated with specific alterations in the metabolic pattern of milk. Indeed, 29 out of the 40 detected metabolites (e.g., fumarate, methanol, hippurate, and 3-hydroxybutyrate) had significantly different levels of expression depending on the sampling period. These findings represent the first step toward the development of more specific trials to screen for milk metabolites, which have been shown to vary significantly according to grazing, in the view of future traceability applications. Particularly in this perspective, results of  the present study may contribute to the development of rapid cow-side tests to discriminate milk from indoor and pasture-grazed animals.