J. Dairy Sci.

101:11040–11051
https://doi.org/10.3168/jds.2018-14980
© 2018, The Authors. Published by FASS Inc. and Elsevier Inc. on behalf of the American Dairy Science Association®.
This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Diversity of dairy goat lactation curves in France

M. Arnal,*†1 C. Robert-Granié,* and H. Larroque*
*GenPhySE, Université de Toulouse, INRA, INPT, ENVT, 31326 Castanet-Tolosan, France
†Institut de l’élevage (Idele), 31326 Castanet-Tolosan Cedex, France

ABSTRACT health, as well as the time–course scale of milk produc-

tion. Indeed, a high level of production at the peak can
A high level of production at the peak of lactation lead to an energy deficit for the animal if it is unable
may be associated with animal health disorders, high to ingest enough food to compensate. Goats, like dairy
feeding costs, and reduced milk supply throughout the cattle, will then begin drawing on their body reserves,
year. The objective of this study was to typologize the which may lead to metabolic and reproductive disorders
lactation curves in French dairy goats and analyze the (Gipson and Grossman, 1990). A goat in energy deficit
influence of environmental and genetic factors on these has to be fed expensive supplement concentrates and
curves. The data set consisted of 2,231,720 monthly high nutrition value forage, whereas a persistent goat
test-day records of 213,534 French Saanen and Alpine can be fed lower-quality forage and fewer concentrates
goats recorded between September 2008 and June 2012. (Sölkner and Fuchs, 1987), which is an asset to farm
First, principal component analysis classified the shape systems aiming to reduce inputs. The ability to main-
of the lactation curves into 3 principal components: the tain lactation in a partially season-cycled production
first component accounted for milk yield level through- context meets the processing industry demand for more
out lactation, the second component accounted for lac- evenly spread milk production, in line with the needs
tation persistency, and the third component accounted of the market (Kearney, 2010), while also avoiding the
for milk yield in mid-lactation. Then, from the principal need to freeze cheeses, which can change their flavor
component scores, the lactations were clustered into 5 (Van Hekken et al., 2005).
different groups. Most lactations had a similar shape Research on shape of the lactation curve has mainly
to the mean curve, except 30% of the lactations that addressed milk persistency, trying to find indicator cri-
fell into 3 clusters that had a high production level at teria of the slope of the curve after the peak of lactation
the peak and then a different persistency according to (Gengler, 1996), or modeling of overall curve shape.
cluster. Estimated breeding value for milk yield and Some authors have focused on finding mathematical
home region of breeding were the factors most related models that describe the biological processes of milk
to lactation production level. Month of kidding, breed, production by mammary gland cells (Pollott, 2000;
and gestation stage had the biggest effect on persis- Elvira et al., 2013a). Others have proposed a statisti-
tency. Month of kidding was the factor most strongly cal approach by modeling the shape of the curve from
linked to mid-lactation production. A herd effect was elementary controls such as random regression test-
observed on all 3 principal components. day models including genetic and nongenetic effects
Key words: French dairy goat, test-day milk record, (Menéndez-Buxadera et al., 2010; Mucha et al., 2014;
classification, lactation curve Brito et al., 2017) or by using methods such as principal
component analysis (PCA) to sum up the information
INTRODUCTION contained in the data (Carta et al., 2014). Principal
component analysis is a dimension-reduction tool that
A goat’s level of milk production is often evaluated reduces a large set of variables to a small set contain-
by its daily average yield or by its total production ing most of the original information. Macciotta et al.
throughout lactation. However, neither approach fac- (2006) used a PCA based on rank of bovine test-days
tors in the change in milk production according to DIM where the first principal component (PC) summarized
(i.e., lactation curve). Shape of the lactation curve can the averaged milk yield level and the second component
affect an animal’s dietary needs, and consequently its summarized persistency. van der Werf et al. (1998) and
Druet et al. (2003) performed data dimension reduction
on a covariance genetic matrix derived from a random
Received April 26, 2018.
Accepted July 23, 2018. regression test-day model and also obtained 2 PC with
1
Corresponding author: mathieu.arnal@​idele​.fr the same biological interpretation.

11040
 LACTATION CURVES IN FRENCH DAIRY GOATS 11041

Bouloc (1991) showed in French dairy goats that lac- French Saanen and Alpine goats in 910 French herds
tation curve shapes may depend on breed and genetics spread throughout France.
of the animal. Menéndez-Buxadera et al. (2010) showed As described by Pulina et al. (2018), in France goats
in Murciano-Granadina goats that persistency is heri- production systems present a great diversity according
table (0.20). Curve shape may also vary with breeding to feeding (grazing or not, concentrate quantity), repro-
and physiological factors such as parity, age at kidding, duction (out of season or not), transformation of milk
kidding season, dry-period length, gestation stage, and (cheese maker/deliverer), herd size, or even breeding
herd environment factors (Gipson and Grossman, 1990; goal (milk yield/composition). However, the database
Bouloc, 1991). used recording animal milk productions and did not
To our knowledge, no recent study has been done on contain detailed information on production system of
lactation curve shapes using test-day records in French herds. Almost 60% of the lactations retained were on
dairy goats. The latest study, performed 27 yr ago, was Alpine goats, and mainly (77.5%) in the north-west of
by Bouloc (1991), who classified protein yield curves France. The final data set contained 77% of lactations
based on different input features: lactation duration, in parity 1 to 3, with other lactations up to parity 8
average daily protein yield, persistency of the protein (Table 1). Lactations lasted 288 d on average, for an av-
yield calculated as the ratio of d 50–d 200-period pro- eraged total production of 964.1 kg of milk (minimum:
tein yield to d 31–d 60-period protein yield. 132.5 kg, maximum: 2,615 kg, SD: 253).
Here we typologized lactation curve shapes based on Goats were measured on average 7 times per lacta-
PCA scores taking into account test-day lactation stage tion (minimum: 4, maximum: 11, SD: 1), at an average
in the 2 main French dairy goat breeds (Alpine and measure-to-measure interval of 39 d (minimum: 10,
Saanen), using a recent large data set of test-day milk maximum: 158, SD: 10.9). Goats were not all measured
production records from the French genetic national at the same DIM, and production varied from 0.5 to
database. We went on to study the influence of different 11.4 kg per test-day record.
genetic and environmental factors on these curves. For each lactation, parturition date served to calcu-
late age at kidding and define kidding month. For goats
in their second lactation or more, “dry-period length”
MATERIALS AND METHODS
was defined as the length of the period between the end
Data of the previous lactation and parturition date. Breeders
give the end date of the previous lactation; otherwise,
The data were extracted from the genetic database it is a deducted date by adding to the date of the last
used for the national French genetic evaluations. France test day 14 d (or 28 d if the interval with the following
has been running a goat milk yield recording scheme lactation is greater than 56 d). Gestation stage at 300
since 1989 that applies the standard A, D, AT, AZ, or DIM was calculated for goats having subsequent lacta-
CZ methods registered by ICAR (2017) with sampling tion.
at 4- or 5-wk intervals. The data retained here was
collected from the A, D, and AZ protocols with precise Analysis
measurement of daily production, taking into account
amount of milk at the 2 daily milkings. About 60% of All statistical analyses were performed using R
the overall database comes from these protocols. The software (RCore Team, R Foundation for Statistical
weights of milk produced by a goat in the morning Computing, 2014, R: A Language and Environment for
and evening milkings were added together to obtain the Statistical Computing, Version 3.0.2, Vienna, Austria,
weight of daily milk. http:​/​www​.R​-project​.org).
Lactations had to have a conventional duration be- The general approach of the analyses was to perform
tween 180 and 350 d to be conserved. Test-day records (1) a PCA on test-days after smoothing the data, then
before 7 and beyond 300 DIM were removed to have a (2) a classification of the curves, and finally, (3) a study
sufficient number of records per DIM. Within a lacta- of the relationships between environmental and genetic
tion, a goat was required to have at least 4 test-day factors and curve shapes.
records, with the first test day recorded before d 81 of Principal Component Analysis. The PCA was
lactation to assess production at the beginning of lacta- performed on the test-day milk yield data, for each lac-
tion, the peak of lactation occurring around 40 DIM. tation, using the R package fdapace (Dai et al., 2016).
Only goats from herds having at least 30 conserved Fdapace was well-adapted to the structure of our data
lactations per year were kept in the analysis. set: sparse data, irregular measurement stage, irregular
The final data set consisted of 2,231,720 test-day number of measurements per lactation, small number
records collected from 319,975 lactations of 213,534 of repeated measurements available per lactation. First,
Journal of Dairy Science Vol. 101 No. 12, 2018
11042 ARNAL ET AL.

Table 1. Percentage of data for each level of different factors in French fdapace calculated the smoothed mean curve using lo-
dairy goats
cal linear Gaussian kernel regression between 51 equi-
Variable % distant support points (default option) aggregating all
Parity
the measurements together. Then, the raw covariance
 1 32.9 for each curve was calculated, and all covariances were
 2 26.5
 3 17.8
aggregated to generate the sample raw covariance. The
 4 10.9 smoothed covariance was estimated from the off-diago-
 5 6.1
 6 3.3
nal of sample raw covariance. Eigen analyses performed
 7 1.7 on the smoothed covariance matrix gave the estimated
 8 0.8
Breed
eigenfunctions, eigenvalues, and PC scores (Yao et al.,
 Alpine 59.1 2005).
 Saanen 40.9
Milk EBV (kg)
To evaluate the accuracy of the smoothed curves us-
  −349 to −79 25.1 ing local linear regression, the true and predicted milk
  −79 to −33 25.3
  −33 to 13 24.7
yields were compared by studying the coefficient of de-
  13 to 301 24.9 termination (R2). This analysis can be done according
Region1
 NE 17.7
to the number of PC.
 NW 77.5 Classification. Using the PC scores (SPC) of the
 S 4.9
Kidding month
lactations, an unsupervised classification, also called
 1 24.51 cluster analysis, was performed using the R package
 2 27.89
 3 14.24
rmixmod (Langrognet et al., 2016). Rmixmod is a tool
 4 3.50 to fit a mixture model of multivariate components to a
 5 1.13
 6 0.19
data set (Biernacki et al., 2006). A new feature of the
 7 0.03 Rmixmod software is that it considers a parameter-
 8 1.06
 9 9.02
ization of the variance matrix of a cluster through its
 10 8.21 eigenvalue decomposition, leading to many meaningful
 11 6.83
 12 3.40
models for clustering (Celeux and Govaert, 1995). Each
Kidding age2 (mo) variance matrix is decomposed as the product of one
 Other 12
 P4+ 22.8
parameter determining the volume of the cluster, a
 Normal 50.8 matrix of eigenvectors determining its orientation, and
 Young 2.9
 Old 11.6
a diagonal matrix determining its shape (Lebret et al.,
Dry-period length3 (d) 2015). By allowing some but not all of these parameters
 1–60 19.5
 60–90 32.7
and matrices to vary between clusters, parsimonious
 90–120 9.1 and easily interpreted models can be obtained. As SPC
  Other lengths 5.8
  Parity 1 32.9
are continuous variables, we tested 14 different Gauss-
SCC EBV4 (base 100) ian models with different variance–covariance struc-
 101–105 23.8
 105–131 21.2
tures as proposed by Celeux and Govaert (1995), and
 65–97 26.2 we selected the best one based on the Bayesian informa-
 97–101 25
  EBV missing 3.8
tion criterion (BIC). The number of components of the
Gestation stage5 (d) mixture model is generally known or fixed beforehand.
 120–150 5.2
 150–180 1.3
If this number is not known, it is possible to compare
 1–60 5.9 models with different numbers of components and to
 60–90 32.5
 90–120 24
use criteria proposed in the Rmixmod package to make
  Other stages 1.3 the choice. Here, the optimal number of components
  No consecutive lactation 29.9
was researched between 2 and 5 clusters, and the BIC
1

2
French herd region (NW = north-west; S = south; NE = north-east). was used to determine it.
Age at kidding for the first 3 parturitions (young = 9–10 mo at first kidding,
16–22 mo at second kidding, or 22–34 mo at third kidding; normal = 11–13 mo The mixture parameter was estimated through
at first kidding, 23–25 mo at second kidding, or 35–37 mo at third kidding; old maximization of the log-likelihood. In a first step, a
= 14–30 mo at first kidding, 26–56 mo at second kidding, or 38–77 mo at third
kidding; P4+ = parity 4 and more; other = other ages at kidding for the first stochastic expectation maximization algorithm was
3 parturitions). used 100 times with 1,000 iterations each time. Then,
3
Dry-period length before last kidding (d; for goats in second and more lacta- an expectation maximization algorithm was performed
tion).
4
Somatic cell count EBV (base 100): <100: increase in SCC, >100: decrease in 100 times with 1,000 iterations each time. The stochas-
SCC. tic expectation maximization algorithm used random
5
Gestation stage (d) at 300 DIM. drawing at each iteration, which avoided converging to

Journal of Dairy Science Vol. 101 No. 12, 2018

 LACTATION CURVES IN FRENCH DAIRY GOATS 11043

stationary points of the log-likelihood (Biernacki et al., from the residuals of model 1, due to the total confu-
2003). sion between herd number and region.
Relationships Between Environmental and
Genetic Factors and PC. The effects of environmen- RESULTS AND DISCUSSION
tal and genetic factors on the first 3 PC were studied
according to the following model: Principal Component Analysis

The first 3 components represented 99.5% of the

yijklmnopqr = µ + Bi + Rj + N k + Al + M m + Dn + Go + I p original variance: 81% for PC1, 15%, for PC2, and 3.5%

+ SCI q + εijklmnopqr , [model 1] for PC3. It was not considered useful to further study
PC due to the total variance already explained. The R2,
used to evaluate the precision of the fit by local linear
where yijklmnopqr is the dependent variable (SPC) and µ regression, was 0.92 with 3 PC and fell to 0.88 with
is overall mean. The fixed effects tested were Bi, breed only the first 2 PC. Clearly, adjustment was better with
(Saanen or Alpine); Rj, French herd region (NW = all 3 components, especially considering the smoothed
north-west; S = south; NE = north-east); Nk, parity (1 lactation curves of extreme goats after peak lactation
to 8); Al, age at kidding for the 3 first parturitions (5 (Figure 1). Therefore, 3 PC were retained.
levels: 1 = 9–10 mo at first kidding, 16–22 mo at second Plotting extreme lactation curves, with SPC extremes
kidding, or 22–34 mo at third kidding; 2 = 11–13 mo at for one PC and other SPC close to 0, showed that each
first kidding, 23–25 mo at second kidding, or 35–37 mo PC characterized a component of curve shape (Figure
at third kidding; 3 = 14–30 mo at first kidding, 26–56 1). Principal component 1 was an indicator of level
mo at second kidding, or 38–77 mo at third kidding; 4 of milk production throughout lactation, where a low
= parity 4 and more; 5 = other ages at kidding for the SPC1 corresponds to low milk production (mean ± SD:
first 3 parturitions); Mm, kidding month (12 levels: 1 = SPC1: −1.4 ± 144.4, total lactation production (kg):
January, …, 12 = December); Dn, class of dry-period 964.1 ± 253.3; correlation between SPC1 and lactation
length before the last kidding (d) for goats in second or
more lactations (5 levels: 1 = [1–60], 2 = [60–90], 3 =
[90–120], 4 = other lengths, 5 = 1st parity); Go, class
of gestation stage (d) at 300 DIM (7 levels: 1 = [1–60],
2 = [60–90], 3 = [90–120], 4 = [120–150], 5 = [150–180
d], 6 = other stages, 7 = no consecutive lactation); Ip,
class of genetic EBV based on total-lactation milk yield
(kg) from official French genetic evaluation (4 levels: 1
= [−349 to −79], 2 = [−79 to −33], 3 = [−33 to 13],
4 = [13 to 301]); SCIq, class of SCC EBV (base 100):
<100: increased SCC, >100: decreased SCC) from of-
ficial French genetic evaluation (5 levels: 1 = [65–97], 2
= [97–101], 3 = [101–105], 4 = [105–131], 5 = missing
information); and εijklmnopqr are the normally distributed
residuals. Model 1 using SPCx with x = {1, …, X} is
noted MSPCx. Percentages of data for each class of the
different factors are presented in Table 1. Despite the
low percentages of data, particularly for the kidding
month effect in summer, the lowest number of data in
a class was 87.
Each effect of model 1 was successively removed to
study its relative effect using a Fisher test and compar-
ing the R2 of the model.
Least squares means were calculated using propor-
tions of the population as weights to take into account
the unbalanced number of data in each level of effects,
using the R package lsmeans (Lenth, 2016). Figure 1. Local linear regressions (curves) for lactations having a
minimum (min) or a maximum (max) score for 1 principal component
The influence of 2-way interactions between herd and (PC; first, PC1; second, PC2; or third, PC3) and scores close to 0 for
year of lactation (2,145 classes) on SPC was studied the 2 remaining principal components.

Journal of Dairy Science Vol. 101 No. 12, 2018

11044 ARNAL ET AL.

Table 2. Mean (SD) per cluster for the first, second, and third principal component (PC1, PC2, and PC3,
respectively), and proportion of lactations per cluster

Cluster
number PC1 PC2 PC3 Proportion (%)
1 −76.5 (94.4) 107.1 (40.7) −39.5 (17.7) 4.2
2 −121.7 (80.9) 7.2 (40.8) 9.7 (18.2) 39
3 221.9 (120.5) −41 (69.3) −23.8 (30.5) 8.9
4 101.1 (108.7) 64.2 (59.5) 29.4 (26.5) 16.2
5 41.7 (78) −29 (42.1) −4.9 (19) 31.7

production: 1); PC2 was an indicator of the difference BIC criteria were 10,775,614, 10,767,801, 10,763,784,
in milk production between the peak and end of lacta- and 10,760,253 for 2, 3, 4, and 5 clusters, respectively.
tion, where a low SPC2 was equivalent to high milk The best model was then obtained by considering 5
production at the end of the lactation (mean ± SD: clusters. To obtain clusters with sufficient sizes, it was
SPC2: 4.9 ± 61.9, production difference (kg) between not desirable to have a greater number of clusters. On
DIM 40 and 240: 1.11 ± 0.9; correlation between SPC2 the whole data file, all iterative procedures used to
and production difference: 0.99); finally, PC3 was an build the best model were run several times, and the
indicator of milk production after the peak of lactation same clusters were found each time.
at around 120 DIM: a low SPC3 equated to low produc- Table 2 reports mean and standard deviation per
tion around 120 DIM (mean ± SD: SPC3: 3.3 ± 27.2, cluster of the lactation scores for each PC. Figures 2a
production difference between DIM 120 and the sum and 2b represent each lactation according to their PC
between the first and the last DIM (kg): −1.6 ± 0.7; scores and their cluster membership. Each ellipse rep-
correlation between SPC3 and production difference: resents the internal covariance between 2 PC within a
0.97). cluster.
On cow data, some authors have defined persistency Clusters 3 and 4, with the biggest ellipses in Fig-
as a degree of maintenance of peak milk yield (Gross- ures 2a and 2b, had the largest SPC variances, with
man et al., 1999), others as the degree of decline in standard deviations of 120.5 and 108.7, respectively, for
milk yield after the peak, using indicators of the speed PC1. In contrast, clusters 2 and 5, although including
of decline such as the ratio of milk yield in the first 100 proportionally more lactations (39 and 31.7%, respec-
DIM to milk yield in the last 100 DIM (Sölkner and tively), had the lowest SPC variances (80.9 and 78,
Fuchs, 1987). Based on these definitions of persistency, respectively, for PC1, for example).
PC2 and PC3 could both be considered persistency The averages of daily milk production were plotted
indicators, where “complete” persistency could be a for all the data as well as according to cluster number
minimal value for PC2 and a maximal value for PC3. (Figure 3). In each cluster, average daily production
For the rest of the study, we defined persistency as a did not fluctuate between 2 consecutive DIM. Each
difference in milk production between the peak and end shape of average daily production looked like a smooth
of lactation. Consequently, PC2 was associated with function.
persistency and PC3 was associated with milk produc- We estimated the probability of each lactation be-
tion in mid-lactation. longing to a cluster. To challenge the robustness of the
Macciotta et al. (2006) performed a PCA on a corre- clustering, lactations with a cluster probability higher
lation matrix from dairy cow test-day milk yields. They than 0.7 were conserved (n = 161,676 lactations; i.e.,
used not DIM but rank of the test-day during lactation 51% of the data set). In this case, the clusters looked
and obtained 2 PC: the first representing total-lactation well separated (the ellipses did not overlap) and there
production, and the second representing persistency. was no substantial change in the means curve profiles
Carta et al. (2014), with a similar approach, obtained (results not shown).
the same results in dairy sheep. According to BIC, the repeatability of the classifica-
tion results, the smoothness of the mean curve of each
Classification of Curve Shapes cluster and the similarity of the mean curves using a
probability higher than 0.7, the choice of 5 clusters ap-
The best Gaussian mixture model for our data based peared robust.
on the BIC was defined by free volume, equal shape, Principal component 3 (representing only 3.5% of
and equal orientation between the different clusters. the original variance) contributed significantly to the
The number of clusters tested varied from 2 to 5. The diversity of the shape of curves highlighted, which is

Journal of Dairy Science Vol. 101 No. 12, 2018

 LACTATION CURVES IN FRENCH DAIRY GOATS 11045

easily understandable as PC3 was responsible for the mean for PC2 corresponding to a “flat” curve (Figure
flexibility in mid-lactation. Indeed, not taking into ac- 3). The mean curves of 3 others clusters (nearly 30%
count PC3 in the classification would make cluster 1 of the data) were very different from the all-data set
disappear and change the shape of clusters 3 and 4 mean curve. Indeed, cluster 3 (the dark blue curve),
(results not shown). characterizing 8.9% of the data, had the highest mean
The average lactation curve could be decomposed milk production (PC1) and highest persistency (PC2)
into 3 phases according to daily milk yield: (1) an initial but the second-lowest value for PC3. Cluster 1 (the red
rapid increase until the 50th DIM, (2) a stabilization curve), which ranked fourth for milk production level,
for 50 d, and finally, (3) a slow decrease through to the had the lowest persistency and the highest decrease in
end of lactation. mid-lactation, resulting in a very specific lactation peak
Day of peak lactation yield differed according to (Figure 3), and constituted the smallest group with
shape of the lactation curve. The lactation peak, de- 4.2% of the lactations. Cluster 4 (the light blue curve
fined as the day of the maximal daily milk production, in Figure 3), which was characterized by the second-
was before 40 DIM for clusters 1 and 5 but was 10 d highest total milk production level and ranked second
later for clusters with the highest milk production level for nonpersistency and first for maintaining milk in
(clusters 3 and 4). Goats in cluster 4 maintained peak mid-lactation, giving an early-bulging convex curve,
milk production for a far longer time than goats in represented 16.2% of lactations.
cluster 1. On a sample of 146,280 French goats of Saanen
Cluster 2 (the green curve in Figure 3), representing breed (31.8%), Alpine breed (56.1%), and other breeds
39% of lactations, was characterized by the lowest milk (12.1%) milked in 1987, Bouloc (1991) realized a clas-
yield level, intermediate means for PC2 and PC3, and a sification based on 3 criteria: (1) persistency of protein
curve shape similar to the mean curve (the black curve yield calculated as the ratio of d 50–d 200-period pro-
in Figure 3). Cluster 5 (the pink curve in Figure 3), tein yield to d 31–d 60-period protein yield, (2) average
representing 31.7% of lactations, was defined by inter- daily protein yield, and (3) duration of lactation. The
mediate means for PC1 and PC3 and the second-lowest author identified 8 distinct clusters and was able to

Figure 2. Plot of the lactations according to (a) the first 2 principal components (PC1 and PC2), and (b) the first and third principal com-
ponents (PC1 and PC3).

Journal of Dairy Science Vol. 101 No. 12, 2018

11046 ARNAL ET AL.

Table 3. Coefficient of determination (R2) of model 1 where the lactation scores of each of the 3 principal
components (PC) are explained by all the environmental factors, and of model 1 successively reduced by 1
environmental factor

Item PC1 PC2 PC3

Model 1 0.402 0.225 0.135
Model 1 without parity 0.399 0.215 0.134
Model 1 without kidding age 0.399 0.222 0.133
Model 1 without kidding month 0.395 0.199 0.034
Model 1 without dry-period length 0.400 0.222 0.135
Model 1 without gestation stage 0.395 0.208 0.130
Model 1 without breed 0.394 0.201 0.135
Model 1 without milk EBV 0.161 0.224 0.135
Model 1 without SCS EBV 0.401 0.219 0.135
Model 1 without region 0.385 0.223 0.134

highlight: (1) lactations with atypical forms (i.e., low ters 3 and 5 (high production and high persistency),
production during the first weeks, which returned to probably due to lower levels of genetics and manage-
an average level thereafter), (2) persistent curves with ment 27 yr ago, nor the profiles of clusters 4 and 1,
a difference of 1.5 kg of milk yield between the peak whereas the atypical curves found by Bouloc (1991)
and end of lactation, (3) curves with marked peaks and were not found here, probably due to the use of differ-
a difference of more than 2.0 kg in milk yield between ent phenotypes for the classification.
the peak and end of lactation. In our data set, only the
mean lactation curve of cluster 2 (low milk yield level Relationship Between Environmental and Genetic
and medium persistency) was similar to what Bouloc Factors and PC
found. Bouloc’s study did not find the profiles of clus-
Model 1 was used to estimate the relative influence of
each environmental and genetic factor on each PC. For
each PC, Table 3 presents the R2 of model 1 and the
R2 of model 1 iterative reduced by one of the factors.
The Fisher’s tests between model 1 and the different
sub-models were all significant (P < 10−5) except for
breed on PC3 (P = 0.08).
The R2 of model 1 was 0.402 for PC1. The biggest
decrease in R2 (R2 = 0.161) occurred when EBV of
total milk yield was removed from model 1, meaning
that it was the factor most related to PC1, far ahead of
region (R2 = 0.385) and breed (R2 = 0.394). For PC2,
the R2 of the model 1 was 0.225. Principal component 2
was not strongly related to environmental and genetic
factors, as the biggest association was found for kidding
month with R2 decreased to 0.199, followed by breed
(R2 = 0.201), gestation stage (R2 = 0.208), parity (R2
= 0.215), and SCS EBV (R2 = 0.219). The R2 of model
1 for PC3 was 0.135. The only factor strongly linked to
this PC was kidding month, with R2 reduced to 0.034.
Figures 4, 5, 6, 7, and 8 chart the effects on PC1 and
PC2 of factors having the greatest effects.
Milk EBV and Region. Clearly, PC1 score in-
creased with EBV milk yield classes (Figure 4). Fur-
thermore, the correlation between milk EBV and PC1
was +0.51, meaning that animals with the highest ge-
netic breeding value for total milk yield had the highest
Figure 3. Average lactation curves for all data (in black) and ac-
cording to cluster number (each point corresponds to the average milk level of daily milk yield. Milk EBV was not linked with
output of the animals). PC2, meaning that genetic selection on total-lactation

Journal of Dairy Science Vol. 101 No. 12, 2018

 LACTATION CURVES IN FRENCH DAIRY GOATS 11047

milk production was not linked to lactation persistency.

The region factor showed a gradient of increasing pro-
duction from south < north-east < north-west France.
Bouloc (1991) observed differences in milk produc-
tion between regions: goats from Pays-de-la-Loire
(north-west in our study) were more numerous in high-
producing groups whereas goats from Midi-Pyrenees
(South in our study) were relatively more numerous in
lower-producing groups.
Gestation Stage, Dry-Period Length, and
Breed. Goats had more persistent lactations when their
gestation was less advanced at 300 DIM and when they
were dried for a short period lower than 90 d before
their last kidding. Saanen goats were more productive
and more persistent than Alpines (Figure 5).
Like the current study, Bouloc (1991) also observed
that Saanen goats were more numerous in the groups of
persistent lactations. Furthermore, Rupp et al. (2011)
observed in the French goat population that Saanens
were more productive than Alpines. Like the current
study, Knight and Wilde (1988) reported that gestation
causes a decrease in milk yield compared with nonges-

Figure 5. Least squares means on the first 2 principal components

(PC1 and PC2) of the environmental factors strongly linked to the
second principal component: gestation stage (d) at 300 DIM (A: 1–60,
B: 60–90, C: 90–120, D: 120–150, E: 150–180), dry-period length (d)
(A: 1–60, B: 60–90, C: 90–120), and breed (Saanen or Alpine).

tating does. Although their study counted a limited

number of goats, they found that the influence of gesta-
tion was significant after 8 wk of gestation. Salama et al.
(2005), on goats milked once a day, also found that the
gestation caused a decrease of milk yield on gestating
goats compared with goats in extended lactation. In the
dairy cow, Druet et al. (2003) also showed a decrease of
milk yield according to the progress of gestation stage.
The effect of dry-period length on lactation curve has
been understudied. Knight and Wilde (1988) specified
that an absence of dry period cause a higher peak and
lower persistency. Caja et al.(2006) studied the effect
of dry-period length on 17 goats and found that no-
dry-period goats produced less than dry-period goats.
Atashi et al. (2013) studied the influence of the dry
period length on lactation curve shape, in dairy cows.
As in our study, they noticed the best persistency, the
smallest peak, and the higher milk lactation production
for the cow dried previously during a short time. In
Figure 4. Least squares means on the first 2 principal components contrary, they observed also that animals previously
(PC1 and PC2) of the environmental factors most strongly linked to dried for a long period produced higher milk yield at
the first principal component: milk EBV (kg) from the French genetic the peak than animals dried during an average period.
evaluation model (A: −349 to −79, B: −79 to −33, C: −33 to 13, D:
13 to 301), and home region of the herd (N-West = north-west; S = Similarly, Elvira et al. (2013b) showed that the longer
south; N-East = north-east). dry periods were not associated with high milk yield

Journal of Dairy Science Vol. 101 No. 12, 2018

11048 ARNAL ET AL.

6). From first to fourth parity, the goats produced more

milk on average but their lactations were less persistent.
Particularly, this highlights the important difference in
lactation production potential between goats in their
first parity and others parities. After the fourth par-
ity, level of milk production decreased more and more
whereas lactations were a little less persistent than in
the fourth lactation.
León et al. (2012) analyzed the production of 38,000
Murciano-Granadina breed goats following their first 5
lactations. Like in the current study, they showed that
persistency was greater in first lactation than second
lactation, which was in turn greater than third lacta-
tion and more.
In the French goat population, Boichard et al. (1989)
showed that age at kidding had a positive near-linear
relationship with milk yield. Bouloc (1991) found the
same effect of age at kidding on shape of the lactation
curve as the current study: goats kidding before 12 mo
had a higher lactation persistency. French dairy cattle
showed a similar pattern (Leclerc et al., 2008).
Kidding Month. Lactations were less persistent in
goats kidding during the spring but more persistent in
goats kidding in autumn (Figure 7). Note that the end
Figure 6. Least squares means on the first 2 principal compo-
nents (PC1 and PC2) of the environmental factors linked to the first
2 principal components: parity (1 to 8) and age at kidding (A: 9–10
mo at first kidding, 16–22 mo at second kidding or 22–34 mo at third
kidding; B: 11–13 mo at first kidding, 23–25 mo at second kidding or
35–37 mo at third kidding, C: 14–30 mo at first kidding, 26–56 mo at
second kidding or 38–77 mo at third kidding).

productions throughout the lactation, on their sheep

data.
SCS EBV. Goats had a more persistent lactation
when they had a high SCS EBV (less SCC). The rela-
tionship between high SCC index (indicator of geneti-
cally low SCC levels) and persistency in goats remains
unknown. One hypothesis could be that goats with
genetically low SCC levels are more resistant to mam-
mary infections and thus better able to maintain high
milk production. Indeed, Leitner et al. (2004) reported
that infected goats produced less milk than the nonin-
fected goats.
In bovines, Appuhamy et al. (2009) showed a negative
genetic correlation (−0.24) between persistency and a
mastitis occurring in the second half of lactation. Cole
and Null (2009) found a negative correlation between
SCS persistency and milk persistency. Note that Nayeri
et al. (2017) found a common QTL between persistency
and mastitis.
Age at Kidding and Parity. Goats that kidded Figure 7. Least squares means on the first 2 principal compo-
youngest (9–10 mo at first kidding, 16–22 mo at second nents (PC1 and PC2) of the environmental factor the most linked to
the second principal component: kidding month (A = January; B =
kidding, or 22–34 mo at third kidding) had a lower February; C = March; D = April; E = May; F = June; G = July; H =
milk production levels but higher persistency (Figure August; I = September; J = October; K = November; L = December).

Journal of Dairy Science Vol. 101 No. 12, 2018

 LACTATION CURVES IN FRENCH DAIRY GOATS 11049

where 84% of the goats did not pasture (Caillat et al.,

2016). Therefore, the spring grass did not explain the
persistency or the effect on PC3. However, abundant
literature is available on the effect of duration of the
day, which has a significant effect on goat lactation
curves (Chemineau et al., 2007; Russo et al., 2013).
These studies showed that long exposure to artificial
light (16 h per day) has a big effect on milk yield
throughout the lactation stage (up to +33% in France;
Delouis and Mirman, 1984). Here, we found this same
effect of photoperiod, independent of other environ-
mental factors, on a large data set of goats raised in
commercial herds.
A kidding in January seems interesting (high lacta-
tion yield, average persistency, and high production
around DIM120) but has to be put in relation to the
milk price, which fluctuates a lot seasonally, and to
justify a kidding out of season.
Herd Effect. The influence of herd was studied
on the residuals of model 1. The remaining variance
explained by the herd effect was relatively significant,
as illustrated by the R2 of 0.397, 0.178, and 0.162 for
residuals of MSPC1, MSPC2, and MSPC3, respectively. The
dispersion of herd effects on residuals was large, ranging
Figure 8. Least squares means on the third principal component from −243.8 to 314.5 for PC1, from −90.43 to 101 for
(PC3) of kidding month, which is the environmental factor most linked
to the PC3.
PC2, and from −35.82 to 43.21 for PC3. Bouloc (1991)
demonstrated that herd effect had the biggest effect,
in comparison to the other factors, on milk production
and persistency.
of the more persistent lactation occurred when days When the herd effect was replaced by the estimated
were long (at the beginning of summer), whereas the herd effect on total-lactation milk yield from the official
end of the less persistent lactation occurred when days genetic evaluation (4 classes), R2 was 0.329 for PC1 and
were short (at the beginning of winter). 0 for PC2 and PC3. This showed that the herd effect
The effect of kidding month on PC3 showed broad from the official genetic evaluation well explained the
variability, with a maximal value in February and a average of all-herd milk production level. However, as
minimum value in August (Figure 8). Principal compo- expected, the herd effect did not take into account the
nent 3 represented the production at approximately 120 effect of the herd on persistency and production at 120
DIM, and variations of the least squares means followed DIM.
the evolution of the photoperiod duration with 120 d Lactation length was not studied as a factor influenc-
difference. For example, a goat that kids in February ing lactation curves because of its links with breeder
will be at its 120 DIM in June when photoperiod is decisions (according to the potential of the animal, the
maximal, whereas a goat that kids in August will be at day of insemination, and breeder strategies). León et al.
its 120 DIM in December when photoperiod is minimal. (2012) specified that number of kids affected the pro-
The goat producing more during the longer days duction and shape of the lactation curve. A goat with
will yield more milk from the middle to the end of more kids had more milk in the beginning of the lacta-
the lactation according to kidding month. The influ- tion and was less persistent. Zumbach et al. (2008), like
ence of kidding month on persistency had already been Fernández et al. (2002), used the “number of born kids”
highlighted in studies by Bouloc (1991), Montaldo et al. factor in their genetic evaluation model. This factor
(1997), and León et al. (2012). In these 3 studies, the was not tested here as the information was unavailable.
greater persistency in goats kidding between October Likewise, we did not have more precise data on the
and January was explained by the thriving spring grass, farm regarding, for example, feeding system, system of
which occurs after the peak of lactation and maintains sales, breeding goal, or climatic conditions, to better
milk production. Here, 80% of the data were recorded understand the large influence of herd on the different
in the north-west (Loire-Atlantique, Poitou-Charente) SPC values.
Journal of Dairy Science Vol. 101 No. 12, 2018
11050 ARNAL ET AL.

CONCLUSIONS Biernacki, C., G. Celeux, and G. Govaert. 2003. Choosing starting

values for the EM algorithm for getting the highest likelihood in
multivariate Gaussian mixture models. Comput. Stat. Data Anal.
This study, based on a big recent data set, typolo- 41:561–575.
gized lactation curve shapes in the current French goat Biernacki, C., G. Celeux, G. Govaert, and F. Langrognet. 2006. Model-
population. It demonstrated a very real diversity of based cluster and discriminant analysis with the MIXMOD soft-
ware. Comput. Stat. Data Anal. 51:587–600.
lactation curve shapes in France: 30% of the lactations Boichard, D., N. Bouloc, G. Ricordeau, A. Piacere, and F. Barillet.
with a shape very different from the mean. This great 1989. Genetic parameters for first lactation dairy traits in the Al-
diversity of curves is mainly due to the global level of pine and Saanen goat breeds. Genet. Sel. Evol. 21:205.
Bouloc, N. 1991. Analyse de la forme de la courbe de lactation. Ap-
milk production and persistency of the lactation, but plication à l’étude des modalités d’allègement du contrôle laitier et
not only this. Indeed, we showed that the curve bends de prévision précoce de la production dans l’espèce caprine. PhD
in mid-lactation, at 120 DIM, which causes particular thesis, Institut national agronomique (INA) Paris Grignon, Paris,
France.
lactation curve shapes tied to duration of the photo- Brito, L. F., F. G. Silva, H. R. Oliveira, N. O. Souza, G. C. Caetano,
period. Furthermore, this trait could be viewed as a E. V. Costa, G. R. O. Menezes, A. L. P. Melo, M. T. Rodrigues,
complementary indicator of milk persistency. Lactation and R. A. Torres. 2017. Modelling lactation curves of dairy goats
by fitting random regression models using Legendre polynomials
persistency plays a key role in the discrimination of or B-splines. Can. J. Anim. Sci. 98:73–83.
lactation curve types, highlighting a certain variability. Caillat, H., B. Ranger, E. Bruneteau, C. Paraud, and R. Delargarde.
However, none of the environmental and genetic fac- 2016. Use of grazing in a dairy goat farm to design sustainable
production systems in France. Page 794 in Book of Abstract of
tors tested in our study had a strong effect on this the 67th Annual Meeting of the European Federation of Animal
trait. Some well-known factors could not be tested here, 67. Annual Meeting of the European Association for Animal Pro-
which may explain these results. This trait also has a duction (EAAP). Wageningen Academic Publishers, Wageningen
(Pays-Bas), Belfast, GBR.
genetic variability that we were unable to factor here. Caja, G., A. A. K. Salama, and X. Such. 2006. Omitting the dry-off
Other approaches will have to be tested to investigate period negatively affects colostrum and milk yield in dairy goats.
the genetic variability of the lactation curves. The use J. Dairy Sci. 89:4220–4228.
Carta, A., S. Casu, M. G. Usai, and S. Salaris. 2014. Heritability of
of random regression models taking into account en- persistency traits and their genetic correlations with milk yield
vironmental and genetic effects that vary throughout and udder morphology in dairy sheep. Page 39th ICAR Session,
lactation will have to be explored. The present study, Berlin (Germany). https:​/​/​www​.icar​.org/​wp​-content/​uploads/​
2016/​07/​Berlin​-2014​-Carta​-Casu​-Usai​-Heritability​-of​-persistency​
showing the effect of some environmental factors on the -traits​-and​-their​-genetic​-correlations​.pdf.
lactation curve shape, will help us to model the envi- Celeux, G., and G. Govaert. 1995. Gaussian parsimonious clustering
ronmental part of the random regression model. Once models. Pattern Recognit. 28:781–793.
Chemineau, P., B. Malpaux, J. P. Brillard, and A. Fostier. 2007. Sea-
the genetic evaluation is carried out, the correlations sonality of reproduction and production in farm fishes, birds and
between lactation shape curve and current or future se- mammals. Animal 1:419–432.
lected characters (on milk quality, longevity, or health) Cole, J. B., and D. J. Null. 2009. Genetic evaluation of lactation per-
sistency for five breeds of dairy cattle. J. Dairy Sci. 92:2248–2258.
will have to be evaluated to ensure the relevance of the Dai, X., P. Z. Hadjipantelis, H. Ji, H.-G. Mueller, and J.-L. Wang.
selection. The effects of the environmental factors could 2016. Fdapace: Functional Data Analysis and Empirical Dynam-
also be used in management advice in the dairy goat ics. R Package Version 0.2.5. Accessed Jan. 18, 2017. https:​/​/​cran​
.r​-project​.org/​web/​packages/​fdapace/​index​.html.
husbandry. Delouis, C., and B. Mirman. 1984. Influence de la durée quotidienne
d’éclairement sur la production laitière de la chèvre. Compte-Ren-
du 9e Journ. Rech. Ovine Caprine 352–360.
ACKNOWLEDGMENTS Druet, T., F. Jaffrézic, D. Boichard, and V. Ducrocq. 2003. Model-
ing lactation curves and estimation of genetic parameters for first
This project receives financial support from Min- lactation test-day records of French Holstein cows. J. Dairy Sci.
istry of Agriculture (Compte d’affectation Spécial au 86:2480–2490.
Elvira, L., F. Hernandez, P. Cuesta, S. Cano, J.-V. Gonzalez-Martin,
Développement Agricole et Rural, CASDAR, Rustic). and S. Astiz. 2013a. Accurate mathematical models to describe the
The first author received financial support from APIS- lactation curve of Lacaune dairy sheep under intensive manage-
GENE (Paris, France) and Association nationale de la ment. Animal 7:1044–1052.
Elvira, L., F. Hernandez, P. Cuesta, S. Cano, J.-V. Gonzalez-Martin,
recherche et de la technologie (ANRT, Paris, France). and S. Astiz. 2013b. Factors affecting the lactation curves of in-
tensively managed sheep based on a clustering approach. J. Dairy
Res. 80:439–447.
REFERENCES Fernández, C., A. Sánchez, and C. Garcés. 2002. Modeling the lacta-
tion curve for test-day milk yield in Murciano-Granadina goats.
Appuhamy, J. A., B. G. Cassell, and J. B. Cole. 2009. Phenotypic and
Small Rumin. Res. 46:29–41.
genetic relationships of common health disorders with milk and fat
Gengler, N. 1996. Persistency of lactation yields: A review. Interbull
yield persistencies from producer-recorded health data and test-
Bull. 12:87–96.
day yields. J. Dairy Sci. 92:1785–1795.
Gipson, T. A., and M. Grossman. 1990. Lactation curves in dairy
Atashi, H., M. J. Zamiri, and M. Dadpasand. 2013. Association be-
goats: A review. Small Rumin. Res. 3:383–396.
tween dry period length and lactation performance, lactation
Grossman, M., S. M. Hartz, and W. J. Koops. 1999. Persistency of
curve, calf birth weight, and dystocia in Holstein dairy cows in
lactation yield: A novel approach. J. Dairy Sci. 82:2192–2197.
Iran. J. Dairy Sci. 96:3632–3638.

Journal of Dairy Science Vol. 101 No. 12, 2018

 LACTATION CURVES IN FRENCH DAIRY GOATS 11051
ICAR. 2017. ICAR guidelines. Accessed Mar. 8, 2018. https:​/​/​www​ Nayeri, S., M. Sargolzaei, M. K. Abo-Ismail, S. Miller, F. Schenkel, S.
.icar​.org/​index​.php/​icar​-recording​-guidelines/​. S. Moore, and P. Stothard. 2017. Genome-wide association study
Kearney, J. 2010. Food consumption trends and drivers. Philos. Trans. for lactation persistency, female fertility, longevity, and lifetime
R. Soc. B Biol. Sci. 365:2793–2807. profit index traits in Holstein dairy cattle. J. Dairy Sci. 100:1246–
Knight, C. H., and C. J. Wilde. 1988. Milk production in concurrently 1258.
pregnant and lactating goats mated out of season. J. Dairy Res. Pollott, G. E. 2000. A biological approach to lactation curve analysis
55:487–493. for milk yield. J. Dairy Sci. 83:2448–2458.
Langrognet, F., R. Lebret, C. Poli, S. Iovleff, B. Auder, C. Biernacki, Pulina, G., M. J. Milán, M. P. Lavín, A. Theodoridis, E. Morin,
G. Celeux, and G. Govaert. 2016. Rmixmod: Supervised, Unsu- J. Capote, D. L. Thomas, A. H. D. Francesconi, and G. Caja.
pervised, Semi-Supervised Classification with MIXture MODel- 2018. Invited review: Current production trends, farm structures,
ling (Interface of MIXMOD Software). R Package Version 2.1.1. and economics of the dairy sheep and goat sectors. J. Dairy Sci.
Accessed Jan. 18, 2017. https:​/​/​cran​.r​-project​.org/​web/​packages/​ 101:6715–6729.
Rmixmod/​index​.html. Rupp, R., V. Clément, A. Piacere, C. Robert-Granié, and E. Manfredi.
Lebret, R., S. Iovleff, F. Langrognet, C. Biernacki, G. Celeux, and G. 2011. Genetic parameters for milk somatic cell score and relation-
Govaert. 2015. Rmixmod: The r package of the model-based un- ship with production and udder type traits in dairy Alpine and
supervised, supervised and semi-supervised classification mixmod Saanen primiparous goats. J. Dairy Sci. 94:3629–3634.
library. J. Stat. Softw. 67:241–270. Russo, V. M., A. W. N. Cameron, F. R. Dunshea, A. J. Tilbrook,
Leclerc, H., D. Duclos, A. Barbat, T. Druet, and V. Ducrocq. 2008. and B. J. Leury. 2013. Artificially extending photoperiod improves
Environmental effects on lactation curves included in a test-day milk yield in dairy goats and is most effective in late lactation.
model genetic evaluation. Animal 2:344–353. Small Rumin. Res. 113:179–186.
Leitner, G., U. Merin, and N. Silanikove. 2004. Changes in milk com- Salama, A. A., G. Caja, X. Such, R. Casals, and E. Albanell. 2005.
position as affected by subclinical mastitis in goats. J. Dairy Sci. Effect of pregnancy and extended lactation on milk production in
87:1719–1726. dairy goats milked once daily. J. Dairy Sci. 88:3894–3904.
Lenth, R. V. 2016. Least-Squares Means: The R Package lsmeans. J. Sölkner, J., and W. Fuchs. 1987. A comparison of different measures
Stat. Softw. 69:1–33. of persistency with special respect to variation of test-day milk
León, J. M., N. P. P. Macciotta, L. T. Gama, C. Barba, and J. V. Del- yields. Livest. Prod. Sci. 16:305–319.
gado. 2012. Characterization of the lactation curve in Murciano- van der Werf, J. H. J., M. E. Goddard, and K. Meyer. 1998. The use
Granadina dairy goats. Small Rumin. Res. 107:76–84. of covariance functions and random regressions for genetic evalu-
Macciotta, N. P. P., D. Vicario, and A. Cappio-Borlino. 2006. Use of ation of milk production based on test day records. J. Dairy Sci.
multivariate analysis to extract latent variables related to level of 81:3300–3308.
production and lactation persistency in dairy cattle. J. Dairy Sci. Van Hekken, D. L., M. H. Tunick, and Y. W. Park. 2005. Effect of
89:3188–3194. frozen storage on the proteolytic and rheological properties of soft
Menéndez-Buxadera, A., A. Molina, F. Arrebola, M. J. Gil, and J. caprine milk cheese. J. Dairy Sci. 88:1966–1972. https:​/​/​doi​.org/​
M. Serradilla. 2010. Random regression analysis of milk yield and 10​.3168/​jds​.S0022​-0302(05)72872​-1.
milk composition in the first and second lactations of Murciano- Yao, F., H.-G. Müller, and J.-L. Wang. 2005. Functional data analysis
Granadina goats. J. Dairy Sci. 93:2718–2726. for sparse longitudinal data. J. Am. Stat. Assoc. 100:577–590.
Montaldo, H., A. Almanza, and A. Juárez. 1997. Genetic group, age Zumbach, B., S. Tsuruta, I. Misztal, and K. J. Peters. 2008. Use of a
and season effects on lactation curve shape in goats. Small Rumin. test day model for dairy goat milk yield across lactations in Ger-
Res. 24:195–202. many. J. Anim. Breed. Genet. 125:160–167.
Mucha, S., R. Mrode, M. Coffey, and J. Conington. 2014. Estimation
of genetic parameters for milk yield across lactations in mixed-
breed dairy goats. J. Dairy Sci. 97:2455–2461.

Journal of Dairy Science Vol. 101 No. 12, 2018