Thailand is located in the tropics between latitudes 5°37′ N to 20°27′ N and longitudes 97°22′ E to 105°37′ E. The total area is about 513,000 km². The weather conditions in Thailand are divided into three seasons, namely winter (mid-October to mid-February), summer (mid-February to mid-May) and rainy season (mid-May to mid-October) with an average temperature and relative humidity of 26°C and 74%, 29°C and 71%, and 28°C and 76%, respectively. The annual rainfall in most areas of the country is from 1,200 mm to 1,600 mm [7]. The three Thai government dairy cattle research farms where the data were obtained from were located in two main regions of Thailand (north and north-eastern). One farm in the north is located in Chiang Mai province (latitude 18°35′41.9″ N and longitude 98°51′53.1″ E) and the other two farms in the north-east are located in Nakhon Ratchasima (latitude 14°40′09.9″ N and longitude 101°26′46.7″ E) and Sakon Nakhon province (latitude 17°09′52.9″ N and longitude 104°02′07.8″ E) as shown in Figure 1.

The data from the three Thai government dairy cattle farms mentioned above were provided by the Bureau of Animal Husbandry and Genetic Improvement (BAHGI), Department of Livestock Development (DLD) in Thailand. The two farms located in Chiang Mai province and Nakhon Ratchasima province raised both Holstein-Friesian (HF) and crossbred HF while the third farm located in Sakon Nakhon province raised only crossbred HF. Upgrading local or Zebu breed with HF semen or natural mating with HF bulls has been used to improve productivity and maintain tropical insect and disease tolerance in these herds.

All animals in the three farms were raised under the same guidelines given by BAHGI. Nevertheless, some aspects were different such as feeding and health management because of the differences of locations, weather conditions, feed resources, farm machinery and disease prevalence. Fresh ruzi (Brachiaria ruziziensis) and Napier grass (Pennisetum purpureum) were fed during the rainy season while in the dry season, ruzi silage and hay were fed. In addition, fresh corn (Zea mays) and corn silage were used in Chiang Mai farm and Nakhon Ratchasima farm in some particular years (nutrient component of roughage and concentrate feed using in dairy cattle farm is described in Pangmao et al [2].

The data comprised of 25,789 test-day milk yield and milk composition records of 1,468 cattle from Lactation 1 to 3 of HF and crossbred HF dairy cattle that calved between 1990 and 2015. The total number of sires, dams and individual animals in the pedigree file were 287, 1,237 and 4,753 respectively. The data consist of cow number, birth date, calving date, drying date, parity, lactation length, monthly or test-day records of milk yield, fat and protein percentage (FP and PP). The records of cows with less than three test-day milk yield records were excluded. The animals were classified into five breed groups (BG) based on the percentage of HF blood (HFB) as 1 (HFB≤75), 2 (75<HFB≤87.5), 3 (87.5<HFB≤93.75), 4 (93.75<HFB≤99.99) and 5 (HFB = 100). Calving months were grouped into three seasons, namely winter (November to February), summer (March to June) and rainy (July to October).

The Wood model [4] specifies the expected milk yield on day t of lactation, and has the following form, the shape depending on the three parameters:

where W(t) is the theoretical or expected milk yield at day t, k is a scalar factor and equal to log_ea, b is related to the rate of increase prior to the peak and c is related to the rate of decrease after the peak. The expected milk yield for cow-lactation i on days-in-milk t can be fitted as the following nonlinear mixed model:

where ɛ_ij is the random effect associated with records. Random effects for each cow’s lactation (i) are indicated as deviations (K_i, B_i, C_i) from the fixed effect means:

and are assumed to have a multivariate normal distribution as follows:

The parameters (k, b, c) were chosen over (a, b, c), as similar to the analyses by Hall [8], based on the preliminary fixed effects model, the former showed a better approximation to a multivariate normal distribution than did the latter (data not shown), a requirement for the nonlinear mixed model specification.

Cumulative milk yield to day T (day 305) was then obtained as MILKW=∫0TW(t)dt. Fitting of the Wood model as a nonlinear mixed model was conducted using the nlme package in R [9] and calculation of cumulative milk yield through use of the pgamma function in R. The parameters derived from Wood model were used to calculate lactation curve characteristics namely peak milk yield (PMILK = e^k(b/c)^be^−b), days to peak milk yield (DPMILK = b/c) and persistency (S = c^−(b+1)) [4] (PSMILK = log_e(S) = log_e(c^−(b+1)) = −(b+1)log_ec.

The 305-day milk yield (MILKTI) calculated by using the TIM [10] with adjustment for the last term as follows:

where M₁, M₂, …, M_n are the weights (kg), given to one decimal place, of the milk yield in the 24 hours of the recording day, I₁, I₂, …, I_n−1 are the intervals, in days, between recording dates, I₀ is the interval, in days, between the lactation period start date and the first recording date and I_n is the interval, in days, between the last recording date and the end of the lactation period. The FP and PP were averaged over the lactation and used for subsequent analysis.

To obtain estimates of genetic parameters, univariate linear mixed models were fitted to each of the seven traits (MILKW, MILKTI, PMILK, DPMILK, PSMILK, FP, and PP), with a separate model for each of the three lactations. The fixed effects of each model comprised herd (Herd), BG, year and season of calving (YSOC); and the random effects were animal linked to a pedigree, and residual error. The form of each model was

where y_ijkl is an observation of the trait on animal l; μ is the overall mean; Herd_i is the fixed effect of herd (level, 1 to 3); YSOC_j is the fixed effect of year and season of calving (level, 1 to 78); BG_k is the fixed effect of BG (level, 1–5); Animl is the random animal effect; and e_ijkl is the random residual effect. Heritability was estimated using REML estimates for h2=σA2/(σA2+σe2), where σA2 and σe2 are the additive genetic and residual variances. Estimated breeding values for each trait were also produced from the fitted univariate models.

Following this a repeatability model was fitted to all three lactations simultaneously, via the following model:

where all terms are as defined previously, with the addition of Lact_m, the fixed effect of lactation m (m = 1, 2, 3); and Pe_l, the random ‘permanent environment’ term for animal l. For this model, heritability was estimated as h2=σA2/(σA2+σpe2+σe2) where σpe2 is the permanent environment variance, and repeatability as r=(σA2+σpe2)/(σA2+σpe2+σe2).

The set of univariable models at each lactation was ex tended to a series of pairwise bivariate models, for each pair of traits using the same terms as the univariate model, to allow estimation of genetic and phenotypic correlations between traits in the same lactation and in the same trait between two different lactations. Model fitting of these univariate and bivariate models was conducted using ASReml-R [11].

The estimated additive genetic variance, residual variance and heritability for milk yield, lactation characteristic traits and milk composition of lactations 1 to 3 are shown in Table 2. In general, the heritability estimates of all traits from first parity cows were higher than in parity 2 and 3 except for FP that were lower (0.08, 0.25, and 0.07, respectively). This was due to a combination of greater genetic variances and smaller residual variances of traits in the first lactation. The heritability estimates of MILKTI in the first and third parity were similar (0.19). The heritability estimates of MILKTI and MILKW were similar (0.21 and 0.19) in lactation 1 while the heritability of MILKW was lower compared to MILKTI in lactation 2 and 3 (0.01 vs 0.12 and 0.08 vs 0.19, respectively). Lactation curve trait heritability estimates from cows in lactation 1 were higher than in lactation 2 and 3 whereas heritability of DPMILK and PSMILK were very low in lactation 2 and were zero or no detectable additive genetic variability in lactation 3. The heritability estimates of FP in lactation 1 and 3 was lower (0.08 and 0.07) than in lactation 2 although heritability of PP in lactation 1 was higher (0.45) compared to lactation 2 and 3 (both 0.22). The overall heritability using the repeatability model across three lactations was low to medium, ranging from 0.14 to 0.31 for all the traits while very low for DPMILK and PSMILK (0.03 and 0.04, respectively) as shown in Table 3. The overall repeatability for lactations 1 to 3 was low for most of the traits at between 0.22 and 0.34 (Table 3). MILKTI has the highest repeatability (0.45) while DPMILK and PSMILK has lowest repeatability (0.03 and 0.04, respectively).

The genetic correlation estimates between milk yield and lactation curve traits were moderate to high (0.57 to 0.99) except for the genetic correlation estimates between PMILK and DPMILK, and between PMILK and PSMILK (0.32 and 0.33) as shown in Table 4. The genetic correlations of milk composition with other traits were low and negative. PP has negative genetic correlations with other traits (−0.46 to −0.06) except for the genetic correlation with FP (0.34). The phenotypic correlations between MILKW, MILKTI and PMILK were in the high range from 0.70 to 0.92 while between DPMILK and PSMILK with other traits were negative and low (−0.07 to 0.25) except between DPMILK and PSMILK (0.93). The phenotypic correlations of FP and PP with other traits were negative and low (−0.11 to 0.11) except between FP and PP which was moderate (0.29).

The genetic correlation estimates between milk yield and lactation curve traits were high to very high (0.64 to 0.99) except for the genetic correlations between DPMILK and PSMILK with MILKTI (0.21 and 0.32 respectively) and between PSMILK and PMILK (0.23) as shown in Table 5. The genetic correlation estimates between FP and other traits were low to moderate (0.14 to 0.29) while between PP and other traits were negative (−0.46 to −0.06) except between PP and FP (0.34). The phenotypic correlation between MILKW with MILKTI, PMILK, PSMILK and PP were high ranging from 0.71 to 0.93 except between MILKW with DPMILK and FP were low (0.16 and 0.09, respectively) while the phenotypic correlation between MILKTI with other lactation curve characteristic traits were negative to moderate (−0.06 to 0.28). The phenotypic correlation between FP and PP with other traits were negative to low (−0.10 to 0.09) except for the phenotypic correlation between PP with FP and MILKW (0.28 and 0.93, respectively).

Table 6 shows the estimated genetic and phenotypic correlations of milk yield, lactation curve traits and milk composition of the cows in lactation 3. The genetic correlation between MILKW and MILKTI is moderate (0.45). The genetic correlation of PMILK, DPMILK, and PSMILK with MILKW and MILKTI are high (0.68 to 0.99) but moderate between MILKW and DPMILK (0.46) and highly negative genetic correlation between MILKW and PSMILK (−0.99). Most of the genetic correlation estimates of FP with other traits are moderate to high (0.44 to 0.99) except with DPMILK (−0.09) while mostly genetic correlations of PP with other traits are negative (−0.82 to −0.14) but with PMILK is 0.25. The phenotypic correlation between MILKW and MILKTI is 0.71. The phenotypic correlation estimates between MILKW and MILKTI with DPMILK and PSMILK are low to moderate ranging from 0.02 to 0.28 but between MILKW and MILKTI with PMILK are high (0.70 and 0.78). Phenotypic correlations between PMILK, DPMILK and PSMILK are negative low (−0.06 and −0.10) but between DPMILK and PSMILK is high (0.93). Phenotypic correlations of FP and PP with other traits are low (−0.04 to 0.08) except between FP and PP (0.25).

The genetic correlation estimates of milk yield, lactation characteristics and milk composition between lactations 1 to 3 are shown in Table 7. Most of the genetic correlation estimates of the traits are high (0.75 to 0.99). However, the genetic correlations of DPMILK, PSMILK, and FP could not be estimated between lactations 1 and 3 and lactations 2 and 3. The phenotypic correlations for milk yield, lactation characteristics and milk composition of Thai dairy cattle between lactations 1, 2, and 3 are shown in Table 8. The phenotypic correlation of all traits is low for most of the traits.

Boxplots of cow EBVs over year of birth of all traits in lactation 1 are shown in Figures 2 to 8. None of the traits showed a systematic improvement over this period and did not show consistent trends. The EBV trend for most traits improved a little from 1982 to 2000, after that the EBV trend fluctuated and decreased during the later years. Boxplots of sire EBVs over the year of birth of all traits in lactation 1 are in the Supplementary document as Figures S1–S7. The sire EBVs of all traits showed the same inconsistent trend over years as the cow, however, there was an improvement in the trend in the last 10 years compared to 10 years before for most traits. The average sire EBV of MILKW in the 10 years before 2002 was −42.7 kg compared to +100.8 kg in the last 10 years after 2002.

The heritability of MILKW was highest in lactation 1 compared to lactation 2 and 3 (0.21 vs 0.01 and 0.08, respectively) while for MILKTI, heritability is similar between lactation 1 and 3 but lower in lactation 2 (0.19, 0.12, and 0.19, respectively). The low heritability of MILKW in lactation 2 was due to a very low additive genetic variance estimate with a relatively larger standard error, which may have been due to data limitations (small sample size) rather than reflecting a true biological result, which would be consistent with the results of Hammami et al [15] who reported similar heritability estimates across the first three lactations. The heritability estimates of MILKW and MILKTI in lactation 1 is comparable to those reported by Mohammed et al [16] (0.24) and Boonkum and Duangjinda [17] (0.23). However, the heritability of MILKW and MILKTI in lactation 1 in this herd is lower than other studies (0.35: König et al [18], 0.43: Seangjun et al [19], 0.34: Sarakul et al [20]). The lower milk yield heritability in this study compared with other studies may be due to a small additive genetic variance and/or a high residual variance, suggesting this trait was highly affected by the environmental factors such as farm and feed management, hot and humid tropical environment. In addition, the different size of the data set and the models used for analysis also may have an effect on estimation of variance components and hence heritability. Nonetheless, milk yield calculated from both methods show the possibility of improving by a selection program. The heritability estimates of lactation curve traits (PMILK, DPMILK, and PSMILK) were low to moderate within the range from 0.00 to 0.23 for lactation 1 to 3. The heritability of PMILK and DPMILK in lactation 1 was similar to that reported by Chegini et al [21] (0.23 vs 0.26 and 0.10 vs 0.10) although PSMILK was higher (0.14 vs 0.05). In general, the heritability of lactation curve traits for all lactations was quite low except the moderate heritability of PMILK in lactation 1 and 3 (0.23 and 0.17, respectively) which means PMILK can be improved by selective breeding while DPMILK and PSMILK trait need an improvement of environmental management. The heritability of fat percentage in lactations 1 and 3 were lower (0.08 and 0.07) than lactation 2 (0.24). The heritability of FP in lactation 1 was lower than most other studies e.g. Welper and Freeman [22] (0.51), Boujenane [23] (0.39), Harris and Pryce [24] (0.48), Kim et al [25] (0.41) and Koonawootrittriron et al [26] (0.22). The heritability of PP in lactation 1 was higher than lactation 2 and 3 (0.45 vs 0.22 for both) and was in the range reported in earlier studies [22] (0.55), [24] (0.52) and [25] (0.43).

The genetic correlation between MILKW and MILKTI was high in lactations 1 and 2 (0.98 and 0.99) but medium in lactation 3 (0.45). The genetic correlation of MILKW with PMILK, DPMILK, and PSMILK for cows in lactation 1 (0.94, 0.57, and 0.57) were comparable with genetic correlations of 305-day milk yield with peak yield, days in milk at peak yield and persistency reported by Chegini et al [21] (0.97, 0.52, and 0.44, respectively) which also used the Wood model for calculation of milk yield and lactation curve traits in Iranian Holstein cows. However, the genetic correlation of MILKTI with lactation curve traits for cows in lactation 1 was high, ranging from 0.70 to 0.86. The high correlation of MILKW and MILKTI with PMILK in lactation 1 was similar to other studies [27–29]. Therefore, the selection of high peak milk yield would tend to improve milk production in these herds although DPMILK and PSMILK would not improve much due to low genetic correlation of both traits with PMILK (0.32 and 0.33). The genetic correlation of MILKW and MILKTI with FP in lactation 1 was low and positive (0.24 and 0.23) and with PP was medium and negative (−0.40 and −0.46) while many studies reported the negative correlation between milk yield with FP and PP [22,24,30]. Selection for high milk production might slightly decrease protein percentage.

The genetic correlations of the same trait between lactations (Table 7) were high for all the traits, ranging from 0.75 to 0.99, although the correlations for DPMILK, PSMILK and FP between lactation 1 and 3, and lactation 2 and 3, could not be estimated as the models did not converge. The high genetic correlation estimates of all traits between lactations suggested that the selection of animals for first lactation curve traits in the herd will improve traits in second and third lactations as well, although the phenotypic correlations between lactations for most of the traits were low and negative for PS.

The genetic trend of sire and cow EBVs over year of birth for all traits showed inconsistent trends in these three herds. In cows there was some decline in EBVs in recent years. However, the sire EBVs showed a slight improvement in the last 10 years after 2002. The use of these bulls may have led to the genetic gain of cows after 2012. However, this will require the analysis of data beyond 2012 to assess the genetic improvement in the cows. Overall, the results on genetic trend suggested that the selection and use of breeding bulls has not been effective. This is in spite of the fact that the selection of bulls has been made based on EBVs. There could be a number of factors which might have impacted the expected genetic gains in these herds. The limitation of financial support, inconsistent animal management and feeding practices might be the issues causing the problems in these herds. In addition, the limitation of selection, semen used, an unintentional culling of high producing cow due to health problems, and an ineffective breeding plan might also impact the genetic progress in these herds. To improve the breeding program, these issues should be examined while establishing a long-term plan by the Department of Livestock Development (DLD) or Ministry of Agriculture and Cooperatives (MOAC) to improve the production status in these herds, in conjunction with the research farm management. Based on the moderate heritability estimates of milk production traits, the selective breeding is expected to result in effective genetic progress in these herds.

In Thailand, most of the data are collected from farms under the support of the DLD and the Dairy Farming Promotion Organization of Thailand and are analysed separately by these two main organizations, although some data analysis is from other sources such as private and educational farms. These organizations publish annual sire summaries and which are mostly used by the contributing farm members. However, to help speed up and improve milk production in the country, Thailand should establish a central organization for managing and analysing the data of dairy cattle across Thailand. In addition, data on a wider variety of traits should be collected and analysed in the form of selection indices which farmers can consider and choose based on their own farm needs.