Bunching behaviour in housed dairy cows at higher ambient temperatures

Bunching behavior in cattle may occur for several reasons including enabling social interactions, a response to stress or danger, or due to shared interest in resources such as feeding or watering areas. There is evidence in pasture grazed cattle that bunching may occur more frequently at higher ambient temperatures, possibly due to sharing of fly-load or to seek shade from the direct sun under heat stress conditions. Here we demonstrate how bunching behavior is associated with higher ambient temperatures in a barn-housed UK dairy herd. A real-time local positioning system (RTLS) was used, as part of a precision livestock farming (PLF) approach, to track the spatial position and activity of a commercial dairy herd (c100 cows) in a freestall barn continuously at high temporal resolution for 4 mo between August and November 2014. Bunching was determined using 4 different spatial measures determined on an hourly basis: herd full and core range size, mean herd inter-cow distance (ICD), and mean herd nearest neighbor distance (NND). For hourly mean ambient temperatures above 20°C, the herd showed higher bunching behavior with increasing ambient temperature (i.e., reduced full and core range size, ICD, and NND). Aggregated space-use intensity was found to positively correlate with localized variations in temperature across the barn (as measured by animal mounted sensors), but the level of correlation decreased at higher ambient barn temperatures. Bunch-ing behavior may increase localized temperatures experienced by individuals and hence may be a maladaptive behavioral response in housed dairy cattle, which are known to suffer heat stress at higher temperatures. Our study is the first to use high-resolution positional data to provide evidence of associations between bunching behavior and higher ambient temperatures for a barn-housed dairy herd in a temperate region (UK). Further studies are needed to explore the exact mechanisms for this response to inform both welfare and production management.


INTRODUCTION
Climate change is a pressing global issue which is increasingly concerning to farmers as changes in environmental conditions can directly impact animal physiology, possibly leading to heat stress.Heat stress occurs when excess body heat cannot be dissipated adequately to maintain thermal equilibrium, usually due to a combination of external factors such as elevated ambient temperature, high humidity, and lack of ventilation (Wang et al., 2020).Indoor-housed livestock are at high risk for heat stress given they are typically housed in confined areas.Under heat stress, cows produce less milk (Bohmanova et al., 2007;André et al., 2011;Bernabucci et al., 2014), and have impaired fertility (Ray et al., 1992;Rensis and Scaramuzzi, 2003;Sammad et al., 2020), leading to both economic and health and welfare concerns.
Physiological changes linked to heat stress are often accompanied by behavioral responses including a reduction in overall activity, to conserve energy and minimize heat production (Cook et al., 2007).In contrast, animals may modify their behavior in cold conditions by increasing activity to generate body heat (Rabaiotti and Woodroffe, 2019;Herbut et al., 2021).Other behavioral responses to heat stress include decreased lying time, particularly above 19°C to 21°C, and increased standing time which serves to increase body surface area exposed for heat loss (Cook et al., 2007;King et al., 2016;Tresoldi et al., 2019).Additionally, cattle reduce their feeding behavior at high temperatures, thought to serve to reduce metabolic heat production (Bouraoui et al., 2002), which is likely to lower nutrient intake and may contribute to reduced conception rates under high temperatures (Marai et al., 2002;Morton et al., 2007;Sammad et al., 2020).This is accompanied by an increase in drinking, to cool individuals down, and an increase in competitive events at resources (Pereyra et al., 2010;McDonald et al., 2020;Tsai et al., 2020).Interestingly, subordinate cows appear to avoid drinkers during the hottest and most competitive hours of the day (McDonald et al., 2020).It is therefore important to consider that behavioral responses to increased heat load may be dependent on the individual animal, housing and the time of day.
Bunching behavior (clustering, spatial aggregation) refers to the close grouping of individual animals within a herd.There is evidence that livestock may bunch when their environment is warmer than typical ambient temperatures including indoor (Mader et al., 2002;Erbez et al., 2012;Javorová et al., 2014) and outdoor-housed cattle (Lefcourt and Schmidtmann, 1989).The reason why cattle may bunch in warm conditions is not well understood, but one possibility is that bunching may occur around shaded areas as a strategy to reduce heat load; tree shading can reduce the body temperature of cattle at pasture (Kendall et al., 2006;Giro et al., 2019).This theory fits less well for indoor housed cattle where they are predominantly blocked from direct sunlight, however shafts of sunlight may still project through panels and gaps in the roofing or walls.It is also possible that cattle group around water at high temperatures since drinking is increased to decrease body temperature.Studies in temperate regions have shown evidence for bunching at high temperatures above a particular threshold (Lefcourt and Schmidtmann, 1989;Erbez et al., 2012;Javorová et al., 2014).In Czech dairy barns, cows were found to bunch more with increasing ambient air temperature above a threshold of 19-20°C (Erbez et al., 2012;Javorová et al., 2014).However, a limitation of these earlier studies is that spatial data collection was not obtained using an automated location system; cow locations were estimated using either video recordings or photographs at 15-min intervals (Erbez et al., 2012;Javorová et al., 2014).
Bunching may also occur at low ambient temperatures, a known response throughout the animal kingdom to help animals thermoregulate (Huynh et al., 2005;Gilbert et al., 2010;Ozella et al., 2020).Similarly, animals may bunch in certain locations due to reasons other than temperature.For example, cows may bunch around areas with good air flow to avoid insects which gather around moist organic material, particularly during warm conditions.Moreover, bunching may help distribute fly burdens; individuals position their tails on the outside of a group for protection through leg kicking and tail swishing (Schmidtmann and Berkebile, 1985;Ashmawy et al., 2019).Likewise, defensive aggregation is a known response of animals to an external threat (Foster and Treherne, 1981;Mooring and Hart, 1992;Bowen et al., 2013).In cattle, this instinctive response may activate in the presence of heightened fly burdens, or during heat stress when flight responses are impeded, leading to increased vigilance.However, if stocking densities within a barn are high, or if significant bunching occurs, body heat within the herd may contribute to an overall increase in the localized temperature, further increasing the internal body temperature of individuals in a potentially negative feedback cycle (Schütz et al., 2010;Polsky and Keyserlingk, 2017).Temperature has indeed been found to increase in dense areas such as milking parlors and holdings areas (Papez and Kic, 2015;Celozzi et al., 2020).If bunching is determined to be a maladaptive response of cattle to increasing temperature, this could have significant implications for the health, welfare and profitability of herds (André et al., 2011;Galán et al., 2018;Ekine-Dzivenu et al., 2020).This would be particularly concerning given the growing demand for increased yields and herd sizes across the dairy industry, alongside the projected increase in heat stress caused by climate change (Ji et al., 2020).
Direct behavioral observations by human observers are time-consuming and may be subjective.However, automated spatial positioning sensor systems used as part of Precision Livestock Farming (PLF) approaches to provide continuous tracking of animals and their environment, offer a promising solution to help better measure and understand behavior (Gygax et al., 2007;Tullo et al., 2016).Automated sensors have the potential to collect high-resolution spatial data which can be analyzed through time series data (Cagnacci et al., 2010) to determine individual time budgets, spaceuse patterns, social interactions, and detect behavioral patterns and changes linked to health issues (Barker et al., 2018;Vázquez Diosdado et al., 2018;Chopra et al., 2020).In this study, we use a local positioning sensor system (LPS) to continuously monitor a housed dairy herd over 4 mo to explore how bunching behavior, defined using 4 different metrics, may be associated with higher ambient temperatures.

Animals and housing
A Holstein-Friesian high-yielding (305-d milk yield of 10,909 L) indoor-housed dairy cow herd was monitored continuously on a commercial farm in southeast England, from 1st August to 30th November 2014.The total number of different individual animals within the herd over this period was n = 127, and a minimum and maximum number of unique cows present per day during the study period was n = 86 and n = 111 respectively.The herd was continually housed in a closed freestall barn with open gable ends containing 98 useable cubicles, at a high stocking density (feed space = 0.53m to 0.68m per cow, lying space = 0.88 to 1.14 cubicles per cow) with free access to the feeding passage and to water troughs.A layout of the barn is shown in Figure 1.The cows were milked thrice daily (at approximately 5am, 1pm and 9pm), with milk yield recorded monthly, and cows were fed a total mixed ration during the first milking (which was subsequently pushed up several times later in the day).Analyses of spatial tracking data for the same herd, in different contexts and over different time periods, were previously reported in (Barker et al., 2018;Vázquez Diosdado et al., 2018;Chopra et al., 2020).

Local positioning system and data pre-processing
Each cow was equipped with a combined real -time local positioning sensor, temperature sensor and accelerometer (see Appendix 1), Omnisense 500, mounted on a weighted neck collar to ensure stability.The Omnisense positioning sensors form a localized radio network which triangulates signals (at 0.1 Hz) between all (fixed and cow-equipped) sensors to determine the spatial position of each individual animal.Fixed sensors were strategically placed around the barn to improve the triangulation network and to fix the absolute spatial position of each sensor.In a separate study in the same barn environment, this sensor system showed a 50% circular error of probability (CEP) measurement of 1.07m (i.e., 50% of all recorded locations lay within a 1.07m radius of the mean location of static sensors), which compares favorably to the commercially advertised CEP specification of 1m, and mean distance error of 2.66m for a static sensor, and 1.90m and 2.90m respectively for an animal-mounted sensor (Barker et al., 2018).
All data processing and analysis was conducted in R for Windows with RStudio (RStudio Team, 2020; R Core Team, 2022).Extended interruption occurred on 15 of the study days due to the system malfunc-tioning and resetting part-way through the day; these days were excluded from the subsequent analysis.The number of cows in the barn on a given day across the study period (not including removed days) varied from n = 86 to 111, and a total of 88,576,716 location data points were collected from these cows.
For the analysis, we excluded hours when most cows were in the milking parlor or collecting yard (05:00-07:59, 12:00-14:59, 20:00-22:59), since their behavior was constrained by farm staff during these periods.The sensor system reset at midnight each day for data upload and hence times between 23:00 and 00:59 were also excluded.A total of 51,630,512 data points remained after removing milking and system reset periods (below, this is referred to as 'original data').
Subsequent data pre-processing consisted of 1) removing location data outside of a 3m buffer of the functional zones due to minor positional inaccuracies (15.88% of original data; 8,201,449 data points removed), 2) removing nonsensical positional data (e.g., sensors "stuck" in exactly the same or a similar, point location for multiple consecutive time points, often shortly after the system reset at the end of each day; 1.89% of original data; 975,291 data points removed), 3) smoothing using a simple moving average (SMA) with a 2-sided window size of 15 data points (0.32% of original data removed as 7 data points were removed from the start and end of the location time series for each sensor per day; 166,555 data points), and finally, 4) plots of trajectories were visualized manually and any further nonsensical data was removed (e.g., where the trajectory of a cow was located within an unrealistically small area for an extended period of time); these instances were deemed as biologically implausible given these cows were not subject to any farm management actions during these times (2.17% of original data; 1,119,142 data points removed).In total, 79.74% of the original data remained after these stages, leaving a total of 41,167,975 data points for the main analysis.Further details of all 4 data pre-processing steps are outlined in Supplementary Material 1 of Chopra et al. (2020) where a similar approach was used.

Environmental data
Wall-mounted temperature sensors automatically recorded the ambient barn temperature (BT) continuously throughout the study period; a mean hourly measurement was calculated across all sensors (n = 26 sensors, where n = 22 recorded every 8 s and n = 4 recorded every hour) (refer to Appendix 2 for further details on temperature sensors and data collected).The sensors were placed at approximately equal distances down the length (x-axis) of the barn at a height of approximately 2.5-3m (Figure 1).
Further temperature data were also obtained from the cow-equipped positional sensors (ST) and from a local weather station (LT), Harold Hill, located approximately 11km from the study site.The ST sensors were housed in an insulated box and were located close to the body of each animal, leading to a buffering effect when comparing the ST ranges to the BT range (Figure A1 B-C in Appendix 2).Meanwhile, LT was used only as a comparator baseline to check temperature trends; LT was weakly correlated to both BT and ST (Figure A1 in Appendix 2), which is unsurprising given the short distance from the study site.

Bunching metrics
There is no formal definition of 'bunching' behavior in the context of PLF, although a range of different metrics are available to measure spatial aggregation, spatial clustering, or social proximity from positional data.Previous studies have used direct observations or basic nearest neighbor aggregation indices to measure bunching, the latter of which would provide a distorted score if cows tend to form dyad pairs (Lefcourt and Schmidtmann, 1989;Mader et al., 2002;Erbez et al., 2012;Ashmawy et al., 2019).Hence, we consider 4 complementary 'bunching metrics', which directly measure the level of spatial separation (i.e., the inverse of bunching), each measured on an hourly basis across the full herd: range size (core and full), mean inter-cow distance (ICD), and mean nearest-neighbor distance (NND).

Range size
Range size was calculated by overlaying a virtual grid (1.5m x 1.5m = 2.25m 2 cells) over the barn map (Figure 1A).At each timestep, location data were used to assign each individual cow within the herd to a given cell (or rejected if it lay outside the main barn region, see Section 2.2).The number of individual data points assigned to each cell across the full herd and over a full hour (360-time steps at 0.1 Hz) were counted giving a final hourly total for each virtual cell within the barn.Following the standard methodology outlined in Vázquez Diosdado et al. ( 2018), the highest density cells cumulatively adding to 50% or 95% were then used to create an hourly utility distribution (space-use intensity map, Figure 1B-C) corresponding to the 'core range' and 'full range' respectively.The number of unique virtual cells included in each hourly core range and full range distribution are then defined as the range size (CR and FR respectively, measured in units of 2.25m 2 virtual cells).

Inter-cow distance
At a given time-step, the herd mean inter-cow distance is the mean distance (in meters) between each of the (N 2 -N)/2 possible dyad pairs across a herd with N individuals.An hourly value (which we define as ICD, measured in meters) is then calculated as the arithmetic mean of all (360) values recorded over that hour (assuming 0.1 Hz sampling frequency): where (i, j) is a given dyad pair and i ≠ j.

Nearest-neighbor distance
At a given time-step, the nearest neighbor distance for a given individual is simply the distance (in meters) to its closest neighbor in the herd (i.e., the smallest inter-cow distance when considering all dyad pairs involving that individual).A mean value can then be determined across the full herd on an hourly basis (which we define as NND, measured in meters): where j indexes all possible dyad pairs for i across the N individuals in the herd.
Although the 4 bunching metrics defined above measure different aspects of space-use and social proximity, they will all decrease as spatial separation decreases and spatial clustering ('bunching') increases.Range size is independent (in a statistical sense) to ICD and NND, but as it is aggregated at an hourly level from the outset, it will not capture short-time variations in proximity between individuals.Nevertheless, a smaller (core or full) range size indicates that the herd as a whole are using less of the barn and are hence sharing a smaller space (Figure 1B, C).NND and ICD are loosely dependent but capture different aspects of herd bunching.ICD measures distances across the full herd and hence can be skewed by individuals that are far apart or in separate subgroups (e.g., a herd split into 2 dis-tinct and tightly bunched subgroups could still have a high ICD value, even though the separate groups would be considered to be bunching).Meanwhile, NND only accounts for the distance between animals nearest each other and hence may not capture bunching between multiple individuals (e.g., a herd split into spatially separated dyad pairs could have a very low NND, but would have high ICD, and wouldn't usually be considered as herd-level bunching).By measuring the above 4 metrics simultaneously we aim to accurately capture a wide range of herd dynamics that may be indicative of bunching behavior.

Segmentation and breakpoints
Motivated by findings from previous studies where behavior was observed to change above a certain threshold (Erbez et al., 2012;Javorová et al., 2014), we use segmentation to explore whether bunching behavior is affected by ambient temperature.A segmented relationship describes changes in a trend at specific values, with a point of change referred to as a breakpoint (BP) (Muggeo, 2003); we consider 2 complementary approaches: Hypothesis-driven. motivated by, and consistent with, previous studies on dairy cow behavior (Chen et al., 2016;Tresoldi et al., 2019), we hypothesize that changes in herd behavior may occur above 20°C.To test this hypothesis, we define a breakpoint (BP) at 20°C.Linear regression models were then created for data <20°C and ≥20°C, for each bunching metric, using the 'lm' function in base R. Chow tests were performed to test the significance of the hypothesized BP of 20°C, for all bunching metrics.For this, the function 'sctest' from the R package 'strucchange' was used (Muggeo, 2003(Muggeo, , 2008)).
Data-driven. the R package 'segmented' was used to identify the locations of forced BPs (Muggeo, 2008).For each bunching metric, the original linear model was inputted into the 'segmented' function, specifying one or 2 BP(s) in BT.Linear regression models were then created using the data-driven BP (i.e., 2 models where one BPs identified: BT < BP and BT ≥ BP; or 3 models where 2 BPs identified: BT < BP1, BP1 ≤ BT < BP2 and BT ≥ BP2) and the results analyzed.Chow tests were performed to test the significance of these BPs.
For both approaches, default R settings were used in all linear models, which exclude observations containing NA values and uses QR factorization, an extensively used method for solving least squares problems.Model assumptions were considered before drawing conclusions: heteroskedasticity through residual vs fitted plots and NVC tests, and residual normality through QQ-plots and the Shapiro-Wilks test (see Appendix 3).
The impact of time of day on the relationship between bunching behavior and temperature was explored to consider the presence of diurnal patterns, differences between the day (morning and/or afternoon) and night, and to consider that milking was a deterministic event at certain points of the day.To do this, the hourly mean data were divided into 3 intervals between milking times: night (01:00 to 04:00; no values where BT ≥20°C across the study period), morning (08:00 to 11:00; n = 33 h where BT ≥20°C) and afternoon/evening (15:00 to 19:00; n = 142 h where BT ≥20°C).Linear regression models were then created for data where BT ≥20°C, for each bunching metric and each time interval, using the lm() function in base R.

Bunching in relation to barn location and localized temperature
Bunching and barn location.Space-use intensity distributions were determined (as described in Section 2.4.1) for those hours corresponding to the lowest 10% of values (n = 127 h) for each individual bunching metric (CR, FR, ICD and NND; i.e., corresponding to high levels of bunching for each metric separately).A combined space-use distribution was also determined by including the union of all hours identified in the previous step across all bunching metrics (n = 321 h).These aggregated space-use distributions corresponding to periods of high bunching were then compared with the aggregated space-use distribution determined across the entire study period to determine whether animals were distributed differently across the barn during bunching events compared with non-bunching periods.
Space-use intensity and temperature.The animal-mounted sensors were used to determine a mean temperature across the entire barn, which we define as ST (see Appendix 2).By matching the spatial location and sensor-measured temperature at each time point for all animals across the herd (i.e., by adapting the method described in Section 2.4.1), it is also possible to determine a spatially-resolved mean hourly temperature for each 2.25m 2 [1.5m x 1.5m] virtual cell of the barn, which we define as ST XY .Hence, we can create a spatial (heat)map of the variations in temperature across the barn over different time periods (e.g., hourly, monthly).It should be stressed that both ST (the mean sensor measured temperature across the whole barn; i.e., the mean of all ST XY values for a given hour) and ST XY are not true measures of either the ambient or localized temperature(s) within the barn; as the sensor is mounted directly on a neck collar (albeit in an insu-Chopra et al.: Dairy cow bunching at higher temperatures lated housing) it will be affected by the heat generated by the animal itself.Nevertheless, relative differences in ST XY across the barn offer a useful measure of how the temperature experienced by each animal, and by the herd as a whole may vary across the barn.
Comparison of spatial distributions: the Bhattacharya coefficient.To compare space-use distributions over different periods or to compare variation in space-use against temperature across the barn, we use the Bhattacharya coefficient (BC) which provides a quantitative measure to assess the level of similarity between 2-dimensional matrices: where P x and Q x are the (normalized) distributions being compared.The resulting BC ranges from 0 to 1, where 1 indicates complete overlap between the 2 distributions, while 0 indicates complete dissimilarity.
The BC only gives a single measure of similarity between 2 distributions; to test if this similarity is significant, it is necessary to undertake a randomization test.Our procedure involved calculating the mean BC for each hour of our study (comparing temperature and space use distributions).The same number of BC values were generated by randomly shuffling the order of hourly space use and temperature distribution matrices, and the mean was determined for these randomized BC values.This randomization process was repeated (n = 10,000 iterations), and the number of times, k, the randomized mean BC was less than the original mean BC from the (non-randomized) data was used to directly determine a p-value (where P = k / 10000).
To compare the spatial distribution for (union of all) hours of highest bunching (see Section 2.6.1) with the spatial distribution over the entire study period, we employed a similar randomization approach.First, we computed the BC for each of the highest bunching hours (n = 321 h) when compared with the entire study period.The mean BC was then calculated.Subsequently, we randomly selected the same number of hours (n = 321) across the study duration and computed the BC for each of these hours in relation to the space-use distribution during the entire study period.This process was repeated 10,000 times, recording the instances where the mean BC of the randomized data exceeded the original mean BC.Finally, the count of such occurrences was divided by the total number of iterations to derive the p-value.A similar randomization procedure was also undertaken for each bunching metric separately (with n = 127 h in each case; see Appendix 4).

Trends in temperature and bunching
Hourly ambient barn temperature (BT) peaked during the afternoon, and dropped during the night, as expected (Figure 2A).The coolest median hourly BT was recorded in November (11.77°C; n = 260 h) and the hottest in August (18.330°C;n = 374 h), with mean for September and October as 18.11°C (n = 276 h) and 15.11°C (n = 359 h) respectively.There is a significant difference in mean hourly temperature between months (Friedman chi-squared = 229.48,P < 2.2e-16), with August and September significantly hotter than October (Wilcoxon rank sum (W) = 30552 and 23677, P < 2.2e-16, with Bonferroni Correction) and November as significantly cooler than all other months (W = 4927, 3696, 77063for August, September and October respectively, P < 2.2e-16, with Bonferroni Correction) (Figure 2A).It should be noted that the mean hourly barn temperature across days reaches above 20°C for some hours during the afternoon in August and September, but not during October and November (Figure 2A).More specifically, Figure 2F shows mean hourly barn temperature over time for the entire study period; most of the values which lie above 20°C occur in August and September (93%, n = 320 of 344 h).
Visualization reveals few clear patterns in the hourly bunching metrics over time, other than a clear decrease in bunching during milking periods; this is expected since cows are located within the restricted space of the milking parlor during these periods (and hence we remove these times from the subsequent segmentation analysis) (Figure 2B-E; Figure 2G-J).

Bunching behavior above and below 20°C.
Above or equal to 20°C, a clear negative trend is shown between all 4 bunching metrics and barn temperature.Although the breakpoint is non-significant for 3 of the 4 bunching metrics, these trends are significant for core range (CR), full range (FR), and inter-cow distance (ICD) (Table 1; Figure 3E-G), but non-significant for NND (Table 1; Figure 3H).Below 20°C, the relationships between barn temperature and FR, ICD and NND are also negative and significant, but are non-significant for CR, e.g., for CR: < 20°C, correlation coefficient (e) = −0.15 and for ≥20°C e = −2.08,and respectively (Table 1; Figure 3A-C).The direction of the relationship between NND and barn temperature changes from positive to negative either side of the 20°C breakpoint (Table 2; Figure 3D).Qualitatively similar results were obtained when considering relative humidity from a lo-Chopra et al.: Dairy cow bunching at higher temperatures cal weather station, through the Temperature Humidity Index (THI) (Appendix 5).
For August, at a given barn temperature, the range size values are generally lower in comparison to September (Figure 3I-L).For example, at 20.1°C (rounded to nearest one decimal point), CR = 60 to 65 (n = 2 h) whereas CR = 65 to 78 for September (n = 4 h).The pattern is the opposite for ICD and NND (Figure 3M-P).This can be explained by the number of cows present in the barn, consistently lower in August (n = 85 to 86 cows) compared with in September (n = 92 to 96 cows).Moreover, there are negative trends between barn temperature and the bunching metrics ≥20°C when considering August and September separately (Table 1; Figure 3I-P).
Behavior above and below data-driven breakpoints.Forcing a data-driven method to find the location of barn temperature breakpoints gives values close to 20°C for all the bunching metrics, except for FR (Table 3; Figure 4).The (significant) negative trends using all 4 bunching metrics also hold above these forced breakpoints (Table 2; Figure 4), giving qualitatively similar results to using a hypothesized BP of 20°C (Table 1, Figure 3).Qualitatively similar results were found using forced breakpoints for THI (Table A3 in Appendix 5).
Effect of time of day for breakpoint of 20°C.Barn temperature did not reach ≥20°C during the night, and most data points ≥20°C were recorded during the afternoon/evening (n = 142 and 33 data points, during the afternoon/evening and the morning respectively).Above (or equal to) 20°C, the negative relationship between barn temperature and the range sizes and ICD, change from non-significant during the morning to significant during the afternoon and evening (Figure 5A-C The relationship between barn temperature and NND remains non-significant and negative across the 2 time periods (Figure 5D; morning: e = −0.05,SE = 0.03, t-value = −1.53,P = 0.14; afternoon/evening: NND: e = −0.01,SE = 0.01, t-value = −0.66,P = 0.51).The highest instances of bunching are mostly during September and August (n = 8 of 9 h; Table 3), which is consistent with results from the segmentation analysis as these are also the hottest months of the study period (refer to Figure 2).The median location of the herd is x = 31.77m and y = 4.52 m across the entire study period (white cross in Figure 6A) and x = 32.01 m and y = 4.63 m during instances of high bunching (white cross in Figure 6B), indicating very little difference in median location between bunching and non-bunching periods.Both coordinates are close to the milking parlor (29.40 m ≤ x ≥42.25 m and (−8.42 m ≤ y ≥ −1.62 m) between the high-density cells in the FZ and the right-hand cubicles (refer to Figure 1).The case is similar when considering the highest instances of bunching for each bunching metric separately (see Appendix 4).

Bunching and barn location
A mean BC = 0.903 was calculated when comparing space-use distribution during the union of specific hours of high bunching (n = 321) to aggregated spaceuse over the entire study period (Figure 6A-B).The randomization procedure (comparing space use across n = 321 randomly selected hours to the entire study period) shows that, although it is still relatively high, the BC during bunching periods is significantly lower than expected by chance (range of mean BC = 0.914 to 0.921 over 10,000 iterations, P < 0.0001), indicating more variation in the spatial distribution during bunching periods when compared with the rest of the study period.Similar results are obtained when considering each bunching metric separately (Appendix 3).
Bunching and temperature variation across the barn On a monthly average basis, ST XY can be seen to vary across the barn with the central areas closest to the feeding passage and exit from the milking Chopra et al.: Dairy cow bunching at higher temperatures Table 1.Outputs from linear regression assessing the relationship between barn temperature and bunching, below and above the breakpoint of 20°C (n = 170 data points), for the full data set ('full') and for August and September alone (O, A and S, respectively), excluding outliers (n = 5 data points).Chow test outputs testing the significance of the breakpoint are also shown.Each test is conducted for each bunching metric (CR = core range, FR = full range, ICD = inter-cow distance, NND = nearest neighbor distance).Periods during which cows were in the milking parlor or collecting yard (05:00-07:59, 12:00-14:59, 20:00-22:59), and periods during which the sensor reset ( 23 parlor consistently hotter across all months than the gable ends (left and right sides) of the barn (Figure 6C-F).Visual comparisons of Figure 6A and Figure 6C-F suggest that, on average and across the study period, those areas of the barn with highest space-use intensity are similar to those areas where ST XY is higher.
This finding is further highlighted when comparing spatial variation in ST XY with space-use intensity across the barn across the entire study, the hourly Bhattacharya coefficient (BC) varied from 0.775 to 0.927, with medians of 0.877 for August,0.900for September, 0.898 for October, and 0.897 for November (all values exclude milking periods).Illustrative examples of the spatial distribution highlight that the coefficient of variation (CV) is very low for temperature (ST XY ), with a median hourly CV of 3.65% (range = 1.41% to 9.41%), indicating low variation across the barn (Figure 6G-H).In contrast, space-use intensity is highly variable across the barn with a median hourly CV of 98.88% (range = 78.13% to 218.38%); this suggests that it is variation in space-use intensity across the barn that is driving differences in the BC rather than temperature variation.
Over the full study period, the mean BC between the hourly space-use intensity and temperature distributions across the barn was 0.888; randomization suggests this is highly significant (range of mean BC = 0.857 to 0.863 for randomly paired space-use and temperature distributions; P < 0.0001), suggesting that there is significant high correlation between space-use intensity and the localized temperature within the barn as measured by animal-mounted sensors.
There is a significant negative relationship between the hourly (temperature and space-use distribution) BC and hourly BT, and this negative relationship is stronger ≥20C compared with <20C (Figure 7E; overall: e = −0.001,SE = 0.0001, t-value = −7.20,P = 1.02e-12where BT <20: e = −0.0004,SE = 0.0002, tvalue = −2.21,P = 0.03; where BT ≥20C: e = −0.007,SE = 0.001, t-value = −7.21,P = 1.23e-11).The fact that BC decreases at higher temperatures is most likely due to bunching (which we previously observed is more likely at higher temperatures); bunching will cause higher variation in space-use intensity across the barn, which in turn will lead to a lower BC since temperature is consistently more uniform across the barn.

DISCUSSION
This study is the first to use automated sensors, continuously over 4 mo as part of a PLF approach, to highlight associations between bunching behavior and higher ambient barn temperatures in an indoor dairy herd in a temperate region (UK).Four different measures of bunching behavior were considered, accounting for various aspects of space-use and social proximity: core and full range size (CR and FR respectively), inter-cow distance (ICD) and nearest neighbor distance (NND).In addition, our novel methodology, which uses combined animal-mounted sensors that measure both location and localized temperature at high spatiotemporal resolution, enables us for the first time to map and compare space-use and localized temperature variations across the barn over time.Our results highlight the importance of understanding animal space-use behavior within the context of farm management under the increased threat of rising temperatures due to climate change.

Bunching behavior increases at higher ambient barn temperatures
Herd bunching behavior increased with increasing barn temperature, and this relationship is significant above or equal to 20°C (Table 1; Figure 3) for CR, FR, and ICD (noting that a decrease in these metrics corresponds to increased bunching).For NND, a negative relationship was similarly found above or equal to 20°C, although this is non-significant, and this relationship did not hold below the 20°C threshold (Table 1; Figure 3).The increase in bunching with temperature ≥20°C Chopra et al.: Dairy cow bunching at higher temperatures Table 2. Data-driven barn temperature breakpoints (BP) and corresponding results assessing the linear relationship between barn temperature and each bunching metric: CR = core range, FR = full range, ICD = inter-cow distance, NND = nearest neighbor distance, on the commercial farm in Essex.Significant p-values are in bold.Outliers have been excluded, based on linear model assumption testing (n = 5).Periods during which cows were in the milking parlor or collecting yard (05:00-07:59, 12:00-14:59, 20:00-22:59) or times when the sensor reset marked with dashed vertical lines), for core range (CR; 50%; measured in virtual cells of 2.25m 2 ), full range (FR; 95%; measured in virtual cells of 2.25m 2 ), inter-cow distance (ICD; m) and nearest neighbor distance (NND; m), respectively.Outliers (n = 5 data points) are excluded, based on model assumption testing (refer to Table 2).A single point represents an average per hour and points are colored by month: August = purple, September = green, October = gray and November = black.Periods during which cows were in the milking parlor or collecting yard (05:00-07:59, 12:00-14:59, 20:00-22:59) or periods when the sensor reset (23:00 to 00:59) were removed.
held across the morning to the evening, but this trend was stronger during the afternoon/evening (15:00 to 19:00, n = 33) compared with the morning (08:00 to 11:00, n = 142) for the CR, FR, and ICD, potentially partly attributed to the smaller sample size available for the morning analysis (Figure 5).
The main finding that cows significantly increase their bunching behavior above warmer than average ambient temperatures (20°C) (Table 1, Figure 3), supports previous observational studies, where cows were found to crowd above but not below an ambient temperature of 19°C to 20°C (Erbez et al., 2012;Javorová et al., 2014).In these earlier studies, bunching was measured by counting the number of individuals in barn segments every 15 min (Erbez et al., 2012;Javorová et al., 2014).Our findings are consistent with this earlier work but provide more detailed understanding of bunching behavior in indoor dairy herds; our data were derived from higher-frequency recordings and by using 4 different metrics of bunching (FR, CR, ICD, NND) we can con-sider various aspects of space-use and social proximity.Specifically, our results indicate that the herd used less space as a whole at high barn temperatures (FR, CR), and not only may dyads (NND) position closer together at higher temperatures, but multiple individuals may, too (ICD) (Table 1, Figure 3).Conversely, bunching has been found to decrease with ambient temperature ≤30°C in temperate and arid climates (Ashmawy et al., 2019), but this earlier study only categorized bunching using a binary method (bunching or not bunching) weekly, twice per day.Another previous study used a more detailed method to measure bunching in a dairy herd in a temperate climate, observing and categorizing individuals into one of 4 nearest neighbor groups, with recordings taken every 20 min for 10 d, then one day per week for 6 weeks (Lefcourt and Schmidtmann, 1989).Although bunching appeared to increase with THI (using an ambient temperature measure), a specific temperature threshold where behavior qualitatively changed was not reported, which was potentially at- tributed to the study's consistently high daily ambient temperatures recorded i.e., weekly least mean squares of daily ambient temperature between 21°C and 28°C (Lefcourt and Schmidtmann, 1989).
The increase in bunching ≥20°C, in terms of 3 of our 4 metrics (FR, CR, ICD), is found to be clearer during the afternoon/evening (15:00 to 19:00; maximum = 27.77°C and median = 18.25°C) than during the morning (08:00 to 11:00; maximum = 25.25°C and median = 16.42°C)(Figure 5).This is consistent with previous studies whereby cows were found to crowd during the afternoon, where temperatures were consistently recorded above 22°C, but not during the night/ early morning, where temperatures did not exceed 11°C (Erbez et al., 2012;Javorová et al., 2014).Alterations to the spatial positioning of the herd during specific Chopra et al.: Dairy cow bunching at higher temperatures periods of the day are likely to disrupt normal behaviors, such as feeding if stocking rates at the feed face are increased (Hill et al., 2009).Such changes have the potential to disproportionately affect cows with health issues; for example, lame cows have been reported to moderate their feeding times to avoid competition at the feed face (Blackie et al., 2011;Barker et al., 2018).The lack of bunching we found overnight (01:00 to 04:00), as temperatures remained below 20°C, suggests these periods are important for cows to undertake activities which may have been avoided during the day (Figure 5).The rise of tropical nights in the UK (Hanlon et al., 2021) may further diminish the ability of cows to adapt to changing climatic conditions and to catch up on important behaviors, such as lying (Schütz et al., 2008).
Our hypothesized temperature breakpoint of 20°C was motivated by previous studies on dairy cows that reported physiological or behavioral changes at or around this threshold (Berman et al., 1985;Nonaka et al., 2008;Dikmen and Hansen, 2009;Hut et al., 2022).Results from using an alternative approach with purely data-driven forced breakpoints support our main conclusions (Table 1, Figure 3).These forced breakpoints were found to be close to 20°C for 3 of the 4 bunching metrics (CR, ICD, NND; Table 2, Figure 4), and the negative trend above these data-driven thresholds is consistent with results from the 20°C breakpoint (Tables 1-2, Figure 3-4).

Other factors that may affect bunching behavior
Many studies have considered the Temperature-Humidity Index (THI), which considers both the effects of relative humidity and temperature, in the context of behavioral responses of cattle to heat stress e.g., (Tresoldi et al., 2019;Tsai et al., 2020).Relative humidity was not recorded within the barn so data were obtained from a local weather station, and results show that the relationships between the bunching metrics and THI were largely qualitatively similar to those for barn temperature alone (Appendix 5).Local or regional measures of relative humidity from local weather stations are unlikely to reflect the fine scale microclimate changes experienced by cows in different areas of the barn at different times.Hence, further investigation should consider finer-scale humidity readings recorded within the barn, as well as other environmental factors (e.g., airflow, methane, or carbon dioxide levels, etc.) to better understand other possible factors contributing to bunching behavior.
Cows may alter their spatial positioning during the day due to a preference for specific areas of the barn at high ambient temperatures.Possible factors which might have resulted in the presence of microclimates in the study barn (i.e., temperature, humidity and air flow) were 1) mechanical ventilation provided by fans in the middle and right of the barn 2) water troughs located in the center cross passage and 3) sunlight radiating through the wooden cladding on the left gable end side (west-southwest) of the shed (due to the position of the parlor, no sunlight would come through the right side of the shed).Cows may have detected microclimatic variations in the barn due to these factors, as found in previous studies i.e., at high temperatures, cattle are known to seek shade and avoid sunlight (Bennett et al., 1985;Kendall et al., 2006;Schütz et al., 2010), but may be drawn to areas with water to cool down and good air flow to avoid flies.These factors may also contribute to increased activity at higher temperatures (see Appendix 1); for instance, cows may be more active as a group at higher ambient temperatures to seek preferable microclimates or to deter fly burdens through behavior such as tail swishing (Ashmawy et al., 2019).However, activity and ambient temperature are not entirely independent; an increase in cow activity may increase the ambient temperature in the barn, given activity generates body heat, although the extent to which is unclear.It is not clear to what extent each of these microclimatic factors affect their bunching behavior, but it can be hypothesized that cows are either drawn to preferable microclimates or driven away from adverse microclimates when temperatures increase above 20°C.
We found no evidence that the average location and spatial distribution of cows was different during bunching periods when compared with the aggregated spatial distribution across the entire study (Figure 6A-B).However, there was more variation within the observed spatial distribution when cows were bunching (median hourly CV of 98.88%; range = 78.13% to 218.38%; Figure 7C-D).In general, cows appeared to prefer positioning in the right-hand cubicles over the left and lower cubicles in the non-feeding zone, across the entire study duration and during hours of high bunching (Figure 6A-B).The right-hand cubicles were closer to the collecting yard/milking parlor which would have allowed cows to relieve their high milk yield sooner compared with if they were positioned elsewhere in the barn (the study group consisted of high-yielding cows).It has previously been observed that higher parity cows and cows with more days in milk spent greater time in cubicle areas closer to a collecting yard (Vázquez Diosdado et al., 2018).Similarly, (Churakov et al., 2021) found that higher-parity cows occupied cubicles closer to a milking area.As discussed above, this apparent preference for the right-hand cubicles could also be due to other possible factors such as better ventilation in these areas.(Churakov et al., 2021) also reported higher Chopra et al.: Dairy cow bunching at higher temperatures occupancy near feeding areas, which is consistent with the high-space use in the feeding zone observed in our study, although this is likely driven by the requirement for feed intake rather than any other factors.
As expected, there was a gradual decrease in both ambient temperature and localized temperatures within the barn over the duration of the study period (August to November; Figure 2F, Figure 6C-F).Localized temperature throughout the barn (ST XY , measured by animal mounted sensors) showed relatively low spatial variation (median hourly CV of 3.65%; range = 1.41% to 9.41%; Figure 6C-F), although visual inspection shows that in general the hottest temperatures were recorded in the passageways and the feeding zone, which is perhaps not surprising as this is where the highest density of cows is typically observed (Figure 6A-B) (Papez and Kic, 2015;Celozzi et al., 2020).There was a significant correlation (Bhattacharyya coefficient, BC = 0.888; P < 0.0001 in randomization test) between localized temperature and space-use intensity across the barn over the study duration.However, the observed BC was negatively correlated with ambient barn temperature (Figure 7E), and this was significant ≥20°C.This result can be explained by the fact that bunching leads to higher spatial variation (higher CV) and hence a lower BC (in comparison to non-bunching periods), and as previously discussed, bunching is more likely to be observed at higher ambient temperatures.Nevertheless, further investigation is required to fully understand the detailed spatiotemporal dynamics of the barn environment and to identify the specific factors driving space-use and bunching behavior, to help improve housing designs.

Implications for herd management
Due to the high metabolic demands of dairy cows their bodies radiate significant heat into their local environment to maintain a stable core body temperature (Spiers, 2003).Bunching of dairy cattle, whether due to management practices such as collection for milking or behavioral choices, has the potential to exacerbate hot and humid conditions capable of increasing heat load i.e., given temperatures are known to increase in areas of high stocking densities (Papez and Kic, 2015;Celozzi et al., 2020), which is consistent with our results in Figure 6.This has implications for welfare, as cows under heat stress undergo physiological changes such as an increase in respiration rate and panting (Osei-Amponsah et al., 2020;Yan et al., 2021) and production losses such as reduced milk yield (Gantner et al., 2011;Gorniak et al., 2014;Wildridge et al., 2018).Moving forward, heat stress is becoming an increasing concern as demand in the dairy sector increases (Bayvel et al., 2012;Cembalo et al., 2016;Miglior et al., 2017) alongside heightened impacts of climate change.In regions such as the UK cattle have historically had limited experience of such high temperatures and hence little time to adapt (Fodor et al., 2018;Hanlon et al., 2021).Further study is required to better understand the relationship between bunching and temperature, alongside impacts on production, considering appropriate heat stress thresholds at which to alert farmers.
Precision Livestock Farming approaches provide opportunities to study and understand animal behavior, allowing us to collect and analyze vast amounts of data in real-time, and to identify patterns and associations that would be difficult to detect using traditional methods such as manual observations (Berckmans, 2014).The findings of this study demonstrate the potential of PLF approaches to uncover previously unknown links between environmental factors and animal behavior (e.g., the overlap between space-use intensity and localized temperature, Figure 6), with implications for health, welfare, and productivity.In this study, dairy cows were found to increase their bunching behavior above average ambient temperatures, and possible mitigation strategies could focus on addressing heat loads, such as through increasing access to shade (Mitlöhner et al., 2001), providing sprinklers and improving ventilation (Fournel et al., 2017), and reducing stocking densities (Celozzi et al., 2020).
Our results have highlighted the association between bunching behavior and ambient temperatures in a large housed dairy herd.Using automated sensors and continuous monitoring at high spatiotemporal resolution that provides an unprecedented level of detail, we provide evidence of how bunching behavior of a dairy cow herd increases above 20°C, as measured by a range of different metrics of spatial clustering and social proximity (Table 1).Given that temperatures are predicted to increase under climate change, there is a need for further exploration and understanding of the drivers of bunching behavior to explore potential mitigation strategies for heat stress (Hanlon et al., 2021), particularly if bunching behavior is itself contributing to the formation of localized higher temperatures.Our study underscores the importance of considering behavioral responses in relation to changing environmental conditions, so that we can implement effective strategies to protect animal health, welfare, and productivity within agricultural settings.
Figure 1.(A) Barn layout of the commercial farm in Essex showing the positions of sensors recording barn temperature (total of n = 26; purple triangles).The area shown in gray (upper barn area) was not used by this group of cows during the study period.The barn is divided into the feeding zone and the non-feeding zone where cubicles were located, and P marks the passageway.The upper wall of the barn (i.e., corresponding to y = 30 m in the plot) is orientated facing north-north-west.(B-C) Space-use distributions showing examples of high bunching (B), where core range (CR, 50%, solid gray boundary) = 46 virtual cells of 2.25m 2 , full range (FR, 95%, dashed gray boundary) = 166 virtual cells of 2.25m 2 , and ICD and NND = 12.78m and 2.66m respectively (for 18/09/2014 at 11:00-12:00; barn temperature (BT) = 20.91°C), and low bunching (C), where CR = 88 virtual cells of 2.25m 2 , FR = 268 virtual cells of 2.25m 2 , and ICD and NND = 20.87m and 3.61m respectively (for 12/10/2014 at 18:00-19:00; BT = 13.10°C).(D-E) Examples of 'snapshots' of point location data for all individuals in the herd shown for specific time points within each of these hours (D: 18/09/2014 at 11:01:32 and E: 12/10/2014 at 18:00:02).
Visualization of space-use over the entire study duration suggests that the cows preferred to position in the right-hand cubicles which are nearer to the exit of the milking parlor (26.90 m ≤ x ≥52.85 m and 0.80 m ≤ y ≥10.50 m and 25.10 m ≤ x ≥52.85 m and −1.62 m ≤ y ≥0.80 m) over the left-hand cubicles (2.55 m ≤ x ≥23.10 m and 0.80 m ≤ y ≥10.50 m) in the NFZ (Figure 6A; refer to Figure 1 for barn layout).

Figure 3 .
Figure 3. Linear regression plots showing the relationship between bunching and barn temperature, below and above (or equal to) the breakpoint of 20°C (≥20°C: n = 170 data points).Plots are shown for each bunching metric: core range (CR; 50%; measured in virtual cells of 2.25m 2 ), full range (FR; 95%; measured in virtual cells of 2.25m 2 ), inter-cow distance (ICD; m) and nearest neighbor distance (NND; m), below the breakpoint of 20°C (A-D respectively) and above or equal to the breakpoint (E-H).A single point represents a mean average per hour and points are colored by month: August = red, September = orange, October = pink and November = purple.Insets showing the patterns for August and September alone where ≥20°C are also shown (I-J for CR, K-L for FR, M-N for ICD, and O-P for NND, respectively).Significant relationships are marked with an asterisk (*), and the corresponding p-values are given.Outliers (n = 5 data points) are excluded based on model assumption testing as shown in Appendix 3. Periods during which cows were in the milking parlor or collecting yard (05:00-07:59, 12:00-14:59, 20:00-22:59) or periods when the sensor reset (23:00 to 00:59) were removed.

Figure 4 .
Figure 4. Linear regression plots showing the relationship between bunching and barn temperature, below and above (or equal to) the forced breakpoints.The plots are shown for each bunching metric where the breakpoints are 22°C, 12°C, 20°C and 22°C (nearest decimal place;marked with dashed vertical lines), for core range (CR; 50%; measured in virtual cells of 2.25m 2 ), full range (FR; 95%; measured in virtual cells of 2.25m 2 ), inter-cow distance (ICD; m) and nearest neighbor distance (NND; m), respectively.Outliers (n = 5 data points) are excluded, based on model assumption testing (refer to Table2).A single point represents an average per hour and points are colored by month: August = purple, September = green, October = gray and November = black.Periods during which cows were in the milking parlor or collecting yard (05:00-07:59, 12:00-14:59, 20:00-22:59) or periods when the sensor reset (23:00 to 00:59) were removed.
Figure 6.(A) Aggregated space-use distribution for the full herd over the entire study duration, (B) during the highest hours of bunching defined as the union of the lowest 10% of any bunching metric (n = 321); (C-H) aggregated sensor temperature distribution (per 2.25m 2 virtual cell) for the full herd over the entire study duration for each month: (C) August (mean ST XY = 23.56°C,range ST XY = 21.63°C to 26.77°C), (D) September (mean ST XY = 23.52°C,range ST XY = 20.71°C to 26.42°C), (E) October (mean ST XY = 22.15°C, range ST XY = 19.23°C to 23.53°C) and (F) November 2014 (mean ST XY = 19.75°C,range ST XY = 17.38°C to 20.73°C), and for hours with the (G) lowest barn temperature (20/11/2014; mean hourly barn temperature = 8.29°C; mean ST XY (ST per cell) = 17.06°C, range ST XY = 15.39°C to 19.40°C) and (H) highest barn temperature (07/08/2014; mean hourly barn temperature = 22.58°C, mean ST XY = 26.24°C,range ST XY = 24.03°C to 29.61°C).In (A) and (B), the median position is marked with a white cross.In (C) to (H), sensor temperatures were calculated using recordings from cow-equipped sensors (n = 127 different cows, with a minimum of n = 86 cows and maximum of n = 111 cows present per day across the study period).A color gradient represents the average per 2.25m 2 cell, with lighter shades of blue indicating areas of low density or space-use, whereas darker shades of blue indicating areas of high density in (A) and (B), whereas yellow cells representing low ST XY from 17°C and red cells representing high ST XY of 32°C in (C) to (H).In all subplots, the highest density 2.25m 2 cells cumulatively adding to 50% and 95% are shown within a solid gray boundary and a dashed gray boundary, respectively.Milking times are not included in any of the calculations.The non-feeding zone is where 1.62 m ≤ y ≤10.5 m, −1.6 m ≤ x ≤58.6 m and the feeding zone is where 10.5 m ≤ y ≤17.2 m.

Figure 7 .
Figure 7. (A and C) The lowest Bhattacharyya coefficient (BC) observed between temperature and space-use intensity across the barn (12/08/2014 at 17:00 (BC = 0.76, barn temperature (BT) = 19.52°C,coefficient of variation for A and C = 2.22% and 98.08% respectively, range of ST XY (sensor temperature per 2.25m 2 virtual cell) = 25.50°C to 28.95°C for (A) and range of density = 3.5e −5 to 4.0e −5 for (C).(B and D) The highest BC between temperature and space-use intensity across the barn observed during the study period (08/10/2014 at 18:00 (BC = 0.92; BT = 15.77°C,coefficient of variation for B and D = 3.27% and 88.49%, range of ST XY = 21.54°C to 25.23°C to (B) and range of density = 2.90 e-5 and 1.30e −2 for (D).(E) The relationship between BT and BC over the full study period; and (F) hourly BC across the study period, both divided by month.The highest density 2.25m 2 cells cumulatively adding to 50% and 95% in A-D (core and full range in both C and D) are shown within a solid gray boundary and a dashed gray boundary, respectively.ST XY in A-B was calculated using recordings from cow-equipped sensors (n = 127 different cows, with a minimum of n = 86 cows and maximum of n = 111 cows present per day across the study period).E excludes milking periods and sensor reset whereas F includes these times.
Chopra et al.:Dairy cow bunching at higher temperatures