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

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


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.

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

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

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

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

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).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 p L-) that produced 155.6 kg less over the lactation than a medium-level cluster (Y p M-).
-a medium persistency cluster with the highest total milk production level that gathered 22.6% of the primiparous (Y p H).
-the highest persistency cluster with a low total milk production level that gathered 14.1% of the primiparous (Y p L+).
Peer Community Journal, Vol. 4 (2024), article e85 https://doi.org/10.24072/pcjournal.449Multiparous 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 m M+) that maintained 20.4 % more the production than a low persistency cluster (Y m M-).
-the highest total milk production level cluster with a medium persistency (Y m H) that gathered 34.6 % of the population.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.

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).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 p L-) that averaged 10.0 kg less at kidding than a high-level cluster (W p H-).Those profiles had a higher BW210 than BWk.-the highest depletion cluster with a high BW level at kidding (W p H+) that gathered 31.4% of the population.Despite having the highest repletion speed, this cluster presented a lower BW210 than BWk due to the high level of depletion, that is not totally compensated at 210 days of lactation.
Peer Community Journal, Vol. 4 (2024), article e85 https://doi.org/10.24072/pcjournal.449Multiparous 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 m L-) that averaged 17.6 kg less at kidding than a high-level cluster (W m H-).For these clusters BW210 was lower than BWk.-the highest depletion cluster with a high BW level at kidding (W m H+) that gathered 26.6% of the multiparous.Despite having the highest repletion speed, this profile presented a lower BW210 than BWk due to the high level of depletion, that is not totally compensated at 210 days of lactation.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 Hcluster.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_Lmin and BCS_Lk) 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 Peer Community Journal, Vol. 4 (2024), article e85 https://doi.org/10.24072/pcjournal.449total 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).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 (STM+) that averaged 0.7 points less at kidding than a highlevel cluster (STH+).STM+ 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 (STM).STM cluster presented the lowest repletion speed.For lumbar BCS, Pradel's Alpine goats were more represented in the LUH+ cluster, whereas Grignon's Alpine goats were more represented in the LUM cluster.Grignon's Saanen goats were more represented in the LUM 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 STM cluster, whereas Grignon's Alpine goats were more represented in the STH+ cluster.Grignon's Saanen goats were more represented in the STH+ cluster.Primiparous represented between 30 % and 35 % of the population in each profile for sternal BCS.See Appendix B section 3 for more details.

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

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 p L-with STM+ and W p H+ with STM accounted for the highest proportion of goats with 18.8 % of the population, followed by the association W p L-with STM with 17.9 % of the population.The association W p H-with STM+ 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 STM accounted for the highest proportion of goats with 18.6 % of the population, followed by the association W m Lwith STM+ with 14.2 % of the population.The association W m H-with STM+ 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.

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 LUM+ with STM+ accounted for the highest proportion of goats with 21.4 % of the population, followed by the association  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 LUM+ with STM+ accounted for the highest proportion of goats with 20.0 % of the population, followed by the association LUH+ with STM with 18.6 % of the goats.The association LUM+ with STH+ 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.

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

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² 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² 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 STH+ profile remained in this cluster in successive lactations, while the other part switched to other clusters.More than half of the goats in the STM+ and STM 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.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 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;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 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 (BWk) than at the beginning of the subsequent gestation (BW210).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 BWk was higher than BW210 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 p H-, W m H).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 (LUM+, LUH+) and sternal (STM+, STH+ ) BCS.Then, we found low depletion clusters for lumbar (LUM) and sternal (STM) 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 LUM+ , 0.3 points for LUH+,and 0.5 points for STM+ and STH+.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 LUM and STH+ 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).

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² test was significant for associations but the Cramer's V showed weak to moderate values (globally less than 0.4) (Kotrlik 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 Peer Community Journal, Vol. 4 (2024), article e85 https://doi.org/10.24072/pcjournal.449productive 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).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 STH+ cluster remained in this cluster in successive lactations.More than half of the goats in the STM+ and STM 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 MYpeak, 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.

Figure 1 -
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: MYpeak = highest milk yield value; MY210 = 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 = (MY250-MY150 )/MY150) x100.

Figure 3 -
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: BWk = body weight at kidding; BWmin = minimum body weight; BW210 = body weight at 210 days; Dep_speedkà30 : Body weight depletion speed between kidding and 30 days = (BW30 -BWk)/ 30; Rep_speed180à210: Body weight repletion speed between 180 and 210 days = (BW210 -BW180)/ 30.
Abbreviations: MY = milk yield; BW = body weight; BCS = body condition score

Figure 4 -
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) (MYpeak = highest milk yield value; MY210 = 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 = (MY250-MY150/MY150) x100; Y p L+= Low milk yield and high persistency cluster for primiparous; Y p L-= Low milk yield and low persistency cluster for primiparous; Y p M-= Medium milk yield and low persistency cluster for primiparous; Y p H = High milk yield and medium persistency cluster for primiparous; Y m M+= Medium milk yield and high persistency cluster for multiparous; Y m M-= Medium milk yield and low persistency cluster for multiparous; Y m H = High milk yield and medium persistency cluster for multiparous).

Figure 5 -
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) (BWk = body weight at kidding; BWmin = minimum body weight; BW210 = body weight at 210 days; Dep_speedkà30: Body weight depletion speed between kidding and 30 days = (BW30 -BWk )/ 30; Rep_speed180à210: Body weight repletion speed between 180 and 210 days = (BW210 -BW180 )/ 30; W p L-= Low body weight and low depletion cluster in primiparous; W p H+ = High body weight and high depletion cluster in primiparous; W p H-= High body weight and low depletion cluster in primiparous; W m L-= Low body weight and low depletion cluster in multiparous; W m H+= High body weight and high depletion cluster in multiparous; W m H-= High body weight and low depletion cluster in multiparous).
2 Y m M+= Medium milk yield and high persistency cluster; Y m M-= Medium milk yield and a low persistency cluster; Y m H = High milk yield and a medium persistency cluster; W m L-= Low body weight and low depletion cluster; W m H+= High body weight and high depletion profile cluster; W m H-= High body weight and low depletion cluster.
2 W m L-= Low body weight and low depletion cluster; W m H+= High body weight and high depletion cluster; W m H-= High body weight and low depletion cluster; STM+ =Medium sternal body condition score and depletion cluster; STM =Medium sternal body condition score and low depletion cluster; STH+ =High sternal body condition score and depletion cluster.

Figure 7 -
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+= Low milk yield and high persistency cluster for primiparous; Y p L-= Low milk yield and low persistency cluster for primiparous; Y p M-= Medium milk yield and low persistency cluster for primiparous; Y p H = High milk yield and medium persistency cluster for primiparous; Y m M+ = Medium milk yield and high persistency cluster for multiparous; Y m M-= Medium milk yield and low persistency cluster for multiparous; Y m H = High milk yield and medium persistency cluster for multiparous).

Figure 8 -
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-= Low body weight and low depletion cluster; W p H+= High body weight and high depletion cluster; W p H-= High body weight and low depletion cluster; W m L-= Low body weight and low depletion cluster; W m H+ = High body weight and high depletion cluster; W m H-= High body weight and low depletion cluster).

Figure 9 -
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 (STM+ =Medium sternal body condition score and depletion profile; STM =Medium sternal body condition score and low depletion profile; STH+ = High sternal body condition score and depletion profile).

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.

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

Table 4a -
Statistical description of synthetic indicators for MY clusters in primiparous goats.Y p L+ = Low milk yield and high persistency cluster; Y p L-= Low milk yield and low persistency cluster; Y p M-= Medium milk yield and low persistency cluster; Y p H = High milk yield and medium persistency cluster. 3

Table 4b -
Statistical description of synthetic indicators for MY clusters in multiparous goats.
a-c Means with superscripts differ significantly by row.1 SumMY = sum of daily milk yield values over 250 days; MYpeak = highest milk yield value; MY210 = milk yield value at 210 days; Peak time = time of the highest milk yield value; Persistency = (MY250-MY150/MY150) x100. 2 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). 3Y m M+ = Medium milk yield and high persistency cluster; Y m M-= Medium milk yield and low persistency cluster; Y m H = High milk yield and medium persistency cluster.

Table 5a -
Statistical description of synthetic indicators for BW clusters in primiparous goats.W p L-= Low body weight and low depletion cluster; W p H+= High body weight and high depletion cluster; W p H-= High body weight and low depletion cluster. 3

Table 5b -
Statistical description of synthetic indicators for BW clusters in multiparous goats.
-a group of depletion clusters with two different lumbar BCS level at kidding (68.7 % of the population): a medium level cluster (LUM+) that averaged 0.4 points less at kidding than a highlevel cluster (LUH+).LUM+ 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 (LUM).LUM cluster presented the same repletion speed than LUH+.

Table 6 -
Statistical description of synthetic indicators for lumbar BCS clusters in goats.
2 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).3LUM+ = Medium lumbar body condition score and depletion cluster; LUM = Medium lumbar body condition score and low depletion cluster; LUH+ = High lumbar body condition score and depletion cluster

Table 7 -
Statistical description synthetic indicators for sternal BCS clusters in goats.

Table 8a -
Contingency table displaying the frequency of individual primiparous lactations affected to MY and BW clusters (see section 2 for clustering methodology).

Table 8b -
Contingency table displaying the frequency of individual multiparous lactations affected to MY and BW clusters (see section 2 for clustering methodology).

Table 9a -
Contingency table displaying the frequency of individual primiparous lactations affected to BW and sternal BCS clusters (see section 2 for clustering methodology).
1 % = proportion of lactations among the 448 primiparous goats. 2 W p L-= Low body weight and low depletion cluster; W p H+= High body weight and high depletion cluster; W p H-= High body weight and low depletion cluster; STM+ =Medium sternal body condition score and depletion cluster; STM =Medium sternal body condition score and low depletion cluster; STH+ =High sternal body condition score and depletion cluster.

Table 9b -
Contingency table displaying the frequency of individual multiparous lactations affected to BW and BCS sternal clusters (see section 2 for clustering methodology).
LUH+ with STM with 19.5 % of the goats.The association LUM with STM+ 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 10a -
Contingency table displaying the frequency of individual primiparous lactations affected to BCS lumbar and BCS sternal clusters (see section 2 for clustering methodology).% = proportion of lactations among the 374 primiparous goats. 2 LUM+ = Medium lumbar body condition score and depletion cluster; LUM = Medium lumbar body condition score and low depletion cluster; LUH+ =High lumbar body condition score and depletion cluster; STM+ =Medium sternal body condition score and depletion cluster; STM = Medium sternal body condition score and low depletion cluster; STH+ =High sternal body condition score and depletion cluster. 1

Table 10b -
Contingency table displaying the frequency of individual multiparous lactations affected to lumbar and sternal BCS clusters (see section 2 for clustering methodology).Medium lumbar body condition score and low depletion cluster; LUH+ =High lumbar body condition score and depletion cluster; STM+ =Medium sternal body condition score and depletion cluster; STM = Medium sternal body condition score and low depletion cluster; STH+ =High sternal body condition score and depletion cluster.
1 % = proportion of lactations among the 740 multiparous goats. 2 LUM+ = Medium lumbar body condition score and depletion cluster; LUM = 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 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.