Influence of environmental factors and parity on milk yield dynamics in barn-housed dairy cattle

We investigated the effects of environmental factors on average daily milk yield and day-to-day variation in milk yield of barn-housed Scottish dairy cows milked with an automated milking system. An incomplete Wood gamma function was fitted to derive parameters describing the milk yield curve including initial milk yield, inclining slope, declining slope, peak milk yield, time of peak, persistency (time in which the cow maintains high yield beyond the peak), and predicted total lactation milk yield (PTLMY). Lactation curves were fitted using generalized linear mixed models incorporating the above parameters (initial milk yield, inclining and declining slopes) and both the indoor and outdoor weather variables (temperature, humidity, and temperature-humidity index) as fixed effects. There was a higher initial milk yield and PTLMY in multiparous cows, but the incline slope parameter and persis-tency were greatest in primiparous cows. Primiparous cows took 54 d longer to attain a peak yield (mean ± standard error) of 34.25 ± 0.58 kg than multiparous (47.3 ± 0.45 kg); however, multiparous cows yielded 2,209 kg more PTLMY. The best models incorporated 2-d lagged minimum temperature. However, effect of temperature was minimal (primiparous decreased milk yield by 0.006 kg/d and multiparous by 0.001 kg/d for each degree increase in temperature). Both primiparous and multiparous cows significantly decreased in day-to-day variation in milk yield as temperature increased (primiparous cows decreased 0.05 kg/d for every degree increase in 2-d lagged minimum temperature indoors, which was greater than the effect in multiparous cows of 0.008 kg/d). Though the model estimates for both indoor and outdoor were different, a similar pattern of the average daily milk yield and day-to-day variation in milk yield and milk yield’s dependence on environmental factors


INTRODUCTION
Dairy business success depends on milk yield and quality.Milk yield and the duration of the lactation of a dairy cow is not only influenced by its genetic capacity, but also environmental factors.Climate change is a concerning issue faced by farmers in many regions of the world, especially in the animal production industry.Though dairy cattle are able to cope with a wide range of climatic conditions, sustainability in the performance of these animals is potentially challenged by climate change, especially as there has been selection for larger breeds such as Holstein (Bos taurus).Meteorological conditions such as ambient temperature, radiant energy, photoperiod, relative humidity (RH), and wind speed contribute to the degree of heat stress events that occurs for the cow (Hammami et al., 2013).Moreover, meteorological conditions, mainly temperature, have a direct influence on biological functions; ultimately, they have a negative influence, including, but not limited to, DMI, rumination, milk production, and reproductive performance in cows (West, 2003;Nardone et al., 2010;Polsky and von Keyserlingk, 2017).Climate change is characterized by higher temperatures and extreme weather events, which may have a direct negative effect on dairy milk production.At high temperatures and humid climates, animals encounter challenges to balance the metabolic heat production and dissipate enough heat to their surroundings.
Heat stress occurs when the animal experiences a condition or state whereby it cannot dissipate enough heat to maintain thermal balance (Bernabucci et al., 2014).Dairy cows use behavioral and physiological strategies to cope with heat stress, including reduced feed intake and taking more water to decrease metabolic heat production, and ultimately milk yield is reduced (Herbut et al., 2018).Therefore, it is important for the dairy cow to be within its thermoneutral zone, which is a range of temperatures that the cow doesn't have to increase energy expenditure to maintain internal body temperature stable.The comfortable environmental temperature for the dairy cow is between 5°C (lower critical temperature) and 25°C (upper critical temperature) (Kadzere et al., 2002).However, the upper critical temperature varies depending on the stage of lactation as the lactation curve dynamics change, degree of acclimatization, RH, and milk production rate in animals (Aharoni et al., 2003;Fuquay, 1981).Therefore, a bioclimatic index, the temperature-humidity index (THI), has been used as a tool to quantify these thresholds.A THI above 72 has been reported as the threshold to induce heat stress in the tropics (Ravagnolo and Misztal, 2000;Kadzere et al., 2002), whereas in temperate zones, milk yield in high yielding cows can be affected at lower a THI value of 60 (Brügemann et al., 2012;Gorniak et al., 2014).
There is some evidence that cows that are exposed to even moderate temperature coupled with high humidity are also exposed to heat stress because it influences the rate of heat loss through evaporative cooling.In the United Kingdom, 90% of the dairy herd is composed of Holstein-Friesians, which have high body temperature due to high milk production (West, 2003;Chebel et al., 2004).There is some evidence that animals in European countries will encounter heat stress, as climate change will bring hotter summers and an elevated risk of extreme weather events (Novak et al., 2009;Hill and Wall, 2017).The ambient temperatures in the United Kingdom are projected to increase up to about 3.5°C by the 21st century (Foskolos and Moorby, 2018).In southwest England, United Kingdom, an average annual milk loss of 2.4% (~170 kg/cow) due to heat stress is projected per annum (Foskolos and Moorby, 2018).By the year 2080, southern Scotland temperatures are predicted to increase in summer by 4.3°C coupled with a slight reduction in humidity ranging between 0 and 5% (Jenkins et al., 2009).
Daily milk yield is not static throughout lactation but forms a very distinctive lactation curve.Temperature and humidity are by far the most critical factors having a major negative effect on livestock productivity (Rojas-Downing et al., 2017), welfare, and health, but it is not clear how they influence each aspect of the typical cow lactation curve.Therefore, as climate change increases Scottish temperatures, it is important to understand the effects it may have on milk production, particularly the dynamics of the lactation curve, to be better able to mitigate its effects on milk production.

Study Location
The study was carried out in a commercial family-run dairy farm, located near Inverurie in Aberdeenshire, Scotland, United Kingdom.The climate is classified as temperate oceanic climate with an annual average temperature of 8.1°C and rainfall of 762 mm (Climatedata.org,2019).

Animals and Housing
All cows were kept in a loose-housing system (barn) on concrete slatted floors with sand-bedded stalls all year round with no access to pasture.All lactating cows were milked by a voluntary robotic milking system supplied by Lely (Astronaut A4, Lely Industries N.V.).Cows were able to decide their own activity (time spent lying, standing, eating and drinking) and visits to the milking system.The building had walls of slatted wood with large doors at each end and an open ridge in the roof to facilitate airflow.The herd consisted of 7 breeds, namely Jersey, Aberdeen Angus, Holstein-Friesian, Holstein-Friesian X, Ayrshire, Belgian Blue, and Swedish Red X, with the majority being Holstein-Friesian (75%).Only Holstein-Friesian and their crosses were used in this study.All cow breeds were fed the same feed, with the amount varying depending on their production stage.Cows evaluated in this study were born from 2008 to 2015 and calved from 2010 to 2019.Cows ranged between their first and eighth lactations, and they were all managed (management and feeding routines) together in the same way.Each cow was identified by Lely QWES-HR tags (electronic transponder) and numbered neck collars, which were recognized in the automatic milking system (AMS).The transponder ensured that the details of all milking and trafficking events were electronically recorded by the Time for Cows (T4C) herd management software, version 3.12.37.603.These collar transponders collected data on individual performance (milk yield, fat, protein content), activity, and rumination time.
The barn was composed of 3 pens, and cows only had access to the AMS associated with their pen.The holding capacity per pen was 110 cows on average.The pens were cleaned by the automated alley scrapers twice a day and additional sand (bedding) was added daily to the cleaned stalls as needed by staff.Cows shared water troughs placed in each pen.Each pen had self-grooming brushes positioned in the crossover zones connecting the main walkways to enhance their grooming time (Mandel et al., 2016).Cows were dried off 2 mo before calving, and gestation lasted on average 283 d.After parturition, the newborn calves were removed within 12 to 18 h following birth, and the lactating cow was prepared for milking.The group structure was dynamic due to year-round calving, with cows entering and leaving depending on calving and drying off dates, in addition to sale or culling.

Feeding and Milking Management
Cows were fed roughage twice daily at 0900 and 1400h at a feed barrier so that it was available ad libitum.In addition, they received pelleted concentrate feed during milking as a reward, depending on their expected level of daily milk production, and also had access to a salt lick stone.Cows were milked at least twice per day with a minimum milking interval time of 2.5 h between visits, and most of the cows were limited to 5 to 6 milkings per day.Some exceptional cases had unrestricted access due to indications of mastitis.

Animal Data Sources
The T4C software was associated with the robotic AMS, and it stored milk production records and cow profile for each milking visit.The downloaded initial data set consisted of 217,406 historical daily milk yield records (observations) of 356 lactating dairy cows between the years 2010 and 2019.Records of lactations ≥5 were excluded due to limited and missing data.All analyses in the current observational study were focused on the daily observations and recordings rather than individual visits.The data set was edited to include records limited to a maximum of 305 DIM, and data points beyond 305 DIM were dropped.In the data set, lactations were analyzed for parity categorized as either primiparous (first lactation) or multiparous (≥2 lactations) to explore how it affected milk production and the characteristics of the lactation curve (Nasri et al., 2008).Parity was coded as either primiparous or multiparous because these 2 parity groups showed different lactation curve shapes, with parity 1 group having a much flatter curve than 2+; the clear distinc-tion of these curves is known between parity 1 and ≥2 but not between parity 2 and 3+ (Andersen et al., 2011;Jingar et al., 2014).The daily milk yield records were matched with the weather data [indoor (microclimate) and outdoor] from the same day.The records for reproductive status were obtained as the separate file from the T4C database for each individual cow, and the records contained the insemination date and day.The cows were inseminated between 60 and 80 d following calving, and pregnancy was later confirmed.

Weather Data Sources
Meteorological Climatic Data.Data on local weather conditions including daily ambient temperature (°C) and humidity (%), over the study period were downloaded from the website www .wunderground.combased on the data from the nearest weather station, which was at Aberdeen Airport, Dyce (57°12′18″N, 02°10′36″W).These data consisted of the 24-h summaries (average, minimum, and maximum values).The weather station was situated approximately 32 km from the farm.

Indoor (Barn Microclimate).
A HOBO microstation (H21-USB data logger) was used in the barn in August 2018 on the walkway directly above the cow pens (4 m).The micro-station recorded temperature, humidity, and light intensity within the barn every 10 min.Due to limited historical weather data within the barn, the daily (24-h summaries) temperature and humidity data from the micro-station and Dyce weather station on the same date were compiled to allow reconstruction of the microclimate within the barn for the previous years, given the known outside temperature and humidity from the weather station.The relationship between the barn interior and exterior data was matched within 10 min and plotted for both temperature and RH; the plots revealed 2 inflection points at low and high data points in all cases.Therefore, a segmented linear model was used to account for the nonlinear nature of the relationship between barn microclimate and outdoor environmental weather data (minimum, average, and maximum values).The models were then used to predict temperature for all historical dates in the study.Daily THI (minimum, average, and maximum) both indoor (microclimate) and outdoor were calculated using the equation by National Research Council (1971) as follows: where T db is the daily temperature (°C).
The weather data were matched with the daily milk yield records from the same day.Weather can have a delayed effect on the biological processes because of the amount of time the ruminant animal takes to fully consume feed, digest, and metabolize the nutrients (Gauly et al., 2013;Hill and Wall, 2017), and the effect of the weather depends on the duration over which the animal experiences it (West, 2003;Renaudeau et al., 2012).Therefore, we explored the relationship between the milk yield and weather variables (temperature, RH, and THI) on the day of and 1, 2, and 3 d before the dependent variable measurements (test day) were determined.We called these effects lagged effects.For example, −3-d minimum temperature was the minimum temperature taken 3 d before the current daily milk yield.We created the lagged (−3, −2 and −1) variables for the weather data, including temperature (minimum, maximum, and average values), RH, and THI for both indoor (microclimate) and outdoor from the test-day data.

Statistical Analysis
The reduced data set after exclusions contained 125,460 daily milk records of lactating cows.The descriptive statistics of the data set can be found in Tables 1 and 2. Modeling the environmental influence on milk production was carried out using statistical computing language R (https: / / www .r-project .org/ ) in Rstudio version 1.2.5033 and LME function in the NLME package (Pinheiro et al., 2018).Statistical analysis was carried out to model milk yield dynamics using the incomplete Wood gamma model (Wood, 1967(Wood, , 1969)), which relates milk yield to DIM using an exponentially truncated power law as follows: where Y t is the average daily milk yield (dependent variable, kg/d) on the nth day (DIM), a is the scaling factor associated with the initial milk yield (kg/d) at the beginning of the lactation (after calving), b is the parameter that informs the rate of increase in milk yield up to the peak, c is the parameter that informs the rate of decrease in milk yield after the peak, and e is the neper number.The parameters b and c have positive values specific to each typical lactation curve (Wood, 1967).Thus, cows that yielded negative lactation curve parameters (a, b, and c) were excluded from the analyses to avoid atypical lactation curves.Therefore, the final data set contained 125,460 test-day milk yield records of 163 lactating cows and 449 lactation cycles.These model parameters were estimated for each cow per lactation curve by the non-least squares procedure using the nlsList function (nlme package; Pinheiro et al., 2018).After derivation, the estimated model parameters became the fixed effects.These were then used to model the lactation shape curve traits for specific lactation number of each cow, namely, DIM at peak yield (DIMpeak = b/c, days), peak milk yield and persistency [S = −(b + 1) × ln(c)], which is the rate at which milk yield declines after the peak (Tekerli et al., 2000).The area under the lactation curve represented predicted total lactation milk yield (PTLMY) up to 305 DIM.

How Does Weather Influence the Parameters of the Lactation Curve?
To determine the effects of weather variables on the average daily milk yield, first, the incomplete gamma Wood function was transformed logarithmically (below) into linear form (fixed effect equation; Wood, 1969) because the nonlinear regressions do not guarantee convergence.When fitted to daily milk yields records per lactation where DIM ≤305 d, the equation is as follows: We used a general linear mixed-effects model (GLMM) to fit the model to the final data set [ln (yield)] using Equation [3].We then generated a total of 6 possible initial models (GLMM) that explained the lactation curve and fitted to the logarithmically transformed Wood model (Equation [3]).The predictors of interest (fixed effects) included incline slope parameter b, decline slope parameter c, parity (primiparous vs. multiparous cows), and insemination time (day), as well as the interaction terms between the parameters and parity.To account for the interindividual variation, random intercept-only models were fitted with the random effects of the animal identity within each lactation.The models were fitted with the temporal first-order autocorrelation (autoregressive model with a lag of 1 d) structure to account for the non-independence of the milk yield values between DIM.The random components of the models were assumed to be normally distributed with zero mean and variance.The parameters for separate linear mixed-effects models were estimated using the maximum likelihood technique.We hypothesized that insemination day might influence the peak milk yield.Therefore, to account the effect of insemination on the milk yield, the insemination day variable, which corresponds to the insemination day for each cow cycle, was created to determine if it did influence milk yield and incorporated into the model.We used the lowest Akaike information criterion (AIC) value to select the best model.
To determine the effects of weather variables on the lactation curve dynamics, we incorporated the continuous effect of both indoor (within the barn) and outdoor (meteorological) daily minimum, maximum, and average values of temperature, humidity, and THI and the lagged effects as the fixed effect into the selected GLMM initial model.We then generated a total of 288 possible GLMM models comprising 144 for the indoor (barn) and 144 for the outdoor weather data based on the possible variations between the combined parameters used in the models and were compared using the maximum likelihood technique.All these models were generated from the initial model selected.The best model that explained the variation in milk yield was selected based on the lowest AIC score.Two final GLMM models that best explained our data were selected for the indoor and outdoor weather effect.The 2 final models were then re-fitted by restricted maximum likelihood technique to minimize bias of the parameter estimates.The coefficients of the final models were re-fitted to the Wood model for the prediction estimation at different values for both indoor and outdoor weather effects.Figure 1 illustrates the model fitting and selection procedure for the study.
How Does Weather Influence the Day-to-Day Variations in Milk Yield?The Wood model was fitted to predict the average daily milk yield (kg/d) specific to the individual cow for each day of lactation by using the estimated model parameters (a, b, and c) and DIM.Model residuals were generated from the difference between the observed and predicted daily milk yield.Analysis of the model residuals was then used to study the effect of daily weather variation on daily milk production.Studying the effect of the weather on daily milk yield through model residuals removes any confounding effect of the DIM (Ehrlich, 2013).
To determine the effect of weather on the day-to-day variation in milk yield, a total of 72 multiple linear regression models were generated based on the possible variations between the combined independent variables used in the model and fitted using the lm function with the model residuals (deviations) as the dependent variable.Fixed effects were weather variables [indoor and outdoor minimum, maximum, and average temperature; RH; THI; and the lagged weather effects (1, 2, and 3 d) and parity (primiparous vs. multiparous cows) or weather variable × parity interaction terms].We fitted the following model: where Y is the response variable (day-to-day variation in milk yield, kg/d), μ is the overall mean for each cow on the given DIM, W is the single weather variable or weather variable plus weather variable × parity interaction term, and ε is the residual error term.Models were compared using AIC, and models (indoor and outdoor) that best explained the day-to-day variations in milk yield were chosen based on the lowest AIC value.The model residuals were assessed for normality and homogeneity of variance in R using diagnostic plots (histogram, residual plots against the fitted values, and quantile-quantile plots).In the present study, all the model results are reported as estimates (β) with associated standard errors, t-values (df), and the P-values.Results were considered significant at P < 0.05.value, was chosen.The first aim was to investigate the weather effects on the lactation curve dynamics over 305 DIM.A total of 6 initial models were fitted to the Wood model.In the second step, weather variables were incorporated as covariates in the model selected in the step 1 process.The second aim was to investigate the weather effects on the day-to-day variation in the average daily milk yield.A Wood model was fitted to the daily milk yield records, and model residuals were generated (the difference between predicted and observed milk yield values).These model sets, including the weather variables and parity as covariates, were used in the model.GLMM = general linear mixed-effects model; ML = maximum likelihood; REML = restricted maximum likelihood.T4C = herd management system linked to the Lely automatic milking system.

Microclimatic Versus Meteorological Conditions During the Experimental Periods
The mean seasonal weather conditions for both indoors (microclimate) and outdoors are presented in Table 3.During the study period, the absolute lowest and highest outdoor temperature was −8.00 and 27.00°C, respectively, although the mean outdoor temperature conditions in winter and summer seasons ranged from 0.89 ± 0.02 to 7.05 ± 0.02°C and 10.07 ± 0.01 to 17.53 ± 0.02°C, respectively.A greater mean indoor value of RH was observed (83.95 ± 0.01%) compared with the outdoor (79.31 ± 0.03%), whereas daily mean indoor THI was higher with mean of 51.31 ± 0.02, which was 2.73 points higher than mean outdoor THI (48.58 ± 0.02).Greater mean ambient temperature, RH, and THI values were observed within the barn than outdoors.Mean ambient temperature in winter and summer was 2.38 and 1.52°C higher within the barn compared with outside temperature, respectively.Also, RH in winter was 4% higher inside the barn, whereas mean THI in summer was 2.36 units greater within the barn compared with the outdoor THI values.In the current study, we found that the indoor minimum and maximum temperature, RH, and THI values ranged between −2.6 and 33.4°C, 32.8 and 96.9%, and 29.4 and 89.4,respectively.Moreover, the outdoor minimum and maximum weather variables, temperature, RH, and THI ranged between 0.89 and 17.53°C, 30.09 and 95.15%, and 38.83 and 63.39, respectively.

Seasonal Effect on the Actual Average Daily Milk Yield
There was a slight but significant difference in the actual average daily milk yield (kg/d) between seasons (F = 3.488, P = 0.015), accounting for the individual differences.Least square averages showed that cows produced 0.207 ± 0.071 kg/d (P = 0.003) higher actual average daily milk yield in winter (35.733

Lactation Curve Parameters and Traits
The average estimated lactation curve parameters and traits including initial milk yield (a), incline slope parameter (b), decline slope parameter (c), peak milk yield (Ymax, kg/d), DIM at peak yield (DIMpeak, d), persistency (S), and PTLMY (kg) and were calculated from the incomplete Wood gamma model fitted to each lactation per cow.The results are presented in Table 4.Our study indicated that multiparous cows had greater initial milk yield (22.7 ± 0.40 vs. 12.7 ± 0.44 kg/d) and higher Ymax average values (47.3 ± 7.76 vs. 34.1 ± 7.19 kg/d) compared with the primiparous cows.Despite multiparous cows having the highest parameter a, greater peak milk yield, and earliest day of peak production (b/c), parameter b (0.30 ± 0.17 vs. 0.25 ± 0.09) and lactation persistency (7.78 ± 0.73) were greater in primiparous cows.Primiparous cows took longer to attain the peak yield (116.6 ± 3.80 d) compared with multiparous cows (60.3 ± 1.08 d).Moreover, primiparous produced 2,209 kg less PTLMY compared with the multiparous cows.

Weather Effect on the Lactation Curve Dynamics in Lactating Dairy Cattle
The GLMM analysis was carried out to investigate the effect of weather on the lactation curve dynamics.Weather factors had significant effects on milk yield.Log milk yield was best explained by the 2-d lagged minimum temperature (indoor = −2d.minTb;outdoor = −2d.minT).There was no influence of maximum temperature, humidity, and THI values.Two final models were selected, one for each of the indoor (barn) and outdoor weather data based on the lowest AIC score.The results of both models demonstrated that the effects of the 2-d lagged minimum ambient temperature on the log milk yield measurements improved the fit of the model compared with that of the current day measures or 1-d lag or 3-d lagged climatic variables (temperature, RH, and THI).
Both retained models (indoor and outdoor) had similar significant fixed covariates, which were 2 threeway interaction terms [b × parity (primiparous, multiparous) × 2-d lagged minimum temperature and c × parity × 2-d lagged minimum temperature] and a 2-way interaction (b × parity).The summary of the models' (indoor and outdoor) coefficients for the explanatory variables are demonstrated in Table 5, and the model estimates were re-fitted to the Wood model to make predictions (Figures 2 and 3).
Outdoor Weather Influence on Lactation Curve Dynamics.The −2d.minT ranged between −8 and +17°C.The model estimates for the outdoor weather effect on the log milk yield were slightly different from that of the indoor conditions; however, a similar trend of milk output and its dependence on the environmental factors was observed for both primiparous and multiparous cows (Table 5, Figure 3).Primiparous cow milk production was most affected by the change in 2-d lagged minimum temperature (for both indoor or outdoor) compared to multiparous cows, but the effect was minimal.
Indoor (Barn Microclimate) Weather Conditions' Influence on Lactation Curve Shape.During the study period, the indoor minimum temperature ranged between −2.7 and 19.6°C and was 2.6°C greater than the outside minimum air temperature.Our results indicated that the main effect of the incline slope, decline slope parameters, and −2d.minTb on the log daily milk yield were significant depending on the parity group (primiparous vs. multiparous).The model results are presented in Table 5.
Our model also revealed that log daily milk yield was also influenced by 2 significant 3-way interaction terms between b × parity (primiparous, multiparous) × −2d.minTb (F = 5.392, P = 0.020) and c × parity × −2d.minTb (F = 8.156, P = 0.004).This revealed that the average log daily milk yield and parameters b and c were influenced by −2d.minTb, but this effect depended on the parity group.Again, we found that it was primiparous cows that had the greatest effect of parameter b on the average log milk yield, but the effect was minimal at only 0.002 kg/d (0.2%) for each degree (°C) increase in the −2d.minTb (t = 2.322, P = 0.020).Moreover, the association between parameter c and −2d.minTb when cows were multiparous had no significant effect on milk yield after the peak was reached (β = 1.3 × 10 −6 , t = 0.329, P = 0.742).However, when cows were primiparous, there appeared to be a small but significant decline, but it was only 1.83 × 10 −5 kg/d in milk yield for each degree (°C) increase in the −2d.minTb (t = −2.856,P = 0.004).
Overall, we had no evidence that −2d.minTb affected average daily log milk yield and parameters b and c of the multiparous cows.However, the decrease in the average log milk yield in both primiparous and multiparous cows was associated with an increase in −2d.minTb (Figure 2).

Weather Effects on the Day-to-Day Variation in Milk Yield
The goal of this analysis was to investigate the effect of weather (indoor and outdoor) on the day-to-day variation in milk yield (kg/d) rather than the effect on the overall aspects of the lactation curve dynamics.Day-to-day variation in milk yield was significantly influenced by the effect of lagged 2 minimum tempera-ture depending on the parity group.The results of both retained models are presented in Table 6.

DISCUSSION
In this study, we investigated the effect of weather on the lactation curve dynamics in Scottish lactating dairy cattle.To understand the milk yield dynamics, we fitted an incomplete Wood gamma function with 3 parameters to the historical milk yield records to derive parameters describing the milk yield curve and incorporated the weather elements.
There has been much effort invested in investigating the effect of climate change on animal productivity and animal welfare.A better understanding of the effects of indoor and outdoor weather events will potentially permit farmers to develop better management systems (Gauly et al., 2013).In the United Kingdom, climate change is also a challenge, especially on dairy cow performance on a pasture-based system (Foskolos and Moorby, 2018).The present study investigated the  1 Intercept = multiparous cows was the reference category; −2d.minT/b = minimum temperature (for both indoor, −2d.minTb, and outdoor, −2d.minT) 2 d before milk yield measurements were taken (°C); −2d.minT = minimum outdoor temperature 2 d before milk yield measurements were taken (°C); −2d.minTb = minimum indoor temperature 2 d before milk yield measurements were taken (°C).
effects of weather on the lactation curve in dairy cows housed indoors in Scotland.Similar to other areas of Northern Europe, indoor temperatures in this study were 3 to 5°C higher than the outdoor temperature; however, RH differed depending on the external temperature (Seedorf et al., 1998).Microclimatic temperature increase can also be influenced by the stocking density, solar radiation, and wind speed (Haider et al., 2017).Several studies report the upper critical temperature of Holstein ranges from 25 to 26°C (Berman et al., 1985;Kadzere et al., 2002;West, 2003), whereas Angrecka and Herbut (2015) report a lower critical temperature of −6.7°C.However, there is a lot of inconsistency in terms of the lower critical temperature below which the cows begin to experience cold stress.
For instance, Kadzere et al. (2002) reported a range of −37 to −16°C for cows to produce a daily milk yield of 30 kg, whereas Brouček et al. (1991) found a reduction of 2 kg of milk yield in cows yielding 15 kg/d on average at a temperature below −10°C.In the current study, the upper critical temperature of +25°C outdoor was observed for only 2 d per year, and the lower temperature of −8.0°C was observed once a year.According to West (2003), a thermoneutral zone ranging between −0.5 and +20°C is acceptable for dairy cows, hence its effect is small to induce milk production and behavioral changes among cows.Dairy cows tend to lie more at low temperatures compared with temperatures >20°C (Vaculíková et al., 2017); ultimately, feed intake is reduced, which affects milk yield.In contrast, the study of Brouček et al. (1991) reported a decline in average milk yield as the minimum temperature approached as low as −19°C in heifers coupled with an increase in feed intake.Dairy cows can acclimatize to changeable temperatures and humidity levels (Kadzere et al., 2002).Several factors influence the thermoneutral zone range (lower and upper critical temperature) of the dairy cow including but not limited to breed, age, previous temperature acclimatization, production level, and feed intake (NRC, 1981;Kadzere et al., 2002;Polsky and von Keyserlingk, 2017).Our data show that, at present, cows in Scotland do not regularly experience ambient temperatures above their upper critical temperature, and hence are not often heat-stressed, and variations in ambient temperature are unlikely to have a significant negative effect on their milk production.Therefore, this could lead to a minimal income loss of approximately £2.36 ($3.26) per cow per lactation based on the current milk farm-gate price of £0.30 ($0.41) per liter (Department for Environment, Food and Rural Affairs, 2021).

Lactation Parameters and Traits
In the present study, multiparous cows had the highest initial milk yield (a), greater peak milk yield (Ymax), total milk yield and earliest day of peak production (b/c); however, persistency (S) and parameter b were greatest in primiparous cows.The difference observed in the initial milk yield and Ymax may be due to the fact that heifers (primiparous cows) are bred at ≤24 mo old, and they are still growing at this stage.The nutrients obtained from feed have to be partitioned between growth and milk production, causing competition for the available resources (Wathes et al., 2007).Milk production, in general, improves with the advancement in lactation number (Vijayakumar et al., 2017) and an increased number of mammary epithelial cells (Herve et al., 2016).The lower initial milk yield in the primiparous cows and increase in the multiparous cows has also been reported in the previous research (Atashi et al., 2013;Jingar et al., 2014;Masía et al., 2020;Sitkowska et al., 2020).
On average, primiparous cows had a greater b parameter than multiparous cows, and this indicated that primiparous cows had a slower increase in milk production to peak after calving.These results are in agreement with the findings of other authors (Bangar and Verm, 2017;Sitkowska et al., 2020).For instance, Sitkowska et al. (2020) found greater b parameter in the first lactation than the second lactation in Polish Holstein-Friesians cows using AMS data; however, Masía et al. (2020) reported similar parameter b values of 0.23 and 0.24 for the primiparous and multiparous cows, respectively.
Multiparous cows produced 2,208.9kg (PTLMY) of more milk across the 305 d than primiparous cows.This is also in agreement with earlier studies (Bagnato and Oltenacu, 1994) that indicated that cows of different parity but same age produce different milk levels.Milk yield increased with the advancement in parity or lactation number (Lee and Kim, 2006), and it also correlates with maturation and increased body size, and mammary gland development with an increasing number of the secretory cells (Wood et al., 1980;Nwosu et al., 2019).In addition, the mammary gland grows together with recurring pregnancies and lactations.
Despite the higher PTLMY in multiparous cows, greater lactation persistency was observed in primiparous cows.The current finding was similar to that reported by Ehrlich (2013).The persistency values of 7.78 ± 0.06 and 6.86 ± 0.02 were found for primiparous and multiparous, respectively.These values are in line with that reported by Dematawewa et al. (2007); Cole and Null (2009); Atashi and Asaadi (2019) but greater than that of Kellogg et al. (1977), who reported values of 3.13, 2.29, and 2.22 for the first-, second-, and thirdlactation Holstein cows, respectively.Lower persistency (greater value of parameter c) observed in multiparous cows is attributed by the age factor because as an animal gets older, mammary alveolar cell regression occurs at a rapid rate, leading to a reduction in milk production (Jingar et al., 2014).It has been reported that high persistency is associated with the slow rate of decline after the peak yield, whereas low persistent lactation is related to high rate of decline due to reduced feed intake (Sölkner and Fuchs, 1987).Efficient feed utilization and less metabolic stress at the peak milk yield have a major influence on the high persistency in lactation observed in primiparous cows (Cole and Null, 2009).Moreover, the higher persistency may be caused by the relatively lower peak yield observed in primiparous cows in the present study, resulting in a flatter lactation curve.These cows additionally tend to have more negative energy balance postpartum than do multiparous cows and may fail to show estrus until the energy balance is more favorable (Reksen et al., 1999).

Indoor and Outdoor Weather Effects on Lactation Curve Dynamics and Day-to-Day Variation in Milk Yield
The influence of the indoor and outdoor minimum, average, and maximum temperatures; RH; THI; and lagged effects (up to 3 d before) were examined with respect to milk yield on both the whole lactation curve and also individual day-to-day variation.When investigating the whole curve dynamics, the 2 models for both −2d.minTb and −2d.minT provided the best fit of the milk yield data, and thus had the greatest influence.The lagged effects of environmental conditions on milk yield were confirmed by the findings of previous studies (Bouraoui et al., 2002;West, 2003;Bertocchi et al., 2014) that observed that the greatest effect of weather variables occurred up to 3 d before test day on the milk yield and DMI.These findings reflect the amount of time the ruminant animal takes to fully consume feed and digest and metabolize the nutrients (West, 2003;Gauly et al., 2013;Hill and Wall, 2017).It seems that unlike previous studies, maximum temperature, humidity, and THI did not explain milk production, and that minimum temperature was much more influential.Our results, however, are discrepant with the study of Gorniak et al. (2014), who noted negative effects on milk yield with THI >60 in temperate conditions.Similarly, recent studies in the United Kingdom found that milk yield and composition decreased with increased THI (Dunn et al., 2014;Hill and Wall, 2015), depending on whether the cows were housed indoor (barn) and outdoor on pasture (Hill and Wall, 2015).The study of Hill and Wall (2015) reported a decrease of 0.02 kg/d in milk yield for the indoor-housed dairy cows in Scotland.In the United Kingdom, the robust heatwaves were reported in 2003 and 2006; however, the studies of Dunn et al. (2014) and Hill and Wall (2015) did not reveal a significant milk yield decrease as a result of heat stress, and they stated an average of 1 d of heat stress each year.This is not surprising because Fodor et al. (2018) also reported a relatively low heat stress occurrence rate in UK dairy cows in the 2010s, with the predicted average annual milk yield loss of 1 kg/cow in the north of the United Kingdom, indicating no evidence of milk yield penalty with the predicted average of 0.6 d with the THI above 70.Moreover, the study of Fodor et al. (2018) indicated that the average annual milk yield loss in the 2090s in the United Kingdom is predicted to be relatively low (2.4%), including the southeast England region.This is the region that is predicted to have an average of 51 d with THI exceeding 70.
A limitation of these previous studies is that the influence of barn microclimate on milk production was not explored.Including barn microclimate would have provided a better understanding of how dairy cow performance is influenced by immediate environment in the temperate maritime climate region.Barn microclimate has been noted to have a significant effect on animal welfare and milk yield and quality in Holstein dairy cows (Vaculíková et al., 2017).
The present study, however, found similar patterns for both the indoor and the outdoor weather effect on the lactation curve dynamics.Both models retained had slightly different estimates but indicated that the effect of weather on milk yield was influenced by if the cow was primiparous or multiparous.Our models indicated that an increase in both −2d.minTb and −2d.minT had an influence on the lactation curve traits (i.e., parameters a, b, and c) and average daily milk yield but depended on the parity of the cows.Surprisingly, the effect of temperature change was extremely low, having minimal effect on the day-to-day milk production of the cows.A decrease of 0.6% for each degree (°C) increase in the −2d.minTb on milk production was observed in primiparous cows, which, if we were to consider the projected increase in temperature of 4.3°C coupled with a slight decrease in humidity ranging between 0 and 5% by 2080 (Jenkins et al., 2009) caused by climate change, would only amount to a 7.87-kg reduction in milk production over a 305-d lactation.The reduction in primiparous cows, although significant, was not biologically important.In multiparous cows, there was no effect at all.
As well as examining the influence of changes in weather conditions on the whole lactation curve, the present study also investigated the effects of weather on the day-to-day variation in milk yield by determining if changes in temperature caused deviations from the predicted lactation curve generated using the Wood's model.However, several factors influence milk production including but not limited to parity, stage of lactation, milking interval (Everett and Wadell, 1970;Syrstad, 1977;Quist et al., 2008), and environmental factors such as temperature (Čandek-Potokar et al., 2006).Understanding the day-to-day variations in milk yield measurements are crucial to monitor herd management improvement (Čandek-Potokar et al., 2006).Similarly, the main findings of this study revealed that day-to-day variation in milk yield was also negatively influenced by both the minimum indoor and outdoor temperature 2 d before the test day.The magnitude of the effect depended on the parity group (primiparous or multiparous cows).These findings are in line with that of the study of Kabuga and Sarpong (1991) but contradict the findings of Gabbi et al. (2017), who highlighted a reduction in milk yield in both primiparous and multiparous cows at maximum temperature.On average, primiparous cows yielded 0.29 kg/d more day-to-day variation in milk yield than did multiparous cows (0.10 kg/d).The less day-to-day deviation in milk yield in multiparous cows is an indication that these cows were less influenced by temperature than primiparous cows.
The findings of this study showed that both primiparous and multiparous experienced a decrease of 0.04 ± 0.005 kg/d and 0.01 ± 0.003 kg/d, respectively, for every 1°C increase in −2d.minTb (illustrated in Figure 4A).This indicated that the barn microclimate affects milk production (Vaculíková et al., 2017).Despite a myriad number of publications reporting a decrease in milk yield at minimum temperatures (MacDonald and Bell, 1958;Brouček et al., 1991), others observed 1239 similar findings as to our study whereby the daily milk yield significantly declined with increasing temperatures (Barash et al., 2001;Yano et al., 2014).
As for the effect of the outdoor weather on the dayto-day variation in milk yield, the same trend as the barn was found.Interestingly, the minimal changes in milk yield and day-to-day variation in milk yield observed in multiparous cows indicated that the cows were able to maintain their body temperature within the −2d.minTb and −2d.minT range that we observed in our study, which enabled them to sustain their performance.The primiparous cows were not so robust to change.During milk synthesis, metabolic heat production increases, which negatively influences feed intake, consequently reducing subsequent milk production level (Kadzere et al., 2002;Yano et al., 2014).Multiparous cows tend to be heavier than primiparous cows, and larger animals have more difficulty in dissipating heat due to their small surface to body ratio (Speakman and Król, 2011).However, this might also be an advantage when the daily temperature fluctuations are below the lower critical temperature.Overall, both analyses of lactation curve dynamics and daily variations showed that predicted milk yield decreased with an increase in both −2d.minTb and −2d.minT for both primiparous and multiparous cows, but at a low level that may not affect farming practices.

CONCLUSIONS
In this study, we tested the effects of indoor (barn) and meteorological temperature, RH, and THI on milk production.Although the model estimates for both indoor and outdoor were slightly different, a similar pattern of milk output and its dependence on environmental factors was observed for both primiparous and multiparous cows.The statistical results indicated that increases in 2-d lagged minimum temperature had a negative effect on milk yield.Primiparous cows were more affected than multiparous cows.However, the effect size was small.The results of the current study indicated that it is important to also make predictions based on the microclimatic conditions, in addition to the meteorological predictions, to better understand the effect that weather has on milk yield in dairy cattle in Scotland.Our results showed that the farmer should implement a suitable strategy or solution to protect the animals, particularly when the microclimatic temperatures are low.As much as the management strategies to deal with heat stress have been investigated, the results of the present study indicated that it is also important for the farmer or breeder to manipulate the microclimate to enhance animal performance during the period of low temperatures.

Figure 1 .
Figure 1.Flow diagram of the statistical model selection procedure.The best model, based on the lowest Akaike information criterion (AIC)value, was chosen.The first aim was to investigate the weather effects on the lactation curve dynamics over 305 DIM.A total of 6 initial models were fitted to the Wood model.In the second step, weather variables were incorporated as covariates in the model selected in the step 1 process.The second aim was to investigate the weather effects on the day-to-day variation in the average daily milk yield.A Wood model was fitted to the daily milk yield records, and model residuals were generated (the difference between predicted and observed milk yield values).These model sets, including the weather variables and parity as covariates, were used in the model.GLMM = general linear mixed-effects model; ML = maximum likelihood; REML = restricted maximum likelihood.T4C = herd management system linked to the Lely automatic milking system.

1a
= initial milk yield at the beginning of the lactation in kg (after calving); b = increasing slope parameter up to the peak milk yield; c = declining slope parameter after the peak; ADMY = average daily milk yield; DIMpeak = days to the peak (b/c); Ymax = peak milk yield [a(b/c) b e −b ]; persistency = rate at which milk yield declines after the peak [−(b +1) × ln(c)]; PTLMY = predicted total lactation milk yield over 305 d.
Marumo et al.: ENVIRONMENTAL FACTORS ON PRODUCTION IN DAIRY COWS Journal of Dairy Science Vol.105 No. 2, 2022 1233 Marumo et al.: ENVIRONMENTAL FACTORS ON PRODUCTION IN DAIRY COWS Table 5.The best general linear mixed-effects models summary results of the weather effects of each parameter on the average log milk yield (kg/d) of both model 1 (indoor) and 2 (outdoor) in 163 Holstein cows (n = 125,460) during the experimental study period (2010cows was the reference category; b = inclining slope parameter; −2d.minT/b = minimum temperature (for both indoor, −2d.minTb, and outdoor, −2d.minT) 2 d before milk yield measurements were taken (°C); −2d.minT = minimum outdoor temperature 2 d before milk yield measurements were taken (°C); −2d.minTb = minimum indoor temperature 2 d before milk yield measurements were taken (°C); c = declining slope parameter.2 Estimate = regression estimate.

Figure 2 .
Figure 2. Predicted average daily milk yield (kg/d) at different temperature levels of the indoor (barn) minimum temperature for primiparous and multiparous cows.The black dots represent the actual observed daily milk yield.

Figure 3 .
Figure 3. Predicted average daily milk yield (kg/d) at different levels of the outdoor minimum temperature for primiparous and multiparous cows.The black dots represent the actual observed daily milk yield.

Figure 4 .
Figure 4.The relationship between day-to-day deviation in milk yield (kg) and 2-d lagged (A) indoor and (B) outdoor minimum temperature by parity group (illustrated by line colors; green and orange as primiparous and multiparous cows, respectively).The regression lines were generated using a linear model from the ggscatter function in the ggpubr package(Kassambara, 2020).

Table 1 .
Marumo et al.: ENVIRONMENTAL FACTORS ON PRODUCTION IN DAIRY COWS Descriptive statistics of the milk yield in relation to lactation number (n = 125,460 d of milking data)

Table 2 .
Descriptive statistics of the milk yield in relation to the parity group (n = 125,460 d of milking data) 1 N = number of daily observations of milk yield.

Table 3 .
± 0.481 kg/d) than in autumn (35.526 ± 0.481 kg/d) but did not differ from summer and spring yields.No significant difference was found between actual average daily milk yield in summer (35.533 ± 0.481 kg/d) and spring (35.648 ± 0.481 kg/d).et al.: ENVIRONMENTAL FACTORS ON PRODUCTION IN DAIRY COWS Descriptive statistics of the mean seasonal weather conditions (2010-2019) indoors (microclimate) and outdoors during the study period; values are expressed as mean ± SE 1 Summer = June-August; autumn = September-November; winter = December-February; spring = March-May.2 N = 125,460.

Table 4 .
Incomplete Wood gamma model parameters estimated (means ± SE) values for the typical lactation curve parameters and traits for primiparous and multiparous cows(n = 125,460)