Relationship among production traits, somatic cell score and temperature–humidity index in the Italian Mediterranean Buffalo

Abstract The temperature–humidity index (THI) has been commonly used to analyse heat stress in dairy cattle, but little is known about its effects on buffaloes. In this study, daily milk yield (MY), fat percentage (FP), protein percentage (PP) and somatic cell count (SCC) data from 808 buffalo cows plus environmental temperature and relative humidity were used to investigate the consequence of heat stress. Two mixed models were used to evaluate the impact of THI on MY, FP, PP and log transformed SCC (SCS). The effect of THI was significant for PP, FP and SCS, whereas its interaction with parity was statistically significant for PP and SCS. The relationship between PP and FP and THI was positive but of different magnitude according to the parity. When THI was below 62, an unfavourable effect was observed, especially in primiparous buffalo cows. A significant interaction between SCS and THI across parities was also observed. The effect of THI on MY across parities was not definite but overall a favourable relationship was observed. Our findings depict a susceptibility of buffaloes to low values of THI, suggesting an optimal THI range for water buffaloes between 59 and 63, although some deleterious effects were observed in primiparous buffaloes at THI values lower than 62. Additional investigations are needed to better elucidate the influence of THI on buffalo species. HIGHLIGHTS The overall effect of THI on buffalo diverges from what commonly observed in dairy cattle Cold stress affects milk and udder health in buffaloes The effect of THI on buffaloes’ performance depends on parity, with a larger susceptibility in primiparous than pluriparous buffalo cows Udder health in buffaloes, evaluated using somatic cell count, is also affected by THI


Introduction
Heat stress can be defined as a condition that occurs when an animal is not able to dissipate an adequate amount of heat to maintain its body thermal balance (Bernabucci et al. 2014). This circumstance has detrimental effects on feed intake, growth efficiency, production and reproduction of dairy species (Hansen 2009;Rhoads et al. 2009;Pragna et al. 2016;Bagath et al. 2019;Summer et al. 2019). Moreover, these effects are expected to increase due to progressive global warming and increased temperature. In fact, in recent years, higher temperatures and humidity, heat waves, and solar flares have led to an important climatic change and global warming causing economic losses in both milk and meat industry (Wankar et al. 2021). The Intergovernmental Panel on Climate Change (IPCC) has confirmed a rise in global temperature by 1.5 C from pre-industrial period and predicted future rise by 0.3 to 4.8 C by the end of twenty-first century (IPCC 2014(IPCC , 2018. To study heat stress in livestock, the temperaturehumidity index (THI), a single value representing the combined effects of air temperature and humidity, is a bioclimatic parameter, which is commonly used to evaluate the degree of heat stress in dairy cattle (Bohmanova et al. 2007;Herbut et al. 2018). Many authors have evaluated the effects of THI and heat stress in dairy cows highlighting a decrease in milk yield and changes in milk's composition (Bernabucci et al. 2010;Nasr and El-Tarabany 2017;Maggiolino et al. 2020).
On the contrary, few studies have evaluated the effects of heat stress on milk yield and quality in buffaloes (Bubalus bubalis) (Upadhyay et al. 2007;Choudhary and Sirohi 2019;Costa et al. 2020c). According to the literature, heat stress seems to severely affect welfare and productivity of buffaloes in tropical environment (Koga et al. 2004;Marai and Haeeb 2010), especially due to some anatomical peculiarities such as few sweat glands, dark skin and rare hairs on the body (Das et al. 1999;Mishra 2021). Nevertheless, buffaloes are mainly distributed in hot and humid climates of tropics and sub-tropics areas. The presence of numerous melanin pigments in the epidermis probably prevents ultraviolet rays from penetrating through the dermis of the skin to the lower tissues and is supposed to be an important anatomical feature against heat stress (Shafie 1985). Furthermore, this peculiarity is also supported by sebaceous glands. These glands secrete sebum, which melts during hot weather and becomes glossier to reflect heat rays, thus relieving the animal from the excessive external heat load (Shafie and El-Khair 1970). Although buffaloes are perfectly suited to their environment, they exhibit signs of great distress when exposed to direct solar radiation or when working in the sun during the warm season. Buffalo body temperature is slightly lower than that of cattle and its skin is usually black and heat absorbent and is only sparsely protected by hair. Furthermore, buffalo skin has a sweat glands density six times lower than cattle skin, making heat dissipation by sweating inefficient (Marai et al. 2009;Marai and Haeeb 2010). In any case, scarce information is available on the effect of THI in buffalo reared in intensive conditions, like in Italy.
Therefore, the aim of this study was to elucidate the result of the variation of THI on milk yield, quality and SCC in primiparous and pluriparous buffaloes throughout the lactation.

Farm management and animals
Animal welfare approval was not needed for this study as data came from pre-existing databases. Records were provided by one commercial buffalo farm located in Cerignola, Foggia (41.2656 N, 15.8936 E) in the southeastern region of Italy. Originally, data were collected on 808 buffalo cows and included 21,644 records. Data recording was between January 2018 and February 2019. Data from December 2018 was not available because the official milk recording service of the Italian Breeders Association was not operating. Information pertaining to each individual included: animal ID, birth date, calving date, parity order, stage of lactation (SOL), milk yield, fat percentage, protein percentage and somatic cell count. Two sources of data were available: A) monthly milk yields and somatic cell counts provided by the official milk recording service of the Italian Breeders Association, and B) daily milk production recorded directly from the herd milking unit one, three and five days before the official milk recording date, respectively. Then, two sets of data were created. The first dataset (D1) included only monthly milk contents and somatic cell counts from the official milk recording service. The original number of records per individual within lactation ranged from a minimum of 1 to a maximum of 11. The second dataset (D2) was created aggregating daily milk production from both official milk recording (i.e. daily milk yield recorded once per month) and herd milking unit (i.e. daily milk yield recorded one, three and five days before the official milk recording date). In this dataset, the original number of records per individual within lactation ranged from a minimum of 12 to a maximum of 40. Indeed, up to a maximum of 4 monthly records per buffalo cow was available. The rationale behind using two different data sources was to maximise the number of available information within trait.

Traits
Production traits such as milk yield (MY, kg/d), protein percentage (PP), fat percentage (FP) and somatic cell count (SCC) were used for this study. In order to cope with non-normality SCC was log-transformed into somatic cell score (SCS) using the formula proposed by (Ali and Shook 1980): The SOL was expressed as 11 classes based on 30 days in milk (DIM) intervals (class 1 from 4 to 30; class 2 from 31 to 60; class 3 from 61 to 90; class 4 from 91 to 120; class 5 from 121 to 150; class 6 from 151 to 180; class 7 from 181 to 210; class 8 from 211 to 240; class 9 from 241 to 270; class 10 from 271 to 300 and class 11 from 300 to 365). Parity was grouped into 5 classes, where 5þ parity included animals that were in their 5th or greater parity (maximum parity number ¼ 8). Months of calving were grouped into 4 season classes namely Autumn (calving from September to November), Winter (calving from January to February, December was not considered because official recording was not available), Spring (calving from March to May) and Summer (calving from June to August). A minimum of 5 records were required per buffalo cow within parity. Additionally, MY, PP and FP outside the range of mean ± 3 standard deviations (SD) were excluded. The resulting final data set used for analysis consisted of 20,396 test-day records observed on 797 buffalo cows.
Descriptive statistics per each trait by parity class, after editing are in Table 1.

Meteorological-information
Environmental temperature and relative humidity (ET and RH, respectively) were recorded from January 2018 to February 2019 at a weather station located 1 km from the commercial farm.
Temperature-humidity index was calculated according to the equation from other authors (Vitali 2009): where ET is the environmental temperature expressed in degrees Celsius and RH is the relative humidity expressed as a fraction of the unit. The ð1:8 Â ET þ 32Þ term is used for the conversion from degree Celsius to Fahrenheit.
The changes of monthly average values of THI by year of recording (2018, 2019) can be observed in Figure 1.

Statistical analyses
This study was made of two parts, according to the data source used, namely D1 (i.e. official monthly milk recording) or D2 (i.e. official monthly milk recording plus milk recorded using the herd milking unit).
Data from D1 were used to elucidate the relationship between PP, FP, SCS and THI while data from D2 were used to investigate the effect of THI on MY.
When using D1, the following general mixed model (1) was used:  where y ijklmn is the dependent variable at test-day (PP, FP, or SCS), l is the overall mean, YS i is the fixed effect for year-season of calving (i ¼ 1, … ,8), where years of calving were 2018 and 2019 and season of calving was defined as Winter (from January to February), Spring (from March to May), Summer (from June to August) and Autumn (from September to November); ðTHI Â ParÞ jl is the fixed effect for the interaction between THI class (j ¼ 1, … ,48) where THI was 43,44,47,51,52,61,63,69,74 and 75 and parity (l ¼ 1, … ,5), ðSOL Â ParÞ kl is the fixed effect for the interaction between SOL (k ¼ 1, … ,11) and parity (l ¼ 1, … ,5), LR is a linear regression on SCS (if y ¼ PP or FP) or on MY (if y ¼ SCS), anim m is the random effect for the buffalo cow and e ijklmn is the random residual error. When using D2, the following mixed model (2) was used: where y ijklmn is the dependent variable test-day MY, l is the overall mean, YS i is the fixed effect for year-season of calving (i ¼ 1, … ,8) where years of calving were 2018 and 2019 and where season of calving was defined as Winter (from December to February), Spring (from March to May), Summer (from June to August) and Autumn (from September to November); ðTHI Â ParÞ jl is the fixed effect for the interaction between THI class (j ¼ 1, … ,48) where THI was and parity (l ¼ 1, … ,5), ðSOL Â ParÞ kl is the fixed effect for the interaction between SOL (k ¼ 1, … ,10) and parity (l ¼ 1, … ,5), anim m is the random effect for the buffalo cow and e ijklmn is the random residual error.
Data preparation and editing, and all statistical analysis were performed using the R programming environment v.3.6.1 (R Core Team 2018). The R package nlme (Pinheiro et al. 2020) was used to fit mixed models. For all models, hypothesis testing and leastsquare means for fixed effects were performed and calculated using the ANOVA and the lsmeans functions from the stats (R Core Team 2018) and emmeans (Lenth et al. 2021) R packages, respectively. Plots were created using package ggplot2 (Wickham 2016). All differences were considered to be significant when p < .05.

Results
The average observed trend for MY across SOL and parity can be observed in Figure 2. As expected, pluriparous buffalo cows had higher daily milk yields than primiparous until approximately 200 days in milk (DIM). Milk yield at peak was on average at 46 6 9 days and ranged from a minimum of 12:162:6 kg/d for primiparous to a maximum of 15:463:8 kg/d for 3rd-parity buffalo cows.
The average observed trend for FP, PP and SCS across SOL and parity can be observed in Figure  3(A-C), respectively. Fat content increased across lactation and, on average, was higher for primiparous than pluriparous. The trend for protein content was specular to milk yield with a nadir around 80 DIM and highest values after 200 DIM. SCS increased steadily for all parities until 150 DIM and afterwards the rate of increase changed according to parity.
Hypothesis testing for FP, PP and SCS (model 1) are summarised in Table 2.
The interaction between year and season of calving (YS) was significant for FP (p < .05) and extremely significant for SCS (p < .001). Interestingly, no significant effect was observed for PP. THI was significant for all traits and its interaction with parity was also significant for PP and SCS, respectively. Least square means for the effect of THI on FP are in Figure 4.
Results for FP suggest an unfavourable effect of THI (i.e. a decrease in FP) when the latter is above 69 and lower than 47. However, a moderate variability can be observed suggesting that the effect of a variation in THI might be more important than a unique, absolute threshold. Model 1 also included the regression on MY on FP, but its effect was not statistically significant.
THI affected also PP, but an interesting and significant interaction between THI and parity was observed. The relationship between FP and THI across parities can be observed in Figure 5(A).
Overall, the relationship between PP and THI was positive ( Figure 5(B)), even if of different magnitude according to the parity. A common drop when THI is below 62 can be observed, especially in primiparous buffalo cows.
A significant interaction with THI across parities was observed also for SCS (Figure 6(B)), which decreased until THI 69 and then rose as before. A quadratic regression was also fitted, suggesting that the decrease was larger for pupiparous than for primiparous.
The interaction between SOL and parity was also significant for FP and PP but not for SCS. An increasing trend was observed for FP across DIM while PP pattern was specular to MY, with a drop around MY peak and an increase afterwards (results not shown). A significant effect of SCS on FP was also observed.
Hypothesis testing for MY from Model 2 is summarised in Table 3.
All effects included in the models were significant. We observed an unfavourable effect on daily milk yield for buffaloes calving in summer and autumn season of year 2018.
The relationship between THI in class and MY across parity can be observed by plotting their relative least square means ( Figure 6). For clarity, a regression line has been added, quantifying the MY variation per THI unit across parities.
The effect of THI on MY across parities is not definite but overall a favourable relationship can be observed. Indeed, fitting a linear regression between MY least square means and THI, the estimated MY  increase per THI unit ranged from a minimum of 0.022 kg/THI in parity 2 to a maximum of 0.037 kg/THI in parity 1 ( Figure 6). Moreover, even if it was not possible to significantly detect the THI value at which MY starts to decline, our results suggest a negative effect of THI on MY when the former is below 60, i.e. a reverse result if compared to what commonly observed in dairy cattle.
Finally, least square means for the parity effect on MY and their relative 95% confidence interval can be observed in Table 4. Primiparous buffalo cows had lowest and significantly different least square means.

Discussion
Because of its adaptability to hot climate and tropical environmental conditions, a great interest has been paid to buffalo species. In the last years, several studies have focussed on the influence of THI on animal fitness. However, only few researches are available for buffalo species, in particular regarding the sensitivity to heat stress and its effect on animal welfare and milk production (Ahmad and Tariq 2010;Yadav et al. 2016;Purohit et al. 2020;Mishra 2021). Milk yield and composition observed in this trial are in agreement with those reported in previous studies, carried out in both Italian Mediterranean Costa, De Marchi, Battisti et al. 2020) and Murrah (Cer on-Muñoz et al. 2002) buffaloes. MY showed an increase until approximately 50 days postpartum, and a gradual decrease in the following months. The highest MY was recorded in third parity buffaloes, as observed in other studies (Costa, Negrini, De Marchi et al., 2020), while FP and PP mean values showed an opposite trend, with greater values in primiparous. This is probably due to the higher milk production recorded in pluriparous animals, compared to primiparous counterparts. The fat content increased across lactation while the protein content had an opposite trend compared to MY with an initial reduction until 80 DIM and a subsequent increase. Interestingly, SCS increased steadily throughout the lactation and parity: also these findings are in agreement with previous studies carried out either in buffalo (Costa, Negrini, De Marchi et al. 2020) and dairy cattle (Yang et al. 2013;Nasr and El-Tarabany 2017;Gonc¸alves et al. 2018). Reasons for the observed pattern can also be found in the health of the mammary gland of primiparous compared to pluriparous cows (McDonald 1968;Berry and Meaney 2005;McCarthy et al. 2007;Olde Riekerink et al. 2008;Walsh et al. 2007;Gonc¸alves et al. 2018), confirming the sensitivity of the animals to inflammations. The high values of Figure 6. The relationship between daily milk yield (MY) and somatic cell score (SCS) and temperature-humidity index (THI) across parities.  SCS recorded at the end of lactation are probably due to the physical, mechanical and biological stress of the gland, as well as the higher probability of inflammation by germs and other agents.
The main purpose of this study was to evaluate the influence of THI on production and qualitative traits of milk in buffaloes, in latitudes of the southeast of Italy and in intensive breeding conditions. It has been proposed that the most suitable condition for buffaloes is an environmental temperature ranging between 13 and 18 C, together with 55-65% relative humidity and 5-8 km/h wind velocity (Williamson and Payne 1978). However, according to other authors (Crudeli 2011) the thermal comfort zone (or thermo-neutral zone -TNZ), which can be defined as a range of environmental temperature in which homeothermic animals survive without any expenditure of energy to maintain homeostasis of the body (Mishra 2021), was estimated at 21 C. For most farm animals TNZ is recognised between 4 and 25 C (Habeeb 2018) and the potential efficiency of livestock species, including buffaloes, is highest in this condition. If the environmental temperature exceeds TNZ, higher and lower critical temperatures expose the livestock to thermal stress conditions. Indeed, any deviation of the environmental temperature that goes beyond the upper and lower critical temperature may result in heat (Mishra 2021;Bharati et al., 2017) or cold stress (Hoelscher 2001;Mader 2003). In our study, the highest average THI values were recorded between June and September with a maximum in July, whereas the lowest values were recorded between January and February of both recording years (2018 and 2019) and with a minimum in February and January of 2018 and 2019, respectively.
MY recorded in a previous study carried out on Mediterranean buffaloes (Costa, Negrini, De Marchi et al. 2020) showed small variability in buffaloes calved throughout the year, although the highest production was recorded in animals that calved in March compared to their counterparts in other seasons. One possible explanation may be found in the most favourable climatic condition, which has a positive effect on both feeding and management. In our study we observed an unfavourable effect on MY when THI is below 59, suggesting that low temperature and humidity may have a detrimental effect on milk production. On the contrary, when the THI is higher than 60, MY is more stable. These results actually differ from what reported especially in dairy cattle (Hill and Wall 2015;Lambertz et al. 2014;Pragna et al. 2016;Liu et al. 2019;Ekine-Dzivenu et al. 2020). Indeed, in dairy cattle a negative relationship between MY and THI was observed when the latter increases (Nasr and El-Tarabany 2017;Pragna et al. 2016;Ekine-Dzivenu et al. 2020). The observed buffalo ability to cope successfully with higher values of THI (i.e. no decrease related to an increase in THI) might be probably due to both the tropical origin of the species (Minervino et al. 2020;Zhang et al. 2020) and a progressive adaptation after the stress period. As suggested in cattle (Ekine-Dzivenu et al. 2020), when an individual is exposed to prolonged heat load it can acclimatise to the heat stress through the process of acclimatising homeostasis. This process is characterised by a decrease in growth hormone, catecholamine, glucocorticoid, T3 and T4 levels eventually resulting in a reduction in basal metabolic rate and heat production.
Nevertheless, our most interesting finding was the susceptibility of buffaloes to low values of THI. This interesting phenomenon can be explained by some anatomical features of buffalo, such as a hairless and thicker skin and by the tropical origin of the species. In particular, when buffaloes undergo cold stress, they need to utilise energy reserves to maintain thermal homeostasis, leading to the detriment of milk production.
A not regular effect of THI on MY across parities has been observed but overall a favourable relationship can be observed. In particular, it was interesting to observe how the MY in first parity is negatively affected by low temperatures compared to other parities. That's probably due to both the tropical origin of the species, as mentioned above (Minervino et al. 2020;Zhang et al. 2020), and because primiparous are more affected by the negative energy balance post partum and fail to adapt their body temperature and as a result there is a loss of milk yield (Janovick and Drackley 2010;Yehia et al. 2020). To our knowledge, there aren't studies that have observed this, as opposed to others where, instead, there is evidence of greater resistance of the animals of first parity than pluriparous to heat stress (Bernabucci et al. 2015;Maggiolino et al. 2020).
Regarding milk composition, significant changes in FP and PP were observed in relation to THI, especially for THI values higher than 62. This is in disagreement with previous studies carried out in cattle, in which a reduction in fat and protein percentages have been highlighted as a consequence of heat stress (Arieli et al. 2004;Rejeb et al. 2012;Bernabucci et al. 2015;Nasr and El-Tarabany 2017). It cannot be ruled out that the reduction of MY may be responsible for this increase, due to solid concentration even if the fitted model actually included an adjustment for the production level.
A significant effect of the season of calving was also observed on SCS. Higher values were recorded in animals calving from June to November (summerautumn), with an opposite trend to MY. Interestingly, higher SCS values have been recorded at THI values <51 and >69. This finding may suggest an effect of either heat and cold stress on mammary gland health, eventually affecting its susceptibility to mastitis as reported in dairy cattle (Alhussien and Dang 2018). The SCS increase with high THI values agrees with previous results from both cattle (Smith et al. 2013;Nasr and El-Tarabany 2017) and buffaloes (Singh and Ludri 2001). However, some authors (Singh and Ludri 2001) reported a reduction in somatic cells during cold weather. In a recent trial by Lovarelli et al. (2020), dairy cattle maintained at low THI values showed higher lying time, compared to the same animals maintained in thermoneutrality or in heat stress (Lovarelli et al. 2020). Furthermore, a higher lying behaviour, particularly during the post-milking period (Watters et al. 2013), has been associated with high risk of experiencing SCC increase. It cannot be ruled out that a similar behaviour may be present also in buffalo, but its tropical origin might suggest a higher sensitivity to cold stress.

Conclusions
This is the first report which investigated the effect of THI on milk yield and somatic cell score of the Italian Mediterranean buffalo. Previous researches have only suggested an upper limit beyond which buffaloes show heat stress. No information is available so far for the possible stress caused by low THI in this species. According to our results, it can be suggested that an optimal THI range for water buffaloes at our latitudes could be between 59 and 63, although THI values lower than 62 were responsible for some unfavourable effect, particularly in primiparous. Thus, it is worth pointing out that the magnitude of the THI effect largely depends on the trait observed and the parity. Primiparous buffalo cows were more affected than pluriparous, especially for fat or protein content and somatic cell score. Nevertheless, further studies including a larger amount of data, possibly from different farms, are needed to better elucidate the influence of THI on buffalo welfare and production in intensive farming conditions, like in the south of Italy.

Ethical approval
Animal welfare approval was not needed for this study because data came from pre-existing databases.

Disclosure statement
No potential conflict of interest was reported by the authors.

Data availability statement
The data that support the findings of this study are available from the corresponding author [SB], upon reasonable request.