Defining and evaluating heat stress thresholds in different dairy cow production systems

The aim of this study was to assess the impact of heat stress in dairy cows on test-day records for production traits and somatic cell score (SCS) in the state of Lower Saxony, Germany. Three different production systems were defined: A production system characterized by intensive crop production (=indoor housing), a pasture based production system, and a maritime region. Heat stress was assessed by two temperature-humidity indices (THI) modelled as random regression coefficients in an analysis of variance: One (THIBo) was defined as an average of hourly THI, calculated from hourly recorded temperatures and humidities, the other (THIRa) was based on daily maximal temperature and daily minimal humidity. In all production systems, THIBo=60 and THIRa=70 were identified as general thresholds denoting a substantial decline in test-day milk yield. For daily fat and protein percentage, no universally valid thresholds were identified. In contrast for SCS, especially in the maritime region, heat stress as well as cold stress thresholds were found. Regression analysis was used to study the change in test day milk yield in response to THI of those THI ranges with an obvious decline in milk yield. Regression coefficients were −0.08 kg/THIBo and −0.16 kg/THIRa for the crop production system, −0.17 kg/THIBo and −0.23 kg/THIRa for the pasture based system, and −0.26 kg/THIBo and −0.47kg/THIRa for the maritime region. Based on statistical information criteria, identified thresholds for THIBo should be given preference over THIRa when applying genetic studies on heat stress in German Holstein cows.


Introduction
In the course of climate change, it is assumed that also in regions traditionally characterized by less extreme climate conditions, cows will be faced with temperatures beyond their »comfort zone« (IPCC 2007).Such extreme temperatures hamper evaporative loss, especially in high yielding dairy cows (e.g.Kadzere et al. 2002, Collier et al. 2006).Many studies from so called 'hot weather areas' in the world found a decline in milk production and a depression of reproductive performance in dairy cows due to heat stress (e.g.Ravagnolo & Misztal 2002, Bohmanova et al. 2007).
Heat stress is influenced by air temperature, humidity, air movement, solar radiation, and precipitation (Igono et al. 1992).Several indices were developed to describe the impact of heat stress on production and functional traits.One of them, the temperature-humidity index (THI), is based only on temperature and relative humidity, which were weighted differently in different climatic zones (Thom 1959, Yousef 1985).Bohmanova et al. (2007) compared several THI formulae for different regions in the USA based on the decline in daily milk yield.In addition, they identified differences in thresholds for the onset of heat stress for different combinations of regions and indices (Bohmanova et al. 2007).Further possible indices like the comprehensive climate index (CCI, Mader et al. 2010) consist of temperatures adjusted for solar radiation, wind speed, and relative humidity.Hence, such an index covers a wider range of environmental conditions, but index complexity complicates analyses of variances using random regression methodology.Furthermore, solar radiation and wind speed are of minor importance for cows kept in insulated barns, implying that there is additional uncertainty, if detailed information for management and housing conditions of farms is missing.
Consequently, most heat stress studies utilised THI formulae, but only a limited number of these studies have been conducted in Europe.For example in Slovakia in 2003, Broucek et al. (2007) found 80 days with mean values for THI above a defined THI threshold of value 72.This specific year 2003 was characterized by an extremely hot summer in East Central Europe.Reiczigel et al. (2009) showed that the number of heat stress days per year (THI >68) has increased countrywide from 5 to 17 days in Hungary in the past 30 years.From the economic point of view, especially for Germany as one of the main exporters of breeding stock beside France and The Netherlands (Peters & Meyn 2005), heat stress in Holstein dairy cows is of increasing concern.In 2011, the total sales of breeding stock included 72 000 heads, and 63% were exported to countries in the tropics or subtropics characterized by high temperatures, e.g.Algier, Marocco or Tunisia (DHV 2012).Bényei et al. (2001) summarized the problems animals have to overcome when moving from temperate to tropical climates: transport stress, change of ambient temperature, change in seasonal day length patterns, change in diet, and new management routines.A moderate antagonistic relationship between general production and heat tolerance on the genetic scale was verified by Ravagnolo & Misztal (2000).They concluded that heat tolerance will further decrease, if selection strategies only focus on productivity.The approach of improving heat tolerance by crossing high-performance purebreds with locally adapted breeds such as Bos indicus (Zebu) was of limited success, because purebred imported Holsteins nevertheless had a higher production level for milk and for protein yield (McDowell et al. 1996).Worldwide highest milk yield of purebred Holstein cows in Israel is an example for successful adaptation and the overall superiority of the Holstein breed, even when heat stress occurs (Aharoni et al. 1999).
However, additional gain is possible when implementing breeding strategies focussing on heat stress tolerance.The development of breeding strategies implies knowledge about the onset of heat stress and the availability of heat stress thresholds as an essential prerequisite.Consequently, the aim of this study is to evaluate the impact of THI on test-day records for production traits and somatic cell score (SCS) for different regions characterized by different production systems in the state of Lower Saxony, Germany.Bohmanova et al. (2007) impressively showed the necessity for evaluating different heat stress indicators in different regions in the USA, and the same applies when focussing on circumstances in Germany.

Data
Data from Holstein cows from calving year 2009 included test-day records for the traits milk yield, fat percentage, protein percentage and somatic cell score.SCS was calculated by taking the logarithm of somatic cell count (SCC): (1) The data set used for statistical analyses comprised 49,397 test-day records from a region with focus on crop production (=indoor housing), 29 704 test-day records from a pasture based region, and 63 658 test-day records from the maritime region at the coast (Figure 1).These regions were pre-determined for scientific studies as typical model regions of Lower Saxony for all research groups involved in the joint research project KLIFF (climate impact and adaption research in Lower Saxony).The crop production region is characterized by keeping cows constantly indoors, whereas the two other production systems can be described as variants of high input grazing systems.The latter two regions differ in climate conditions, e.g.substantial cooling effects due to sea breeze in the maritime region.Descriptive statistics of all traits for the analysed production systems are given in Table 1.Mean test-day milk yield was highest in the crop production region (30.2 kg), but combined with lowest fat percentages (3.92 %).Somatic cell score was slightly higher in the pasture based region (2.66 vs. 2.63).
Only larger herds comprising at least 50 records per test-day were used for statistical analyses to ensure sufficient size of contemporary groups.Lactation was divided into five stages according to Huth (1995).Stage one (<2 weeks) is characterized by a steep increase of milk yield.It is followed by the peak yield (stage two: 2-11 weeks).In the next stages, milk yield is decreasing linearly at different rates: steep slope (stage three: 11-20 weeks), decelerated (stage four: 20-33 weeks), and again steeper towards the end (stage five: >33 weeks).After eliminating observations outside the 95 % confidence interval of a regression curve modelled with Legendre polynomials of order four within lactation stage, a cow had to fulfil the requirement of at least five test-day records per lactation.Production data was merged with hourly temperature and humidity information recorded three days before the respective test date.Weather data was recorded in five different weather stations.Distances between farms and weather stations were calculated using longitudes and latitudes converted into the Gauss-Krueger coordinate system, and applying the simple Pythagoras triangle.The maximum distance between a farm and a weather station was 60 km.

Heat stress indicators
Heat stress indicators were two different THI formulae, which are generally used when identifying the genetic component of heat stress.One (THI Bo ) was the daily THI as suggested by Bohmanova et al. (2005) for the humid to semiarid climate in the USA.This is the mean of hourly THI calculated with the formula of NRC (1971) using hourly measurements for dry bulb temperature (T) and relative humidity (RH): The other (THI Ra ) based on the formula of NOAA (1979) which was found to fit best to data from public weather stations in Georgia in the study by Ravagnolo & Misztal (2000).This formula used the maximum of the daily temperature (T max ) and the minimum of the daily relative humidity (RH min ):

Statistical models
A linear mixed model including fixed effects and random regressions for different THI, i.e.THI Bo and THI Ra with Legendre polynomials up to order four, was applied to assess the impact of heat stress on test-day records.The procedure Proc Mixed as implemented in SAS 9.2 (SAS Institute, Cary, NC, USA) was used for the statistical analyses.In matrix notation, the statistical model was: where y represents the vector of observed test-day data.Vector b with incidence matrix X included the fixed effects herd-year-season of calving, lactation stage and lactation number, as well as regression coefficients.Vector a includes the random cow effect with the corresponding incidence matrix Z, and e was the vector of random residual effects.Specific heat stress thresholds within production systems were selected by visual inspections, i.e. studying the decline in production level with increasing THI.Regression analysis was used to study the change in test day milk yield and test day SCS in response to THI of those THI ranges where an obvious decline in milk yield was observed.
For the identification of the most appropriate heat stress indicator, i.e. the comparison of THI Bo and THI Ra , residual variances and AIC (Akaike Information Criterion, Akaike 1973) values have been used as evaluation criteria.

Characterization of production systems using heat stress indicators
General descriptive statistics (means and standard deviations) for monthly THI Bo and THI Ra in different production systems are depicted in Figure 2.An interesting finding from the meteorological point of view is that high summer temperatures are compensated by low values for humidities in the maritime region at the coast.Consequently, compared to the other production systems, monthly THI Ra are lowest in the maritime region from March to September.Theses extreme values mostly occurred in the months of July and August.Only 0.1 % of all test-day records are from days with THI Bo >70 in the pasture based production system, and no record corresponds to those extreme THI Bo in the other two production systems.For THI Ra , we consistently found the lowest percentage of records in extreme classes in the maritime region.For example, 30 % of test-day records in the crop production region, 25 % of test-day records in the pasture region, and only 17 % of records from the maritime region correspond to THI Ra >70.Despite differences in absolute values for THI Bo and THI Ra , simple Pearson correlations between both indices within the same production system were throughout higher than 0.97.

Evaluation of THI Bo
The discussion and illustration in Figure 3 is focussed on lactation stage three according to Huth (1995) referring to the relatively long period from 11 to 20 weeks after calving.The general pattern was similar, though less pronounced in the other four lactation stages (not shown).
Least square means from analysis of variances with random regressions on THI Bo revealed a substantial decline in daily milk yield for THI Bo >60.THI Bo =60 corresponds to a temperature of 16 °C, and is relatively independent from humidity.This value could be considered as a general threshold for the upper point of the thermoneutral zone.Different ranges for the thermoneutral zone of cows in different regions are given in the literature: −5 to 25 °C (Hahn 1999), 18-28 °C and 40-60 % RH (Aceves et al. 1987) and 0-16 °C (Yousef 1985).The latter corresponds well to our findings, i.e. daily milk yield increases up to a value of THI Bo =30, then stagnates, and decreases beyond a threshold value of THI Bo >60.
The pasture based region and the region of crop production show similar patterns in terms of decreasing daily milk yield due to heat stress.The identified threshold was THI Bo =62 for both regions, respectively.In the maritime region, THI Bo =60 was identified as a threshold, and in this case, the steepest decrease within a production system was found.Schierenbeck et al. (2011) identified potential contract herds for progeny testing, characterized by superior management within Lower Saxony, and most of them were located in the region corresponding to the pasture based systems of the present study.Accordingly, anticipated effects of heat stress in the pasture based system may be compensated by superior herd management.Superior management in the pasture based system can also explain the higher heat stress threshold compared to the maritime region characterized by higher wind speed and cooling effects from the adjacent sea.Also in the region characterized by crop production, cows have more comfortable husbandry conditions which overcome detrimental heat stress impact.Cows are kept indoors in modern and partly insulated buildings.The explanation for lower susceptibility to heat stress in the crop production system is also valid low values of THI Bo which indicate cold stress.The generally higher production level in the region of crop production does not necessarily correspond to a stronger decrease in daily milk yield due to heat stress.This disagrees with a study by Berman (2005) conducted in Israelian Holstein cows, describing a strong association between production level and decline of daily milk yield due to heat stress.As one specific example, the author found that the threshold temperature decreased by 5 °C when milk production increased from 35 to 45 kg/day.The overall threshold for heat stress of THI Bo =60 as identified for milk yield in the present study differs from those reported in studies from the USA.For example, Bohmanova et al. (2007) used the R-square value from regression analysis, and they reported that heat stress in dairy cows mostly occurs beyond values of THI Bo =70.Similar thresholds were shown by Ravagnolo et al. (2000), who compared values of THI calculated from mean and maximum temperatures.In the US, Dikmen & Hansen (2009) suggested that the difference between their THI Bo threshold of value 78.2 (Florida) and those of Bohmanova et al. (2007), i.e. value 72 in Georgia, and value 74 in Arizona, could be due to better adapted cows, or special housing features.Another reason could be the usage of afternoon measurements instead of daily means.This fact also might have caused that the upper critical temperature for the occurrence of hyperthermia in Florida was 28.4 °C (Dikmen & Hansen 2009), in contrast to measurements of Berman et al. (1985) in Israel (upper critical temperature: 25 to 26 °C).Substantially lower thresholds in our study for THI Bo could be attributed to the reduced adaptability of German Holstein cows to heat stress scenarios.Dikmen et al. (2009) measured vaginal temperatures in dairy cows in intervals of 15 min over the whole day.They explained higher vaginal temperatures of Swedish Red compared to U.S. Holsteins with the development of the Swedish Red in Northern Europe and limited selection under North American heat stress conditions.
In contrast to milk yield, no specific heat stress threshold for fat percentage could be identified.For the interpretation of results, the pronounced genetic antagonism between milk yield and fat percentage should be kept in mind.Fat percentages decreased steadily with increasing THI, without an obvious THI threshold.Surprisingly, there is even a little increase or stagnation in daily fat percentage starting from THI Bo =67 in all regions.Potentially, the substantial reduction in milk yield discussed above leads to higher concentrations of milk contents, and thus unchanged or higher fat percentages during times of heat stress.This finding is most obvious in the maritime region which revealed the strongest decrease in milk yield.In the pasture based system, where cows are more exposed to radiation, and no cooling exists, daily fat percentage remains on a constant level.Protein percentage is continuously decreasing with increasing THI Bo .
Daily SCS shows different patterns for alteration in THI Bo in different production systems.In contrast to the pasture based production system and the maritime region, SCS is lowest for low values of THI Bo in the region characterized by crop production.In the latter region, we observed only a marginal increase in SCS for values of THI Bo >60.Following results from our studies, heat stress in terms of undesired effects on SCS only occurred in the maritime region.In this region, also indications for 'cold stress' were identified.SCS is not only an indicator trait for mastitis (Philipsson et al. 1995), but also reflects a cow's immune response to general stress situations (Coffey et al. 1986).In Arizona, Shathele (2009) observed the highest incidence of clinical mastitis for most of the bacterial isolates in cold months with a peak of coliform bacteria for temperatures below 21 °C after the rainy and muddy periods.In contrast, Elvinger et al. (1991) and Broaddus et al. (2001) found the highest SCC and highest values of infectious bacteria, respectively, when heat stress occurs.Hence, heat stress can have differential impact on SCC which needs further clarification.
In order to derive outcomes per THI Bo unit, linear regression coefficients were calculated within production systems.Generally for this type of analysis, THI values ranging between the identified threshold and the maximal THI value were considered.Accordingly, regression analyses for milk yield were carried out for the THI Bo range of values 62-69, 62-71 and 60-70, for the crop production region, pasture, and coast, respectively.Regression coefficients were −0.08 kg/THI Bo for the crop production region, −0.17 kg/THI Bo for the pasture-based system, and −0.26 kg/THI Bo for the coast.These differences between different regions situated within close proximity and thus having similar macroclimates indicate that cow's milk production is also sensitive to regional climatic factors.Barash et al. (2001) reported that average milk production was reduced by 0.38 kg/°C and protein by −0.01 kg/°C in Israeli Holstein dairy cows.Reiczigel et al. (2009) observed that two of six indices were able to indicate milk production losses (1.5 to 2 litres per cow per day) due to heat stress.In the study by Bohmanova et al. (2007) the different rates of milk decline per THI unit ranged between −0.40 kg and −0.27 kg in Athens (USA), and from −0.59 kg to −0.23 kg in Phoenix (USA).Bouraoui et al. (2002) found a negative correlation (r=−0.76) between daily THI and milk yield.For THI values of 69 and higher, milk yield decreased by 0.41 kg/THI.Lower milk contents were identified for the summer season (milk fat: 3.24 vs. 3.58 %; and milk protein: 2.88 vs. 2.96 %).

Evaluation of THI Ra
THI Ra , which is calculated using maximal daily temperatures and humidities, revealed values for heat stress thresholds beyond THI Ra =70 (Figure 3).This was shown in many other previous heat stress studies (e.g.Ravagnolo & Misztal 2000).
Compared to THI Bo , test day milk yield is more sensitive to extreme values of THI Ra .The threshold indicating cold stress in all production systems, i.e. leading to a decline in daily milk yield, was a value of THI Ra =40.Thresholds for heat stress for daily milk yield were THI Ra =67 in the maritime region, THI Ra =73 for the pasture based system, and THI Ra =74 for the region characterised by crop production.The patterns of daily milk yield by THI Ra are very similar to those modelled in dependency of THI Bo .Differences for the observed traits when comparing the three production systems were highest for high values of THI Ra , and for low values of THI Bo .Daily protein and fat percentage are increasing with increasing THI Ra beyond a threshold of value 70 for the coast, and beyond a threshold of value 75 for pasture and crop production.Again, this is likely due to the antagonistic nature of milk yield and milk contents.Also the shape of curves for SCS modelled in dependency of THI Ra is similar to those modelled for THI Bo .However, when using THI Ra , there is an extreme increase of SCS due to heat stress in the maritime region, and a more pronounced reaction to cold stress in the crop production system.
Also for THI Ra , linear regression coefficients were calculated for observed ranges of heat stress.Intervals were: values of THI Ra ranging between 74 to 84 for the crop production region, values between 73 to 83 for the pasture based system, and values between 67-82 for the maritime region.The regression coefficients for the decline in milk yield per THI Ra unit were −0.16 kg/THI Ra , −0.23 kg/THI Ra and −0.47 kg/THI Ra for the crop production region, pasture and coast, respectively.In a dataset with different regions of the USA, Freitas et al. (2006) found a decline in milk yield varying between −0.15 and −0.36 kg/THI Ra .Ravagnolo et al. (2000) found a decline for test-day production traits beyond a value of THI Ra =72, i.e. −0.2 kg/unit THI Ra , −0.009 kg/unit THI Ra and −0.012 kg/unit THI Ra for milk yield, protein yield, and fat yield, respectively.
In summary, also in Germany, both heat stress indicators, i.e.THI Bo based on average daily values for temperature and humidity, as well as THI Ra , which is based on maximum temperatures and minimal humidities, can be used to identify heat stress thresholds for testday observations.When deciding for one specific THI for on-going studies that focus on the genetic component of heat stress, we recommend to use THI Bo for statistical modelling.For all traits and all analyses, residual components used as evaluation criterion were lower for THI Bo compared to THI Ra .For example for test-day milk yield, differences between residual variances for models firstly run using THI Ra , and secondly run using THI Bo were 1.16 kg 2 , 1.59 kg 2 , and 0.07 kg 2 , for the crop production region, pasture and coast, respectively.Also AIC values were substantially lower (=preferable) for models using THI Bo .
In conclusion, in accordance with countries characterized by high temperatures and humidities, i.e. regions in the USA, South America, or in the tropics, heat stress thresholds have been identified for daily milk yield in different production systems within the region of Lower Saxony in Germany.For both heat stress indicators used in the present study, heat stress thresholds were lower than in other international studies.An explanation could be the lower level of adaption of German Holstein cows to increasing outside temperatures and humidities.Both heat and cold stress had a negative impact on daily SCS.To summarize, heat stress is of increasing relevance for different dairy cattle production systems in Germany, especially for pasture based production systems.The identified thresholds provide the essential pre-requisites for identification of genetic components of heat stress as done by Ravagnolo & Misztal (2000) and Aguilar et al. (2010) in the USA.

Figure 1
Figure 1 Map of production systems as used in the present study

Figure 2
Figure 2 Means and standard deviations for THI Bo and THI Ra by month and production system (white:  crop production, dark grey:  pasture, and light grey:  coast)

Figure 3 3
Figure 3 3 Least square means and belonging standard errors for daily milk yield, fat percentages, protein percentages, and somatic cell score (lactation stage 11-22 weeks) in dependency of THI Bo and THI Ra for different production systems (white:  crop production, dark grey:  pasture, and light grey:  coast)

Table 1
Descriptive statistics for test-day milk kg, fat %, protein %, and SCS in different production systems SCS: somatic cell score, SD: standard deviation

Table 2
Bohmanova 2007 test-day records belonging to classes of THI Bo ≥60, ≥65, and ≥70; and classes of THI Ra ≥70, ≥75, and ≥80Table2shows the low percentages of records for extreme values of THI Bo >60, >65 or >70 compared to studies conducted in North America (e.g.Bohmanova 2007).Also in Table2, percentages for values of THI Ra >70, >75 and >80 are shown.Only 17 % and 4 % of test-day records in Lower Saxony correspond to such extreme THI Bo >60, and THI Bo >65, respectively.