The dairy goat sector faces many challenges, such as animals with reduced longevity (Palhière et al., 2018 I. Palhière; C. Oget; R. Rupp Functional longevity is heritable and controlled by a major gene in French dairy goats, Proceedings of the World Congress on Genetics Applied to Livestock Production, 11-16 February 2018 (2018), p. 165 | DOI) 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 D. E. Bauman; B. W. Currie Partitioning of Nutrients During Pregnancy and Lactation: A Review of Mechanisms Involving Homeostasis and Homeorhesis, Journal of Dairy Science, Volume 63 (1980) no. 9, pp. 1514-1529 | DOI; Friggens et al., 2017 N. Friggens; F. Blanc; D. Berry; L. Puillet Review: Deciphering animal robustness. A synthesis to facilitate its use in livestock breeding and management, Animal, Volume 11 (2017) no. 12, pp. 2237-2251 | DOI). 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 J. E. Pryce; M. P. Coffey; G. Simm The relationship between body condition score and reproductive performance, Journal Of Dairy Science, Volume 84 (2001) no. 6, pp. 1508-1515 | DOI; Roche et al., 2009 J. R. Roche; N. Friggens; J. K. Kay; M. W. Fisher; K. J. Stafford; D. P. Berry Invited review: Body condition score and its association with dairy cow productivity, health, and welfare., Journal of dairy science, Volume 92 (2009) no. 12, pp. 5769-5801 | DOI; Friggens et al., 2010 N. C. Friggens; C. Disenhaus; H. V. Petit Nutritional sub-fertility in the dairy cow: towards improved reproductive management through a better biological understanding, Animal : an international journal of animal bioscience, Volume 4 (2010) no. 7, pp. 1197-1213 | DOI). 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 N. C. Friggens Body lipid reserves and the reproductive cycle: towards a better understanding, Livestock Production Science, Volume 83 (2003) no. 52nd Annual EAAP Meeting, pp. 219-236 | DOI).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 P. D. P. Wood Algebraic Model of the Lactation Curve in Cattle, Nature, Volume 216 (1967) no. 5111, pp. 164-165 | DOI; Cobby & Le Du, 1978 J. M. Cobby; Y. L. Le Du On fitting curves to lactation data, Animal Production, Volume 26 (1978) no. 2, pp. 127-133 | DOI; Dhanoa, 1981 M. S. Dhanoa A note on an alternative form of the lactation model of Wood, Animal Science, Volume 32 (1981) no. 3, pp. 349-351 | DOI; Wilmink, 1987 J. Wilmink Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation, Livestock Production Science, Volume 16 (1987) no. 4, pp. 335-348 | DOI) and some wanted to have models based on a biological framework (Dijkstra et al., 1997 J. Dijkstra; J. France; M. Dhanoa; J. Maas; M. Hanigan; A. Rook; D. Beever A Model to Describe Growth Patterns of the Mammary Gland During Pregnancy and Lactation, Journal of Dairy Science, Volume 80 (1997) no. 10, pp. 2340-2354 | DOI; Friggens et al., 1999 N. Friggens; G. Emmans; R. Veerkamp On the use of simple ratios between lactation curve coefficients to describe parity effects on milk production, Livestock Production Science, Volume 62 (1999) no. 1, pp. 1-13 | DOI; Pollott, 2000 G. E. Pollott A biological approach to lactation curve analysis for milk yield, Journal of Dairy Science, Volume 83 (2000) no. 11, pp. 2448-2458 | DOI). 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 A. Ben Abdelkrim; L. Puillet; P. Gomes; O. Martin Lactation curve model with explicit representation of perturbations as a phenotyping tool for dairy livestock precision farming, Animal, Volume 15 (2020) no. 1, p. 100074 | DOI). 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 M. Arnal; C. Robert-Granie; H. Larroque Diversity of dairy goat lactation curves in France, Journal of Dairy Science, Volume 101 (2018) no. 12, pp. 11040-11051 | DOI). Studies on modelling the shape of BW or BCS curves through lactation (Macé et al., 2023 T. Macé; E. González-García; G. Kövér; D. Hazard; M. Taghipoor PhenoBR: a model to phenotype body condition dynamics in meat sheep, Animal, Volume 17 (2023) no. 6, p. 100845 | DOI) are less frequent. Some mathematical functions with an exponential approach (Sauvant et al., 2012 D. Sauvant; S. Giger-Reverdin; F. Meschy; L. Puillet; P. Schmidely Actualisation des recommandations alimentaires pour les chèvres laitières, Productions Animales, Volume 25 (2012) no. 3, pp. 259-276 | DOI) or a random regression approach (Berry et al., 2003 D. P. Berry; F. Buckley; P. Dillon; R. D. Evans; M. Rath; R. F. Veerkamp Estimation of genotype×environment interactions, in a grass-based system, for milk yield, body condition score, and body weight using random regression models, Livestock Production Science, Volume 83 (2003) no. 2-3, pp. 191-203 | DOI) 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.

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).

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 P. Morand-Fehr; J. Hervieu Apprécier l'état corporel des chèvres. Intérêt et méthode, La chèvre, Volume 231 (1999), pp. 22-25 | DOI). 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).

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).

Models were selected according to data frequency.

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

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 P. D. P. Wood Algebraic Model of the Lactation Curve in Cattle, Nature, Volume 216 (1967) no. 5111, pp. 164-165 | DOI) 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.

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.

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.

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

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 J. W. Tukey Exploratory Data Analysis (Book Section), Explor Data Anal (1977), 768, pp. 61-100 | DOI) 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.

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).

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 N. Gafsi; O. Martin; F. Bidan; B. Grimard; L. Puillet Diversity of performance patterns in dairy goats: multi-scale analysis of milk yield, body condition score and body weight trajectories, Zenodo, 2023 | DOI). 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.

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.

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.

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, 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.

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.

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

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).

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

Primiparous clusters were characterized by:

Multiparous clusters were characterized by:

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^p_M- 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^m_H cluster, whereas Grignon’s Alpine goats were less represented in this cluster. Grignon’s Saanen goats were more represented in the Y^m_M+ cluster. See Appendix B section 1 for more details.

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).

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

Primiparous clusters were characterized by:

Multiparous clusters were characterized by:

For primiparous goats, Pradel’s Alpine goats were more represented in the W^p_L- and W^p_H+ clusters. Grignon’s Alpine goats were more represented in the W^p_L- cluster. Grignon’s Saanen goats were more represented in the W^p_H- cluster. For multiparous goats, Pradel’s Alpine goats and Grignon’s Alpine goats were more represented in the W^m_L- cluster. Grignon’s Saanen goats were more represented in the W^m_H- cluster. See Appendix B section 2 for more details.

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_min and BCS_L_k) 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_S_min and BCS_S₂₁₀) 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).

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

Lumbar BCS clusters were characterized by:

Sternal BCS profiles were characterized by:

For lumbar BCS, Pradel’s Alpine goats were more represented in the LU_H+ cluster, whereas Grignon’s Alpine goats were more represented in the LU_M cluster. Grignon’s Saanen goats were more represented in the LU_M 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_M cluster, whereas Grignon’s Alpine goats were more represented in the ST_H+ cluster. Grignon’s Saanen goats were more represented in the ST_H+ cluster. Primiparous represented between 30 % and 35 % of the population in each profile for sternal BCS. See Appendix B section 3 for more details.

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- with W^p_L- accounted for the highest proportion of goats with 17.8 % of the population, followed by the associations Y^p_L- with W^p_L- and Y^p_M- with W^p_H+ with 13.9 % of the population. The association Y^p_L+ with W^p_H+ 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.

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- accounted for the highest proportion of goats with 18.6% of the population. The association Y^m_M+ with W^m_H+ 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.

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.

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^p_L- with ST_M+ and W^p_H+ with ST_M accounted for the highest proportion of goats with 18.8 % of the population, followed by the association W^p_L- with ST_M with 17.9 % of the population. The association W^p_H- 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.

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^m_L- with ST_M accounted for the highest proportion of goats with 18.6 % of the population, followed by the association W^m_L- with ST_M+ with 14.2 % of the population. The association W^m_H- 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.

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

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_H+ with ST_M with 19.5 % of the goats. The association LU_M with ST_M+ 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.

For multiparous goats, the association between lumbar and sternal BCS clusters is shown in Table 10b. The Chi² test was significant (P<0.001) with a Cramer’s V of 0.35. The association LU_M+ with ST_M+ accounted for the highest proportion of goats with 20.0 % of the population, followed by the association LU_H+ with ST_M with 18.6 % of the goats. The association LU_M+ with ST_H+ 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.

Individual lactation transition in MY curves between successive lactations is shown in Figure 7. Between parity 1 to 4, the Chi² 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^p_H cluster remained the most productive ones in parity 2 (Y^m_H), while the other half switched to other clusters. More than half of the goats in the two lowest productive clusters (Y^p_L- and Y^p_L+ ) switched to the Y^m_M+ cluster. Goats in the Y^p_M+ cluster were almost equally distributed among the clusters in parity 2. For multiparous goats, more than two third of the goats in the Y^m_H cluster remained in this cluster in successive lactations. The proportion of goats that 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.

Individual lactation transition in BW curves between successive lactations is shown in Figure 8. Between parity 1 and 4, the Chi² 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.

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² 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.

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.

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 T. A. Gipson; M. Grossman Lactation curves in dairy goats: a review, Small Ruminant Research, Volume 3 (1990) no. 4, pp. 383-396 | DOI). 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.

This can be explained by a lower level of development of the mammary gland (Safayi et al., 2010 S. Safayi; P. K. Theil; V. S. Elbrønd; L. Hou; M. Engbæk; J. V. Nørgaard; K. Sejrsen; M. O. Nielsen Mammary remodeling in primiparous and multiparous dairy goats during lactation, Journal of Dairy Science, Volume 93 (2010) no. 4, pp. 1478-1490 | DOI).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^p_H, Y^m_H) 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 T. A. Gipson; M. Grossman Lactation curves in dairy goats: a review, Small Ruminant Research, Volume 3 (1990) no. 4, pp. 383-396 | DOI; Rupp et al., 2011 R. Rupp; V. Clément; A. Piacere; C. Robert-Granié; E. Manfredi Genetic parameters for milk somatic cell score and relationship with production and udder type traits in dairy Alpine and Saanen primiparous goats, Journal of Dairy Science, Volume 94 (2011) no. 7, pp. 3629-3634 | DOI; Arnal et al., 2018 M. Arnal; C. Robert-Granie; H. Larroque Diversity of dairy goat lactation curves in France, Journal of Dairy Science, Volume 101 (2018) no. 12, pp. 11040-11051 | DOI), 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.

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^p_L-, W^m_L-, W^p_H-, W^m_H-) and high depletion clusters (W^p_H+, W^m_H+). 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_k) than at the beginning of the subsequent gestation (BW₂₁₀). 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_k was higher than BW₂₁₀ at the end of lactation. BW is easy to measure on farm to monitor animals, especially to quantify energy balance (Thorup et al., 2012 V. M. Thorup; D. Edwards; N. C. Friggens On-farm estimation of energy balance in dairy cows using only frequent body weight measurements and body condition score, Journal of Dairy Science, Volume 95 (2012) no. 4, pp. 1784-1793 | DOI). 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^p_{H- ,} W^m_H). They were generally heavier than Alpine goats (Sauvant et al., 2012 D. Sauvant; S. Giger-Reverdin; F. Meschy; L. Puillet; P. Schmidely Actualisation des recommandations alimentaires pour les chèvres laitières, Productions Animales, Volume 25 (2012) no. 3, pp. 259-276 | DOI).

For lumbar and sternal BCS we found three clusters for all parities. First for all parities, we found high depletion clusters for lumbar (LU_M+, LU_H+) and sternal (ST_M+, ST_H+ ) BCS. Then, we found low depletion clusters for lumbar (LU_M) and sternal (ST_M) 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_M+ , 0.3 points for LU_H+,and 0.5 points for ST_M+ and ST_H+. 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 personel carrying out the BCS evaluation (although differences in herd management, or a random distribution linked to the clustering approach cannot be excluded).

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² test was significant for associations but the Cramer’s V showed weak to moderate values (globally less than 0.4) (Kotrlik et al., 2011 J. W. Kotrlik; J. C. Atherton; H. A. Williams; M. Khata Jabor Reporting and Interpreting Effect Size in Quantitative Agricultural Education Research, Journal of Agricultural Education, Volume 52 (2011) no. 1, pp. 132-142 | DOI). 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 S. S. Waltner; J. P. McNamara; J. K. Hillers Relationships of Body Condition Score to Production Variables in High Producing Holstein Dairy Cattle, Journal of Dairy Science, Volume 76 (1993) no. 11, pp. 3410-3419 | DOI; Roche et al., 2007 J. Roche; K. Macdonald; C. Burke; J. Lee; D. Berry Associations Among Body Condition Score, Body Weight, and Reproductive Performance in Seasonal-Calving Dairy Cattle, Journal of Dairy Science, Volume 90 (2007) no. 1, pp. 376-391 | DOI), whereas other studies did not find any relationship between these variables (Garnsworthy & Topps, 1982 P. C. Garnsworthy; J. H. Topps The effect of body condition of dairy cows at calving on their food intake and performance when given complete diets, Animal Production, Volume 35 (1982) no. 1, pp. 113-119 | DOI; Garnsworthy & Jones, 1987 P. C. Garnsworthy; G. P. Jones The influence of body condition at calving and dietary protein supply on voluntary food intake and performance in dairy cows, Animal Production, Volume 44 (1987) no. 3, pp. 347-353 | DOI). 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 F. Blanc; F. Bocquier; J. Agabriel; P. D'Hour; Y. Chilliard Adaptive abilities of the females and sustainability of ruminant livestock systems. A review, Animal Research, Volume 55 (2006) no. 6, pp. 489-510 | DOI; Friggens et al., 2017 N. Friggens; F. Blanc; D. Berry; L. Puillet Review: Deciphering animal robustness. A synthesis to facilitate its use in livestock breeding and management, Animal, Volume 11 (2017) no. 12, pp. 2237-2251 | DOI). 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 B. Dumont; L. Puillet; G. Martin; D. Savietto; J. Aubin; S. Ingrand; V. Niderkorn; L. Steinmetz; M. Thomas Incorporating Diversity Into Animal Production Systems Can Increase Their Performance and Strengthen Their Resilience, Frontiers in Sustainable Food Systems, Volume 4 (2020) no. July, pp. 1-15 | DOI). 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² 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^m_H 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 M. Arnal; C. Robert-Granie; H. Larroque Diversity of dairy goat lactation curves in France, Journal of Dairy Science, Volume 101 (2018) no. 12, pp. 11040-11051 | DOI). 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^m_L- ) in parity 2, while W^p_H+ switched to the W^m_L- cluster. Goats in the W^p_H+ 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^m_H- profile remaining in this cluster in successive lactations. Half of the goats in the W^m_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_H+ cluster remained in this cluster in successive lactations. More than half of the goats in the ST_M+ and ST_M 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 T. J. Hayden; C. R. Thomas; I. A. Forsyth Effect of Number of Young Born (Litter Size) on Milk Yield of Goats: Role for Placental Lactogen, Journal of Dairy Science, Volume 62 (1979) no. 1, pp. 53-57 | DOI; Zamuner et al., 2020 F. Zamuner; K. DiGiacomo; A. Cameron; B. Leury Effects of month of kidding, parity number, and litter size on milk yield of commercial dairy goats in Australia, Journal of Dairy Science, Volume 103 (2020) no. 1, pp. 954-964 | DOI). 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 L. Puillet; O. Martin A dynamic model as a tool to describe the variability of lifetime body weight trajectories in livestock females, Journal of Animal Science, Volume 95 (2017) no. 11, pp. 4846-4856 | DOI). 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 I. Palhière; C. Oget; R. Rupp Functional longevity is heritable and controlled by a major gene in French dairy goats, Proceedings of the World Congress on Genetics Applied to Livestock Production, 11-16 February 2018 (2018), p. 165 | DOI),which increases replacement costs.

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.

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 A. Ben Abdelkrim; L. Puillet; P. Gomes; O. Martin Lactation curve model with explicit representation of perturbations as a phenotyping tool for dairy livestock precision farming, Animal, Volume 15 (2020) no. 1, p. 100074 | DOI; Poppe et al., 2020 M. Poppe; R. F. Veerkamp; M. L. van Pelt; H. A. Mulder Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding, Journal of Dairy Science, Volume 103 (2020) no. 2, pp. 1667-1684 | DOI; Adriaens et al., 2021 I. Adriaens; I. van den Brulle; L. D'Anvers; J. M. Statham; K. Geerinckx; S. De Vliegher; S. Piepers; B. Aernouts Milk losses and dynamics during perturbations in dairy cows differ with parity and lactation stage, Journal of Dairy Science, Volume 104 (2021) no. 1, pp. 405-418 | DOI).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 S. Lerch; A. De La Torre; C. Huau; M. Monziols; C. Xavier; L. Louis; Y. Le Cozler; P. Faverdin; P. Lamberton; I. Chery; D. Heimo; C. Loncke; P. Schmidely; J. A. Pires Estimation of dairy goat body composition: A direct calibration and comparison of eight methods, Methods, Volume 186 (2021) no. May 2020, pp. 68-78 | DOI). 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 S. Lerch; A. De La Torre; C. Huau; M. Monziols; C. Xavier; L. Louis; Y. Le Cozler; P. Faverdin; P. Lamberton; I. Chery; D. Heimo; C. Loncke; P. Schmidely; J. A. Pires Estimation of dairy goat body composition: A direct calibration and comparison of eight methods, Methods, Volume 186 (2021) no. May 2020, pp. 68-78 | DOI). 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 I. Adriaens; N. C. Friggens; W. Ouweltjes; H. Scott; B. Aernouts; J. Statham 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, Journal of Dairy Science, Volume 103 (2020) no. 8, pp. 7155-7171 | DOI).

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.

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 A. Atzori Understanding milk and body reserves trajectories and nutrient partitioning in dairy goats through a modelling approach, Peer Community in Animal Science (2024) | DOI).

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).

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 and scripts/code are deposited on Recherche Data Gouv (https://doi.org/10.57745/C1XPQ2; Gafsi, 2023a)

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