Application of infrared thermal technology to assess the level of heat stress and milk yield reduction of cows in tropical smallholder dairy farms

Panting score (PS) is a common research tool used to assess the physiological state of cows exposed to heat stress, but it is subjective. Infrared temperature (IRT), measured by either infrared thermometers or cameras, may be a more objective and reliable alternative. Very few studies thus far have evaluated the associations be-tween PS, IRT, and milk production. We investigated the applicability of IRT compared with PS as a means of assessing heat stress and milk yield reduction in dairy cows in tropical smallholder dairy farms (SDF). In autumn 2017, SDF located across 4 typical dairy regions of Vietnam were each visited once to collect farm (n = 32) and individual cow data (n = 344). For each SDF, heat load index (HLI) inside the cowsheds, an indicator of environmental heat load calculated from ambient temperature, humidity, and wind speed, was measured. For each cow, PS (0 indicates a cow breathing normally, not panting; 4.5 indicates an extremely heat-stressed cow with excessive panting, tongue fully extended, and excessive drooling), IRT of the cow’s body, single-day energy-corrected milk yield (ECM), body weight, and body condition score were measured. Cow genotype, age, lactation number, and days in milk were recorded. The IRT of the cows’ inner vulval lip (IVuT) were measured with an infrared thermometer; and the IRT of the cows’ vulval surface (OVuT), inner tail base surface (ITBT), ocular area, muzzle, armpit area, paralumbar fossa area, fore udder, rear udder, fore hoof, and hind hoof were also measured with an infrared camera. Multivariate mixed-effects models were used to assess the associations between HLI with PS and IRT, and associations between PS and IRT with ECM while accounting for the effects of other cow variables. All IRT correlated positively with PS (Pearson correlation, r = 0.23–0.50). Each unit increase in HLI was associated with increases of 0.07 units in PS and 0.09 to 0.23°C in IRT. Each degree (°C) increase in IVuT, OVuT, and ITBT was associated with decreases of 0.75, 0.87, and 0.70 kg/cow per day in ECM, respectively, whereas PS and other IRT were not significantly associated with ECM. Thus, all IRT showed potential to assess the heat stress level of cows; and IVuT, OVuT, and ITBT, but not PS and other IRT, showed potential to predict ECM reduction in cows during heat stress. First cross (F1) Holstein Brown Swiss and F1 Holstein Jersey showed lower PS and yielded higher ECM than the third backcross (B3) Holstein Zebu (7/8 Holstein + 1/8 Zebu) and pure Holstein. Thus, F1 Holstein Brown Swiss and F1 Holstein Jersey could be more suitable for tropical SDF than B3 Holstein Zebu and pure Holstein.


INTRODUCTION
In dairy cattle, heat stress is associated with a reduction in feed intake, milk production, and reproduction, which in turn cause economic losses and raise welfare concerns (Kadzere et al., 2002;Hansen, 2007;Polsky and von Keyserlingk, 2017). Climate change is predicted to make heat stress a more severe issue, especially in the tropics (Morton, 2007;Henry et al., 2012;Sejian et al., 2015). Heat stress can be assessed using environment-based indicators such as temperature, humidity, a temperature-humidity index (THI), a heat load index (HLI; Herbut et al., 2018), or animal-based indicators such as panting score (PS), body temperature, respiration rate, heart rate, DMI, or behavior (Gaughan et al., 2002;Galán et al., 2018;Wang et al., 2018;Rashamol et al., 2019). Although environment-based indicators can help anticipate the effects of heat stress on cows, they do not reflect the physiological changes experienced by the cow as animal-based indicators do (Liu et al., 2019).
Among animal-based indicators, PS, rectal temperature, and respiration rate are most commonly used to assess heat stress in cattle (Galán et al., 2018;Liu et al., 2019). Although rectal temperature and respiration rate are general indicators of animal health, PS was specifically developed to assess the heat stress levels of the cow. A PS of 0 indicates normal panting (normal breath, no drooling), and a PS of 4.5 indicates excessive panting (fast breath from the flank, tongue fully extended, excessive drooling, neck extended, and head held down). When environmental THI <68 equivalent or HLI <70, cattle are usually not heat-stressed and have rectal temperatures of 38.5 ± 0.5°C, respiration rates of 20 ± 10 breaths/min, and PS of 0 (Jackson and Cockcroft, 2002;Gaughan et al., 2008;Zimbleman et al., 2009). When THI and HLI are above the heat stress thresholds, PS, rectal temperature, and respiration rate of the cattle increases (Regan and Richardson, 1938;Gaughan et al., 2008). Measuring PS is simple and noninvasive, without the need to restrain animals or buy equipment. Therefore, PS has been widely studied and applied to dairy and beef cattle in both experimental and commercial conditions (Mader et al., 2006;Gaughan et al., 2009;Alfonzo et al., 2016;Unruh et al., 2017). Panting score appears to be suitable for smallholder dairy farms (SDF) where restraining facilities are limited. However, scoring PS is subjective, and its accuracy depends on the experience of the observers. Measurements of rectal temperature and respiration rate, in contrast, are objective. However, these measurements are quite complicated and timeconsuming, and the measurement of rectal temperature is quite invasive because animals usually need to be restrained. These disadvantages limit the application of these methods under commercial conditions. As a result, rectal temperature and respiration rate are often used to assess heat stress during heat stress experiments but are rarely used on commercial farms.
Infrared thermometers and infrared thermal cameras, which both apply infrared thermal technology, can be alternative devices to the traditional rectal thermometer to measure the temperature of cattle. Both are noninvasive devices that can measure the infrared temperature (IRT) of the surface of cows' bodies with little or no direct contact with the animals (Hoffmann et al., 2013;Tattersall, 2016). Compared with PS, IRT is much more objective. When comparing infrared thermometer and infrared thermal camera, the infrared thermometer is less expensive and easier to use, but users need to stand close to or maintain contact with animals. However, infrared thermal cameras are expensive, require expertise to use, and require time to process the captured images (Unruh et al., 2017). The advantages of infrared thermal cameras are their capacity to record moving animals, their ability to measure temperatures in many locations simultaneously, and the potential to capture images automatically (Hoffmann et al., 2013). A few studies have shown the potential application of infrared thermal cameras for recognizing heat stress in beef heifers (Unruh et al., 2017), Nellore cattle , purebred lactating Holstein cows, Girolando cows (Holstein × Gir crossbreeds; Daltro et al., 2017), and Jersey heifers (Salles et al., 2016). These studies have reported positive correlations between IRT of the animals' body surface and PS, physiological parameters of the cows, and THI of the environment. However, these studies were research projects, in which confounding variables were managed. To our knowledge, IRT technology has not been studied within a production system in SDF, where breeds, management, and microclimate are diverse. Also, studies assessing the associations of PS and IRT with cows' milk production remain very scarce.
Therefore, this study aimed to investigate the possibility of using PS and IRT as tools for assessing the heat stress status and productivity response of the cows in tropical SDF while accounting for variations in breeds of cows. To do this, we investigated the relationships between changes in the shed microclimate as measured by HLI, PS (visual assessment) and body temperatures (obtained via infrared thermal devices), and cow milk production measured by ECM yield. It was hypothesized that, within an SDF, (1) increased HLI will be associated with increased PS and IRT of the cows, and (2) increased PS and IRT, in turn, will be associated with decreased milk production.

MATERIALS AND METHODS
This study was approved by the University of Queensland Human Research Ethics Committee, approval number 2016001815 (December 13, 2016); and the University of Queensland Animal Ethics Unit, approval numbers ANRFA/SVS/565/16/VIETNAM (January 13, 2017) and SVS/010/18/VIETNAM (March 26, 2018). Informed consent was obtained from all subjects involved in the study.

Farm Visit
This study was conducted on 32 SDF located in 4 typical dairy regions of Vietnam: a southern lowland (SL), a southern highland (SH), a northern lowland (NL), and a northern highland region (NH; Bang et al., 2021b). These 4 regions were selected to represent Vietnamese smallholder dairy regions because of their Bang et al.: INFRARED ASSESSMENT OF HEAT STRESS ON TROPICAL FARMS high population of dairy cows and the diverse weather conditions between the provinces. During a period between August 24 and October 7, 2017, each SDF was visited for an afternoon and again the next morning. At each visit, weather data and individual lactating cow data (n = 344) were obtained.
Infrared temperatures of the cows were obtained within this study, and cowshed microclimate data and other individual cow data, including PS, milk yield, age, DIM, lactation, BW, and BCS, were obtained from previous studies conducted concurrently on the same SDF (Bang et al., 2021a(Bang et al., ,b,c, 2022. All measurements of the cows in a given SDF and the microclimate data within the cowshed of that SDF were gathered on the same date.

Microclimate Data
Microclimate data were obtained from Bang et al. (2021c). Briefly, the microclimate parameters inside the cowshed of each SDF were measured at 0600, 0800, 1000, 1100, 1400, 1600, and 1800 h by holding a Kestrel 5400 Heat Stress Tracker (Nielsen-Kellerman) at about 1.8 m above the floor, and standing as nearly as possible in the middle of the cowshed. The microclimate parameters measured by the device included wind speed (WS, m/s; accuracy 3% reading value, resolution 0.1 m/s), ambient or dry bulb temperature (AT, °C; accuracy 0.5°C, resolution 0.1°C), relative humidity (RH, %; accuracy 2% RH, resolution 0.1% RH), black globe temperature (GT, °C; accuracy 1.4°C, resolution 0.1°C), and dew point temperature (°C; accuracy 1.9°C, resolution 0.1°C).
The THI (units) was calculated using the equation of Yousef (1985): where AT = dry bulb temperature (°C) and Tdp = dew point temperature (°C).
The THI and HLI measurements were then averaged to obtain average daytime THI and HLI at each SDF.

Individual Cow Data
All surveyed SDF had either freestall or tiestall cowsheds, where cows were kept inside the cowsheds all days. Therefore, all cows were under shade when all individual cow data were measured.
Infrared Temperature Measurements. To facilitate easy identification and the recording of the data of individual cows, cow identification was written with a large black marker pen on size A5 paper, and this was attached with nontoxic, water-soluble Thien Long paper glue (Thien Long Group Corporation) to the left side of each cow (Figure 1 a).
The average temperature (°C) of the inner vulval lip (IVuT) of each cow was the average of morning and afternoon temperatures of the inner vulval lip, which were measured between 0600 and 0700 h and between 1500 and 1600 h, respectively, using a Combi Infrared Thermometer (Combi Corporation). The inner vulva was chosen to measure temperature because the tail usually obscures this area, so it is usually clean, and the temperature in this area is less affected by environmental variables such as direct sunlight or wind. When obtaining the temperature, the thermometer probe was placed between the vulval lips and pointed to the vulval wall at a position of less than 1 cm from outside the vulval lips, and the measurements were triplicated (Figure 1 c).
Infrared thermal images of all cows in each SDF were taken in an afternoon only, and it often took 2 to 5 min to capture all necessary thermal images of each cow. An infrared camera (FLIR E50 IR camera, FLIR Systems Inc.) was used to capture infrared thermal images of the cows at 10 targeted body areas ( Figure  1 d-l). The specifications of the FLIR E50 were as follows: resolution of 240 × 180-pixels, thermal sensitivity of 0.05°C, accuracy of ±2% of reading, and temperature range of −20 to 650°C. The cow images were taken from a distance of 1.0 to 1.5 m (Montanholi et al., 2015) and within 20 min after afternoon milking (usually between 1530 and 1730 h) when cows were usually standing and eating in shaded areas inside the cowsheds. The reflected apparent temperature was set at 20°C, and the emissivity was set at 0.98, which is the emissivity of human and mammal skins (Kastberger and Stachl, 2003;Hoffmann et al., 2013;Sathiyabarathi et al., 2016). The RH and environmental temperature were obtained using a Kestrel 5400 Heat Stress Tracker (Nielsen-Kellerman) when the first thermal image was captured. The captured images were then used to determine the IRT (°C) of the outer vulval surface (OVuT), rear udder (RUdT), inner tail base surface (ITBT), ocular area (EyeT), muzzle, armpit, paralumbar fossa area, fore udder (FUdT), fore hoof, and hind hoof (HHoT; Figure 1 d-l). When analyzing the images, a circle was drawn around the eyes, and rectangles were drawn around other regions. Although maximum, minimum, and average IRT were obtained from each shape, only the maximum IRT data were used in the subsequent statistical analysis.
Panting Score. The PS data of the cows were obtained from Bang et al. (2021b). Briefly, the day PS was the average of morning and afternoon PS, which were scored between 0500 and 0600 h and between 1400 and 1500 h, respectively. A PS scale from 0 to 4.5 was used (0 indicates a cow breathing normally, not panting; 4.5 indicates excessive panting with fast breath from the flank, tongue fully extended, excessive drooling, neck extended, and head held down; Gaughan et al., 2009;Figure 1 b).
Milk Yield and Other General Cow Data. Milk yield, age, DIM, lactation number, BW, and BCS of each cow were obtained from Bang et al. (2021b). Briefly, the single-day milk yield of each cow was obtained by adding morning and afternoon milk yields obtained directly after each milking. The morning and afternoon milk outputs of each cow were also sampled to analyze milk fat and protein concentrations at the Nutrition Laboratory, Vietnam National University of Agriculture (Hanoi, Vietnam). Single-day milk yield, fat concentrations, and protein concentrations were used to calculate single-day ECM (3,138 kJ/kg of ECM) using the equation of Tyrrell and Reid (1965).
The heart girth of each cow was measured using a tape measure and was used to estimate cow BW using the equation of Goopy et al. (2018). Cows' BCS were obtained by the methods described in Bang et al. (2021b). Body condition score (5-point scale: 1 = lean, 5 = obese) was determined independently by 2 team members and averaged, using the guidelines of Edmonson et al. (1989). The means (±SD) obtained for BW and BCS were 479 ± 74 kg and 2.83 ± 0.50, respectively.
Age, number of lactations, and DIM of each cow were obtained by checking farmers' record books. When farmers did not have record books or did not have complete data in the record book, the missing data were , and measurements of infrared temperature (c-l) of smallholder farm dairy cows in Vietnam. When analyzing thermal images with the FLIR infrared camera (FLIR E50 IR camera, FLIR Systems Inc.), within each shape, the solid red triangle shows the hottest point (maximum temperature), and the solid blue triangle shows the coolest point (minimum temperature). obtained by asking farmers directly. The mean (±SD) obtained for age (yr), number of lactations, and DIM of those cows were 4.5 ± 1.7 yr, 2.3 ± 1.4 lactations, and 189 ± 121 d, respectively.

Statistical Analyses
All statistics were calculated using the base and additional packages of R software (R Core Team, 2018). Cows were the experimental units in all analyses. All individual cow data were measured at the animal level, and THI and HLI were measured at the farm level. Thus, when analyzing the data using cows as the experimental units, the cows within an SDF were assigned the same HLI and THI values as measured in that SDF.
Descriptive Statistics. Descriptive statistics for quantitative variables were calculated for each region using the R package psych (Revelle, 2019).

Statistical Comparisons and Correlations.
A normality test showed that all quantitative variables in this study were normally distributed; thus, one-way ANOVA tests followed by Tukey-Kramer tests (P < 0.05) using the R package agricolae (de Mendiburu, 2019) were applied to compare the means of variables between regions. Pearson correlation coefficients and corresponding significance levels between variables were calculated using the R package psycho (Makowski, 2018).
Multivariate Regression Analysis. A total of 24 multivariate linear mixed effect models were fitted using a restricted maximum likelihood approach, to evaluate the effects of the independent variables on the dependent variables, using the R package lmerTest (Kuznetsova et al., 2017). Farm identification was fitted as a random effect in all models to account for possible clustering. Four farms with fewer than 5 lactating cows per farm (a total of 12 cows) were removed to ensure the quality of the model. The matrix notation describing the model was as follows: where y was the vector of dependent variables, β was the vector of fixed effects, c was the vector of the random farm effect 28 farms; c N I c ∼ 0 2 , , σ ( )       e was the vector of residual random farm effects e N I e ∼ 0 2 , , σ ( )       and X and W were the incidence matrices of the fixed effects and random environmental effects, respectively. Random farm effects and residual effects were assumed to be independently distributed. The parameters of the model σ c 2 and σ e 2 were estimated by the restricted maximum likelihood method.
The first 12 models were used to evaluate the effect of environmental heat load measured by HLI on PS and 11 IRT measurements. In these models, HLI, region, breed, age, lactation, DIM, BW, BCS, and ECM of the cows were fitted as fixed effects (β), and each of the PS and 11 IRT measurements was an outcome variable (y). Although both THI and HLI were indicators of environmental heat load, the calculation of THI included only AT and RH, whereas the calculation of HLI included AT, RH, and WS (Yousef, 1985;Gaughan et al., 2008). Thus, we have chosen HLI rather than THI to be included as a fixed effect in building models.
The other 12 models were used to evaluate the effects of PS and 11 IRT measurements on ECM. In these models, ECM was fitted as an outcome variable (y), whereas region (4 regions), cow breed, age, lactation number, DIM, BW, BCS, and each of the PS and IRT measurements were fitted as fixed effects (β).
When fitting models, Satterthwaite's approximation method was used to calculate degrees of freedom and P-values for both F-statistics and t-statistics. To avoid multicollinearity, the independent variables with the variance inflation factor >5 were removed first (Kock and Lynn, 2012). Based on this approach, the variable "age" was removed from all models. Then, the manual backward elimination technique was applied to remove one by one the variables with F-statistic P-values > 0.1. To choose the final models, models were compared using likelihood ratio tests for the significance of dropped independent variables. Distributions of standard residuals of the final models were plotted to confirm that normality and homoscedasticity assumptions were met (Barberg et al., 2007). The results of the final chosen models are presented. Marginal R 2 of the final models, indicating the percentage of the total model variance explained by the fixed effects, was calculated as the ratio of fixed-effect variance to the total variance, which is the sum of fixed-effect, random-effect, and residual variances (Bartoń, 2019). Conditional R 2 of the final models, indicating the percentage of the model vari-ance explained by both fixed and random effects, was calculated as the ratio of the sum of the fixed-and random-effect variances to the total variance (Bartoń, 2019).
After choosing the final models, whenever the effects of the categorical variables, either breed or region, were significant, the least squares means of different levels of the categorical variables were compared.

Correlations Between Main Variables
The THI, HLI, ECM, PS, and IRT measurements of cows in 4 regions are presented in Table 1. The overall means (±SE) for THI and HLI were 79.4 ± 1.9 and 86.4 ± 3.3, respectively. The means of THI and HLI in SL and NL were significantly higher than those in SH and NH (P < 0.05). The overall mean (±SE) of ECM was 15.7 kg/cow per day. The overall mean (±SE) of PS was 1.3 ± 0.2. Except for the means of EyeT, RUdT, ITBT, and OVuT, which were similar across the 4 regions, the means of PS and all other IRT measurements were highest in SL, followed by NL and NH, and lowest in SH (P < 0.03).

Associations of Heat Load Index with Inner Infrared Vulva Temperature, Outer Infrared Vulva Temperature, and Panting Score
Results of multivariate linear mixed effect models showed that PS and all IRT measurements except  HHoT had significant and positive associations with HLI (P ≤ 0.008; Table 3, Appendix Table A1). Each unit increase in HLI was associated with an increase of 0.07 unit in PS and increases from 0.09°C to 0.23°C in IRT measurements except HHoT when the other predictors were held fixed (Table 3, Appendix  Table A1). We found a suggestive association of HHoT with HLI (P = 0.070; Appendix Table A1). The breed of the cows showed significant associations only with PS (P < 0.001) but not with any IRT (P > 0.1; Table  3).
The changes in confidence intervals (95% CI) of PS, IVuT, OVuT, and ITBT over the range of HLI after adjusting for the effects of all other factors in the models are presented in Figure 2. When HLI changed from the lowest value (74.03) to the highest value (97.27), the changing patterns of IRT measurements were similar to the changing pattern of PS ( Figure 2).
It is necessary to consider the thresholds to interpret the meaning of IRT measurements when assessing the heat stress status of the cows. Gaughan et al. (2008) suggested that, based on PS, cows can be classified into the following categories: normal (PS <0.4), slightly heat-stressed (PS = 0.4 to <0.8), moderately heatstressed (PS = 0.8 to <1.2), and highly heat-stressed (PS ≥1.2). Using the PS thresholds suggested by Gaughan et al. (2008) and the prediction model for PS in Table 3, the thresholds for HLI can be calculated. The calculated thresholds for HLI then can be used to estimate thresholds for IRT using other prediction models for IRT, as in Table 3 and Appendix Table  A1. For example, when using the prediction model for PS (first column of Table 3) and holding other predictors of the model fixed, a PS value of 0.4 (the first threshold of Gaughan et al., 2008) gave an HLI threshold of (0.4 + 4.71)/0.07 = 73 units. Using the prediction model for IVuT (second column of Table  3) and holding other predictors of the model fixed, the HLI threshold of 73 gave an IVuT threshold of 28.93 + 0.10 × 73 = 36.2°C. Other thresholds for HLI and IRT were calculated using the same methods, and Table 4 presents the thresholds calculated for some selected measurements. When using a measurement, the heat stress level of a cow can be categorized as follows: normal when that measurement < threshold 1, slightly heat-stressed when threshold 1 ≤ measurement < threshold 2, moderately heat-stressed when threshold 2 ≤ measurement < threshold 3, and highly heat-stressed when measurement ≥ threshold 3. For example, using the thresholds for IVuT, SDF cows can be considered normal when IVuT <36.2°C, slightly heat-stressed when 36.2°C ≤ IVuT <36.8°C, moderately heat-stressed when 36.8°C ≤ IVuT <37.4°C, and highly heat-stressed when IVuT ≥37.4°C.  Table 1. *Significant at P ≤ 0.05; **significant at P ≤ 0.01; ***significant at P ≤ 0.001.
The least squares means of PS, IVuT, OVuT, and ITBT across regions after adjusting for the effects of all other factors in the models are presented in Figure 3. The corresponding results for other IRT measurements are not presented because they showed the same patterns as the presented IRT measurements. Geographical regions showed significant associations with all PS, IVuT, OVuT, and ITBT ( Figure 3). Consistently, least squares means of PS, IVuT, OVuT, and ITBT were highest in NH and lowest in NL (P < 0.05). The least squares means of PS, IVuT, OVuT, and ITBT in SL were all similar to those in SH, although the pattern of PS was different from the patterns of IRT measurements.
The least squares means of PS for the crossbreeds are shown in Figure 4. In all crossbreed pairs (HOL_BSW, HOL_JER, and HOL_ZEB), an overall trend was that F1HOL crossbreeds tended to have the lowest PS, and the greater the percentage of HOL, the greater the PS. Comparing crossbreeds, the least squared mean of PS of F1HOL_BSW (0.81), B1HOL_BSW (1.08), F1HOL_JER (1.11), and B1HOL_ZEB (1.12) were lower than those of B3HOL_ZEB (1.31) and pure HOL (1.33) (P < 0.05).

Associations of Panting Score, Inner Vulva Infrared Temperature, Outer Vulva Infrared Temperature, and Inner Tail Base Infrared Temperature with Energy-Corrected Milk
Results of multivariate linear mixed effect models showed that only IVuT, OVuT, and ITBT had significant associations (P ≤ 0.035) with ECM, whereas PS and all other IRT measurements did not ( Table 5). The effect patterns of IVuT, OVuT, and ITBT on ECM were very similar to each other. Each degree increase in IVuT, OVuT, and ITBT was significantly associated with decreases of 0.75, 0.87, and 0.70 kg of ECM/cow per day, respectively. The regression coefficients of the effects of PS and other IRT measurements on ECM were positive for PS (0.32, P = 0.593), close to zero for HHoT (0.00, P = 0.978) and paralumbar fossa area temperature (0.07, P = 0.783), and negative for the rest of the IRT (ranging from −0.45 for EyeT to −0.13  for fore hoof temperature, P > 0.113). However, these associations were not significant, and therefore are not presented. The effects of DIM, BW, region, and breed on ECM were consistent across the models (Table 5). Comparing models, the model including the variable OVuT had the highest explanatory power (conditional R 2 = 51%). Therefore, this model was used to draw the effect plots that compare the least squares of ECM across the regions and between crossbreeds ( Figure 5). After accounting for effects of OVuT, DIM, BW, and breed, cows in NH yielded highest ECM (19.7 kg/cow per day), followed by cows in NL (16.8 kg/cow per day) and SH (15.3 kg/cow per day), and cows in SL yielded lowest ECM (14.7 kg/cow per day; P < 0.05; Figure  5a).
When comparing breeds (Figure 5b), the overall trend was that the least squares means ECM of F1 crossbreeds tended to be higher than that of B1 crossbreeds. We found that HOL_BSW tended to yield the highest ECM, followed by HOL_JER, whereas HOL_ZEB and pure HOL tended to yield the least ECM. After accounting for effects of OVuT, DIM, BW and region, the least squares means ECM of F1HOL_BSW (20.9 kg/cow per day) and F1HOL_JER (17.6 kg/cow per day) were highest and significantly higher than those of B1HOL_ZEB (14.6 kg/cow per day) and B3HOL_ZEB (15.2 kg/cow per day; P < 0.05). The least squares  means ECM of F1HOL_BSW was also higher than that of pure HOL (15.7 kg/cow per day).

Pearson Correlations Between Temperature-Humidity Index, Panting Score, Infrared Temperature Measurements, and Energy-Corrected Milk Yield
In the current study, significant negative associations between PS and all IRT measurements with HLI and THI were found, as well as significant negative correlations of ECM with THI, HLI, PS, OVuT, EyeT, and FUdT (P < 0.05). These negative correlations were expected because the mean HLI (86.4 units) and mean THI (79.4 units) between 0600 h and 1800 h in the current study were higher than the heat stress thresholds of THI (68 units; Zimbleman et al., 2009) and HLI (70 units;Gaughan et al., 2008). When heat stress thresholds are exceeded, the body temperature and respiration rate of cattle increases, feed intake decreases, and milk yields decrease (Regan and Richardson, 1938;Könyves et al., 2017;Habeeb et al., 2018;Liu et al., 2019). Moreover, in the current study, all IRT measurements were significantly and positively correlated with PS, which is a heat stress indicator commonly applied in monitoring heat stress in dairy and beef cattle (Mader et al., 2006;Gaughan et al., 2009;Alfonzo et al., 2016;Unruh et al., 2017).
The positive correlations between PS and IRT measurements with HLI and THI in the current study were also consistent with previous studies (Bouraoui et al., 2002;Gaughan et al., 2008;Daltro et al., 2017;Unruh et al., 2017). For example, a study in Brazil using pure HOL and crossbreeds of HOL and Gir cows showed that EyeT, FUdT, and RUdT in the afternoon correlated significantly with PS (0.37, 0.54, and 0.41, respectively; P < 0.01; Daltro et al., 2017). The same authors also reported positive correlations between EyeT, FUdT, and RUdT with THI, but these correlations were not significant (P > 0.01; Daltro et al., 2017). In a study by Bouraoui et al. (2002), THI was positively correlated with respiration rate (r = 0.89), heart rate (r = 0.88), rectal temperature (r = 0.85), and cortisol (r = 0.31). Gaughan et al. (2008) reported positive correlations between HLI and PS (r = 0.93, P < 0.001) and between THI and PS (r = 0.61, P < 0.001). The correlation between HLI and PS (r = 0.56) and between THI and PS (r = 0.54) in the current study were lower than those reported by Gaughan et al. (2008). These differences could be because Gaughan et al. (2008) examined the variation of PS within cows over a wide range of HLI and THI, whereas the current study examined it over a more restricted range. The  . Plot comparing least squares means (black dots) with confidence interval (red arrows) and range (blue shaded areas) of panting score (PS) between different dairy cow crossbreeds on smallholder dairy farms in Vietnam. H, Z, J, and B stand for Holstein, Zebu, Jersey, and Brown Swiss, respectively; and F1, B1, B2, B3, and PHF indicate that the proportion of Holstein in a cow was >1/4 to < 3/4, ≥3/4 to <7/8, ≥7/8 to <15/16, ≥15/16 to <31/32, and ≥31/32, respectively. Least squares means with different letters were significantly different from each other (P < 0.05). negative correlations between ECM with THI and HLI in the current study are in the same direction as the results of previous studies, which reported that each unit increased in THI was associated with a decline in milk yield of 0.08 to 0.41kg/cow per day (du Preez et al., 1990;Ravagnolo et al., 2000;Bouraoui et al., 2002;Könyves et al., 2017).

Applicability of Panting Score and Infrared Technology in Monitoring Heat Stress
Statistical models that describe the associations between HLI with PS and IRT measurements show that each unit increase in HLI is associated with an increase of 0.07 unit in PS and increases ranging from 0.09°C to 0.23°C in IRT measurements. The positive regression coefficient of PS on HLI found in the current study (0.07) is consistent with and slightly higher than the result of Veissier et al. (2018), who reported that within the HLI range of 50.8 to 78.8, each unit increase in HLI was associated with a 0.05-unit increase in PS in grazing dairy cows. The significant positive associations between all IRT measurements and HLI and the significant correlations between IRT measurements and PS, as discussed earlier, suggest that IRT measurements at 11 different positions on the cow's body, determined by using either an infrared thermometer or an infrared thermal camera, could potentially be used to assess the heat stress level of the cows in SDF. Table 4 was constructed to serve as a guideline for interpreting the heat stress level of cows based on HLI and IRT measurements. In Table 4, the thresholds for HLI (73.0, 78.7, and 84.4) were calculated from the thresholds for PS suggested by Gaughan et al. (2008). Animal-based indicators of heat stress were quite close to the HLI thresholds of 70, 77, and 86 suggested by the same authors. This showed the consistency of HLI as an environment-based indicator of heat stress for cattle in different production conditions. However, it should be noted that the thresholds for IRT, as presented in Table 4, were calculated from the thresholds for HLI, which is an environment-based indicator of heat stress, not an animal-based indicator of heat stress. Thus, further studies using animal-based indicators of heat stress such as rectal temperature, respiration rate, or biological markers are recommended to validate the thresholds for IRT identified in the current study.
Results of models describing associations between PS and IRT measurements with ECM indicated that some IRT measurements, including IVuT, OVuT, and ITBT, had significant associations with ECM, but PS did not. The reason for PS not having significant associations with ECM might be because PS is based on respiratory rate, deepness of panting, and degree of drooling (Gaughan et al., 2008), which are all among the adaptive responses of cattle in dealing with high environmental heat load and maintaining their body temperature within a small range. By contrast, IRT measures the heat load retained on the cow's body surface after all adaptive responses are activated.
Data from the current study suggest that IVuT, OVuT, and ITBT are good indicators to assess the heat stress status of the dairy cows in SDF, and they may also be potential traits to use in the selection of heattolerant dairy cattle. This is because IVuT, OVuT, and ITBT were positively associated with HLI and negatively associated with ECM. Desirable cows for selection could be the cows that can maintain low IRT in a hot environment. However, further randomized experiments are necessary to validate the associations between these IRT measurements and ECM. Experiments are also required to estimate levels of genetic variation in these IRT measurements and their genetic relationships with other important production and welfare Bang et al.: INFRARED ASSESSMENT OF HEAT STRESS ON TROPICAL FARMS Figure 5. Effect plots comparing least squares means (black dots) with confidence interval (red arrows) and range (blue shaded areas) of ECM across regions (a) and between crossbreeds (b) of smallholder farm dairy cows in Vietnam. Regions: SL = southern lowland, SH = southern highland, NL = northern lowland, NH = northern highland. Breeds: H, Z, J, and B stand for Holstein, Zebu, Jersey, and Brown Swiss, respectively; and F1, B1, B2, B3, and PHF indicate that the proportion of Holstein in a cow was >1/4 to < 3/4, ≥3/4 to <7/8, ≥7/8 to <15/16, ≥15/16 to <31/32, and ≥31/32, respectively. Least squares means with different letters were significantly different from each other (P < 0.05).

Suitable Dairy Crossbreeds for Tropical Smallholder Dairy Farms
The overall trend shown in Figure 4 suggests that higher PS was associated with a higher percentage of HOL in crossbreeds, which is consistent with the results of other authors (Alfonzo et al., 2016;Dalcin et al., 2016;Daltro et al., 2017). For example, in the same climatic conditions (THI between 74.2 and 84.6 units), heart rate (46.6 beats/min), rectal temperature (39.0°C), respiration rate (52.9 breaths/min), and PS (0.38) of F1HOL_GIR were within the normal ranges of 40 to 60 beats/min for heart rate, 38.3 to 39.3°C for rectal temperature, 23 to 40 breaths/min for respiration rate, and 0 to 0.4 for PS, whereas these physiological parameters of B1HOL_GIR and pure HOL were all higher than the upper threshold of the normal ranges and notably higher than the corresponding measurements of F1HOL_GIR (Alfonzo et al., 2016). Similarly, a study by Daltro et al. (2017) showed that at an ambient temperature range of 20.7°C to 37.9°C and RH nearly 95%, the rectal temperature (40.84°C), respiratory frequency (111.36 breaths/min), and heart rate (99.22 beats/min) of pure HOL were higher than those of F1HOL_GIR and B1HOL_GIR. The lower PS of F1HOL_BSW, B1HOL_BSW, F1HOL_JER, and B1HOL_ZEB compared with those of B3HOL_ZEB and pure HOL indicated that, in terms of heat stress measured by PS, F1 and B1 crossbreeds of HOL, especially with BSW, would be more heat-tolerant of the tropical conditions of Vietnam than B3 or pure HOL.
For milk production, the trend that F1 crossbreeds of HOL with other breeds tended to yield higher ECM than F1 crossbreeds of HOL and pure HOL might reflect the effect of heterosis, which was also observed in other studies (Tadesse and Dessie, 2003;Madalena et al., 2012;Coffey et al., 2016). In addition, the tendency of HOL_BSW crossbreeds and, to a lesser extent, HOL_JER crossbreeds to yield higher ECM compared with HOL_ZEB crossbreeds and pure HOL, as observed in the current study, are in agreement with other studies, which showed that in terms of milk production BSW, to a lesser extent JER, and their crosses with HOL, were more tolerant of tropical conditions than pure HOL (Johnson and Vanjonack, 1976;West et al., 2003;El-Tarabany et al., 2017). Johnson and Vanjonack (1976) reported that at an ambient temperature of 34°C and RH of 80%, the milk yields of HOL, BSW, and JER cows were 41%, 56%, and 76%, respectively, of their normal milk yields at an ambient temperature of 24°C and RH of 38%. Similarly, a study in a subtropical region of Egypt by El-Tarabany et al. (2017) reported that when THI changed from low to high levels, the reduction in daily milk yield and peak milk yield of pure BSW and F1HOL_BSW was much smaller than that of pure HOL. The higher ECM of F1HOL_BSW compared with B1HOL_ZEB, B3HOL_ZEB, and pure HOL, and the higher ECM of F1HOL_JER compared with B1HOL_ZEB and B3HOL_ZEB, suggested that, in terms of milk production, F1HOL_BSW and F1HOL_JER would be more suitable for Vietnamese SDF than B3HOL_ZEB and pure HOL. However, it should be noted that, in the current study, the numbers of F1HOL_BSW (4 cows) and F1HOL_JER (20 cows) were small. Therefore, the results should be interpreted with care, and further controlled experiments are recommended to validate the results of the current study.

CONCLUSIONS
In the context of SDF, all IRT showed potential to predict the level of heat stress in cows. Three IRT measurement sites, IVuT, OVuT, and ITBT, exceeded PS as a predictor of ECM reduction in cows during heat stress. Based on these IRT measures of heat stress, crossbreed cows, F1HOL_BSW and F1HOL_JER, could be more suitable for Vietnamese SDF than pure HOL and B3HOL_ZEB cows.  Variable cow age was excluded from all models due to variance inflation factor >5. The independent variables that were excluded from each final model due to lack of significant effect (P > 0.1) were ECM yield (kg) and the variables with a dash in the P column.