Interactions among breed, farm intensiveness and cow productivity on predicted enteric methane emissions at the population level

Abstract

The milk Fourier-transform infra-red spectrometry (FTIRS) can be used to predict the enteric methane emissions (EME) at population level. In this study, the variability in FTIRS predicted EME traits due to the breed of cow, farm, the production level of individual cows within herds, and their interactions were evaluated. A dataset obtained from milk recordings, which covered four breeds (Holstein, Brown Swiss, Simmental and Alpine Grey), 6,430 herds, 115,819 cows, and 1,759,706 test-day milk/spectra records was used. The herds were stratified into 5 production levels considering their average daily milk energy production; individual cows within herd were similarly stratified considering their individual production levels. The EME traits were predicted directly from milk FTIR spectra or indirectly from six informative milk FA predicted from milk spectra. The statistical model included, separately for each trait and method, breed, herd intensiveness level, cow production level, and their interactions, year, month, parity, and lactation stage. The direct and indirect methods yielded similar results in predicting CH₄ yield per kg of DMI and CH₄ intensity per kg of fat-protein corrected milk. The indirect method was reliable in predicting daily EME production per cow, the indirect one did not. EME was affected by the breed x herd production intensiveness interaction, and to a lesser degree, by the breed x cow production level interaction. A better understanding of the complex interactions influencing EME in dairy herds was achieved. This would be useful for the genetic improvement, the environmental certification of farms, and for setting prices in milk payment schemes.

HIGHLIGHTS

• Causes of variation of the FTIR predicted EME were studied at population level using data from 6,430 herds and 1,759,706 milk spectra records

• Reliable daily EME estimates were achieved from 6 six informative fatty acids predicted from milk spectra

• EME was influenced by the interaction breed × level of herd intensiveness and, to a lesser extent, by the individual production level

Keywords:

• Carbon footprint
• enteric methane emission
• global warming
• cow milk
• infra-red spectra

Introduction

International treaties aimed at reducing the rate of increase in global warming over the coming decades have been signed by the majority of the world’s countries. The major interventions concern reducing dependency on fossil fuels, increasing the soil/forest carbon sink, and limiting CH₄ emissions from domestic ruminants and landfill (Moumen et al. 2016). Enteric CH₄ emissions (EME) are the main source of greenhouse gases from agriculture. The livestock sectors of several countries are currently reducing their yearly EME, and are actively contributing to contrast global climate warming by increasing production efficiency through genetic improvement, better nutrition and management practices, and improved animal health and welfare.

The gold standard method for EME quantification remains the calculation of material and energy balances using respiration chambers (Hristov et al. 2018), but it is not feasible at the population or large nucleus levels (de Haas et al. 2017). Proxies of daily CH₄ emissions are based on analyses of exhaled gases sampled around the nostrils of animals as they access automatic feeding stations or automatic milking parlours (Difford et al. 2020). Other proxies of CH₄ emissions are based on the milk fatty acid (FA) profiles, considering their relationships with rumen fermentation (van Engelen et al. 2015; van Gastelen and Dijkstra 2016; Bougouin et al. 2019). Both methodologies were used in large surveys carried out on commercial farms, although they are expensive and labour intensive.

Predictive equations obtained from selected informative milk FA analysed using gas chromatography (van Lingen et al. 2014) and used on commercial dairy farms (1,158 cows from 85 farms) yielded results that matched expectations with regards to the effects of dairy system, parity and lactation stage (Bittante et al. 2018), and the selected FA used and the EME traits obtained were heritable (Bittante and Cecchinato 2020). Among the equations published the two proposed by van Lingen et al. (2014) for predicting CH₄ yield (g CH₄/DMI) and CH₄ intensity (g CH₄/kg CM, fat and protein corrected milk), were obtained from a meta-analysis that encompassed 30 different diets (eight experiments in two countries) comprising not only different feedstuffs and different forage:concentrate ratios, but also certain fat supplements sometimes used in commercial dairy farms.

To reduce costs and labor, and extend the possibility of EME prediction to population level, equations based on Fourier-transform infra-red (FTIR) spectra of milk were proposed (Negussie et al. 2017). The results were variable (Lassen et al. 2012; van Gastelen and Dijkstra 2016; Vanlierde et al. 2016) and raised debates over their scientific and practical value (Shetty et al. 2017; Wang and Bovenhuis 2019).

FTIRS predictions of EME were used in large surveys to quantify the effects of individual herds and of the parity and lactation stage of individual cows, and were found to have a genetic basis. FTIR spectrometry produced satisfactory results to predict the informative FA needed to use van Lingen et al.'s (2014) equations for EME prediction (indirect approach: milk spectra → milk FA → EME) (Bittante and Cipolat-Gotet 2018). A direct approach, where EME was predicted directly from milk spectra was also proposed, and both direct and indirect predictions were found to have a genetic background (Bittante et al. 2020).

However, FTIR predictions of EME have not yet been studied at the population level with different breeds reared in herds under different dairy farming systems, and with animals stratified for different potential milk yields. A scientific question concerns the effect of daily milk yield on EME traits. It is important to understand whether differences in average milk yield among different herds and differences in milk yield among different cows within herd affect EME traits; whether these two effects are of the same magnitude, direction and pattern, and whether the differences are similar across different breeds or if there is some interaction. Aim of this research was to study FTIRS predictions of informative FA and EME traits at the population level based on a large dataset of milk infra-red spectra obtained during routine milk recording operations and to compare different dairy and dual-purpose breeds across different herd intensiveness and cow productivity levels.

Material and methods

The experimental design

The study was based on a dataset obtained from the milk recording systems of Italian dairy populations and did not include any direct handling of animals to take measurements or samples, and did not require authorisation by the University Ethical Committee. It drew on the entire historical database of milk recording and artificial insemination outcomes of the Breeders Federation of Alto Adige/Südtirol (Bolzano, Italy) in the northeastern Italy. This database was used in previous studies on the relationships between milk FTIR spectra and the fertility of dairy cows (Toledo-Alvarado et al. 2017, 2018a, 2018b, 2021), and on the prediction of the cheese-making ability and urea content of milk at the population level (Bittante et al. 2021; Bittante 2022).

Milk recording data and editing

We extracted all the data from cows of the four most common breeds in the province: two specialised dairy breeds, Holstein and Brown Swiss; and two dual-purpose breeds, Simmental and Alpine Grey. Multi-breed herds, i.e. with cows of two or more breeds, about one third of total, were treated as two or more single-breed herds.

A total of 1,898,994 test-day records were extracted. Only herds with more than five cows, and cows with more than five records were retained. After editing, 6,430 herds, 115,819 cows, 291,129 lactations and 1,759,706 test days/milk spectra remained for data analyses. Parity orders were grouped into first (95,049), second (74,976), third (55,780), fourth (33,578), and fifth and subsequent (31,746).

Test day milk recording data included daily milk yield. A milk sample with preservative added was transported to the laboratory for composition analysis.

FTIR spectra of milk samples

The milk samples collected were analysed with a MilkoScanTM FT+ 6000 (Foss A/S) on the basis of their FTIR spectra. The instrument and its operations are validated and certified according to ISO 9622:2013/IDF 141:2013 (ISO 2013) recently updated to ISO 21543:2020/IDF 201:2020 (ISO 2013, 2020; Niermöller and Holroyd 2019). The complete spectrum of every milk sample analysed was stored in the experimental database. A total of 1,060 absorbance values were recorded for each milk sample covering the infra-red wavenumbers ranging from 5,000/cm (corresponding to a wavelength of 2.0 μm) in the near-infra-red subdivision of the infra-red area, through the mid-infra-red to wave number 930/cm (corresponding to a wavelength of 10.8 μm) in the far infra-red subdivision.

EME traits and informative FA

For indirect prediction of the EME traits we used the equations developed by van Lingen et al. (2014).

Methane yield per unit of dry matter intake: $${\text{CH}}_4/\text{DMI}\mathrm{ }\left(\mathrm{g}/\text{kg}\right)\mathrm{ }=\mathrm{ }23.39\mathrm{ }+\mathrm{ }9.74\times 16:0iso–\mathrm{ }1.06\times 18:1t10+t11\mathrm{ }–\mathrm{ }1.75\times 18:2c9,c12$$ where 16:0iso is iso-palmitic acid, 18:t10,+t11 is the sum of the isooleic and vaccenic acids, and 18:2c9,c12 is the α-linoleic acid of milk, all expressed as % of the sum of all milk FAs.

Methane intensity per unit of fat-protein corrected milk produced: $${\text{CH}}_4/\text{CM}\mathrm{ }\left(\mathrm{g}/\text{kg}\right)\mathrm{ }=\mathrm{ }21.13\mathrm{ }–\mathrm{ }1.38\times 4:0\mathrm{ }+\mathrm{ }8.53\times 16:0iso–\mathrm{ }0.22\times 18:1c9\mathrm{ }–\mathrm{ }0.59\times 18:1t10+t11$$ where 4:0 is butyric acid, and 18:1c9 is oleic acid, all expressed as % of the sum of all milk FAs.

In addition, daily CH₄ production per cow was obtained from CH₄ intensity.

Daily CH₄ production per cow: $${\text{dCH}}_4\left(\mathrm{g}/\mathrm{d}\right)\mathrm{ }={\mathrm{ }\text{CH}}_4/\text{CM}\mathrm{ }\times \mathrm{ }\text{dMY}$$ where dMY is the daily production of milk per cow corrected for fat and protein.

FTIR spectrometry of informative FAs and EME traits

Two datasets of EME traits predicted using the equations developed by Bittante and Cipolat-Gotet (2018) were compiled:

• Methane yield, CH₄ intensity, and daily CH₄ production were indirectly predicted (Indirect-IR) from the FTIR-predicted informative FA and the daily milk yields measured during milk recording using the aforementioned equations;

• Methane yield, CH₄intensity, and daily CH₄ production were directly predicted (Direct-IR) from FTIR spectra using prediction equations of previous papers (Bittante and Cipolat-Gotet, 2018).

Prior to data analysis, each single wavenumber of the spectra was centred and standardised to a null mean and a unit sample variance. Mahalanobis distances were calculated using the Mahalanobis function in the R software (R Core Team 2015; R Foundation for Statistical Computing, Vienna, Austria). The spectra were centred to a null mean and standardised to a unit sample variance. The spectra were not subjected to any other mathematical pre-treatment.

Separate models were fitted for all the EME traits and for the informative FA. A Bayesian model (the BayesB model) was implemented in the BGLR package of the R software (Pérez and De Los Campos 2014). Details of the chemometric procedure used are described in Ferragina et al. (2015), while details of the development of the FTIR equations for predicting informative FA and EME traits are described in Bittante and Cipolat-Gotet (2018). The accuracy of the prediction equations was assessed by a training-testing procedure.

Accuracy of the FTIR prediction equations

A random training data set (80% of the total records) was used to build the equation, and a testing data set (20% of the total) was used as validation; the procedure was repeated 10 times for each trait. The coefficients of correlation between the measured and predicted traits (r) and the root mean square error of cross-validation (RMSE_CV) of the predicted traits (Bittante and Cipolat-Gotet 2018) were:

Informative FA for indirectly predicting EME traits:

• 4:0 butyric acid, r = 0.69; RMSE_CV = 0.66%;

• 16:0iso iso-palmitic acid, r = 0.55; RMSE_CV = 0.08%;

• 18:1t10 isooleic acid, r = 0.65; RMSE_CV = 0.07%;

• 18:1t11 vaccenic acid, r = 0.66; RMSE_CV = 0.28%;

• 18:1c9 oleic acid, r = 0.82; RMSE_CV = 1.85%;

• 18:2c9,c12 linoleic acid, r = 0.70; RMSE_CV = 0.44%;

Directly predicted EME traits:

• CH₄/DMI, CH₄ yield (g/kg DM), r = 0.70; RMSE_CV = 1.18 g/kg;

• CH₄/CM, CH₄ intensity (g/kg CM), r = 0.75; RMSE_CV = 1.17 g/kg;

• dCH₄, daily CH₄ production (g/d/cow), r = 0.60; RMSE_CV = 86.0 g/d.

Stratification of herd intensiveness and cow productivity levels

Within breed, the herd average amount of milk energy produced daily by lactating cows was used as an indicator trait of herd ‘intensiveness level’. This is because it represents the cow’s net energy requirement for milk production as well as the animal’s major – and more variable – metabolic load. Unlike daily milk yield, energy output does not include the weight of water and ashes, but does include fat, protein and lactose in terms of their energy content. The herds were stratified into 5 categories of herd intensiveness level (HL) determined according to the average daily milk energy output (dMEO, MJ/d) yielded by all the lactating cows in the herd, as described in detail in a previous study using the same dataset (Bittante et al. 2021). Briefly, the net energy content (NE_L) of milk was estimated by applying the following equation, recommended by NRC (2001): $${\text{NE}}_{\mathrm{L}}\left(\text{Mcal}/\text{kg}\right)\mathrm{ }=\mathrm{ }0.0929\mathrm{ }\times \mathrm{ }\text{fat},\%\mathrm{ }+\mathrm{ }0.0547\mathrm{ }\times \mathrm{ }\text{protein},\%\mathrm{ }+\mathrm{ }0.0395\mathrm{ }\times \mathrm{ }\text{lactose},\%$$ where NE_L is the energy content of one kg milk.

The least squares means (LSM) of the dMEO of each herd were calculated separately for each breed after correcting the data for parity, class of days in milk, and year-season. The herd solutions were used to classify them into five dMEO levels: HL-1, herds < −1.5σ; HL-2, −1.5σ to −0.5σ; HL-3, −0.5σ to +0.5σ; HL-4, +0.5σ to +1.5σ; HL-5, > +1.5σ; where σ is the standard deviation of the herds’ LSM of dMEO within breed.

Similarly, within each breed and herd, the average milk energy produced daily by individual lactating cows was used as an indicator trait of the animal’s ‘productivity level’, i.e. the metabolic effort it requires for milk production. The model included the effects of parity, class of days in milk, and year-season. The milk production levels of the cows (CL) were stratified into five classes (CL-A to CL-E) according to the solutions of individual cows within HL.

The number of herds, cows and spectra available by breed, HL and CL are shown in Figure 1.

Figure 1. Number of herds, cows and spectra collected by breed, herd intensiveness level class (HL) and cow production level class (CL). Five classes of daily milk energy output (dMEO, MJ/d) for each factor: HL-1 and CL-A <-1.5σ; HL-2 and CL-B -1.5σ to -0.5σ; HL-3 and CL-C -0.5σ to +0.5σ; HL-4 and CL-D + 0.5σ to +1.5σ; HL-5 and CL-E > +1.5σ.

Statistical analysis

The informative FA and EME traits were analysed using the following linear model: , $$y_{\mathit{ijklmnop}}=\mu +P_i+D_j+{Y_k+\mathrm{ }M}_l+{HL}_m+{CL}_n+B_o+{HL}_m\times B_o+{CL}_n\times B_o+e_{\mathit{ijklmnop}}$$, where $y_{\mathit{ijklmnop}}$ is the response on the trait (FA and EME traits); µ is the general mean; $P_i$ is the parity order (i = 1st, 2nd, 3rd, 4th or ≥ 5th); $D_j$ is the category of DIM (j = 12 groups of 30 days each, with the last category open, max DIM = 450 d); $Y_k$ is the effect of Year (k = years 2010 to 2017); $M_l$ is the effect of month (l = January to December); HL_l is the Herd Level class (l = 1 to 5); CL_n is the Cow Level class (n = 1 to 5); $B_o$ is the effect of breed (o = Holstein, Brown Swiss, Simmental, and Alpine Grey); ${HL}_m\times B_o$ is the effect of the interaction between Herd Level l and breed o; ${CL}_n\times B_o$ is the effect of the interaction between Cow Level n and breed o; and $e_{\mathit{ijklmnop}}$ is the random residual $\sim \mathrm{ }\mathit{NIID}\mathrm{ }(0,{\sigma }_e^2).$ The models were fitted using the lm and aov functions in R. Details of the statistical models are given in a previous study on the same dataset (Bittante et al. 2021).

Results and discussion

Main sources of variation in milk yield and informative milk FA used for indirect EME predictions

The results of the analysis of variance carried out on the daily milk yield and the six FTIR-predicted informative FA included in the equations used for indirect prediction of EME traits are given in Table 1. Due to the very large amount of records, all the effects included in the statistical model were highly significant for every trait (p < 0.001). The relative importance of the various factors for each trait can be assessed from the size of the F-values. In accordance with the aims of this study, only the effects of breed, HL, CL, and their interactions are reported and discussed.

Table 1. Analysis of variance of milk yield and of the predicted informative fatty acids content: F-values of the factors included in the model of analysis (all the effects are highly significant for p < 0.001).

Trait	DF	Milk yield, kg/d	Predicted informative fatty acids, %:
			4:0	16:0iso	18:1t10	18:1t11	18:1c9	18:2c9c12
Mean	–	23.20	3.40	0.32	0.29	1.21	18.13	2.02
SD	–	7.40	0.59	0.06	0.11	0.43	2.54	0.44
Year of calving	7	19,078	3,952	2,799	790	809	2,709	3,104
Month of the year	11	1,061	308	593	659	355	784	560
Parity	4	41,237	687	6,558	3,568	5,891	6,497	2,438
Stage of lactation	11	130,812	29,666	17,690	895	3,982	26,647	8,661
Breed	3	372,766	966	433	240	1,147	176	159
Herd level	4	246,646	554	4,110	1,135	5,283	3,967	8,064
Cow level	4	131,360	1,222	933	300	3,708	2,723	357
Herd level by breed	12	214	153	268	262	519	60	427
Cow level by breed	12	170	46	27	51	54	27	28
RMSE	1,759,637	4.0	0.54	0.05	0.10	0.41	2.23	0.41

DF: degrees of freedom; SD: standard deviation.

The average daily milk yield (23.2 ± 7.4 kg/d) is lower than, but with a similar variability to the values (24.5 ± 7.9 kg/d) reported in previous research by Bittante and Cipolat-Gotet (2018) on direct and indirect FTIR predictions of EME traits, and in Lingen et al.'s (2014) research on fatty acid prediction of EME traits (28.9 ± 6.4 kg/d). The major source of variation in daily milk yield was breed of cow and HL, followed by CL and stage of lactation.

While oleic and linoleic acid are almost always included in the FA profiles of bovine milk, the other informative FA are found only when detailed FA profiles are analysed. The average values and SDs of the six FTIR-predicted informative FA considered here are similar to those found in previous research on detailed FA profiles of bovine milk analysed by gas chromatography (Heck et al. 2009; Pegolo et al. 2016). It should be noted that it is usually only the major milk FA that are predicted by FTIR, as those present in milk in small amounts generally have exhibit a modest level of prediction accuracy (Eskildsen et al. 2014; Ferrand-Calmels et al. 2014; Rovere et al. 2021).

The most important sources of variation in informative FA and in daily milk yield had different rankings. In neither case was breed of cow among the most important factors, despite being always significant. This is consistent with the results of previous analyses of milk FA, which revealed minor differences in the profiles across different breeds, with the exception of Jersey breed, not included here (Maurice-Van Eijndhoven et al. 2011; Vanbergue et al. 2017; Poulsen et al. 2012). Similar modest differences were found using FTIR predictions (Soyeurt et al. 2006; Maurice-Van Eijndhoven et al. 2013; Gottardo et al. 2017).

Stage of lactation was the most important source of variation for almost all the FA considered, except for the two 18:1t isomers (t10 and t11), which were more affected by parity (important also for oleic and iso-palmitic acids). The two 18:1t isomers and linoleic acid were also greatly affected by HL. The variation due to CL was less important than that due to HL (except in the case of butyric acid). The interactions between breed and both HL and CL, although less important than the main factors, were always highly significant. The LSMs of the major effects of breed of cow, HL and CL are therefore not shown in the figures, whereas the LSMs of the interactions are shown and discussed in the relevant chapters.

Main sources of variation in EME traits

The results of the analysis of variance carried out on the EME traits predicted directly by FTIR spectra (Direct-IR) and indirectly through prediction of the informative FA(Indirect-IR) are given in Table 2. Also, for EME traits, due to the large amount of data available, all the effects included in the statistical model were highly significant (p < 0.001) for every trait. In accordance with the aims of this study, only the effects of breed of cow, HL, CL, and their interactions are shown and discussed.

Table 2. Analysis of variance of the enteric CH₄ yield, intensity and production predicted directly from milk FTIR spectra or through prediction of informative fatty acids: F-values of the factors included in the model of analysis (all the effects are highly significant for p < 0.001).

Trait	DF	CH4 yield, g/kg DMI:		CH4 intensity, g/kg milk:		Daily CH4 production, g/d:
		Direct-IR	Indirect-IR	Direct-IR	Indirect-IR	Direct-IR	Indirect-IR
Mean	–	21.3	21.3	14.2	14.2	357	333
SD	–	1.3	1.1	1.3	1.2	77	100
Year of calving	7	5,030	2,957	6,830	6,266	6,249	10,575
Month of the year	11	145	196	91	360	1,193	726
Parity	4	1,917	2,148	1,944	3,515	9,051	43,591
Stage of lactation	11	19,233	18,016	46,397	65,960	17,553	39,948
Breed	3	1,108	918	1,168	1,442	686	343,342
Herd level	4	3,537	4,506	196	114	6,378	278,162
Cow level	4	368	148	25	21	4,350	139,400
Herd level by breed	12	614	158	64	97	677	405
Cow level by breed	12	34	18	11	31	57	323
RMSE	1,759,637	1.2	1	1	1	70	58

DF: degrees of freedom; SD: standard deviation.

The means and SDs of Direct-IR and Indirect-IR CH₄ yield were very similar, as were those of Direct-IR and Indirect-IR CH₄ intensity. Note that all these predictions were based on FTIR spectra information, the difference being that in the case of Indirect-IR the spectra were used to predict the informative FA that, in turn, were used to estimate the EME traits. The mean values of these two traits (CH₄ yield and intensity) in both versions (Direct-IR and Indirect-IR) were very similar but the SDs were lower than those of the population used for calibrating the FTIR equations (21.3 ± 1.6 g/kg DMI, and 14.2 ± 1.8 g/kg CM; Bittante and Cipolat-Gotet 2018), and those used for defining the equations based on FA (van Lingen et al. 2014).

Methane yield and intensity were both particularly affected by stage of lactation, but also by year of calving and HL.

The results for daily CH₄ production differ from the previously cited EME traits because the Direct-IR and Indirect-IR predictions yielded very different results and were affected by different sources of variability. It is worth noting that while the Direct-IR daily CH₄ production was predicted directly and solely from FTIR spectra, the reference values used for calibrating the prediction equation were obtained by multiplying the daily milk yield by the CH₄ intensity per kg of fat-protein corrected milk of each cow (Bittante and Cipolat-Gotet 2018). Chemometric methods are required to extrapolate from the milk FTIR spectra (which constitute substantially qualitative information) the amount of milk produced daily by the cows and the amount of CH₄ released into the atmosphere per kg of milk. Using FTIR spectrometry to predict quantitative information is expected to be difficult, which explains the fact that Direct-IR CH₄ production was the EME trait prediction with the lowest correlation with the reference data (r = 0.60).

Indirect-IR daily CH₄ production was, in contrast, obtained from the cows’ daily milk yield measured during milk recording (quantitative information), whereas FTIR spectra alone were used to predict CH₄ intensity per kg of corrected milk (r = 0.75) (qualitative information).

This explains the very different means and SDs obtained for Direct-IR and Indirect-IR daily CH₄ production and for the ANOVA. That the ranking of the major sources of variation for Indirect-IR CH₄ production is much more similar to daily milk yield (Table 1) than to the other EME traits (Table 2) was expected. Breed was the most important factor affecting Indirect-IR CH₄ production (and daily milk yield), followed by HL then CL.

Effects of herd intensiveness level and cow productivity levels by breed of cow on informative milk FA

In our database, breed of cow was the most important source of variability in milk yield, which ranged on average from 26.8 kg/d in the Holstein cows to 15.6 kg/d in the Alpine Grey cows, with Brown Swiss and Simmental intermediate (Figure 2). Note that this study was carried out in an Alpine region, where there is a variety of dairy farming systems, ranging from the very traditional to the modern (Schiavon et al. 2019; Berton et al. 2020). The differences in the average milk yields of the various breeds therefore reflect not only different genetic backgrounds but also different farming systems. Previous research on the same area which compared six breeds of cow reared in mixed herds (Stocco et al. 2018) found a difference within herd between Holstein and Alpine Grey cows of 8.0 kg/d, as opposed to the 11.2 kg/d found here across herds. However, the variability among farms was also very large within breed, which is why we included HL in the statistical model. Note that this is a within-breed classification factor and that it is not corrected for the different environmental conditions in which the different breeds are raised. To make the significance of the HL within breed clearer, we plotted all the LSMs of the investigated traits (y axis) against the actual daily milk yield of each breed and the Herd or CLs (Figure 3).

Figure 2. Average daily milk production by breed of cows (kg/d).

Figure 3. Milk fat, butyric acid (4:0) and iso-palmitic acid (16:0 iso), in % of total fatty acids (LSM ± CI), by breed, herd intensiveness class level or cow production class level plotted against actual cow yield (kg/d).

The first plot in Figure 3 reports the 20 LSMs of the milk fat content (in %) of each of the four breeds and at each of the five HLs. Here, milk fat serves to compare the pattern of total lipid content of the milk with the proportion of the informative FA in the same fat. Large differences between the four breeds in terms of milk fat content were evident (compare the medians, i.e. the third value of each curve). All breeds exhibited similar but not identical patterns of increasing HL. All patterns increased curvilinearly with small differences between the first three HLs and increasingly large differences moving to the fourth and especially the fifth HLs. The increase from first to fifth HL was greater for the Holstein breed than for the Brown Swiss, with the two dual-purpose breeds intermediate (the interaction between breed and HL was highly significant, data not shown). No overlap in the relationships between milk fat content and daily milk yield in the different breeds was found. However, the Alpine Greys and Holsteins had a common trajectory, which could be interpreted as the difference in fat content of these two breeds being dependent mainly on their different milk yields. Simmental, and especially Brown Swiss cows had different parallel trajectories.

The effect of CL within breed had similarities with the HL effect, but also some differences (second plot of Figure 3). The rankings and the median values of the four breeds were obviously the same, as the two plots are based on the same dataset. The increase in milk fat with increasing CL was, in this case, almost linear and the difference between the first and last classes of each breed was smaller than for HL. The fact that both HL and CL were favourably associated with the increase in milk fat content was because these classes were based not simply on milk yield, but on the daily output of net energy. So, the classes represent an estimation of the ability of the cows (herds) to produce valuable nutrients, and not simply water excretion, and also quantify the cows’ metabolic burden for lactation.

The informative FA exhibited patterns that were very different from each other and also from total milk fat content. The butyric acid content of milk fat (Figure 3, 3rd and 4th plots) increased curvilinearly for all breeds at the CL and for Holsteins and Brown Swiss at the HL, but in an opposite pattern to that of milk fat, increasing rapidly from the lower to the central classes then much less so to the higher classes. The curves in part overlap because the average butyric acid content differed significantly in the four breeds, although the differences were not relevant. Note that butyric acid was included in the equation for estimating CH₄ intensity per kg milk with a negative sign (van Lingen et al. 2014), which means that the increase in this fatty acid with increasing Herd and CLs results in a favourable effect of herd intensiveness and cow production, reducing the EME of lactating cows per kg of milk produced.

Iso-palmitic acid (Figure 3rd, 5th and 6th plots), in contrast, decreased with increasing Herd and CLs. The median values of the breeds were similar, but the decrease was much greater at the HL than at the CL, and in Holsteins and Simmentals than in Brown Swiss and Alpine Greys. Note that iso-palmitic acid is the only one of the six informative FA with positive correlations with both CH₄ yield and intensity (van Lingen et al. 2014). This means that the decrease in the proportion of iso-palmitic acid in milk at increasing Herd and CLs along with the increase in butyric acid leads to a decrease in CH₄ emitted per unit of DMI and per unit of milk produced. A positive association between iso-palmitic acid in milk and EME was also found by Rico et al. (2016), and positive associations between palmitic acid and EME were found by several other authors (Chilliard et al. 2009; Dijkstra et al. 2011), as was the case for many saturated FA and especially those derived by mammary de novo synthesis (Chilliard et al. 2009).

Isooleic acid (18:1 trans10) increased with HL, steeply in the case of Alpine Grey herds, moderately in Brown Swiss and Simmental herds, and very modestly and curvilinearly in Holstein herds (1st and 2nd plots in Figure 4). The effect of CL on this fatty acid was smaller and almost the opposite (decreasing, especially in Holsteins and Simmentals). This milk fatty acid was also correlated with a decreased CH₄ yield and intensity (van Lingen et al. 2014), like several other isomers of oleic acid (Chilliard et al. 2009; Mohammed et al. 2011). It is worth noting that isooleic acid had a favourable effect on the carbon footprint at increasing levels of dairy farm intensiveness, but an unfavourable effect at increasing cow production levels within farm (Figure 4, 1st and 2nd plots).

Figure 4. Milk isooleic acid (18:1 trans10), vaccenic (18:1 trans11), oleic (18:1 cis9), and linoleic acid (18:2 cis9,cis12), in % of total fatty acids (LSM ± CI), by breed, herd intensiveness class level or cow production class level plotted against actual cow yield (kg/d).

Vaccenic acids exhibited very different patterns among the HLs (Figure 4, 3rd plot), with the extremes represented by the small linear increase in Alpine Grey herds and the strong curvilinear decrease in Holstein herds. The effect of CL (Figure 4, 4th plot), on the other hand, was very similar in the four breeds with relevant linear decreases in all cases. The negative association between vaccenic acid and EME traits (van Lingen et al. 2014) was also reported by other authors (Dijkstra et al. 2011; Mohammed et al. 2011; van Gastelen et al. 2018). In addition, the variation in vaccenic acid contributed to an unfavourable effect on EME at increasing HL and CLs.

The major fatty acid considered, 18:1 cis9 (oleic acid), decreased at increasing Herd and CLs (more in the former and less in the latter case) (Figure 4, 5th and 6th plots). Linoleic acid exhibited an almost opposite pattern to oleic acid, with a strong curvilinear increase for HL and much smaller curvilinear increases for CL (Figure 4, 7th and 8th plots). As both FA were found to have a negative relationship with the EME of dairy cows (van Lingen et al. 2014), they had opposite effects on CH₄ emission: oleic acid unfavourable and linoleic acid favourable at increasing HL and CL. Several authors have found oleic acid to have a tempering effect on EME, whereas the effect of linoleic acid was less consistent and was found to be negatively (Mohammed et al. 2011) and positively (van Gastelen et al. 2018) correlated with the EME of dairy cows.

The complexity of the interactions between breed, farm intensiveness level and cow productivity, and the different weights assigned to the informative FA in the equations for estimating CH₄ yield and CH₄ intensity may have the potential to capture the differences in feeding regimes, metabolic burdens and the management and environmental conditions of cows and herds. At least in part, the greater differences observed for increasing HL, compared to CL, would be likely due to different feeding practices applied to different breeds and different farms within breeds in terms of forage to concentrate ratio and feed ingredients, that could modify the rumen fermentations towards different proportion of volatile fatty acids.

Effects of herd intensiveness level and cow productivity level by breed of cow on EME traits

EME traits, whether derived directly or indirectly from informative FA and milk FTIR spectra, were also characterised by different patterns in different breeds at increasing HL and CL. The effects on EME of interactions between the cow’s genotype and farming systems were observed in previous studies (Ross et al. 2014).

The effect of HL on CH₄ yield was always curvilinear, decreasing particularly at higher HLs, and stronger in Holstein (and Simmental) than in the other breeds (Figure 5). Comparing the two methods of estimating EME, Direct-IR (EME predicted from milk FTIR spectra) and Indirect-IR (EME predicted from informative FA predicted from milk FTIR spectra), we find they had similar curved shapes, similar rankings of breeds, but the HL effect had a broader range in the former case than in the latter. The effect of CL on CH₄ yield was much less evident than the effect of HL, and the difference between the two methods of EME prediction was modest. It is evident that the variation in dMEO (daily milk net energy output; MJ/d) of the different herds reflected above all differences in the compositions of the diets and DMI.

Figure 5. Predicted enteric CH₄ yield (LSM ± CI, CH₄/DMI, g/kg dry matter intake) predicted directly from milk FTIR spectra or indirectly from predictive informative milk fatty acids by breed, herd intensiveness class level or cow production class level plotted against actual cow yield (kg/d).

At the CL, the animals belong to the same herd and were often fed the same diets (especially if total mixed rations are used). In this case, the only differences arise from the different DMI of cows of different productivity levels and from feed selection behaviour. Bear in mind that stratification of the cows was based on the average output of their entire productive life and therefore did not reflect differences in different lactation stages. In more traditional farms, an amount of concentrate was provided in the milking parlour or, in the case of tie stalls, at the manger (sometimes using automatic feeders). In any case, the same compound feed was generally provided to all animals and the average difference between individual cows over the entire lactation period was moderate.

Regarding CH₄ intensity per kg of fat-protein corrected milk, the prediction methods yielded similar patterns with respect to the interaction between the HL and CL and the breed of cow (Figure 6). The Alpine Greys had the highest average values and the greatest decrease in CH₄ intensity at increasing Herd and CLs. In contrast, Holsteins had the lowest emissions and were the least affected by HL and CL. The Simmental and Brown Swiss herds and cows were intermediate with respect to the other breeds.

Figure 6. Predicted enteric CH₄ intensity (LSM ± CI, g/kg of corrected milk) predicted directly from milk FTIR spectra or indirectly from predictive informative milk fatty acids by breed, herd intensiveness class level or cow production class level plotted against actual cow yield (kg/d).

Lastly, the results for daily CH₄ production (Figure 7) clearly show the inability of the Direct-IR method to capture the expected increase in CH₄ emissions with increasing HL and CL, and hence the need for further information in addition to the milk FTIR spectra (van Gastelen and Dijkstra 2016). In fact, the patterns were flat and there was overlap among the different breeds. In the case of the Indirect-IR method, there was a strong linear relation between milk yield and CH₄ emission. The lines of the four breeds were substantially on the same trajectory. The median values of the different breeds were very different, in accordance with their different average milk yields. It appeared that at the same actual milk yield the CH₄ emissions of Holstein herds and cows were slightly lower than those of the other breeds. This was in accordance with the lower CH₄ intensity per kg corrected milk (Figure 6), but also with the lower average energy content of Holstein milk. These results could be useful for better understanding the determinants and complex interactions resulting in CH₄ emissions in commercial dairy herds.

Figure 7. Predicted daily enteric CH₄ production (LSM ± CI, dCH₄, g/d) predicted directly from milk FTIR spectra or indirectly from predictive informative milk fatty acids by breed, herd intensiveness class level or cow production class level plotted against actual cow yield (kg/d).

How could FTIR EME predictions be used to mitigate the dairy sector’s carbon footprint?

The availability of a tool for rapidly and cheaply predicting EME traits of dairy cows at the population level has a threefold potential:

1. First, spectra of individual milk samples taken during milk recording can be used for the genetic improvement of dairy populations.

2. Second, spectra of bulk milk samples collected by dairy cooperatives from member farms for milk payment schemes can be used to determine possible premiums/penalties on the price of the milk in accordance with the environmental impact of each farm.

3. Third, almost all the milk produced for commercial purposes could be certified using predictions obtained from FTIR spectra of individual or bulk milk samples.

In all cases, an initial problem concerns the selection of the most useful EME trait to attain the desired target at the dairy chain level. At the commercial level, the objective would be to certify and apply penalties according to the amount of CH₄ [or better, the global temperature potential to consider other greenhouse gases and the short lifetime of methane in the atmosphere (Mazzetto et al. 2022] emitted for every kg of milk (or dairy product) produced. The most useful trait in these cases is CH₄ intensity.

Genetic improvement is a matter of debate (Knapp et al. 2014; Hristov et al. 2018; Beauchemin et al. 2020). The majority of publications in this area are concerned with daily CH₄ production (Negussie et al. 2018). As genetic improvement concerns the animal, the objective should be to minimise the cow’s daily EME.

The use of FTIR spectrometry to develop prediction equations for estimating daily CH₄ production from selected informative milk FA was investigated by Kandel et al. (2015) using a procedure similar to the one used in this study. They estimated daily CH₄ production from the detailed fatty acid profiles obtained from 602 milk samples analysed by gas chromatography, to which they applied four of the equations developed by Chilliard et al. (2009) based on a study (8 cows, 4 diets, 32 balances) exploring the correlations between milk fatty acid profiles obtained by gas chromatography and daily CH₄ production predicted using sulphur hexafluoride as tracer. They then regressed the predicted daily CH₄ outputs on the FTIR spectra of the same milk samples to obtain four FTIR-based prediction equations. In the next step they applied these (Direct-IR) equations to a population database of 604,028 milk spectra (70,872 cows), then estimated the genetic parameters of the predicted EME traits. They concluded that ‘the study showed the feasibility of using MIR spectrometry results to predict fatty acid-derived CH₄ indicator traits developed in the literature’, and that ‘the estimated genetic parameters of these traits suggested a potential phenotypic and genetic variability of the daily quantity of CH₄ eructed by Holstein dairy cows’. In a later research using the same database (Kandel et al. 2017), they found that the phenotypic and genetic correlations between FTIR-predicted daily CH₄ production and daily production of milk, fat and protein were generally modest and sometimes negative. This contradicts the frequent observation that higher daily milk yields go hand-in-hand with higher nutrient requirements and feed intake, hence greater CH₄ production (von Soosten et al. 2020). The positive aspect of increasing daily milk yield is that CH₄ production tends to increase less than proportionally, hence the CH₄ intensity per kg of milk is expected to decrease. This means that the selection objective of dairy cow populations should not be to reduce daily CH₄ production per animal (correlating with lower milk yields, and the rearing of more cows), but rather to reduce CH₄ intensity, mainly by increasing production efficiency.

It is also important that when planning a strategy for mitigating the environmental impact of a given sector a holistic approach should be taken, as one action aimed at reducing a source of greenhouse gas could cause an increase in other emissions or a worsening of other impact categories (Gerber et al. 2013; Lynch 2019). Moreover, even when considering EME alone, it is essential to evaluate the correlations with fertility, longevity, and other welfare and health traits of the cows, as it has been shown that variations in these traits may have a more important effect on the sector’s EME than the variation expected from the genetic improvement of production traits (Van Middelaar et al. 2015; Llonch et al. 2017; von Soosten et al. 2020; Beauchemin et al. 2020).

Conclusions

The enteric CH₄ emissions of dairy cows predicted at the population level from milk FTIR spectra and informative fatty acid proportions vary according to the predicted EME trait considered, the methodology used, the breed of cow, herd intensiveness, the productivity of individual cows within herd, and their interactions. Direct-IR (milk FTIR spectra → EME) and Indirect-IR (milk FTIR spectra → informative milk FA → EME) yield similar results when used to predict CH₄ yield per kg of DMI and CH₄ intensity per kg of fat-protein corrected milk. However, Direct-IR was found to be unreliable for predicting daily CH₄ production in g/d, which could be more accurately estimated by multiplying the recorded daily corrected milk yield by the CH₄ intensity per kg milk. Different breeds were differently affected by the level of herd intensiveness and, to a lesser degree, by the individual cow’s production level. These results offer a better understanding of the determinants and complex interactions that give rise to the CH₄ emissions of commercial dairy herds, and could be useful for the genetic improvement of dairy cows, environmental certification of farms, and perhaps also for determining new milk payment schemes.

Ethical approval

The study was based on a dataset obtained from the milk recording systems of Italian dairy populations and did not include any direct handling of animals to take measurements or samples, and did not require authorisation by the University Ethical Committee.

Acknowledgments

We thank the Associazione Provinciale delle Organizzazioni Zootecniche Altoatesine/Vereinigung der Südtiroler Tierzuchtverbände and the Federazione Latterie Alto Adige/Sennereiver- band Südtirol (Bolzano/Bozen, Italy) for providing the data.

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper

Data availability statement

The results and analyses presented in this paper freely available upon request.

Additional information

Funding

The authors thank the DAFNAE department of the University of Padua that approved the project 'Studio delle problematiche relative all’utilizzo della spettroscopia nel vicino e medio infrarosso per la valutazione della qualità dei prodotti di origine animale’ [DOR1970230, DOR2070528, DOR2129201, DOR2254442] and for the financial support.