Genetic parameters of milk and lactation curve traits of dairy cattle from research farms in Thailand

Objective This study was aimed to estimate the genetic parameters, including genetic and phenotypic correlations, of milk yield, lactation curve traits and milk composition of Thai dairy cattle from three government research farms. Methods The data of 25,789 test-day milk yield and milk composition records of 1,468 cattle from lactation 1 to 3 of Holstein Friesian (HF) and crossbred HF dairy cattle calved between 1990 and 2015 from three government research farms in Thailand were analysed. 305-day milk yield was estimated by the Wood model and a test interval method. The Wood model was used for estimating cumulative 305-day milk yield, peak milk yield, days to peak milk yield and persistency. Genetic parameters were estimated using linear mixed models with herd, breed group, year and season of calving as fixed effects, and animals linked to a pedigree as random effects, together with a residual error. Univariate models were used to estimate variance components, heritability, estimated breeding values (EBVs) and repeatability of each trait, while pairwise bivariate models were used to estimate covariance components and correlations between traits in the same lactation and in the same trait across lactations. Results The heritability of 305-day milk yield, peak milk yield and protein percentage have moderate to high estimates ranging from 0.19 to 0.45 while days to peak milk yield, persistency and fat percentage have low heritability ranging from 0.08 to 0.14 in lactation 1 cows. Further, heritability of most traits considered was higher in lactation 1 compared with lactations 2 and 3. For cows in lactation 1, high genetic correlations were found between 305-day milk yield and peak milk yield (0.86±0.07) and days to peak milk yield and persistency (0.99±0.02) while estimates of genetic correlations between the remaining traits were imprecise due to the high standard errors. The genetic correlations within the traits across lactation were high. There was no consistent trend of EBVs for most traits in the first lactation over the study period. Conclusion Both the Wood model and test interval method can be used for milk yield estimates in these herds. However, the Wood model has advantages over the test interval method as it can be fitted using fewer test-day records and the estimated model parameters can be used to derive estimates of other lactation curve parameters. Milk yield, peak milk yield and protein percentage can be improved by a selection and mating program while days to peak milk yield, persistency and fat percentage can be improved by including into a selection index.


INTRODUCTION
Milk production in Thailand in 2018 was 1.25 million tons [1], most of which is used for consumption within the country. However, the current milk production is not enough to meet the domestic demand. Hence, milk production is the main focus in the dairy cattle breeding improvement programs in Thailand. Crossbreeding local Zebu cattle with exotic Bos taurus cattle, mainly Holstein Friesian (HF), is the main approach used to improve milk production in the country. Artificial insemination with HF semen for upgrading has been the main method of herd improvement. The dairy system is generally driven by smallscale farmers and cooperative organizations which have support from the Thai government [2]. In developing countries where there are often less data being recorded, genetic evaluation using the computationally simpler 305-day lactation yield record in an animal model is still commonly used.
Many methods have been used for calculating 305-day milk yield. The test interval method (TIM) is one of the standard methods approved by the International Committee for Animal Recording [3]. Milk production estimation by TIM is calculated as the area under the lactation curve up to a specific day, e.g. day 305, using a simple trapezoidal numerical integration method. A typical lactation curve involves a rapid increase after calving until it reaches the peak early in the lactation, followed by a decline in milk yield until the end of milk production. The lactation curve can be an important management tool as a summary of milk production in dairy cattle. Understanding the lactation curve and factors that influence its shape can help farmers to manage their herd effectively in terms of feed allocation, reproductive and health management which can be optimised to increase the milk yield. A number of mathematical models have been used to describe lactation curves. The Wood model [4] is one of the models that has been widely used due to the relatively few parameters to be estimated (three) and variety of model shapes [5]. It is more robust compared to other models when fitting lactation curve to irregular and infrequent test-day sampling regimes [6]. The inclusion of suitable and accurate estimates of genetic parameters of lactation curve characteristics into a breeding improvement program can be used to improve milk production. This study used data from three dairy cattle research centres managed under the supervision of the Department of Livestock Development of Thailand. The aim of these centres is to study various aspects in dairy cattle raising in Thailand such as nutrition, farm management and breed improvement in order to gain and extend knowledge to smallholder farmers. The objective of this study was to estimate the genetic parameters, including heritabilities, genetic and phenotypic correlations, genetic trends of estimated breeding values (EBVs) of milk yield, lactation curve traits and milk composition of Thai dairy cattle in these three government research farms.

Geography and climate conditions
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 2 . The weather conditions in Thai-land 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.

Description of data and herd management
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].

Data and statistical analysis
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 e a, 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 cowlactation 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: 150 ki = k+Ki, bi = b+Bi, ci = c+Ci 151 152 and are assumed to have a multivariate normal distribution as follows:

156
The parameters (k, b, c) were chosen over (a, b, c), as similar to the analy 157 the preliminary fixed effects model, the former showed a better approximation 158 distribution than did the latter (data not shown), a requirement for the 159 specification.

160
Cumulative milk yield to day T (day 305) was then obtained as

166
The 305-day milk yield (MILKTI) calculated by using the test interval met 167 for the last term as follows:

168
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 150 ki = k+Ki, bi = b+Bi, ci = c+Ci 151 152 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 157 the preliminary fixed effects model, the former showed a better approximation to a multivariate normal 158 distribution than did the latter (data not shown), a requirement for the nonlinear mixed model 159 specification.

160
Cumulative milk yield to day T (day 305) was then obtained as  The 305-day milk yield (MILKTI) calculated by using the test interval method [10] with adjustment 167 for the last term as follows:
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) b e -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 e c.
The 305-day milk yield (MILKTI) calculated by using the TIM [10] with adjustment for the last term as follows:    where M 1 , M 2 , …, M n are the weights (kg), given to one decimal place, of the milk yield in the 24 hours of the recording day, I 1 , I 2 , …, I n-1 are the intervals, in days, between recording dates, I 0 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 l l; μ is the overall mean; Herdi is the fixed effect of r and season of calving (level, 1-78); BGk is the fixed ndom animal effect; and eijkl is the random residual timates for 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 where is the permanent environment variance, and repeatability as . The set of univariable models at each lactation was extended 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].

RESULTS
The descriptive summary of milk yield, lactation characteristic traits and milk composition of lactations 1 to 3 is shown in Table 1. Milk yield from both models (MILKW and MILKTI) was highest in lactation 2 (3,685 and 3,682 kg) while lowest in lactation 1 (3,348 and 3,379 kg). Peak milk yield (PMILK) was lowest for lactation 1 at 14.95 kg while for lactation 2 and 3 were 17.52 kg and 17.49 kg, respectively. The cows in lactation 1 took longer to attain peak yield (DPMILK) than the cows in lactation 2 or 3 (49 days vs 39 days and 41 days), however, the persistency (PSMILK) of lactation 1 cows was higher than cattle in lactation 2 and 3 (6.68 vs 6.43 and 6.43). For milk composition, fat percentage (FP) of the cows in lactation 1 was less than lactations 2 and 3 (3.55% vs 3.64% and 3.64%) while protein percentage (PP) was the same for all lactation (3.04%).

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

Genetic and phenotypic correlations among traits in lactation 1
The genetic correlation estimates between milk yield and

Genetic and phenotypic correlations among traits in lactation 2
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

Genetic and phenotypic correlation of each trait across lactations
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.

DISCUSSION
This study reported the genetic parameters of milk yield, lactation curve traits and milk composition from three government dairy farms estimated using data from the first three lactations of cows. For first lactation cows, milk yield and peak milk yield have moderate heritability, protein percentage has high heritability which can be improved by a selection scheme while days to peak milk yield, persistency and fat percentage have low heritability. The similarity of estimated milk yield and high genetic correlation between both methods suggested that either method can be used for a genetic improvement program in these herds. However, the Wood model has advantages over the test-interval method by requiring fewer test-day records, particularly if strategically selected [12]. In addition, lactation curve traits  such as persistency can be estimated easily by the Wood model. However, even with the Wood model the estimates of heritability of milk yield and lactation curve traits in lactations 2 and 3 had high standard errors compared to the heritability estimates. Further study including more data from lactations 2 and 3 would be required to obtain more precise estimates. The positive genetic correlations between cumulative 305-day milk yield and peak milk yield, and between days to peak milk yield and persistency within a lactation implies that selection to improve peak milk yield would also improve milk yield, and that selection for improved days to peak milk yield would improve persistency as well. In addition, the high genetic correlations of all traits between lactations suggested that selection of favourable animals in the first lactation would also improve the corresponding trait in second and third lactations, and with these correlations being so high, may indicate they are essentially the same trait. Both cumulative 305-day milk yields (MILKW and MILKTI) were similar within all the three lactations. The  lowest mean MILKW was for lactation 1 compared with lactations 2 and 3 (3,348 kg vs 3,685 kg and 3,564 kg, respectively) and this ranking is in agreement with Hossein-Zadeh [13] (9,186 kg vs 10,386 kg and 10,000 kg, respectively) which may be due to partition of nutrition for growth and milk production in lactation 1 cows. The highest PS in lactation 1 was in agreement with Gengler [14] who reported higher persistency in first lactation than the other lactations. In the present study, FP increased over the three lactations while PP was steady across all three lactations.

Heritability
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 re-  sults 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).

Genetic and phenotypic correlations
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][28][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.

Genetic trend
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 longterm 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.

CONFLICT OF INTEREST
We certify that there is no conflict of interest with any financial organization regarding the material discussed in the manuscript.

FUNDING
The authors received no financial support for this article.

ACKNOWLEDGMENTS
The data for this work were provided by the Department of Livestock Development, Thailand. We would like to sincerely thank the Royal Thai Government Scholarships for SP's PhD studies.