Assessment of milk metabolites as biomarkers for predicting feed efficiency in dairy sheep

Estimating feed efficiency (FE) in dairy sheep is challenging due to the high cost of systems that measure individual feed intake. Identifying proxies that can serve as effective predictors of FE could make it possible to introduce FE into breeding programs. Here, 39 Assaf ewes in first lactation were evaluated regarding their FE by 2 metrics, residual feed intake (RFI) and feed conversion ratio (FCR). The ewes were classified into high, medium and low groups for each metric. Milk samples of the 39 ewes were subjected to untargeted metabolomics analysis. The complete milk metabolomic signature was used to discriminate the FE groups using partial least squares discriminant analysis. A total of 41 and 26 features were selected as the most relevant features for the discrimination of RFI and FCR groups, respectively. The predictive ability when utilizing the complete milk metabolomic signature and the reduced data sets were investigated using 4 machine-learning algorithms and a multivariate regression method. The Orthogonal Partial Least Square algorithm outper-formed other ML algorithms for the FCR prediction in the scenarios using the complete milk metabolite signature (r2 = 0.62 ± 0.06) and the 26 selected features (0.62 ± 0.15). Regarding RFI predictions, the scenarios using the 41 selected features outperformed the scenario with the complete milk metabolite signature, where the Multilayer feedforward artificial neural network (r 2 = 0.18 ± 0.14) and extreme gradient boosting (r 2 = 0.17 ± 0.15) outperformed other algorithms. The functionality of the selected metabolites implied that the metabolism of glucose, galactose, fructose, sphingolipids, amino acids, insulin, and thyroid hormones was at play. Compared with the use of traditional methods, practical applications of these biomarkers might simplify and reduce costs in selecting feed-efficient ewes.


INTRODUCTION
Feed efficiency (FE) is an essential phenotype in animal production as it represents the capacity of an animal to convert feed into a product.Identifying more efficient animals might also potentially reduce production costs and environmental impact simultaneously with maximization of profitability (Connor, 2015;Madilindi et al., 2022).The FE in sheep is more difficult to predict and measure than in cattle due to a wider range of farming technologies, feeding systems, geographical areas and populations (Cannas et al., 2019).Two commonly studied measures of feed efficiency (FE) are the feed conversion ratio (FCR) and the residual feed intake (RFI).The first one, FCR, is calculated by dividing feed intake by weight gain, and the RFI is the difference between actual and predicted feed intake required to achieve a certain rate of gain (Koch et al., 1963).Similarly, in dairy animals, these measures are estimated based mainly on milk production (Connor et al., 2013;Barrio et al., 2023;Toral et al., 2023).
The difficulty and the costs of calculating FE on a large scale to be used in breeding programs is a challenge for dairy sheep; however, this information is crucial for the industry.However, gathering data for FE estimation in dairy breeds poses challenges, requiring the exploration of alternative approaches.A potential solution is to establish a reference population and utilize biomarkers to effectively train predictive models for estimating FE in additional animals.
Milk metabolites have the potential to serve as valuable biomarkers, as they represent the end products of a wide range of biological processes and physiological mechanisms associated with FE.Recently, metabolomic signatures in plasma and milk were described as effectively discriminating Assaf lactating ewes that were divergent for FE (Toral et al., 2023).While the Assessment of milk metabolites as biomarkers for predicting feed efficiency in dairy sheep H. Marina, 1 J. J. Arranz, 1 * A. Suárez-Vega, 1 R. Pelayo, 1 B. Gutiérrez-Gil, 1 P. G. Toral, 2 G. Hervás, 2 P. Frutos, 2 and P. A. S. Fonseca 1 potential applications of these findings are promising, it is equally important to explore a methodology that enables the prediction of actual FE.Developing such an approach would be an intriguing alternative, as it would provide a more comprehensive and inclusive assessment of FE, allowing for more accurate predictions that would be applicable across diverse populations and enabling individual animals to be ranked by those values.Therefore, identifying a set of metabolites that can be easily measured at the locus using quick and costeffective methods, such as mid-infrared spectroscopy, and demonstrating an efficient predictive performance for FE metrics holds great potential for the dairy sheep industry.
Recently, the application of machine learning (ML) models, such as partial least squares discriminant analysis (PLS-DA), has proven effective in leveraging rumen and plasma metabolomic data to differentiate animals with high and low RFI values (Touitou et al., 2022).Supervised and semi supervised methods using genomic and metabolomic data have shown promising results in predicting RFI in pigs and cattle, with promising levels of predictive accuracy (Yao et al., 2016;Martin et al., 2021;Piles et al., 2021), which indicates the potential for improving FE estimation in other livestock species.
In light of the above, this study aimed to achieve the following objectives: 1) Assess the effectiveness of utilizing complete metabolomic signatures obtained from milk samples of Assaf ewes in predicting both RFI and FCR; 2) identify a subset of milk features that demonstrate comparable or superior predictive performance for RFI and FCR compared with the complete metabolomic signature; and 3) annotate the selected features to enhance the understanding of the underlying biological mechanisms involved in feed efficiency (FE) predictions, specifically in dairy sheep.

Ethics Statements
All experimental procedures were approved by the Research Ethics Committees of the Instituto de Ganadería de Montaña, the Spanish National Research Council (CSIC), and the Junta de Castilla y León (Spain), following procedures described in Spanish and European Union legislation (Royal Decree 53/2013 and Council Directive 2010/63/EU).

Sampling and feed efficiency estimations
The studied animal population involved 40 Assaf ewes in their first lactation.These ewes underwent estrus synchronization to minimize lactation stage variations.
The ewes were individually housed in pens and milked twice daily using a 1x10 stall milking parlor (DeLaval).The ewes had ad libitum access to a total mixed ration (TMR) composed of alfalfa hay and concentrate in a 50:50 forage: concentrate ratio.Detailed information regarding diet formulation and animal management protocols can be found in Barrio et al. (2023).
After a 3-week adaptation period, the dry matter intake (DMI) and milk yield and composition of the ewes were measured individually and daily for 3 weeks.At this stage, all the ewes were in the second month of lactation (60.59 ± 7.93 d in milk).The daily feed intake was estimated by weighing the dry matter offered and refused.In addition, the morning and evening milk production of each ewe was also measured to calculate the daily milk yield.The protein, fat, and lactose content of the milk were determined as described by Barrio et al. (2023).Furthermore, body weight (BW) changes were calculated for each ewe by recording the weights on 2 consecutive days at the beginning and end of the assay.One ewe was excluded due to acute mastitis during the trial; consequently, the data of 39 Assaf ewes were available.
The FCR and RFI of each ewe were estimated for the evaluated period.The FCR was estimated as the ratio between DMI and energy-corrected milk yield (ECM, L/d), where ECM was calculated using INRA ( 2018 The RFI was estimated using the GLM procedure of the SAS software package (version 9.4; SAS Institute Inc.) and expressed as residuals of the following regression model: where DMI is the mean dry matter intake during the experiment (kg/d); μ is the intercept; ECM is the energy-corrected milk yield (kg/d); MBW is the mean metabolic body weight during the experiment (BW 0.75 ; kg); BWC is the mean metabolic body weight change (BWC 0.75 ; kg); and RFI is the residual.

Metabolomics
Sample processing.A 50 mL milk sample was collected from each of the 39 ewes for untargeted metabolomics on the last day of the trial.Each sample was aliquoted into 200 μL portions and stored at −80°C.Each aliquot was mixed with 1 mL of methanol-chloroform (3:1) in an Eppendorf tube.The mixture was vortexed at 3000 rpm for 1 min and then sonicated in an ice-water bath for 2 min.Following this, centrifugal precipitation was performed at 21,000 g and 5°C for 10 min.The resulting transparent supernatant was collected for further processing.
Reversed-phase liquid chromatography-mass spectrometry (LC-MS).Reversed-phase LC-MS was performed to detect features ranging from polar to less polar compounds and lipids.For each sample, 200 μL of the supernatant was transferred to an LC injection vial.Subsequently, 5 μL aliquots of the supernatant from each sample were injected into a Waters BEH reversephase LC column.High-resolution LC-MS analysis was conducted using a Thermo Ultimate 3000 UHPLC system coupled to a Thermo LTQ-Orbitrap Velos Pro mass spectrometer equipped with an electrospray ionization source.The column flow rate was set to 400 μL/min, and the column temperature was maintained at 55°C.Each sample underwent 2 rounds of sample solution injections, one for positive ion detection and another for negative ion detection.For compound detection and relative quantitation, the MS instrument was operated in the Fourier transform MS detection mode at a mass resolution of 60,000 FWHM @ m/z (mass-to-charge ratio) 400 and within a mass range of m/z 50 to 1800.To ensure the quality of LC-MS detection during data acquisition, a quality control (QC) sample was prepared by pooling samples from 10 randomly selected test samples.The QC sample was injected in parallel with the test samples in completely random order during data acquisition.

Hydrophilic interaction LC-MS (HILIC-MS).
HILIC-MS was performed to detect highly polar features.For each sample, 600 μL of the supernatant was mixed with 150 μL of water and 200 μL of chloroform.The mixture was vortex-mixed at 3,000 rpm for 1 min, followed by centrifugation for 5 min to achieve clarification.Next, 500 μL of the upper aqueous phase was collected and dried in an LC injection microvial using a nitrogen gas flow.The residue was reconstituted in 600 μL of 70% acetonitrile for the subsequent LC-MS runs.Subsequently, 5 μL aliquots from each sample solution were injected into a Waters HILIC column.The highresolution LC-MS analysis was conducted using the same equipment and a similar procedure as described above.In this case, the column flow rate was set to 300 μL/min, and the column temperature was maintained at 40°C.The mass detection range for this analysis was m/z 50 to 1000.
Data processing.Four high-resolution LC-MS data sets were obtained, and each data set was processed using the R package XCMS (Smith et al., 2006).The processing steps included peak detection, peak grouping, peak alignments, and retention time shift correction.The resulting data set included pairs of MS m/z and LC retention time (RT, min) for the detected compound features across all the samples within each data set.This procedure revealed 3760 feature signals.Those for which less than 75% of the data were available from the ewes were filtered out, resulting in a final data set composed of 3749 features.As untargeted metabolomics data are compositional in nature (Sisk-Hackworth and Kelley, 2020), the feature signals were normalized using a centered log-ratio (CLR) transformation.The processed data were implemented in subsequent statistical analyses.

Analysis of the discrimination between high and low FE groups and selection of features with better discriminating profiles
The ewes were classified into high, medium, and low FE groups for RFI (H-RFI, M-RFI, L-RFI) and FCR (H-FCR, M-FCR, L-FCR) based on the distribution of each metric.For the animal classification in the High, Medium and Low FE groups, first, the distribution of the individual values for RFI and FCR was analyzed.Subsequently, the rank of the animals was compared for both metrics and the 10 animals from both extremes of the distribution (10 High FE and 10 Low FE) that did not show big discrepancies between RFI and FCR metrics were selected to compose the High and Low groups of each metric.The rest of the animals were grouped in the respective Medium FE groups.The mean of each group was compared using an ANOVA test using R v.4.2.0, where the significant differences among all the groups were defined based on a p-value < 0.05.After this step, Tukey's post-hoc analysis was performed to identify the significant differences between each pair of groups.Additionally, the mean DMI was compared between each group for RFI and FCR using the same approach.All p values were adjusted to a false discovery rate (FDR) threshold of 0.05.The R package mixOmics (Rohart et al., 2017) was used to perform a partial least squares discriminant analysis (PLS-DA) to evaluate the potential of selected features to discriminate the high, medium, and low FE groups for both RFI and FCR.The area under the curve (AUC) was calculated to estimate the discrimination potential of PLS-DA among the 3 FE groups.The first 2 principal components were plotted with centroids to define the groups identified by the PLS-DA using the plotIndiv function in the mixOmics package.Additionally, the performance of PLS-DA for the RFI and FCR groups was evaluated based on the mean error rate and Q 2 for the first 5 principal components through crossvalidation using 10-folds.
Once the ewes were assigned to each group, the variable importance in the projection (VIP) was estimated to compute the influence on the features of every predictor in the PLS-DA model.Features with a VIP > 2 were selected individually from the outputs of the discriminant analysis for RFI and FCR.

Annotation of selected features and metabolic pathway enrichment analysis
The features obtained through the discriminant analysis for RFI and FCR groups were matched against the available information in MetaboAnalyst 5.0 (Pang et al., 2022) by comparing their m/z values to 3 different databases: BioCyc (Karp et al., 2018), KEGG (Kanehisa and Goto, 2000), and MNF (a manually curated database, originally from the mummichog package (https: / / shuzhao -li .github.io/mummichog .org/index .html),that integrates KEGG, BiGG, and Edinburgh model data).Additionally, the selected features were matched against the information available in MassBank (Horai et al., 2010).During the identification process, a maximum error of 5 ppm was permitted for all database searches, and a maximum difference of 0.1 between the m/z of the selected features and potential matches was allowed.Only potential matches with the same ion form (positive-ion or negative-ion modes) were considered for the selected features.All the potential matches were investigated regarding their presence in milk, and those previously identified as milk-related metabolites were retained for functional analysis.Functional enrichment analysis was performed, including the potential matches for the features exclusively identified for RFI and FCR and those shared between both metrics in the IMPaLA software (Kamburov et al., 2011).The enriched pathways were defined using an FDR threshold of 0.05.

Prediction of RFI and FCR using machine learning algorithms
Predictive models were built using the CLR-transformed values for each of all 3749 detected features (RFI_all and FCR_all) and for the subsets of features selected in the PLS-DA analysis for RFI and FCR with VIP > 2 (RFI_VIP and FCR_VIP).The ML models were built based on a multilayer feedforward artificial neural network (deep) and random forest (RF) using the R package h2o (Candel, Arno, Viraj Parmar, Erin LeDell, 2016), extreme gradient boosting (xgboost) using the R package xgboost (Chen and Guestrin, 2016), and support vector machine (SVM) using the caret R package (Kuhn, 2008).For the deep and RF algorithms, a random-discrete search composed of 5-fold cross-validation was performed during the hyperparameter tuning stage.The hyperparameter tuning stage for SVM and XGBoost was performed using an extensive search.The tuning stage for xgboost used 1,000 as maximum boosting interactions and 10 interactions as early stopping (stop if there is no improvement for 10 consecutive trees).SVM hyperparameter tuning was performed using the caret package with the "repeatedcv" method with 10-folds and 5 repeats.The grids explored in each ML algorithm are presented in Supplemental Table 1 (https: / / zenodo .org/record/ 8154651).
The effect of the reduced sample size in the current experiment was evaluated by different compositions of training and test sets.In total, 100 models were constructed for each ML algorithm evaluated here.Each model was built randomly by assigning 39 ewes to train (2/3, 26 animals) and testing (1/3, 13 animals) samples.It is important to note that in all 100 iterations of the random assignments, the training and testing samples were not perfectly equal.The model evaluation was performed as described by Martin et al. (2021).The precision of each model generated was determined using the squared Pearson correlation (r 2 ).The accuracy was evaluated by the root mean squared error (RMSE) and the concordance correlation coefficient (CCC) between the predicted and actual values of RFI and FCR.The best model among the 100 models generated in each scenario was defined based on the maximum ratio between r 2 and RMSE values.
To compare the performance of ML models with traditional regression models, the RFI and FCR values were predicted using the complete milk metabolomic signature and the reduced feature data sets by applying an orthogonal partial least squares (OPLS) model.Similar to the workflow for ML models, 100 iterations of random sample assignments in the training and test sets were performed.The optimal number of orthogonal components (from 1 to 9) was determined within each iteration based on the maximum cumulative Q2Y metric.The models were fitted using the ropls package (Thevenot, 2016) in R. The predictive performance of the OPL models was also evaluated based on RMSE and r 2 .The best model for the RFI and FCR predictions using the complete milk metabolomic signature and the reduced feature data sets was also defined based on the maximum ratio between r 2 and RMSE values.
A summary of the workflow applied for the prediction of RFI and FCR values is shown in Figure 1.

Distribution of residual feed intake and feed conversion ratio among groups
The RFI and FCR values for the 39 ewes are available in Supplemental Table 2 (https: / / zenodo .org/record/ 8154651).For RFI, 10, 19, and 10 ewes were classified in the H-RFI, M-RFI and L-RFI groups, respectively.The classification based on the FCR values resulted in 10, 18, and 11 ewes in the H-FCR, M-FCR, and L-FCR groups, respectively.Significant differences (p value < 0.05) were observed among the groups within RFI and FCR (Table 1), indicating that the division of animals into distinct groups based on the values obtained for each FE metric was effective in creating meaningful distinctions in terms of RFI and FCR.A low and nonsignificant Pearson correlation was observed between the RFI and FCR values (r = 0.28, p value = 0.08).Additionally, no significant differences (FDR > 0.05) were observed in the pairwise group comparisons for RFI and FCR regarding DMI.

Metabolomics
The features composing the complete milk metabolomic signature identified in the 4 high-resolution LC-MS data sets evaluated here and the respective m/z and LC retention time (RT, min) are available in Supplemental Table 3 (https: / / zenodo .org/record/ 8154651).In total, 3469 features were retained after the processing and filtering steps for the statistical and functional analyses.

Discriminant analysis and feature selection
The PLS-DA results indicated a high potential of the complete milk metabolomic signature to discriminate the high, medium, and low FE groups for both RFI and FCR metrics (Figure 2).The AUCs obtained for the discrimination among the H-RFI, L-RFI, and M- RFI groups and the other groups were 0.969, 0.951, and 0.992, respectively.In addition, for FCR, the AUCs obtained for the discrimination among H-FCR, L-FCR, and M-FCR groups and the other groups were 0.976, 0.987, and 0.981, respectively.The mean error rate and Q 2 metrics for the RFI groups were 0.621 and 0.379, respectively, while those for the FCR groups were 0.426 and 0.574, respectively.Therefore, all 3 groups were well discriminated for RFI and FCR.For the reduced set of features, 41 and 26 features with VIP > 2 were selected for RFI and FCR, respectively (Supplemental Table 4, https: / / zenodo .org/record/ 8154651).

Predictive accuracy of residual feed intake and feed conversion ratio using milk features
In general, the predictions for FCR outperformed the predictions for RFI.Higher means of r 2 and CCC were observed for FCR predictions in all scenarios (Tables 2 and 3).On the other hand, a larger mean RMSE was obtained for FCR in all scenarios.However, RMSE should not be directly compared between the metrics, as RFI and FCR have different distributions and scales.In the scenario where all the features were used for prediction, the highest r 2 for FCR was obtained using OPLS (0.62 ± 0.06), while that for RFI was obtained using RF (0.08 ± 0.10).Regarding RMSE, the smallest mean for FCR was obtained using SVM (0.14 ± 0.03) and that for RFI using xgboost (0.14 ± 0.02).Regarding the CCC values, in all the scenarios where all the features were used for RFI prediction, means close to zero were obtained, while in the RFI predictions, using the 41 features selected based on VIP > 2 resulted in substantially higher CCC means (Table 2).The exception was the scenario where OPLS was used to predict RFI with the subset of 41 features, which also resulted in a mean close to zero (−0.04).The 41 features selected for RFI (VIP > 2) provided better predictions (higher r 2 ) for RFI in all the models, excluding OPLS.The highest r 2 mean for RFI using the 41 features was obtained using the deep model (0.18 ± 0.14).Similarly, for FCR, only SVM was not outperformed by the predictions performed using the 26 features selected in the PLS-DA analysis for FCR groups.The mean r 2 and CCC for all features using SVM were 0.58 ± 0.15 and 0.67 ± 0.10, while using the 26 features, the r 2 and CCC means were 0.51 ± 0.17 and 0.61 ± 0.11, respectively.The predictions performed using OPLS in the reduced data set performed exactly (for r 2 , CCC, and RMSE) as the complete milk metabolomic signature.The highest r 2 means obtained using ML models for FCR predictions based on the 26 features was with the xgboost algorithm (0.55 ± 0.15).On the other hand, higher CCC means were obtained for RF and OPLS, 0.68 ± 0.11 and 0.68 ± 0.10, respectively.Therefore, in some scenarios, it was possible to note only a slightly worse performance for the FCR predictions using the selected subset of features when r 2 , CCC, and RMSE were evaluated.
An important variation was observed across the 100 iterations regarding the predictive performance (Tables 2 and 3).These results reinforce the potential effect of different compositions of training and test sets in reduced data sets.The best model across the 100 iterations was defined for each scenario (RFI_all, FCR_all, RFI_VIP, FCR_VIP) based on the maximum ratio between r 2 and RMSE.
The best overall model for RFI prediction using the complete milk metabolomic signature was obtained through the RF algorithm (r 2 = 0.68, RMSE = 0.20).For the scenario where FCR values were predicted using the complete milk metabolomic signature, the best overall model was obtained using SVM (r 2 = 0.79, RMSE = 0.07).Regarding the scenarios of RFI and FCR predictions using the respective reduced feature data sets, the models were obtained using xgboost (r 2 = 0.53, RMSE = 0.07) and OPLS (r 2 = 0.83, RMSE = 0.09), respectively.The full description of all the models tested for each ML algorithm (train and test composition, hyperparameters, r 2 , and RMSE) is available in Supplemental Table 5 (https: / / zenodo .org/record/ 8154651).

Annotation of selected features and pathway enrichment analysis
The potential matches for the 41 and 26 features with VIP > 2 selected for RFI and FCR during the discriminant analysis are available in Supplemental Table 6 (https: / / zenodo .org/record/ 8154651).It is important to highlight that some features presented more than one candidate-matched metabolite (Figure 3).In total, 29 and 15 features showed a potential match during the annotation process out of the 41 and 26 features selected for RFI and FCR, respectively.The 29 RFI features were assigned to 123 potential matches, while the 15 FCR features were assigned to 76 potential me- tabolites (Figure 3).Twenty-three potential matches were shared between RFI and FCR.Consequently, 100 and 53 potential matches were exclusively identified for the features selected for RFI and FCR, respectively.
The enriched pathways and the potential corresponding match of the metabolites identified exclusively for RFI, FCR, and both metrics are available in Supplemental Table 7 (https: / / zenodo .org/record/ 8154651).Figure 4 shows the distribution of the enriched pathways regarding the FDR values, the richness factor (number of annotated metabolites for a specific pathway divided by the total number of metabolites assigned to this pathway in the database), and the overlapping of enriched pathways identified for the potential matches of RFI, FCR and both metrics.
Regarding the enriched pathways identified for each matched metabolite data set, 15 and 14 enriched pathways were identified exclusively for the RFI and FCR data sets, respectively.On the other hand, 71 enriched pathways were shared between RFI-and FCR-matched The ROC curves show the area under the curve (AUC) obtained for discriminating between one group and the others using the first 2 principal components.The X and Y axes in the left-hand side plots represent the percentage of the variance explained by the first and second components identified in the PLS-DA analysis.For right-hand side plots, the x-axis shows the 100-specificity (100%) while the y-axis shows the sensitivity (%) obtained in the discriminant analysis for the 3 feed efficiency groups.
metabolites in different contexts: 5 enriched pathways shared among the 3 data sets (RFI, FCR, and shared); 26 exclusively enriched for the matched metabolites shared between RFI and FCR; 10 shared between the metabolites exclusively matched for RFI or FCR; 23 shared between the metabolites exclusively matched for FCR and those matched metabolites shared between RFI and FCR; and 7 shared between the metabolites exclusively matched for RFI and those matched metabolites shared between RFI and FCR.
The shared enriched pathways for the candidate metabolites identified for RFI and FCR (or for the commonly matched metabolites between both) are mainly related to glycolysis, gluconeogenesis, galactose metabolism, and fructose metabolism.The enriched pathways exclusively identified for the RFI-matched metabolites were mainly related to lactose and amino acid metabolism (emphasis on lysine).For FCR, the enriched pathways were associated with fructose metabolism, insulin action, thyroid hormone synthesis, sphingolipid de novo biosynthesis, and energy metabolism.

DISCUSSION
This study aimed to contribute valuable insights into the potential application of milk metabolomics for accurate FE estimation and to advance the understanding of the biological factors influencing FE in dairy sheep.Here, the predictive ability of the complete milk metabolomic signature and the reduced data sets selected for RFI and FCR were investigated using 4 ML algorithms and a multivariate regression method.

Marina et al.: Milk Metabolites and Feed Efficiency in Sheep
Table 2. Summary statistics (mean ± standard deviation (sd), minimum (min), and maximum (max)) for the squared Pearson correlation (r 2 ), concordance correlation coefficient (CCC), and root mean squared error (RMSE) values obtained for the predictions of residual feed intake (RFI) using the complete milk metabolomic signature (all) and the reduced feature data sets selected based on Variable Importance in the Projection (VIP) obtained in the PLS-DA analysis The metabolomic signature in the milk of divergent animals for FE has the potential to provide informative biomarkers.Indeed, the PLS-DA analysis resulted in an AUC greater than 0.96 for the discrimination of the RFI (H-RFI, L-RFI, and M-RFI) and FCR (H-FCR, L-FCR, and M-FCR) groups.Similar discriminant potentials between FE groups were obtained using serum and milk features for divergent FE Assaf ewes in a previous study (Toral et al., 2023) and superior to those obtained using serum metabolites in crossbreed meat sheep divergent for RFI (Goldansaz et al., 2020).Additionally, the PLS-DA analysis allowed the identification of features with higher importance for the discrimination of FE groups for both RFI and FCR metrics.In the present study, the highest mean r 2 was obtained using the OPLS model for both feature data sets.The results of the predictions of FCR showed similar performance of the complete feature data set and the reduced data set (composed of 26 features).Consequently, the results obtained here suggest the reduced data sets can be used effectively to predict RFI and FCR.In all the scenarios, the predictive performance for FCR outperformed the predictions of RFI.There are several hypotheses that could explain these results.For example, the RFI is the error term of a regression model used to predict DMI (Wu et al., 2021).Suboptimal adjustment of this model could result in less accurate RFI metrics and, consequently, impact the RFI predictions using the milk features.Here, 39 ewes were used to estimate RFI, which could be considered a moderate sample size.Therefore, larger sample sizes for RFI predictions would increase the accuracy of RFI estimations and help improve predictive performance using milk features.Another hypothesis relies on the fact that FCR is the ratio between DMI and ECM.This ratio could represent a more direct metric related to milk production when compared with the RFI due to the high genetic correlation between FCR and milk yield in dairy cattle (Oldenbroek, 1989).RFI could be interpreted as a more complex metric representing different energy allocations in the animal, displaying a low genetic correlation with milk production traits (Veerkamp et al., 1995).Nevertheless, it is important to mention that using FCR in breeding programs might be considered inappropriate because it is a ratio, which might result in unpredictable selection if the ratio components are disproportionate (Connor et al., 2012).
The prediction of individual FE metrics instead of the classification of animals in divergent groups is difficult, as reflected in the limited studies conducting the same methodology.In dairy cows, milk components and milk yield for the prediction of RFI results in r 2 values of 0.01 to 0.06 (Martin et al., 2021).In the current study, the highest mean of r 2 outperformed the abovementioned results using all milk features and the reduced feature set.Additionally, our results indicated that the reduced set of features (41 features with VIP > 2) selected for RFI outperformed the complete set of features for the predictions in all ML algorithms.
Regarding the effect of the reduced sample size on ML prediction accuracy, the current experiment evaluated different training and test set compositions (different groups of animals assigned to the training and test sets in each iteration).In omics studies, where reduced sample sizes are usually observed, providing a unique combination of hyperparameters and features to be applied in new samples is not a simple task.The results obtained here suggested large variability across the models generated for each scenario.This variability could be explained by the reduced sample size or other concepts such as data leakage (use of information in the training set, which is not expected to be available in the testing step) (Kaufman et al., 2011).In the current study, we observed similar levels of variability across the 100 iterations for the data sets composed of all milk features and those previously selected for RFI and FCR in the PLS-DA analysis, which could reduce the probability of data leakage.An in-depth analysis of the model structure and composition (features and samples) might be useful to clarify the applicability of such a model.The exact cause of this variability and outperformance of some scenarios can be identified by applying methodologies such as explainable artificial intelligence (XAI).Therefore, applying such methodologies would help improve the understanding of the relationship between the samples, which consist of the training and test sets with the features that most contribute to the overall model performance.Furthermore, leveraging XAI techniques makes it possible to improve the accuracy of predictions and facilitate the selection of more efficient biomarkers, thereby optimizing livestock breeding strategies and promoting advancements in the field.
The annotation process for features is not a simple task either and still requires detailed work for each potential match for the list of candidate features (Alseekh et al., 2021).Successful annotation of potential matches was observed for 29 and 15 features out of the 41 and 26 features selected for RFI and FCR, respectively.Additional studies are needed to validate the potential matches for each feature selected here for the predictions of RFI and FCR and to identify potential matches for those not annotated in the current study.
The functional profile of the potential matches for the selected features for RFI and FCR was evaluated regarding the enrichment of metabolic pathways.These were classified into 3 groups: enriched pathways shared between RFI and FCR, enriched pathways exclusive for RFI, and enriched pathways exclusive for FCR.
Sugar metabolism is strongly represented in the enriched pathways shared between RFI and FCR candidate features.Interestingly, in sheep, the feeding level and nutrient availability were associated with glucose uptake and utilization through alterations in the expression of genes involved in glucose metabolism, such as GLUT1, GLUT3, SGLT1, β-(1,4)GAT1 and LALBA (Tsiplakou et al., 2015).Glucose forms part of the main milk disaccharide, lactose (Glänze, 2022).Here, several features were associated with glucose metabolism, such as D-glucose 6-phosphate, D-glucose 1-phosphate, β-Dglucose 6-phosphate, D-fructose 6-phosphate, L-malate, L-lactate, L-valine, L-isoleucine, and L-lysine, among others.Pathways associated with galactose metabolism were also frequently observed as enriched for RFI and FCR candidate metabolites (Figure 4 and Supplemental Table 7, https: / / zenodo .org/record/ 8154651).Galactose has been suggested in cows as an interesting biomarker for metabolic status and energy balance (Pires et al., 2022), underscoring our results.Several candidate metabolites were identified as associated with galactose metabolism, such as D-galactose, UDP-α-D-galactose, D-mannose, lactose, galactinol, myo-inositol, N-acetyl-D-glucosamine, and D-tagatose 6-phosphate.Fructose metabolism was also strongly represented in the enriched pathways shared between RFI and FCR.Fructose metabolism was previously found to be enriched in the data sets of metabolites differentially abundant in the plasma of dairy cows during parturition and in the metagenomic data differentially abundant in the rumen of divergent cows for RFI (Luo et al., 2019;Xie et al., 2022).The metabolites β-D-fructose 6-phosphate, Dmannose 6-phosphate, and D-mannose 1-phosphate are among those associated with fructose metabolism and annotated for both RFI-and FCR-selected features.In addition, it is worth highlighting the presence of the enriched pathways "Sphingolipid Metabolism" and "Metabolism of proteins," which are associated with lipids and protein production.The sphingolipid, specifically the sphingolipid ceramide, is proposed as a potential antagonist of insulin-stimulated glucose during lactation (McFadden and Rico, 2019).Among the candidate metabolites associated with "Sphingolipid Metabolism," it is interesting to highlight galactosylceramide, ethanolamine phosphate, PC(20:3(8Z,11Z,14Z)/18:1(9Z)), sphingomyelin, dehydrosphinganine, and sphingosine.
The pathways exclusively enriched for RFI are mainly associated with amino acid metabolism ("Lysine degradation," "Lysine degradation I (saccharopine pathway)," and "Metabolism of amino acids and derivatives") and translocation of substrates across the membrane.Similar pathways were observed to be enriched in the study by Toral et al. (2023).Therefore, the levels of lysine or features associated with lysine metabolism could be attractive proxies for FE in dairy sheep.Among the RFI features associated with amino acid metabolism are L-malate, 3,4-dihydroxy-L-phenylalanine, L-2-aminoadipate, creatinine, L-isoleucine, L-glutamate 5-semialdehyde, L-valine, L-glutamate, 2-oxoadipate, N6-(L-1,3-dicarboxypropyl)-L-lysine, carnitine, betaine, N-acetyl-L-glutamic acid, L-lysine, and L-leucine.Among these features, carnitine plays a crucial role in energy metabolism via mitochondrial β-oxidation (De Vivo and Tein, 1990) and in health aspects in ruminants (Ringseis et al., 2018).Elevated carnitine concentrations were observed in low-RFI beef steers (higher-efficiency animals) (Clemmons et al., 2017).The increased carnitine concentration might be associated with more efficient nutrient utilization (Carlson et al., 2007;Murali et al., 2015), and supplementation with carnitine results in higher plasma glucose concentrations in cattle and sheep (Greenwood et al., 2001;Çetin et al., 2003).
For pathways exclusive to FCR, it is relevant to highlight the enriched pathways associated with insulin action, thyroid hormone synthesis, sphingolipid de novo biosynthesis, and energy metabolism.The features associated with sphingolipid de novo biosynthesis (NADP+, phytosphingosine, sphingosine, and 3-dehydrosphinganine) were also associated with the "sphingolipid metabolism" pathway, which was observed to be enriched for both RFI and FCR.In rams, RFI was significantly associated with the animal response to stress (cortisol levels) events expressed as insulin-induced hypoglycemia (Knott et al., 2010).Additionally, the divergent selection for insulin-like growth factor-I (IGF) suggests a minimal effect on RFI in beef cattle (Lancaster et al., 2008), and polymorphisms in the IGF2 gene were associated with average daily gain and FCR in beef cattle (Sherman et al., 2008).Interestingly, exposure to heat stress in cows decreased milk yield and resulted in higher insulin secretion in lactation periods (Itoh et al., 1998).Therefore, these results reinforce the impact of insulin action over FE and the putative role of the features identified here (D-erythrose 4-phosphate, Dfructose 6-phosphate, and ribose 5-phosphate) as associated with this pathway for the prediction of FCR.D-Erythrose 4-phosphate was previously identified as differentially abundant in the ruminal fluid between beef steers showing high and low average daily gain (Artegoitia et al., 2017).Consequently, this underscores the potential use of these features as a proxy for FCR in dairy ewes.The "thyroid hormone synthesis" pathway was also exclusively related to the FCR features NADP+, glucose 6-phosphate, and ribose 5-phosphate.Thyroid hormones have a galactopoietic action and help in the development of the mammary gland during lactation (Capuco et al., 2008).In the absence of thyroid hormones, growth and differentiation of the mammary epithelium are reduced (Vonderhaar and Greco, 1979).In addition, the thyroid hormone synthesis pathway was previously identified in beef cattle as a key component of the pleiotropic effect observed among feed efficiency, reproduction, and production-related traits (Fonseca et al., 2018).Therefore, features associated with the regulation of thyroid hormone synthesis might be interesting proxies for the prediction of FE metrics, such as FCR.
The selected features for RFI and FCR are enrolled with biological processes related to milk production and/or composition.However, no differences were Marina et al.: Milk Metabolites and Feed Efficiency in Sheep found in DMI between the animals classified in the high, medium and low groups during the trial.Consequently, the significant differences observed for RFI and FCR among the groups should be explained by the productivity of these animals.Indeed, this is reinforced by the results obtained in Toral et al. (2023) and Barrio et al. (2023), where higher FE animals were more productive than lower FE animals.It is important to mention that Barrio et al. (2023) analyzed the same data set used in the current study.Therefore, this suggests an important role of the selected features in the production component of FE metrics.

CONCLUSIONS
In this study, a global metabolomic analysis of Assaf ewe milk led to the identification of reduced subsets of features that effectively predicted RFI and FCR.Notably, FCR predictions consistently outperformed RFI predictions across all scenarios, with similar results observed using various ML models.The selected features are associated with metabolic pathways regulating the metabolism of glucose, galactose, fructose, sphingolipids, amino acids, insulin, and thyroid hormones, which directly impact milk production and yield in sheep and other species.These findings highlight the potential of milk biomarkers for predicting FE metrics in dairy sheep, offering a more accessible and cost-effective approach to selecting feed-efficient animals.The direct relationship between the potential functional role of the selected features and milk yield and composition underscores this outcome.The results in the current study showed no significant differences in DMI between the FE groups for FCR and RFI.Consequently, the differences observed for FE between the groups should be explained by productivity.Therefore, implementing the use of these biomarkers could help in the identification of ewes with higher potential feed efficiency.
) equations for sheep.ECM = kg/d of milk yield × [(0.0071 × g/kg of milk fat) + (0.0043 × g/kg of milk protein) + 0.2224 Figure 1.Workflow applied for the predictions of residual feed intake (RFI) and feed conversion ratio (FCR) based on the milk metabolomic signature.

Figure 2 .
Figure2.Discriminating performance of high, medium, and low groups for RFI (A) and FCR (B) obtained by applying partial least squares discriminant analysis (PLS-DA).The ROC curves show the area under the curve (AUC) obtained for discriminating between one group and the others using the first 2 principal components.The X and Y axes in the left-hand side plots represent the percentage of the variance explained by the first and second components identified in the PLS-DA analysis.For right-hand side plots, the x-axis shows the 100-specificity (100%) while the y-axis shows the sensitivity (%) obtained in the discriminant analysis for the 3 feed efficiency groups.

Figure 3 .Figure 4 .
Figure 3. Annotation results for the 41 and 26 features selected for residual feed intake (RFI) and feed conversion ratio (FCR), respectively.The red cells in the heatmap show a putative match between the m/z of the selected features and those available in the annotation databases.The heatmap was divided into blocks based on m/z values to improve the visualization: 91.03917-217.01706(A), 220.11845-520.50812(B), and 567.06116-812.54852(C).The green and purple bars over the heatmaps indicate whether the feature was selected for FCR or RFI.The number of features annotated and shared between FCR and RFI is shown in the Venn diagram (D).
Marina et al.: Milk Metabolites and Feed Efficiency in Sheep

Table 1 .
Marina et al.: Milk Metabolites and FeedEfficiency in Sheep Comparison of the mean (±standard deviation) among the high, medium and low groups of animals for residual feed intake (RFI) and feed conversion ratio (FCR)

Table 3 .
Summary statistics (mean ± standard deviation (sd), minimum (min), and maximum (max)) for the squared Pearson correlation (r 2 ), concordance correlation coefficient (CCC), and root mean squared error (RMSE) values obtained for the predictions of feed conversion ratio (FCR) using the complete milk metabolomic signature (all) and the reduced feature data sets selected based on Variable Importance in the Projection (VIP) obtained in the PLS-DA analysis Deep: Multi-layer feedforward artificial neural network, RF: Random Forest, SVM: Support Vector Machine, xgboost: Extreme Gradient Boosting, OPLS: Orthogonal Partial Least Square.
Marina et al.: Milk Metabolites and Feed Efficiency in Sheep