# Multivariate analysis of milk metabolite measures shows potential for deriving new resilience phenotypes

M. Ithurbide, H. Wang, T. Fassier, Z. Li, J. Pires, T. Larsen, J. Cao, R. Rupp, Nicolas N.C. Friggens

► To cite this version:

M. Ithurbide, H. Wang, T. Fassier, Z. Li, J. Pires, et al. Multivariate analysis of milk metabolite measures shows potential for deriving new resilience phenotypes. Journal of Dairy Science, 2023, 106 (11), pp.8072 - 8086. 10.3168/jds.2023-23332. hal-04317387

HAL Id: hal-04317387

https://hal.inrae.fr/hal-04317387

Submitted on 1 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multivariate analysis of milk metabolite measures shows potential for deriving new resilience phenotypes

M. Ithurbide,<fnref n="1" />* H. Wang,<fnref n="2" /> T. Fassier,<fnref n="3" /> Z. Li,<fnref n="2" /> J. Pires,<fnref n="4" /> T. Larsen,<fnref n="5" /> J. Cao,<fnref n="2" /> R. Rupp,<fnref n="1" /> and N. C. Friggens<fnref n="6" />

## ABSTRACT

In a context of growing interest in breeding more resilient animals, a noninvasive indicator of resilience would be very valuable. We hypothesized that the time-course of concentrations of several milk metabolites through a short-term underfeeding challenge could reflect the variation of resilience mechanisms to such a challenge. We submitted 138 one-year-old primiparous goats, selected for extreme functional longevity (i.e., productive longevity corrected for milk yield [60 low longevity line goats and 78 high longevity line goats]), to a 2-d underfeeding challenge during early lactation. We measured the concentration of 13 milk metabolites and the activity of 1 enzyme during prechallenge, challenge, and recovery periods. Functional principal component analysis summarized the trends of milk metabolite concentration over time efficiently without preliminary assumptions concerning the shapes of the curves. We first ran a supervised prediction of the longevity line of the goats based on the milk metabolite curves. The partial least square analysis could not predict the longevity line accurately. We thus decided to explore the large overall variability of milk metabolite curves with an unsupervised clustering. The large year × facility effect on the metabolite concentrations was precorrected for. This resulted in 3 clusters of goats defined by different metabolic responses to underfeeding. The cluster that showed higher β-hydroxybutyrate, cholesterol, and triacylglycerols increase during the underfeeding challenge was associated with poorer survival compared with the other 2 clusters. These results suggest that multivariate analysis of noninvasive milk measures show potential for deriving new resilience phenotypes.

**Key words:** resilience, milk metabolites, multivariate modeling, dairy goats

## INTRODUCTION

Today, there is growing interest in selecting for resilience, as livestock are expected to face increasingly harsh environmental and climatic conditions. Animal resilience is defined as the ability to overcome short-term environmental disturbances and quickly return to a predisturbance state (Colditz and Hine, 2016). In this context, resilience can be seen as an underlying component of longevity since it corresponds to the ability to cope with and recover from challenges to allow the animal to carry on its productive life (Friggens et al., 2017; Scheffer et al., 2018). Longevity corresponds to true longevity (all culling reasons) and functional longevity that includes all culling reasons, except productivity (Sasaki, 2013). Several studies estimated heritability of functional longevity to be around 10% in cattle and goats (Castañeda-Bustos et al., 2017; Nayeri et al., 2017; Palhière et al., 2018). Ithurbide et al. (2022) showed that selection on functional longevity in a commercial population of dairy goats translated into significant differences in longevity and resilience-related traits such as better mammary health and lower body fat mobilization during the beginning of the first lactation for goats selected for longer functional longevity. Selection seems to be possible; however, improvements are expected to be slow due to low heritability. This low heritability could be explained by the fact that longevity is a multifactorial trait, that is, there are factors other than resilience contributing to longevity, and strong genetic × environmental interactions can be involved (Tsartsianidou et al., 2021). Thus, there is a need to find more direct resilience indicators. Being less multifactorial, more direct resilience indicators could have a higher heritability than functional longevity, and allow a more efficient selection and for instance

<footnote n="1">¹GenPhySE, Université de Toulouse, INRAE, Castanet Tolosan, France 31326</footnote>

<footnote n="2">²Department of Statistics and Actuarial Science, Simon Fraser University, Burnaby BC, Canada V5A 1S6</footnote>

<footnote n="3">³Domaine de Bourges, INRAE, Osmoy, France 78910</footnote>

<footnote n="4">⁴INRAE, Université Clermont Auvergne, Vetagro Sup, UMR Herbivores, Saint-Genès-Champanelle, France 63122</footnote>

<footnote n="5">⁵Department of Animal Science, Aarhus University, 8830 Tjele, Denmark</footnote>

<footnote n="6">⁶UMR 0791 Modélisation Systémique Appliquée aux Ruminants, INRAE, AgroParisTech, Université Paris-Saclay, 75005 Paris, France</footnote>
<table><caption>**Table 1.** Distribution of the 138 goats within the 2 divergent lines selected on high longevity (High_LGV) or low longevity (Low_LGV) bred at INRAE facilities of P3R Bourges and Mosar Paris that underwent the underfeeding challenge during early lactation (36.7 DIM ± 6.2 SD)</caption><thead><tr><th rowspan="2">Item</th><th colspan="5">Year of the underfeeding challenge/INRAE facility</th></tr><tr><th>2020/<br>P3R Bourges</th><th>2021/<br>P3R Bourges</th><th>2021/Paris</th><th>2022/Paris</th><th>Total</th></tr></thead><tbody><tr><td>Low_LGV</td><td>15</td><td>14</td><td>17</td><td>14</td><td>60</td></tr><tr><td>High_LGV</td><td>18</td><td>22</td><td>17</td><td>21</td><td>78</td></tr><tr><td>Total</td><td>33</td><td>36</td><td>34</td><td>35</td><td>138</td></tr></tbody></table>

select animals for longevity at an early stage of productive life.

We hypothesized that the metabolic response to short-term feed restriction could provide information about some genetic characteristics of goat resilience. The objective of this study is to explore the existence of underlying resilience components within the time-course of 13 milk metabolites and 1 enzyme activity during an underfeeding challenge imposed on 2 divergent lines of goats for functional longevity. We propose a new statistical approach to model and explore multi-variate longitudinal data.

## MATERIALS AND METHODS

The experiment was carried out in agreement with French National Regulations for the humane care and use of animals for research purposes. Animals were bred at 2 experimental INRAE Farms: P3R Bourges (UE0332, La Sapinière, Osmoy, France, license to carry out animal experiments: C18-174-01) and Experimental Installation, UMR MoSAR (Route de la Ferme, Thiverval-Grenon, France) close to Paris (license to carry out animal experiments: A 78 615 1002). This article followed the STROBE-Vet guidelines (O'Connor et al., 2016). All procedures performed on animals were approved by the Ethics Committee on Animal Experimentation and the French Ministry of Higher Education, Research and Innovation (APAFIS#8613-2017012013585646 V4 and APAFIS#24314-2019120915403741).

### Animals

Following the method developed by Palhière et al. (2018) and described by Ithurbide et al. (2022), we created 2 functional longevity lines of Alpine goats. Since 2017, we have run the genetic evaluation for functional longevity over 8,787 Alpine AI bucks based on the productive longevity of their daughters (time difference between first kidding and culling) corrected for milk yield. We selected the 16 bucks who had the highest EBV and the 19 bucks who had the lowest EBV among the

whole AI buck population to find the low longevity line (Low_LGV) and high longevity line (High_LGV), respectively. From 2019 to 2022, 138 goats were bred: 60 Low_LGV goats and 78 High_LGV goats. Among them, 69 were bred in the INRAE P3R Bourges facility and 69 in the INRAE Paris facility (Table 1). Within each facility, Low_LGV and High_LGV goats were housed in common pens.

Farm management and animal monitoring in the P3R Bourges facility are described in Ithurbide et al. (2022). Briefly, goats were not culled for milk production reasons, which allows a clean assessment of functional longevity. For farm management reasons, low-producing goats could not be kept on farm in the Paris facility. Thus, survival data were not available in the latter facility. Weight, chest size, and height were measured every month in both facilities during the first year of life. The milk yield (MY), milk fat content (MFC), milk protein content (MPC), and SCS were measured every month during lactation. Moreover, the EBV of the goats sires for functional longevity, MY, MFC, MPC, and SCS were estimated and provided from the national genetic evaluation procedure.

### Underfeeding Challenge

A total of 138 one-year-old primiparous dairy goats were exposed to a 2-d underfeeding challenge during early lactation (35.5 DIM ± 5.6 SD). The design of the underfeeding challenge followed the protocol described in detail in Friggens et al. (2016); briefly, the challenge consisted of a 2-d, straw-only feeding. Milk samples were collected for 4 d prechallenge, throughout the challenge period, and for 4 d following the challenge. From parturition and for 2 wk postchallenge, animals received a standard lactation diet. At P3R Bourges, the goats received a ration based on Lucerne hay offered in collective troughs, complemented with concentrate that was dispensed by automatic concentrate feeders and in the milking parlor. At Paris, the lactation diet was offered as a TMR containing on a DM basis: 20% concentrate, 24% hay, 29% Lucerne, 27% beet pulp,
and 1% mineral and vitamin supplement (as described by M. Gindri [Université Paris-Saclay, INRAE, AgroParisTech, UMR Modélisation Systémique Appliquée aux Ruminants, Palaiseau, France], M. Ithurbide, J. Pires, R. Rupp, L. Puillet [Université Paris-Saclay, INRAE, AgroParisTech, UMR Modélisation Systémique Appliquée aux Ruminants, Palaiseau, France], and N. C. Friggens; unpublished data). Forage and water were offered ad libitum. All goats were milked twice a day. During prechallenge, challenge, and recovery periods, respectively, 3, 2, and 4 milk samples were collected during morning milking. Fixed standard volume were taken after mixing the total production in the milking jar. The concentrations of 13 milk metabolites and 1 enzyme were measured: glucose-6-phosphate (Glu6P), glucose (Glu), galactose (Gal), BHB, isocitrate, glutamate, NH₂ groups, lactate dehydrogenase (LDH), urea, choline, malate, urate, triacylglycerols (TAG), and cholesterol (Chol). Each goat had 13 milk metabolites and 1 enzyme curve with data points at d −7, −4, −1, 0, 1, 2, 3, 4, 5, and 6 for Bourges and every day from d −4 to 12 in the Paris facility. Day 0 was the last morning milking before the underfeeding challenge that started the same day.

Milk urea was analyzed with a FIAstar 5000 Analyzer (Foss Tecator AB, Höganäs, Sweden) using flow injection analysis (Nielsen et al., 2005). Enzymatic-fluorometric methods were used to analyze TAG and minor milk constituents: LDH activity (Larsen, 2005), BHB (Larsen and Nielsen, 2005), urate (Larsen and Moyes, 2010), TAG (Larsen et al., 2011), Chol (Larsen, 2012), isocitrate (Larsen, 2014), Glu, and Glu6P (Larsen, 2015). The Gal in milk was analyzed by an analog procedure to Glu, using β-galactose dehydrogenase (EC 1.1.1.48) to start the fluorometric determination. Moreover, weight, MY, milk composition (MFC, MPC), and udder health indicator (SCS) were measured the same days as the milk samples in both facilities.

## Statistical Analysis

**Milk Metabolite Curve Modeling with Functional Principal Component Analysis.** All statistical analyses were done in the R statistical environment (https://www.r-project.org/). To model the individual metabolite concentration curves, we used the functional data analysis smoothing method described by Ramsay and Silverman (2005). We used a spline interpolation, that is, a piece-wise interpolation that joins several low degree polynomial functions at knots (predetermined time points along the time-series of data). We used natural cubic splines (i.e., a piece-wise cubic polynomial that is a continuous when differentiated twice), fixing a minimum degree of the polynomial at 5. The

degree of smoothing of the spline was controlled by a roughness penalty.

Three goats exhibited outlier metabolic trajectories with BHB concentrations in milk above 3 SD. The recordings of these goats were excluded. In addition, 10 implausible data points were excluded from the analysis (10 out of 27,594 data points). None of these points belonged to the underfeeding period (d 0 to 2) and each belonged to different goats and metabolite curve. As such, removing these points did not distort the general shape of the curves. Figure 1 shows the smoothed curves of the 13 milk metabolites and 1 enzyme from one randomly selected goat.

**Correction for the Year Facility Effect With Functional Regression.** To minimize the impact of nongenetic factors, such as the global environment, on the metabolic response to the underfeeding challenge, we accounted for the facility × year effect by running a functional regression analysis:

$$X_i(t) = \beta_0(t) + \beta_1(t) \times Year\_i + \xi_i(t),$$

where $X_i(t)$ was the milk metabolite curve for the ith goat, $\beta_0(t)$ was the intercept function, $\beta_1(t)$ was a function of time corresponding to the regression coefficient associated with the year-facility effect, $Year\_i$ was a dummy variable corresponding to the year-facility of study of the ith goat, and $\xi_i(t)$ was the residual term. Using a functional regression coefficient allows a correction of the year-facility effect. The corrected individual curves were then estimated as

$$X_{i\_correct}(t) = \beta_0(t) + \xi_i(t).$$

Figure 2 sets out a summary of the milk metabolite curves modeling steps and shows smoothed curves of isocitrate milk concentration before and after correction for year-facility effect. Note that, at Paris, animals were reared from weaning until mid-gestation on 2 different diets, but these were balanced and equally distributed between years and were also equally distributed between clusters in the present analyses, and consequently were ignored.

**Functional Principal Component Analysis.** We characterized milk metabolite curves upon challenge using a functional PCA (FPCA) for each year-facility-corrected milk metabolite (Yao et al., 2005) using the R package “FDA.” Functional PCA is a statistical method for investigating the dominant modes of variation of a functional data set. It allows the time-related variation to be captured in a small number of principal components (see Figure 2). In other words, FPCA decomposes a set of random function $X_j(t)$ from the jth metabolite in the following representation:
<figure><img src="image_2.png" /><figcaption>Figure 1. The smoothed curves of 13 milk metabolites and 1 enzyme from one randomly selected goat. Day 0 corresponds to the beginning of the 2-d underfeeding challenge. Points correspond to observed values. The 13 milk metabolites and 1 enzyme are glucose-6-phosphate (Glu6P, µM), glucose (Glu, µM), galactose (Gal, µM), BHB (µM), isocitrate (µM), glutamate (µM), NH₂ (glutamate micro equivalent), lactate dehydrogenase (LDH, IU), urea (mM), choline (mM), malate (µM), urate (µM), triacylglycerols (TAG, mM), cholesterol (Chol, µM). The 138 goats belonged to 2 divergent lines selected on high longevity or low longevity bred at INRAE facilities of P3R Bourges and Mosar Paris.</figcaption></figure>

$$X_j(t) = \sum_{k=1}^{\infty} \gamma_{jk} \omega_{jk}(t),$$

where $\omega_{jk}$ are orthogonal functions across k knots, that is, the functional principal components (FPC) that are common to all goats, and $\gamma_{jk}$ are the FPC scores that characterize individual curves.

The first step to this decomposition was to estimate the functional principal components $\omega_{jk}$. Let $C_j(s,t)$ be the covariance function of $X_j(t)$, and it corresponds to a self-adjoint and positive semi-definite operator $C_j:L2(\tau) \to L2(\tau)$, where L2 is a Hilbert space. Specifically, L2(τ) denotes the space of square-integrable functions defined on a compact interval τ, equipped with the inner product and norm that make it a Hilbert space. The FPC $\omega_{jk}(t)$ satisfy the following eigen equation:

$$C_j \omega_{jk} = \rho_{jk} \omega_{jk},$$

where $\rho_{jk}$ are the eigenvalues of $C_j$ and $C_j$ gives the following integral transform:

$$C_j \omega_{jk} = \int_{\mathbb{R}} C_j(s,t) \omega_{jk}(s) ds.$$

To obtain the FPC, we could solve the eigen equations for $k = 1, ..., K$ for a fixed K. Equivalently, the solution fits the maximization problem of

$$\max Q_0(\omega_{jk}) = \max \langle \omega_{jk}, C_j \omega_{jk} \rangle,$$

subject to the constraints of $||\omega_{jk}|| = 1$ and $\langle \omega_{jk}, \omega_{jk'} \rangle = 0$ for $k' < k$. Given a set of observed trajectories $x1, ..., xn$ we have the observed covariance function:
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Schematic diagram showing the different stages of analysis of one milk metabolite curves data set (here isocitrate is shown as example). In the raw curves plot (1), smoothed curves plot (2), and curves corrected for year × facility effect (3), each red line corresponds to one goat. The bold blue line corresponds to one randomly chosen goat. The functional principal components (PC) of the functional principal component analysis (PCA) for isocitrate are plotted in 4.1. and the corresponding scores for the randomly chosen goat are shown in 4.2. The functional PC (FPC) scores of the 13 milk metabolites and 1 enzyme are then used to predict the longevity line of the goats with a sparse partial least square discriminant analysis (sPLS-DA; 5A.) and classify the goats within clusters with the same overall metabolic response to underfeeding challenge (5B.).</figcaption></figure>

$$ \hat{C}_j(s,t) = \frac{1}{n} \sum_{i=1}^{n} x_i(s) x_i(t) $$

with its associated operator $\hat{C}_j$. For obtaining estimates $\hat{\omega}_{jk}$, we solve the maximization or eigen equation problem by replacing $C(s,t)$ with its empirical version $\hat{C}(st)$.

The second step was the calculation of the individual FPC scores as the projections of $X_j(t)$ onto the FPC through the inner product:

$$ \gamma_{jk} = X_j, \omega_{jk} = \int_{\tau} X_j(t) \omega_{jk}(t) dt. $$

We chose the minimum number of components that explained at least 90% of the variability. A small number of FPC scores thus characterized each individual milk metabolite curve. The Figure 2.4. represents the 3 FPC for isocitrate and the corresponding scores for one given goat. Note that by construction, the mean values of the FPC equal 0. The first FPC plot was positive across the whole period of time, increasing slightly during challenge. The chosen goat has a negative first score for isocitrate (1isocitrate), indicating that this goat has a low overall isocitrate concentration. The individual curve plotted in Figure 2.3. confirms that. Likewise, we can interpret the negative second isocitrate score value as a high isocitrate concentration compared with pre-
challenge and recovery concentrations, and the negative third score for isocitrate as a low postchallenge/prechallenge concentration ratio compared with other goats.

**Prediction of the Longevity Lines of Goats Based on the Metabolite Curve Model.** For the supervised clustering we used sparse partial least square discriminant analysis (sPLS-DA) to evaluate the ability of milk metabolite curves to distinguish the longevity lines of the 138 goats. This is a linear multi-variate model which performs classification tasks and is able to predict the class of new samples (R package "MixOmics"; Lê Cao et al., 2011). The method integrated a continuous data matrix comprising the individual FPC of the 13 milk metabolites and 1 enzyme and enzyme and a categorical outcome variable: the line of the goat (High_LGV vs. Low_LGV). The sPLS-DA seeks the components that best separate the sample groups, and also selects variables that best discriminate between groups using lasso penalization. We chose the number of components using cross validation on a non-sparse model (comprising all variables) and then tuned the number of variables to select on each component using lasso selection. We assessed the final performance of the model using a 5-fold cross-validation.

## Unsupervised Clustering of the Milk Metabolite Curves

**Hierarchical Clustering.** We ran an unsupervised hierarchical clustering on all the FPC (Dash et al., 2003, R package "FactoMineR"), to define the different types of metabolic responses to the underfeeding challenge for each goat, and independently of the longevity line. First, a 5-dimensions PCA was run on the FPC. The PCA hierarchical clustering starts by treating each goat as a separate cluster. Then, it recursively executes the following steps: (1) identify the 2 closest clusters; (2) merge the 2 closest clusters. This process continues until all the clusters are merged together. The final number of clusters was automatically chosen based on the inertia gain (i.e., finding a minimum number of clusters allowing a low intracluster variability and a high intercluster variability).

To easily understand the differences between clusters we compared the milk metabolite curves between clusters using a permutation test (Ramsay and Silverman, 2005; Sirski, 2012). The test begins by taking the absolute value of a t-test-type statistic at each point along the curve:

$$F(t) = \frac{[\overline{X1}(t) - \overline{X1}(t)]^2}{\sum_i [X1i(t) - \overline{X1}(t)]^2 + \sum_i [X1i(t) - \overline{X2}(t)]^2}$$

Then it uses a permutation test to assess significance, by randomly reordering the curves and recalculating the test statistic with the new groups of curves. We used the default setting of 200 random reorderings. One main advantage of the permutation test is that, unlike parametric tests, it does not assume theoretical probability distributions.

**Data Used to Compare Clusters.** The R package 'survival' was used to compare lifespan between the clusters of goats resulting from the unsupervised clustering. Survival analysis was performed using a Cox model (Cox, 1972). Because of the lack of survival data in Paris data set, the survival analysis was thus run over a sub data set of the 3 clusters, excluding the Paris data.

Analysis of variance tests were used to compare the following data. Weight, chest size, and height at the withers were measured at 6 mo old. The EBVs of the sires of the goats for functional longevity and milk performances were also compared between clusters. The weight and milk performance curves (MY, MFC, MPC, ratio of fat content to protein content [F:P ratio], and SCS) during the underfeeding challenge were compared between clusters following the exact same methodology as the milk metabolite curves: spline interpolation, correction for year × site effect, and permutation test.

# RESULTS

## Modeling of the Individual Milk Metabolite Curves with Functional PCA

The smoothed curves of the 13 milk metabolites and 1 enzyme of one randomly selected goat are presented in the Figure 1, and Table 2 shows the distribution of the milk metabolite concentrations during the whole period of sampling among the 138 goats. After smoothing, the general shape of the curve was preserved and the bounce that occurred during the challenge was correctly fitted. Linear regression correctly corrected for the year-facility effect, leading to similar mean curves per year-facility for each metabolite. Between 2 and 4 functional components were necessary to explain 90% of variation for each metabolite. The interpretation of those principal components should be made as follows (Figure 2): using the isocitrate components as an example, the first component (PC1) of isocitrate roughly corresponds to a flat line over the whole period, a goat with a higher than average milk isocitrate concentration over the whole period will have a proportionally high score FPC on this component. The second component (PC2) shows a positive flat line before challenge and a shift to negative value during the challenge: a goat with higher than average concentration of isocitrate be-
<table><caption>**Table 2.** Concentrations of 13 milk metabolites and 1 enzyme collected during 10 morning milkings among 138 goats that underwent an underfeeding challenge during early lactation; 2 samples were taken during underfeeding challenge of 2 d, 4 d before and 4 d after<fnref n="1" /></caption><thead><tr><th rowspan="2">Item<sup>2</sup></th><th colspan="2">Prechallenge<br>(n = 552 samples)</th><th colspan="2">Challenge<br>(n = 276 samples)</th><th colspan="2">Postchallenge<br>(n = 552 samples)</th></tr><tr><th>Mean</th><th>(SD)</th><th>Mean</th><th>(SD)</th><th>Mean</th><th>(SD)</th></tr></thead><tbody><tr><td>BHB (μM)</td><td>27.3</td><td>(±8.8)</td><td>39.8</td><td>(±25.1)</td><td>25.0</td><td>(±7.5)</td></tr><tr><td>Chol (μM)</td><td>192.0</td><td>(±74.4)</td><td>480.9</td><td>(±192.7)</td><td>258.6</td><td>(±144.3)</td></tr><tr><td>Choline (mM)</td><td>1.5</td><td>(±0.9)</td><td>3.3</td><td>(±1.3)</td><td>1.5</td><td>(±1.0)</td></tr><tr><td>Gal (μM)</td><td>65.9</td><td>(±22.4)</td><td>86.5</td><td>(±36.2)</td><td>66.3</td><td>(±29.9)</td></tr><tr><td>Glu (μM)</td><td>226.8</td><td>(±120.6)</td><td>117.2</td><td>(±57.2)</td><td>269.8</td><td>(±161.1)</td></tr><tr><td>Glu6P (μM)</td><td>60.8</td><td>(±19.0)</td><td>40.9</td><td>(±20.9)</td><td>55.6</td><td>(±20.7)</td></tr><tr><td>Glutamate (μM)</td><td>250.0</td><td>(±95.3)</td><td>103.3</td><td>(±43.3)</td><td>283.2</td><td>(±157.8)</td></tr><tr><td>Isocitrate (μM)</td><td>162.2</td><td>(±45.5)</td><td>239.3</td><td>(±76.5)</td><td>118.7</td><td>(±42.2)</td></tr><tr><td>LDH (IU)</td><td>10.2</td><td>(±4.7)</td><td>45.5</td><td>(±27.8)</td><td>13.9</td><td>(±14.0)</td></tr><tr><td>Malate (μM)</td><td>98.0</td><td>(±51.2)</td><td>38.0</td><td>(±20.5)</td><td>78.6</td><td>(±34.1)</td></tr><tr><td>NH<sub>2</sub> (μEq)</td><td>1,687.1</td><td>(±377.1)</td><td>1,356.5</td><td>(±300.6)</td><td>1,893.5</td><td>(±414.9)</td></tr><tr><td>TAG (mM)</td><td>44.7</td><td>(±15.5)</td><td>83.6</td><td>(±31.3)</td><td>44.4</td><td>(±20.6)</td></tr><tr><td>Urate (μM)</td><td>58.9</td><td>(±31.7)</td><td>116.5</td><td>(±70.7)</td><td>76.2</td><td>(±57.7)</td></tr><tr><td>Urea (mM)</td><td>6.7</td><td>(±3.7)</td><td>6.5</td><td>(±2.1)</td><td>4.8</td><td>(±3.1)</td></tr><tr><td>MY (kg)</td><td>3.1</td><td>(±0.5)</td><td>1.7</td><td>(±0.6)</td><td>2.5</td><td>(±0.6)</td></tr><tr><td>MFC (g/kg)</td><td>40.5</td><td>(±6.2)</td><td>65.9</td><td>(±15.5)</td><td>38.7</td><td>(±10.7)</td></tr><tr><td>MPC (g/kg)</td><td>33.2</td><td>(±2.5)</td><td>34.8</td><td>(±3.9)</td><td>32.8</td><td>(±2.8)</td></tr><tr><td>F:P ratio</td><td>1.3</td><td>(±0.2)</td><td>1.9</td><td>(±0.4)</td><td>1.2</td><td>(±0.3)</td></tr><tr><td>SCS</td><td>4.9</td><td>(±1.8)</td><td>5.6</td><td>(±3.0)</td><td>6.1</td><td>(±1.8)</td></tr></tbody></table>

fore challenge and a lower than average concentration postchallenge will get a high score on this component. In total 48 FPC were attributed to each goat to characterize the variation of the 13 milk metabolites and 1 enzyme through the underfeeding challenge. The FPC were then used to compare the milk metabolite variations between lines.

## Supervised Clustering to Compare the Milk Metabolite Curves Between the 2 Longevity Lines of Goats

The optimal number of components in the sPLS-DA to discriminate the 2 longevity lines of goats was 1. The lasso penalization selected 13 variables for this component (Figure 3). Chief among these were Gal (third FPCA component), glutamate and urea (respectively first and second FPCA component). The balanced error rate (i.e., the percentage of misclassifications) estimated overall for the model was 49.5%. It was 61.3% and 37.7%, respectively, for Low_LGV and High_LGV lines. Thus, the milk metabolite curves during an underfeeding challenge could not predict the longevity line of the goats coming from the 4 year-sites of the experiment. However, if the analysis was run on P3R Bourges and Paris separately, the balanced error (BER) rate was respectively 44 and 37%. Figure 3 shows the contributions of the selected FPC scores to the prediction of

the longevity line within Paris and P3R Bourges data sets. When the analysis was made on each of the 4 year-facilities separately, the BER ranged from 30% (Paris 2022) to 39% (P3R Bourges 2021).

## Unsupervised Clustering on the Milk Metabolite Curves

**Description of the Clusters.** The correlation circle of the PCA applied on the 48 FPC of the 138 goats is shown in Figure 4A. Three clusters were identified: cluster 1, 2, and 3, respectively, gathered 36, 53, and 49 goats (Figure 4B). Distribution of goats from the 4 year-facility combinations did not differ along clusters nor between longevity lines and the number of kids per kidding (*P* > 0.70). The mean milk metabolite curves per cluster and permutation test result are shown in Figure 5. The permutation test over the milk metabolites curves between clusters indicates which milk metabolites were significantly different between clusters. The metabolite curves that were significantly different between clusters are BHB, Chol, choline, Glu, Glu6P, glutamate, LDH, malate, NH₂, and TAG (permutation test, 5% critical value). Except for TAG and Chol, the differences were significant only after the beginning of the feed restriction. For Chol and TAG the values were significantly higher for cluster 2 before d −2.

<footnote n="1">¹The 138 goats belonged to 2 divergent lines selected on high longevity or low longevity bred at INRAE facilities of P3R Bourges and Mosar Paris.</footnote>

<footnote n="2">²Glucose-6-phosphate (Glu6P), glucose (Glu), galactose (Gal), triacylglycerols (TAG), cholesterol (Chol), and lactate dehydrogenase (LDH). Milk performance: daily milk yield (MY), ratio of fat content to protein content (F:P ratio), fat content (MFC), protein content (MPC).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Loading plot from the sparse partial least square discriminant analysis (sPLS-DA) to discriminate the 2 lines selected on functional longevity, based on 13 milk metabolites and 1 enzyme curve through a 2-d underfeeding challenge run separately on goats from the P3R Bourges facility (left, n = 69) and Paris facility (right, n = 69). Colors indicate the longevity line in which the median is maximum for each functional principal component score: red = low longevity (Low_LGV), blue = high longevity (High_LGV). Variable names indicate the functional principal component score of the metabolite that was selected. Glu6P = glucose-6-phosphate; Glu = glucose; Gal = galactose; TAG = triacylglycerols; Chol = cholesterol; LDH = lactate dehydrogenase.</figcaption></figure>

Cluster 1 was mainly characterized by lower Glu, malate, and glutamate during the recovery period; higher Gal and Chol during the recovery period; higher BHB, TAG, and choline during challenge; and lower NH₂ and Glu6P during challenge and early recovery. The cluster 2 was mainly characterized by lower TAG, choline, BHB, LDH, and Chol during challenge. The cluster 3 was mainly intermediate between clusters 1 and 2 with the exception of a higher LDH during challenge.

**Comparison of Resilience-Related Features Between Clusters.** The Kaplan–Meier survival curves of the 3 clusters are displayed in Figure 6. The Cox analysis shows poorer survival of goats belonging to the cluster 1 relative to both clusters 2 and 3 (*P* = 0.04, hazard ratio = 2.63 and *P* = 0.02, hazard ratio = 3.70, respectively, Table 3). Note that a Cox analysis comparing cluster 1 relative to the rest of the goats (i.e., cluster 2 and 3 merged) shows a more significant effect (hazard ratio = 2.97, *P* = 0.009). Table 4 presents ANOVA results between clusters. No significant difference could be seen in sire's EBV (longevity, MY, milk components, and SCS), morphology at 6 mo old, and DIM at the beginning of the challenge (*P* > 0.05). Figure 7 shows the mean curves of the milk performance and weight within the 3 clusters through a 2-d underfeeding challenge. Cluster 1 showed higher F:P ratio and MFC, during and after challenge, as well as

higher MPC during challenge and higher SCS after challenge (permutation test, 5% critical value).

## DISCUSSION

### Context

This study presents an innovative design using longevity lines exposed to a short-term challenge with repeated measures of multiple milk metabolites. We hypothesized that the metabolic responses to a short-term feeding restriction would characterize a resilience mechanism that has an impact on goat survival within herd. Repeated measurements over time were of great value in understanding the temporal aspect of resilience (Döring et al., 2015). Moreover, animal resilience is a complex trait as it involves many interconnected physiological regulations and metabolic pathways. Novel data analysis methods of 13 milk metabolites and 1 enzyme concentration over time allowed us to both grasp the time varying aspect of the process and some of its complexity. Several studies report the modeling of a physiological response to short-term perturbation (Sadoul et al., 2015; Friggens et al., 2016). Those models made strong assumptions concerning the shape of the curves to decipher the different components of the reaction (prechallenge baseline, response, recovery). To deal with the complexity of the metabolic pathways that we explored, and reduce the number of assumptions made a priori, we used spline interpolations as they were flexible and do not make a priori assumptions regarding curve shapes. Both the sparsity of the time points and the heterogeneity of variance between days (a sharp difference occurred during the 2 d of underfeeding challenge) made it difficult to settle on a proper roughness penalty that would be strong enough to prevent boundary effects and flexible enough to capture the bounce during challenge. That is why we used natural cubic splines, fixing a minimum degree of the polynomial at 5. The natural cubic spline is considerably “stiffer” than a polynomial in the sense that it has less tendency to oscillate between data points. Imposing a minimum complexity via the natural cubic spline allowed both a small boundary effect and a good fitting of the sharp increases and decreases during the underfeeding challenge. Moreover, Friggens et al. (2016) showed an interesting variability in the reaction to the challenge but also a strong correlation between the presupposed components of the reaction, suggesting redundancy among them. This is why we decided to use FPCA, which allowed an efficient dimension reduction because each principal component is orthogonal to the others, avoiding any redundancy. In a sense, one can see FPCA as an alternative piece-wise modeling
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Hierarchical clustering on the principal component (PC) analysis of 138 goats described with 48 functional principal component scores (FPC) of 13 milk metabolites and 1 enzyme curve through a 2-d underfeeding challenge. The 138 goats belonged to 2 divergent lines selected on high longevity or low longevity bred at INRAE facilities of P3R Bourges and Mosar Paris. (A) Correlation circle of the principal component analysis applied on the 48 FPC of the 138 goats. (B) Individual plot of the principal component analysis. The colors correspond to the 3 clusters determined by the hierarchical clustering.</figcaption></figure>

because the individual curves can be estimated as the linear combination of the functional principal components weighted by the FPC scores (Figure 2), but with automatically optimized components rather than presupposed components.

## Findings of Supervised Clustering

The prediction of the longevity line of goats by the sPLS-DA was associated with a 49.5% error rate, showing that no discrimination of the genetic line was possible by this approach. There may be several explanations for this high error rate. First, as previously stated, the selection for functional longevity might lead to a large intra-line variability (i.e., many factors affecting the longevity could be selected). Ithurbide et al. (2022) showed that the high longevity line of goats had higher BW and lower F:P ratio in milk at the beginning of the first lactation, suggesting that the better survival of the High_LGV line was linked with lower body fat mobilization. However, the nature of the 2-d underfeeding challenge we imposed in the present study does not exactly mimic the challenges that can be naturally undergone during the beginning of the first lactation. The differences in the metabolic reaction to the early lactation-related energy deficit and a negative energy balance induced by feed restriction have been investigated in dairy cows (Gross and Bruckmaier, 2015).

Moreover, functional longevity is a complex trait, and selection for better longevity can result in animals with different kind of resilience or robustness mecha-

nism (resilience or resistance to diseases for example). The present study only explored one aspect of the resilience: the energy metabolism. This diversity of the possible underlying components of longevity reduces the statistical power of the analysis (some goats could be considered High_LGV because they have good genetic value for disease resilience despite low energy metabolism resilience).

Despite the finding that the 2-d underfeeding challenge we used was shown to induce acute metabolic and production deviations (Friggens et al., 2016), resilience and longevity may reflect a broader range of (short- and long-term) coping mechanisms to a diversity of challenges such as heat waves, behavioral stress, and infectious diseases (Friggens et al., 2022). We found a large batch effect between the 4 year × facility combinations of the study. That finding is corroborated by several studies that showed large farm to farm variability in either the proportion of variance explained or in the panel of dynamic features that best predicted resilience (Adriaens et al., 2020; Krogh et al., 2020; Poppe et al., 2020). We chose to apply a linear functional regression to deal with this batch effect. This functional linear regression relied on the hypothesis that the difference we observe between years was not due to resilience-related differences. This correction was necessary to run an unsupervised clustering, but not for the sPLS-DA. We decided to present the result of the sPLS-DA based on FPCA over the milk metabolite curves corrected for year-facility effect to compare the conclusions of the 2
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Mean curves of the milk 13 milk metabolites and 1 enzyme within the 3 clusters identified by unsupervised clustering in 138 goats through a 2-d underfeeding challenge. These curves are corrected for the year × facility effect with a functional linear regression. The red area indicates the time period during which the variables are significantly different between clusters (permutation test, *P* < 0.05). Chol = cholesterol; Gal = galactose; Glu = glucose; Glu6P = glucose-6-phosphate; LDH = lactate dehydrogenase; TAG = triacylglycerols.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Plot of Kaplan-Meier curve showing survival against time in 69 goats of the 3 clusters identified by unsupervised clustering. The unsupervised clustering defined 3 overall metabolic responses to a 2-d underfeeding challenge based on milk metabolite curves. The 69 goats belonged to 2 divergent lines selected on high longevity or low longevity bred at INRAE facilities of P3R Bourges. The Cox model analysis showed significantly poorer survival of the cluster 1 over cluster 2 and 3 (*P* < 0.05).</figcaption></figure>

approaches based on the same FPC scores. However, the sPLS-DA based on the noncorrected curves did not result in better prediction (results not shown). The prediction run separately on the 2 years of experiments in Paris and the 2 years in Bourges showed better results (respectively 37% and 44% BER). This suggests a possible interaction between the longevity line and the environment, as well as the importance to further study those effects. The housing and the staff differed between the 2 facilities. As explained in the Materials and Methods section, the diets were different between facilities. Feed quality can also vary between years due to the prevailing weather and other factors. The 2 farms had performance levels similar to commercial farms, as described in Ithurbide et al. (2022).

<table><caption>**Table 3.** Hazard ratios (HR) with 95% lower and upper CI from Cox hazard model for culling data in 69 goats of the 3 clusters identified by unsupervised clustering<fnref n="1" /></caption><thead><tr><th>Risk factor</th><th>HR</th><th>CI</th><th>P-value</th></tr></thead><tbody><tr><td>Cluster 1 vs. 2</td><td>2.63</td><td>1.07</td><td>6.25</td><td>0.04</td></tr><tr><td>Cluster 1 vs. 3</td><td>3.70</td><td>1.25</td><td>11.11</td><td>0.02</td></tr><tr><td>Cluster 2 vs. 3</td><td>0.70</td><td>0.24</td><td>2.06</td><td>—</td></tr></tbody></table>

The previous points highlighted that, even if the selection for functional longevity implied differences for several resilience traits (Ithurbide et al., 2022), it did not result in 2 strictly different metabolic responses to the underfeeding challenge but rather to a large variability of response that overlapped between the 2 longevity lines of goats. This led us to explore the diversity of responses to the challenge without any preliminary hypothesis on the level of resilience of the goats (i.e., without taking into account line through the unsupervised clustering of the metabolic responses to the underfeeding challenge).

## Findings of Unsupervised Clustering

The unsupervised clustering based on the FPC scores of all the 13 metabolites and the activity of 1 enzyme was a powerful method to explore the diversity of metabolism responses to underfeeding. This analysis defined 3 clusters of metabolic response to the underfeeding challenge. The survival of goats of cluster 1 was lower than cluster 2 and 3, with an estimated hazard ratio equal to 2.97 (P = 0.009); that is, at any age of life, a goat from the cluster 1 had 2.97 times more risk of being culled than other goats (Cox model analysis).

<table><caption>**Table 4.** The ANOVA testing (LSM) for the difference between the 3 clusters identified by unsupervised clustering in 138 goats through a 2-d underfeeding challenge<fnref n="1" /></caption><thead><tr><th>Item<sup>2</sup></th><th>Number of records</th><th>Cluster 1</th><th>Cluster 2</th><th>Cluster 3</th><th>P-value</th></tr></thead><tbody><tr><td>EBV<sub>sire</sub> for functional longevity (d)</td><td>124</td><td>-1.82</td><td>23.59</td><td>6.30</td><td>—</td></tr><tr><td>EBV<sub>sire</sub>-MPC (g/L)</td><td>124</td><td>0.50</td><td>0.02</td><td>0.39</td><td>—</td></tr><tr><td>EBV<sub>sire</sub>-MFC (g/L)</td><td>124</td><td>-0.04</td><td>0.00</td><td>0.13</td><td>—</td></tr><tr><td>EBV<sub>sire</sub>-SCS</td><td>124</td><td>100.7</td><td>100.6</td><td>98.4</td><td>—</td></tr><tr><td>EBV<sub>sire</sub>-MY (kg)</td><td>124</td><td>15.06</td><td>14.43</td><td>4.75</td><td>—</td></tr><tr><td>Weight_6mo (kg)</td><td>137</td><td>34.0</td><td>33.0</td><td>33.8</td><td>—</td></tr><tr><td>Chest_6mo (cm)</td><td>136</td><td>71.4</td><td>71.2</td><td>71.0</td><td>—</td></tr><tr><td>Height_6mo (cm)</td><td>137</td><td>66.2</td><td>66.7</td><td>66.7</td><td>—</td></tr><tr><td>DIM</td><td>138</td><td>34.7</td><td>35.8</td><td>35.9</td><td>—</td></tr></tbody></table>

<footnote n="1">¹The unsupervised clustering defined 3 overall metabolic responses to a 2-d underfeeding challenge based on milk metabolite curves. The 69 goats belonged to 2 divergent lines selected on high longevity or low longevity bred at INRAE facilities of P3R Bourges.</footnote>

<footnote n="1">¹The ANOVA included the year × site effect. The 138 goats belonged to 2 divergent lines selected on high longevity or low longevity bred at INRAE facilities of P3R Bourges and Mosar Paris.</footnote>

<footnote n="2">²Estimated breeding values of the goats' fathers for functional longevity (EBV<sub>sire</sub> for functional longevity), milk protein content (EBV<sub>sire</sub>-MPC), milk fat content (EBV<sub>sire</sub>-MFC), SCS (EBV<sub>sire</sub>-SCS), and milk yield (EBV<sub>sire</sub>-MY). Morphology at 6 mo of age: weight (Weight_6mo), chest size (Chest_6mo), height at the withers (Height_6mo). Days in milk at the beginning of the challenge (DIM).</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Mean curves of the milk performance and weight within the 3 clusters identified by unsupervised clustering in 138 goats through a 2-d underfeeding challenge: daily milk yield (MY), ratio of fat content to protein content (F:P ratio), fat content (MFC), protein content (MPC), SCS, and weight. These curves are corrected for the year × facility effect with a functional linear regression. The red area indicates the time period during which the variables are significantly different between clusters (permutation test, *P* < 0.05). The 138 goats belonged to 2 divergent lines selected on high longevity or low longevity bred at INRAE facilities of P3R Bourges and Mosar Paris.</figcaption></figure>

It should be noted that survival records were only available for goats in the Bourges facility.

Interestingly, the cluster that was associated with the lowest survival (cluster 1) had the highest milk TAG, choline, Chol, and BHB concentrations during challenge and recovery periods (Figure 6). Milk Chol and choline are shown to be associated with milk TAG (Billa et al., 2020). Increased MFC and BHB during the underfeeding challenge suggests higher body fat mobilization in cluster 1 (Bjerre-Harpøth et al., 2012; Pires et al., 2022). A possible interpretation is that a high body fat mobilization during short-term feed restriction is linked with lower resilience mechanism. This was confirmed by the higher F:P ratio and higher MFC of the cluster 1 during challenge (Figure 7). Interestingly, Ithurbide et al. (2022) showed that Low_LGV goats had higher F:P

ratio during early lactation, indicating a link between resilience and body fat mobilization.

Cluster 1 was also defined by lower Glu and Glu6P during the recovery period and from the beginning of challenge, respectively. Milk G6P is synthesized in the mammary gland from Glu and is a precursor for NADPH via the pentose phosphate pathway that provides reduction equivalents for preventing oxidative stress and also for reductive biosynthesis (Garnsworthy et al., 2006). Several studies report an increased G6P milk concentration during feed restriction (Chaiyabutr et al., 1981; Faulkner and Peaker, 1982; Larsen et al., 2016; Billa et al., 2020). Zachut et al. (2016) suggested that the increase in milk Glu6P concentrations observed at the onset of lactation may be due to activation of the pentose phosphate pathway in mammary epithelial
cells. The Glu6P increase would meet the NADPH requirements for the attenuation of cellular oxidative stress during periods of increased fatty acid oxidation. The lower Glu6P among cluster 1 goats might indicate a lower ability to mitigate oxidative stress. Surprisingly, cluster 1 also presented higher Gal concentration during the recovery period. Similarly to Glu6P, Gal is synthesized in the mammary gland from Glu. That might indicate that Glu is preferably used for Gal synthesis rather than Glu6P among cluster 1 animals, which might increase oxidative stress. Interestingly, Ben Abdelkrim et al. (2023) found that Glu and BHB milk concentrations were part of the most informative milk components for determining membership of clusters of milk metabolite curves through a 2-d underfeeding challenge in late-lactating dairy goats.

The glutamate, malate, and NH₂ decrease during challenge tended to be greater in cluster 1 with a slower increase during recovery. Overall, cluster 1 corresponds to goats that have stronger modifications of milk metabolite concentrations during challenge. The idea that a better resilience is associated with smaller metabolic variations is explored in several articles, for example, lower variation and autocorrelation of the daily MY (Poppe et al., 2020) or the relative height of the MY maximum compared with the MY in late lactation (Arnal et al., 2019). The comparison of the milk composition and SCS between clusters showed that cluster 1 had significantly higher SCS during the recovery period. Milk SCS in goats is an indicator of inflammation and bacterial mastitis (Paape et al., 2001; Luengo et al., 2004; Moroni et al., 2005). Interestingly, cluster 1 also showed higher LDH concentration around d 4 after the beginning of the challenge. Endogenous LDH in milk originates mainly from somatic cells, leucocytes, and invading microorganisms (Larsen, 2005) and is an indicator of inflammation (Krogh et al., 2020). Increased LDH during feed restriction could be partly explained by cell damage of mammary tissue during the challenge period and was also reported in Ben Abdelkrim et al. (2023). Inflammation imposes a metabolic burden because it requires glucose and other limiting nutrients in ruminants, and may explain decreased concentrations of glucogenic milk metabolites concomitant with increased SCS in cluster 1 (Bouvier-Muller et al., 2018; Kvidera et al., 2017). Our study suggests that the goats of the cluster 1 were characterized by lower resilience mechanisms, related both to energy metabolism and the inflammatory system.

## CONCLUSIONS

This study presented the curves of 13 milk metabolites and 1 enzyme through an underfeeding challenge

among 138 early lactating primiparous goats selected for extreme functional longevity. A novel functional PCA approach was used to model the milk metabolite curves, allowing us to address the dynamic and multifactorial patterns of the responses. The approach did not discriminate the 2 longevity lines, highlighting a large variability within lines. Unsupervised clustering of such profiling, however, showed distinct metabolite curves associated with length of productive life in the flock. Moreover, we found that cholesterol, Glu6P, Glu, TAG, and BHB were the most discriminating metabolites for the cluster. These results confirm that multivariate analysis of noninvasive milk measures shows potential for deriving new resilience phenotypes.

## ACKNOWLEDGMENTS

The authors thank the staff of the INRAE experimental unit P3R at Bourges and the Experimental Installation of UMR MoSAR at Paris Grignon, France, for technical support in producing, raising, and monitoring the animals. This work was performed in the framework of the EU H2020 SMARTER project (grant no. 772787, Castanet Tolosan, France) and the national APIS-GENE (Paris, France) funded projects Active-Goat and Resilait. M. Ithurbide acknowledges EU H2020 (SMARTER) and APIS-GENE (Resilait) for her PhD grant. The authors have not stated any conflicts of interest.

## REFERENCES

Adriaens, I., N. C. Friggens, W. Ouweltjes, H. Scott, B. Aernouts, and J. Statham. 2020. Productive life span and resilience rank can be predicted from on-farm first-parity sensor time series but not using a common equation across farms. J. Dairy Sci. 103:7155-7171. https://doi.org/10.3168/jds.2019-17826.

Arnal, M., H. Larroque, H. Leclerc, V. Ducrocq, and C. Robert-Granié. 2019. Genetic parameters for first lactation dairy traits in the Alpine and Saanen goat breeds using a random regression test-day model. Genet. Sel. Evol. 51:43. https://doi.org/10.1186/s12711-019-0485-3.

Ben Abdelkrim, A., M. Ithurbide, T. Larsen, P. Schmidely, and N. C. Friggens. 2023. Milk metabolites can characterize individual differences in animal resilience to a nutritional challenge in lactating dairy goats. Animal 17:100727. https://doi.org/10.1016/j.animal.2023.100727.

Billa, P. A., Y. Faulconnier, T. Larsen, C. Leroux, and J. Pires. 2020. Milk metabolites as noninvasive indicators of nutritional status of mid-lactation Holstein and Montbéliarde cows. J. Dairy Sci. 103:3133-3146. https://doi.org/10.3168/jds.2019-17466.

Bjerre-Harpøth, V., N. C. Friggens, V. M. Thorup, T. Larsen, B. M. Damgaard, K. L. Ingvartsen, and K. M. Moyes. 2012. Metabolic and production profiles of dairy cows in response to decreased nutrient density to increase physiological imbalance at different stages of lactation. J. Dairy Sci. 95:2362-2380. https://doi.org/10.3168/jds.2011-4419.

Bouvier-Muller, J., C. Allain, F. Enjalbert, Y. Farizon, D. Portes, G. Foucras, and R. Rupp. 2018. Somatic cell count-based selection reduces susceptibility to energy shortage during early lactation
in a sheep model. J. Dairy Sci. 101:2248-2259. https://doi.org/10.3168/jds.2017-13479.

Castañeda-Bustos, V. J., H. H. Montaldo, M. Valencia-Posadas, L. Shepard, S. Pérez-Elizalde, O. Hernández-Mendo, and G. Torres-Hernández. 2017. Linear and nonlinear genetic relationships between type traits and productive life in US dairy goats. J. Dairy Sci. 100:1232-1245. https://doi.org/10.3168/jds.2016-11313.

Chaiyabutr, N., A. Faulkner, and M. Peaker. 1981. Changes in the concentrations of the minor constituents of goat's milk during starvation and on refeeding of the lactating animal and their relationship to mammary gland metabolism. Br. J. Nutr. 45:149-157. https://doi.org/10.1079/BJN19810087.

Colditz, I. G., and B. C. Hine. 2016. Resilience in farm animals: Biology, management, breeding and implications for animal welfare. Anim. Prod. Sci. 56:1961-1983. https://doi.org/10.1071/AN15297.

Cox, D. R. 1972. Regression models and life-tables. J. R. Stat. Soc. B 34:187-220.

Dash, M., H. Liu, P. Scheuermann, and K. L. Tan. 2003. Fast hierarchical clustering and its validation. Data Knowl. Eng. 44:109-138. https://doi.org/10.1016/S0169-023X(02)00138-6.

Döring, T. F., A. Vieweger, M. Pautasso, M. Vaarst, M. Finckh, and M. Wolfe. 2015. Resilience as a universal criterion of health. J. Sci. Food Agric. 95:455-465. https://doi.org/10.1002/jsfa.6539.

Faulkner, A., and M. Peaker. 1982. Reviews of the progress of dairy science: Secretion of citrate into milk. J. Dairy Res. 49:159-169. https://doi.org/10.1017/S00220299000224X.

Friggens, N. C., I. Adriaens, R. Boré, G. Cozzi, J. Jurquet, C. Kamphuis, F. Leiber, I. Lora, T. Sakowski, J. Statham, and Y. De Haas. 2022. Resilience: Reference measures based on longer-term consequences are needed to unlock the potential of precision livestock farming technologies for quantifying this trait. Peer Community Journal 2:e38. https://doi.org/10.2407/pcjournal.136.

Friggens, N. C., F. Blanc, D. P. Berry, and L. Puillet. 2017. Review: Deciphering animal robustness. A synthesis to facilitate its use in livestock breeding and management. Animal 11:1-15. https://doi.org/10.1017/S175173111700088X.

Friggens, N. C., C. Duvaux-Ponter, M. P. Etienne, T. Mary-Huard, and P. Schmidely. 2016. Characterizing individual differences in animal responses to a nutritional challenge: Toward improved robustness measures. J. Dairy Sci. 99:2704-2718. https://doi.org/10.3168/jds.2015-10162.

Garnsworthy, P. C., L. L. Masson, A. Lock, and T. Mottram. 2006. Variation of milk citrate with stage of lactation and de novo fatty acid synthesis in dairy cows. J. Dairy Sci. 89:1604-1612. https://doi.org/10.3168/jds.S0022-0302(06)72227-5.

Gross, J. J., and R. M. Bruckmaier. 2015. Repeatability of metabolic responses to a nutrient deficiency in early and mid lactation and implications for robustness of dairy cows. J. Dairy Sci. 98:8634-8643. https://doi.org/10.3168/jds.2014-9246.

Ithurbide, M., C. Huau, I. Palhière, T. Fassier, N. C. Friggens, and R. Rupp. 2022. Selection on functional longevity in a commercial population of dairy goats translates into significant differences in longevity in a common farm environment. J. Dairy Sci. 105:4289-4300. https://doi.org/10.3168/jds.2021-21222.

Krogh, M. A., M. Hostens, M. Salavati, C. Grelet, M. T. Sorensen, D. C. Wathes, C. P. Ferris, C. Marchitelli, F. Signorelli, F. Napoli-tano, F. Becker, T. Larsen, E. Matthews, F. Carter, A. Vanlierde, G. Opsomer, N. Gengler, F. Dehareng, M. A. Crowe, K. L. Ingvar-tsen, and L. Foldager. 2020. Between- and within-herd variation in blood and milk biomarkers in Holstein cows in early lactation. Animal 14:1067-1075. https://doi.org/10.1017/S1751731119002659.

Kvidera, S. K., E. A. Horst, M. Abuajamieh, E. J. Mayorga, M. V. S. Fernandez, and L. H. Baumgard. 2017. Glucose requirements of an activated immune system in lactating Holstein cows. J. Dairy Sci. 100:2360-2374. https://doi.org/10.3168/jds.2016-12001.

Larsen, T. 2005. Determination of lactate dehydrogenase (LDH) activity in milk by a fluorometric assay. J. Dairy Res. 72:209-216. https://doi.org/10.1017/S0022029905000865.

Larsen, T. 2012. Enzymatic-fluorometric quantification of cholesterol in bovine milk. Food Chem. 135:1261-1267. https://doi.org/10.1016/j.foodchem.2012.05.061.

Larsen, T. 2014. Fluorometric determination of free and total isocitrate in bovine milk. J. Dairy Sci. 97:7498-7504. https://doi.org/10.3168/jds.2014-8018.

Larsen, T. 2015. Fluorometric determination of free glucose and glucose 6-phosphate in cows' milk and other opaque matrices. Food Chem. 166:283-286. https://doi.org/10.1016/j.foodchem.2014.06.017.

Larsen, T., L. Alstrup, and M. R. Weisbjerg. 2016. Minor milk constituents are affected by protein concentration and forage digestibility in the feed ration. J. Dairy Res. 83:12-19. https://doi.org/10.1017/S0022029915000692.

Larsen, T., M. K. Larsen, and N. C. Friggens. 2011. Enzymatic and fluorometric determination of triacylglycerols in cow milk and other opaque matrices. Food Chem. 125:1110-1115. https://doi.org/10.1016/j.foodchem.2010.09.104.

Larsen, T., and K. M. Moyes. 2010. Fluorometric determination of uric acid in bovine milk. J. Dairy Res. 77:438-444. https://doi.org/10.1017/S0022029910000580.

Larsen, T., and N. I. Nielsen. 2005. Fluorometric determination of β-hydroxybutyrate in milk and blood plasma. J. Dairy Sci. 88:2004-2009. https://doi.org/10.3168/jds.S0022-0302(05)72876-9.

Lé Cao, K.-A., S. Boitard, and P. Besse. 2011. Sparse PLS discriminant analysis: Biologically relevant feature selection and graphical displays for multiclass problems. BMC Bioinformatics 12:253. https://doi.org/10.1186/1471-2105-12-253.

Luengo, C., A. Sánchez, J. C. Corrales, C. Fernández, and A. Contreras. 2004. Influence of intramammary infection and non-infection factors on somatic cell counts in dairy goats. J. Dairy Res. 71:169-174. https://doi.org/10.1017/S0022029904000019.

Moroni, P., G. Pisoni, G. Ruffo, and P. Boettcher. 2005. Risk factors for intramammary infections and relationship with somatic-cell counts in Italian dairy goats. Prev. Vet. Med. 69:163-173. https://doi.org/10.1016/j.prevetmed.2004.10.013.

Nayeri, S., M. Sargolzaei, M. K. Abo-Ismail, S. Miller, F. Schenkel, S. S. Moore, and P. Stothard. 2017. Genome-wide association study for lactation persistency, female fertility, longevity, and lifetime profit index traits in Holstein dairy cattle. J. Dairy Sci. 100:1246-1258. https://doi.org/10.3168/jds.2016-11770.

Nielsen, N. I., T. Larsen, M. Bjerring, and K. L. Ingvartsen. 2005. Quarter health, milking interval, and sampling time during milking affect the concentration of milk constituents. J. Dairy Sci. 88:3186-3200. https://doi.org/10.3168/jds.S0022-0302(05)73002-2.

O'Connor, A. M., J. Sargeant, I. R. Dohoo, H. N. Erb, M. Cevallos, M. Egger, A. Ersbøll, S. W. Martin, L. R. Nielsen, D. L. Pearl, D. Pfeiffer, J. Sánchez, M. E. Torrence, H. Vigre, C. Waldner, and M. P. Ward. 2016. Explanation and elaboration document for the STROBE-Vet statement: Strengthening the reporting of observational studies in epidemiology-veterinary extension. J. Vet. Intern. Med. 30:1896-1928. https://doi.org/10.1111/jvim.14592.

Paape, M., B. Poutrel, A. Contreras, J. C. Marco, and A. Capuco. 2001. Milk somatic cells and lactation in small ruminants. J. Dairy Sci. 84:E237-E244. https://doi.org/10.3168/jds.S0022-0302(01)70223-8.

Palhière, I., C. Oget, and R. Rupp. 2018. Functional longevity is heritable and controlled by a major gene in French dairy goats. Page 1130 in 11th World Congress on Genetics Applied to Livestock Production, Auckland, New Zealand.

Pires, J. A., T. Larsen, and C. Leroux. 2022. Milk metabolites and fatty acids as noninvasive biomarkers of metabolic status and energy balance in early-lactation cows. J. Dairy Sci. 105:201-220. https://doi.org/10.3168/jds.2021-20465.

Poppe, M., R. F. Veerkamp, M. L. van Pelt, and H. A. Mulder. 2020. Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding. J. Dairy Sci. 103:1667-1684. https://doi.org/10.3168/jds.2019-17290.

Ramsay, J., and B. W. Silverman. 2005. Functional Data Analysis. 2nd ed. Springer Series in Statistics. Springer-Verlag, New York, NY.

Sadoul, B., O. Martin, P. Prunet, and N. C. Friggens. 2015. On the use of a simple physical system analogy to study robustness features in
animal sciences. PLoS One 10:e0137333. https://doi.org/10.1371/journal.pone.0137333.

Sasaki, O. 2013. Estimation of genetic parameters for longevity traits in dairy cattle: A review with focus on the characteristics of analytical models. Anim. Sci. J. 84:449–460. https://doi.org/10.1111/asj.12066.

Scheffer, M., J. E. Bolhuis, D. Borsboom, T. G. Buchman, S. M. W. Gijzel, D. Goulson, J. E. Kammenga, B. Kemp, I. A. van de Lemput, S. Levin, C. M. Martin, R. J. F. Melis, E. H. van Nes, L. M. Romero, and M. G. M. Olde Rikkert. 2018. Quantifying resilience of humans and other animals. Proc. Natl. Acad. Sci. USA 115:11883–11890. https://doi.org/10.1073/pnas.1810630115.

Sirski, M. 2012. On the statistical analysis of functional data arising from designed experiments. PhD Diss. Faculty of Graduate Studies, University of Manitoba, Winnipeg, Manitoba.

Tsartsianidou, V., V. V. Kapsona, E. Sánchez-Molano, Z. Basdagianni, M. J. Carabaño, D. Chatziplis, G. Arsenos, A. Triantafyllidis, and G. Banos. 2021. Understanding the seasonality of performance resilience to climate volatility in Mediterranean dairy sheep. Sci. Rep. 11:1889. https://doi.org/10.1038/s41598-021-81461-8.

Yao, F., H.-G. Müller, and J.-L. Wang. 2005. Functional data analysis for sparse longitudinal data. J. Am. Stat. Assoc. 100:577–590. https://doi.org/10.1198/016214504000001745.

Zachut, M., G. Kra, Y. Portnik, F. Shapiro, and N. Silanikove. 2016. Milk glucose-6-phosphate dehydrogenase activity and glucose-6-phosphate are associated with oxidative stress and serve as indicators of energy balance in dairy cows. RSC Advances 6:65412–65417. https://doi.org/10.1039/C6RA11924G.

ORCIDS

M. Ithurbide https://orcid.org/0000-0002-5987-3194

Z. Li https://orcid.org/0000-0002-7434-5947

J. Pires https://orcid.org/0000-0002-3773-9293

T. Larsen https://orcid.org/0000-0003-3058-5938

R. Rupp https://orcid.org/0000-0003-3375-5816