The TD records of Thai dairy cows from the first three lactations, with calving between 1993 and 2013, were obtained from the dairy database of the Bureau of Biotechnology in Livestock Production (BBLP), DLD. To obtain data sets with a consistent size, the following records were included: age at calving was restricted to 18 to 48, 30 to 60, and 41 to 75 months for the first, second and third lactations respectively; an interval of between 5 and 35 days from parturition to the first TD; daily milk yields between 2 and 40 kg; at least 5 TD records per lactation; a minimum of 150 days in milk (DIM). Furthermore, all three traits (milk, fat, protein yield) were required for each TD. All third lactation cows were required to have first and second lactation records. Likewise, the second lactation cows were required to have first lactation records. Additionally, the sires for the cows in the dataset were identified.

After applying these criteria, a total of 168,996, 63,388, and 27,145 TD records (milk, fat and protein yield) of the first, second and third lactations, respectively, measured in different calendar months within herds from 29,230 lactations, which were daughters of 1,116 sires belonging to 124 herds, were left to be analyzed. Ancestors of cows in the final data set were traced back in pedigree as far as the parents were known. Furthermore, the numbers of TD records for milk, fat and protein yields were not equal because fat and protein yields were missing in some TD due to technical reasons.

The breed groups were classified into three groups using the percentage of Holstein Friesian (HF) blood level: <87.5% HF, 87.50% HF to 93.75% HF and >93.75% HF. Because dairy cattle in Thailand used up-grading system to increase milk yield and expected <87.5% HF group provide low milk yield, 87.50% HF to 93.75% HF group provide medium milk yield and >93.75% HF group provide high milk yield. The months of milk testing were divided into 12 months, as follows: January to December. Calving seasons were classified into three groups: winter (November to February), summer (March to June) and rainy season (July to October). Details of the final data file are given in Table 1.

Data were analyzed with a 3-trait, 3-lactation RR-TDM. The matrix notation of the model is

where y = a vector of milk, fat and protein yields; b = a vector of the fixed effects: herd-year-month of testing, breed-season of calving-month in tested milk group, and linear and quadratic age at calving; a = a vector of RR coefficients for animal genetic (AG) effect; p = a vector of RR coefficients for permanent environmental (PE) effect; e = a vector of residual (R) effects; Q = a matrix of five modified Legendre polynomials (constant, linear, quadratic, cubic, and quartic), as defined by Gengler et al [6]; and X, Z, and W = incidence matrices relating observations to various effects. The residual effects were taken to be independently distributed and with constant variance along DIM. Modeling of the herd-year-month of testing effect as a fixed effect across lactations was found to be advantageous in populations with small herd sizes [10]. Druet et al [3] also found that the use of fixed classes assured the best fit compared with parametric curves (Legendre polynomials, Ali-Schaeffer curve and Wilmink curve). Due to their large number of parameters, fixed classes allow for more flexibility than parametric curves. In addition, any record in the parametric curve will influence the whole curve. In contrast, the influence of the data is local in the case of fixed classes. Also, classes of month of tested milk can be cross-classified with other effects with the potential to influence lactation shapes. The covariance structure of the model is:

where A = the numerator relationship matrix; ⊗ = a Kronecker product function; G₀ and P₀ = 45×45 for (co)variance matrices of RR coefficients for AG and PE effects, respectively; I = an identity matrix (the number of cows with records); and R = 9×9 for diagonal matrix of residual variance matrices corresponding to each trait with elements on lactation n assumed to be constant throughout the lactation.

Figure 1A to 1C shows trajectories of the milk, fat and protein yields by DIM for each lactation of Thai dairy cattle. The curves were based on a 15-day moving average. The average peak curves for milk yields of the first, second and third lactations were 43, 43, and 28 DIM, with milk yields of 15.84± 4.25, 17.56±4.85, and 18.14±4.96 kg, respectively. There was a trend for the peak DIM decreasing and the peak yield increasing as the lactation increased, which is similar to the lactation curves that were previously found in the literature, but lower than those presented by Miglior et al [12]. Differences from the previous studies could be the breed (purebred vs crossbred) and population differences. The first and second lactation curves took longer to reach the peak yield than the third lactation, indicating that it is more persistent than the others. Curves for fat and protein yield had the same patterns as those for milk yield.

The estimation of the daily AG, PE, and R variances in the first three lactations for milk, fat and protein yields are shown in Table 2. Generally, all variances increased with lactation for all yield traits. Similar trends were reported in other studies [4,5,12–16]. The AG variances in the second and third lactations were consistently higher than the PE and R variances for all traits, unlike in the first lactation. The difference between AG variances in the second and third lactations was small for all productions. However, its difference increased greatly from the first to second lactation. The difference in AG variances between primiparous and multiparous cows indicated that the expression of their genetic potential was higher in multiparous than primiparous cows.

The pattern of AG variances by DIM across lactation for milk, fat and protein yields are shown in Figure 2. The AG variance for milk yield had the same pattern in all parities, which was high at DIM 80 after parturition, and then decreased exponentially. On the other hand, AG variance for fat and protein yields had a different trend as milk yield, being large at the beginning, small in the middle, and moderate at the end of lactation. Figure 3 shows the pattern of PE variances for all traits. The PE variances for milk yield in all lactations were highest at the beginning of lactation and subsequently decreased until the end of lactation, whereas PE variance for fat and protein yield were relatively high at the periphery of lactation.

The trends in the AG and PE variance estimates for all yields throughout lactation obtained in this study were comparable to the trends found by Miglior et al [12] and de Roos et al [17], in which RR-TDM was also applied and Legendre polynomials were used to describe the random curves of Chinese Holsteins and Dutch dairy cattle, respectively. However, the AG variance patterns in our study are different from those of Hammami et al [15], who reported higher estimates of AG variances for milk yield in the beginning and end than in the middle of lactation. Our patterns also differ from those of Muir et al [4] and Zavadilova et al [13], who found that AG and PE variances increased with progressive lactations. The decrease in AG variability after the onset of the first lactation was the same as that from a study of Holstein cows [18]. Misztal et al [19] reported that the level and pattern of daily milk yield variances obtained by random regression models were heterogeneous. Pool et al [20] reported that the shape of the variance curves across lactation could be modeled with sufficient accuracy by using a third-order polynomial for the genetic part, but a fourth-order Legendre polynomial was needed for the PE part. López-Romero and Carabaño [21] also reported that smaller order polynomials could be more suitable for AG than for PE. The differences in variance shape through lactations were caused by using various orders of Legendre polynomials, which described random curves for the AG and PE effects. Besides this, the R variance was not assumed to be constant during lactation, as in our study. A constant R variance was assumed, which might somehow influence the estimation of the PE variance [22], but Olori et al [23] found that keeping a constant R variance produced a bias in the residual term in early lactation, without significantly affecting the estimates of the other variance components.

For higher estimates of AG and PE variance in early lactation, we can discuss with postpartum condition. The transition period and the early lactation period are characterized by the mobilization of fat body reserves to cover the energy demands and include various metabolic and endocrine adaptations. Recent research has suggested considerable variation among high selected dairy cows for the release of energy fuels from stored adipose tissues [24]. Tamminga et al [25] observed great differences in fat mobilization among cows during the first 8 weeks of lactation, ranging from 8 to 57 kg of body fat. Some cows can cope with metabolic stress very well, whereas others do not [26]. This ability to deal with metabolic disturbances apparently also depends on important individual cow factors and results from genetic variance. Postpartum management resulted in higher PE variance. The TD genetic variances for the production yield of multiparous cows were also higher than those of primiparous cows. In general, the increase of genetic variations among lactations could reflect differences in the genetic ability of cows to produce milk and the milk composition. We speculate that animals in later lactation expressed their genetic ability differently. One possible reason might be due to metabolic differences, inducing a differential nutrient partitioning between adult and first calving cows [27]. Moreover, primiparous cows need to divert nutrients for growth in order to achieve a mature size.

The average daily and 305-d yield heritability estimates for production yield traits in the first three lactations are shown in Table 3. Estimates of the average daily and 305-d heritability for yield traits ranged from 0.33 to 0.48 and 0.46 to 0.64, respectively. They likely increased with lactation. Average daily heritabilities for milk yield (0.36 to 0.48) were generally higher than those for fat yield (0.33 to 0.44) and close to those for protein yield (0.37 to 0.48). These values were significantly higher than those presented in Polish black and whites [14], Italian Holsteins [4], Tunisian Holsteins [15], Chinese Holsteins [12], and Australian Holsteins [16]. Heritability estimates of 305-d milk and fat yields were slightly lower than protein yield. The estimates of all production yields in the first lactation were similar to the results of de Roos et al [17] and Konstantinov et al [16]. de Roos et al [17] found that heritability estimates for all yields in the first three lactations, using a random regression model with Legendre polynomials, were 0.53, 0.53, and 0.54 for milk yield, 0.51, 0.50, and 0.54 for fat yield, and 0.46, 0.45, and 0.47 for protein yield, respectively. Konstantinov et al [16] found that the values of 305-d heritabilities in the first lactation were 0.44, 0.48, and 0.39 for milk, fat and protein yield, respectively.

In general, the level and pattern of milk yield heritabilities obtained from RR models are sensitive to the model applied; Misztal et al [19] and other recent studies have confirmed this fact. Nevertheless, large estimates of AG variances and heritabilities are associated with high milk production levels [3, 4,6,17]. Low AG and heritability estimates have been reported for populations with low to medium production levels [6,14]

Heritabilities of TD production yields along the lactation trajectory for selected DIM in the first three lactations are also illustrated in Figure 4. The heritability estimates ranged from 0.29 to 0.44, 0.38 to 0.55, and 0.34 to 0.56 for milk yields; from 0.30 to 0.46, 0.38 to 0.51, and 0.37 to 0.52 for fat yields; and from 0.34 to 0.48, 0.44 to 0.53, and 0.44 to 0.56 for protein yields in the first, second and third lactations, respectively. The heritabilities in multiparous cows were higher than in primiparous cows, with a similar shaped curve for each production trait. However, daily milk yield heritability patterns were different from daily fat and protein yields.

Estimated heritabilities for daily milk yield were moderate at the early stage, increased to the highest heritability at the middle stages after parturition, and then gradually decreased until the late stage of lactation, whereas estimated heritabilities for fat and protein yield were relatively high at the periphery of lactation. The heritability curves for all production yield traits were similar to those for Dutch dairy cattle reported by de Roos et al [17]. However, these were different from other reports [4,12,15]. The differences in the heritability estimates for all production yields among the studies could be due to breed differences, the different types of models, and the effects included within the model. In addition, the feeding system and heat stress in Thailand can explain the opposite heritability shapes for production yields. Dairy farms in this study were distributed throughout the region of Thailand. Roughage sources for dairy cattle were by-products from agriculture, whole corn and grass, grass and rice straw in the north, central, west and east, respectively. A separate feeding diet was operated in most of Thailand. Under this feeding system, the farmer used different high-concentrate rations based on the quantity and quality of roughages that were readily available. Moreover, dairy cattle were kept under an open house with high heat-humidity weather for years. Therefore, dairy cattle tried to adapt to production under the environment of Thailand in terms of naturally-produced forage. We might speculate that these factors indicate that dairy cattle in tropical regions express their genetic potential differently.

The AG and PE correlations of 305-d production yields in the first three lactations are shown in Table 4. The estimated correlations varied depending on the yield traits, and there were all positive. In general, correlations between consecutive lactations were higher compared with those for lactations that were further apart. The AG correlations obtained between the yield traits in the first and second lactation (0.48 to 0.81) were the largest among all genetic correlations. The AG correlations between the second and third lactation ranged from 0.39 to 0.63, whereas correlations between the first and third lactation ranged from 0.41 to 0.72. Overall, genetic correlations (305-d yield) among production yields within lactations (0.44 to 0.69) were higher than those across lactations (0.36 to 0.68). The values ranged from 0.56 to 0.60 between milk and fat yields, from 0.58 to 0.69 between milk and protein yields, and from 0.44 to 0.52 between fat and protein yields. The PE correlations (305-d production yields) between yield traits within lactations (0.57 to 0.68) were also higher than the genetic correlation across lactations (0.03 to 0.19).

Results from our study were in accordance with those of Muir et al [4] and Hammami et al [15], who used a multiple-trait RR model in Italian Holsteins and Tunisian Holsteins, respectively, but with lower estimated values. However, a larger genetic correlation between milk and protein yield than between milk and fat yield was also reported by Jakobsen et al [2].

The AG and PE correlations of production yields on different DIM within the first three lactations are shown in Table 5. In general, genetic correlations of each production yield during the different stages of lactation were all positive. The correlations between yields on days that were close together were higher compared to those for days that were further apart. The AG correlation of milk yields over the first three lactations ranged from 0.46 to 0.55 for early, from 0.73 to 0.95 for mid- and from 0.52 to 0.86 for late lactation. The AG correlation of fat yield and protein yield followed the same trend as milk yield, ranging from 0.06 to 0.25 for early, from 0.47 to 0.89 for mid- and from 0.00 to 0.34 for late lactation. There were high correlations of production yields between 65 DIM and 245 DIM, ranging from 0.73 to 0.95, from 0.51 to 0.89 and from 0.45 to 0.89 for the milk, fat and protein yields, respectively. The high correlations of production yields imply that selection to improve production traits in the mid stage of lactation was more rapidly achieved than in the other stages of lactation.

Estimates of AG correlations between production yields on the same DIM in the first three lactations are shown in Figure 5. For all traits, the largest genetic correlations occurred between the first and second lactation, and the lowest was observed between the second and third lactation. The shapes of correlations across DIM showed a similar pattern for all traits. Generally, they increased from early lactation, fluctuated in mid-lactation and decreased during late lactation. For milk and protein yield, the correlations between the same DIM in the consecutive lactations were below 0.94 at the beginning of lactation, between 0.70 and 0.93 in the middle of lactation, and below 0.74 at the end of lactation. However, for fat yields, the genetic correlation between the first and third lactations were moderate, but they were smaller, 0.70 over the whole trajectory of lactations.

Among various lactations, similar shapes of correlation curves at the same DIM were also reported by Strabel and Jamrozik [14] and Hammami et al [15]. Berry et al [18] found a genotype by environment interaction for body condition score (BCS), implying that genes that influence BCS may differ according to nutritional (i.e. concentrate feeding level, grazing severity, and silage quality) or environmental management. Our study found that genetic correlations between all yields were moderate in distant lactations. We might speculate that our study indicated that production yields would be under the influence of a similar gene set. However, gene sets were limited in expression under the natural diversity in Thailand in later lactations. Furthermore, this study found that genetic correlations between fat yields were low in both distant lactations and consecutive DIM. Hammami et al [15] reported that high temperature and the lack of quality forage had an influence on the decline of fat more strongly than that of milk and protein. We might speculate that these factors limited the expression of identical gene sets in each DIM and lactation.

It should be noted that, in this study, there were no negative estimates between any of the selected DIM. Negative genetic correlations for selected DIM were reported in several studies [5,13]. The problem of selecting the best random regression model and related co-variance components is not trivial and has been discussed in several studies [3,21]. On one hand, it comes from the fact that different countries use different functions to describe the random curves. Although Legendre polynomials have become a standard for this part of the model, there are differences in their order between different countries. For example, the fourth order is used in Canada [28] and the fifth order is used in the United Kingdom [29]. Overall, in this study, estimates of the genetic parameters were in good agreement with the literature values.

Changes in the average EBVs for different 305-day production yields of Thai dairy cattle in the first, second and third lactations, against birth year from 2001 to 2011, are presented in Figure 6; the corresponding linear random regression coefficients are given in Table 6. Generally, the annual genetic trends of production yield were positive and increased as the lactation number increased. Milk yield of the third lactation was 10.83 and 9.96 kg higher than that of the first and second lactations, respectively, whereas the fat and protein yields of all lactations increased slightly. These trends demonstrate the effectiveness of selection for the improvement of milk, fat and protein yields. These results were consistent with a previous report [30], in which genetic trends for Thai dairy cattle were positive for milk, fat and protein yields, whereas those of milk component (fat and protein percentage) were close to zero or negative. The negative or zero genetic trends for fat and protein percentage are likely the result of a major emphasis on milk yield, with farmers neglecting fat and protein percentage in sire selection at the farm level over the past years. This condition could cause a correlated response for fat and protein percentages as a result of the selection for milk because of the probable negative correlation between milk yield and milk components [31]. Consequently, it could cause a slightly positive genetic trend in fat and protein yields. It is, therefore, necessary that, apart from milk yield, fat and protein percentage are also included in the breeding goals in order to optimize the genetic improvement of Thai dairy cattle.