Ten non-lactating Holstein-Friesian dairy cows (7.1 ± 0.6 years old; 781 ± 16 kg body weight; 7.5 ± 0.2 mean ± SE months pregnant) were selected from the University's commercial dry cow herd. All cows were experienced with strip-grazing pasture but naïve to the VF neckbands and stimuli and had not been used in previous experiments.

Cows were habituated to the VF neckbands for 3 days. On day one, all cows were weighed (Gallagher G02601 TWI weigh scale, Gallagher Group Ltd., Hamilton, New Zealand) and spray marked with numbers on each side of the flank for identification purposes (Tell tail, GEA Fil Ltd., Mount Maunganui, New Zealand). While restrained in the cattle crush each cow was fitted with a neckband. Cows were then moved into home pen 1 (Figure 1) and remained for 3 days to habituate to the VF neckbands, which remained inactive during this time (Lomax et al., 2019; Colusso et al., 2020, 2021). During the habituation period the cows had ad libitum access to water and shade and were fed an amount of oaten hay (Dry matter 83%) to meet their energy requirements for maintenance and pregnancy (ARC, 1980; Moran, 2005), as determined by the following equation (Corbett et al., 1987):

Pre-commercial prototype neckbands (eShepherd™, Agersens Pty Ltd., Melbourne, VIC, Australia) which have been previously described in similar studies (Lomax et al., 2019; Colusso et al., 2021; Langworthy et al., 2021) were used in this experiment. The VF experimental prototype neckbands consist of a strap with a counterweight (total weight approximately 1.4 kg) and an electronic unit (approximately 725 g and 17 cm L × 12 cm W × 13 cm H) positioned on the top of each animal's neck, with two electrodes that contact the skin on the right side of the neck. The VF neckbands use a GPS module within the unit to monitor cow movement and provide real-time data on cow location, heading and speed. A VF boundary separating the IZ from the EZs is specified using GPS coordinates and which is transmitted to the unit using a radio frequency link (base station). As an animal reaches the VF boundary, the unit will emit a distinctive AT within the cow's hearing range (values are commercial in confidence). If the cow remains at that location (stops) or turns away, no EP is emitted. If the animal continues to move through the VF boundary into the EZ, the unit delivers a short, sharp pulse (EP) in the kilovolt range (values are commercial in confidence). This sequence of an AT followed by the EP is repeated if the animal continues through the VF and into the EZ. A grazing function is implemented to account for the natural grazing behavior that can mimic a cow correctly responding to an AT (slow movement forward with pauses). As the cow slowly reaches the VF the function will deliver a maximum of three consecutive AT's, followed by an EP if the cow continues its forward grazing (Langworthy et al., 2021). If a cow was moving above a specified velocity (values are commercial in confidence) a safeguard function would suspend the neckband from delivery any stimuli (Langworthy et al., 2021). The date, time, GPS location, and “stimuli event” which included where the cow was located in relation to the IZ and details of stimuli delivery, were recorded for later download from the unit.

On day four the cows were removed from the home pen and trained to the VF system and stimuli over a 6-day period on a 1.2 ha paddock of irrigated annual ryegrass (Lolium multiforum), as per Lomax et al. (2019) and Colusso et al. (2021). The paddock was then divided into two above maintenance allocations areas that were set at 3,000 kg DM/ha, with allocation 1 (A1) 4,356 m² and allocation 2 (A2) 3,417 m². The remainder of the paddock with 667 m² of grazable pasture was not offered. A water trough was provided that was accessible by all cows across allocations (Figure 2).

On day one of training, VF1 was activated in the morning to provide cows with access to allocation 1 (A1) for a total of 3 days. At the end of day three the VF1 was deactivated and VF2 was activated to provide the cows access to the residual of A1 in addition to the fresh pasture in A2 for three more days. At the end of day six the VF neckbands were deactivated, and the cows were returned to the home pen for 24 days to allow for pasture regrowth. During this period, cows were alternated between home pens 1–3 to avoid soil compaction and fed as per the acclimatization period.

Strip grazing was conducted on the 1.2 ha paddock of annual ryegrass, which was divided into nine areas that consisted of eight pasture allocations and an un-grazed area (Figure 3). The allocations were determined by pasture mass on offer, using the pre-calibrated RPM set to meet the daily maintenance energy requirement (111.3 ± 1.3 MJ DM/cow, equating to 11.1 ± 0.1 kg DM/cow) for the 10 cows. The pre-graze target of each allocation averaged at 2,324 ± 81 kg DM/ha with the post-graze residual target averaging at 1,649 ± 48 kg DM/ha. A second post-graze reading was taken only on the hold-off day 5 (allocation 4) and day 10 (allocation 8) averaging at 1,572 ± 21 kg DM/ha. The pre- and post-readings consisted of sampling across the allocation in a zig-zag transect with 55 reading points for greater estimation of the available pasture.

The VF was activated at the grazing front to contain cows from future fresh pasture allocations, and either a physical fence (physical external fencing of the paddock, Day 1) or a non-electrified polywire tape (days 2–10, not shown in figure) that was moved up with each allocation to prevent cows from accessing previous sections of the paddock. On average each allocation offered an area of 2,403 m² that included an average of 1,223 m² of fresh pasture and 1,180 m² of residual from the previous day. The small residual that was included in each allocation area was to account for any differences between individual cow VF and/or for any GPS drift. One moveable water trough was shifted forward along the fixed fence with each new allocation offered, its position displayed in Figure 3.

All 10 cows were strip grazed as a group over a 10-day period. On every day except for days 5 and 10, the VF was shifted forward in the morning (07:30–08:00 h) to provide cows with a fresh pasture allocation. On days 5 and 10 the VF was not shifted, and the cows were only offered the previous day's residual pasture for 24 h.

Date and time stamped GPS location (HH:MM:SS) and associated stimuli details (AT and EP) were recorded by the neckband. Neckbands were removed after the training period for data download and remained off until 2 days prior to strip grazing when they were placed back on the cows again. The neckbands were then removed at the end of the 10-day period of strip grazing for data collection. The process involved removing each neckband in the cattle crush to access the SD card. A logged text file was downloaded from the SD cards to produce an Excel^® file (version 1908, 2016) and stacked and processed by cows and day as reported in Colusso et al. (2020, 2021). Briefly, the stimuli data was filtered to determine the total number of AT and EP each cow received per day. The grazing function ATs were included in the total sum. There was a total of 66 h of missing data points from a total of 2,400 hours, which reflects only 2.8% of SD data that was not captured.

ArcMap v10.5 (ESRI, Redlands, CA) was used to determine the location of a cow relative to the VF as to ascertain if/when a cow crossed over the VF. In order to achieve this, it was necessary to determine the location of each individual cows VF on each day. This was determined by plotting the spatial coordinates of all ATs in Excel^® for each cow per day and the furthest left and right GPS location (relative to the paddock entry point, Figure 1) were selected to represent the VF location for that cow and this was then repeated to give a total of 20 latitude and longitudes points for each day. All points were then averaged to produce a generalized location (latitude/longitude) of the left and right side of the VF which when joined together represented the VF for the entire herd on that day. In addition, the standard deviation (StdDev) of either side of the VF was also determined. This process was repeated for each of the 10 days of the study. The average and StdDev latitudinal and longitudinal locations were converted to Universal Transverse Mercator (UTM) Zone 56 S using Earth Point ^®. The four variables were then imported into the below formula to determine the difference:

where aVF = averaged virtual fence coordinate (eastings or northings), and sdVF = corresponding standard deviations of coordinates of virtual fence.

The difference (in m) between the average and StdDev VF was then used to create a buffer zone (BZ) to account for variation between individual neckband VF locations. In Google Earth Pro ^® (Version 7.3.3.7786) the VF points were used to draw a polygon of the IZ for each day. A second polygon was then drawn to account for the remaining area of the paddock and thus the EZ. The IZ and EZ polygons were then imported into ArcMap and joined to produce a single polygon of the paddock. The converted GPS data from the neckbands was imported and clipped to the paddock polygon to remove any spurious data points (“cleaned”). Spurious data points were identified as being outside the physical boundary of the 1.2 ha paddock, which the cows were unable to access. On average 3.6% ± 0.63 (ranging from 1.7 to 8.8%) of spurious GPS points were removed in this process.

The cleaned GPS datapoints were then determined to be in either the IZ or EZ polygons and the distance from the VF was calculated. This data was then exported to Excel^® and the meter difference between the average and StdDev VF for each day was applied to the distance column and those data points were re-assigned to the BZ (mean ± SE, 8.4 m ± 1.1). This created three locations a data point could be located within: IZ, EZ, or BZ.

The locational data were processed into bouts, defined as the duration of time a cow spent in a particular zone. These individual time intervals were calculated from the raw time data by the time between zone changes. The bout was designed to allow hourly analysis of cow location, speed, and distance into the EZ. If there were multiple consecutive time stamps found within the same zone, they were collapsed into a single bout. The total daily duration (24 h periods) each animal spent in each of the three zones across the 24 h period was also analyzed. The total duration was calculated by summing up each individual bout by zone within each day for each cow.

Density plots were created in R using the “smoothScatter” function. This function created a 2-D density plot for all cows, each day by working out the frequency of points in a local neighborhood, following kernel density smoothing. The density values were then plotted using a standard graphic plot.

Statistical analyses were conducted in RStudio© (v1.2.5019) (Rstudio, 2015), an integrated development environment for R (v4.0.2) (R Core Team, 2020). Two data sets were used in the analysis: (1) stimuli count data of the AT, EP, and proportion of EP/AT, and the (2) locational data.

Stimuli count data were analyzed using two negative binomial generalized linear mixed models (GLMM) using the “lme4” package (Bates et al., 2015) with AT count and EP count as the response variables. Within these models, the fixed effects included day (1–10) with Cow ID as the random effect. A negative binomial model was used to account for the over dispersion of the count data. The proportion of EP to AT, i.e., the number of ATs resulting in EPs being delivered, was analyzed using a logistic GLMM with the fixed effects of Day (1–10) and Cow ID as the random effect. The cow standard deviation estimate of random effects (Cow ID) and its corresponding 95% confidence interval (CI) from the three models was used to assess the variation between cows in stimuli delivery.

All the locational data were analyzed with linear mixed models (LMM) using the “lme4” package. The duration of bouts and the speed within each zone was analyzed with two LMM with the fixed effects of zone (IZ, EZ, BZ), day (1–10) and time of day (0–23 h) including interactions, with a random effect of zone nested within cow ID. Nested random effect allowed between-cow differences to be explored both in overall bout duration, but also bout durations of individual cows across the three zones. Analysis was performed on log-transformed data to account for the positively skewed data and residual variance. For the speed LMM a constant of 0.05 m/s was added prior to transformation due to stationary bouts. The bout physical location was then plotted for each cow per day using a standard graphic plot.

Daily total duration in the IZ was analyzed by day (1–10), with cow ID as a random effect. The IZ was considered due to the data constraint, with majority of the observed time being within the IZ. To analyse distance into the EZ, a subset of the EZ data was used, with the fixed effects of day (1–10), time of day (0–23 h) with cow ID as a random effect. Analysis was performed on log-transformed data to account for the positively skewed data and unstable residual variance.

All model-based means for fixed effects were obtained using the “emmeans” package (Lenth, 2021) and pairwise comparisons using the cld function in the “multcomp” package (Hothorn et al., 2008).

The spatial utilization of areas throughout the IZ and spatial coordinate data were manipulated using the “sp” (Pebesma and Bivand, 2005; Bivand et al., 2013) and “geosphere” (Hijmans, 2019) packages. The GPS data have been analyzed on a per day basis, and records were filtered to only those within the inclusion area of the day. Each inclusion area was divided into a grid of 0.1 s (longitude) × 0.1 s (latitude), resulting in grid pixels of 1.63 × 3.09 m. This resulted in between 198 and 317 pixels per study period, depending on the inclusion area size. The number of observations within each pixel were calculated over each 1-day study period. Next, these counts were converted into proportions of observations within each of the n pixels for that study day, p_i, i = 1, …, n. From these pixel proportions, the Shannon diversity index (H) was calculated as

Values of H increase with increasing number of pixels, but more importantly with more even distribution across the pixels. It may be considered as a measure of variability for categorical data, like the role of the standard deviation for quantitative data. The maximum diversity occurs when all p_i are equal, i.e., all equal to 1/n, and this results in H_max = log_en. As the number of grid pixels differ across the 10 study days (i.e., across the 10 inclusion areas), it is preferable to compare spatial utilization using evenness, defined as E = H/H_max. A value of E = 1 would indicate the cows use all parts of the inclusion area equally. Bootstrap standard errors of H and hence E were calculated by performing 1,000 re-samples form a multinomial distribution [N = total number of observations p = (p₁, …, p_n)] and re-calculating H and E at each sample, with the SE being the SD of these 1,000 values.

The mean ± SE estimated pasture consumption on hold-off days was 165 ± 0.0 kg DM/ha as compared to 675 ± 60.5 kg DM/ha for the grazing days.

There was no significant difference in the number of ATs delivered on the hold-off days as compared to the grazing days. There was a difference in the number of ATs delivered across days (P = 0.002, Table 1). Less ATs were delivered on Day 1 as compared to Days 2–9. The StdDev of the random effect of cow was identified with a CI of 0.35–0.89. The CI range was above 0 identifying significant variation between individual cows in the number of ATs delivered across days.

There was no difference in the number of EPs delivered on the hold-off days as compared to grazing days. However, there was a difference between days (P = 0.002, Table 1). The lowest number of EPs were delivered on Day 1 as compared to Days 2–10. The StdDev of the random effect of cow was identified with a CI 0.60–1.58 indicating variation between individual cows in the number of EPs delivered across days.

There was no difference in the proportion of EP/AT between hold-off days and grazing days. There was a difference across days (P = 0.005, Table 1). Day 1 had lower EP/AT as compared to Days 4, 5, 8, 9, and 10. The StdDev of the random effect of cow was identified with a CI of 0.28–0.84, indicating variation between individual cows in EP/AT across days.

Density plots of the herd location for each day of the experiment were created, to visualize where cows spent most of their time in each allocation. We have included plots from grazing Day 4 and the subsequent hold-off Day 5 to illustrate the similarity in herd density across days on the same pasture allocation (Figure 4). The darker the color represents a greater density of data points.

Cow location by hour was plotted for every cow across all days using the bout per hour model. The GPS location by hour is visually presented as the location of one individual (cow 3) within the same allocation across Days 4 and 5 (Figure 5).

The evenness values, used to evaluate spatial utilization of the IZ, were high across all days with an average of E = 0.941 and SE = 0.0003, with variation due to the large numbers of observations. A value of E = 1 indicates cows used all parts of the IZ equally.

There was no significant difference in the time spent in the EZ on the hold-off days as compared to the grazing days. The daily duration in the EZ across the experiment was low, with cows spending an average of 20 ± 5 min (range 4–44 min) in the EZ within each 24-h allocation (Figure 6). Time spent within the BZ changed significantly across the study (P < 0.001). Time spent in the BZ decreased from Day 6 in conjunction with increased time in the IZ (Figure 6). Cows spent more time in the IZ on the hold-off days as compared to the grazing days (P = 0.05). There were differences between individual animals in time spent in the IZ (average 21.5 ± 0.377 h, minimum 19.0 h and maximum 22.8 h, P < 0.001).

There were significant differences in bout duration (time between zone changes) across zones and days (zone × day P < 0.001). There was an increase in bouts in the IZ across the study period (Figure 7A), with more consecutive time spent in the IZ, and less transition into other zones. There was also a significant zone × hour interaction (P < 0.001), with peaks in time spent in the IZ between 01:00 and 02:00 and 17:00 and 18:00 h but without any corresponding peaks in activity in the other two zones (Figure 7B). There was significant variation in bout duration between cows (P < 0.001).

The distance that cows traveled into the EZ reduced across the 10 days of the experiment (P < 0.001, Figure 8A), with cows moving <8 m into the EZ from day 5. Distance traveled into the EZ varied across hours of the day (P < 0.001), with cows moving a greater distance into the EZ in the 2-h period from 06:00 to 08:00 h (Figure 8B).

There was a significant zone × hour interaction for speed of cow (P < 0.001), with speed above 0.1 m/s between 07:00 and 18:00 h. The greatest speed occurred in the EZ, followed by the BZ and then the IZ (Figure 9). There was overall variation in speed between cows (P < 0.001).