# Understanding milk and body reserves trajectories and nutrient partitioning in dairy goats through a modelling approach

Alberto Atzori

► To cite this version:

Alberto Atzori. Understanding milk and body reserves trajectories and nutrient partitioning in dairy goats through a modelling approach. 2024, 10.24072/pci.animsci.100234. hal-04711577

HAL Id: hal-04711577

https://hal.inrae.fr/hal-04711577v1

Submitted on 27 Sep 2024

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
# Peer Community Journal

Section: Animal Science

Research article

Published
2024-09-04

**Cite as**
Nicolas Gafsi, Olivier Martin, Fabrice Bidan, Bénédicte Grimard and Laurence Puillet (2024)
*Diversity of performance patterns in dairy goats: multi-scale analysis of the lactation curves of milk yield, body condition score and body weight*, Peer Community Journal, 4: e85.

**Correspondence**
nicolas.gafsi@vet-alfort.fr

**Peer-review**
Peer reviewed and
recommended by
PCI Animal Science,
https://doi.org/10.24072/pci animsci.100234

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

This article is licensed
under the Creative Commons
Attribution 4.0 License.

## Diversity of performance patterns in dairy goats: multi-scale analysis of the lactation curves of milk yield, body condition score and body weight

Nicolas Gafsi<sup>D,1,2</sup>, Olivier Martin<sup>D,1</sup>, Fabrice Bidan<sup>D,2</sup>, Bénédicte Grimard<sup>D,3</sup>, and Laurence Puillet<sup>D,1</sup>

Volume 4 (2024), article e85

https://doi.org/10.24072/pcjournal.449

## Abstract

In the dairy goat sector, reduced longevity is a key issue leading to higher replacement rates in the herd and a poor dilution of doe-rearing costs. There is a need to better understand the determinants of lifetime performance. Thus, the general objective of this work was to analyze the phenotypic variability of lifetime trajectories (milk yield (MY), body weight (BW) and body condition score (BCS)) through a 3-step approach: (1) characterize individual phenotypic lactation curves, (2) explore the associations between MY, BW and BCS curves at the lactation scale and (3) assess the diversity of phenotypic curves over successive lactations. Routine data from two experimental farms: Le Pradel (Dataset 1, Ardeche department, France) and MoSAR experimental farm (Dataset 2, Yvelines department, France) were used. Dataset 1 included 793 Alpine goats from 1996 to 2020. Dataset 2 included 339 Alpine and 310 Saanen goats from 2006 to 2022. Weekly MY records (Dataset 1) and daily MY records (Dataset 2) were fitted using a lactation model with explicit representation of perturbations. Monthly BW records (Dataset 1) and BCS records (Dataset 1&2) were fitted using the Grossman multiphasic model. Daily BW records (Dataset 2) were fitted using a weight model. Each individual lactation curve modeled for MY, BW and BCS was thus summarized by synthetic indicators of level and dynamics. Principal component analysis was performed on the MY, BW and BCS indicators separately, and clusters of phenotypic curves identified. At the lactation scale, associations between MY, BW and BCS clusters were evaluated by contingency tables with a chi-square test. Lifetime-scale bar plots were used to display cluster changes throughout parities. For MY curves, 4 and 3 clusters were found for primiparous and multiparous goats respectively. For BW, lumbar and sternal BCS curves, 3 clusters were found for all parities. At the lactation scale, no major association was found among phenotypic curves suggesting a diversity of energy partitioning strategies between life functions. At the lifetime scale, change among clusters occurred primarily between first and second lactation, whereas a pattern of stable cluster membership appeared for multiparous goats. Further analyses are needed to include reproductive performance in analyzing lifetime performance clusters, to better identify clusters or combinations of clusters at risk for culling.

<footnote n="1">¹Université Paris-Saclay, INRAE, AgroParisTech, UMR Modélisation Systémique Appliquée aux Ruminants, 91120, Palaiseau, France, ²Institut de l'Elevage, 75595 Paris, France, ³Université Paris-Saclay, INRAE, ENVA UMR BREED, 78350, Jouy-en-Josas, France</footnote>
Introduction

The dairy goat sector faces many challenges, such as animals with reduced longevity (Palhière et al., 2018) and high replacement costs. In the future design of livestock farming, breeding and managing robust animals is on the agenda of many research programs. One of the key elements of robustness is to consider goats as a biological system in which productive functions (e.g., lactation, growth, reproduction, etc...) dynamically interact through complex mechanisms involving nutrient partitioning (Bauman & Currie, 1980; Friggens et al., 2017). Nutrient partitioning implies that energy cannot be maximized across all productive functions and therefore some functions are given priority over others, especially to support some physiological stages (e.g. lactation). Thus, individual variability in performance could reveal different nutrient partitioning strategies. A first important aspect to explain changes in nutrient partitioning is the succession of reproductive cycles throughout life. This modifies priorities among functions to support a given physiological stage (e.g., gestation, lactation). In addition to these homeorhetic drivers, priorities can be modified by various aspects of the farming system environment. For instance, it is well documented that genetic selection for milk production has altered priorities among functions in dairy cattle leading to health and reproductive disorders (Pryce et al., 2001; Roche et al., 2009; Friggens et al., 2010). Indeed, high genetic merit for milk has led to energy partitioning in favor of lactation over other biological functions. It is also known that priorities can be modified to cope with nutritional constraints. For instance, most of female mammals will not invest energy in pregnancy during feed shortage (Friggens, 2003). As a central function to support lactation and as a buffer for variation in nutritional environment, body reserves play a central role in energy partitioning among productive functions.

Assessing the diversity of phenotypic lactation curves reflecting productive functions (e.g., milk yield (MY) and body reserves (body weight (BW), body condition score (BCS)) is a way to understand interactions among biological functions and potential trade-offs between them. With time series data based on more frequent measures (e.g., MY, BW, BCS...), the use of mathematical models can provide information about individual phenotypic lactation curves and their variability. Models can be used to transform raw data into biologically meaningful information. Over the past decades, authors have proposed mathematical models to capture the shape of the lactation curve (Wood, 1967; Cobby & Le Du, 1978; Dhanoa, 1981; Wilmink, 1987) and some wanted to have models based on a biological framework (Dijkstra et al., 1997; Friggens et al., 1999; Pollott, 2000). With more frequent data, a recent model was developed to characterize the lactation curve with an explicit representation of perturbations (Ben Abdelkrim et al., 2020). This model allows a better estimation of the lactation potential for a given animal. Having an estimation of the potential lactation curve can help to identify those goats that need improved feeding management (Arnal et al., 2018). Studies on modelling the shape of BW or BCS curves through lactation (Macé et al., 2023) are less frequent. Some mathematical functions with an exponential approach (Sauvant et al., 2012) or a random regression approach (Berry et al., 2003) have been used. In dairy cows, Ollion et al. (2016) developed a method to characterize trade-offs among biological functions. This method was based on principal component analysis (PCA) followed by agglomerative hierarchical classification (AHC) using MY curves, BCS curves and reproductive performance.

Studying the diversity of lactation curve sequences on a lifetime scale opens the perspective to look at potential changes in priorities among functions across the lifespan, and thus to see how early lifetime performance can impact the subsequent productive lifetime. Understanding the career diversity within a herd would allow the development of management strategies adapted to different curve sequence types, thereby favoring animal longevity. To our knowledge, no studies in dairy goats have used models to compare milk, body weight and body condition dynamics at a lactation scale or at a lifetime scale. In this study, we hypothesized that a multi-scale approach (lactation and lifetime scale) based on phenotypic curves would bring insights on energy partitioning strategies among biological functions.
The general objective of this work was to analyze the variability of lifetime phenotypic trajectories
through a 3-step approach: (1) characterize individual lactation curves, (2) explore the associations
between MY, BW and BCS curves at the lactation scale and (3) assess the diversity of phenotypic curves
over successive lactations.

# Material and methods

Ethics approval

This paper did not require animal experimentation approval because the datasets came from routine
data recorded on farm. The two farms, housed their animals in conditions that fully complied with the
current regulations on animal housing (directive 98/58/CE).

Datasets

Dataset 1 (1996-2020).

Data came from the experimental farm Le Pradel (agricultural high school Olivier de Serres) located in the French department Ardeche (44° 34' 58.4364" N; 4° 29' 53.2068" E). The data set contained 2,460 lactations from 793 Alpine goats including 93,965 weekly milk records, 28,099 monthly BW records and 26,271 monthly BCS records. Over this period, goats were milked twice daily, and the recorded milk yield value was the sum of the two milkings. BW was measured once a month on a weighing balance. BCS was evaluated at lumbar and sternal regions on a 0 to 5 scale (Morand-Fehr & Hervieu, 1999). Le Pradel farm had a seasonal system with a kidding period between January and February. During the breeding season in August, inseminated goats received a hormonal treatment. Males were introduced 18 days after artificial insemination (AI). Males stayed until October to mate the goats that returned to heat after AI and those that were not inseminated. Goats produced milk from January to November-December. All lactations retained for milk records had a first record less than 30 days after kidding, a last record after 240 days in milk and had less than 30 days interval between two records. All lactations retained for BW and BCS records had a first record less than 17 days after kidding, a last record after 240 days, more than 8 records per lactation and less than 100 days interval between two records. Lactations lasted on average 289.6 ± 28.5 days. The final dataset 1 concerned 2,271 lactations for milk records, 1,935 lactations for BW records and 1,851 lactations for BCS records (Table 1).

Dataset 2 (2006-2022).

Data came from the MoSAR experimental farm (INRAE-AgroParisTech) located in the French department of Yvelines (48° 50' 31.4801" N; 1° 56' 56.5843" E). The data set contained 1,608 lactations from 339 Alpine and 310 Saanen goats including 396,814 daily milk records, 252,725 daily BW records and 11,525 monthly BCS records. The farm has a rotary parlor with an automatic weighing platform, goats were milked and weighed twice a day. The recorded value for milk was the sum of the two milkings. The recorded value for BW was an average of the two measurements. BCS was assessed as the same way as in dataset 1. The MoSAR experimental farm had a seasonal system with a kidding period between January and February. During the breeding season in August, all goats received a hormonal treatment. Selected goats were inseminated after treatment on a fixed date in August. For the goats that were naturally mated, a male was introduced in small groups of 10-12 goats over 6-7 days. Goats produced milk from January to November-December. All lactations retained for milk records in our dataset had a first record less than 5 days after kidding, a last record after 240 days in milk and had less than 30 days interval between two records. All lactations retained for BW and BCS records had a first record less than 20 days after kidding, a last record after 240 days, more than 8 records per lactation and less than 80 days interval between two records. Lactations lasted on average 280.1± 35.1 days. The final dataset 2
concerned 1,256 lactations for milk records, 1,299 lactations for BW records and 381 lactations for BCS records (Table 1).

<table><caption>**Table 1** - Lactation selection criteria for milk yield, body weight and body condition score records with parity and breed distribution for dataset 1 and 2.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Milk yield</th><th colspan="2">Body weight</th><th colspan="2">Body condition score</th></tr><tr><th>Dataset 1</th><th>Dataset 2</th><th>Dataset 1</th><th>Dataset 2</th><th>Dataset 1</th><th>Dataset 2</th></tr></thead><tbody><tr><td>Lactation stage</td><td>First record</td><td>&lt;30</td><td>&lt;5</td><td>&lt;17</td><td>&lt;20</td><td>&lt;17</td><td>&lt;20</td></tr><tr><td>(d)</td><td>Last record</td><td>&gt;240</td><td>&gt;240</td><td>&gt;240</td><td>&gt;240</td><td>&gt;240</td><td>&gt;240</td></tr><tr><td>Interval between records (d)</td><td></td><td>&lt;30</td><td>&lt;30</td><td>&lt;100</td><td>&lt;80</td><td>&lt;100</td><td>&lt;80</td></tr><tr><td>Record per lactation</td><td></td><td>/</td><td>/</td><td>&gt;=8</td><td>&gt;=8</td><td>&gt;=8</td><td>&gt;=8</td></tr><tr><td>Parity</td><td>Primiparous</td><td>671</td><td>520</td><td>606</td><td>499</td><td>549</td><td>143</td></tr><tr><td></td><td>Multiparous</td><td>1,600</td><td>736</td><td>1,329</td><td>800</td><td>1,302</td><td>238</td></tr><tr><td>Breed</td><td>Alpine</td><td>2,271</td><td>716</td><td>1,935</td><td>742</td><td>1,851</td><td>191</td></tr><tr><td></td><td>Saanen</td><td>0</td><td>540</td><td>0</td><td>557</td><td>0</td><td>190</td></tr><tr><td>Total</td><td></td><td>2,271</td><td>1,256</td><td>1,935</td><td>1,299</td><td>1,851</td><td>381</td></tr></tbody></table>

## Models of individual phenotypic lactation curves

Models were selected according to data frequency.

### Lactation curve fitting of both daily and monthly data (dataset 1 and 2)

The perturbed lactation model proposed by Ben Abdelkrim et al. (2020) was fitted to the MY time-series data (Figure 1).

<figure><img src="image_3.png" /><figcaption>**Figure 1** - Example of daily milk records fitted using the model proposed by Ben Abdelkrim et al. (2020) with empty white circles representing raw data, black bold lines representing the unperturbed lactation model (ULM), black thin lines representing the perturbed lactation model (PLM), and grey dotted lines representing the theoretical Wood model. The ULM trajectory was summarized with synthetic indicators: MY<sub>peak</sub> = highest milk yield value; MY<sub>210</sub> = milk yield value at 210 days; SumMY = sum of daily milk yield values over 250 days; Peak time = time of the highest milk yield value; Persistency = (MY<sub>250</sub>-MY<sub>150</sub>)/MY<sub>150</sub> x100.</figcaption></figure>

This model was designed to decompose lactation dynamics into two components: a theoretical unperturbed lactation curve, and the perturbations in milk yield. This approach was selected to characterize lactation curves corrected for perturbations because it captures a proxy of the lactation
potential. The model used for the unperturbed lactation was a modified version of the Wood model (Wood, 1967) integrating a late lactation decrease. The model was fitted in Scilab (Version 6.1.1, www.scilab.org) using an updated version (Martin, unpublished/personnal communication) of the fitting protocol described in Ben Abdelkrim et al. (2020). For further details about the model and the fitting procedure see Appendix A, section 1.

### BCS curve fitting of monthly data (dataset 1 and 2)

The triphasic model proposed by Grossman et al. (1999) was fitted to monthly BCS time-series data (Figure 2). This model was designed to decompose body condition dynamics into three parts: a depletion phase, a plateau phase and a repletion phase. This model allows characterization of curves with less frequent data (at least five records were needed). The model was fitted using RStudio (version 2023.06.01). For further details about the model and the fitting procedure see Appendix A, section 3.

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

**Figure 2** - Example of monthly sternal body condition records fitted using the model proposed by Grossman et al. (1999) with empty white circles representing raw data, black straight lines representing the fitted curve. This fitted curve was summarized with synthetic indicators: BCS_S<sub>k</sub> = sternal BCS at kidding; BCS_S<sub>min</sub> = minimum sternal BCS; BCS_S<sub>210</sub> = sternal BCS at 210 days; Dep_speed_S<sub>k→30</sub> : sternal BCS depletion speed between kidding and 30 days = (BCS_S<sub>30</sub> - BCS_S<sub>k</sub>)/ 30; Rep_speed_S<sub>180→210</sub>: sternal BCS repletion speed between 180 and 210 days = (BCS_S<sub>210</sub> - BCS_S<sub>180</sub>)/ 30.

### BW curve fitting of daily data

The unperturbed weight model proposed by Martin & Ben Abdelkrim, (2019) was fitted to the daily BW time-series data from dataset 2 (Figure 3). This model was designed to decompose the BW dynamics during a lactation into a sequence of depletion/repletion of BW. This model was built to be flexible and to capture various shapes of BW curves. The model was fitted using RStudio (version 2023.06.01). For further details about the model and the fitting procedure see Appendix A section 2.
<figure><img src="image_5.png" /><figcaption>**Figure 3** - Example of daily body weight records fitted using the model proposed by Martin & Ben Abdelkrim, (2019) with empty white circles representing raw data, black straight lines representing the fitted trajectory. This fitted trajectory was summarized with synthetic indicators: BW<sub>k</sub> = body weight at kidding; BW<sub>min</sub> = minimum body weight; BW<sub>210</sub> = body weight at 210 days; Dep_speed<sub>k→30</sub>: Body weight depletion speed between kidding and 30 days = (BW<sub>30</sub> - BW<sub>k</sub>)/30; Rep_speed<sub>180→210</sub>: Body weight repletion speed between 180 and 210 days = (BW<sub>210</sub> - BW<sub>180</sub>)/30.</figcaption></figure>

The same fitting procedure used for BCS was used to fit monthly BW data from dataset 1.

## Fitting convergence

Non-convergence of the fitting procedure occurred in situations where the model was irrelevant to describe data. Non-convergence of the fitting procedure accounted for 0 % of lactations of the datasets for MY, 3 % of lactations of the datasets for BW and 30 % of lactations of the datasets for lumbar BCS and 22 % of lactations of the datasets for sternal BCS. Modelled curves with extreme features were removed using the Tukey's rule (Tukey, 1977) applied to estimates of model parameters and root mean square error (RMSE) (exclusion of values above the third quartile plus three times the interquartile range). Loss associated to extreme features accounted for 3 % of lactations of the datasets for MY, 7 % of lactations of the datasets for BW and 6 % of lactations of the datasets for lumbar and sternal BCS.

## Synthetic indicators to describe fitted individual phenotypic lactation curves

Finally, we used synthetic indicators derived from modelled curves to describe lactation, BW and BCS curves during lactation. Two types of indicators were used: level indicators were considered to characterize performance at specific times and dynamic indicators were considered to characterize temporal changes in performance (Table 2).

## Statistical analysis

### Clusters of phenotypic curves at lactation scale

All statistical analyses were performed using RStudio (version 2023.06.01). Data and scripts can be found in the repository linked to this manuscript (Gafsi et al., 2023). To characterize groups of individual phenotypic curves, principal component analysis (PCA) was performed on the synthetic indicators of MY, BW and BCS separately. The number of principal components (PC) was based on the cumulative variance. To choose the number of PCs at least 75 % of the total variance was needed. PCA was followed by an agglomerative hierarchical clustering (AHC) based on the retained number of PCs for each of MY, BW and BCS, using Ward's linkage procedure.
<table><caption>**Table 2** - Description of the set of synthetic indicators used to describe fitted individual phenotypic curve for milk yield (MY), body weight (BW) and body condition score (BCS).</caption><thead><tr><th>Curve</th><th>Type<sup>2</sup></th><th>Indicator</th><th>Description and calculation</th><th>Unit</th></tr></thead><tbody><tr><td rowspan="5">Milk production</td><td>L</td><td>SumMY</td><td>Total milk produced between 0 and 250 days of lactation, calculated as the sum of daily milk yield values</td><td>kg</td></tr><tr><td>L</td><td>MY<sub>peak</sub></td><td>Highest daily milk yield reached during lactation</td><td>kg/d</td></tr><tr><td>L</td><td>MY<sub>210</sub></td><td>Daily milk yield at 210 days of lactation</td><td>kg/d</td></tr><tr><td>D</td><td>Peak time</td><td>Lactation time at which the maximum milk yield value is reached</td><td>d</td></tr><tr><td>D</td><td>Persistency</td><td>Rate of decrease of milk production between 150 and 250 days of lactation: (MY<sub>250</sub>-MY<sub>150</sub>/MY<sub>150</sub>) x100</td><td>%</td></tr><tr><td rowspan="5">Body weight</td><td>L</td><td>BW<sub>k</sub></td><td>Daily body weight at kidding</td><td>kg</td></tr><tr><td>L</td><td>BW<sub>min</sub></td><td>Minimum daily body weight reached during lactation</td><td>kg</td></tr><tr><td>L</td><td>BW<sub>210</sub></td><td>Daily body weight value at 210 days of lactation</td><td>kg</td></tr><tr><td>D</td><td>Dep_speed<sub>k</sub><br/>≥30</td><td>Speed of body weight depletion between 0 and 30 days of lactation, calculated as: (BW<sub>30</sub> - BW<sub>k</sub>)/ 30</td><td>kg/d</td></tr><tr><td>D</td><td>Rep_speed<sub>18</sub><br/>0-≥210</td><td>Speed of body weight depletion between 180 and 210 days of lactation, calculated as: (BW<sub>210</sub> - BW<sub>180</sub>)/ 30</td><td>kg/d</td></tr><tr><td rowspan="5">Lumbar or sternal<br/>body condition<br/>score<sup>1</sup></td><td>L</td><td>BCS_X<sub>k</sub></td><td>Lumbar/sternal BCS at kidding</td><td>[0-5] scale</td></tr><tr><td>L</td><td>BCS_X<sub>min</sub></td><td>Minimum lumbar/sternal BCS reached during lactation</td><td>[0-5] scale</td></tr><tr><td>L</td><td>BCS_X<sub>210</sub></td><td>Lumbar/sternal BCS at 210 days of lactation</td><td>[0-5] scale</td></tr><tr><td>D</td><td>Dep_speed_X<sub>k</sub><br/>≥30</td><td>Speed of lumbar/sternal BCS depletion between 0 and 30 days of lactation, calculated as: (BCS_X<sub>30</sub> - BCS_X<sub>k</sub>)/ 30</td><td>[0-5] scale/d</td></tr><tr><td>D</td><td>Rep_speed_X<sub>180-≥210</sub></td><td>Speed of lumbar/sternal BCS depletion between 180 and 210 days of lactation, calculated as: (BCS_X<sub>210</sub> - BCS_X<sub>180</sub>)/ 30</td><td>[0-5] scale/d</td></tr></tbody></table>

Ward's method is a hierarchical procedure that iteratively merges groups of individuals represented by points in a Euclidean space resulting in the smallest increase in the sum of within-group sums of squares. This clustering method produces groups that minimize intra-group dispersion and maximize inter-group dispersion at each binary fusion. Preliminary analysis was conducted including the farming systems, breed, and parities all together. Breed and farming systems did not play a strong role on cluster characterization. Parity played a strong role in cluster characterization for MY and BW. So, we performed a clustering by parity (primiparous vs. multiparous) for MY and BW, whereas we performed a single clustering for all parities together for BCS. The optimal number of clusters was based on the higher relative loss of inertia criteria. Differences between clusters for each synthetic indicator were assessed using a one-way ANOVA followed by a Tukey test.

### At lactation scale, contingency tables between clusters of phenotypic curves

To assess the associations between MY, BW and BCS curves at the lactation scale, we produced two-way contingency tables. After clustering, each lactation was assigned a MY, BW or BCS cluster. A contingency table summarized the conditional frequencies of two clusters (e.g., MY and BW clusters). It was used to assess if a cluster membership for a given phenotypic curve was associated to a particular cluster membership for another phenotypic curve, i.e. it showed how these two clusters were dependent on each other. MY, BW and BCS records concerned different numbers of lactations, so each contingency table (e.g., MY with BW or MY with lumbar BCS) considered different sub-populations. Chi-squared tests were performed to assess for associations, between phenotypic curves. Cramer's V test was performed on significant associations to evaluate the strength of the associations. Cramer's V values ranged from 0 to 1. Values close to 1 indicate a strong association, whereas values close to 0 indicate a weak association.

### At lifetime scale, changes in cluster composition for each parity

To assess the diversity of phenotypic lactation curves at lifetime scale, we produced bar plots of the composition of each cluster for parity *n* in terms of clusters in the next parity *n+1*. With this visual display,

<footnote n="1">¹X stands for lumbar (L) or sternal (S).</footnote>

<footnote n="2">²L = level; D = dynamic.</footnote>
it is possible to characterize if goat's assignment to a cluster is stable across parities (reflecting goats with a stable type of lactation curve across parities) or if assignment to a cluster varies across parities (reflecting goats with various dynamics during their lifetime). Chi-squared tests were performed to assess for associations between lactation curves. Cramer's V test was performed on significant associations to evaluate the strength of the associations.

# Results

## Goodness-of-fit

For the two data sets, the RMSE averaged 5.0% ± 1.9% of the average MY per lactation, 2.7% ± 1.0% of the average BW per lactation, 3.6% ± 1.6% of the average lumbar BCS per lactation, and 3.1% ± 1.3% of the average sternal BCS per lactation.

## Phenotypic lactation curves characterization

For all clusters of MY, BW and BCS, a detailed description of cluster names is given in table 3.

<table><caption>**Table 3** - Detailed description of cluster names for MY, BW and BCS. The upper script letter describes the phenotype (Y: milk yield; W: body weight; LU: lumbar BCS and ST: sternal BCS). The superscript describes if the cluster is for primiparous (p) or multiparous (m) goats. The subscript describes the key feature of the cluster (letter for level and plus or minus sign for dynamics).</caption><thead><tr><th>Phenotypic curve</th><th>Cluster</th><th></th></tr><tr><th></th><th>Primiparous</th><th>Multiparous</th></tr></thead><tbody><tr><td>MY<sup>1</sup></td><td>Y<sup>P<sub>L</sub></sup> = Low MY and high persistency.<br>Y<sup>P<sub>L</sub></sup> = Low MY and low persistency.<br>Y<sup>P<sub>M</sub></sup> = Medium MY and low persistency.<br>Y<sup>P<sub>H</sub></sup> = High MY and medium persistency.</td><td>Y<sup>m<sub>M+</sub></sup> = Medium MY and high persistency.<br>Y<sup>m<sub>M</sub></sup> = Medium MY and low persistency.<br>Y<sup>m<sub>H</sub></sup> = High MY and medium persistency.</td></tr><tr><td>BW</td><td>W<sup>P<sub>L</sub></sup> = Low BW and low depletion.<br>W<sup>P<sub>H</sub></sup> = High BW and high depletion.<br>W<sup>P<sub>H</sub></sup> = High BW and low depletion.</td><td>W<sup>m<sub>L</sub></sup> = Low BW and low depletion.<br>W<sup>m<sub>H</sub></sup> = High BW and high depletion.<br>W<sup>m<sub>H</sub></sup> = High BW and low depletion.</td></tr><tr><td>Lumbar BCS</td><td colspan="2">All parities<br>LU<sup>M+</sup> = Medium lumbar BCS and depletion.<br>LU<sup>M</sup> = Medium lumbar BCS and low depletion.<br>LU<sup>H+</sup> = High lumbar BCS and depletion.</td></tr><tr><td>Sternal BCS</td><td colspan="2">ST<sup>M+</sup> = Medium sternal BCS and depletion.<br>ST<sup>M</sup> = Medium sternal BCS and low depletion.<br>ST<sup>H+</sup> = High sternal BCS and depletion.</td></tr></tbody></table>

## Clusters of MY lactation curves

The first two PCs accounted for 83.5% of the total variance for primiparous goats and 81.6% for multiparous goats. The first PC captured the total amount of milk produced during the lactation and accounted for 53.1% of the total variance for primiparous goats and 50.7% for multiparous goats. The second PC captured the persistency and peak time of the lactation curve and accounted for 30.4% of the total variance for primiparous goats and 30.8% for multiparous goats. Based on the highest loss of inertia, four clusters were retained for primiparous goats, and three clusters were retained for multiparous goats (Figure 4).

<footnote n="1">¹Abbreviations: MY = milk yield; BW = body weight; BCS = body condition score</footnote>
<figure><img src="image_6.png" /></figure>

**Figure 4** - PCA and clusters of milk yield synthetic indicators in primiparous (a) and multiparous (b) goats with grey circles representing raw data, lines representing the mean cluster and dotted lines representing a paragon cluster (i.e., the most representative goat in the cluster) (MY<sub>peak</sub> = highest milk yield value; MY<sub>210</sub> = milk yield value at 210 days; SumMY= sum of daily milk yield values over 250 days; Peak time = time of the highest milk yield value; Persistency = (MY<sub>250</sub>-MY<sub>150</sub>)/MY<sub>150</sub>) x100; Y<sup>P</sup><sub>L+</sub>= Low milk yield and high persistency cluster for primiparous; Y<sup>P</sup><sub>L-</sub>= Low milk yield and low persistency cluster for primiparous; Y<sup>P</sup><sub>M-</sub> = Medium milk yield and low persistency cluster for primiparous; Y<sup>P</sup><sub>M+</sub> = High milk yield and medium persistency cluster for primiparous; Y<sup>m</sup><sub>M+</sub>= Medium milk yield and high persistency cluster for multiparous; Y<sup>m</sup><sub>M-</sub> = Medium milk yield and low persistency cluster for multiparous; Y<sup>m</sup><sub>H</sub> = High milk yield and medium persistency cluster for multiparous).

Full details for each cluster are given in tables 4a and 4b.

Primiparous clusters were characterized by:

- a group of low persistency clusters with two different total milk production levels (63.3% of the primiparous): a low-level cluster (Y<sup>P</sup><sub>L</sub>) that produced 155.6 kg less over the lactation than a medium-level cluster (Y<sup>P</sup><sub>M</sub>).

- a medium persistency cluster with the highest total milk production level that gathered 22.6% of the primiparous (Y<sup>P</sup><sub>H</sub>).

- the highest persistency cluster with a low total milk production level that gathered 14.1% of the primiparous (Y<sup>P</sup><sub>L+</sub>).
<table><caption>**Table 4a** - Statistical description of synthetic indicators for MY clusters in primiparous goats.</caption><thead><tr><th>Indicator</th><th>Y<sup>P<sub>L</sub></sup><sub>n = 273</sub></th><th>Y<sup>P<sub>L</sub></sup><sub>n = 163</sub></th><th>Y<sup>P<sub>M</sub></sup><sub>n = 459</sub></th><th>Y<sup>P<sub>H</sub></sup><sub>n = 262</sub></th><th>Pooled SE</th><th>p-value<sup>2</sup></th></tr></thead><tbody><tr><td>SumMY<sup>1</sup></td><td>629.1<sup>a</sup></td><td>675.5<sup>b</sup></td><td>784.7<sup>c</sup></td><td>925.4<sup>d</sup></td><td>67.2</td><td>***</td></tr><tr><td>MY<sub>peak</sub></td><td>3.0<sup>a</sup></td><td>3.0<sup>a</sup></td><td>3.7<sup>b</sup></td><td>4.2<sup>c</sup></td><td>0.4</td><td>***</td></tr><tr><td>MY<sub>210</sub></td><td>2.1<sup>a</sup></td><td>2.7<sup>b</sup></td><td>2.6<sup>b</sup></td><td>3.4<sup>c</sup></td><td>0.3</td><td>***</td></tr><tr><td>Peak time</td><td>47.4<sup>a</sup></td><td>106.0<sup>b</sup></td><td>49.8<sup>a</sup></td><td>71.4<sup>c</sup></td><td>26.7</td><td>***</td></tr><tr><td>Persistency</td><td>-36.7<sup>a</sup></td><td>-19.2<sup>b</sup></td><td>-35.2<sup>a</sup></td><td>-27.2<sup>c</sup></td><td>10.9</td><td>***</td></tr></tbody></table>

a-d Means with superscripts differ significantly by row.

¹SumMY = sum of daily milk yield values over 250 days; MY<sub>peak</sub> = highest milk yield value; MY<sub>210</sub> = milk yield value at 210 days; Peak time = time of the highest milk yield value; Persistency = (MY<sub>250</sub> - MY<sub>150</sub>)/MY<sub>150</sub>) x100.

²p-value resulting from Tukey's test assessing the significance of differences between profiles for each variable. NS (p<0.1), *(p<0.05); and ***(p≤0.001).

³Y<sup>P<sub>L</sub></sup> = Low milk yield and high persistency cluster; Y<sup>P<sub>L</sub></sup> = Low milk yield and low persistency cluster; Y<sup>P<sub>M</sub></sup> = Medium milk yield and low persistency cluster; Y<sup>P<sub>H</sub></sup> = High milk yield and medium persistency cluster.

Multiparous clusters were characterized by:

- a group of medium total milk production levels with two different persistency (65.4 % of the multiparous): a high persistency cluster (Y<sup>mM+</sup>) that maintained 20.4 % more the production than a low persistency cluster (Y<sup>mM-</sup>).

- the highest total milk production level cluster with a medium persistency (Y<sup>mH</sup>) that gathered 34.6 % of the population.

<table><caption>**Table 4b** - Statistical description of synthetic indicators for MY clusters in multiparous goats.</caption><thead><tr><th>Indicator</th><th>Y<sup>mM+</sup><sub>n = 741</sub></th><th>Y<sup>mM-</sup><sub>n = 740</sub></th><th>Y<sup>mH</sup><sub>n = 783</sub></th><th>Pooled SE</th><th>p-value<sup>2</sup></th></tr></thead><tbody><tr><td>SumMY<sup>1</sup></td><td>911.4<sup>a</sup></td><td>940.9<sup>b</sup></td><td>1,212.4<sup>c</sup></td><td>111.5</td><td>***</td></tr><tr><td>MY<sub>peak</sub></td><td>4.1<sup>a</sup></td><td>4.7<sup>b</sup></td><td>5.7<sup>c</sup></td><td>0.6</td><td>***</td></tr><tr><td>MY<sub>210</sub></td><td>3.4<sup>a</sup></td><td>2.9<sup>b</sup></td><td>4.1<sup>c</sup></td><td>0.5</td><td>***</td></tr><tr><td>Peak time</td><td>71.1<sup>a</sup></td><td>38.1<sup>b</sup></td><td>58.3<sup>c</sup></td><td>27.4</td><td>***</td></tr><tr><td>Persistency</td><td>-25.9<sup>a</sup></td><td>-46.3<sup>b</sup></td><td>-36.3<sup>c</sup></td><td>12.4</td><td>***</td></tr></tbody></table>

a-c Means with superscripts differ significantly by row.

¹SumMY = sum of daily milk yield values over 250 days; MY<sub>peak</sub> = highest milk yield value; MY<sub>210</sub> = milk yield value at 210 days; Peak time = time of the highest milk yield value; Persistency = (MY<sub>250</sub> - MY<sub>150</sub>)/MY<sub>150</sub>) x100.

²p-value resulting from Tukey's test assessing the significance of differences between profiles for each variable. NS (p<0.1), *(p<0.05); and ***(p≤0.001).

³Y<sup>mM+</sup> = Medium milk yield and high persistency cluster; Y<sup>mM-</sup> = Medium milk yield and low persistency cluster; Y<sup>mH</sup> = High milk yield and medium persistency cluster.

With respect to farm, for primiparous goats, Pradel's Alpine goats, Grignon's Alpine goats, and Grignon's Saanen goats were more represented in the Y<sup>P<sub>M</sub></sup> cluster because this is the cluster with the highest number of goats overall. For multiparous goats, Pradel's Alpine goats were more represented in the Y<sup>mH</sup> cluster, whereas Grignon's Alpine goats were less represented in this cluster. Grignon's Saanen goats were more represented in the Y<sup>mM+</sup> cluster. See Appendix B section 1 for more details.

### Clusters of BW lactation curves

The first two PCs accounted for 77.4% of the total variance for primiparous goats and 79.4% for multiparous goats. The first PC represented the level of BW at different times of lactation and accounted for 52.8% of the total variance for primiparous goats and 56.9% for multiparous goats. The second PC represented the BW speed loss in the 30 days after kidding and accounted for 24.7% of the total variance for primiparous goats and 22.5% for multiparous goats. Three clusters were retained for each parity group due to the highest loss of inertia (Figure 5).
<figure><img src="image_7.png" /></figure>

**Figure 5** - PCA and clusters of body weight synthetic indicators in primiparous (a) and multiparous (b) goats with grey circles representing raw data, lines representing the mean cluster and dotted lines a paragon cluster (i.e., the most representative goat in the cluster) (BW<sub>k</sub> = body weight at kidding; BW<sub>min</sub> = minimum body weight; BW<sub>210</sub> = body weight at 210 days; Dep_speed<sub>k→30</sub>: Body weight depletion speed between kidding and 30 days = (BW<sub>30</sub> - BW<sub>k</sub>)/30; Rep_speed<sub>180→210</sub>: Body weight repletion speed between 180 and 210 days = (BW<sub>210</sub> - BW<sub>180</sub>)/30; W<sup>m</sup><sub>L</sub> = Low body weight and low depletion cluster in primiparous; W<sup>P</sup><sub>H+</sub> = High body weight and high depletion cluster in primiparous; W<sup>mL</sup> = Low body weight and low depletion cluster in multiparous; W<sup>mH+</sup> = High body weight and high depletion cluster in multiparous; W<sup>mH</sup> = High body weight and low depletion cluster in multiparous).

Full details for each cluster are given in tables 5a and 5b.

Primiparous clusters were characterized by:

- a group of low depletion clusters with two different BW level at kidding (68.6% of the primiparous): a low-level cluster (W<sup>P</sup><sub>L</sub>) that averaged 10.0 kg less at kidding than a high-level cluster (W<sup>P</sup><sub>H</sub>). Those profiles had a higher BW<sub>210</sub> than BW<sub>k</sub>.

- the highest depletion cluster with a high BW level at kidding (W<sup>P</sup><sub>H+</sub>) that gathered 31.4% of the population. Despite having the highest repletion speed, this cluster presented a lower BW<sub>210</sub> than BW<sub>k</sub> due to the high level of depletion, that is not totally compensated at 210 days of lactation.
<table><caption>**Table 5a** - Statistical description of synthetic indicators for BW clusters in primiparous goats.</caption><thead><tr><th>Indicator</th><th>W<sup>P<sub>L</sub></sup><sup>3</sup><br>n = 418</th><th>W<sup>P<sub>H+</sub></sup><br>n = 312</th><th>W<sup>P<sub>H+</sub></sup><br>n = 264</th><th>Pooled SE</th><th>p-value<sup>2</sup></th></tr></thead><tbody><tr><td>BW<sub>k</sub><sup>1</sup></td><td>47.7<sup>a</sup></td><td>54.3<sup>b</sup></td><td>57.7<sup>c</sup></td><td>4.0</td><td>***</td></tr><tr><td>BW<sub>min</sub></td><td>45.2<sup>a</sup></td><td>47.6<sup>b</sup></td><td>55.6<sup>c</sup></td><td>3.5</td><td>***</td></tr><tr><td>BW<sub>210</sub></td><td>49.5<sup>a</sup></td><td>52.9<sup>b</sup></td><td>61.5<sup>c</sup></td><td>4.3</td><td>***</td></tr><tr><td>Dep_speed<sub>k→30</sub></td><td>-0.05<sup>a</sup></td><td>-0.17<sup>b</sup></td><td>-0.03<sup>c</sup></td><td>0.07</td><td>***</td></tr><tr><td>Rep_speed<sub>180→210</sub></td><td>0.04<sup>a</sup></td><td>0.06<sup>b</sup></td><td>0.05<sup>c</sup></td><td>0.03</td><td>***</td></tr></tbody></table>

<sup>a-c</sup> Means with superscripts differ significantly by row.

¹ BW<sub>k</sub> = body weight at kidding; BW<sub>min</sub> = minimum body weight; BW<sub>210</sub> = body weight at 210 days; Dep_speed<sub>k→30</sub> = (BW<sub>30</sub> - BW<sub>k</sub>)/30; Rep_speed<sub>180→210</sub> = (BW<sub>210</sub> - BW<sub>180</sub>)/30.

² p-value resulting from Tukey's test assessing the significance of differences between profiles for each variable. NS (p<0.1), *(p<0.05); and ***(p≤0.001).

³ W<sub>P<sub>L</sub></sub> = Low body weight and low depletion cluster; W<sub>P<sub>H+</sub></sub> = High body weight and high depletion cluster; W<sub>P<sub>H+</sub></sub> = High body weight and low depletion cluster.

Multiparous clusters were characterized by:

- a group of low depletion clusters with two different BW level at kidding (73.4 % of the multiparous): a low-level (W<sup>m<sub>L</sub></sup>) that averaged 17.6 kg less at kidding than a high-level cluster (W<sup>m<sub>H+</sub></sup>). For these clusters BW<sub>210</sub> was lower than BW<sub>k</sub>.

- the highest depletion cluster with a high BW level at kidding (W<sup>m<sub>H+</sub></sup>) that gathered 26.6% of the multiparous. Despite having the highest repletion speed, this profile presented a lower BW<sub>210</sub> than BW<sub>k</sub> due to the high level of depletion, that is not totally compensated at 210 days of lactation.

<table><caption>**Table 5b** - Statistical description of synthetic indicators for BW clusters in multiparous goats.</caption><thead><tr><th>Indicator</th><th>W<sup>m<sub>L</sub></sup><sup>3</sup><br>n = 835</th><th>W<sup>m<sub>H+</sub></sup><br>n = 513</th><th>W<sup>m<sub>H+</sub></sup><br>n = 583</th><th>Pooled SE</th><th>p-value<sup>2</sup></th></tr></thead><tbody><tr><td>BW<sub>k</sub><sup>1</sup></td><td>64.1<sup>a</sup></td><td>78.1<sup>b</sup></td><td>81.7<sup>c</sup></td><td>6.2</td><td>***</td></tr><tr><td>BW<sub>min</sub></td><td>57.8<sup>a</sup></td><td>65.5<sup>b</sup></td><td>74.2<sup>c</sup></td><td>5.2</td><td>***</td></tr><tr><td>BW<sub>210</sub></td><td>61.1<sup>a</sup></td><td>69.4<sup>b</sup></td><td>76.2<sup>c</sup></td><td>5.4</td><td>***</td></tr><tr><td>Dep_speed<sub>k→30</sub></td><td>-0.14<sup>a</sup></td><td>-0.35<sup>b</sup></td><td>-0.14<sup>a</sup></td><td>0.12</td><td>***</td></tr><tr><td>Rep_speed<sub>180→210</sub></td><td>0.04<sup>a</sup></td><td>0.04<sup>b</sup></td><td>0.01<sup>c</sup></td><td>0.03</td><td>***</td></tr></tbody></table>

<sup>a-c</sup> Means with superscripts differ significantly by row.

¹ BW<sub>k</sub> = body weight at kidding; BW<sub>min</sub> = minimum body weight; BW<sub>210</sub> = body weight at 210 days; Dep_speed<sub>k→30</sub> = (BW<sub>30</sub> - BW<sub>k</sub>)/30; Rep_speed<sub>180→210</sub> = (BW<sub>210</sub> - BW<sub>180</sub>)/30.

² p-value resulting from Tukey's test assessing the significance of differences between profiles for each variable. NS (p<0.1), *(p<0.05); and ***(p≤0.001).

³ W<sub>P<sub>L</sub></sub> = Low body weight and low depletion cluster; W<sub>P<sub>H+</sub></sub> = High body weight and high depletion cluster; W<sub>P<sub>H+</sub></sub> = High body weight and low depletion cluster.

For primiparous goats, Pradel's Alpine goats were more represented in the W<sub>P<sub>L</sub></sub> and W<sub>P<sub>H+</sub></sub> clusters. Grignon's Alpine goats were more represented in the W<sub>P<sub>L</sub></sub> cluster. Grignon's Saanen goats were more represented in the W<sub>P<sub>H+</sub></sub> cluster. For multiparous goats, Pradel's Alpine goats and Grignon's Alpine goats were more represented in the W<sub>P<sub>L</sub></sub> cluster. Grignon's Saanen goats were more represented in the W<sub>P<sub>H+</sub></sub> cluster. See Appendix B section 2 for more details.

### Clusters of BCS lactation curves

For lumbar and sternal BCS, clusters were built all parities together. For lumbar BCS, the first two PCs accounted for 75.8% of the total variance. The first PC represented levels of lumbar score at different times of the lactation (BCS_L<sub>min</sub> and BCS_L<sub>k</sub>) and accounted for 46.9% of the total variance. The second PC represented the lumbar BCS speed loss in the 30 days after kidding and accounted for 28.9% of the
total variance. Three clusters were retained due to the highest loss of inertia. For sternal BCS, the first two PCs represented 78.6% of the total variance. The first PC represented levels of sternal score at different times of the lactation (BCS_Smin and BCS_S210) and accounted for 50.7% of the total variance. The second PC represented the sternal BCS speed loss in the 30 days after kidding and accounted for 27.9% of the total variance. Three clusters were retained due to the highest loss of inertia (Figure 6).

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

**Figure 6** - PCA and clusters of lumbar (a) and sternal (b) body condition score synthetic indicators with grey circles representing raw data, lines representing the mean cluster and dotted lines a paragon cluster (i.e., the most representative goat in the cluster) (BCS_Lk = lumbar BCS at kidding; BCS_Lmin = minimum lumbar BCS; BCS_L210 = lumbar BCS at 210 days; Dep_speed_Lk→30: lumbar BCS depletion speed between kidding and 30 days = (BCS_L30 - BCS_Lk)/ 30; Rep_speed_L180→210: lumbar BCS repletion speed between 180 and 210 days = (BCS_L210 - BCS_L180)/ 30;BCS_Sk = sternal BCS at kidding; BCS_Smin = minimum sternal BCS; BCS_S210 = sternal BCS at 210 days; Dep_speed_Sk→30: sternal BCS depletion speed between kidding and 30 days = (BCS_S30 - BCS_Sk) / 30; Rep_speed_S180→210: sternal BCS repletion speed between 180 and 210 days = (BCS_S210 - BCS_S180) / 30; LU_M+ = Medium lumbar body condition score and depletion cluster; LU_M = Medium lumbar body condition score and low depletion cluster; LU_H+ = High lumbar body condition score and depletion cluster; ST_M+ = Medium sternal body condition score and depletion cluster; ST_M = Medium sternal body condition score and low depletion cluster; ST_H+ = High sternal body condition score and depletion cluster).

Full details for each cluster are given in tables 6 and 7.

Lumbar BCS clusters were characterized by:
- a group of depletion clusters with two different lumbar BCS level at kidding (68.7 % of the population): a medium level cluster (LU<sub>M+</sub>) that averaged 0.4 points less at kidding than a high-level cluster (LU<sub>H+</sub>). LU<sub>M+</sub> profile presented the highest repletion speed and the lowest minimum lumbar BCS value.

- the lowest depletion cluster with a medium lumbar BCS level at kidding that gathered 31.3% of the population (LU<sub>M</sub>). LU<sub>M</sub> cluster presented the same repletion speed than LU<sub>H+</sub>.

<table><caption>Table 6 - Statistical description of synthetic indicators for lumbar BCS clusters in goats.</caption><thead><tr><th>Indicator</th><th>LU<sub>M+</sub><sup>3</sup><br>n = 437</th><th>LU<sub>M</sub><br>n = 459</th><th>LU<sub>H+</sub><br>n = 572</th><th>Pooled SE</th><th>p-value<sup>2</sup></th></tr></thead><tbody><tr><td>BCS_L<sub>k</sub><sup>1</sup></td><td>2.5<sup>a</sup></td><td>2.4<sup>b</sup></td><td>2.9<sup>c</sup></td><td>0.2</td><td>***</td></tr><tr><td>BCS_L<sub>min</sub></td><td>2.1<sup>a</sup></td><td>2.3<sup>b</sup></td><td>2.6<sup>c</sup></td><td>0.2</td><td>***</td></tr><tr><td>BCS_L<sub>210</sub></td><td>2.3<sup>a</sup></td><td>2.5<sup>b</sup></td><td>2.7<sup>c</sup></td><td>0.2</td><td>***</td></tr><tr><td>Dep_speed_L<sub>k→30</sub></td><td>-0.009<sup>a</sup></td><td>0.002<sup>b</sup></td><td>-0.006<sup>c</sup></td><td>0.005</td><td>***</td></tr><tr><td>Rep_speed_L<sub>180→210</sub></td><td>0.002<sup>a</sup></td><td>0.001<sup>b</sup></td><td>0.001<sup>b</sup></td><td>0.001</td><td>***</td></tr></tbody></table>

a-c Means with superscripts differ significantly by row.

¹ BCS_L<sub>k</sub> = lumbar BCS at kidding; BCS_L<sub>min</sub> = minimum lumbar BCS; BCS_L<sub>210</sub> = lumbar BCS at 210 days; Dep_speed_L<sub>k→30</sub> = (BCS_L<sub>30</sub> - BCS_L<sub>k</sub>)/30; Rep_speed_L<sub>180→210</sub> = (BCS_L<sub>210</sub> - BCS_L<sub>180</sub>)/30.

² p-value resulting from Tukey's test assessing the significance of differences between profiles for each variable. NS (p<0.1), *(p<0.05); and ***(p≤0.001).

³ LU<sub>M+</sub> = Medium lumbar body condition score and depletion cluster; LU<sub>M</sub> = Medium lumbar body condition score and low depletion cluster; LU<sub>H+</sub> = High lumbar body condition score and depletion cluster

Sternal BCS profiles were characterized by:

- a group of depletion clusters with two different sternal BCS level at kidding (56.5 % of the population): a medium-level cluster (ST<sub>M+</sub>) that averaged 0.7 points less at kidding than a high-level cluster (ST<sub>H+</sub>). ST<sub>M+</sub> cluster presented the lowest minimum sternal BCS. These clusters presented the highest and the same repletion speed.

- the lowest depletion cluster with a medium sternal BCS level at kidding that gathered 43.5 % of the population (ST<sub>M</sub>). ST<sub>M</sub> cluster presented the lowest repletion speed.

<table><caption>Table 7 - Statistical description synthetic indicators for sternal BCS clusters in goats.</caption><thead><tr><th>Indicator</th><th>ST<sub>M+</sub><sup>3</sup><br>n = 489</th><th>ST<sub>M</sub><br>n = 708</th><th>ST<sub>H+</sub><br>n = 433</th><th>Pooled SE</th><th>p-value<sup>2</sup></th></tr></thead><tbody><tr><td>BCS_S<sub>k</sub><sup>1</sup></td><td>3.0<sup>a</sup></td><td>3.1<sup>b</sup></td><td>3.7<sup>c</sup></td><td>0.2</td><td>***</td></tr><tr><td>BCS_S<sub>min</sub></td><td>2.5<sup>a</sup></td><td>2.9<sup>b</sup></td><td>3.2<sup>c</sup></td><td>0.2</td><td>***</td></tr><tr><td>BCS_S<sub>210</sub></td><td>2.6<sup>a</sup></td><td>3.0<sup>b</sup></td><td>3.4<sup>c</sup></td><td>0.2</td><td>***</td></tr><tr><td>Dep_speed_S<sub>k→30</sub></td><td>-0.010<sup>a</sup></td><td>-0.003<sup>b</sup></td><td>-0.010<sup>a</sup></td><td>0.006</td><td>***</td></tr><tr><td>Rep_speed_S<sub>180→210</sub></td><td>0.0020<sup>a</sup></td><td>0.0004<sup>b</sup></td><td>0.0020<sup>a</sup></td><td>0.001</td><td>***</td></tr></tbody></table>

a-c Means with superscripts differ significantly by row.

¹ BCS_S<sub>k</sub> = sternal BCS at kidding; BCS_S<sub>min</sub> = minimum sternal BCS; BCS_S<sub>210</sub> = sternal BCS at 210 days; Dep_speed_S<sub>k→30</sub> = (BCS_S<sub>30</sub> - BCS_S<sub>k</sub>)/30; Rep_speed_S<sub>180→210</sub> = (BCS_S<sub>210</sub> - BCS_S<sub>180</sub>)/30.

² p-value resulting from Tukey's test assessing the significance of differences between profiles for each variable. NS (p<0.1), *(p<0.05); and ***(p≤0.001).

³ ST<sub>M+</sub> = Medium sternal body condition score and depletion cluster; ST<sub>M</sub> = Medium sternal body condition score and low depletion cluster; ST<sub>H+</sub> = High sternal body condition score and depletion cluster.

For lumbar BCS, Pradel's Alpine goats were more represented in the LU<sub>H+</sub> cluster, whereas Grignon's Alpine goats were more represented in the LU<sub>M</sub> cluster. Grignon's Saanen goats were more represented in the LU<sub>M</sub> cluster. Primiparous represented between 30 % and 38 % of the population in each profile for lumbar BCS. For sternal BCS, Pradel's Alpine goats were more represented in the ST<sub>M</sub> cluster, whereas Grignon's Alpine goats were more represented in the ST<sub>H+</sub> cluster. Grignon's Saanen goats were more
represented in the ST<sub>H+</sub> cluster. Primiparous represented between 30 % and 35 % of the population in each profile for sternal BCS. See Appendix B section 3 for more details.

## Diversity of phenotypic curves at lactation scale

### Associations between MY curves and BW curves

In this section, the association between MY and BW is presented. For primiparous goats, the association between MY and BW clusters is shown in Table 8a. The Chi² test was significant (P<0.001) with a Cramer's V of 0.17. The association $Y_{P-M}^m$ with $W_{P-L}^m$ accounted for the highest proportion of goats with 17.8 % of the population, followed by the associations $Y_{P-L}^m$ with $W_{P-L}^m$ and $Y_{P-M}^m$ with $W_{P-H}^m$ with 13.9 % of the population. The association $Y_{P-L}^m$ with $W_{P-H}^m$ accounted for the lowest proportion of goats with 2.8 % of the population. The remaining 51.6% of the population was almost equally distributed among the clusters.

<table><caption>**Table 8a** - Contingency table displaying the frequency of individual primiparous lactations affected to MY and BW clusters (see section 2 for clustering methodology).</caption><thead><tr><th rowspan="2">Milk yield profile</th><th colspan="2">W<sup>P-L</sup></th><th colspan="4">Body weight profile<br>W<sup>P-H</sup></th><th colspan="2">Total</th></tr><tr><th>n</th><th>%¹</th><th>n</th><th>%</th><th>n</th><th>%</th><th>n</th><th>%</th></tr></thead><tbody><tr><td>Y<sup>P-L</sup></td><td>124</td><td>13.9</td><td>45</td><td>5.0</td><td>43</td><td>4.8</td><td>212</td><td>23.7</td></tr><tr><td>Y<sup>P-L+</sup></td><td>44</td><td>4.9</td><td>25</td><td>2.8</td><td>44</td><td>4.9</td><td>113</td><td>12.7</td></tr><tr><td>Y<sup>P-M</sup></td><td>159</td><td>17.8</td><td>124</td><td>13.9</td><td>77</td><td>8.6</td><td>360</td><td>40.3</td></tr><tr><td>Y<sup>P-H</sup></td><td>60</td><td>6.7</td><td>86</td><td>9.6</td><td>62</td><td>6.9</td><td>208</td><td>23.3</td></tr><tr><td>Total</td><td>387</td><td>43.3</td><td>280</td><td>31.4</td><td>226</td><td>25.3</td><td>893</td><td>100.0</td></tr></tbody></table>

1 % = proportion of lactations among the 893 primiparous goats.

2 YpL+= Low milk yield and high persistency cluster; YpL-= Low milk yield and low persistency cluster; YpM- = Medium milk yield and low persistency cluster; YpH = High milk yield and medium persistency cluster; WpL- = Low body weight and low depletion cluster; WpH+= High body weight and high depletion cluster; WpH-= High body weight and low depletion cluster.

For multiparous goats, the association between MY and BW clusters is shown in Table 8b. The Chi² test was significant (P<0.001) with a Cramer's V of 0.17. The association $Y_{M+}^m$ with $W_{M-L}^m$ accounted for the highest proportion of goats with 18.6% of the population. The association $Y_{M+}^m$ with $W_{M-H}^m$ accounted for the lowest proportion of goats with 5.6% of the population. The remaining 75.8% of the population was almost equally distributed among the clusters.

<table><caption>**Table 8b** - Contingency table displaying the frequency of individual multiparous lactations affected to MY and BW clusters (see section 2 for clustering methodology).</caption><thead><tr><th rowspan="2">Milk yield profile</th><th colspan="2">W<sup>M-L</sup></th><th colspan="4">Body weight profile<br>W<sup>M-H</sup></th><th colspan="2">Total</th></tr><tr><th>n</th><th>%¹</th><th>n</th><th>%</th><th>n</th><th>%</th><th>n</th><th>%</th></tr></thead><tbody><tr><td>Y<sup>M+</sup></td><td>313</td><td>18.6</td><td>95</td><td>5.6</td><td>145</td><td>8.6</td><td>553</td><td>32.8</td></tr><tr><td>Y<sup>M-M</sup></td><td>242</td><td>14.4</td><td>166</td><td>9.9</td><td>140</td><td>8.3</td><td>548</td><td>32.5</td></tr><tr><td>Y<sup>M-H</sup></td><td>169</td><td>10.0</td><td>200</td><td>11.9</td><td>215</td><td>12.8</td><td>584</td><td>34.7</td></tr><tr><td>Total</td><td>724</td><td>43.0</td><td>461</td><td>27.4</td><td>500</td><td>29.7</td><td>1,685</td><td>100.0</td></tr></tbody></table>

1 % = proportion of lactations among the 1,685 multiparous goats.

2 $Y_{M+}^m$ = Medium milk yield and high persistency cluster; $Y_{M-M}^m$ = Medium milk yield and a low persistency cluster; $Y_{M-H}^m$ = High milk yield and a medium persistency cluster; $W_{M-L}^m$ = Low body weight and low depletion cluster; $W_{M-H}^m$ = High body weight and high depletion cluster; $W_{M-H}^m$ = High body weight and low depletion cluster.

The conclusions were the same for the associations between MY and lumbar BCS curves and for the associations between MY and sternal BCS curves see Appendix C section 1 and 2.
### Associations between BW and sternal BCS curves

In this section, the association between BW and sternal BCS is presented. For primiparous goats, the association between BW and sternal BCS clusters is shown in Table 9a. The Chi² test was significant (P<0.001) with a Cramer's V of 0.25. The association $W_{PL}^m$ with $ST_M$ and $W_{PH}^m$ with $ST_M$ accounted for the highest proportion of goats with 18.8% of the population, followed by the association $W_{PL}^m$ with $ST_M$ with 17.9% of the population. The association $W_{PH}^m$ with $ST_M$ accounted for the lowest proportion of goats with 1.6% of the population. The remaining 42.9% of the population was almost equally distributed among the clusters.

<table><caption>**Table 9a** - Contingency table displaying the frequency of individual primiparous lactations affected to BW and sternal BCS clusters (see section 2 for clustering methodology).</caption><thead><tr><th rowspan="2">Body weight profile</th><th colspan="4">Sternal BCS profile</th><th colspan="4">Total</th></tr><tr><th>n</th><th>ST<sub>M+</sub><sup>2</sup></th><th>%<sup>1</sup></th><th>n</th><th>ST<sub>M</sub></th><th>%</th><th>n</th><th>ST<sub>H+</sub></th><th>%</th></tr></thead><tbody><tr><td>W<sub>PL</sub><sup>2</sup></td><td>84</td><td>18.8</td><td></td><td>80</td><td>17.9</td><td></td><td>24</td><td>5.4</td><td>188</td><td>42.0</td></tr><tr><td>W<sub>PH</sub><sup>+</sup></td><td>75</td><td>16.7</td><td></td><td>84</td><td>18.8</td><td></td><td>29</td><td>6.5</td><td>188</td><td>42.0</td></tr><tr><td>W<sub>PH</sub><sup>-</sup></td><td>7</td><td>1.6</td><td></td><td>29</td><td>6.5</td><td></td><td>36</td><td>8.0</td><td>72</td><td>16.0</td></tr><tr><td>Total</td><td>166</td><td>37.1</td><td></td><td>193</td><td>43.1</td><td></td><td>89</td><td>19.9</td><td>448</td><td>100.0</td></tr></tbody></table>

¹% = proportion of lactations among the 448 primiparous goats.

²$W_{PL}^m$ = Low body weight and low depletion cluster; $W_{PH}^m$ = High body weight and high depletion cluster; $W_{PH}^m$ = High body weight and low depletion cluster; $ST_M$ = Medium sternal body condition score and depletion cluster; $ST_M$ = Medium sternal body condition score and low depletion cluster; $ST_H$ = High sternal body condition score and depletion cluster.

For multiparous goats, the association between BW and sternal BCS clusters is shown in Table 9b. The Chi² test was significant (P<0.001) with a Cramer's V of 0.18. The association $W_{PL}^m$ with $ST_M$ accounted for the highest proportion of goats with 18.6% of the population, followed by the association $W_{PL}^m$ with $ST_M$ with 14.2% of the population. The association $W_{PH}^m$ with $ST_M$ accounted for the lowest proportion of goats with 2.8% of the population. The remaining 64.4% of the population was almost equally distributed among the clusters.

<table><caption>**Table 9b** - Contingency table displaying the frequency of individual multiparous lactations affected to BW and BCS sternal clusters (see section 2 for clustering methodology).</caption><thead><tr><th rowspan="2">Body weight profile</th><th colspan="4">Sternal BCS profile</th><th colspan="4">Total</th></tr><tr><th>n</th><th>ST<sub>M+</sub><sup>2</sup></th><th>%<sup>1</sup></th><th>n</th><th>ST<sub>M</sub></th><th>%</th><th>n</th><th>ST<sub>H+</sub></th><th>%</th></tr></thead><tbody><tr><td>W<sub>PL</sub><sup>2</sup></td><td>139</td><td>14.2</td><td></td><td>182</td><td>18.6</td><td></td><td>74</td><td>7.6</td><td>395</td><td>40.5</td></tr><tr><td>W<sub>PH</sub><sup>+</sup></td><td>120</td><td>12.3</td><td></td><td>132</td><td>13.5</td><td></td><td>93</td><td>9.5</td><td>345</td><td>35.3</td></tr><tr><td>W<sub>PH</sub><sup>-</sup></td><td>27</td><td>2.8</td><td></td><td>115</td><td>11.8</td><td></td><td>94</td><td>9.6</td><td>236</td><td>24.2</td></tr><tr><td>Total</td><td>286</td><td>29.3</td><td></td><td>429</td><td>44.0</td><td></td><td>261</td><td>26.7</td><td>976</td><td>100.0</td></tr></tbody></table>

¹% = proportion of lactations among the 976 multiparous goats.

²$W_{PL}^m$ = Low body weight and low depletion cluster; $W_{PH}^m$ = High body weight and high depletion cluster; $W_{PH}^m$ = High body weight and low depletion cluster; $ST_M$ = Medium sternal body condition score and depletion cluster; $ST_M$ = Medium sternal body condition score and low depletion cluster; $ST_H$ = High sternal body condition score and depletion cluster.

The conclusions were the same for the associations between BW and lumbar BCS curves see Appendix C section 3.

### Association between lumbar and sternal BCS curves

For primiparous goats, the association between lumbar and sternal BCS clusters is shown in Table 10a. The Chi² test was significant (P<0.001) with a Cramer's V of 0.27. The association $LU_{M+}$ with $ST_M$ accounted for the highest proportion of goats with 21.4% of the population, followed by the association
LU<sub>H+</sub> with ST<sub>M+</sub> with 19.5 % of the goats. The association LU<sub>M+</sub> with ST<sub>M+</sub> accounted for the lowest proportion of goats with 6.7 % of the population. The remaining 52.4 % of the population was almost equally distributed among the clusters.

<table><caption>**Table 10a** - Contingency table displaying the frequency of individual primiparous lactations affected to BCS lumbar and BCS sternal clusters (see section 2 for clustering methodology).</caption><thead><tr><th rowspan="2">Lumbar BCS profile</th><th colspan="2">ST<sub>M+</sub><sup>2</sup></th><th colspan="2">Sternal BCS profile</th><th colspan="2">ST<sub>H+</sub></th><th colspan="2">Total</th></tr><tr><th>n</th><th>%<sup>1</sup></th><th>n</th><th>%</th><th>n</th><th>%</th><th>n</th><th>%</th></tr></thead><tbody><tr><td>LU<sub>M+</sub><sup>2</sup></td><td>80</td><td>21.4</td><td>31</td><td>8.3</td><td>28</td><td>7.5</td><td>139</td><td>37.2</td></tr><tr><td>LU<sub>M</sub></td><td>25</td><td>6.7</td><td>51</td><td>13.6</td><td>30</td><td>8.0</td><td>106</td><td>28.3</td></tr><tr><td>LU<sub>H+</sub></td><td>27</td><td>7.2</td><td>73</td><td>19.5</td><td>29</td><td>7.8</td><td>129</td><td>34.5</td></tr><tr><td>Total</td><td>132</td><td>35.3</td><td>155</td><td>41.4</td><td>87</td><td>23.3</td><td>374</td><td>100.0</td></tr></tbody></table>

For multiparous goats, the association between lumbar and sternal BCS clusters is shown in Table 10b. The Chi<fnref n="2" /> test was significant (P<0.001) with a Cramer's V of 0.35. The association LU<sub>M+</sub> with ST<sub>M+</sub> accounted for the highest proportion of goats with 20.0 % of the population, followed by the association LU<sub>H+</sub> with ST<sub>M</sub> with 18.6 % of the goats. The association LU<sub>M+</sub> with ST<sub>H+</sub> accounted for the lowest proportion of goats with 4.1 % of the population. The remaining 57.3 % of the population was almost equally distributed among the clusters.

<table><caption>**Table 10b** - Contingency table displaying the frequency of individual multiparous lactations affected to lumbar and sternal BCS clusters (see section 2 for clustering methodology).</caption><thead><tr><th rowspan="2">Lumbar BCS profile</th><th colspan="2">ST<sub>M+</sub><sup>2</sup></th><th colspan="2">Sternal BCS profile</th><th colspan="2">ST<sub>H+</sub></th><th colspan="2">Total</th></tr><tr><th>n</th><th>%<sup>1</sup></th><th>n</th><th>%</th><th>n</th><th>%</th><th>n</th><th>%</th></tr></thead><tbody><tr><td>LU<sub>M+</sub><sup>2</sup></td><td>148</td><td>20.0</td><td>52</td><td>7.0</td><td>30</td><td>4.1</td><td>230</td><td>31.1</td></tr><tr><td>LU<sub>M</sub></td><td>59</td><td>8.0</td><td>108</td><td>14.6</td><td>45</td><td>6.1</td><td>212</td><td>28.6</td></tr><tr><td>LU<sub>H+</sub></td><td>36</td><td>4.9</td><td>138</td><td>18.6</td><td>124</td><td>16.8</td><td>298</td><td>40.3</td></tr><tr><td>Total</td><td>243</td><td>32.8</td><td>298</td><td>40.3</td><td>199</td><td>26.9</td><td>740</td><td>100.0</td></tr></tbody></table>

## Diversity of phenotypic lactation curves at lifetime scale

### *MY lactation curves throughout parities*

Individual lactation transition in MY curves between successive lactations is shown in Figure 7. Between parity 1 to 4, the Chi<fnref n="2" /> test was significant (P<0.001) with a Cramer's V ranging from 0.27 to 0.32. For primiparous goats, almost half of the goats in the Y<sup>P</sup><sub>H</sub> cluster remained the most productive ones in parity 2 (Y<sup>m</sup><sub>H</sub>), while the other half switched to other clusters. More than half of the goats in the two lowest productive clusters (Y<sup>P</sup><sub>L</sub> and Y<sup>P</sup><sub>L+</sub>) switched to the Y<sup>m</sup><sub>M+</sub> cluster. Goats in the Y<sup>P</sup><sub>M+</sub> cluster were almost equally distributed among the clusters in parity 2. For multiparous goats, more than two third of the goats in the Y<sup>m</sup><sub>H</sub> cluster remained in this cluster in successive lactations. The proportion of goats that

<footnote n="1">¹% = proportion of lactations among the 374 primiparous goats.</footnote>

<footnote n="2">² LU<sub>M+</sub> = Medium lumbar body condition score and depletion cluster; LU<sub>M</sub> = Medium lumbar body condition score and low depletion cluster; LU<sub>H+</sub> =High lumbar body condition score and depletion cluster; ST<sub>M+</sub> =Medium sternal body condition score and depletion cluster; ST<sub>M</sub> = Medium sternal body condition score and low depletion cluster; ST<sub>H+</sub> =High sternal body condition score and depletion cluster.</footnote>

<footnote n="1">¹% = proportion of lactations among the 740 multiparous goats.</footnote>

<footnote n="2">² LU<sub>M+</sub> = Medium lumbar body condition score and depletion cluster; LU<sub>M</sub> = Medium lumbar body condition score and low depletion cluster; LU<sub>H+</sub> =High lumbar body condition score and depletion cluster; ST<sub>M+</sub> =Medium sternal body condition score and depletion cluster; ST<sub>M</sub> = Medium sternal body condition score and low depletion cluster; ST<sub>H+</sub> =High sternal body condition score and depletion cluster.</footnote>
remained in the $Y^{m}_{M}$-cluster in successive lactations increased with parity. Goats in the $Y^{m}_{M+}$ cluster were almost equally distributed among the clusters in successive lactations.

### BW lactation curves throughout parities

Individual lactation transition in BW curves between successive lactations is shown in Figure 8. Between parity 1 and 4, the $Chi^2$ test was significant ($P<0.001$) with a Cramer's V ranging from 0.41 to 0.44. For primiparous, goats in the $W^{P}_{H-}$ cluster switched clusters in parity 2. More than 80% of the goats in the $W^{P}_{H+}$ and in the $W^{P}_{L-}$ clusters switched to the $W^{m}_{L-}$ cluster in parity 2. For multiparous, more than two third of the goats in the $W^{m}_{H-}$ cluster remained in this cluster in successive lactations. Half of the goats in the $W^{m}_{H+}$ cluster remained in this cluster, while the other half switched clusters in successive lactations. Half of the goats in the $W^{m}_{L-}$ cluster remained in this cluster, while the other half switched clusters in successive lactations.

### BCS lactation curves throughout parities

Only sternal BCS is presented here. Individual lactation transition in sternal BCS curves between successive lactations is shown in Figure 9. Between parity 1 and 4, $Chi^2$ test was significant ($P<0.001$) with a Cramer's V ranging from 0.35 to 0.49. For primiparous, more than half of the goats in the three clusters remained in their cluster in parity 2, while the other part switched to other clusters. For multiparous, more than three quarters of the goats in the $ST_{H+}$ profile remained in this cluster in successive lactations, while the other part switched to other clusters. More than half of the goats in the $ST_{M+}$ and $ST_{M-}$ profile remained in their cluster in successive lactations, while the other part switched to other clusters.

## Discussion

Our dataset is relatively large and the frequency of measurement of the different variables is high. However, it reflected only the management of two farms. The observations made here are a starting point for a better understanding of the relationships between milk production, body condition score and body weight in goats, but will need to be confirmed in various systems. In addition, it will be necessary to add reproductive performance, which is also considered when making decisions about culling.

The first objective of this work was to characterize the diversity of phenotypic curves of performance (MY, BW, BCS) at the lactation scale.

### MY curves

For MY curves we found four clusters for primiparous goats and three clusters for multiparous goats. Parity had a strong effect on the scale of the lactation curve. Over the lactation, primiparous had lower total milk yield than multiparous goats (Gipson & Grossman, 1990). Parity also affected the shape of the lactation curve. For all parities, some clusters presented the same shape characterized by a low persistency with different milk production levels ($Y^{P}_{L-}, Y^{P}_{M-}, Y^{m}_{M-}$). These clusters, in terms of shape, were close to the mean curve of cluster 2, which represented the most common shape of lactation observed by Arnal et al. (2018) over the French dairy goat population. This cluster 2 represented 39% of the French dairy goat population, characterized by a marked peak and a medium persistency, i.e. a low persistency for our study because Arnal et al. had an additional atypical cluster with a very low persistency. For primiparous goats, one cluster combined a low level of milk with the highest persistency over the whole population ($Y^{P}_{L+}$). This is consistent with observations made by Gipson & Grossman, (1990), where persistency was the highest in primiparous goats and decreased when parity increased.
<figure><img src="image_9.png" /></figure>

Parity_1

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

Parity_2

Parity_3

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

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

Parity_3

Parity_4

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

**Figure 7** - Barplots displaying the frequency of goats affected to a MY cluster between (a) parity 1 and 2, (b) parity 2 and 3, (c) parity 3 and 4 ($Y_{P_L}^m$= Low milk yield and high persistency cluster for primiparous; $Y_{P_L}^p$= Low milk yield and low persistency cluster for primiparous; $Y_{P_M}^m$= Medium milk yield and low persistency cluster for primiparous; $Y_{P_M}^p$= High milk yield and medium persistency cluster for primiparous; $Y_{P_M}^m$= Medium milk yield and high persistency cluster for multiparous; $Y_{P_M}^p$= Medium milk yield and low persistency cluster for multiparous; $Y_{P_H}^m$= High milk yield and medium persistency cluster for multiparous).
<figure><img src="image_14.png" /></figure>

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

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

**Figure 8** - Barplots displaying the frequency of goats affected to a BW cluster between (a) parity 1 and 2, (b) parity 2 and 3, (c) parity 3 and 4 ($W_{P-L.}^P$ = Low body weight and low depletion cluster; $W_{P-H+}^P$= High body weight and high depletion cluster; $W_{H-L.}^m$ = Low body weight and low depletion cluster; $W_{H+}^m$ = High body weight and high depletion cluster; $W_{H-}^m$= High body weight and low depletion cluster).
<figure><img src="image_17.png" /><figcaption>**Figure 9** - Barplots displaying the frequency of goats affected to a sternal BCS cluster between (a) parity 1 and 2, (b) parity 2 and 3 and (c) parity 3 and 4 (ST<sub>M+</sub> = Medium sternal body condition score and depletion profile; ST<sub>M</sub> = Medium sternal body condition score and low depletion profile; ST<sub>H+</sub> = High sternal body condition score and depletion profile).</figcaption></figure>
This can be explained by a lower level of development of the mammary gland (Safayi et al., 2010). This shape of lactation curve was also observed in the study of Arnal et al. (2018). However, persistency and MY are not always negatively correlated, because for all parities we observed that the highest productive clusters were those with a medium persistency (Y<sup>p</sup><sub>H</sub>, Y<sup>m</sup><sub>H</sub>) rather than those with the lowest persistency. This result is close to the finding of Arnal et al. (2018) that their highest productive cluster had a high persistency. We can hypothesize that better fed goats produce more milk and are better able to maintain that production.

Despite differences between Saanen and Alpine goats being reported in the literature (Gipson & Grossman, 1990; Rupp et al., 2011; Arnal et al., 2018), breed did not have a significant effect on the scale or the shape of the lactation curves in the present study, regardless of parity on the two farms. It should be noted that the two breeds were only present on one farm (Grignon) and thus had the same feeding and management environment.

## Lactation curves of BW

For BW curves we found three clusters for primiparous and multiparous goats. Parity and breed played a strong role on the scale of BW curves. As expected, primiparous goats were lighter than multiparous ones. For all parities, we found low depletion clusters (W<sup>p</sup><sub>L</sub>, W<sup>m</sup><sub>L</sub>, W<sup>p</sup><sub>H</sub>, W<sup>m</sup><sub>H</sub>) and high depletion clusters (W<sup>p</sup><sub>H+</sub>, W<sup>m</sup><sub>H+</sub>). The low depletion clusters had the same shape but differed in terms of level. Only the high depletion clusters differed in terms of shape from the other clusters. However, the depletion speed was lower in primiparous goats than in multiparous goats. Indeed, for primiparous goats, the difference between kidding and the minimum of BW averaged 3.7 kg, while for multiparous this difference averaged 8.3 kg. These results are consistent with what Sauvant et al. (2012) observed when they modelled the BW curve by parity. They observed that primiparous were lighter and lost less BW (4.0 kg on average) than multiparous goats (7.3 kg on average). To our knowledge, little work has been done to characterize BW curves in dairy goats. Our work can be compared to the study of Macé et al. (2019) in meat sheep. They analyzed BW longitudinal data in 1146 ewes to characterize curves over multiple production cycles. Most of their curves had the same shape but differed in terms of level.

All multiparous clusters had a higher BW at kidding (BW<sub>k</sub>) than at the beginning of the subsequent gestation (BW<sub>210</sub>). In contrast, for primiparous clusters the opposite was mainly true, indicating that primiparous goats were still growing in first lactation. For multiparous and for all clusters BW<sub>k</sub> was higher than BW<sub>210</sub> at the end of lactation. BW is easy to measure on farm to monitor animals, especially to quantify energy balance (Thorup et al., 2012). However, BW measures also include digestive content, growth, gravid uterus and body reserves. Therefore, BW measures alone are not consistent enough to quantify body reserve changes. They need to be analyzed with BCS to better understand body reserves dynamics.

A breed effect was observed for BW curves: Saanen goats were more represented in the high-level clusters for all parities (W<sup>p</sup><sub>H-</sub>, W<sup>m</sup><sub>H-</sub>). They were generally heavier than Alpine goats (Sauvant et al., 2012).

## Lactation curves of BCS

For lumbar and sternal BCS we found three clusters for all parities. First for all parities, we found high depletion clusters for lumbar (LU<sub>M+</sub>, LU<sub>H+</sub>) and sternal (ST<sub>M+</sub>, ST<sub>H+</sub>) BCS. Then, we found low depletion clusters for lumbar (LU<sub>M</sub>) and sternal (ST<sub>M</sub>) BCS. High depletion clusters presented the same shape but differed in terms of level. Only the low depletion clusters differed in terms of shape. These results are also consistent with the observations of Macé et al. (2019). They found the same shape but differing in terms of level. Moreover, for the high depletion clusters, the variation between kidding and the minimum of BCS averaged 0.4 points for LU<sub>M+</sub>, 0.3 points for LU<sub>H+</sub>, and 0.5 points for ST<sub>M+</sub> and ST<sub>H+</sub>. Our values, especially for sternal BCS, are lower but close to those described in the French feeding system (Inra, 2018).
Parity did not significantly affect BCS curves. Indeed, primiparous goats represented a third of the whole population in each of the clusters. Breed did not significantly affect BCS curves. We observed only a farm effect on BCS because Grignon's Alpine and Saanen goats were more represented in the $LU_M$ and $ST_{H+}$ clusters. This can probably be explained by differences in the personnel carrying out the BCS evaluation (although differences in herd management, or a random distribution linked to the clustering approach cannot be excluded).

## A great diversity of associations among biological functions

The second objective of this work was to assess the diversity of associations among the different phenotypic curves. We investigated whether one phenotypic curve was associated with another. At the lactation scale, the $Chi^2$ test was significant for associations but the Cramer's V showed weak to moderate values (globally less than 0.4) (Kotrlík et al., 2011). This lack of strong associations among lactation curves of MY, BW and BCS suggests there exists a relatively large diversity of energy partitioning strategies among individuals. Associations among MY, BW and BCS were well-studied in dairy cows. Some studies showed a positive correlation between pre-calving BCS and milk production (Waltner et al., 1993; Roche et al., 2007), whereas other studies did not find any relationship between these variables (Garnsworthy & Topps, 1982; Garnsworthy & Jones, 1987). More recently, Ollion et al. (2016) assessed the diversity of trade-offs between milk production, body reserves and reproduction in early lactation dairy cows. They showed four different trade-off profiles according to a priority given to a biological function. The first trade-off profile represented cows giving priority to lactation instead of reproduction, the second trade-off profile represented cows giving priority to reproduction instead of lactation. The third trade-off profile represented cows with poor performances in all functions, and the last trade-off profile represented cows with no trade-off among functions. All of these approaches considered correlations between traits at one time point and not at the whole lactation scale. Moreover, these performance traits were evaluated at the beginning of the lactation where cows exhibited a negative energy balance allowing energy partitioning in favor of milk over body reserves. Another possible explanation for the lack of strong associations found in our study is that trade-off between life functions, and therefore correlations between traits, are well expressed when animals face feed shortage (Blanc et al., 2006; Friggens et al., 2017). Our data came from two experimental farms where we can assume that animals are well managed and not so constrained in terms of nutrition.

The diversity of associations among biological functions found in the present study could also suggest a great diversity in intake at the individual level. However, we are not able today to accurately evaluate individual intake. Furthermore, BCS was used as a proxy to evaluate body reserve dynamics. As a subjective evaluation of body reserves, BCS was probably not accurate enough to capture a relationship with MY. MY fat would have been important to consider, because at equal MY fat could vary a lot. However, this variable was not considered in our study because although some data were available, the frequency was low (less than one measurement per month). On the one hand, this diversity of biological profiles can be seen as a potential resource to improve farming system resilience (Dumont et al., 2020). On the other hand, this diversity raises questions about feeding systems that assumed a relationship between a BW and a MY curve. There is a need to better quantify body reserves contribution in terms of energy to goat's requirements (Inra, 2018). These findings question management strategies that are based on the average animal, i.e; ignoring cluster types. A perspective can be to adapt management strategies to the diversity of individual profiles in terms of phenotypic curves and then better match animal's requirements.

The final objective of this work was to assess the diversity of phenotypic curves at the lifetime scale. For each phenotypic trait, the $Chi^2$ test was significant. Cramer's V test showed lower values for MY than for BW and BCS suggesting stronger associations for BW and BCS. For MY curves, we saw for primiparous goats that almost half of the goats in the $Y_{P_H}$ cluster remained in this cluster in parity 2, while the lowest
productive goats ($Y_{P-L}$ and $Y_{P+L}$) switched cluster in parity 2. For multiparous, we observed a more stable pattern of cluster membership with two thirds of the goats in the $Y_{P+L}$ cluster remaining in this cluster in successive lactations. Usually, milk production increases from first to fourth parity. After the fourth parity, the level of milk production decreases (Arnal et al., 2018). However, with genetic improvement, we can make the hypothesis that some goats can reach their milk potential earlier. Goats that stayed in the highest productive clusters could be animals that have reached their milk potential relatively early. Goats that are changing clusters could be the ones that have not reached their potential early. For BW curves across parities, we saw that for primiparous goats, most of the goats in the $W_{P-L}$ remained in the lowest BW cluster ($W_{P+L}$) in parity 2, while $W_{P+L}$ switched to the $W_{P-L}$ cluster. Goats in the $W_{P+L}$ presented the highest depletion speed, so they were not able to recover from the intense depletion and remained in the lowest cluster in parity 2. For multiparous goats, we also observed a pattern of cluster membership, with more than three-quarters of the goats in the $W_{P+L}$ profile remaining in this cluster in successive lactations. Half of the goats in the $W_{P-L}$ remained in this cluster in successive lactations. For sternal BCS curves across parities, we saw that more than half of the primiparous goats in the three clusters remained in their cluster in parity 2. For multiparous goats, we observed that three-quarters of the goats in the $ST_{P+L}$ cluster remained in this cluster in successive lactations. More than half of the goats in the $ST_{P+L}$ and $ST_{P-L}$ remained in their cluster in successive lactations. These observations on BW and BCS over successive lactations, are consistent with what Macé et al. (2019) observed in meat sheep. They observed one-third up to half of ewes remaining in the same profile during successive cycles of production. They supposed that changes in profile distribution could be linked to litter size that can play a role in body weight depletion. We did not consider the prolificacy of the goats in our study. This information was missing for 35% of the animals. When the number of kids was known (single kid for 33% of the goats, two kids for 51%, three kids and more for 16%), we did not find any relationship with our clusters. However, it is an information to consider in further analysis as it is described to be a factor related to milk production (Hayden et al., 1979; Zamuner et al., 2020). These results highlighted the importance of a lifetime approach to better understand potential changes in priorities among functions and see how an early lifetime performance can impact the whole productive lifespan (Puillet & Martin, 2017). Lifetime and longevity approaches are increasingly being studied because in France from 1991 to 2011, the female productive life decreased by 346 days, which led to an average productive lifespan of 2.7 years per goat (Palhière et al., 2018), which increases replacement costs.

## A methodology to analyze trade-off between phenotypic curves with heterogeneous data frequency

This methodology was built to analyze the trade-off between phenotypic lactation curves based on longitudinal data with different frequencies. We used models adapted to the data frequency to better characterize our curves. However, this approach implied the creation of synthetic indicators to have the same baseline for phenotypic curves characterized by different models. For MY curves, synthetic indicators were simple to find, because we used common indicators to summarize a lactation curve with level and dynamic indicators such as the $MY_{peak}$, Peak time and Persistency. However, because the BW (dataset 1) and BCS data were less frequent (datasets 1 and 2) less elaborate models were used. This then meant that a more simple set of summary indicators was used to characterize these curves, which may not be as informative as those for MY. With heterogeneity in frequencies, it is difficult to use the same models to capture phenotypic curves. Differences in frequencies could lead to use simple models with parameters that are not always biologically meaningful. Or it may lead to the use of more complex models that deal with problems of parameters identifiability. It is important to find a way to use biologically meaningful parameters from different models as inputs for a clustering approach. This approach with model parameters will help to summarize the phenotypic curves without considering synthetic indicators.
## Further development and potential use of on-farm record for managing animal

With development of on-farm automatic measuring technologies, more frequent data for MY or BW are becoming available. Some authors developed methods to characterize new indicators such as the deviation of milk production from a theoretical potential production (Ben Abdelkrim et al., 2020; Poppe et al., 2020; Adriaens et al., 2021). Intense and rapid MY or BW losses might be used as indicators of disease or metabolic disorders. Being able to identify these animals is of great interest for farming management. In our study, we used specific models that dissociated the effects of perturbation from a theoretical unperturbed curve. To characterize phenotypic curves, we focused only on unperturbed curves, which represented the potential production that an animal could have in a non-perturbed environment. With unperturbed MY and BW curves we saw a diversity of associations, it would be of future interest to also consider perturbations. The extent to which there are common perturbations on MY and BW curves may be informative. This approach has been used in dairy cows where Ben Abdelkrim et al. (2021) identified common perturbations in MY and BW. Using perturbations in a phenotypic curve analysis could help to select animals that better cope with their environment.

Data acquisition for BCS is more complicated in goats than in dairy cows. Manual BCS evaluation provided satisfactory results but is still a subjective method depending on the operator (Lerch et al., 2021). Recent studies have shown that new methods such as 3-dimension imaging did not provide satisfactory estimators of body composition and further developments may be needed to develop a robust phenotyping tool (Lerch et al., 2021). For all parities, BCS curves were well discriminated one month after kidding and stayed constant over the whole lactation. This observation suggests that BCS measures frequency can be reduced to key periods (kidding period, two months before breeding period, dry-off). This paper is the first step of a study that will include reproductive success in the analysis. Including reproduction outcome will help to predict fertility according to phenotypic curves for a given lactation. This analysis will be conducted also on the lifetime scale to look for potential unfavorable clusters. These further analyses will clarify this diversity of phenotypic curves and will provide metrics to better manage at-risk animals in terms of reproduction (e.g., finding the best periods to monitor at-risk animals). In the dairy goat sector, extended lactations became an alternative farming management to reduce culling and give another chance for a goat to reproduce. Being able to make early decisions on reproductive management, can be of economic interest and may increase sustainability (Adriaens et al., 2020).

## Conclusion

With a multi-scale approach on MY, BW and BCS time-series data, it was possible to characterize the diversity of associations between phenotypic lactation curves related to milk production and the use of body reserves. For each of MY, BW and BCS, the lactation curves clustered into 4 (MY) or 3 (BW, BCS) clusters. The diversity of associations at the lactation scale between clusters suggests a diversity of energy partitioning strategies among goats, which may provide different adaptive responses to environmental perturbations. Our results challenge mainstream management strategies that are based on average animal profiles. Rather, considering diversity of performance profiles can be a way to better adapt management to individuals or groups of individuals to improve their robustness. At the lifetime scale, change among clusters are more pronounced between first and second lactation, while a stable pattern of cluster membership appears for multiparous goats. Indeed, more than two thirds of the highest clusters for each phenotypic curve remained in these clusters in successive lactations. To further identify some clusters or combination of clusters that are at risk of culling, a first perspective of this study is to combine reproductive performance with MY, BCS and BW curves and then provide metrics to better manage animals at risk of culling.
## Acknowledgements

We gratefully acknowledge the teams at the INRAE UMR 791 Modélisation Systémique Appliquée aux Ruminants (Ile-De-France, France) experimental farm and Le Pradel experimental farm (Ardeche, France) for the care of the animals and their work to provide the different dataset. We gratefully acknowledge Maxime Legris from l'Institut de l'Elevage (Ile-De-France, France) for his help on statistical support for this paper. We gratefully acknowledge Nicolas Friggens for his help to revise the English of the manuscript.

Preprint version 3 of this article has been peer-reviewed and recommended by Peer Community In Animal Science (https://doi.org/10.24072/pci.animsci.100234; Atzori, 2024).

## Funding

This work was funded by the French Livestock Institute and was supported by ANRT (Association Nationale de la Recherche et de la Technologie) with a CIFRE doctoral fellowship (2020/1547). This doctoral fellowship is a part of a research axis of a Technical Joint Unit called "Sustainable goat systems of tomorrow" (UMT SC3D).

## Conflict of interest disclosure

The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article. Laurence Puillet is a recommender of PCI Animal Science.

## Data, scripts, code, and supplementary information availability

Data and scripts/code are deposited on Recherche Data Gouv (https://doi.org/10.57745/C1XPQ2; Gafsi, 2023a)

## Appendices

Appendices are deposited on Zenodo (https://doi.org/10.5281/zenodo.10090604; Gafsi, 2023b)

## References

Adriaens I, Friggens NC, Ouweltjes W, Scott H, Aernouts B, Statham J. 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 2020;103:7155-71. https://doi.org/10.3168/jds.2019-17826.

Adriaens I, van den Brulle I, D'Anvers L, Statham JME, Geerinckx K, De Vliegher S, et al. Milk losses and dynamics during perturbations in dairy cows differ with parity and lactation stage. J Dairy Sci 2021;104:405-18. https://doi.org/10.3168/jds.2020-19195.

Arnal M, Robert-Granie C, Larroque H. Diversity of dairy goat lactation curves in France. J Dairy Sci 2018;101:11040-51. https://doi.org/10.3168/JDS.2018-14980.

Atzori, A. (2024) Understanding milk and body reserves trajectories and nutrient partitioning in dairy goats through a modelling approach. Peer Community in Animal Science, 100234. https://doi.org/10.24072/pci.animsci.100234

Bauman DE, Bruce Currie W. Partitioning of Nutrients During Pregnancy and Lactation: A Review of Mechanisms Involving Homeostasis and Homeorhesis. J Dairy Sci 1980;63:1514-29. https://doi.org/10.3168/jds.S0022-0302(80)83111-0.
Ben Abdelkrim A, Puillet L, Gomes P, Martin O. Lactation curve model with explicit representation of perturbations as a phenotyping tool for dairy livestock precision farming. Animal 2020;15:100074. https://doi.org/10.1016/j.animal.2020.100074.

Ben Abdelkrim A, Tribout T, Martin O, Boichard D, Ducrocq V, Friggens NC. Exploring simultaneous perturbation profiles in milk yield and body weight reveals a diversity of animal responses and new opportunities to identify resilience proxies. J Dairy Sci 2021;104:459-70. https://doi.org/10.3168/jds.2020-18537.

Berry DP, Buckley F, Dillon P, Evans RD, Rath M, Veerkamp RF. Estimation of genotype×environment interactions, in a grass-based system, for milk yield, body condition score, and body weight using random regression models. Livest Prod Sci 2003;83:191-203. https://doi.org/10.1016/S0301-6226(03)00098-8.

Blanc F, Bocquier F, Agabriel J, D'Hour P, Chilliard Y. Adaptive abilities of the females and sustainability of ruminant livestock systems. A review. Anim Res 2006;55:489-510. https://doi.org/10.1051/animres:2006040.

Cobby JM, Le Du YLP. On fitting curves to lactation data. Anim Prod 1978;26:127-33. https://doi.org/10.1017/S0003356100039532.

Dhanoa MS. A note on an alternative form of the lactation model of Wood. Anim Prod 1981;32:349-51. https://doi.org/10.1017/S0003356100027276.

Dijkstra J, France J, Dhanoa MS, Maas JA, Hanigan MD, Rook AJ, et al. A Model to Describe Growth Patterns of the Mammary Gland during Pregnancy and Lactation. J Dairy Sci 1997;80:2340-54. https://doi.org/10.3168/jds.S0022-0302(97)76185-X.

Dumont B, Puillet L, Martin G, Savietto D, Aubin J, Ingrand S, et al. Incorporating Diversity Into Animal Production Systems Can Increase Their Performance and Strengthen Their Resilience. Front Sustain Food Syst 2020;4:1-15. https://doi.org/10.3389/fsufs.2020.00109.

Friggens NC, Blanc F, Berry DP, Puillet L. Review: Deciphering animal robustness. A synthesis to facilitate its use in livestock breeding and management. Animal 2017;11:2237-51. https://doi.org/10.1017/S175173111700088X.

Friggens NC, Disenhaus C, Petit H V. Nutritional sub-fertility in the dairy cow: towards improved reproductive management through a better biological understanding. Animal 2010;4:1197-213. https://doi.org/10.1017/S1751731109991601.

Friggens NC, Emmans GC, Veerkamp RF. On the use of simple ratios between lactation curve coefficients to describe parity effects on milk production a. vol. 62. 1999, https://doi.org/10.1016/S0301-6226(99)00110-4.

Friggens NC. Body lipid reserves and the reproductive cycle: towards a better understanding. Livest Prod Sci 2003;83:219-36. https://doi.org/10.1016/S0301-6226(03)00111-8.

Gafsi N, Martin O, Bidan F, Grimard B, Puillet L. Diversity of performance patterns in dairy goats: multi-scale analysis of milk yield, body condition score and body weight trajectories_Data_and_R_code. 2023a. Recherche Data Gouv, V1 https://doi.org/10.57745/C1XPQ2.

Gafsi N, Martin O, Bidan F, Grimard B, Puillet L. Diversity of performance patterns in dairy goats: multi-scale analysis of milk yield, body condition score and body weight trajectories. 2023b. Zenodo. https://doi.org/10.5281/zenodo.10090605

Garnsworthy PC, Jones GP. The influence of body condition at calving and dietary protein supply on voluntary food intake and performance in dairy cows. Anim Prod 1987;44:347-53. https://doi.org/10.1017/S0003356100012277.

Garnsworthy PC, Topps JH. The effect of body condition of dairy cows at calving on their food intake and performance when given complete diets. Anim Prod 1982;35:113-9. https://doi.org/10.1017/S0003356100000878.

Gipson TA, Grossman M. Lactation curves in dairy goats: a review. Small Rumin Res 1990;3:383-96. https://doi.org/10.1016/0921-4488(90)90019-3.
Grossman M, Hartz SM, Koops WJ. Persistency of lactation yield: A novel approach. J Dairy Sci 1999;82:2192-7. https://doi.org/10.3168/jds.S0022-0302(99)75464-0.

Hayden TJ, Thomas CR, Forsyth IA. Effect of Number of Young Born (Litter Size) on Milk Yield of Goats: Role for Placental Lactogen. J Dairy Sci 1979;62:53-7. https://doi.org/10.3168/jds.S0022-0302(79)83201-4.

Inra. INRA feeding system for ruminants. 2018:728 p. https://doi.org/10.3920/978-90-8686-292-4.

Kotrlik JW, Atherton JC, Williams HA, Khata Jabor M. Reporting and Interpreting Effect Size in Quantitative Agricultural Education Research. J Agric Educ 2011;52:132-42. https://doi.org/10.5032/jae.2011.01132.

Lerch S, De La Torre A, Huau C, Monziols M, Xavier C, Louis L, et al. Estimation of dairy goat body composition: A direct calibration and comparison of eight methods. Methods 2021;186:68-78. https://doi.org/10.1016/j.ymeth.2020.06.014.

Macé T, González-García E, Carrière F, Douls S, Foulquié D, Robert-Granié C, et al. Intra-flock variability in the body reserve dynamics of meat sheep by analyzing BW and body condition score variations over multiple production cycles. Animal 2019;13:1986-98. https://doi.org/10.1017/S175173111800352X.

Macé T, González-García E, Kövér G, Hazard D, Taghipoor M. PhenoBR: a model to phenotype body condition dynamics in meat sheep. Animal 2023;17:100845. https://doi.org/10.1016/j.animal.2023.100845.

Martin O, Ben Abdelkrim A. Detecting pertubations in dairy cows liveweight trajectories, 70th EAAP Annual Meeting, Ghent, Belgium (hal-03190731) 2019. https://hal.inrae.fr/hal-03190731

Morand-Fehr P, Hervieu J. Apprécier l'état corporel des chèvres. Intérêt et méthode. La Chèvre 1999;231:22-5.

Ollion E, Ingrand S, Delaby L, Trommenschlager JM, Colette-Leurent S, Blanc F. Assessing the diversity of trade-offs between life functions in early lactation dairy cows. Livest Sci 2016;183:98-107. https://doi.org/10.1016/j.LIVSCI.2015.11.016.

Palhière I, Oget C, Rupp R. Functional longevity is heritable and controlled by a major gene in French dairy goats. Proc World Congr Genet Appl to Livest Prod 11-16 Febr 2018 2018:165.

Pollott GE. A biological approach to lactation curve analysis for milk yield. J Dairy Sci 2000;83:2448-58. https://doi.org/10.3168/jds.S0022-0302(00)75136-8.

Poppe M, Veerkamp RF, van Pelt ML, Mulder HA. Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding. J Dairy Sci 2020;103:1667-84. https://doi.org/10.3168/jds.2019-17290.

Pryce JE, Coffey MP, Simm G. The relationship between body condition score and reproductive performance. J Dairy Sci 2001;84:1508-15. https://doi.org/10.3168/jds.S0022-0302(01)70184-1.

Puillet L, Martin O. A dynamic model as a tool to describe the variability of lifetime body weight trajectories in livestock females. J Anim Sci 2017;95:4846-56. https://doi.org/10.2527/jas2017.1803.

Roche JR, Friggens NC, Kay JK, Fisher MW, Stafford KJ, Berry DP. Invited review: Body condition score and its association with dairy cow productivity, health, and welfare. J Dairy Sci 2009;92:5769-801. https://doi.org/10.3168/jds.2009-2431.

Roche JR, Macdonald KA, Burke CR, Lee JM, Berry DP. Associations Among Body Condition Score, Body Weight, and Reproductive Performance in Seasonal-Calving Dairy Cattle. J Dairy Sci 2007;90:376-91. https://doi.org/10.3168/jds.S0022-0302(07)72639-5.

Rupp R, Clément V, Piacere A, Robert-Granié C, Manfredi E. Genetic parameters for milk somatic cell score and relationship with production and udder type traits in dairy Alpine and Saanen primiparous goats. J Dairy Sci 2011;94:3629-34. https://doi.org/10.3168/jds.2010-3694.

Safayi S, Theil PK, Elbrønd VS, Hou L, Engbæk M, Nørgaard J V., et al. Mammary remodeling in primiparous and multiparous dairy goats during lactation. J Dairy Sci 2010;93:1478-90. https://doi.org/10.3168/jds.2009-2422.
Sauvant D, Giger-Reverdin S, Meschy F, Puillet L, Schmidely P. Actualisation des recommandations alimentaires pour les chèvres laitières. Prod Anim 2012;25:259-76. https://doi.org/10.20870/productions-animales.2012.25.3.3214.

Thorup VM, Edwards D, Friggens NNC. On-farm estimation of energy balance in dairy cows using only frequent body weight measurements and body condition score. J Dairy Sci 2012;95:1784-93. https://doi.org/10.3168/jds.2011-4631.

Tukey JW. Exploratory Data Analysis (Book Section). Explor Data Anal 1977:61-100.

Waltner SS, McNamara JP, Hillers JK. Relationships of Body Condition Score to Production Variables in High Producing Holstein Dairy Cattle. J Dairy Sci 1993;76:3410-9. https://doi.org/10.3168/jds.S0022-0302(93)77679-1.

Wilmink JBM. Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation. Livest Prod Sci 1987;16:335-48. https://doi.org/10.1016/0301-6226(87)90003-0.

Wood PDP. Algebraic model of the lactation curve in cattle. Nature 1967;216:164-5. https://doi.org/10.1038/216164a0.

Zamuner F, DiGiacomo K, Cameron AWN, Leury BJ. Effects of month of kidding, parity number, and litter size on milk yield of commercial dairy goats in Australia. J Dairy Sci 2020;103:954-64. https://doi.org/10.3168/jds.2019-17051.