Genetic selection in dairy cattle breeding programs is nowadays increasingly relying on test-day (TD) models in many countries using TD records instead of aggregated 305-d yield records [1]. Advantages of the utilization of TD models compared to 305-d lactation yields models are: i) The TD record is not dependent on the accuracy of 305-d extension methods; ii) The ability of the TD to account for environmental effects that are specific to each TD; and iii) The shape of the lactation curve describing genetic variation during milking periods may be obtained from TD models within lactations as well as across lactations; and the breeding values can be estimated for any day of the lactation period [2,3]. A distinct economic advantage of the use of TD models in breeding programs, compared to only using records of complete lactations is to reduce the cost of recording dairy cattle performances [1]. Furthermore, TD models provide 4% to 8% more accurate genetic evaluations in dairy cattle over genetic evaluation from 305-d models [4]. Numerous studies have investigated the estimates of genetic parameters from TD milk traits using various random regression models (RRM) [5–7]. In the repeatability model, all TD records are considered to be a measure of the same trait [8,9]. In multiple trait models, lactation is considered as a separate trait [10,11]. Variance components in a random regression TD model were usually estimated by Legendre polynomials. However, different orders of Legendre polynomials for additive genetic and permanent environmental effects were used. Nevertheless, the choice of order can be based on practical considerations, accounting for technicalities such as the complexity of the used models, number of parameters, and the time required running the resolution for parameter estimates. The disadvantage of RRM is related to computing requirements to store all TD and time of computing compared with 305-d yields [6]. Pool et al [8] reported that for using a TD approach on a national genetic evaluation scale, the size of equation system can exceed hundreds of millions of equations. However, computer capacity is continually improving, and the ability to use RRM might be resolved [10]. Holstein cows (around 200,000 cows) account for more than a half of the Tunisian cow population. The adoption of a random regression TD model as the official genetic evaluation for production traits in Tunisia, as implemented by most Interbull countries, may be an alternative method to enhance the genetic improvement of the Holstein population. The first step in implementing a routine genetic evaluation using a TD model is to estimate necessary genetic parameters [6]. The objective of this study was to estimate genetic parameters of milk, fat, and protein yields within and across lactations in Tunisian Holsteins using a random regression TD model.

Data used in this study were obtained from the official milk recording maintained by the Center for Genetic Improvement of the Livestock and Pasture Office. The original data were TD of milk, fat, and protein yields from the first 3 lactations of Tunisian Holsteins collected from 1997 to 2014 and included a total of 960,000 records. Data were edited on the following requirements: Day in milk ranged from 5 to 305 days. Reasonable ages of calving in a specific lactation were applied. The latter were those of cows that had their tests between ages of 20 and 48, 34 and 65, and between 45 and 85 months for the first, second, and third calving, respectively. Lactations were required to have at least 5 TD records for milk, fat, and protein yields. Herds with fewer than 5 cows per herd-year of calving were discarded. Only records from the first 3 parities that had data for all production traits on a given TD were kept. Within cow, if lactation p (p = 2 to 3) was present, lactations p-1 and p-2 when applicable were also present. There were 90,327 animals in the pedigree file. Animals with missing pedigree information were discarded. The edited data structure and the mean for milk, fat, and protein yields and other variables are shown in Table 1. The model for trait r (milk, fat, and protein) in lactation p (first, second, and third) was

Where, y_ijmprt are records on trait r of cow m in lactation p on days in milk (DIM) t, CG_i is the effect of contemporary group defined as herd test date recording; AS_j is the fixed effect of the jth subclass of age at calving-season of calving (nested within parity). Four seasons were defined (September to November, December to February, March to May, and June to July); z_tq = is the qth Legendre polynomial corresponding to day t of lactation; a_mpq = random additive genetic coefficients of cow m corresponding to polynomial q of parity p; pe_mpq = random permanent environmental coefficients of cow m corresponding to polynomial q of parity p; e_ijmprt was the residual effect for each observation. The choice of the third-order polynomials used in this study was inferred from the study by Hammami et al [12] who reported that the constant, linear, and quadratic polynomials were highly related to the first three eigenvalues, and explained more than 95% of the variance components for all three random effects.

In matrix notation, the model can be written as

Where, y = the vector of TD record for milk, fat, and protein yields for the first three lactations; b = vector of fixed effects; a = vector of random regression coefficients for genetic additive animal effects; p = vector of random regression coefficients for permanent environmental effects; e = vector of residual effects. X, Z, and W = incidence matrices that relate observations to their respective effects. The phenotypic covariance matrix V of the observations is given by.

Where G, and P are the random regression (co)variances matrices for the genetic and permanent environmental effects, respectively. I = identity matrix and ⊗ is the Kronecker product. A is the additive genetic covariance matrix among all animals, and R is the diagonal matrix of the residual variance. Residual variance was assumed to be constant within DIM intervals. Variance components for additive genetic and permanent environmental random regression were performed by the Bayesian method using the GIBBS3F90 program [13]. Posterior means of the parameters of interest were calculated using 100,000 samples after discarding the 20,000 iterations as a burn-in period. Convergence of Gibbs chains was determined based on inspection of plots of realizations of selected parameters. The genetic variance matrix for all DIM was obtained as G* = ΨGΨ′ where G* is a 301 by 301 genetic (co)variance matrix for all DIMs ranging from 5 to 305 days. Ψ is a 301 by 36 matrix with the values of the 12 coefficients of the 3 order Legendre polynomial for each DIM from 5 to 305 days. The same operation was applied to P. Genetic parameters for complete 305-d lactations for all traits can be computed by using G, and P (co)variance matrices and 305-d vectors of Legendre polynomials Ψ₃₀₅. Vector of Legendre polynomials was obtained by summing up the four coefficients from 5 to 305 d for genetic and permanent environmental effects. As an example, the genetic additive variance for 305 d milk, fat, and protein yield was defined as: Ψ305GΨ305′. The same process was applied to permanent environmental variances for 305-d milk yields. Heritability for day i(hi2) was calculated as:

Where, σa(i)2 and σp(i)2 are the additive genetic and permanent variances for day i. Heritability for 305-day yield h(305d)2 was calculated as:

Where, σa(305)2, σp(305)2 are the additive genetic and permanent environmental variances for 305 day for each trait. Genetic correlations between traits i and j were calculated as:

Where G_{(i, j)} are the genetic (co)variance matrix for i and j.

The estimated additive genetic, permanent environmental variances, and heritabilities for TD milk, fat, and protein yields in lactations 1, 2, and 3 are given in Tables 2 to 4. The genetic variances and permanent environmental variances were high at the beginning and end of lactation and lower in between. Additive genetic variances of daily milk, fat, and protein yields tended to be more stable over DIM in the 3 parities compared with permanent environmental variances. For all traits in all lactations, the permanent environmental variance estimates were consistently higher than the genetic additive variance. On the other hand, the permanent environmental variances for milk, fat, and protein yields were more heterogeneous across lactations associated with high difference between the peripheries and the middle part of lactation. For TD milk, fat, and protein yields, the pattern of the additive genetic variances was relatively similar in all lactations, decreasing from the beginning of lactation until around DIM 185 and increased afterward. Heritability of milk yield was lowest in mid-lactation and higher in early and late lactation, especially in lactations 2 and 3. Daily heritability of milk, fat, and protein yields tended to increase as parity increases. Heritabilities for daily milk, fat, and protein yields were higher in the third lactation, and ranged from 0.095 to 0.287, 0.069 to 0.258, and 0.087 to 0.285, respectively. Heritability trends for daily milk, fat, and proteins yields were similar to trends for additive genetic variance but were less extreme at the start and end of the curve because of higher permanent environmental variances. Estimates of heritabilities and genetic correlations of 305-d yields and their standard errors are given in Table 5. Heritability estimates ranged from 0.14 in the first lactation to 0.2 in the third lactation. Estimates of heritability for milk were higher than for fat and protein yields. Genetic correlations among lactations for yield traits were relatively high and ranged from 0.78±0.01 for milk to 0.82±0.03 for protein, between the first and second parities, respective estimates were from 0.73±0.03 to 0.80±0.04 between the first and third parities, and from 0.82±0.02 to 0.84±0.04 between the second and third parities. Genetic correlations between milk, fat, and protein yields within lactation were moderate to high, ranging from 0.74±0.02 to 0.8±0.06 between milk and fat yield, from 0.8±0.03 to 0.85±0.06 between fat and protein yields, and from 0.83±0.02 to 0.88±0.05 between milk and protein yields.

A concave shaped trend for the additive genetic and permanent environmental variances; that is, greatest values at the peripheries of the lactation, was reported by Naserkheil et al [3], Pool et al [8], and Pereira et al [14]. However, the shapes obtained in this study were different from those reported by Druet et al [6] and Mayeres et al [15] and both using Legendre Polynomials. In fact, Druet et al [6] observed that the additive genetic variance was highest in mid-lactation and lowest at the beginning and the end of lactation. The genetic additive variance was small in comparison with the permanent environmental variance for milk, fat, and protein yields, this is in agreement with results obtained by Caccamo et al [11], and Druet et al [16]. It was also observed that the genetic and permanent environmental variances were smaller in the first lactation compared with the second and third lactation. Similar pattern was obtained by Miglior et al [17] using a random regression TD model in Italian and Chinese Holsteins. Permanent environmental variances were larger than additive genetic variances (25.85 vs 4.9 kg2 around 5 DIM for milk yield in first parity) indicating that environmental effects had a higher impact on the variation of milk yields among cows than genetic effects. Pool et al [8] noticed that the permanent environmental effect needs to be modeled by higher orders of Legendre polynomials than the additive genetic effect. Hammami et al [12] using fourth-order Legendre Polynomial reported that permanent environmental variances in the Tunisian cow population were 65% larger than that in Luxembourgian population in a joint analysis of both cow populations. Limited feed resources and harsh environmental factors are more likely to occur in Tunisian climatic conditions; where cows are managed under conditions completely different from those where their sires were selected, might be a possible explanation for these high estimates of the permanent variances [18]. Considerable variation in the pattern of the variance components have been reported in the literature probably due to differences in genotypes, climatic conditions, herd management levels among countries and production systems. Several researchers have considered each country as a distinct environment such as Weigel et al [19] and Ojango et al [20], whereas other authors argued that some herds in different countries could be considered to be in the same environment [21]. Gebreyohannes et al [22] reported that the genetic parameter estimates with random regression TD models are influenced by covariance structure of additive genetic, residual, and permanent effects, and especially by the regression functions. Misztal et al [23] argued that the results based on random regression were very heterogeneous. Heritability of daily milk, fat, and protein yields were low and ranged from 0.1 to 0.28 for all parts of lactations. Estimates of heritability for milk were higher than for fat and protein yields. This result is in agreement with results of Strabel and Jamrozik [24]. Several studies reported that high values of heritability at the peripheries of the lactation can be explained by the difficulties in modeling the corresponding variances because of the biological process generated at the onset and termination of the lactation [22]. Other studies advanced that the unreasonably high heritability estimates at the beginning and the end of lactation could be explained by a lack of information to model the variability, particularly nearing the end of lactation [14,24]. The use of herd by year of calving as a random effect improves modeling the trajectory of variance components and heritability curves [24]. The range of daily heritabilities in this study was from 0.1 to 0.28. Current estimates are lower than those reported by De Ross et al [5], Druet et al [6], Pereira et al [14], and Miglior et al [17] who found heritabilities higher than 0.3 for milk, fat, and protein yields in some parts of the lactation trajectory. Heritabilities for 305-d milk, fat, and protein yields from the current study were comparable with those obtained by Ben Gara et al [25]. Hammami et al [12] also reported a low heritability (around 0.16) of 305-d yield for first lactation in Tunisian Holsteins using test-day random regression sire model. Likewise, heritability estimates were in the same range as found by Ojango et al [20] in Kenyan Holsteins. In contrast, results from the current study on heritability estimates are markedly smaller than those found in other previously published results on Holsteins using multiple-trait multiple lactation models [5,7,17]. Several authors explained low heritability estimates in some occasions by limited production levels [18,20,25]. Ojango et al [20] in the same study found heritability values (0.26 vs 0.45 in the Kenyan vs UK Holstein populations, respectively) associated with two 305-day milk yield (4,557 vs 7,674 kg in the Kenyan vs UK Holstein populations, respectively). Nevertheless, other authors attributed the low heritability values to the environmental factors especially to heat stress conditions [12,26,27]. Hammami et al [27] argued that animals with the high genetic merit in cold environments do not necessarily have the high genetic merit in warm stressful environments. In Tunisia, the climate varies from arid in the South to humid in the North, and characterized by hot summers coupled with high humidity [18]. The feeding system is unbalanced and rations are based mainly on concentrates. The forage is characterized by poor quality, and high rate in indigestible cellulosic constituents that could be possible causes of the lower milk yield. All these factors could lead to decreases in production performances and increase in the incidence of health troubles such as acidosis at the herd level. Within-parity genetic correlations obtained in the present study were relatively high among milk, fat, and protein yields (>0.73) especially in the third parity. Genetic correlations between milk and protein yields were the highest suggesting the possibility to introduce protein yield in breeding objectives without any pressure on milk yield. Estimates of genetic correlations between second and third lactations were larger than between the first and second lactations, which confirm the results obtained by Muir et al [7]. However, Hammami et al [28], using a comparable model with smaller Tunisian data and including RR coefficients of the herd by year of calving, reported genetic correlations estimates of 0.64 to 0.86 between the first and second lactation, and of 0.6 to 0.81 between the second and third lactation. Reasons for differences in estimates are not apparent, but can be explained by the evolution of the Tunisian management system in recent years. Genetic parameters for milk, fat, and protein yields for Tunisian Holstein dairy cows were estimated using a random regression TD model. Results are similar to previous reports where a comparable model to that used in the current study was applied. Genetic parameter estimates suggest that the adoption of a random regression TD model as the official genetic evaluation for production traits in Tunisia is feasible.