Introduction
The Nublang is a dual-purpose breed (milk and draft) of Bos indicus cattle. It is a genetically unique breed (Dorji et al., Reference Dorji, Hanotte, Arbenz, Rege and Roder2003; Sharma et al., Reference Sharma, Maitra, Singh and Tantia2013) that evolved from crossing of humpless cattle (Bos taurus) from Tibet and humped cattle (B. indicus) from India (Dorji et al., Reference Dorji, Kaneda, Lin, Takahashi, Oyama, Sasazaki, Yamamoto, Kawamoto and Mannen2010). They are highly adapted to harsh and varied environments ranging from sub-tropical to cool temperate climate and between an altitude range of 200–3 000 m asl. They are also relatively tolerant to some local ecto-parasites and diseases (DAD-IS, 2014).
Nublang cattle are central to household nutrition and subsistence agriculture prevalent in remote and rural areas of Bhutan. They are found throughout the country with the most concentrated populations in the districts of Haa, Chukha and Samtse. They are also in the Indian states of Sikkim and Darjeeling district of West Bengal where they are known as Siri. In the past, their number symbolized the wealth of the family. They used to be sacrificed during annual household rituals and construction of new houses. They are also worshipped by Hindu communities in Bhutan as a Goddess of Wealth (Arbenz and Tshering, Reference Arbenz and Tshering2000).
The milk yield of Nublang cattle is low and for this reason, they were traditionally crossbred with Mithun (Bos frontalis) to produce Jatsham (F1 female progeny) for better milk yield (Dorji, Tshering and Rai, Reference Dorji, Tshering, Rai, Dhawan and Koyyana2009). From the 1960s onwards they have also been extensively crossed with high yielding B. taurus breeds (Jersey and Brown Swiss cattle). Consequently, the quality and population size of purebred Nublang has declined (RNRRC Jakar, 2008). Currently the purebred Nublang population constitutes <35 percent of total cattle population of 298 916 in the country (DoL, 2013). However, the current population is classified as “breed not at risk” (DAD-IS, 2014).
Genetic parameter estimates are specific to a trait, population and environment where measured and it is not appropriate to generalize them to other populations. Although heritability (h 2) and repeatability (r) estimates for milk production traits have been extensively reported for temperate and tropical dairy breeds, they are lacking for Nublang cattle. This paucity is likely due to the lack of a functional herd recording system for this breed.
Therefore, this study assessed the phenotypic variation and genetic parameters for milk yield in Nublang cattle to gain insights into the potential for improving milk yield through selective breeding, towards conservation and sustainable utilization of the breed.
Materials and methods
Data and milk recording
Milk yields and pedigree records of Nublang cows from the National Nublang Breeding Farm, Tashiyangphu, Bhutan were used for the study. The farm was established in 1994 for conservation of purebred Nublang. The animals were maintained on kikuyu grass (Pennisetum clandestinum) pasture in the summer and fed its hay and silage during winter. Supplementary feeding of 1 kg of concentrate feed mixture in summer and 2 kg during the winter were provided to the milking cows irrespective of level of production. Both artificial insemination and natural service were used on the farm and calving occurred year round.
The final clean data consisted of 2 052 test day (TD) milk yields of 66 milking cows over a period from 1997 to 2013. These cows were sired by 18 bulls: ten sires with less than three daughters and eight sires with three to ten daughters. TD milk yields were recorded monthly on the 30th day of each month except for the 28th in February. Cows were milked in morning and evening and, immediately prior to milking, their calves were allowed to suckle from one udder quarter to initiate the milk let-down reflex and removed from their mother. TD milk yield was the sum of morning and evening milk yields. The average number of cows lactating per year was 15.
Total lactation milk yield (LMY) for each cow was estimated from TD milk yields and the interval in days between subsequent milk recordings using the Test Interval Method (ICAR, 2009). Overall there were 261 lactation records (average of four lactations per cow) available considering only lactation with TD records of four and above for calculation of LMY:
$$\eqalign{{\rm LMY} & = {I_0}{M_1} + \displaystyle{{{I_1}^\ast({M_1} + {M_2})} \over 2} \cr &\quad + \displaystyle{{{I_2}^\ast({M_2} + {M_3})} \over 2} \ldots \ldots \displaystyle{{{I_{n - 1}}^\ast({M_{n - 1}} + {M_n})} \over 2} + {I_n}{M_n}}$$
where
LMY is the estimated total milk yield for the lactation period of (I 0 + I 1 + I 3 + … + I n ); I 1, I 2, I n−1 are the intervals, in days, between TD recording dates; I 0 is the interval, in days, between the lactation period start date and the first recording date; I n is the interval, in days, between the last recording date and the end of the lactation period.
M 1, M 2, M n are milk yields in litre, given to one decimal place of the milk yielded in the 24 h of the recording day. That is, M 1 is the milk yield for TD1, M 2 for TD2, etc.
Data analysis
Phenotypic analyses of the data were done using univariate GLM in Statistical package for Social Sciences (SPSS, 2007). The final statistical model used for this analysis was;
$${\rm DM}{{\rm Y}_{ijklm}} = \mu + {\rm Pa}{{\rm r}_i} + {\rm TD}{n_j} + {\rm Sea}\_\,{r_k} + {\rm Yea}{{\rm r}_i} + {e_{ijklm}}$$
where
DMY ijklm , the observation of daily milk yield on the mth cow that on ith parity, jth TDn, kth season and lth year; μ the overall mean; Par i the ith parity subclasses (i = 1 to ≥ 6); TDn j the jth month in milk (j = 1 to ≥9); Sea_r k the kth season of milk recording (k = spring, March–May; summer, June–August; autumn, September–November; and winter, December–February).
Year l the lth year (i = 1997–2013); e ijklm the residual associated with y ijklm .
Where there were few observations within a factor level, these levels were combined: this included combining parity7, parity8 and parity9 with parity6 as parity≥6 and TD9, TD10 and TD11 into TD ≥ 9.
The factors fitted for LMY were number of TD records per lactation (NTD) with six levels (4 to ≥9), as well as parity and season of calving (sea_c) with levels as described for DMY. Wherever a factor had a significant effect, a pairwise comparison of level means was carried out using least significant difference (LSD) (SPSS, 2007).
Estimation of variance components and genetic parameters
The variance components of DMY and LMY were estimated using an Average Information Restricted Maximum Likelihood mixed model procedure with ASReml statistical package version 3 (Gilmour et al., Reference Gilmour, Cullis, Welham and Thompson2009) using the following univariate linear repeated measure animal model;
$${\bf y} = {\bf Xb} + {{\bf Z}_1}{\bf a} + {{\bf Z}_{\bf 2}}{\bf p} + {\bf e}$$
where
y the phenotypic repeat measures of the trait (DMY, LMY) where TD milk yields were treated as repeated measures of DMY; b the vector of fixed effects; a the vector of random additive genetic effects for all animals, assumed to be distributed as N (0,2 g A) where A is the pedigree relationship matrix and σ2 g is the additive genetic variance; p the vector of permanent environmental effects (PE) for all animals with milk records, and; e the vector of residual effects.
The matrices X and Z are design matrices that associate observations to particular fixed or random effects, respectively.
The final model for DMY included fixed effects (parity, TD n and sea_r) and random effects of additive genetic effects, permanent environment and residuals. Similarly the fixed effects for LMY included parity, NTD and sea_c, and random effects as for DMY. The variance components from these models were used for estimation of heritability and repeatability using standard formula from Falconer and Mackay (Reference Falconer and Mackay1996) in ASReml software package.
Results
The mean DMY of Nublang cows in the farm was 2.12 ± 0.7 litres (N = 2 052, range 0.3–5.0 litres). The mean LMY was 519.2 ± 151.4 litres (N = 261, range 115–881.7 litres) with an average lactation length of 239.6 days. The variation in these traits, as indicated by the coefficient of variations (DMY 34 percent, LMY 29.1 percent), is high. Parity, season and stage of lactation (TDn) significantly affected DMY (P < 0.05). For example, least-squares mean DMY in the first parity was 1.84 ± 0.04 litres and increased thereafter to peak in the fourth parity (2.37 ± 0.04 litres) (Table 1). Mean DMY was highest in summer (2.40 ± 0.03) and lowest in winter (1.81 ± 0.04) (P < 0.05). For stage of lactation, (as measured by the TDn), DMY peaked on TD2 (2.54 ± 0.04) and then declined gradually up to TD ≥ 9 (1.70 ± 0.07) (Figure 1). The effect of interaction between parity and TD n , parity and year, and season and year were significant (P < 0.05).
LSM Daily milk yield of Nublang cows by stages of lactation.
Least-squares means and standard error (SE) of DMY and LMY of Nublang cows.
N, no. of observations; NC, no. of cows; DMY, daily milk yield; LMY, lactation milk yield.
*Season = season during milk recording for DMY (Sea_r); season of calving for LMY (Sea_c).
a,b,c,dLeast-squares mean within a column group with different letters differ significantly (P < 0.05).
Total lactation milk yields were significantly affected by parity and NTD (P < 0.05). It was highest in the fourth parity (572.5 ± 20.90 litres) and lowest in the first parity (411.6 ± 19.68 litres) (Table 1). The LMY of spring calvers (456.3 ± 14.70 litres) was significantly low compared to other seasons. As expected, NTD and LMY are highly correlated because an increasing number of TD milk yields per lactation corresponds to a longer lactation and therefore a higher LMY (Figure 2). The effect of interaction between parity and NTD was significant (P < 0.05) among other factors.
LSM Lactation milk yield of Nublang cows by number of test records in a lactation.
The h 2 estimates for DMY and LMY were 0.22 ± 0.16 and 0.13 ± 0.20, respectively. The repeatability of DMY and LMY were 0.44 ± 0.05 and 0.49 ± 0.08, respectively.
Discussion
Milk yield
The milk yield of Nublang in this study is low compared with its crossbreds with Jersey (4.9 litres) and Brown Swiss cattle (4.3 litres) in the country (Tamang and Dorji, in press). Further, it is lower than some of the other B. indicus tropical dairy breeds; Sahiwal (1 200–1 800 kg) (Bajwa et al., Reference Bajwa, Khan, Khan and Gondal2004) and Tharparkar (2 179 litres) (Gahlot, Pant and Kachawaha, Reference Gahlot, Pant and Kachawaha2000) but comparable to the dual purpose Ongole breed of cattle (675 litres) (Gaur, Kaushik and Garg, Reference Gaur, Kaushik and Garg2002). Within the Nublang breed, the milk yield in this study was similar to that found in Tamang and Perkins (Reference Tamang and Perkins2005) (2.2 litres). However, there is wide variation in reported milk yields of Nublang; 1.5 litres (Arebenz and Tshering, Reference Arbenz and Tshering2000), (2.5 litres) Phanchung and Roden (Reference Phanchung and Roden1996) and 5.0 litres (Dorji et al., Reference Dorji, Tshering, Dorji, Dhendup, Tshering and Chettri2005) which is similar to the range observed across cows in the present study. It is likely that the differences between studies may be attributed to the differences in animal management, environment, sampling methodology and genetic merit of the animals. Overall, the milk yield in Nublang cattle is low yet highly variable reflecting a similar pattern seen in other B. indicus breeds in general (Syrstad, Reference Syrstad1993).
Within a lactation, peak yield occurred around TD2 which is within 2 months from date of calving. In other B. indicus breeds, peak milk yield was observed in about 1 month post calving in Tharparkar (Balbir, Pankaj, and Pathodiya, Reference Balbir, Pankaj and Pathodiya2011) and about 7 weeks into milk in Sahiwal (Dongre and Gandhi, Reference Dongre and Gandhi2014). The increase in milk yield from first to fourth parity in this study is in line with a number of studies (Rekik et al. Reference Rekik, Ben Gara, Ben Hamouda and Hammami2003; Gradiz et al. Reference Gradiz, Alvarado, Kahi and Hirooka2009; Amasaib et al. Reference Amasaib, Fadel-Elseed, Mahala and Fadlelmoula2011; Jingar et al. Reference Jingar, Mehla, Singh and Roy2014), which according to Gradiz et al. (Reference Gradiz, Alvarado, Kahi and Hirooka2009), is a consequence of the increasing development of the udder, increase in body size, increase in feed intake and possibly improved partitioning of energy into milk production.
The seasonal variation in milk yields with higher yields in summer in this study is in agreement with Rekik et al. (Reference Rekik, Ben Gara, Ben Hamouda and Hammami2003), where animals were mainly maintained on pasture in summer and conserved fodder in winter. Such seasonal variation may be capitalized upon to produce more milk by planning breeding programs to coincide with availability of high-quality fodder in summer to increase milk production. Scheduling mating programmes when cows are in good condition and most likely to become pregnant should be considered. Although breeding occurred all year round in this study, the number of calvings in spring (95) was higher than other seasons (45–69). This may be due to better conception rates in summer due to availability of quality fodder and better nutrition.
Genetic parameters
The h 2 estimates for milk yield in this study suggest a moderate additive genetic effect on milk production, although standard errors were large. They suggest that genetic progress can be made by selection for these traits. The estimates were within the range for other B. indicus cattle (0.05–0.40) reviewed by Rehman et al. (Reference Rehman, Khan, Rehman and Hassan Faiz2013). The low accuracy of the estimates may be attributed to smaller datasets as in Eid et al. (Reference Eid, Elsheikh, Ibrahm and Yousif2012) and Bilal et al. (Reference Bilal, Khan, Bajwa and Shafiq2008) and inaccurate recording of milk yield or errors in pedigree records (Javed, Mohiuddin and Ahktar, Reference Javed, Mohiuddin and Ahktar2001). However, for better accuracy of estimates, analysis of a larger dataset is recommended.
The medium repeatability estimates in this study are within the range reported for Sahiwal cattle (0.13–0.50) reviewed by Rehman et al. (Reference Rehman, Khan, Rehman and Hassan Faiz2013). This suggests that information from first parity could be used as a reasonable prediction of future milk yields, thereby hastening the process of selection and genetic gain. This is useful information for bull selection where breeding value for milk yield would generally be evaluated through the time-consuming process of progeny testing.
Conclusions
The results of the present study indicate potential for genetic improvement of milk production within the Nublang breed. However, for better accuracy of heritability estimates, computation using a larger dataset of Nublang cattle is recommended in future. Performance recording of Nublang cows in villages in addition to the government farm is suggested for obtaining adequately large datasets.
Acknowledgements
The authors would like to thank the Department of Livestock, Ministry of Agriculture and Forests, Bhutan for permission to work on the Nublang dataset and Integrated Livestock Crop Conservation Project (ILCCP), Bhutan for financially supporting this study. We are also indebted to anonymous reviewers for useful suggestions and editing to improve the manuscript.
Statement of interest
The authors do not have any conflict of interest that will influence the judgment and potentiality of being biased. Further, the authors did not have any financial arrangements or connections pertinent to submitted manuscript.
