Milk metabolome reveals pyrimidine and its degradation products as the discriminant markers of different corn silage-based nutritional strategies

The purpose of this study was to evaluate the effect of 6 different feeding systems (based on corn silage as the main ingredient) on the chemical composition of milk and to highlight the potential of untargeted metabolomics to find discriminant marker compounds of different nutritional strategies. Interestingly, the multivariate statistical analysis discriminated milk samples mainly according to the high-moisture ear corn (HMC) included in the diet formulation. Overall, the most discriminant compounds, identified as a function of the HMC, belonged to AA (10 compounds), peptides (71 compounds), pyrimidines (38 compounds), purines (15 compounds), and pyridines (14 compounds). The discriminant milk metabolites were found to significantly explain the metabolic pathways of pyrimidines and vitamin B 6 . Interestingly, pathway analyses revealed that the inclusion of HMC in the diet formulation strongly affected the pyrimidine metabolism in milk, deter-mining a significant up-accumulation of pyrimidine degradation products, such as 3-ureidopropionic acid, 3-ureidoisobutyric acid, and 3-aminoisobutyric acid. Also, some pyrimidine intermediates (such as l-aspartic acid, N-carbamoyl-l-aspartic acid, and orotic acid) were found to possess a high discrimination degree. Additionally, our findings suggested that the inclusion of alfalfa silage in the diet formulation was potentially correlated with the vitamin B 6 metabolism in milk, being 4-pyridoxic acid (a pyridoxal phosphate degradation product) the most significant and up-accumulated compound. Taken together, the accumulation trends of different marker compounds revealed that both pyrimidine intermediates and degradation products are potential marker compounds of HMC-based diets, likely involving a complex metabolism of microbial nitrogen based on total splanchnic fluxes from the rumen to mammary gland in dairy cows. Also, our findings highlight the potential of untargeted metabolomics in both foodomics and foodomics-based studies dairy products.


INTRODUCTION
In Northern Italy, corn silage represents up to 90% of the forage component in the lactating cow diet, because of the soil fertility, favorable climate for corn silage, and its high DM yield potential per hectare (Borreani et al., 2013;Bellingeri et al., 2019). Accordingly, most dairy farms in this area adopted nutritional strategies based on TMR, in which corn silage is the main ingredient of the dairy ration, in combination or not with other grass, legume, or grass-legume silages. The complexity of the role of nutrition in the dairy farm can be exploited deepening the study on the entire transformation process of nutrients in milk but adopting different approaches. Looking at the quantitative relationships amount of input used versus the produced outputs, a particular focus can be oriented to economics and efficiency targets to improve profitability and sustainability of production managerial choices (Atzori et al., 2013;Bellingeri et al., 2020;Atzori et al., 2021;Gallo et al., 2022). Otherwise, focusing on qualitative aspects of relationships among diet, feces, and milk components is particularly useful to improve knowledge of the biological complexity and the most relevant biomarkers of efficient nutrient use.
Over the past few years, the application of metabolomics in milk and other dairy products has grown, particularly concerning health and economic benefits (Zhu et al., 2021). Metabolomics represents the study of low molecular weight organic molecules (<2,000 Da) characterizing a target cell, tissue, organ, or organism. Metabolites involved in several chemical processes can Milk metabolome reveals pyrimidine and its degradation products as the discriminant markers of different corn silage-based nutritional strategies G. Rocchetti, 1 * F. Ghilardelli, 1 E. Carboni, 1 A. S. Atzori, 2 F. Masoero, 1 and A. Gallo 1 be measured by large-scale and high-throughput techniques, including nuclear magnetic resonance and highresolution mass spectrometry (HRMS; Rocchetti and O'Callaghan, 2021). Therefore, the recent application of targeted and untargeted metabolomics to the field of dairy science allowed to select some biomarker compounds related to animal health, production, authentication, and contributors to milk techno-functional properties.
Milk is a complex biological fluid produced by the mammary glands of mammals. As the end products of gene expression in mammary epithelial cells, peripheral blood, and microorganisms, the final composition of milk will be different according to both inherent and external factors (Zhu et al., 2021). However, the feeding system is recognized as the main factor affecting the metabolic processes from nutrient intake to mammary gland release (Lamanna et al., 2011). Among the most impacting inherent factors, it's important to list the cow breeds and the lactation period (Xu et al., 2018(Xu et al., , 2020, whereas seasonal variation, feeding systems, herd management, geographical origin, cow's health condition, processing, storage, and technological conditions represent the most studied external factors, potentially affecting the milk metabolome. Regarding the external factors, the effect of different feeding systems on the chemical profile of milk represents one of the most recently studied topics (O'Callaghan et al., 2018;Magan et al., 2019;Rocchetti et al., 2020Bellassi et al., 2021;Yanibada et al., 2021). Also, the combination of botanical origin and conservation method of forage is known to modify the milk composition, including fatty acids, fat-soluble vitamins, N-compounds, and organic acids (O'Callaghan et al., 2016). Recently, Ran et al. (2021) investigated the effect of partial substitution of corn silage with sweet sorghum silage in the diet of lactating cows, showing that the substitution altered the relative abundances of some predominant rumen bacteria. However, these changes had little effect on ruminal fermentation and milk yield and composition. Also, Boivin et al. (2013) evaluated the effect of grain and forage fractions of corn silage on milk production and composition in dairy cows, highlighting the importance of energy and nutrients supplied by the grain fraction of corn silage to support milk yield and fat composition.
In a previous work by Sun et al. (2015), a metabolomic approach (based on GC-TOF-MS), was exploited to find potential common and exclusive biomarkers in 4 different biofluids such as rumen fluid, milk, serum, and urine, and considering nutritional strategies based on corn stover versus alfalfa hay. In particular, the authors demonstrated that some AA, such as Phe and Tyr, may be used more in the liver in the alfalfa-fed cows than in the animals fed corn stover. An HRMS approach was also used by Craige Trenerry et al. (2013) for the targeted characterization of milk lipidic fraction from dairy cows fed different diets, while Riuzzi et al. (2021) recently performed the authentication of foragebased milk by mid-level data fusion of direct analysis in real-time-HRMS signatures, providing 25 potential biomarker compounds in milk.
To the best of our knowledge, there is still limited and fragmented information on the application of milk metabolomics to evaluate potential correlations with different feeding strategies based on corn silage as the main ingredient of the ration. Therefore, in this work, we evaluated the potential of untargeted metabolomics, based on ultra-high-pressure liquid chromatography (UHPLC) coupled with HRMS, to discriminate the chemical profile of bulk milk samples according to different nutritional strategies. Milk metabolomics may be used to capture most of the small molecular compounds in milk (as an easily collected and stored biofluid) and to identify those that significantly changed among different dietary treatments. This approach can discern valid biomarkers and discriminant metabolic pathways to better understand the complex mechanisms involving the feeding system and milk composition.

MATERIALS AND METHODS
No human or animal subjects were used, so this analysis did not require approval by an Institutional Animal Care and Use Committee or Institutional Review Board.

Collection of Milk Samples and Classification Criteria
In a previous work by our research group (Gallo et al., 2022), an observational study was carried out to verify the influence of different nutritional corn silagebased strategies on efficient use of dietary nutrients, fecal fermentation profile, and profitability in a cohort of intensive dairy farms (i.e., a total of 66 commercial dairy farms of a specific area of Italy, namely Po valley). Overall, farms had a single diet for lactating cows, a single diet for dry cows, and 2 diets for heifers (from weaning to first calving), except for 2 dairy farms that divided the lactating dairy cows into 2 groups, being early-(<150 DIM) and late-lactation (>150 DIM) cows. Interestingly, an unsupervised hierarchical cluster analysis was found to successfully discriminate the farm feeding choices in 6 main nutritional strategies, based on high use of (1) high-moisture ear corn (HMC) and legume silage (cluster 1), (2) compound feed (cluster 2), (3) corn and soy meals (cluster 3), (4) HMC and soy meal (cluster 4), (5) cornmeal and protein compound feed (cluster 5), and (6) HMC and protein compound feed strategies (cluster 6).
The raw bulk tank milk samples (500 mL) were taken between January and June 2018 (Gallo et al., 2022), with this period being the best of the year in terms of highest productive performances and less influence due to seasonal changes on herd composition, milk quality, and silage bunker. The milk samples were collected before the day of the visit by the milk tank of 36 farms containing both afternoon and successive morning milking and then stored at 4°C. The 36 farms represented a random subsample of the cluster identified from the original 66 farms described by Gallo et al. (2022) and more details regarding the clusters of considered nutritional strategies are reported in Supplemental Table S1 (https: / / www .doi .org/ 10 .17632/ ggpyyzkd4p .1; Rocchetti, 2022a). These farms can be considered representative of the intensive dairy farm system in Italy, raising Holstein-Friesian cows, housed in freestall barns, fed TMR, and with no access to pasture.

Extraction of Milk Metabolites
The extraction of milk metabolites was carried out as previously reported in . The collected bulk milk samples (a total of 108 samples, considering 3 biological replicates for each milk sample) were skimmed by centrifugation at 4,500 × g for 10 min at 4°C and then frozen at −80°C for further processing. Thereafter, samples were thawed and thoroughly vortex mixed. An aliquot of 2 mL was extracted with 14 mL of acetonitrile (liquid chromatography-MS grade, Sigma-Aldrich) added with 3% formic acid (vol/vol), mixed by vortexing for 2 min, and processed with ultrasounds for 5 min. Milk samples were then centrifuged at 12,000 × g for 15 min at 4°C to remove large biomolecules (such as proteins). The supernatants were filtered through 0.22-μm cellulose syringe filters in UHPLC vials until further untargeted metabolomic profiling.

Metabolomics Profiling Based on HRMS
The untargeted UHPLC-HRMS analysis was done using a Q Exactive Focus Hybrid Quadrupole-Orbitrap Mass Spectrometer (Thermo Scientific) coupled to a Vanquish UHPLC pump and equipped with heated electrospray ionization-II probe (Thermo Scientific; . The chromatographic separation was based on a water-acetonitrile (both liquid chromatography-MS grade, from Sigma-Aldrich) gradient elu-tion (6-94% acetonitrile in 35 min), using 0.1% formic acid as phase modifier, and using an Agilent Zorbax Eclipse Plus C18 column (50 × 2.1 mm, 1.8 μm). The flow rate was 200 μL/min, and the injection volume was 6 μL, using a full scan MS analysis (100-1,200 m/z range) with a positive ionization mode and a mass resolution of 70,000 at m/z 200. The automatic gain control target and the maximum injection time were 1 × 10 6 and 200 ms, respectively. Also, pooled quality control samples were randomly injected and analyzed in a data-dependent (Top N = 3) MS/MS mode, with full scan mass resolution reduced to 17,500 at m/z 200, with an AGC target value of 1 × 10 5 , maximum injection time of 100 ms, and isolation window of 1.0 m/z, respectively. The Top N ions were fragmented using a stepped Normalized Collisional Energy (i.e., 10, 20, 40 eV). The heated electrospray ionization parameters were the following: sheath gas flow 40 arbitrary units, auxiliary gas flow 20 arbitrary units, spray voltage 3.5 kV, capillary temperature 320°C. Before data collection, the mass spectrometer was calibrated using Pierce positive ion calibration solution (Thermo Fisher Scientific). The collected .RAW data were further processed using the MS-DIAL software (version 4.70; http: / / prime .psc .riken .jp/ compms/ msdial/ main .html). Automatic peak finding, LOWESS normalization, and annotation via spectral matching against the database Bovine Metabolome Database, Phenol-Explorer (considering the intestinal metabolites), and Mass Bank of North America were performed. The mass range 100-1,200 m/z was searched for features with a minimum peak height of 10,000 counts per second. The MS and MS/MS tolerance for peak centroiding was set to 0.01 and 0.05 Da, respectively. Retention time information was excluded from the calculation of the total score. Accurate mass tolerance for identification was 0.01 Da for MS and 0.05 Da for MS/MS. The identification step was based on mass accuracy, isotopic pattern, and spectral matching. The selected identification criteria were used to calculate a total identification score. The total identification score cut off was >50%, considering the most common heated electrospray ionization + adducts. Finally, gap filling (using the peak finder algorithm) was used to fill in missing peaks, considering a 5 ppm tolerance for m/z values.

Multivariate Statistics and Pathway Analysis
The metabolomics-based data were elaborated for multivariate statistical modeling using 3 different software, namely Mass Profiler Professional (Agilent Technologies), MetaboAnalyst, and SIMCA 13 (Umetrics), as reported in previous works (Rocchetti et al., 2020. After data normalization, both unsupervised and supervised multivariate statistics were carried out. The unsupervised approach was based on hierarchical cluster analysis, while the orthogonal projections to latent structures discriminant analysis (OPLS-DA) was used as a supervised tool. Additionally, the OPLS-DA model validation parameters (goodness-of-fit R 2 Y together with goodness-of-prediction Q 2 Y) were inspected, considering a Q 2 Y prediction ability of >0.5 as the acceptability threshold. Thereafter, the OPLS-DA model produced was inspected for outliers and permutation testing (n > 100) was used to exclude model over-fitting. The importance of each milk metabolite for discrimination purposes (i.e., when considering the different nutritional clusters) was then calculated according to the variable selection method variable importance in projection (VIP), using as the minimum significant threshold a score >0.8. Thereafter, volcano plots were produced for the pairwise comparisons of clusters 2, 3, and 5 versus HMC-based group (combining clusters 1, 4, and 6), coupling FC analysis (cutoff value >1.2) and ANOVA (P-value < 0.05). Finally, the online tool MetaboAnalyst was used to inspect both the most important pathways represented by the metabolites annotated (using as pathway library: Bos taurus, Kyoto Encyclopedia of Genes and Genomes, KEGG), and then to provide a metabolite set enrichment analysis based on the discriminant and significant class of metabolites and metabolic pathways outlined by the multivariate statistical approach used.

Discrimination of Milk Samples According to the Different Feeding Strategies
In this work, the screening based on the untargeted metabolomics approach allowed for the putative annotation of 2,954 mass features. Therefore, to reduce data complexity and to provide a data set of compounds to be used for the following multivariate statistical analyses, 3 liquid chromatography-MS comprehensive databases (when considering the electrospray ionization positive polarity) were exploited, namely Bovine Metabolome Database, Phenol-Explorer, and Mass Bank of North America. Following this strategy as a first step, the isotopic identification of 697 milk metabolites was achieved, also recording 88 structurally confirmed compounds, based on the tandem MS information provided by quality control samples (i.e., annotated by tandem MS). Overall, among the compounds identified, we found a great abundance of AA, peptides, carboxylic acids and derivatives, feed-related compounds (such as derivatives of hippuric acid and phenolic metabolites), nucleic acids derivatives, and pyridine derivatives. A detailed list reporting all the compounds annotated according to this level 2 (i.e., putative identification and structural confirmation of the most abundant mass features) can be found in Supplemental Table S2 (https: / / www .doi .org/ 10 .17632/ 9vtb2dwtkj .1; Rocchetti, 2022b).
The next step was the application of a supervised statistical analysis, namely OPLS-DA to maximize the explained covariance on the first latent variables (a small subset of hidden variables in the X-data to predict the response variables), while the remaining latent variables capture variance in the predictors which is orthogonal (i.e., statistically uncorrelated with the response variables). The corresponding score plot, built considering as classification parameters the 6 clusters related to the nutritional strategies (i.e., different feeding systems; Gallo et al., 2022), is reported in Figure 1. The OPLS-DA score plot was characterized by significant goodness parameters, being the R 2 Y (goodness of fitting) = 0.953 and the Q 2 (goodness-of-prediction) = 0.606. In addition, the score plot revealed interesting information paving the way for additional considerations. First, milk samples were separated into 3 different groups, namely cluster 5 (labeled as protein compound feed strategy), cluster 2 (labeled as compound feed strategy), and cluster 3 (labeled as corn and soy meals strategy). Interestingly, milk samples belonging to cluster 1 (labeled as HMC and legume silage), 4 (labeled as HMC and soy meal), and 6 (labeled as HMC and protein compound feed strategies) were found to cluster together thus revealing a very similar metabolomic profile. Therefore, according to the observational study previously published by our research group (Gallo et al., 2022), the OPLS-DA provided clear evidence of the potential effect of HMC-based nutritional strategies on the metabolomic profile of bulk milk samples under investigation. Accordingly, the next part of this work was based on the pairwise comparisons between HMCbased milk samples (combining clusters 1, 4, and 6) versus the remaining clusters (i.e., 2, 3, and 5) to better understand the most important and discriminant variables driving the hyperspace separation observed in Figure 1.

Discriminant Compounds and Enriched Chemical Classes
The most discriminant variables (i.e., VIP marker compounds) driving the class separation observed in Figure 1 were then extrapolated and reported in Supplemental Table S2. Overall, a total of 322 milk metabolites showed a VIP score >0.8 (i.e., the minimum cutoff for discriminant ability), while 211 compounds (65.5%) were characterized by a VIP score >1 (maximum degree for discriminant ability). The most discriminant and numerically represented classes were peptides and derivatives (71 compounds), followed by pyrimidine derivatives (38 compounds), purine derivatives (15 compounds), pyridine derivatives (14 compounds), and essential or nonessential AA (10 compounds). Accordingly, the enrichment analysis (based on the most represented classes) confirmed that pyrimidine derivatives were the most significantly (P < 0.05) enriched class with the higher number of metabolic intermediates and when compared with the other classes (Figure 2A and 2B). Those milk metabolites presenting the higher VIP score within each most discriminant class were then reported in Table 1. As can be observed, the most of discriminant markers (such as guanosine diphosphate, l-aspartic acid, thymidine, and 4-pyridoxic acid) discriminated HMC-based milk samples from the other 3 clusters under investigation.
Thereafter, to point out both common and exclusive VIP metabolites for the pairwise comparisons against the HMC-based cluster, a Venn diagram was inspected and reported in Figure 3. Overall, the 37.3% of the discriminant metabolites were shared between the different conditions, while the pairwise comparison "HMCbased versus cluster 3" was found to provide the highest number of exclusive discriminant compounds (13.8%), followed by "HMC-based versus cluster 2" (9.6%), and "HMC-based versus cluster 5" (6.2%). Among the most discriminant milk metabolites characterizing the different comparisons, we found N1-acetylspermine (VIP score = 2.70; HMC-based versus cluster 2), glycylaspartate (VIP score = 2.12; HMC-based versus cluster 3), and l-cysteinylglycine disulfide (VIP score = 2.12; HMC-based versus cluster 5). The low number of exclusive discriminant compounds when considering the 3 different comparisons (Figure 3) is coherent with the high impact of the HMC factor previously described. Rocchetti et al.: CORN SILAGE-BASED STRATEGIES ON MILK METABOLOME Figure 1. Orthogonal projection to latent structures discriminant analysis considering the milk samples classified according to the nutritional strategy of dairy cows. Cluster 1: labeled as high-moisture ear corn (HMC) and legume silage strategy, in which HMC was used; cluster 2: labeled as compound feed strategy, in which HMC was not used; cluster 3: labeled as corn and soy meal strategy, in which HMC was not used; cluster 4: labeled as HMC and soy meal strategy; cluster 5: labeled as protein compound feed strategy, in which HMC was not used; cluster 6: labeled as HMC and protein compound feed strategy. Cum = cumulative; t[1] = first latent vector; t[2] = second latent vector.

Pathway Analysis and Significant Metabolites
The multivariate statistical analysis (OPLS-DA modeling) clearly showed a potential correlation between the use of HMC (in substitution of cornmeal) and the metabolomic profile of bulk milk samples. Therefore, to find the most represented metabolic pathways (i.e., those most affected by the HMC inclusion in the ration), 3 different pathway analyses were built, considering the VIP discriminant marker compounds of the comparisons against the cluster HMC. Overall, as shown in Figure 4, 3 metabolic pathways were found to possess a high impact for each comparison against the HMC-based cluster, namely pyrimidine metabolism, vitamin B 6 metabolism, alanine, aspartate, and gluta-mate metabolism. As a general consideration, a greater and significant (P < 0.01) impact of the pyrimidine metabolism for each comparison against the HMCbased cluster could be highlighted. Additionally, the 3 metabolome view maps built according to the VIP metabolites annotated by UHPLC-HRMS revealed that other significant and common metabolic pathways were those involving vitamin B 6 metabolism (i.e., pyridoxal derivatives) and alanine, aspartate, and glutamate metabolism (AA). Going into detail, the comparison "HMC-based versus cluster 2" revealed an impact = 0.75 of the pyrimidine metabolism on the metabolomic differences observed, with the 47.4% of discriminant pyrimidine metabolites involved. The impact of pyrimidine metabolism for this comparison was extremely Rocchetti et al.: CORN SILAGE-BASED STRATEGIES ON MILK METABOLOME Figure 2. Enriched discriminant metabolic pathways (A) and discriminant main chemical classes (B) considering those milk metabolites annotated by ultra-high-pressure liquid chromatography, high-resolution mass spectrometry, and having a variable importance in projection score >0.8, as a function of the high-moisture ear corn. significant, recording a P-value = 4.6 × 10 −11 and a false discovery rate (FDR) = 3.8 × 10 −9 . Another significant finding was represented by the great impact of compounds involved in the vitamin B 6 metabolism; this latter, although represented by a lower number of metabolites (i.e., 6/9) was highly impacting (0.87) with a P-value = 1.6 × 10 −5 and FDR = 6.7 × 10 −4 . Looking at the comparison ""HMC-based versus cluster 3," a major impact (0.91) of the pyrimidine metabolism on the metabolomic differences of milk was observed, recording also extremely significant values, being P-value = 8.78 × 10 −16 and the FDR = 7.38 × 10 −14 . In addition, 22 pyridine derivatives (on a total of 38 compounds) were found to discriminate the HMC-based milk samples when compared with cluster 3-related milk samples. Additionally, other significant and impacting pathways were represented by alanine, aspartate, and glutamate metabolism (P-value < 0.01; impact = 0.46), arginine biosynthesis (P-value < 0.05; impact = 0.46), phenylalanine metabolism (P-value < 0.05; impact = 0.36), and vitamin B 6 metabolism (P-value < 0.05; impact = 0.33). Finally, the comparison "HMC-based versus cluster 5" revealed again similar information, with a significantly higher impact of pyrimidine metabolism (P-value = 3.6 × 10 −9 , FDR = 3.06 × 10 −7 , impact = 0.57), followed by alanine, aspartate, glutamate, and arginine metabolisms (P-value < 0.05, impact = 0.46), and vitamin B 6 metabolism (P-value < 0.05, impact = 0.33). Interestingly, this latter comparison was characterized by a numerically lower number of significant pyrimidine derivatives (15/38) when compared with the others. Considering that pyrimidine metabolism was highlighted as the most impacting and significant for each comparison under investigation, the next step of this work was based on the evaluation of the significant variations fold-change (FC) of each pyrimidine metabolic intermediate to demonstrate the potential correlation between the feeding strategy (i.e., based on HMC as a major factor) and the final metabolomic profile of milk. According to the pyrimidine metabolism reported on KEGG for Bos taurus, the pyrimidine intermediates (33 compounds) were extrapolated from the OPLS-DA modeling and reported in Table 2, considering the 3 different comparisons against clusters 2, 3, and 5. Overall, we focused on different pyrimidine metabolic pathways, involving 3 nucleobases, namely cytosine, uracil, and thymine. The list of compounds included several intermediates, such as ureidosuccinic acid (also known as N-carbamoyl-l-aspartate; VIP score = 1.25), l-aspartic acid (VIP score = 1.24), and orotic acid (VIP score = 1.39). In addition, several nucleotides and nucleosides have been observed with different log2 FC values (both up-and down-accumulation).
Interestingly, some of the most important degradation products of thymidine showed positive and significant (P-value < 0.05) log2 FC values, such as 3-Ureidoisobutyric acid and 3-Aminoisobutyric acid (Table 2). In this regard, the highest accumulation for 3-Aminoisobutyric acid was found for the comparison "HMC-based versus cluster 2" (log2 FC = 0.63), followed by "HMCbased versus cluster 5" (log2 FC = 0.34), and "HMCbased vs. cluster 3" (log2 FC = 0.30). Similarly, an up-accumulation trend of 3-Ureidopropionic acid (the degradation product of both cytosine and uracil) was also observed (Table 2), although with lower log2 FC values. Finally, considering the lower but significant effect of vitamin B 6 metabolism highlighted for each comparison, we reported in Table 3 the vitamin B 6 intermediates according to KEGG metabolic pathway for Bos taurus. Overall, the vitamin B 6 intermediates revealed a significant up-accumulation of 4-pyridoxic acid (the degradation product of pyridoxal phosphate) for each comparison under investigation. Interestingly, the highest accumulation of 4-pyridoxic acid was found for the comparison HMC versus cluster 2 (log2FC = 0.18), followed by the comparisons against clusters 3 and 5 (Table 3). Taken together, our findings seem to Rocchetti et al.: CORN SILAGE-BASED STRATEGIES ON MILK METABOLOME Figure 3. Venn diagram considering the common and exclusive variable importance in projection discriminant compounds of the comparisons high-moisture ear corn (HMC) versus clusters 2, 3, and 5. support a direct correlation between the feeding strategies based on HMC and the accumulation of pyrimidine and pyridoxal degradation products in milk.

DISCUSSION
In this work, an untargeted metabolomics approach was used to discriminate bulk milk samples associated with different nutritional strategies based on the large utilization of corn silage. The different nutritional strategies (Supplemental Table S1) were previously grouped and classified according to an unsupervised hierarchical cluster analysis (Gallo et al., 2022) in which the different feed inclusion rate, expressed as a percentage on a DM basis, was used as the independent variable to cluster diets (categorical dependent or target variable). The utilization of this multivariate statistical approach could be particularly useful to identify, without a priori criterion and following a data-driven approach, differ-ent nutritional approaches to be further applied to farm data. Overall, the cluster analysis discriminated the farm feeding strategies in 6 main nutritional groups, namely HMC and legume silage (cluster 1), compound feed (cluster 2), corn and soy meals (cluster 3), HMC and soy meal (cluster 4), cornmeal and protein compound feeds (cluster 5), and HMC and protein compound feed (cluster 6). Interestingly, combining the information from unsupervised hierarchical clustering and supervised OPLS-DA, a successful separation of milk samples, according to the feeding strategies, could be observed. Therefore, untargeted metabolomics was then exploited to provide new insights into the most affected metabolomic pathways and potential discriminant compounds of physiological and nutritional processes.
Overall, the first result that deserves to be highlighted is the hyperspace separation provided by the OPLS-DA score plot and built considering the me-   Each metabolite is reported with its VIP (variable importance in projection) score (cutoff >0.8) and cross-validation SE, together with its log2 fold-change value for each comparison against high-moisture ear corn (HMC)-labeled milk samples. 2 MSMS = tandem mass spectrometry; MS1 isotopic = untargeted mass spectrometry (annotation based on accurate mass and isotopic profile).
3 OPLS-DA = orthogonal projections to latent structures discriminant analysis. Each metabolite is reported together with its log2 fold-change value for each comparison against high-moisture ear corn (HMC)-labeled milk samples.
tion observed in the OPLS-DA score plot (Figure 1) revealed that milk samples belonging to clusters 1, 4, and 6 were characterized by a similar chemical profile. Therefore, considering that this multivariate analysis tends to maximize the explained covariance between the latent variables to provide a clear separation of the observation groups (Worley and Powers, 2013), it is possible to hypothesize a common discriminant factor driving the hyperspace separation of clusters 1, 4, and 6. Interestingly, our previous observational study (Gallo et al., 2022) allowed us to identify this driving factor in the utilization of HMC that exclusively characterized (from 10% up to 18.3% of DM) the diet ingredient formulation of clusters 1, 4, and 6 (Supplemental Table  S1). The potential correlation between the use of HMC in the rations and the milk metabolomic profile represents a novel contribution to the field of dairy science, scarcely explored in scientific literature. The next step of this work was to evaluate the effect of HMC-based feeding strategies on the chemical profile of milk, extrapolating potential marker compounds and discriminant metabolomic pathways. As reported in Figure 2, the enrichment analysis carried out on the VIP marker compounds of the OPLS-DA model allowed us to observe a significant (P-value < 0.05) abundance of metabolites involved in the metabolic pathways of nucleic acid derivatives (purines and pyrimidines), AA (both essential and not essential), peptides, and pyridines (pyridoxal derivatives). Most of these significant classes of compounds discriminated the possible pairwise comparisons, namely "HMC versus cluster 2," "HMC versus cluster 3," and "HMC versus cluster 5" (Table 1). In our previous work (Gallo et al., 2022), we found that the inclusion level of HMC into diets increased the feed efficiency (P-value < 0.05), the kinetics of gas production (P-value < 0.05), crude protein, and starch digestibility of the diets (P-value < 0.05). Additionally, we observed production levels higher than 33 kg/cow per day in nutritional strategies including HMC (clusters 1, 4, and 6) despite differences in other silages or starch and protein sources. Usually, a greater ruminal and total-tract starch digestibility is measured when dairy cows are fed HMC compared with dry corn (Firkins et al., 2001;Ferraretto et al., 2013). This is likely related to the breakdown of the hydrophobic starch-protein matrix surrounding starch granules during ensiling (Gallo et al., 2014), which allows for greater microbial fermentation and enzymatic digestion of starch by ruminants. Although the starch digestibility of HMC can be altered by several factors (Gallo et al., 2022), the HMC has been previously described to possess a greater rumen rate of starch degradation than cornmeal, being the latter harvested in a more advanced stage of maturity (Masoero et al., 2011;Hoffman et al., 2012;Ferraretto et al., 2015), with potential effects on starch degradability and rumen fermentation syncronization. Accordingly, Eun et al. (2014) found that feeding HMC in high-forage diets increased NDF and CP digestibilities, microbial protein synthesis, together with feed and N utilization efficiencies, with a decrease in DMI relative to steam-flaked corn. Consequently, cows fed with HMC may have improved energetic efficiency compared with those fed with steam-flaked corn. However, there is no information in the scientific literature about how feeding HMC affects the milk metabolomic profile and final composition of milk and this work provides a preliminary indication of its effects on a large panel of dairy cow diets.
Taken together, the HRMS-based findings seem to suggest that the inclusion of HMC in the diet formulation had a great effect on the metabolism of nitrogen, with special reference to the accumulation of nucleic acids and their degradation products in milk ( Figures  2 and 4). Nitrogen (N) efficiency in dairy cows is low, and optimization of the diet, with a focus on dietary N in the form of protein, AA, and urea, has only led to minor improvements in the utilization of N in dairy cows (Fujihara and Shem, 2011). Also, the importance of other N-containing compounds such as microbial nucleic acids (i.e., DNA and RNA) in the nutritional physiology of ruminants has so far been scarcely investigated, although they correspond to more than 20% of the total microbial N synthesized in the rumen. Overall, the resynthesized microbial nucleic acids flow into the small intestine where they are digested before subsequent absorption takes place (Stentoft et al., 2015). Regarding specific biomarkers in milk, quantitative analysis of nucleic acid derivatives in dairy cattle research has almost solely focused on purine derivatives in both urine and milk, where the purine degradation products uric acid and allantoin have been proposed as indirect markers of rumen microbial synthesis (Schager et al., 2003). To go into detail, it was postulated that nucleic acids entering the duodenum are almost all of microbial origin. Purine nucleotides are degraded in the intestine and then absorbed. Adenine and guanine are catabolized and excreted via 4 possible routes, namely renal, mammary gland, recycling via saliva, and secretion into the intestines (Schager et al., 2003). The purine derivatives are excreted primarily as allantoin, but also as hypoxanthine, xanthine, and uric acid. Allantoin and uric acid have been related to the high xanthine oxidase activity in the blood and tissues, converting xanthine and hypoxanthine into uric acid before excretion. Interestingly, although some studies are available regarding milk allantoin extraction to estimate microbial protein flow to the duodenum (Stentoft et al., 2015), data concerning potential biomarkers related to other purines and mainly pyrimidines are at present very limited. In our experimental conditions, the purine metabolism was found to possess a nonsignificant effect on the metabolomic perturbations induced by HMC-based diets, while a higher effect of pyrimidine metabolism could be observed (Figure 4). Accordingly, lower up-accumulation values of both uric acid and allantoin were found in the HMC-related milk samples when considering each possible comparison (Table 2). Overall, the rumen inner environment has a great influence on milk quality. However, the correlation between the changes in rumen microorganisms and metabolites together with the udder health and milk metabolic profiles of dairy cows is still in short of knowledge. Our findings suggest a potential co-occurrence network in ruminal physiology and milk function considering the significant enrichment of not only pyrimidine degradation products (such as 3-ureidoisobutyric acid, 3-aminoisobutyric acid, and 3-ureidopropionic acid), but also pyrimidine biosynthetic intermediates (such as l-aspartate, N-carbamoyl-l-aspartate, orotic acid derivatives, and nucleotides; Table 2) required for their de novo synthesis (Löffler et al., 2015). A similar interconnection was recently reported by Wang et al. (2021), investigating milk and ruminal fluid of Holstein cows supplemented with Perilla frutescens leaf extracts. Regarding orotic acid, it represents a not very well studied milk component, and it is an indicator of the metabolic cattle disorder deficiency of uridine monophosphate synthase (Wehrmüller et al., 2008;Zaalberg et al., 2020). Orotic acid has been previously described as a by-product of protein biosynthesis in milk (Tiemeyer et al., 1984). However, the exact function of orotic acid in milk and its effect on calf health, the health of humans consuming milk or dairy products, manufacturing properties of milk, and its potential as an indicator trait, despite attracted the attention of different researchers in the past (Robinson, 1980;Tiemeyer et al., 1984;Saidi and Warthesen, 1989), appears still largely unknown.
Nucleotides, nucleosides, and nucleobases belong to the nonprotein-nitrogen fraction of milk, which are semi-essential dietary nutrients characterizing the milk of mammals and exerting many physiological functions. However, the degradation efficiency greatly depends on the rumen enzymatic activity, reflecting differences in microbial populations between animals and diets, and the rumen passage rate. In this regard, the total amount of DNA and RNA synthesized in the rumen depends largely upon the amount of bacterial growth; therefore, ad hoc studies combining metagenomic pro-filing of rumen microbiota and the metabolomic profile of different biofluids seem to be of great interest for the dairy industry, considering the need to accurately measure microbial protein supply and to better evaluate the shift of rumen microbiota toward the pyrimidine metabolism (as observed in this work by targeting the milk metabolome). According to literature (Kanehisa et al., 2014), the pyrimidine metabolites are degraded during pre-absorption, in the blood, and the hepatic tissue in much the same manner as the purine metabolites. However, data concerning the pyrimidine metabolic pathway in cattle is at present very limited. It is known that the pyrimidine end products (such as β-alanine and β-aminoisobutyric acid) can act as intermediate products in other nitrogen-related metabolisms (such as β-alanine metabolism together with valine, leucine, and isoleucine degradation). This could indicate that the degradation pathways of the pyrimidine metabolites differ from that of the purine metabolites. Generally, how much, and what types of purine and pyrimidine metabolites are absorbed from the intestinal lumen and to what extent they are metabolized across the hepatic tissue are unknown (Stentoft et al., 2015). To the best of our knowledge, the secretion and comprehensive profiling of pyrimidine metabolites into milk have not been studied in dairy cows and the final pyridine composition in milk represents a complex interplay involving total splanchnic fluxes from the rumen to mammary gland. In a previous work by Stentoft et al. (2015), the absorption and intermediary metabolism of purine and pyrimidines in lactating dairy cows was studied. Interestingly, the authors demonstrated that the pyrimidine metabolism differed from the purine metabolism in the splanchnic tissues, with a percentage of the total influx of approximately 50%, thus suggesting that the pyrimidine degradation enzymes in the hepatic tissue are not capable of removing all the pyrimidine intermediates and metabolites entering from the portal drained viscera and the peripheral tissues. This trend likely reflected the fact that much higher levels of pyrimidine nucleotides and degradation products entered the hepatic tissue and that the effectiveness of the hepatic tissue and the body's requirements or tolerance of the pyrimidine metabolites were different from that of the purine metabolites. The possible consequences related to the excess of pyrimidine metabolites are so far unknown; one possible hypothesis, strictly related to the findings observed in this work from milk metabolomics, could be their re-utilization in peripheral tissues or final excretion in biofluids, such as milk. Therefore, this work seems to support the hypothesis that pyrimidine intermediates and degradation products are potential marker compounds of the complex metabolism of mi-crobial nitrogen in dairy cows, and for the first time, an indirect correlation between HMC-based diet (described in the scientific literature to increase the rumen microbial synthesis) and pyrimidine derivatives in milk was outlined.
Finally, the pathway analysis ( Figure 4 and Table  3) also highlighted a significant effect of vitamin B 6 metabolism when considering the comparison "HMCbased" versus "the remaining clusters" (i.e., 2, 3, and 5). As a general consideration, our findings revealed a significant (P-value <0.05) up-accumulation of 4-pyridoxic acid (i.e., the degradation product of pyridoxal phosphate) in HMC-related milk samples. This result is not surprising considering that the HMC-based group also included cluster 1, which was characterized by a high inclusion of alfalfa silage in the rations (i.e., 14.9% of DM). Accordingly, a previous study by Castagnino et al. (2016) demonstrated that diets based on alfalfa silage resulted in greater degradation of vitamin B 6 in the rumen of dairy cows. In particular, the authors stated that changing the nutrient composition of the diets likely affected the microbial population in the rumen and their B vitamin metabolism. In another work, the increase in dietary forage concentration was found to increase the ruminal degradation of vitamin B 6 (Seck et al., 2017). Last, Beaudet et al. (2020) evaluated the addition of extruded linseed in the diet of dairy cows, showing that increasing linseed concentration in isonitrogenous and isoenergetic diets increased vitamin B 6 and folate supply to dairy cows, both with hay and corn silage-based diets. Consequently, further studies are needed to verify how dietary strategies can affect the metabolism of vitamin B 6 . Therefore, the significant effect of vitamin B 6 metabolism (i.e., pyridoxal related compounds) outlined by the pathway analysis in milk could be explained by the inclusion of alfalfa silage in the rations.

CONCLUSIONS
In this work, an untargeted metabolomics approach based on HRMS was used to provide new insights into the potential correlation existing between different nutritional strategies of dairy cows (previously characterized according to the different diet formulations) and the chemical profile of milk. Overall, a combination of different multivariate statistical approaches allowed us to identify a potential effect of HMC on pyrimidine metabolism. The accumulation of both pyrimidine intermediates and their degradation products could be considered in the future as potential markers of HMCbased diets, known to be characterized for greater microbial fermentation and enzymatic digestion of starch by ruminants. Also, our preliminary metabolomic findings suggest a complex metabolism of microbial nitrogen based on splanchnic fluxes from the rumen to mammary gland in dairy cows. Therefore, further works, based on the comprehensive evaluation of other biofluids (such as rumen, plasma, and urine) appear to be of great interest to better elucidate the biochemical and physiological role of pyrimidine metabolites.