Exploring the presence of genotype-by-environment interactions between dairy cow herds milking once-a-day or twice-a-day for the entire lactation

Abstract

The objective was to explore if the regression of phenotypic performance for six milk production traits on the respective estimate of genetic merit for that trait differed by herd milking frequency; variance components for each trait in the two milking frequency environments were also estimated as well as the genetic correlation between the same trait in both environments. The data used included 12,581 lactations from 5,456 cows in 32 spring-calving once-a-day (OAD) milking herds. Each OAD herd was matched with three contemporary twice-a-day (TAD)-milking herds; 35,823 lactations from 15,188 cows in 96 TAD herds were used. Mean yield was 20% (fat yield) to 31% (milk yield) lower in OAD herds. Milk protein concentration was 11% higher in OAD herds, while milk fat concentration was 16% higher in OAD herds. The mean back-transformed somatic cell score (SCS) was 100,390 cells/mL in OAD herds and 72,493 cells/mL in TAD herds. The association between each milk production trait and its respective estimate of genetic merit differed by herd milking frequency; the estimated regression coefficients were larger in TAD for just milk yield and SCS. The genetic correlation between the same trait in OAD versus TAD was all ≥0.73 with the exception of SCS (genetic correlation of 0.48) which suggests some re-ranking of sires between environments. In conclusion, differences in the scale of the genetic variance were evident for both milking frequencies and possible re-ranking was evident for SCS.

Main article text

Introduction

Dairy cows in most jurisdictions are milked, on average, twice-a-day (TAD; Stelwagen et al., 2013). Some, however, milk cows just once-a-day (OAD) for either part of the lactation (Kay et al., 2012; Phyn et al., 2014; Kennedy et al., 2021) or the whole lactation (Lembeye et al., 2016a; Edwards, 2020). Several motivatory factors influence the decision to embark on a strategy of milking cows just OAD, especially in systems where maximising yield per cow is of lesser importance relative to other considerations (Stelwagen et al., 2013). Milking OAD in early lactation is often used as a tactic to reduce the extent and duration of negative energy balance (Kay et al., 2012; Phyn et al., 2014), with an anticipated benefit for cow reproductive performance and metabolic health status (Patton et al., 2006; McNamara et al., 2008). Lifestyle choices or shortage of labour also contribute to the decision to milk cows just OAD (Edwards, 2020). Periods of feed shortages could also force producers to milk cows just OAD as a tactical option to reduce herd feed demand. What is not clear, however, is if the breeding policy of producers operating an OAD-milking system should differ from those milking TAD. The phenomenon of where differences between genotypes of animals (i.e., breeds, strains of animals or indeed individual animals or families) vary between environments is termed a genotype-by-environment (G×E) interaction (Falconer, 1989). Two possible manifestations of G×E exist: (1) re-scaling—when the scale of the difference between genotypes varies between environments without impacting the relative ranking of the genotypes or (2) re-ranking—when the correct ranking of genotypes is dependent on the environment (Lynch & Walsh, 1997).

Extensive research has been undertaken, primarily in New Zealand, on the performance and merits of OAD and TAD milking (Clark et al., 2006; Tucker et al., 2007; Kay et al., 2012) as well as how breeds or animals may rank differently depending on herd milking frequency (Lembeye et al., 2016b). Controlled studies have also been undertaken in France (Rémond et al., 2004; Pomiès et al., 2007), Germany (Schlamberger et al., 2010) and Ireland (Gleeson et al., 2007; Kennedy et al., 2021). Yields of milk, fat and protein are, on average, lower in OAD herds, but milk concentration of both fat and protein tends to be higher in OAD herds (Lembeye et al., 2018; Edwards, 2019). Milk somatic cell count (SCC) is also generally higher, on average, in OAD- versus TAD-milked herds (Lembeye et al., 2018; Edwards, 2019). Whether differences exist in the genetic variance for milk production traits between herds differing in milking frequency has also been documented to vary not only by breed but also by parity (Lembeye et al., 2018). Despite the often-cited differences in variance between herd milk frequency sub-populations, strong genetic correlations between OAD- and TAD-milking herds have been documented for milk production traits (≥0.67; Lembeye et al., 2018); nevertheless, few have actually attempted to estimate the genetic correlations. The difference in genetic variance by herd milking frequency in the presence of strong genetic correlations between milking frequency environments suggests a re-scaling G×E exists in some situations (i.e., particular breeds and parities) but without necessarily a significant re-ranking effect. Whether G×E exists for milk production in Irish herds milked OAD versus TAD, however, has never been explored using a cross-sectional analysis of national data.

Therefore, the objectives of the present study were to (1) describe the performance and breeding policy characteristics of Irish OAD-milked herds and compare them to matched TAD-milked herds of similar geographical location, herd size, age structure and calving system and (2) explore the presence of G×E for milk production traits between OAD- and TAD-milking herds. Results will be useful for producers, breeders and advisors in supporting breeding policy decisions for OAD-milking herds in Ireland.

Materials and methods

All data used in the present study originated from the Irish Cattle Breeding Federation (www.icbf.com) which manages the national database; only data between the years 2019 and 2022, inclusive, were considered. A list of dairy farmers operating an OAD farming system for the entire lactation was obtained from a Teagasc advisor (with farmer consent); no information was available on how many years they had been practicing full lactation OAD milking prior to the study period. The data available from the database included milk production lactation yields (including mean SCC), calving dates, individual animal breed composition, parity number and ancestral information. Individual animal predicted transmitting abilities (PTA) for 305-day milk yield, fat yield, protein yield, milk fat concentration, milk protein concentration and SCC were obtained from the December 2018 national genetic evaluation. The estimated breeding value (EBV) of all animals for each milk production trait was generated as the sum of the respective sire and dam PTA. National genetic evaluations for milk, fat and protein in Ireland are undertaken using an across-breed test-day model where parity 1, 2 and 3+ are treated as separate traits in a multi-trait model. Milk, fat and protein genetic evaluations are run separately. The published PTA is a weighted average of the individual parity PTAs as 0.41:0.33:0.26 for parity, 1, 2, and 3+, respectively. Breed is fitted via genetic groups in the founder population which are allocated to breeds. The genetic evaluations for SCC are undertaken using an animal repeatability model where log_eSCC is the dependent variable.

Contemporary TAD-milking herds

A total of 32 spring-calving herds known to milk cows just OAD throughout the whole lactation and that also participated in milk recording were available. Where possible, all OAD herds had to participate in milk recording for each of the 4 years of the study period but this was relaxed due to the impact of COVID-19 restrictions on milk recording participation. Nonetheless, 69% of the OAD herds milk recorded in all 4 years of the study between 2019 and 2022, inclusive. All herds had >40 calving events in each year.

An algorithm was developed to locate at least three unique matching TAD herds for each of the 32 OAD-milking herds. All TAD herds had to have participated in milk recording for each of the 4 years of the study period. Herds were matched based on the county in Ireland, herd size, average calving date, and proportion of the herd that was primiparous. Herd size of the TAD herds had to be within 10 cows of its matching OAD herd, with the percentage of primiparous cows in the TAD herds having to be within ±5 percentage units of its OAD match; the mean calving date of the TAD herds had to be within 10 days of the mean of the matching OAD herd. The restriction on herd size was relaxed from being within 10 cows to being within 15 cows for 12 of the 32 OAD herds so that at least three TAD herds could be identified which also fulfilled the other matching criteria. A random sample of just three unique TAD herds for each OAD herd was subsequently retained.

Data edits

The logarithm to the base 10 of SCC per lactation was calculated to normalise the distribution and will here forward be referred to as somatic cell score (SCS). Only animals with a recorded sire and dam were retained. The parity number was collapsed into 1, 2, 3, 4, 5+. Because of the number of both Holstein–Friesian crossbreds and Holstein–Jersey crossbreds in the dataset, three separate heterosis coefficients were calculated for each animal: (1) Holstein×Friesian (HO×FR), (2) Holstein×Jersey (HO×JE) and (3) residual heterosis defined as a general heterosis coefficient less the sum of the HO×FR and HO×JE heterosis coefficients. The HO×FR coefficient was calculated as:

$$\text{HO}\text{}\times \text{}\text{FR\hspace{0.17em}=\hspace{0.17em}}{P}_{\text{HO}}^{\text{Sire}}{P}_{\text{FR}}^{\text{Dam}}\text{+}{P}_{\text{FR}}^{\text{Sire}}{P}_{\text{HO}}^{\text{Dam}}$$

where P represents proportion, the subscripts HO and FR represent the proportion of Holstein and Friesian, respectively, and the superscript represents whether the proportion of the breed (subscript) relates to the sire or the dam. The HO×JE heterosis coefficient was calculated the same way except that FR was replaced by JE (i.e., Jersey). The general heterosis coefficient per animal was calculated as:

$$\text{heterosis}=1-{\displaystyle \sum _{i=1}^n{\text{sire}}_i\times {\text{dam}}_i}$$

$$\text{recombination\hspace{0.17em}loss}=1-{\displaystyle \sum _{i=\text{1}}^n\frac{{\text{sire}}_i^2+{\text{dam}}_i^2}{2}}$$

where sire _i and dam _i are the proportion of breed i in the sire and dam, respectively.

Each lactation record was allocated to a within-herd contemporary group based on the date of calving. The approach used to generate a contemporary group was that used for most of the national genetic evaluations in Ireland (Berry et al., 2013). Only contemporary groups with at least seven records were considered further where the difference in calving date between the start and end of the contemporary group was no longer than 30 days. Following all edits, 12,581 lactations from 5,456 cows in 32 OAD herds remained while 35,823 lactations from 15,188 cows in 96 TAD herds remained. The ancestry of each animal was traced back to the founders which, in turn, were allocated to genetic groups based on breed. The pedigree consisted of 104,756 animals.

Statistical analyses

The association between phenotypic milk performance and the EBV for that respective trait was quantified using linear mixed models in SAS (SAS Institute Inc., Cary, NC, USA). Parity within cow was included as a repeated effect in all models with a first-order autoregressive covariance structure with heterogeneous variances assumed among errors per parity. The contemporary group was included in all models as a random effect. Fixed effects considered in the model included parity (categorical), HO×FR heterosis coefficient (continuous), HO×JE heterosis coefficient (continuous), residual heterosis coefficient (continuous), milking frequency (categorical) and EBV (continuous) for the dependent variable being investigated. Two-way interactions between each of the model features and milking frequency were also tested with a particular focus on the two-way interaction between EBV and milking frequency.

Genetic, permanent environmental and residual variance components for each of the six milk production traits were estimated for OAD and TAD herds separately using a repeatability animal linear mixed model in ASREML (Gilmour et al., 2009) based on the model:

$$\begin{array}{l}{Y}_{ijk}{\text{=\hspace{0.17em}CG}}_j\text{}\text{+\hspace{0.17em}}{H}_{\text{HO}\text{}\times \text{}\text{JE}}\text{+\hspace{0.17em}}{H}_{\text{HO}\text{}\times \text{}\text{FR}}\text{+\hspace{0.17em}}{H}_{\text{Residual}}{\text{+\hspace{0.17em}parity}}_k\\ {\text{+\hspace{0.17em}animal}}_i{\text{+\hspace{0.17em}PE}}_i\text{+\hspace{0.17em}}{e}_{ijk}\end{array}$$

where Y_ijk is the milk trait under investigation (milk yield, fat yield, protein yield, fat percent, protein percent, SCS), CG _j is the jth contemporary group of herd-year-season of calving, H _HO×JE and H _HO×FR are the heterosis coefficient between Holstein and Jersey and between Holstein and Friesian, respectively, H _Residual is the general heterosis coefficient less the previous two heterosis coefficients, parity _k is the parity k of the cow i, a_i is the random additive genetic effect of animal i ( $a\hspace{0.17em}~\hspace{0.17em}N(Qg,A{\sigma }_a^2)$ , where Q is a matrix relating a to genetic groups of breed, g is a vector of genetic group means, A is the numerator relationship matrix and ${\sigma }_a^2$ is the additive genetic variance), PE _i is the random permanent environmental effect of animal i across lactations ( $N(0,I{\sigma }_{PE}^2)$ where ${\sigma }_{PE}^2$ is the permanent environmental variance and I is the identity matrix) and e_ijk represents the residual term where $N(0,I{\sigma }_e^2)$ with ${\sigma }_e^2$ representing the residual variance and I an identity matrix. The genetic covariance between the same trait in both milking frequency environments was estimated using bivariate repeatability animal linear mixed models with the same model construction as used in the univariate analyses. No residual covariance was assumed since no animal had a record in both milking frequency environments.

Results

Summary statistics

The raw mean and s.d. for the yield and composition traits as well as SCS for OAD and TAD herds are in Table 1. Mean yield was 20% (fat yield) to 31% (milk yield) lower in OAD herds compared to TAD while little difference was detected in the s.d. of the yield traits between both milking frequency environments. Milk protein concentration was 11% higher in OAD herds, while milk fat concentration was 16% higher in OAD herds relative to their TAD contemporary herds; the observed variation in milk composition was larger in the OAD herds. The mean back-transformed SCC was 100,390 cells/mL in OAD herds and 72,493 cells/mL in TAD herds with similar variability.

Table 1:

Raw 305-day lactation mean (raw s.d. in parenthesis) for cows milk OAD or TAD

Trait	OAD	TAD
Milk yield (kg)	4,580.0 (1,193.0)	6,626.6 (1,371.3)
Fat yield (kg)	220.8 (55.03)	276.1 (54.29)
Protein yield (kg)	183.2 (45.61)	240 (47.01)
Fat percent (%)	4.89 (0.72)	4.20 (0.47)
Protein percent (%)	4.03 (0.30)	3.64 (0.23)
SCS (log10 units)	5.00 (0.384)	4.86 (0.400)

OAD = once-a-day; TAD = twice-a-day.

The mean PTA and breed composition of the sires of all cows in the OAD and TAD herds, weighted by number of daughters in each, are in Table 2. Also included in Table 2 is the mean EBV and breed composition of the cows in both herd types. A difference (P < 0.001) between both milking frequencies existed for all traits for both sires and cows. The mean genetic merit for milk yield for both the sires and cows was lower in the OAD herds compared to the TAD herds, but the opposite was true for fat and protein yield, particularly fat yield. This manifested itself as the mean genetic merit for milk fat and protein percent also being higher in OAD herds. Minimal biological differences existed in genetic merit for SCS of the two milking frequency groups. Jersey bloodlines were more prevalent in the OAD herds.

Table 2:

Mean (s.d. in parenthesis) genetic merit and breed proposition of the sires of the cows and the cows themselves in the OAD- and TAD-milking herds

	Sires		Cows
	OAD	TAD	OAD	TAD
Genetic merit
Milk yield (kg)	−6.2 (192.10)	73.8 (194.15)	−38.7 (254.34)	103.1 (276.12)
Fat yield (kg)	13.4 (7.65)	10.1 (8.10)	20.1 (11.97)	14.5 (12.20)
Protein yield (kg)	8.4 (5.29)	7.9 (6.50)	11.7 (8.19)	11.1 (9.45)
Fat percent (%)	0.23 (0.182)	0.12 (0.139)	0.37 (0.232)	0.18 (0.191)
Protein percent (%)	0.14 (0.079)	0.09 (0.068)	0.22 (0.109)	0.13 (0.097)
SCS (log10 units)	−0.02 (0.090)	−0.03 (0.086)	−0.03 (0.108)	−0.04 (0.105)
Breed proportion
Jersey	0.36 (0.386)	0.06 (0.197)	0.28 (0.240)	0.04 (0.120)
Holstein	0.41 (0.328)	0.63 (0.314)	0.47 (0.235)	0.65 (0.233)
Friesian	0.20 (0.220)	0.29 (0.279)	0.18 (0.175)	0.26 (0.216)

OAD = once-a-day; TAD = twice-a-day.

Of the 1,603 sires with daughters in the entire dataset, 416 (i.e., 26%) had daughters in both OAD and TAD herds; a total of 84 sires had at least 10 daughters in each of OAD and TAD herds. Of the 20,644 individual cows in the dataset, 14,954 had at least one paternal half-sib in the other milking frequency group.

Regression analyses

The association between herd milking frequency and each milk production trait differed (P < 0.05) by parity (Table 3). The difference in yield between the two milking frequencies relative to the TAD yield reduced with each progressing lactation from a 37% to 38% reduction in parity 1 to a 27% to 30% reduction in parity 3 and a 22% to 26% reduction in parity 5+. The heterosis effect for the yield traits was greater in OAD than in TAD (Table 3). The association between EBV and the respective phenotype differed (P < 0.05) by milking frequency for all traits (Table 3). The linear regression of phenotypic milk yield on milk EBV was 0.98 in OAD herds and 1.14 in TAD herds. The trend in the association between EBV and phenotype was opposite for fat and protein yield compared to milk yield in that the average response per unit change in EBV was 15% (protein yield) to 37% (fat yield) greater in OAD herds compared to TAD herds. A similar trend was observed for fat and protein percent where the response per incremental difference in the respective EBV was 25% (protein percent) to 37% (fat percent) greater in OAD compared to TAD herds. A greater reduction in phenotypic SCS was, on average, detected per unit decrease in EBV for SCS in TAD herds; nonetheless, the sign of the observed association was in line with expectations.

Table 3:

Milk production predicted marginal means for parity by milking frequency (assuming value of 0 for all covariates) for OAD- and TAD-milking herds along with the linear regression coefficients of each trait on HO×FR, HO×JE and the respective EBV for the trait under investigation

Trait	Parity					HO×FR	HO×JE	EBV
	1	2	3	4	5+
Milk yield (kg)
OAD	3,205 (46.0)	4,270 (46.4)	4,732 (47)	5,065 (47.9)	5,227 (46.1)	145.3 (68.65)	277.6 (43.05)	0.98 (0.032)
TAD	5,152 (24.7)	6,258 (24.8)	6,744 (25)	7,023 (25.4)	7,042 (23.8)	90.1 (29.29)	19.5 (41.3)	1.14 (0.019)
Fat yield (kg)
OAD	133 (2.2)	179 (2.2)	199 (2.2)	214 (2.2)	220 (2.1)	6.5 (2.95)	17.1 (1.84)	1.29 (0.03)
TAD	209 (1.1)	250 (1.1)	272 (1.1)	283 (1.1)	282 (1.1)	4.0 (1.25)	12.3 (1.75)	0.94 (0.02)
Protein yield (kg)
OAD	111 (1.7)	154 (1.7)	173 (1.8)	185 (1.8)	190 (1.7)	5.1 (2.38)	15.8 (1.48)	1.24 (0.04)
TAD	176 (0.9)	219 (0.9)	237 (0.9)	246 (0.9)	245 (0.9)	3.8 (1.02)	7.5 (1.42)	1.07 (0.02)
Fat percent (%)
OAD	4.60 (0.02)	4.37 (0.02)	4.31 (0.02)	4.30 (0.02)	4.27 (0.02)	−0.093 (0.036)	−0.001 (0.022)	1.50 (0.022)
TAD	4.06 (0.01)	3.97 (0.01)	3.99 (0.01)	3.98 (0.01)	3.98 (0.01)	0.028 (0.015)	0.158 (0.022)	1.10 (0.014)
Protein percent (%)
OAD	3.72 (0.009)	3.70 (0.009)	3.69 (0.009)	3.67 (0.009)	3.63 (0.008)	−0.034 (0.015)	0.042 (0.010)	1.53 (0.021)
TAD	3.43 (0.004)	3.47 (0.004)	3.48 (0.004)	3.48 (0.004)	3.46 (0.004)	0.017 (0.006)	0.105 (0.009)	1.22 (0.013)
SCS (log10 units)
OAD	5.03 (0.01)	4.92 (0.01)	4.99 (0.02)	4.99 (0.02)	5.05 (0.01)	0.061 (0.031)	0.011 (0.019)	0.38 (0.032)
TAD	4.86 (0.01)	4.79 (0.01)	4.84 (0.01)	4.88 (0.01)	5.01 (0.01)	−0.030 (0.013)	−0.018 (0.018)	0.48 (0.019)

EBV = estimated breeding value; HO×FR = Holstein×Friesian heterosis coefficient; HO×JE = Holstein×Jersey heterosis coefficient; OAD = once-a-day; TAD = twice-a-day.

Variance components

Genetic and phenotypic s.d. as well as the heritability and repeatability estimates for the six evaluated milk traits in OAD and TAD herds are in Table 4. The genetic s.d. for milk yield in OAD herds was 88% that of TAD herds equating to a difference of 58.7 kg. In contrast, the genetic s.d. for fat yield was 32% greater in OAD than in TAD herds. There was minimal difference in genetic s.d. for protein yield between the two milking frequencies. Large differences in genetic variance were evident for milk fat and protein percent with the genetic s.d. in OAD herds being 2.73 to 2.95 times that estimated in TAD herds. The genetic s.d. for SCS was lower in OAD herds. Minimal difference between milking frequency environments was evident for either the heritability or the repeatability estimates for any of the milk traits evaluated. With the exception of SCS, the genetic correlation between each trait in OAD herds with its respective trait in TAD herds varied from 0.73 to 0.98 (Table 4); the respective genetic correlation for SCS was 0.48.

Table 4:

Genetic and phenotypic s.d. as well as heritability and repeatability estimates (standard errors in parenthesis) for all milk production traits in herds milked OAD or TAD; also included is the estimated genetic correlation (r _g; standard error in parenthesis) between the same trait in both milk frequency environments

Trait	Genetic s.d.		Phenotypic s.d.		Heritability		Repeatability		r g
	OAD	TAD	OAD	TAD	OAD	TAD	OAD	TAD
Milk yield (kg)	417.9	476.5	782.9	835.3	0.28 (0.03)	0.31 (0.02)	0.60 (0.01)	0.63 (0.01)	0.73 (0.07)
Fat yield (kg)	22.6	17.1	39.7	33.4	0.32 (0.03)	0.25 (0.02)	0.57 (0.01)	0.52 (0.01)	0.86 (0.06)
Protein yield (kg)	16.2	15.6	30.6	27.5	0.28 (0.03)	0.27 (0.02)	0.59 (0.01)	0.60 (0.01)	0.73 (0.07)
Fat percent (%)	0.443	0.162	0.589	0.423	0.57 (0.03)	0.57 (0.02)	0.71 (0.01)	0.72 (0.005)	0.97 (0.02)
Protein percent (%)	0.204	0.069	0.251	0.204	0.66 (0.03)	0.69 (0.02)	0.82 (0.01)	0.80 (0.003)	0.98 (0.02)
SCS (log10 units)	0.120	0.161	0.367	0.360	0.11 (0.02)	0.13 (0.01)	0.39 (0.01)	0.33 (0.01)	0.48 (0.15)

OAD = once-a-day; TAD = twice-a-day.

Discussion

A myriad of different dairy production systems exist, not just globally, but even within the country. Ireland is no exception to the diversity in production systems—for example, seasonal spring-calving, seasonal spring- and autumn-calving, or year-round calving; producers of milk for manufacturing versus liquid consumption; OAD versus TAD. Understanding not only the characteristics of such production systems but also their breeding policy can help provide bespoke advice. Should a G×E interaction exist between production systems, then advice on animal selection should be tailored for the different production systems. While no official statistics exist for the number of Irish herds operating OAD for the entire lactation, they are very much in the minority representing probably <1% of the herd population.

G×E interactions have been the topic of research in animal studies for decades (for reviews see Montaldo, 2001; Dillon et al., 2006; Hammami et al., 2009). In dairy cattle, a plethora of environmental differentiation features have been explored including countries (Nilforooshan & Jorjani, 2022) or groups of countries (Zwald et al., 2003), herd milk yield (Berry et al., 2003; Calus and Veerkamp, 2003), herd mean SCC (Calus & Veerkamp, 2003), herbage quality (Berry et al., 2003) and concentrate feeding level (Berry et al., 2003). Especially when based on within-country analyses of dairy cow populations, re-scaling G×E is far more prevalent than re-ranking G×E (Berry et al., 2003; Calus & Veerkamp, 2003) with strong evidence of the latter rarely being evident. The greater the difference in the environments being compared, the greater the likelihood of detecting a G×E (Dillon et al., 2006). Furthermore, re-scaling can be accounted for in genetic evaluations by adjusting for heterogeneous variances (Meuwissen et al., 1996); heterogeneous variance adjustments for milk production are already undertaken as part of the Irish national genetic evaluation system. Nonetheless, even if just re-scaling exists for individual traits, re-ranking of animals for an overall breeding index which combines these individual traits is possible (Namkoong, 1985); the same is true for efficiency metrics which may be composites of traits (Charagu & Peterson, 1998). Nonetheless, Dickerson (1977) recommended that the genetic correlation between populations is the most informative criterion to quantify the impact of G×E. Irrespective, while some re-ranking of genotypes across environments is inevitable, justification for a separate genetic evaluation or a separate breeding programme only exists when large-scale re-ranking is occurring (among other reasons). In fact, Roberston (1959) stated that a separate breeding programme for different environments is only justifiable if the genetic correlation between the environments is weaker than 0.80. More recently, Mulder et al. (2006) stated that only if the genetic correlation between environments is <0.60, are different breeding programmes needed.

Herd characteristics of OAD versus TAD-milking herds

The mean statistics of the TAD-milking herds in the present study (Tables 1 and 2) should not be taken to represent the population average of all TAD Irish herds since the TAD herds used in the present study were chosen to mimic a matched case–control analysis to the 32 OAD herds. The statistical model used in the present study accounted for differences in age structure which could mask any potential benefit of improved reproductive performance (i.e., longevity) in OAD herds (Patton et al., 2006; McNamara et al., 2008) and its impact on age structure and thus achieving mature yields. Nonetheless, the greater impact of OAD milking on first parity cows observed in the present study has been reported elsewhere both in controlled studies (Rémond & Pomiès, 2005; Clark et al., 2006) and from cross-sectional analyses of large datasets (Lembeye et al., 2016b). In fact, the absolute difference in milk yield between OAD and TAD herds across parities represented approximately two phenotypic s.d. Information on the impact of OAD on milk fat and protein concentration is less well explored globally albeit, because yield data are provided in other studies, it was possible to calculate the impact of full lactation OAD on milk fat and protein concentration. Corroborating the results from the present study, calculations on the national data of OAD versus TAD dairy herds in New Zealand revealed that milk fat concentration and protein concentration were both 5% higher in cows milked OAD (when all breeds and parities were considered); this is lower than the 11–16% higher milk fat and protein concentrations observed in the present study in favour of OAD-milked cows. While the greater concentration of the milk could be partly attributed to a lower dilution effect with less yield in OAD cows, in the present study at least, approximately a quarter of this was due to breeding policy with the sires used in OAD herds having genetically greater milk fat and protein concentration. Nonetheless, the greater milk composition in OAD herds in the present study is consistent with those observed elsewhere in short-term studies (for review see Davis et al., 1999) or full lactation controlled studies (Rémond et al., 2004; Clark et al., 2006).

Because few studies have undertaken a cross-sectional analysis of commercial datasets on OAD and TAD-milking herds, there are limited data on the breeding decisions in both herd types. It should be noted, however, that while genetic merit and breed differences did exist among the different milking frequency environments in the present study, it is simply a statement of what occurred and does not imply that either breeding policy is optimum. This is emphasised by the fact that the genetic merit for SCS in the OAD cows was inferior to that of the TAD cows; inferior udder health (including SCS) is a known risk factor of OAD (Lacy-Hulbert et al., 2005; the present study; Lembeye et al., 2016a). There is a general tendency for SCS to be 2–3% higher in OAD-milked cows compared to TAD-milked cows (present study; Lembeye et al. 2016a). In order to negate the impact of OAD on high SCS, such producers should place more emphasis on SCS (and general udder health) when choosing parents of the next generation. In the present study at least, this was not the case with minimal difference in the SCS PTA of the sires of the cows in OAD versus TAD herds.

Edwards (2019) documented the mean genetic merit for eight different traits in OAD and TAD cows in New Zealand. Consistent with the observations in the present study, the mean genetic merit of the OAD cows for milk yield was lower than TAD cows (Edwards, 2019); in direct contrast to the New Zealand data (Edwards, 2019), however, Irish cows in OAD herds had a higher genetic merit for fat and protein yield as well as a higher (i.e., worse) genetic merit for SCS relative to their TAD contemporaries. The difference in trends of genetic merit between the two milking frequency environments for milk yield versus fat and protein yield is due to the greater emphasis obviously being placed on genetic merit for more concentrated milk composition in the breeding programmes of OAD herds.

Although the present study did not differentiate between purebred or crossbred status (<0.2% of the cows in the present study were ≥87.5% Jersey), the greater usage of Jersey bloodlines in the OAD herds investigated in the present study is consistent with that reported by Edwards (2019) in a cross-sectional analysis of New Zealand data for dairy herds milking OAD or TAD. The greater use of Jersey bloodline is most likely due to their documented greater suitability, as a breed average, to OAD-milking systems (Carruthers et al., 1993; Clark et al., 2006) in that the relative difference in yield between OAD and TAD herds is less for Jersey cows. The existence of more crossbreds in OAD herds in the present study could also have contributed to the lack of a perfect genetic correlation between milking frequency systems since crossbred performance is genetically different from purebred performance due to different genetic backgrounds (Wientjes & Calus, 2017). Breed was not included in the statistical model used in the present study to quantify the association between phenotype and EBV because the national genetic evaluations are across breeds and thus breed is already implicit within the EBV. Nonetheless, in a supplementary analysis, the regression coefficients of phenotype on its respective EBV did not change much when the breed was included in the model signifying that the national genetic evaluations are indeed capturing both within-breed and across-breed effects.

Evidence of genotype-by-environment interactions

With the exception of both milk fat and protein concentration, the variance component analyses failed to detect any real biologically important scale differences in genetic variance between OAD and TAD herds. It should, nonetheless, be noted that the OAD and TAD environments investigated in the present study did not differ solely on milking frequency but also on other associated differences in the production systems such as level of supplementation, stocking rate and other herd management features. Given the high heritability for milk traits, especially milk composition, a strong G×E is not expected as otherwise, the underlying variance due to G×E could partition some of the phenotypic variance away from the genetic variance, contributing to a lower heritability; the lower heritability of SCS in the present study could lend itself to be more susceptible to G×E. The genetic correlation between the same trait in different milking frequency environments was stronger than 0.72 for all traits except for SCS (0.48); nonetheless, the correlation of 0.48 between SCS in both OAD and TAD herds was associated with a relatively large standard error of 0.15 owing mainly to the lower heritability of SCS in both environments (0.11–0.13) relative to the yield traits (0.25–0.32) or milk composition traits (0.57–0.69). In fact, the only other study that has estimated the genetic correlation between milk production traits in OAD versus TAD herds also reported that the weakest correlation between the milk production traits was for SCS (Lembeye et al., 2018); the genetic correlations for SCC in the New Zealand study of 156,074 cows were as low as 0.67 depending on the parity number and breed while the weakest correlation for the yield traits was 0.77 with 83% of the 18 different yield traits evaluated by Lembeye et al. (2018) being >0.85 in strength. Similar to the present study, the standard error associated with genetic correlations was largest for SCS; the heritability of SCS in the New Zealand study (0.08–0.19) was often the lowest of all the reported heritability estimates for the milk performance traits across the different parity×breed combinations investigated. Nonetheless, the genetic correlations reported for SCS in the two milking frequency environments in New Zealand (Lembeye et al., 2018) and Ireland (the present study) were not weaker (i.e., P > 0.05) than the threshold genetic correlation of 0.60 suggested by Mulder et al. (2006) to be of practical significance.

The conclusion from the variance component analyses in the present study was generally corroborated by the regression analyses that regressed the phenotype on its respective EBV and tested if that regression coefficient differed by whether or not the cow was milked OAD or TAD. The expectation for the regression coefficient is one implying that a one unit difference in EBV is expected to translate to a one unit difference in phenotype—in a large cross-sectional analysis of the Irish dairy cow population, Ring et al. (2020) demonstrated how this expectation largely held true for most traits (when expressed at the level of EBV rather than PTA). If the regression coefficient is >1, then phenotypic change per unit change in EBV is, on average, greater than expected. The ratio of the regression coefficient of phenotype on EBV in OAD herds relative to TAD varied from 0.79 (SCS) to 1.37 (fat yield and concentration) indicating that while a statistically significant difference in regression coefficients by milk frequency existed (i.e., a scaling effect), its biological significance was questionable. No other study in dairy cows appears to have investigated the regression of milk phenotype on milk EBV in OAD versus TAD herds. Nonetheless, a linear regression coefficient can be calculated from the estimated correlation between two traits times the ratio of the s.d. of these traits. Using this approach it was possible to calculate the genetic regression of genetic merit for OAD on TAD herds based on the (co)variance components presented by Lembeye et al. (2018) for New Zealand dairy herds milked OAD or TAD. This could then be compared to the results from the present study derived using the same approach. The calculated regression for milk yield in the present study was 0.65 which is very similar to the average of 0.66 calculated from parity 1 and 2 cows across different breeds presented by Lembeye et al. (2018) using New Zealand data. In absolute terms, therefore, daughters of genetically superior bulls for milk yield are likely to suffer more than genetically inferior bulls. The genetic regression coefficient for fat yield in the present study was 1.13 which is higher than the average value of 0.74 based on the information provided by Lembeye et al. (2018) with a regression coefficient >1 only existing in that study for parity 1 Jersey cows. The regression coefficient of 0.83 for protein yield in the present study is very similar to the average value of 0.84 calculated from parity 1 and 2 cows across breeds in New Zealand (Lembeye et al., 2018).

Practical implications

The lack of any real strong evidence of re-ranking G×E between the OAD and TAD herds except for possibly SCS implies that a separate breeding programme for both herd types may not be justified; however, only milk production traits were investigated in the present study and all traits in the breeding programme as well as other features like the level of recording in each environment and the relative economic weight per trait in the breeding objectives need to be firstly determined before any definitive conclusion can be arrived at. Nonetheless, for Ireland at least, the population of OAD herds in Ireland would be too small to support an efficient and effective breeding programme. The expected annual genetic gain in OAD herds from a breeding programme optimised for TAD herds can be deterministically derived having modified the original genetic gain equation of Rendel & Robertson (1950) as:

$$\Delta G_{\text{OAD}}=r_{\text{g}}·i_{\text{TAD}}·r_{\text{TAD}}·{\sigma }_{\text{gOAD}}$$

where ΔG _OAD is genetic gain per generation in the OAD population, r _g is the genetic correlation between the two milking frequency environments, i _TAD is the selection intensity in the TAD herds, r _TAD is the accuracy of selection in the TAD populations, and σ _gOAD is the genetic s.d. of the trait/index in OAD herds.

Nonetheless, it could still be possible to express (national) genetic evaluations on an OAD scale by either a simple post hoc re-scaling of the EBVs generated using predominantly TAD-milking data from the national genetic evaluation or by including OAD data in a multi-trait evaluation for each trait. The accuracy of EBVs would simply be the accuracy of the EBVs from the TAD herds times the genetic correlations reported in Table 4 between the two environments. The benefit of such is expected to be negligible or even counter-productive relative to the additional computational requirements and potential for confusion of how and where to use the OAD and TAD EBVs, especially in herds that operate OAD for only a part of the lactation.

Given that, for some traits, a difference in genetic variance existed between OAD and TAD herds (e.g., milk yield) or indeed the genetic correlation between the same trait in both herd types differed from one (i.e., SCS), whether variance exists in how an animal’s EBV changes across milking frequency environments or, in fact, if this variation was heritable, was explored. The variance of the change between two points is simply the sum of the variance at each point, less twice the covariance between them. The genetic s.d. of the change in animal EBV between OAD and TAD herds for milk yield calculated in the present study was 329 kg; the respective heritability was 0.11. The heritability of the change in performance from OAD to TAD for the other five milk production traits investigated in the present study varied from 0.06 to 0.10. Hence, it could be possible to generate EBVs for sires based on the expected change in their daughters’ milk productivity should they transition from TAD to OAD. Such a trait could be considered in a bespoke OAD breeding objective or at least as a stand-alone trait. This supports the hypothesis of Lembeye et al. (2016a) who speculated that there could indeed be some genetic component to a cow’s ability to tolerate and adapt better to OAD. Having an EBV for this change could potentially be useful for herds that operate a OAD-milking system for a part of the lactation or indeed that might operate other reduced milk frequency systems such as milking three times every 2 days (Edwards et al., 2022). The approach proposed here would, nonetheless, still require a separate genetic evaluation for OAD and TAD production performance treating them as separate traits. Nevertheless, what the calculations from the estimated variance components in the present study do highlight is that considerable (genetic) variability does exist in how the genetic merit of an individual differs depending on if milked OAD or TAD.

The (co)variance components presented in the present study can also be used to estimate the genetic covariance (and by extension the genetic correlation) between performance in TAD herds and the change in EBV from TAD to OAD. Using the (co)variance components estimated in the present study, the genetic correlation between milk yield in TAD herds and the change in yield from TAD to OAD was 0.49 suggesting that animals genetically superior for milk yield in TAD suffered more in yield if only milked OAD. Interestingly, the genetic correlations for fat yield and protein yield were −0.19 and 0.10, respectively.

Conclusions

Comparison with controlled research studies is difficult since controlled studies, by their very name, are controlled for many factors while in reality, changes to the production system may be adopted by OAD herds in response to milking just OAD. Nonetheless, across recent studies, the reduction in milk yield in herds milking OAD for the entire lactation relative to TAD herds is in the region of 21–31% (the present study; Rémond et al., 2004; Clarke et al., 2016; Lembeye et al., 2016c); a reduction of between 17% and 29% is expected for fat or protein yield (Rémond et al., 2004), although the reduction does vary by parity. While G×E does exist for milk production traits both in the present study and elsewhere, the manifestation is generally as a scaling effect rather than a re-ranking effect; the one trait which seems most at risk of a re-ranking effect is SCS. Where re-ranking is expected to exist, the genetic correlation with TAD herds is still very strong, often greater than the required 0.60 to justify a separate breeding programme.

Conflict of interest

The authors declare no conflict of interest.

References

1. Berry DP, Buckley F, Dillon P, Evans RD, Rath M, Veerkamp RF. 2003. Estimation of genotype × environment interactions, in a grass-based system, for milk yield, body condition score, and body weight using random regression models. Livestock Production Science. Vol. 83:191–203. [Cross Ref]

2. Berry DP, Kearney JF, Twomey K, Evans RD. 2013. Genetics of reproductive performance in seasonal calving dairy cattle production systems. Irish Journal of Agricultural and Food Research. Vol. 52:1–16

3. Calus MPL, Veerkamp RF. 2003. Estimation of environmental sensitivity of genetic merit for milk production traits using a random regression model. Journal of Dairy Science. Vol. 86:3756–3764. [Cross Ref]

4. Carruthers R, Davis SR, Bryant AM, Henderson HV, Morris CA, Copeman PJA. 1993. Response of Jersey and Friesian cows to once a day milking and prediction of response based on udder characteristics and milk composition. Journal of Dairy Research. Vol. 60:1–11. [Cross Ref]

5. Charagu P, Peterson R. 1998. Estimates of G × E effects for economic efficiency among daughters of Canadian and New Zealand sires in Canadian and New Zealand dairy herdsProceedings INTERBULL Open Management, 18–19 January; Rotorua, New Zealand. Vol. Vol. 17. INTERBULL. Uppsala, Sweden: p. 100–104

6. Clark DA, Phyn CVC, Tong MJ, Collis SJ, Dalley DE. 2006. A systems comparison of once- versus twice-daily milking of pastured dairy cows. Journal of Dairy Science. Vol. 89:1854–1862. [Cross Ref]

7. Davis SR, Farr VC, Stelwagen K. 1999. Regulation of yield loss and milk composition during once-daily milking: a review. Livestock Production Science. Vol. 59:77–94. [Cross Ref]

8. Dickerson GE. 1977 Genetic-environmental interaction-Background information. United States Department of Agriculture. June. 1977. Texas:

9. Dillon PG, Berry DP, Evans RD, Buckley F, Horan B. 2006. Consequences of genetic selection for increased milk production in European seasonal pasture based systems of milk production. Livestock Production Science. Vol. 99:141–158. [Cross Ref]

10. Edwards JP. 2019. Comparison of milk production and herd characteristics in New Zealand herds milked once or twice a day. Animal Production Science. Vol. 59:570–580. [Cross Ref]

11. Edwards JP. 2020. A comparison of profitability between farms that milk once or twice a day. Animal Production Science. Vol. 60:102–106. [Cross Ref]

12. Edwards JP, McMillan N, Bryant RH, Kuhn-Sherlock N. 2022. Reducing milking frequency from twice each day to three times each two days affected protein but not fat yield in a pasture-based dairy system. Journal of Dairy Science. Vol. 105:4206–4217. [Cross Ref]

13. Falconer DS. 1989. Introduction to Quantitative Genetics. Longman. New York, NY, USA:

14. Gilmour AR, Gogel BJ, Cullis BR, Thompson R. 2009. ASReml User Guide Release 3.0. VSN International Ltd. Hemel Hempstead, UK:

15. Gleeson DE, O’Brien B, Boyle L, Earley B. 2007. Effect of milking frequency and nutritional level on aspects of the health and welfare of dairy cows. Animal. Vol. 1:125–132. [Cross Ref]

16. Hammami H, Rekik B, Gengler N. 2009. Genotype by environment interaction in dairy cattle. Biotechnology Agronomy Society and Environment. Vol. 13:155–164

17. Kay JK, Phyn CVC, Rius AG, Morgan SR, Grala TM, Roche JR. 2012. Once-daily milking during a feed deficit improves energy status in early lactating dairy cows. Proceedings of the New Zealand Society of Animal Production. Vol. 72:129–133

18. Kennedy E, Murphy JP, Delaby L, O’Donovan M. 2021. Effects of short-term once-a-day milking in early lactation on dairy cow performance when managed in a seasonal-calving pasture based system. Livestock Science. Vol. 251:104648. [Cross Ref]

19. Lacy-Hulbert SJ, Dalley DE, Clark DA. 2005. The effects of once a day milking on mastitis and somatic cell count. Proceedings of the New Zealand Society of Animal Production. Vol. 65:137–142

20. Lembeye F, López-Villalobos N, Burke JL, Davis SR. 2016a. Estimation of genetic parameters for milk traits in cows milked once-or twice-daily in New Zealand. Livestock Science. Vol. 185:142–147. [Cross Ref]

21. Lembeye F, López-Villalobos N, Burke JL, Davis SR. 2016b. Breed and heterosis effects for milk yield traits at different production levels, lactation number and milking frequencies. New Zealand Journal of Agricultural Research. Vol. 59:156–164. [Cross Ref]

22. Lembeye F, López-Villalobos N, Burke JL, Davis SR. 2016c. Milk production of Holstein-Friesian, Jersey and crossbred cows milked once-a-day or twice-a-day in New Zealand. New Zealand Journal of Agricultural Research. Vol. 59:50–64. [Cross Ref]

23. Lembeye F, López-Villalobos N, Burke JL, Davis SR. 2018. Genotype-by-milking frequency interaction for milk production traits in New Zealand. New Zealand Journal of Agricultural Research. Vol. 61:92–101. [Cross Ref]

24. Lynch M, Walsh B. 1997. Genotype × environment interactionGenetics and Analysis of Quantitative Traits. Lynch M, Walsh B. Sinauer. Sunderland: p. 657–683

25. McNamara S, Murphy JJ, O’Mara FP, Rath M, Mee JF. 2008. Effect of milking frequency in early lactation on energy metabolism, milk production and reproductive performance of dairy cows. Livestock Science. Vol. 117:70–78. [Cross Ref]

26. Meuwissen THE, De Jong G, Engel B. 1996. Joint estimation of breeding values and heterogeneous variances of large data files. Journal of Dairy Science. Vol. 79:310–316. [Cross Ref]

27. Montaldo HH. 2001. Genotype by environment interactions in livestock breeding programs: A review. Interciencia. Vol. 26:229–235

28. Mulder HA, Veerkamp RF, Ducro BJ, van Arendonk JAM, Bijma P. 2006. Optimization of dairy cattle breeding programs for different environments with genotype by environment interactions. Journal of Dairy Science. Vol. 89:1740–1752

29. Namkoong G. 1985. The influence of composite traits on genotype by environment relations. Theoretical and Applied Genetics. Vol. 70:315–317. [Cross Ref]

30. Nilforooshan MA, Jorjani H. 2022. Invited review: a quarter of a century—International genetic evaluation of dairy sires using MACE methodology. Journal of Dairy Science. Vol. 105:3–21. [Cross Ref]

31. Patton J, Kenny DA, Mee JF, O’Mara FP, Wathes DC, Cook M, Murphy JJ. 2006. Effect of milking frequency and diet on milk production, energy balance, and reproduction in dairy cows. Journal of Dairy Science. Vol. 89:1478–1487. [Cross Ref]

32. Phyn CVC, Kay JK, Rius AG, Morgan SR, Roach CG, Grala TM, Roche JR. 2014. Temporary alterations to postpartum milking frequency affect whole-lactation milk production and the energy status of pasture-grazed dairy cows. Journal of Dairy Science. Vol. 97:6850–6868. [Cross Ref]

33. Pomiès D, Martin B, Chilliard Y, Pradel P, Rémond B. 2007. Once-a-day milking of Holstein and Montbéliarde cows for 7 weeks in mid-lactation. Animal. Vol. 1:1497–1505. [Cross Ref]

34. Rémond B, Pomiés D, Dupont D, Chilliard Y. 2004. Once-a-day milking of multiparous Holstein cows throughout the entire lactation: milk yield and composition, and nutritional status. Animal Research. Vol. 53:201–212. [Cross Ref]

35. Rémond B, Pomiès D. 2005. Once-daily milking of dairy cows: a review of recent French experiments. Animal Research. Vol. 54:427–442. [Cross Ref]

36. Rendel J, Robertson A. 1950. Estimation of genetic gain in milk yield by selection in a closed herd of dairy cattle. Journal of Genetics. Vol. 50:1–8. [Cross Ref]

37. Ring SC, Evans RD, Cromie AR, Berry DP. 2020. Cross-sectional analyses of a national database to determine if superior genetic merit translates to superior dairy cow performance. Journal of Dairy Science. Vol. 104:8076–8093

38. Robertson A. 1959. The Sampling Variance of the Genetic Correlation Coefficient. Biometrics. Vol. 15:469–485

39. Schlamberger G, Wiedemann S, Viturro E, Meyer HHD, Kaske M. 2010. Effects of continuous milking during the dry period or once daily milking in the first 4 weeks of lactation on metabolism and productivity of dairy cows. Journal of Dairy Science. Vol. 93:2471–2485. [Cross Ref]

40. Stelwagen K, Phyn CVC, Davis SR, Guinard-Flament J, Pomiès D, Roche JR, Kay JR. 2013. Invited review: reduced milking frequency: milk production and management implications. Journal of Dairy Science. Vol. 96:3401–3413. [Cross Ref]

41. Tucker CB, Dalley DE, Burke J-LK, Clark DA. 2007. Milking cows once daily influences behavior and udder firmness at peak and mid lactation. Journal of Dairy Science. Vol. 90:1692–1703. [Cross Ref]

42. Wientjes YCJ, Calus MPL. 2017. Board invited review: the purebred-crossbred correlation in pigs: a review of theory, estimates, and implications. Journal of Animal Science. Vol. 95:3467–3478. [Cross Ref]

43. Zwald NR, Weigel KA, Fikse WF, Rekaya R. 2003. Identification of factors that cause genotype by environment interaction between herds of Holstein cattle in seventeen countries. Journal of Dairy Science. Vol. 86:1009–1018. [Cross Ref]