Genetic aspects of Wood's lactation curve parameters in Jersey crossbred cattle using Bayesian approach

Abstract

The study was undertaken to estimate the genetic parameters of lactation curve parameters of Wood's function in Jersey crossbred cattle using the Bayesian approach. Data on 33,906 fortnightly test day milk yields of 1,718 lactation records of Jersey crossbred cows, maintained at the ICAR-National Dairy Research Institute in West Bengal, were collected over a period of 40 years. The lactation curve parameters including ‘a’ (initial milk yield after calving), ‘b’ (ascending slope up to peak yield) and ‘c’ (descending slope after peak yield) and lactation curve traits, peak yield (y_max), time of peak yield (t_max) and persistency of milk yield (P) of individual cow for each lactation were estimated using the incomplete gamma function (Wood's model) by fitting the Gauss–Newton algorithm as an iteration method using PROC NLIN procedure of SAS 9.3. Variance components and genetic parameters of lactation curve parameters/traits were estimated by a repeatability animal model using the Bayesian approach. Estimates of heritabilities were found to be 0.18 ± 0.05, 0.09 ± 0.03 and 0.11 ± 0.04 for parameters ‘a’, ‘b’ and ‘c’, respectively and 0.24 ± 0.05, 0.12 ± 0.04, and 0.15 ± 0.05 for y_max, t_max and P, respectively. Repeatability estimates were 0.31 ± 0.03, 0.21 ± 0.04 and 0.30 ± 0.04 for parameters ‘a’, ‘b’ and ‘c’ respectively and 0.39 ± 0.03, 0.24 ± 0.03 and 0.37 ± 0.03 for y_max, t_max and p, respectively. Genetic correlations among lactation curve parameters/traits ranged from −0.75 to 0.95. Existence of genetic correlations among lactation curve parameters/traits indicated substantial genetic and physiological relationships among lactation curve parameters/traits of Jersey crossbred cattle.

Milk production is a complex physiological process which generally follows a definite trend of milk secretion throughout lactation in dairy animals (Farhangfar and Rowlinson, 2007). Lactation curve, the graphical representation of milk yield against time after calving, represents the biological efficiency of lactating animals in different phases of lactation (Scott et al., 1996). Appropriate knowledge of the shape of the lactation curve of any breed is required to maintain an optimum level of milk production in a dairy farm. Several researchers (Pérochon et al., 1996; Macciotta et al., 2005) postulated that the lactation curve provides valuable information on the features of milk production, which could be used in decision making process for herd management and breeding and ultimately helps to design effective selection strategies for evaluating the genetic potentiality of dairy animals. The biological aspects of the shape of the lactation curve can be determined by models using only three parameters (Rekaya et al., 2000) but to explore the biological efficiency of dairy cows, the incomplete gamma function or Wood's model (Wood, 1967) is well known for its accurate biological interpretation of different phases of lactation in dairy animal (Radjabalizadeh et al., 2022). This model also efficiently estimates the lactation curve traits such as peak milk yield, time of peak yield and persistency of milk yield from its estimated curve parameters (Lopez-Ordaz et al., 2009). Several authors have included the lactation curve traits such as peak milk yield and persistency as selection criteria to improve total milk yield (Ferris et al., 1985; Batra et al., 1987). Lactation curve parameters of Wood's model are subjected to genetic variation, suggesting that genetic information may be used to alter the lactation curve's shape (Varona et al., 1998). Since the shape of lactation curves differs between breeds and species of animals, it is important to study the genetic aspects of the curve parameters which are treated as different traits.

Precise estimations of (co)variance components and genetic parameters of any trait are important in animal breeding programs for more accurate prediction of breeding value, prediction of genetic gain and to achieve the optimum selection response for faster genetic improvement. Several researchers adopted Frequentist as well as Bayesian statistical methods for genetic parameter estimation and statistical inferences. However, Bayesian approaches can provide more precise estimates where there are problems of a finite sample size because any set of data, regardless of how big or small, this approach allows for the drawing of conclusions (Misztal, 2008; Aspilcueta-Borquis et al., 2010). Since the Bayesian approach overcomes the limitations of having a limited sample size by providing an accurate a posteriori distribution for any large or small data set from which interpretations can be inferred, so this statistical method was used in the present study. Scanty literatures are available on the genetic aspects of lactation curve functions in crossbred cattle. Therefore, the present investigation was undertaken to elucidate the genetic variations of Wood's lactation curve parameters and their associations with lactation curve traits in Jersey crossbred cattle using a multi-trait animal model employing the Bayesian method.

Material and methods

Experimental animals and location of study

Lactation data on Jersey crossbred cows, maintained at the Eastern Regional Station of ICAR-National Dairy Research Institute, Kalyani, Nadia, West Bengal, India, were collected and used for the present investigation. The farm is located at latitude of 22.59°N and longitude of 88.29°E with an altitude of 9.75 m above the mean sea level (MSL) in the lower Gangetic basin of West Bengal, India. The climate is sub-tropical and hot-humid with average annual rainfall of 1250 mm and a maximum precipitation (83%) during July to October. The temperature ranges from 13°C to 25°C during winter and 26°C to 37°C during summer with minimum and maximum humidity of 58% and 91% respectively. Hence, dairy animals maintained in this farm are exposed to a wide range of climatic changes round the year. The managemental practices of animals in this herd have been described by Koloi and Mandal (2020). Due to non-existence of multiple Jersey crossbred herds in India, data on lactation traits from multiple herds over this period are not available (Koloi and Mandal, 2020).

Data and statistical analyses

Data pertaining to test day milk yields and pedigree records of Jersey crossbred cows were collected from daily milk yield registers and pedigree registers of the farm, respectively. A total of 33 906 fortnightly test day milk yields spanning the 5^th to 305^th day of lactation from 1718 lactation records were collected over a period of 40 years. A total of 21 test-day milk yields in each lactation was used for the study. The animals with culling, abortion, stillbirth and disease records, abnormal lactations or unknown pedigree were removed from the data sets. Cows having less than 100 d of lactation length and less than 2 kg daily milk yield at the beginning of lactation were discarded. The outliers in the data sets (> mean ± 3SD) were removed to ensure normal distribution of data sets for all the traits before analysis and only sires having at least 3 progenies were considered for genetic analysis. Due to the small number of records in each calving year, the entire duration of study was divided into eight 5-year periods of calving. Based on prevalent climatic conditions of the area, each year of calving was divided into three seasons, winter (November–February), summer (March–June) and rainy (July–October). The parities of animals were classified as 7 groups.

Fitting lactation curve and estimation of lactation curve parameters

The lactation curve was described by incomplete gamma function or Wood's model (Wood, 1967). The lactation curve parameters of Wood's function were derived using the following nonlinear model.

$$y_t = at^be^{{-}ct}$$

where y_t is the daily milk yield at t ^th test day, and ‘a’ (initial milk yield after calving), ‘b’ (ascending slope up to peak yield) and ‘c’ (descending slope after peak yield) are the model parameters that described the biological property and shape of the lactation curves. Animals having negative estimates of ascending slope up to peak yield (‘b’) and/or descending slope after peak yield (‘c’) obtained through Wood's model failed to generate the lactation curve traits and were excluded from further analysis, being considered as having an atypical lactation curve.

The lactation curve parameters (‘a', ‘b' and ‘c’) of incomplete gamma function were estimated for each individual cow for each lactation based on the fortnightly test day milk using PROC NLIN procedure of SAS 9.3 (SAS Institute Inc., 2011) where the Gauss–Newton algorithm was used as an iteration method. From the estimated lactation curve parameters, the lactation curve traits peak milk yield (y _max), time of peak milk yield (t _max) and persistency of milk yield (P) for individual animal for each lactation were calculated as described by Wood (1976) (as detailed in the online Supplemental Table S1).

Estimation of variance components and genetic parameters

Lactation curve parameters/traits generated for each individual cow were subjected to genetic analysis. In the preliminary analysis, the significant fixed effects for each trait were identified by least-squares analysis of variance (Harvey, 1990), so as to be included in the final model. The fixed effects considered were period of calving (8 periods), season of calving (3 seasons), parity of dam (7 parities) and genetic groups of animals (12 groups). Animals calved from 1981 through to 2020 with 38, 103, 139, 195, 184, 224, 224 and 321 calves born in each 5-year period of calving. Calving season discriminated between calves born in winter (n = 476), summer (n = 480) and rainy (n = 472). Parity of the dam ranged from 1 to 7 with 422, 314, 245, 161, 115, 87 and 84 dams recorded under the parity 1, 2, 3, 4, 5, 6 and 7 or above, respectively. Animals of 12 genetic groups had 41 to 228 animals in each group. All these effects except genetic groups were significant (P < 0.05) for all lactation curve parameters and traits and were retained in the final model for genetic parameters estimation. Estimates of variance components and genetic parameters (heritability, repeatability and genetic correlations) for lactation curve parameters and traits of pooled lactations (up to 7^th lactations) were obtained by fitting multi-trait repeatability animal model (detailed in online Supplementary File) using Bayesian method under GIBBS3F90 programs of BLUPF90 family program (Misztal et al., 2018).

The posterior densities of covariance components and marginal posterior distribution for each parameter were obtained by Gibbs sampling. The Markov chain was subjected to 2 × 10⁶ iterations with 6 × 10⁵ rounds burn in period and a thinning interval of 1000 iterations. Hence a total of 1400 effective samples were generated which provide the best convergence of Gibbs chains for all the parameters under study monitored through graphical inspection (trace-plots and histograms). The program also calculates the 95% high posterior density (HPD 95%) interval for all parameters from the individual marginal posteriors. The HDP 95% interval provides a region in which 95% of samples fall into and is a measure of reliability.

Results and discussion

Out of the total of 1718 individual lactation records, 16.88% revealed atypical lactation curve due to negative estimates of parameter ‘b’ and/or ‘c’. Hence, a total of 1428 individual lactation records of Jersey crossbred cows showing typical lactation curve were subjected to further genetic analysis. The occurrence of atypical lactation curve, with absence of the peak yield in lactation in Jersey crossbred cattle in our study were within the range as described by Rekik et al. (2003) and Atashi et al. (2009) in Holstein-Friesian (HF) cattle of different regions. The descriptive statistics of lactation curve parameters/traits under the study are depicted in Table 1. Milk yield increased from calving to the peak milk yield of 11.2 kg at day 50.3, and then decreased gradually until dry off with a persistency of 3.2 for these Jersey crossbred cattle.

Table 1. Descriptive statistics for lactation curve parameters/traits of Jersey crossbred cows

^aa, initial milk yield after calving; b, ascending slope up to peak yield; c, declining slope after peak yield; y _max, peak milk yield; t _max, time of peak yield and P, persistency of milk yield.

Lactation curve parameters

Estimates of posterior means, medians and modes of variance components and genetic parameters of different lactation curve parameters of Wood's function for Jersey crossbred cattle using Bayesian method are shown in Table 2. The efficient convergences for heritability and repeatability estimates of different lactation curve parameters based on density of distribution and trace of iteration process during Bayesian analysis are shown in Supplementary Figures S1 and S2, respectively. The direct heritability estimates for lactation curve parameters, ‘a', ‘b' and ‘c’ were 0.18, 0.09 and 0.11, respectively. Low heritability estimates of curve parameters of Jersey crossbred cattle in the present findings are in agreement with the result of Rao and Sundaresan (1979), Gama et al. (1994) and Yilmaz et al. (2011) in Sahiwal, Holstein and Brown Swiss cattle breeds, respectively. Compared to our study, low estimates of heritability (0.02–0.04) for different lactation curve parameters of Wood function were reported in Holstein cattle (Torshizi, 2016; Saghanezhad et al., 2017; Radjabalizadeh et al., 2022). However, other researchers (Varona et al., 1998; Rekaya et al., 2000) found moderate estimates of heritability for different lactation curve parameters ranging from 0.17 to 0.44 in dairy cattle. These differences could be due to differences in breed, location or the estimation methods. The low direct heritability estimates of the lactation curve parameters, ‘a’, ‘b’ and ‘c’ of Wood's incomplete gamma function indicate that these lactation curve parameters are largely influenced by environmental variation. Due to lack of substantial additive genetic variations between animals for these parameters, direct selection based on lactation curve parameters may not effectively change the shape of the lactation curve of Jersey crossbred cows under the prevailing management system.

Table 2. Posterior means, medians and modes of direct genetic, permanent environmental, residual and phenotypic variance components and genetic parameters for lactation curve parameters of Jersey crossbred cattle

^aσ²_a, additive genetic variance; σ²_pe, permanent environmental variance; σ²_e, residual variance; σ²_p, phenotypic variance; h ², heritability; pe ², permanent environmental effect; r, repeatability; PSD (Posterior standard deviation) and HPD (Highest posterior density)

The proportion of variance associated with permanent environmental effect of animals (pe²) for parameters ‘a', ‘b' and ‘c’ were estimated as 0.13 ± 0.04, 0.12 ± 0.04 and 0.19 ± 0.04, respectively. Repeatability estimates for the different lactation curve parameters (‘a', ‘b' and ‘c’) of animals were moderate in nature, ranging from 0.21 to 0.31 (Table 2). Moderate estimates of repeatability of lactation curve parameters were also reported by Osorio-Arce and Segura-Correa (2005) and Gebreyohannes et al. (2013) in crossbred cows. However, several researchers (Boujenane and Hilal, 2012; Gebreyohannes et al., 2013; Saghanezhad et al., 2017) reported low estimates of repeatability for lactation curve parameters of Wood's function in different cattle breeds. Our moderate repeatability estimates indicate that these parameters were influenced by both permanent and temporary environment, and selection based on early measurement on these parameters could be useful to improve the shape of the lactation curve in the Jersey crossbred cow population.

Lactation curve traits

Table 3 shows the estimates of posterior means, medians and modes of variance components and genetic parameters for lactation curve traits of Jersey crossbred cattle. Convergence for heritability and repeatability estimates of lactation curve traits based on density of distribution and trace of iteration process during Bayesian analysis are shown in Supplementary Figures S3 and S4, respectively. Direct heritability estimates for different lactation curve traits were observed to be 0.24 ± 0.05 for peak milk yield (y _max), 0.12 ± 0.04 for time of peak milk yield (t _max) and 0.15 ± 0.05 for persistency (P) of milk yield of Jersey crossbred cattle using Bayesian method (Table 3). Moderate heritability of peak milk yield was also reported by Saghanezhad et al. (2017) and Atashi et al. (2019) in HF cattle and Yilmaz et al. (2011) in Brown Swiss cattle, which supports our findings. However, higher (Rekik et al., 2006) and lower (Boujenane and Hilal, 2012; Torshizi, 2016; Bakri et al., 2022) estimates of heritability for peak milk yield were also reported in different cattle breeds. Low heritability estimates for time of peak yield and persistency of milk yield were also reported by Rekik et al. (2006), Wasike et al. (2014) and Atashi et al. (2019) which is in agreement with the results of this study. However, comparative higher (Yilmaz et al., 2011) and very low (Saghanezhad et al., 2017; Bakri et al., 2022; Pangmao et al., 2022) estimates of heritability than the present findings have been reported. Low-to-moderate heritability estimates of different lactation curve traits, i.e., y _max, t _max and P in this study suggest the existence of some genetic variations for these traits in Jersey crossbred population and hence slow genetic improvement of these traits could be possible through appropriate selection strategies.

Table 3. Posterior means, medians and modes of direct genetic, permanent environmental, residual and phenotypic variance components and genetic parameters for lactation curve traits of Jersey crossbred cattle

^a See Table 2 for abbreviations.

Estimates of pe² for y _max, t _max and P in this study (Table 3) were 0.15 ± 0.04, 0.12 ± 0.03 and 0.22 ± 0.04, respectively. Compared to the present findings, Wasike et al. (2014) observed the permanent environmental variance as a proportion of phenotypic variance (pe²) for milk yield at peak lactation, days in milk at peak lactation and persistency to be 0.03, 0.08 and 0.25, respectively. in dairy cattle. The repeatability estimates of lactation curve traits of Jersey crossbred cows were estimated as moderate which ranged from 0.24 to 0.39 (Table 3). Various researchers (Saghanezhad et al., 2017; Atashi et al., 2019; Pangmao et al., 2022) obtained moderate estimates of repeatability for peak milk yield in different dairy cattle breeds. However, they observed lower estimates of repeatability for time of peak yield and persistency of milk yield than the present findings, which could be attributed to differences in animal's genetic makeup, variation in data structure of the present herd and difference in analytical methodology. Moderate repeatability estimates of different lactation curve traits (peak milk yield, time of peak yield and persistency of milk yield) of Jersey crossbred cattle in our study indicate the importance of both permanent and temporary environment in observed phenotypic variation of these traits. Hence, selection based on early measurement of these traits in a lactation can improve the persistency and peak milk yield as well as time of peak yield.

Correlations among lactation curve parameters

Table 4 depicts the genetic and phenotypic correlations of lactation curve parameters with lactation curve traits in Jersey crossbred cattle. Genetic correlations between parameters ‘a’ and ‘b’ and between ‘a’ and ‘c’ were negative and moderate in nature, ranging from −0.43 to −0.53. Similarly, negative phenotypic correlations between parameters ‘a’ and ‘b’ (−0.43) and between ‘a’ and ‘c’ (−0.13) were also observed. Further, high and positive correlations (0.82) existed between the curve parameters ‘b’ and ‘c’ at both genetic and phenotypic level (Table 4). Similar to the present findings, negative genetic correlations (Rao and Sundaresan, 1979; Gama et al., 1994; Rekaya et al., 2000) and negative phenotypic correlations (Chegini et al., 2015; Saghanezhad et al., 2017; Radjabalizadeh et al., 2022) were observed between ‘a’ and ‘b’ as well as between ‘a’ and ‘c’ with a wide range of magnitudes in different breeds of cattle. Positive and high genetic correlations between ‘b’ and ‘c’, as observed in the current study, were also reported by Rekaya et al. (2000), Chegini et al. (2015) and Saghanezhad et al. (2017) in dairy cattle. In contrast to our present findings, positive phenotypic correlation between parameter ‘a’ and ‘c’ was reported by Rekaya et al. (2000), Farhangfar and Rowlinson (2007) and Boujenane and Hilal (2012) in dairy cattle. Regarding the associations of ‘b’ and ‘c’, several researchers (Farhangfar and Rowlinson, 2007; Boujenane and Hilal, 2012; Radjabalizadeh et al., 2022) detected strong phenotypic relationships between parameter ‘b’ and parameter ‘c’, which supports our findings. In our study, moderate and negative genetic correlations between parameters ‘a’ and ‘b’ and between ‘a’ and ‘c’ indicate that cows with higher initial yields had a slower rate of incline until peak yield and a slower rate of decline. The strong and positive genetic and phenotypic correlations (0.82) between parameters ‘b’ and ‘c’ in our study implies that Jersey crossbred cows with a higher ascending slope of milk production would be expected to have a faster rate of decline of milk production after peak yield at phenotypic level.

Table 4. Genetic (below diagonal) and phenotypic correlations (above diagonal) among lactation curve parameters and lactation curve traits of Jersey crossbred cows

^a See Table 1 for abbreviations.

Correlations among lactation curve traits

Genetic associations among different lactation curve traits, i.e., y _max, t _max and P were found to be positive and ranged from 0.08 to 0.61 (Table 4). However, weak and negative phenotypic correlations were observed between y _max and t _max (−0.19) and t _max and P (−0.09) in our study. Further, positive and moderately high phenotypic correlation (0.60) was also observed between t_max and P (Table 4). Our findings are in agreement with those of Boujenane and Hilal (2012), Saghanezhad et al. (2017) and Pangmao et al. (2022), who observed positive genetic associations between the lactation curve traits (peak yield, time of peak yield and persistency) in Moroccan Holstein, Holstein crossbred and Iranian Holstein cattle, respectively. Contrary to the present findings, Gama et al. (1994) and Albarrán-Portillo and Pollott (2008) found a negative phenotypic correlation between persistency (P) and time of peak yield (t _max) in dairy cattle. In this study, positive genetic associations among the lactation curve traits revealed that selection of one trait will reward correlated response to other lactation curve traits in Jersey crossbred cattle. Further, high and positive correlation (0.60) of persistency (P) with time of peak yield (t _max) at both genetic and phenotypic level in our study indicate that cows that reach peak yield later during lactation would have higher persistency. Weak antagonistic phenotypic relationships of peak milk yield with persistency and time of peak yield did exist but improvement of peak yield may not affect the persistency and time of peak yield in these Jersey crossbred cattle.

Correlations of curve parameters with curve traits

Association study of lactation curve parameters with lactation curve traits (Table 4) revealed that parameter ‘a’ was negatively correlated with time of peak yield (−0.23), positively with persistency (0.24) and peak yield (0.95) at genetic level. Further, the parameter ‘a’ had negative phenotypic correlations with P (−0.28) and t _max (−0.55) but positive correlation with y _max (0.89). Similarly, high positive genetic correlation of parameter ‘a’ with peak yield (Farhangfar and Rowlinson, 2007; Gebreyohannes et al., 2013), moderate positive genetic correlation with persistency (Gama et al., 1994) and moderate negative genetic correlation with time of peak yield (Rekaya et al., 2000; Gebreyohannes et al., 2013; Chegini et al., 2015) have been reported in different breeds of cattle. Positive phenotypic correlation of parameter ‘a’ with peak milk yield and negative phenotypic correlations of parameter ‘a’ with time of peak yield and persistency, as evidenced from this study, were also reported by Gama et al. (1994), Farhangfar and Rowlinson (2007) and Boujenane and Hilal (2012). On contrary, negative phenotypic correlation of parameter ‘a’ with peak milk yield (Batra et al., 1987) and positive correlation with time of peak milk yield (Rekaya et al., 2000; Gebreyohannes et al., 2013) and persistency (Rao and Sundaresan, 1979) were reported.

Genetic and phenotypic correlations of parameter ‘b’ with y _max (−0.39 and −0.18), t _max (0.50 and 0.46) and P (−0.30 and −0.36) were observed in Jersey crossbred cattle (Table 4). Gama et al. (1994) also found negative genetic associations of parameter ‘b’ with lactation curve traits except time of peak yield which agrees with the present findings. Regarding the associations of parameter ‘c’ with lactation curve traits, it was found that parameter ‘c’ was negatively correlated with all lactation curve traits at genetic and phenotypic level. Genetic correlation estimates of parameter ‘c’ were low for time of peak yield (−0.02), moderate for peak yield (−0.47) and high for persistency (−0.75), whereas the corresponding values for phenotypic correlations with these traits were −0.001, −0.10 and −0.71, respectively (Table 4). Similar to the present result, Rao and Sundaresan (1979), Boujenane and Hilal (2012) and Chegini et al. (2015) obtained moderate to high and negative estimates of genetic correlations of parameter ‘c’ with peak yield and persistency in different cattle breeds. Antagonistic relationships of parameter ‘c’ with time of peak yield (Rekik et al., 2003; Farhangfar and Rowlinson, 2007) and persistency (Rekaya et al., 2000; Farhangfar and Rowlinson, 2007; Boujenane and Hilal, 2012) were reported in various dairy cattle breeds, which were consistent with our findings. In contrast to our study, Tekerli et al. (2000) and Chegini et al. (2015) found positive estimates of phenotypic correlations of parameter ‘c’ with different lactation curve traits. In our study moderate negative genetic correlations of parameter ‘b’ and ‘c’ with peak and persistency indicate that selection against parameter ‘b’ will result in improvement of the shape of lactation curve. However, selection based on parameter ‘b’ could shorten the peak yield time due to high positive relationship of parameter ‘b’ with time of peak yield which may affect the 305 d milk yield of animals. Moderate-to-high negative genetic correlations of parameter ‘c’ with peak yield and persistency indicate that prediction of lactation curve based on parameter ‘c’ can improve lactation persistency and peak yield of Jersey crossbred population at genetic level.

In conclusion, it was observed from the present study that there exist genetic variations of lactation curve parameters/traits of Wood's incomplete gamma function. However, most of the traits under study revealed low magnitude of heritability and a large influence of environmental factors, hence direct genetic selection based on lactation curve parameters/traits will lead to only slow genetic improvement, and it will take longer to change the shape of lactation curve of Jersey crossbred cows. Moderate repeatability estimates of different lactation curve parameters/traits indicated that there is scope to select animals with better performance in early parities. Genetic and phenotypic correlations between lactation curve parameters and lactation curve traits in Jersey crossbred cattle were mostly moderate to high in magnitude in our study. This indicates that there are substantial genetic and physiological relationships between lactation curve parameters and lactation curve traits of Jersey crossbred cattle. These need to be considered in genetic improvement programs of the Jersey crossbred population.