An Untargeted Metabolomic Comparison of Milk Composition from Sheep Kept Under Different Grazing Systems

: This study aimed to evaluate the effects of different feedings on main traits and polar and semi-polar metabolite profiles of ovine milk. The milk metabolome of two groups of Sarda sheep kept under different grazing systems were analyzed by gas chromatography coupled with mass spectrometry (GC-MS) and multivariate statistical analysis (MVA). The results of discriminant analysis indicated that the two groups showed a different metabolite profile, i.e. ES milk samples of sheep kept under Grazing System 1 (GS1) were richer in nucleosides, inositols, hippuric acid, and organic acids, while milk of sheep under Grazing System 2 (GS2) showed higher levels of phosphate. Statistical analysis of milk main traits indicates that fat content was significantly higher in GS1 samples while milk from GS2 sheep had more urea, trans -vaccenic acid, and rumenic acid. MVA studies of the associations between milk main traits and metabolite profile indicated that the latter reflects primarily the long chain fatty acid content, the somatic cell count (SCC), and lactose levels. All together, these results demonstrated that an integrated holistic approach could be applied to deepen knowledge about the effects of feeding on sheep’s milk composition. 1 polyphite natural pasture: Ferula communis, Thapsia garganica, Asphodelus ramosus, Pteridium aquilinum subsp. Aquilinum, Urginea maritima, and Trifolium repens; 2 selected pasture for direct grazing and to obtain forage: 70% Avena sativa, 15% Hordeum vulgare and 15% Trifolium incarnatum; 3 DMI = dry matter intake.


Introduction
Compared to the rest of the world, where bovine milk represents the most common dairy source, the Mediterranean area is characterized by a significant number of sheep flocks designated to milk rather than meat production [1]. Sardinia, Italy, is one of the leading Mediterranean regions in the production of sheep dairy products. Sheep breeding and pecorino cheese manufacturing have strong economic and cultural links with this region. Though most of the herds in Sardinia are still under traditional grazing on natural pastures of native flora [1], pasture with selected plants and diet supplements have been introduced to ameliorate milk quality, the well-being of sheep, and for environmental and economical sustainability [2]. The effects of these feeding systems on milk quality should be studied. It is widely recognized that diet and season play a major role in modulating chemical composition of ruminant's milk. It has been assessed that the main components of sheep's milk appear to be less influenced by the type of farming system (indoors vs. outdoors) than by the feeding system (pasture alone vs. pasture plus feed supplements) [3], and milk urea measure has been successfully proposed to evaluate the energy contents of a diet [4]. Most of the studies on the effects of feeding systems on sheep milk quality are focused on the production of health beneficial fatty acids (FA), such as rumenic acid (C18:2 cis-9, trans-11) and linolenic acid (C18:3, cis-9, 12,15) [5], as well as on the effects of polyphenols on milk quality and animal performances [2,6]. Differences in ovine milk and cheese composition due to barley and corn supplement to dry pasture have been recently studied by Caprioli et al. [7]. They assessed positive effects on milk retinol and alpha-tocopherol contents and cheese volatiles. Although several researches have studied different aspects of ovine milk main composition [1,[8][9][10][11], there is still little information in literature regarding the metabolite profile of sheep's milk [12][13][14], as most of the works focused on cow's milk. One of the most valuable approaches for the investigation of metabolite profiles in a biological system is metabolomics [15]. Through the application of a gas chromatography coupled with mass spectrometry (GC-MS) based metabolomic approach, it has been possible to study different metabolic pathways linked to quality and other aspects of milk production [16][17][18][19], such as mastitis [12], diet [20], and ketosis [21]. In this study, by a GC-MS metabolomics approach, the milk metabolite profiles of sheep bred in Sardinia under different feeding systems were investigated in association with measured milk main traits such as fat, proteins, caseins, lactose, and urea.

Animals and Pasture
A total of seventy milk samples were collected from Sarda breed sheep. Animals belonged to two different flocks kept under two different grazing systems (GS1 and GS2) located in the same geographical area (coordinates 40°7'48.5" N, 8°40'42.9", E and 40°11'15.1" N, 8°41'38.1" E for GS1 and GS2, respectively) and underwent two different feeding systems (Table 1). In particular, one flock (GS2) was housed from 6 p.m. to 5 a.m. and during the day (10 a.m. to 14 p.m.), it had access to natural pasture while was allowed to graze in a selected pasture (Avena sativa 70%, Hordeum vulgare 15%, and Trifolium incarnatum 15%) only for 1-2 h in the afternoon. GS2 flock had ad libitum access to dry forage (from the selected pasture) during night. During the day, the GS1 flock had unlimited access to natural pasture. The natural pasture available for both flocks was composed by Ferula communis, Thapsia garganica, Asphodelus ramosus, Pteridium aquilinum subsp. Aquilinum, Urginea maritima, and Trifolium repens. Twice a day, both flocks received barley and corn grains (50/50 w/w) before milking (500 and 200 g/head/d for GS1 and GS2 groups, respectively), with an average dry matter intake (DMI) of 438 and 175 g/d/head for GS1 and GS2 groups, respectively ( Table 1). The feeding systems were those usually adopted by the farmers, so they did not represent an additional management and manual labor or a stress for the animals. GS1 flock was composed of approximately 100 ewes and the GS2 flock of 120 ewes. The weight of animals was estimated as 43 ± 6 for GS1, and 40 ± 5 for GS2 and a body condition scoring (BCS) of 3.3 ± 0.2 for GS1 and 3.0 ± 0.2 for GS2. Table 1. Description of the two diet systems.

Samples Collection
The morning individual milk samples were collected in spring 2019 from 70 pluriparous Sarda breed sheep in their late lactation period. Milk samples (n = 37 for GS1 and n = 33 for GS2) were collected in sterilized tubes through manual milking. Each milk sample was added with the preservative bronopol (2-bromo-2-nitro-1,3-propanediol) and divided in two aliquots. One was analyzed by the Sardinian Regional Animal Farmers Association (Associazione Regionale Allevatori, ARA, Sardegna) milk laboratory, and the other was stored at −20 °C for GC-MS analysis, which was performed within two days from collection. The milk samples were provided by the ARAS within the International Committee for Animal Recording (ICAR) program.
Milk yield was registered for each ewe at the morning milking and the total amount of milk produced by the two groups was recorded.

GC-MS Analysis
Thawed milk (1 mL) was subjected to ultrasounds for 15 min, and 0.1 mL were extracted following the Folch procedure [23] using 0.375 mL of a methanol and chloroform mixture (2/1, v/v). Samples were vortexed every 15 min up to 1 h, when 0.38 mL of chloroform and 0.09 mL of aqueous KCl 0.2 M solution were subsequently added. Samples were vortexed again and centrifuged for 15 min at 15294 g (Eppendorf 5810R, Milan, Italy). Two hundred µL of the hydrophilic supernatant were dried in glass vials using a nitrogen stream and then derivatized using 0.05 mL of methoxamine chloride dissolved in pyridine at 10 mg/mL, homogenized for 20 s, and kept at room temperature for 17 h. Then, 100 µL of MSTFA were added and samples were vortexed. After 1 h, 600 µL of hexane were added and samples homogenized again before GC-MS analysis. A Hewlett Packard 6850 gas chromatograph, a 5973 mass selective detector, and a 7683B series injector (Agilent Technologies, Palo Alto, CA, USA) were used to analyze samples using helium as the carrier gas at 1.0 mL/min flow. One µL of sample was injected in split-less mode and separated using a 30 m × 0.25 mm × 0.25 µm DB-5MS column (Agilent Technologies, Palo Alto, CA, USA). The temperatures for the inlet, interface, and ion source were 250 °C, 250 °C, and 230 °C, respectively. The oven temperature was programmed to increase from 50 to 230 °C by 5 °C/min over 36 min and kept at this temperature for 2 min. The mass range was set between 50 and 550 m/z using and electron voltage of 70. MSD ChemStation software (Agilent Technologies, Santa Clara, CA, USA) was used to elaborate the data. The GC-MS spectra deconvolution was performed by the AMDIS tool in the NIST08 library. Retention time and relative mass spectra of each compound were compared with the related analytical standards with the aim to identify those compounds previously recognized by consulting the NIST14 library of the National Institute of Standards and Technology (Gaithersburg, MD, USA), and the library developed at the Max Planck Institute of Golm (http://gmd.mpimp-golm.mpg.de/). Chromatograms in the AIA format were then uploaded to the XCMS Online platform [24]. The output of XCMS consisted in a 70 X 1230 matrix where each variable corresponded to the intensity value of an m/z ion at a specific retention time value. GC-MS features that annotated to bronopol with a retention time of 21.84 min were excluded.

Statistical Data Analysis
To obtain descriptive information on the data and to calculate the means and their standard deviations, we performed univariate analysis with the software OriginPro2016 (OriginLab, Northampton, MA, USA). When the two groups of samples were compared, the null hypothesis ("the means are not significantly different among the two sets of samples") was tested by using the unpaired samples two-tailed Student t-test. The output of XCMS pipeline was submitted to MVA, and analysis were performed as implemented in SIMCA-P+ software (version 14.1, Umetrics, Umeå, Sweden). Variables were mean centered and scaled to unit variance. At first, for sample distribution overview we performed a Principal Component Analysis (PCA) [25]. Subsequently, to identify discriminant metabolites between the two milk typologies, we performed an Orthogonal Partial Least-Squares-Discriminant Analysis (OPLS-DA). The obtained variable importance in projection (VIP) scores summarize the contribution of each variable to the model. In this work, GC-MS features showing VIP values greater than 1 underwent manually to a supervised analytes identification. A metabolite was considered significant only when at least two of its most abundant mass fragments and a retention time deviation <0.05 min were found in the list of VIP > 1. For quantification purposes for each metabolite, we considered the intensity of the most abundant mass fragment. Milk metabolite profiles were correlated to the measured main traits by a single-Y Partial Least-Squares (PLS) regression. The quality of the models was evaluated based on the cumulative parameters R 2 Y and Q 2 Y estimated by the default leave-1/7th-out cross validation in the corresponding PLS-DA models, as implemented in SIMCA-P+ program. The models were considered significant only when the difference between R 2 Y and Q 2 Y was <0.50 [26].

Main Traits
The milk yield was higher in the sheep kept under grazing system 2 (GS2 flock, 784 ± 73 g/d) compared to (GS1 flock, 581 ± 64 g/d). The milk content of proteins, caseins, fat, lactose, FA, and other milk traits (expressed as mean ± SD) are reported in Table 2.
Milk of GS1 sheep showed a statistically higher fat content when compared to the GS2 group (5.9 ± 1.0 and 4.8 ± 1.2 g/100 mL, p < 0.001). An average milk fat level of 6-7 g/100 mL was reported for the Sarda breed [9]. It is well known that the milk fat content has a high variability, especially when compared with proteins and lactose levels, being strongly influenced by the feeding composition [3]. In the GS1 milk samples, consistently with their higher levels of fat, saturated and unsaturated long chain (C14-C18) FA linked to diet were found at higher levels, with the exception of n-3 FA. Contents of the de novo synthetized FA produced by the mammary gland (C4:0, C6:0, C8:0, C:10, C12:0) were found similar to those previously reported in literature for the Sarda breed sheep [27]. Differently from all other de novo FA, the concentration of C4:0 was statistically higher in GS1 samples. A different behavior of C4:0 with respect to the other de novo synthetized FA was already observed in sheep's [27] and cow's [28] milk. Genetic studies confirmed the independence of C4:0 from de novo mammary FA synthesis [29].
Among milk FA, due to their health benefits [30], a particular attention has been paid to rumenic acid, trans-vaccenic acid (C18:1 trans-11), and linolenic acid. Levels of trans-vaccenic acid (3.7 ± 0.5 vs. 3.0 ± 0.7 g/100 g of fat in the GS2 and GS1 groups, respectively) and rumenic acid (2.0 ± 0.2 vs. 1.6 ± 0.3 g/100 g of fat in the GS2 and GS1 groups, respectively) were significantly higher in GS2 milk. Milk products are the only food that naturally contain trans-vaccenic acid and rumenic acid produced by the ruminal biohydrogenation and the mammary Δ-9 desaturation pathways [1]. The pasture plays a pivotal role in increasing milk levels of trans-vaccenic acid and rumenic acid [5]. The preferred precursor of trans-vaccenic acid and of rumenic acid is the dietary linoleic acid. Oats, one of the plants in the selected pasture of GS2, is rich in linoleic acid and this fact can explain the higher content of trans-vaccenic acid and rumenic acid in the GS2 group.  A parameter strictly connected with the quality of ewe's diet is the milk urea content, which is representative of urea levels of blood and other body fluids. Milk urea is a breakdown product of proteins and thus is considered a useful management tool to evaluate the metabolism and intake of proteins, minimizing potential negative effects on animals [4]. High milk urea content is associated with a worsening of animal welfare condition and different ruminant's pathologies, such as fertility and subclinical mastitis [31]. High urea content has been related both to the characteristics of the pasture with young sward poor in fiber and rich in highly fermentable proteins and/or non-protein nitrogen. High levels of dietary proteins can cause animal health and reproduction problems as well as the impairment of the milk characteristics and technological properties [4,32,33]. Sheep's milk of GS1 had significantly less urea than GS2 (24 ± 5 vs. 41 ± 8 mg/100 mL). However, both values are within the reported ovine milk urea level ranges [2,9], although some GS2 samples reached values of 60 mg/100 mL, which is a clear sign of unbalanced nutrition [2], while GS1 had values comparable with those found in summer season [9] when pasture is dry and poor. The higher urea content in GS2 can be due to a higher soluble protein content in the selected pasture and lower energy intake due to a probably insufficient grains ration.

Metabolomics
Forty-six polar and semi polar metabolites were identified and reported in Table S1. In Figure 1, the GC-MS chromatograms of milk samples from sheep kept under the two different grazing systems (GS1 and GS2) are shown. The former exhibits, in the range of 33-36 min, besides peaks annotated as lactose, a different chromatographic pattern, when compared to GS2, attributed to disaccharides such maltose and others. This fact is probably linked to a greater consumption of natural polyphite meadow, known to be rich in fiber, and to a higher portion of grains rich in starch in the GS1 feeding system (Table 1). To have an overview of samples similarities and dissimilarities we performed a PCA of GC-MS features. As shown in the score plot reported in Figure S1, milk samples clustered in two groups and GS2 milk samples were more scattered. Clustering indicated that the different feeding systems influenced the metabolite milk profile. In order to find discriminating metabolites between the two groups of sheep's milk samples, we performed a pair-wise OPLS-DA. The discriminant analysis correctly classified milk samples. Score plot is shown in Figure 2, and we report the metabolites that mostly contributed to sample classification in Table 3.  Table 3. Table 3. Pair-wise OPLS-DA a discriminant metabolites between milk from sheep kept under grazing system 1 (GS1) and grazing system 2 (GS2). Further information can be obtained by associating the whole metabolite profile to the main traits, this approach can indicate whether modifications of a single trait are accompanied by modification of the whole metabolite profile. Results of the PLS analysis reported in Table 4 indicate that the metabolite profile is correlated to the content of the long chain FA, among which rumenic acid and trans-vaccenic acid, the most studied indexes of different pastures in the diet [5]. No n-3 FA, neither de novo synthetized FA were found to be significantly linked to the metabolite profiles. Correlations between metabolites and SCC and lactose milk content were already observed [14,16].

Discussion
GS1 milk samples showed higher levels of maltose, inosine, and uridine (Table 3). An increase of maltose and purine/pyrimidine cycle metabolites was observed in the rumen and milk of dairy cows fed with an increasing proportion of diet cereal grains [34][35][36], in agreement, as reported in Table 1, ewes in the GS1 received more concentrate than the GS2 group. Inosine and uridine belong to the non-protein-nitrogen (NPN) fraction of milk. These nucleosides are bioactive components with health beneficial effects [37]. Milk of sheep under GS1 showed higher levels of inositols (myo-inositol and scyllo-inositol, Table 3). Free inositols occur naturally in the environment as a bioactive component of living cells. Milk is a significant source of myo-inositol and for its role in neonatal nutrition, this polyol is often added to infant formulas to prevent a potential deficiency during early neonatal growth stage [38]. As well as urea, myo-inositol is osmolyte and the lower level of this polyol in GS2 can counterbalance their high urea content. Among the discriminant metabolites listed in Table 3, gluconic acid was found upregulated in the GS1 milk samples. Gluconic acid is the product of the glucose oxidation [39]. It has been demonstrated that in the rumen epithelium, butyrate (C4:0) and glucose oxidation mechanisms compete [40], and in the light of this observation we can hypothesize that the GS1 animals with a diet richer in concentrate and therefore in starch have more glucose available as substrate for the oxidizing enzymes, therefore saving more butyrate which is significantly higher in milk samples of the GS1 animals compared to GS2 ( Table 2).
As shown in Table 3, hippuric acid, malic acid, and succinic acid were found in higher levels in GS1 milk samples. The content of hippuric acid in milk has been linked to the presence of polyphenols in forages [41]. Greater levels of hippuric acid in milk from pasture-based feeding systems have previously been reported in milk from goats and cows [19,36,42] and this molecule has been proposed as biomarker of pasture-derived milk [36,42]. In accordance, we found higher levels of hippuric acid in the milk of sheep allowed to graze on natural pasture for a longer time and thus used as biomarker for ovine milk. Malic acid and succinic acid are intermediates in the succinate and propionate pathway of ruminal bacteria for carbohydrate fermentation. Their higher levels in GS1 group could be linked to a greater consumption of natural polyphite meadow rich in fiber and grains when compared to GS2 group (see Table 1).
Compared to the GS1 samples reported in Table 3, milk of the GS2 group showed higher phosphorus content. A relationship between phosphorus and casein milk content is expected since this mineral is bound to casein micelles. Besides protein, lactose content can also affect the phosphorus content of milk since inorganic phosphorus is generated during the formation of lactose in the mammary gland and subsequently secreted into milk [43]. Although significant relationships between milk phosphorus and both protein and lactose contents were found [44], a considerable portion of the phosphorus content of milk remained unexplained, implying that other factors may affect milk phosphorus levels [45]. Feeding management is another important factor that can influence mineral concentration in milk. One source of organic phosphorus is phytate, which contributes to the majority of phosphorus in grains [46]. Ruminants can utilize phosphorus from phytate because ruminal microorganisms are capable of synthesizing phytase enzyme, which can release a phosphate group from the phytate [47]. Ruminal phytate hydrolysis is influenced by the type of grain, processing of feed ingredients, and supplemental exogenous phytase enzyme [48,49]. Lastly, as for the highly discriminant not-annotated metabolites (Unk#1 and Unk#2) characterizing GS2 milk samples, further studies will be carried out by GC-TOFMS to elucidate their chemical structure.

Conclusions
In summary, grazing system 1 led to a higher milk yield, more trans-vaccenic acid, and rumenic acid. Although accompanied by high urea levels, grazing system 2 was associated with higher fat and presence of metabolites linked to the grain supplements.
Our results suggest that the GC-MS metabolite profiles of ovine milk well reflect variations in diet with respect to pasture and concentrate supplementation and is strongly associate with the content of health beneficial FA. Further work with an ampler data set over a longer period of time will be carried out by this approach with the aim of establishing a model to classify milk and dairy products from sheep grazing natural pasture.

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