Introduction
The performance of farm animals is strongly influenced by the thermal environment. Ambient temperatures below or above the thermo-neutral range create stress conditions in animals. The thresholds (both cold and hot environmental challenges) beyond which potential thermal stressors can influence performance and health are influenced by various factors such as animal species, genetics, age or life stage, level of nutrition and prior conditioning (Deshazer et al., Reference Deshazer, Hahn, Xin and Deshazer2009). Among farm animals, the upper critical temperature of dairy cattle is lower than other livestock species (Wathes et al., Reference Wathes, Jones and Webester1983). The approximate thermal-comfort zone for optimum performance of adult cattle is reported to be 5–15 °C (Hahn, Reference Hahn1999). However, significant changes in feed intake or in numerous physiological processes will not occur within the range of 5–25 °C (McDowell, Reference McDowell1972).
In temperate climatic conditions, heat stress in dairy cattle during certain times of the year is also a major concern for the dairy industry as it has detrimental effects on feed intake, growth efficiency, production and reproduction of animals (Ravagnolo et al., Reference Ravagnolo, Misztal and Hoogenboom2000; West et al., Reference West, Mullinix and Bernard2003; Gantner et al., Reference Gantner, Mijic, Kuterovac, Barac and Potocnik2015). In the tropical and sub-tropical climatic conditions found in India, the susceptibility of dairy animals to heat stress would be even higher as the average annual temperature in most parts is 25 °C or higher, that is at or above the thermal-comfort zone for maximum productive performance of these animals (Sirohi and Michaleowa, Reference Sirohi and Michaelowa2007). Studies in Africa (Du Preez et al., Reference Du Preez, Giesecke and Hatiingh1990), Australia (Davison et al., Reference Davison, McGowan, Mayer, Young, Jonsson, Hall, Matschoss, Goodwin, Goughan and Lake1996) and the Indian sub-continent (Zewdu et al., Reference Zewdu, Thombre and Bainwad2014) provide ample empirical evidence for the negative association between cattle productivity and heat stress in the tropics and sub-tropics.
The literature is, however, mostly focused on the physiological responses of cattle of predominantly Bos taurus species or its crosses, with limited evidence on native dairy breeds in the tropical countries, especially the buffalo. Buffalo have been an integral part of animal agriculture for centuries. Of the two generic types, African (Syncerus caffer) and Asian (Bubalus spp.), Asian water buffaloes (Bubalus bubalis) have been heavily domesticated and are now widespread in around 42 countries across the world ranging from Indian subcontinent to parts of south and central America.
The morphological and anatomical characteristics of buffaloes make them well-suited to hot and humid climates of the tropics and sub-tropics. The presence of numerous melanin pigments in the epidermis prevents ultraviolet rays from penetrating through the dermis of the skin to the lower tissues and is an important anatomical feature of the species making them tolerant to heat stress (Shafie, Reference Shafie and Yousef1985). This beneficial characteristic is reinforced by well-developed sebaceous glands, with greater secretion activity than in cattle. The sebum layer secreted by the glands melts during hot weather and becomes glossier to reflect many of the heat rays, thus relieving the animal from the excessive external heat load. Despite this, it is generally believed that buffaloes are less physiologically adaptive to heat stress than cattle, owing to thick black skin colour with sparse hair coat and fewer deeply situated sweat glands (Marai and Haeeb, Reference Marai and Haeeb2010).
Only a few studies carried out in India report the detrimental effect of heat stress on the reproduction of buffaloes (Kaur and Arora, Reference Kaur and Arora1982; Tailor and Nagda, Reference Tailor and Nagda2005) and their productive performance (Upadhyay et al., Reference Upadhyay, Singh, Kumar and Gupta2007). Therefore, since the empirical work on the extent of production losses in buffaloes accounting for variables other than climatic conditions and during different levels of heat stress is limited, the objective of this study was to establish the relationship between weather conditions and production performance of lactating buffaloes.
Materials and methods
Data
Data were obtained from Livestock Research Centre (LRC), National Dairy Research Institute (NDRI) located in Karnal district (latitude 29°41′0″N, longitude 76°59′0″E) of the Indian state Haryana. With mean annual temperature of around 25 °C, and two distinct seasons viz., hot summer and cold winter, the weather conditions in the district are typical of those prevalent in vast expanse of north Indian plains.
The data comprised fortnightly records of milk yield, calving date, age at first calving and stage of lactation for the period 2001–10 of Murrah breed of buffaloes. The data was collected from one single research farm, which in climatic conditions is typical of the home tract of this highly preferred dairy breed which stretches around the southern parts of the Haryana state. All the animals in the LRC were maintained in loose housing barns with no specific shelter management for ameliorating environmental stress. There was joint feeding system on farm, animal ration was replenished in the mangers twice a day, morning and evening
Each buffalo was required to have at least 14 fortnight milk yield records (viz. minimum about 210 d in milk in a lactation) to be part of the analysis for that parity. Each month was divided into two fortnights, 1 to 15th and 16th to the remaining days of the month. After cleaning the data for outlier values, the final data set consisted of 12 763 observations for 695 buffaloes.
The data on weather parameters was sourced from the nearest located weather station (within 5 km. radius), maintained by the sister organization, Central Soil Salinity Research Institute (CSSRI), Karnal. Weather variables included daily maximum, minimum, and average temperature and relative humidity (RH). The relative humidity data were for two time periods, viz., 7.22 am and 2.22 pm representing daily maximum and minimum RH, respectively.
Statistical analysis
Estimates of THI
Temperature humidity index (THI) has been extensively used as an indicator of thermal stress (Ravagnolo et al., Reference Ravagnolo, Misztal and Hoogenboom2000; Gantner et al., Reference Gantner, Mijic, Kuterovac, Barac and Potocnik2015) Based on the availability of maximum and minimum temperature and relative humidity the present study worked out THI based on the following equation (Ravagnolo et al., Reference Ravagnolo, Misztal and Hoogenboom2000; St-Pierre et al., Reference St-Pierre, Cobanov and Schnitkey2003): THI = (9/5 temperature °C + 32) − (11/2 − 11/2 × relative humidity) × (9/5 temperature °C − 26). The mean daily and fortnightly temperature recorded for the period 2001–10 ranged between around 18–30 °C with average temperature of around 24 °C (Table 1). The mean daily relative humidity was marginally high than that of fortnightly humidity. Maximum, minimum and average THI were estimated. Incorporating maximum temperature and minimum relative humidity (RH at 2.22 pm) in the equation gave maximum THI (THI_Max), minimum temperature and maximum relative humidity (at 7.22 am) was used to work out minimum THI (THI_Min), and average THI (THI_Avg) was obtained using average temperature and average relative humidity.
Mean temperature and humidity across the study area, 2001–10
The estimated fortnightly maximum, minimum and average THI based on ensemble mean of temperature and relative humidity of each fortnight for the period 2001–10 provides a record of the fortnightly variation of all the three variants of THI (Fig. 1). The flattened bell shaped curvature of the THI curve during the year is not only typical in tropical countries (Upadhyay et al., Reference Upadhyay, Singh, Kumar and Gupta2007; Zewdu et al., Reference Zewdu, Thombre and Bainwad2014) but also reported in the temperate weather conditions (Ravagnolo et al., Reference Ravagnolo, Misztal and Hoogenboom2000). The curve signifies that with passing fortnights THI increases and reaches at its maximum during 13th and 14th fortnights i.e., in month of July, and then gradually starts to decline.
Pattern of temperature humidity index (THI), based on ensemble mean of temperature and relative humidity for the period 2001–10, Karnal.
Analytical model
The extent of production losses due to heat stress may typically be influenced by the animal management practices (du Preez et al., Reference Du Preez, Giesecke and Hatiingh1990) and other factors such as order and stage of lactation, breed, age (Ravagnolo et al., Reference Ravagnolo, Misztal and Hoogenboom2000; West et al., Reference West, Mullinix and Bernard2003). As mentioned earlier, since the animals were from the same breed and maintained under same feeding and management practices, the breed and management effect was controlled. Taking the cue from Ravagnolo et al. (Reference Ravagnolo, Misztal and Hoogenboom2000), three models were used for studying the effect of thermal stress on production of buffaloes.
Model 1 was defined as,
$$\eqalign{{\rm FM}{\rm Y}_{{\rm ift}}{\rm}={\rm } &{\rm} {\rm \beta} _{\rm 0}{\rm}+{\rm } {\rm \beta} _{\rm 1}{\rm AFCi}+{\rm } {\rm \beta} _{\rm 2}{\rm O}{\rm L}_{{\rm ift}}{\rm}+{\rm } {\rm \beta} _{\rm 3}{\rm O}{\rm L}^{\rm 2}_{{\rm ift}} {\rm}+{\rm } {\rm \beta} _{\rm 4}{\rm SO}{\rm L}_{{\rm ift}}{\rm \_E} \cr &{\rm}+{\rm } {\rm \beta} _{\rm 5}{\rm SO}{\rm L}_{{\rm ift}}{\rm \_M}+{\rm } {\rm \beta} _{\rm 6}{\rm THI\_Ma}{\rm x}_{{\rm ft}}{\rm}+{\rm } {\rm \beta} _{\rm 7}{\rm THI\_Ma}{\rm x}^{\rm 2}_{{\rm ft}} {\rm}+{\rm \varepsilon}} $$where, FMYift represents fortnightly milk yield of ith lactating animal (i = 1, 2, …, 396) in fth fortnight (f = 1, 2,.., 24) in tth year (t = 2001, 2002,…, 2010), AFC is age at first calving of i th animal, OLift representing stage of lactation of ith animal in fth fortnight in the year t, THI_Maxft is maximum THI of fth fortnight in the year t. SOL represent stage of lactation of ith animal in fth fortnight in tth year. Three stage of lactation viz., Early, Middle and Late, were defined for every 90 d interval after calving, starting from zero and were captured through dummy variables as,
$$\eqalign{{\rm SO}{\rm L}_{{\rm ift}}{\rm \_E} = &{\rm 1,}\,{\rm for\ early\ stage}\,({\rm up\ to\ 90\ days\ after\ calving}) \cr =\, &{\rm 0,}\,{\rm otherwise}.} $$
$$\eqalign{{\rm SO}{\rm L}_{{\rm ift}}{\rm \_M} = \,&{\rm 1,\ for\ middle\ stage}\,{\rm (91}-{\rm 180\ days\ after\ calving)} \cr = &\,{\rm 0,\ otherwise}{\rm.}} $$The other two models were the same as Model 1, but minimum THI (THI_Min.) and average THI (THI_Avg.) replaced the maximum THI in Model 2 and Model 3 respectively.
The age at first calving (AFC), stage of lactation and stage of lactation served as control variables. AFC is an important trait influencing not only the milk productivity of the animal per year but also the lifetime productivity of the animal. There is an inverse relationship between AFC and milk productivity of dairy animals (Schutz et al., Reference Schutz, Hansen, Steuernagel and Kuck1990; Niloforooshan and Edriss, Reference Niloforooshan and Edriss2004). Therefore, consideration of AFC as a control variable in regression model became important to remove its effect on productivity while estimating the effect of climatic stress on milk yield of animal. Another important non-environmental trait affecting performance of dairy animals is stage of lactation, which has a quadratic relation with milk yield (Torshizi et al., Reference Torshizi, Aslamenejad, Nassiri and Farhangfar2011). The lactation curve of the dairy animal shows that milk yield increases in the early stage (up to 90 d after calving) and starts declining thereafter. Since the stage of lactation varies over the fortnight, inclusion of this variable was considered very important so that the effect of climatic parameters was not unduly overestimated or underestimated.
The models were estimated using the panel data estimation approach. As each cross-section (dairy animal) was neither in-milk during all the 24 fortnights in a year nor was in the dataset for the entire data duration of 10 years, the data set was an unbalanced panel. Both fixed and random effects model were attempted and based on the Hausman test, the fixed effect model was selected (Baltagi, Reference Baltagi2005). The results of one-way period effect were more robust than cross-section effect model as it captured the changes in overall management condition of the entire herd over the time period. The data set tested positive for the problem of heteroscedasticity that would have led to inefficient regression estimates. Hence, to correct for this problem and get best unbiased estimates of the regression parameters, the generalized least square (GLS) estimation was done (Wooldridge, Reference Wooldridge2002).
Results and discussion
The average milk yield of buffaloes for fortnights during the year is given in Table 2. The estimated milk response functions of buffaloes presented in Table 3 indicate that the control and climatic variables were highly significant with expected sign of regression coefficients in all the three models. The relative magnitude of regression estimates was highest for the maximum THI (Model 1) in comparison to the corresponding coefficients in Model 2 (average THI) and Model 3 (minimum THI). Based on the weather data of Georgia state in south-eastern United States, Ravagnolo et al. (Reference Ravagnolo, Misztal and Hoogenboom2000) also concluded that maximum daily temperature and minimum daily humidity were the most critical variables to quantify heat stress in dairy cows. The studies have indicated that the full impact of climatic variables on production may not be instantaneous but delayed due to a time gap between altered feed intake and utilization of consumed nutrients, or changes in the endocrine status of dairy animals. Hence, decline in milk yield per unit increase in the environmental measure was substantially less when evaluated on same day climatic measures in comparison with climatic variables two days earlier (West et al. Reference West, Mullinix and Bernard2003). However, as the present study models the effect on fortnightly milk yield rather than daily productivity such lag effect of climatic variable would not be discernible.
Descriptive statistics of fortnightly productivity of buffaloes
Milk response function of buffaloes
*Significant P < 0.01.
The marginal effect of THI change on the fortnightly milk yield (δFMY/δTHI) at decadal mean level of the THI in each fortnight indicates that buffaloes begin to experience heat stress due to rise in maximum THI (Fig. 2) from the 5th fortnight (early March) and it lasts up to 21st fortnight (mid-November). The decrease in milk yield continues to escalate from 5th to 13–14th fortnight (July) and thereafter the rate of decline gradually decreases till the 21st fortnight. The maximum THI crosses 72 during this duration in the north Indian plains. However, as night temperature (minimum temperature) is still not very high in the month of March (5–6th fortnight) and in early November, the marginal effect of rise in minimum THI on milk yield is negative only during the 7th to 20th fortnight. The marginal effects based on changes in average THI also show a fall in milk yield of buffaloes for a unit increase in average THI during 14 fortnights stretching from April to October (Fig. 2).
Marginal effect of THI on fortnightly milk yield of buffaloes.
The classification of stress zones for dairy animals based on the level of THI has been done in various ways (Armstrong, Reference Armstrong1994; Huhnke et al., Reference Huhnke, McCowan, Meraz, Harp and Payton2001). Broadly, the THI between 72 to 79 is considered as mild stress, 80–89 as moderate stress and ≥90 as severe stress zone. In the Indian dataset also, it was evident that as maximum THI increased beyond 72, mild stress conditions in buffaloes was reflected in slight decline in milk productivity. The moderate stress zone wherein the negative marginal effect on milk production was ≥0.5 l per unit THI rise corresponds to the maximum THI of above 82. Hence, the rate of fortnightly decline in milk yield was estimated for two classes of THI, discomfort period (THI > 72) and stress period (>82) by fitting the log-linear function:
$$\eqalign{{\rm Log} \lpar {{\rm FM}{\rm Y}_{{\rm ift}}} \rpar =\; &{\rm \beta} _{\rm 0}\,{\rm}+{\rm } \,{\rm \beta} _{\rm 1}{\rm AFCi}\,{\rm}+{\rm } \,{\rm \beta} _{\rm 2}\,{\rm O}{\rm L}_{{\rm ift}}\,{\rm}+{\rm } \,{\rm \beta} _{\rm 3}{\rm O}{\rm L}^{\rm 2}_{{\rm ift}} \,{\rm}\cr &+{\rm } \,{\rm \beta} _{\rm 4}{\rm SO}{\rm L}_{{\rm ift}}{\rm \_E} + \,{\rm \beta} _{\rm 4}{\rm SO}{\rm L}_{{\rm ift}}{\rm \_M}\,{\rm}+{\rm } \,{\rm \beta} _{\rm 5}{\rm TH}{\rm I}_{{\rm ft}}\,{\rm}+{\rm } \,{\rm \varepsilon}} $$where FMYift indicates the fortnightly milk yield of ith animal during the fortnights when THI exceeds the level of 72 (discomfort period) and 82 (stress period) in a year t. Similarly, THIft represents the fortnightly maximum, minimum and average THI value in discomfort period and stress period in a year t. The other variables in equation (2) were same as in equation (1). The percentage decline in productivity was worked out as {exp(β5) − 1} × 100.
The results indicated that when maximum THI was >72 (5–21st fortnight), the rate of productivity decline was 0.74% per unit increase in THI (Fig. 3). The rate of decrease was much higher (1.34%), during 9–18th fortnight when THI max >82. The night time discomfort due to minimum THI going beyond 72 sets in from early June (11th fortnight) and continues till mid- September (17th fortnight), but minimum THI never exceeds the stress level of 82. The rate of decline in milk yield resulting from minimum THI of >72 was 0.37% per unit increase in minimum THI (Fig. 3). Working out the productivity decline with respect to the average THI > 72 and >82 indicated a milk yield loss of 0.54 and 0.71% for a unit rise in average THI, respectively (Fig. 3). Using a similar database on daily basis for the period 1994–2004, Upadhyay et al. (Reference Upadhyay, Singh, Kumar and Gupta2007) studied the susceptibility of Murrah buffaloes to climatic stress in the same locale and reported that a sudden decline in minimum temperature (>3 °C) during winter or increase in maximum temperature (>4 °C) during summer adversely affected buffalo milk production. The extent of decline in milk yield was less at mid lactation stage (5–20%) than in early stage (10–30%). In summer season when THI crossed 80, the lactation period of buffaloes was shortened by several days (3–7 d). Estimating a simple linear regression between milk yield and THI, without controlling for animal specific characteristics and period effect, a decline in milk yield by 0.0553 kg/d for a unit increase in THI beyond threshold level was reported in the study.
Rate of decline (%) in fortnightly milk yield (FMY) in the discomfort and stress period.
Based on the findings of response function and marginal effects, the threshold level of maximum, minimum and average THI for the buffaloes from the dataset in our study was estimated to be 74.37, 61.73 and 68.15 respectively. The threshold level indicates the critical level of THI up to which the animal can tolerate the heat stress and after which a significant drop in productivity takes place. Thus, in buffaloes, although a downturn in milk production was discernible after maximum THI crossed 72, the decline was significant only after THI higher than 74. Other studies also more or less substantiate threshold THI of 74+ in case of buffaloes; effect of thermal stress on milk production performance was seen to be relatively small for THI of 75 (Upadhyay et al., Reference Upadhyay, Singh, Kumar and Gupta2007) and this level was also identified as the threshold THI for pregnancy rate of Murrah buffaloes (Dash et al., Reference Dash, Chakravarty, Sah, Jamuna, Behera, Kashyap and Deshmukh2015).
Thus, the better suitability of buffaloes to hot-humid conditions is also reflected in higher critical THI (by 2–3 points) than purebred Bos taurus cows and their crosses. The thermo-neutral level of THI is influenced by many factors including the adaptation level of animals, farm management or housing specialities, nevertheless, empirical evidence shows that stress level for dairy cows sets in at a substantially lower THI of 68 as reported by Gantner et al. (Reference Gantner, Mijic, Kuterovac, Barac and Potocnik2015) for the first parity Holsteins in eastern Croatia. Bouraoui et al. (Reference Bouraoui, Lahmar, Majdoub, Djemali and Belyea2002) observed THI value 69 as threshold value in a Mediterranean climate and Zewdu et al. (Reference Zewdu, Thombre and Bainwad2014) in their study in Maharashtra (India) reported that crossbred dairy cattle were exposed to heat stress when THI crosses the level 72. In terms of the magnitude of losses, the estimates of productivity decline in cows vary widely in different studies, ranging from 0.2 kg per unit increase in THI (Ravagnolo et al., Reference Ravagnolo, Misztal and Hoogenboom2000) to 0.88 kg per THI unit (West et al., Reference West, Mullinix and Bernard2003). In the case of buffaloes the corresponding daily decline is on the lower side, about 0.052 kg per unit increase in THI. One major reason for the lower absolute magnitude of loss is the much lower productivity of buffaloes in India compared to the productivity of dairy cows in the countries where most of the work on heat stress has been carried out. But, even in comparison to crossbred cows of lower productive potential kept in Indian conditions, the production losses from heat stress in buffaloes are lower (Choudhary, Reference Choudhary2017).
Buffalo have formed an important place in the livestock industry of many developing countries like Bulgaria, Philippines, Malaysia, Vietnam, Brazil and Sri Lanka (Mingala et al., Reference Mingala, Villanueva, Cruz and Presicce2017). Owing to its popularity among dairy farmers this breed has spread to almost all parts of the Indian-subcontinent, providing livelihood security to millions of dairy farmers. The animal is suited to the climatic condition in the tropics and the sub-tropics but we can conclude that its production performance is susceptible to heat stress leading to decline in milk productivity of more than 1% per unit increase in maximum THI over 82 in the vast expanse of the north Indian plains. Such stress conditions prevail in this region for about 5 months starting from early May. The maximum crosses the critical threshold of 74 for a longer stretch beginning mid- March to early November. The maximum temperature and minimum humidity that are captured through maximum THI are the most critical weather parameters causing thermal stress in animals. The climatic conditions in this region are such that not only maximum THI but also minimum THI crosses the critical threshold providing little relief to the animals during the night. With the looming threat of climate change, the stress conditions are likely to increase, causing further losses to the dairy farmers unless suitable animal management options are taken to the field level. Buffalo farmers earning cash income from milk production operate on a small gross margin per litre of milk produced (Sirohi et al., Reference Sirohi, Bardhan and Chand2016), which could be further eroded by heat stress. The estimates of marginal effects of THI on milk productivity provided in this study can be useful for designing appropriate managerial interventions for optimizing the production performance of buffaloes in different months of the year and preventing economic losses due to heat stress.
Acknowledgements
We express our sincere thanks to staff of NICRA unit of ICAR-NDRI, Karnal for providing required data to conduct the present research project.
