Heritability estimates of enteric methane emissions predicted from fatty acid profiles, and their relationships with milk composition, cheese-yield and body size and condition

Abstract

In the present study we estimated the genetic parameters of enteric methane emissions (EME) traits predicted from milk fatty acid profile (FA) and those of their predictors in 1,091 Brown Swiss cows reared on 85 farms in order to assess the potential of using EME-related phenotypes in selective breeding. Univariate and bivariate genetic models were fitted in a Bayesian framework. The means of the marginal posterior distribution of intra-herd heritability ranged from 0.12 for estimated methane production (g/d/cow) to 0.24 for estimated methane yield (g/kg dry matter intake [DMI]), with intermediate values for estimated methane intensity, increasingly higher when expressed per kg of corrected milk (0.13), fresh cheese (0.16), or cheese solids (0.20). Regarding the correlations, the milk quality traits and percentage cheese yields were generally moderately correlated with the estimated EME traits, and were variable in terms of sign. Daily milk and cheese yield traits were, as expected, all highly positively correlated with estimated daily methane production. In contrast, they were negatively correlated with estimated methane yield and intensity, the estimates being large in the case of phenotypic and herd correlations, and low in the case of additive genetic and residual correlations. With the exception of the negative correlations with daily methane production, EME traits exhibited trivial correlations with body size and BCS of cows, which, in turn, were negatively correlated with milk yield. Although the results should be validated on a larger population and different breeds, our study demonstrate the presence of additive genetic variation of EME traits, which could be exploited in breeding programmes for the improvement in both milk production and the ecological footprint of dairy farming.

Highlights

• Enteric methane emissions (EME) of dairy cows can be estimated on the basis of milk fatty acid profile.

• EME exhibited exploitable genetic variation.

• Genetic selection could be preferentially based on predicted methane intensity per kg of milk, or per kg of cheese in countries where milk production is used mainly for cheese-making.

Keywords:

• Enteric methane emissions
• dairy cattle
• genetic parameters
• ecological footprint
• greenhouse gases
• global warming

Introduction

In agriculture, ruminant productions are considered to be the major activities responsible for global greenhouse gas (GHG) emissions, because of enteric methane emissions (EME) (Knapp et al. 2014), to which dairy cows make a significant contribution. Also in Italy EME represent the most important source of GHG of the dairy cows (Pirlo and Carè 2013; Battaglini et al. 2014; Pulina et al. 2017; Lovarelli et al. 2019).

Direct quantification of GHG through the gold standard method, based on respiration chambers, requires facilities, tools, capital investments and knowledge that are available only in some research centres, making it very expensive to directly evaluate EME on a large number of animals.

Analysis of milk fatty acid (FA) profiles and the use of combinations of FAs are considered valuable methods for use in the field, as recently reviewed by Negussie et al. (2017). This is because there are direct links between EME, the type of fermentation in the fore-stomach, the quantities and proportions of volatile FAs in the rumen, the quantities and proportions of FAs absorbed by the intestines, and the transfer and de novo synthesis of FAs in the udder, as shown in many experiments on cow’s feeding (Chilliard et al. 2000).

Interesting studies recently carried out on the relationships between the FA profiles of milk obtained by gas-chromatography (GC), considered to be the gold standard for this type of analysis, and some EME traits have been reviewed by van Gastelen and Dijkstra (2016). These studies were carried out with small numbers of cows in respiration chambers, so the results cannot be easily generalised to field conditions. On the other hand, van Lingen et al. (2014) undertook a meta-analysis of results from research in this area, by combining the data from 4 experiments carried out at the University of Reading (UK) and 4 experiments carried out at Wageningen University (the Netherlands), covering 30 different diets and 146 complete gas balances. They devised two equations for predicting methane yield per unit of DMI (CH₄/DMI, in g/kg), and methane intensity per unit of fat/protein-corrected milk (CH₄/CM, in g/kg), both with acceptable levels of accuracy (R²: 0.54 and 0.47, respectively) considering that the data came from different experiments with different diets and environmental conditions. Similar prediction performances where recently obtained by the same group (van Gastelen et al. 2018) using the results of Dutch experiments for predicting also the daily methane production per cow (dCH₄, in g/d).

In a previous study (Bittante et al. 2018), we applied these equations of van Lingen et al. (2014) to a large dataset of detailed FA profiles of milk samples from a survey on 1,158 Brown Swiss cows from 85 herds from different dairy systems. We also used this information, combined with the cheese yield of individual cows to estimate methane intensity per unit of fresh cheese (CH₄/CY_CURD, g/kg) and per unit of cheese solids (CH₄/CY_SOLIDS, g/kg). In addition, we estimated dCH₄, g/d by multiplying CH₄/CM by daily corrected milk yield (dCMY, kg/d). We were able to characterise the estimated EME traits of different dairy systems and the effects of individual sources of variation (i.e. DIM and parity). No studies on the individual genetic variation of all the aforementioned traits have been conducted so far.

The aims of this study were: i) to estimate the heritabilities of 5 EME traits obtained from the yield, composition and detailed FA profiles of milk, and from individual cheese yields; ii) to test, through a sensitivity analysis, the robustness of these estimates across different dairy systems and feeding, parity, and lactation stage; iii) to estimate the additive genetic, phenotypic, herd, and residual (co)variances among the 5 estimated EME traits, and between these and the phenotypes used for their estimation, and also the body size and condition.

Material and methods

Animals, samples and analyses

This study is part of the Cowability-Cowplus project (Bittante et al. 2015) and details about the animals, herds, milk sampling and analyses, detailed gas-chromatographic fatty acid profiling and individual model-cheese making were described in details in the previous study on the estimation of EME traits from FA profile and on the analysis of their phenotypic sources of variation (Bittante et al. 2018). Briefly, after editing, a total of 1,091 Brown Swiss cows from 85 herds, all enrolled in the Herd Book of the breed, were considered. Pedigree information was provided by the Italian Brown Swiss Cattle Breeders Association (ANARB, Verona, Italy). We considered cows with phenotypic records available for the investigated traits and all known ancestors. Each sampled cow had known ancestors for at least four generations and the pedigree file included 8,845 animals.

The body size and condition of each cow were defined as: measured heart girth, and body weight (BW) and body condition score (BCS), estimated by a trained university technician. The milk yield was sampled once, and two milk subsamples per cow were immediately refrigerated (4 °C) without preservative and used to analyse milk quality traits (50 mL), and to produce individual model cheeses (2,000 mL) at the Cheese-Making Laboratory of the Department of Agronomy, Food, Natural Resources, Animals and Environment (DAFNAE) of the University of Padova (Legnaro, Padova, Italy). Individual cheese yields, expressed in percentages, were measured as the ratio between the weight of the fresh cheese wheel and the weight of milk processed (%CY_CURD), and as the ratio between the weight of the fresh wheel multiplied by its percentage DM and the weight of the milk processed (%CY_SOLIDS). From the daily milk yield (dMY), individual daily productions (kg/d) of fresh cheese (dCY_CURD) and cheese solids (dCY_SOLIDS) were calculated by multiplying dMY by the traits concerned (%CY_CURD or %CY_SOLIDS), and of daily fat and protein corrected milk production (dCMY) according to van Lingen et al. (2014).

Detailed FA profiles of the milk samples were obtained using a ThermoQuest gas chromatograph (Thermo Electron Corp., Waltham, MA) fitted with a flame-ionization detector and a high polar fused-silica capillary column (Chrompack CP-Sil88 Varian, Middelburg, the Netherlands; 100 m, 0.25 mm i.d.; film thickness 0.20 μm) by the method described in detail by Mele et al. (2016). Milk FA composition was expressed as grams per 100 g of total FAs, and, in conformity with van Lingen et al. (2014) equations, we examined the following FAs: 4:0 (butyric acid); 16:0iso (iso-palmitic acid); 18:1t10; 18:1t11 (vaccenic acid); 18:1c9 (oleic acid); and 18:2c9,c12 (linoleic acid).

Definition of EME phenotypes

In the present study, we examined the following estimated EME traits, as defined by de Haas et al. (2017):

• Methane yield: (CH₄/DMI), predicted according to van Lingen et al.’s (2014) equation: $${\text{CH}}_4/\text{DMI}\left(\mathrm{g}/\text{kg}\right)=23.39+9.74\times \text{16:0}iso–1.06\times 18:1t10+t11–1.75\times 18:2c9,c12\left({\mathrm{R}}^2=0.54\right)$$

where 16:0iso is iso-palmitic acid, 18:t10,+t11 is the sum of the 18:1t10 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-milk: (CH₄/CM), predicted according to van Lingen et al.’s (2014) equation: $${\text{CH}}_4/\text{CM}\left(\mathrm{g}/\text{kg}\right)=21.13–1.38\times 4:0+8.53\times 16:0iso–0.22\times 18:1c9–0.59\times 18:1t10+t11\left({\mathrm{R}}^2=0.47\right)$$

where 4:0 is butyric acid, and 18:1c9 is oleic acid, all expressed as % of the sum of all milk FAs.

• Methane production: (dCH₄), calculated as: $${\text{dCH}}_4\left(\mathrm{g}/\mathrm{d}\right)={\text{CH}}_4/\text{CM}\times \text{dCMY}$$

• Methane intensity-cheese: (CH₄/CY_CURD), calculated as: $${\text{CH}}_4/{\text{CY}}_{\text{CURD}}\left(\mathrm{g}/\text{kg}\right)={\text{dCH}}_4/{\text{dCY}}_{\text{CURD}}$$

• Methane intensity-cheese solids: (CH₄/CY_SOLIDS), calculated as: $${\text{CH}}_4/{\text{CY}}_{\text{SOLIDS}}\left(\mathrm{g}/\text{kg}\right)={\text{dCH}}_4/{\text{dCY}}_{\text{SOLIDS}}$$

• Daily DMI of each cow (dDMI_est), indirectly estimated (only to test the consistency of van Lingen et al.’s (2014) equations) as: $${\text{dDMI}}_{\text{est}}\left(\text{kg}/\mathrm{d}\right)={\text{dCH}}_4/\left({\text{CH}}_4/\text{DMI}\right)$$

Statistical analysis

Non-genetic effects, which were included in the mixed models used to estimate the genetic parameters for the estimated EME traits and for the traits used to make the estimations, were identified through a preliminary analysis reported in detail by Bittante et al. (2018). For all traits, the model accounted for the effects of herd/date (85 levels), days in milk (DIM: class 1, <60 d; class 2, 60 to 120 d; class 3, 121 to 180 d; class 4, 181 to 240 d; class 5, 241 to 300 d; class 6, >300 d), and parity (1 to 4 or more).

The genetic determinism of the estimated EME traits and their predictors ($\mathrm{y}$) was investigated by analysing the data with the following mixed model: [1] $$\mathrm{y}=\mathbf{Xb}+{\mathrm{Z}}_1\mathbf{h}+{\mathrm{Z}}_2\mathbf{a}+\mathbf{e}$$[1] where $\mathrm{y}$ is the vector of phenotypic records with dimension n; $\mathbf{X},$ ${\mathrm{Z}}_1,$ and ${\mathrm{Z}}_2$ are the appropriate incidence matrices for systematic effects ($\mathbf{b}$), herd/date effects ($\mathbf{h}$), and polygenic additive genetic effects ($\mathbf{a}$), respectively. In $\mathbf{b}$ were included the effects of days DIM and parity. The priors for $\mathbf{b}\mathrm{ }$ and the variance components were assumed to be flat.

A sensitivity analysis was carried out repeating the genetic analyses on reduced datasets obtained excluding some cow’s category. To test the effect of dairy system and feeding, two datasets were obtained removing all traditional farms with tied cows, and, alternatively, all modern farms using total mixed rations (TMR) with or without use of maize silage. To test the effect of parity, two datasets were obtained removing all primiparous or all old (≥4 parities) cows, respectively. Finally, to test the effect of the stage of lactation, two datasets were obtained removing all cows in early (<60 DIM) or late (>240 DIM) lactation, respectively.

We estimated the genetic, herd-date and residual correlations between the studied variables by conducting a set of bivariate analyses that implemented model 1 in its multivariate version. In this case, the traits involved were assumed to jointly follow a multivariate normal distribution along with the additive genetic, herd/date and residual effects. The corresponding prior distributions for these effects were: $$\mathbf{a}|{\mathbf{G}}_0,\mathrm{ }\mathbf{A}\mathrm{ }\sim \mathrm{ }\mathit{MVN}\left(0,{\mathbf{G}}_0,⨂\mathbf{A}\right),$$ $$\mathbf{h}|{\mathbf{H}}_0,\mathrm{ }\sim \mathrm{ }N\left(0,{\mathbf{H}}_0,⨂{\mathbf{I}}_n\right)\mathrm{ }\mathrm{and}$$ $$\mathbf{e}|{\mathbf{R}}_0,\mathrm{ }\sim \mathrm{ }N\left(0,{\mathbf{R}}_0,⨂{\mathbf{I}}_m\right),$$ where ${\mathbf{G}}_0,\mathrm{ }{\mathbf{H}}_0,\mathrm{ }{\mathbf{R}}_0$ are the corresponding variance-covariance matrices between the traits involved, and $\mathbf{a},\mathrm{ }\mathbf{h},\text{ and }\mathbf{e}$ are vectors of dimension equal to the number of animals in the pedigree (n and m) times the number of traits considered.

Marginal posterior distributions of all unknowns were estimated using the Gibbs Sampling algorithm. The analyses were implemented using the software TM (available at http://cat.toulouse.inra.fr/∼alegarra/). For all analyses, the total number of iterations was 850,000 with a burn-in of 50,000 and a thin interval of 100. The posterior mean was used as the point estimate for all parameters. For heritability estimates, the lower and upper bounds of the highest 95% probability density regions (HPD95%) were obtained from the estimated marginal densities. For the phenotypic, genetic, herd and residual correlations, besides the mean of each marginal posterior distribution, we also estimated the probability of each mean being greater than 0 when the mean is positive, or lower than 0 when the mean is negative (P). We considered all estimates with P greater than 95% as ‘relevant’ correlations.

Intra-herd heritability ($h^2$) was computed as $h^2=\mathrm{ }{\sigma }_{\mathrm{a}}^2/{(\sigma }_{\mathrm{a}}^2+{\sigma }_{\mathrm{e}}^2),$ and the herd/date incidence as $HD=\mathrm{ }{\sigma }_{\mathrm{h}}^2/{(\sigma }_{\mathrm{a}}^2+{\sigma }_{\mathrm{h}}^2+{\sigma }_{\mathrm{e}}^2).$

Results

Mean and standard deviation of the estimated EME traits and of their predictors (informative milk FAs, quality traits, and daily yield traits), and also of cow’s size and condition traits are reported in Table 1. A comprehensive discussion of these traits was reported by Bittante et al. (2018).

Table 1. Descriptive statistics of all the investigated traits.

Item	Mean	SD
Estimated Methane emissions^a
CH4/DMI, g/kg	21.34	1.58
CH4/CM, g/kg	14.26	1.63
dCH4, g/d	358	104
CH4/CYCURD, g/kg	100	16
CH4/CYSOLIDS, /kg	208	31
dDMIest, kg/d	17.1	5.5
Informative milk FAs (%FA)^b
4:0 - Butyric acid	3.45	0.90
16:0iso - Iso-Palmitic acid	0.32	0.09
18:1t10	0.29	0.09
18:1t11 - Vaccenic acid	1.20	0.37
18:1c9 - Oleic acid	18.33	3.19
18:2c9,c12 - Linoleic acid	2.04	0.59
Quality traits^c
Milk fat, %	4.23	0.72
Milk protein, %	3.71	0.43
%CYCURD, %	15.06	1.89
%CYSOLIDS, %	7.23	0.94
Daily yield traits^d
dMY, kg/d	24.5	7.9
dCMY, kg/d	25.6	8.1
dCYCURD, kg/d	3.66	1.17
dCYSOLIDS, k/d	1.76	0.57
Body size and condition
Heart girth, cm	200	8
Est. body weight, kg	668	51
Body condition score	2.97	0.34

^a Estimated Methane Emissions: CH₄/DMI, g/kg: methane yield, emitted per kg DMI; CH₄/CM, g/kg: methane intensity per kg fat and protein corrected milk produced; dCH₄, g/kg: daily methane production per cow; CH₄/CY_CURD, g/kg: methane intensity per kg of fresh cheese produced; CH₄/CY_SOLIDS g/kg: methane intensity per kg of cheese solids produced; dDMI_est, kg/d: estimated daily DMI of cows.

^b Informative milk FAs are the fatty acids included as independent variables in the equations used to estimate the enteric methane emissions (van Lingen et al. 2014).

^c Quality traits: %CY_CURD: wt of fresh cheese as % of processed milk; %CY_SOLIDS: wt of cheese solids as % of processed milk.

^d Daily yield traits: dMY: daily milk yield; dCMY: daily fat and protein corrected milk yield; dCY_CURD: daily production of fresh cheese per cow; dCY_SOLIDS: daily production of cheese solids.

Variance components, and estimates of heritability and herd incidence

Point estimates and features of the marginal posterior densities for the additive genetic variance and heritability ($h^2)\mathrm{ }$ of traits are reported in Table 2. The intra-herd heritability estimates of EME traits was: 0.25 for CH₄/DMI, 0.12 for CH₄/CM, 0.17 for dCH₄, 0.16 for CH₄/CY_CURD, 0.20 for CH₄/CY_SOLIDS and 0.12 for dDMI_est. As the incidence of herd-date variance on total variance ($HD$) was very high for all traits (46 to 71%), the corresponding across-herd heritabilities were lower, varying from 0.06 to 0.09 (data not shown). Heritability estimates for the daily production traits (dMY, dCMY, dCY_CURD and dCY_SOLIDS) varied from 0.12 (dCMY) to 0.15 (dCY_CURD), whereas the value of $HD$ for these traits was around 50%. The quality and technological traits of the milk samples had higher $h^2\mathrm{ }$ values than the yield traits (varying from 0.13 for fat percentage to 0.28 for protein percentage), and a relatively low incidence of $HD.$ The informative FAs used for predicting EME exhibited a high incidence of $HD$ and relatively low values of $h^2,$ varying from 0.07 (butyric acid) to 0.21 (linoleic acid). All $h^2$ estimates are characterised by a relatively low range of HPD95. Heart girth and BW showed intermediate within-herd heritability (0.27 and 0.25, respectively) and BCS slightly higher (0.33); their HD were similar to those of milk quality traits (0.23 to 0.25).

Table 2. Marginal posterior densities of genetic variance (${\sigma }_{\mathrm{a}}^2$), herd-date variance (${\sigma }_{\mathrm{h}}^2$), residual variance (${\sigma }_{\mathrm{e}}^2$), heritability $(h^2)$ and herd-date variance as proportion of total variance $(\mathrm{HD})$ for all the investigated traits.

Item	Variances			$h^2$		$\mathrm{HD}$
	${\sigma }_{\mathrm{a}}^2$	${\sigma }_{\mathrm{h}}^2$	${\sigma }_{\mathrm{e}}^2$	Mean	HPD95^a	Mean
Estimated Methane emissions^b
CH4/DMI, g/kg	0.18	1.78	0.53	0.25	0.10; 0.43	0.71
CH4/CM, g/kg	0.12	1.45	0.94	0.12	0.01; 0.28	0.57
dCH4, g/d	916	4166	4371	0.17	0.02; 0.28	0.44
CH4/CYCURD, g/kg	17	149	88	0.16	0.03; 0.32	0.59
CH4/CYSOLIDS, /kg	82	497	317	0.20	0.05; 0.41	0.55
dDMIest, kg/d	1.7	13.3	11.6	0.12	0.02; 0.28	0.50
Informative milk FAs (%FA)^c
4:0 - Butyric acid	0.024	0.436	0.299	0.07	0.01; 0.20	0.57
16:0iso - Iso-Palmitic acid	0.001	0.004	0.003	0.18	0.07; 0.33	0.53
18:1t10	0.001	0.005	0.004	0.20	0.04; 0.43	0.51
18:1t11 - Vaccenic acid	0.011	0.089	0.041	0.20	0.06; 0.39	0.63
18:1c9 - Oleic acid	0.980	1.719	5.993	0.14	0.02; 0.33	0.20
18:2c9,c12 - Linoleic acid	0.025	0.247	0.091	0.21	0.08; 0.42	0.68
Quality traits^d
Milk fat, %	0.054	0.103	0.356	0.13	0.03; 0.29	0.20
Milk protein, %	0.023	0.023	0.058	0.28	0.12; 0.49	0.22
%CYCURD, %	0.513	0.920	1.568	0.25	0.11; 0.42	0.31
%CYSOLIDS, %	0.143	0.155	0.468	0.23	0.09; 0.41	0.20
Daily yield traits^e
dMY, kg/d	2.85	22.84	19.80	0.13	0.02; 0.27	0.50
dCMY, kg/d	3.35	25.21	24.28	0.12	0.02; 0.24	0.47
dCYCURD, kg/d	0.085	0.628	0.484	0.15	0.01; 0.39	0.52
dCYSOLIDS, k/d	0.021	0.143	0.120	0.15	0.04; 0.29	0.50
Body size and condition
Heart girth, cm	13	15	35	0.27	0.07; 0.51	0.24
Est. body weight, kg	438	601	1340	0.25	0.08; 0.44	0.25
Body condition score	0.026	0.021	0.050	0.33	0.17; 0.51	0.23

^a HPD95: bounds of the 95% high posterior density interval.

^b Estimated Methane Emissions: CH₄/DMI, g/kg: methane yield, emitted per kg DMI; CH₄/CM, g/kg: methane intensity per kg fat and protein corrected milk produced; dCH₄, g/kg: daily methane production per cow; CH₄/CY_CURD, g/kg: methane intensity per kg of fresh cheese produced; CH₄/CY_SOLIDS g/kg: methane intensity per kg of cheese solids produced; dDMI_est, kg/d: estimated daily DMI of cows.

^c Informative milk FAs are the fatty acids included as independent variables in the equations used to estimate the enteric methane emissions (van Lingen et al. 2014).

^d Quality traits: %CY_CURD: wt of fresh cheese as % of processed milk; %CY_SOLIDS: wt of cheese solids as % of processed milk.

^e Daily yield traits: dMY: daily milk yield; dCMY: daily fat and protein corrected milk yield; dCY_CURD: daily production of fresh cheese per cow; dCY_SOLIDS: daily production of cheese solids.

Sensitivity analysis

Table 3 shows the results of the sensitivity analysis on $h^2$ and HD estimates. The only factors affecting the $h^2$ estimates (average $h^2$ of reduced datasets being at least 1.0 SD unit above or below the corresponding average estimate found on the entire dataset, shown in bold) were: the dairy system on CH₄/CM, whose $h^2$ was increased by both the exclusion of traditional farms with tied cows or the exclusion of modern farms using TMR, and the stage of lactation because removing late lactation cows led to a decrease of CH₄/DMI and an increase of the 3 methane intensities. The variations observed for HD effect were much smaller. The only relevant change was observed for dCH₄ when traditional farms with tied cows were removed from the analysed dataset: in this case the herd variability was reduced from 0.44 to 0.30 (Table 2).

Table 3. Sensitivity analysis of the heritability ($h^2$) and herd effect (HD) obtained excluding alternatively from the entire dataset of phenotyped cows those belonging to the extreme dairy systems (traditional farms with tied cows or modern farms using total mixed rations), or parities (primiparous or old cows), or lactation stages (fresh cows and late lactation cows).^a

Item	Entire dataset	All dairy systems excluded		All parities excluded		Entire lactation excluded
		Tied cows	Modern TMR	1st	≥ 4th	DIM <60 d	DIM >240 d
Cows with EME, N.	1,091	753	746	768	817	926	762
CH4/DMI, g/kg
h^2, mean	0.25	0.29	0.31	0.34	0.31	0.25	0.17
h^2, SD	0.08	0.10	0.11	0.12	0.11	0.08	0.08
HD, mean	0.71	0.68	0.69	0.70	0.69	0.73	0.70
CH4/CM, g/kg
h^2, mean	0.12	0.25	0.21	0.14	0.09	0.09	0.29
h^2, SD	0.06	0.11	0.09	0.09	0.04	0.04	0.14
HD, mean	0.57	0.59	0.55	0.54	0.56	0.64	0.54
CH4/CYCURD, g/kg
h^2, mean	0.16	0.17	0.22	0.15	0.18	0.16	0.39
h^2, SD	0.08	0.08	0.09	0.08	0.09	0.07	0.13
HD, mean	0.59	0.56	0.58	0.55	0.58	0.61	0.58
CH4/CYSOLIDS, g/kg
h^2, mean	0.20	0.24	0.30	0.22	0.19	0.15	0.39
h^2, SD	0.09	0.10	0.12	0.10	0.09	0.08	0.12
HD, mean	0.55	0.54	0.54	0.52	0.55	0.59	0.53
dCH4:
h^2, mean	0.17	0.23	0.24	0.21	0.21	0.16	0.15
h^2, SD	0.08	0.10	0.10	0.10	0.11	0.09	0.08
HD, mean	0.44	0.30	0.41	0.46	0.44	0.44	0.46

^a Boldface indicates estimated average $h^2$ values of reduced datasets being at least 1.0 SD unit above or below the corresponding estimated h² value of the entire dataset.

Estimated Methane Emissions (EME): CH₄/DMI, g/kg: methane yield, emitted per kg DMI; CH₄/CM, g/kg: methane intensity per kg fat and protein corrected milk produced; CH₄/CY_CURD, g/kg: methane intensity per kg of fresh cheese produced; CH₄/CY_SOLIDS g/kg: methane intensity per kg of cheese solids produced; dCH₄: daily methane production per cow.

Additive genetic correlations among traits

Table 4, shows the point estimates of the marginal posterior densities of the additive genetic correlations among the estimated EME traits, and between them and the traits used for their estimation. All estimates with P greater than 95%, indicating relevance in a Bayesian framework, are shown in bold. Phenotypic, herd, and residual correlations have been also obtained and are reported in supplementary material (Tables S1, S2 and S3, respectively).

Table 4. Mean of the marginal posterior distribution of additive genetic correlations among the estimated EME traits and between EME traits and informative milk fatty acids (FAs), milk quality, daily yield, and body size traits.^a

Item	Yield	Intensity			Production
	CH4/DMI	CH4/CM	CH4/CYCURD	CH4/CYSOLIDS	dCH4
Estimated EME traits^b
CH4/CM, g/kg	0.53	–
CH4/CYCURD, g/kg	0.78	0.68	–
CH4/CYSOLIDS, /kg	0.76	0.76	0.94	–
dCH4, g/d	0.02	0.23	0.10	0.27	–
Informative milk FAs (%FA)^c
4:0 - Butyric acid	0.03	−0.72	−0.19	0.07	0.04
16:0iso - Iso-Palmitic acid	0.64	0.03	0.45	0.56	−0.49
18:1t10	−0.43	−0.64	−0.39	−0.10	−0.65
18:1t11 - Vaccenic acid	−0.06	−0.48	−0.24	−0.31	−0.86
18:1c9 - Oleic acid	−0.25	−0.89	−0.54	−0.79	−0.70
18:2c9,c12 - Linoleic acid	−0.92	−0.65	−0.78	−0.66	−0.05
Quality traits^d
Milk fat, %	−0.64	0.15	−0.19	−0.24	0.27
Milk protein, %	−0.53	−0.03	−0.35	−0.47	−0.23
%CYCURD, %	−0.81	0.00	−0.64	−0.54	0.26
%CYSOLIDS, %	−0.70	−0.17	−0.49	−0.62	0.04
Daily yield traits^e
dMY, kg/d	0.08	−0.19	−0.03	0.17	0.84
dCMY, kg/d	−0.28	−0.26	−0.28	−0.14	0.95
dCYCURD, kg/d	−0.48	−0.18	−0.50	−0.24	0.91
dCYSOLIDS, k/d	−0.33	−0.16	−0.41	−0.40	0.86
Body size and condition
Heart girth, cm	0.32	0.20	−0.18	0.09	−0.63
Estimated body weight, kg	0.31	0.40	0.26	0.44	−0.65
Body condition score	0.03	0.14	−0.09	−0.17	−0.74

^a Boldface indicates additive genetic correlations with >95% of posterior probability accumulated above 0 (positive estimates) or below 0 (negative estimates).

^b Estimated Methane Emissions: CH₄/DMI, g/kg: methane yield, emitted per kg DMI; CH₄/CM, g/kg: methane intensity per kg fat and protein corrected milk produced; dCH₄, g/kg: daily methane production per cow; CH₄/CY_CURD, g/kg: methane intensity per kg of fresh cheese produced; CH₄/CY_SOLIDS g/kg: methane intensity per kg of cheese solids produced; dDMI_est, kg/d: estimated daily DMI of cows.

^c Informative milk FAs are the fatty acids included as independent variables in the equations used to estimate the enteric methane emissions (van Lingen et al. 2014).

^d Quality traits: %CY_CURD: wt of fresh cheese as % of processed milk; %CY_SOLIDS: wt of cheese solids as % of processed milk.

^e Daily yield traits: dMY: daily milk yield; dCMY: daily fat and protein corrected milk yield; dCY_CURD: daily production of fresh cheese per cow; dCY_SOLIDS: daily production of cheese solids.

Estimated Methane Emissions (EME): CH₄/DMI, g/kg: methane yield, emitted per kg DMI; CH₄/CM, g/kg: methane intensity per kg fat and protein corrected milk produced; CH₄/CY_CURD, g/kg: methane intensity per kg of fresh cheese produced; CH₄/CY_SOLIDS g/kg: methane intensity per kg of cheese solids produced; dCH₄: daily methane production per cow.

In general, the estimated methane yield and the 3 estimated methane intensity traits are all positively correlated with each other from the phenotypic, genetic, herd and residual point of view. Estimated methane production, on the other hand, had much lower correlation coefficients with the other EME traits, and sometime negative correlation coefficients with estimated methane yield.

As expected, the estimated EME traits were correlated with the informative FAs used for their estimation, and with the proper sign (positive for iso-palmitic acid, negative for the others). Vaccenic acid was the only FA almost not correlated with the majority of estimated EME traits.

The correlations between milk quality traits (i.e. composition, and percentage cheese yields) and the estimated EME traits were generally moderate to low and of varying signs. In fact, the phenotypic, genetic, herd and residual correlations were generally lower than 30% in absolute value, except for those between the two %CYs and the methane intensities that they were jointly used to estimate (CH₄/CY_CURD and CH₄/CY_SOLIDS). Important exceptions were the high negative additive genetic correlations between estimated methane yield and all the milk quality traits (−0.53 to −0.81, p > 95%), and between estimated methane intensity per unit of cheese solids and milk protein content (−0.47, p > 95%).

Daily milk and cheese yield traits were all, as expected, highly positively correlated with estimated daily methane production from the phenotypic, genetic, herd and residual point of view. In contrast, they were negatively correlated with the estimated methane yield and intensities.

Lastly, from the genetic point of view, body size and condition traits were weakly correlated with the 4 methane yield and intensity traits, whereas they were negatively correlated with methane production (Table 4).

Discussion

Heritability and herd incidence of informative traits used for estimating EME traits

Heritability estimates of milk traits used to estimate EME traits ranged from 0.07 to 0.28, in line with previous findings obtained from the same database. The genetic parameters of milk and cheese-making traits were reported in detail by Bittante et al. (2013). The milk fatty acid profiles were described and discussed by Pegolo et al. (2016, 2017), and Cecchinato et al. (2019). It is worth noting that none of the FAs van Lingen et al. (2014) included in their equations for predicting EME traits are synthesised in the cow’s udder. Moreover, with the exception of oleic acid, these FAs are present in milk in low proportions and have large phenotypic variability mainly due to the large differences among different herds (with different feeding practices), as reported in Table 2 with the incidence of $HD.$ The intra-herd heritability of these FAs, with the only exception of butyric acid, is of the same order as the milk yield and quality traits. It is worth noting that the heritability estimates in Holsteins were also similar to (iso-oleic and vaccenic acids) or greater than (butyric acid) our estimates (van Engelen et al. 2015), and that $HD$ was very high for the two former FAs (56–60%), and much lower for the latter (16%).

Oleic acid (18:1c9), on the other hand, is present in milk in large proportions, has a much lower coefficient of phenotypic variability, and is less affected by herd than the other FAs (Table 2). As discussed in the previous phenotypic study (Bittante et al. 2018), oleic acid has been frequently found to have an inhibitory effect on rumen fermentations, particularly those producing acetate and butyrate (Chilliard et al. 2000), that are the basis of fat synthesis in the cow’s udder. It is not surprising, then, that the oleic acid proportion in milk is negatively correlated with de novo synthetised FAs. In a multivariate factor analysis this FA is included, with a negative sign, in a latent explanatory factor together with the main de novo FAs (Mele et al. 2016). But oleic acid, during negative energy balance, could be originated also from body fat mobilisation.

Heritability and herd incidence on EME traits estimated from milk FA profile and cheese traits

Methane in the fore-stomach and intestines of ruminants is produced by extensive microbial fermentation of ingested and ruminated feedstuffs. There is no known bovine gene directly involved in enteric methane production, but there is a growing body of research on the complex, reciprocal relationships between the organism of the ruminant host and its gastro-intestinal microbiota (Weimer et al. 2010; Kiaosa-Ard and Zebeli 2013), also known as the animal’s ‘hidden organ’. In light of this, we may speak of ‘indirect heritability’ quantifying the indirect effect of the bovine genome on rumen production of methane by microorganisms.

Some studies have been carried out on the EME genetics of bovine species in different conditions, and with different populations and EME quantification techniques. As recently reviewed by de Haas et al. (2017) and by Brito et al. (2018), no estimates have been made of the genetic parameters of dairy cows for EME traits measured using the gold standard method, respiratory chambers. Moreover, this method, which requires highly sophisticated experimental facilities, is not expected to be used in the near future for genetic studies on large numbers of cows and it cannot be used in commercial herds in field conditions. In light of this the estimates of genetic parameters of EME available, included ours, should be taken with prudency. Two approaches are possible: i) the use of less accurate direct methods; or ii) the use of indirect methods correlated with or calibrated on direct methods (possibly using respiratory chambers).

Direct methods, suitable for field conditions, involve the use for analysing breath/eructate gas composition during milking or feeding are much less accurate and are the object of a strong scientific debate between authors that criticise heavily the use of these methods (Hristov et al. 2018; Huhtanen and Hristov 2018) and authors affirming ‘… that even if measurements are inaccurate, imprecise, or biased, they might provide valuable information for selective breeding.’ (Bovenhuis et al. 2018). Using these systems, Lassen and Løvendahl (2016) estimated a heritability value of 16% for the CH₄/CO₂ ratio of the air, which is the phenotype measured by the apparatus, whereas van Engelen et al.(2018) found for the same trait a repeatability of 14% and a heritability of only 3% (31% and 11%, respectively, after log transformation). Similar questions have been raised for infra-red analysers placed in feeder bin (Wu et al. 2018). Moreover, these methods sometimes yield different results among them (Hristov et al. 2016), and respect to respiration chambers in response to different cow’s diets (Hammond et al. 2016), as recently reviewed by Hristov et al. (2018).

As a comparison between these methods of direct measurement of gas composition of air/breath and an indirect method like that based on milk FA profile and used in this study is not available, future research on this issue is needed.

Among the indirect methods, some involve measuring the cow’s feed intake and assessing the correlation between feed intake and methane production. De Haas et al. (2011) estimated a heritability of predicted methane production of 35%, close to that of residual feed intake (40%). Methane intensity (g/kg CM) was also estimated and found to have very high heritability (58%), but it should be considered that all cows were reared in one experimental farm. Pickering et al. (2015), using similar methods, estimated h² of methane production to be 13%. The assumption that methane yield and methane intensity benefit from an increase in feed efficiency (de Haas et al. 2011; Negussie et al. 2017), was recently demonstrated not corresponding to the experimental results obtained in the rumen, where the most efficient cows showed a different microbial populations and more intense feed fermentation and digestion (Flay et al. 2019).

Moreover, measurement of DMI and body tissue deposition/mobilisation for calculating residual feed intake and feed efficiency is not easy at field level and can explain some inconsistency in relations between feed efficiency and expected EME traits (Negussie et al. 2017).

Predictions of EME based only on simulations based on milk yield and composition, expected body weight, and expected daily feed intake (Yin et al. 2015) depend only on equations used, and could suffer from circularity of arguments. Furthermore, simulations based only on milk yield and body weight are not useful for genetic selection of EME traits because they do not give any information on the actual EME production of a given cow respect to predicted one (residual methane emission), which should be the true objective of selection.

Indirect prediction of methane yield and/or methane intensity from measured milk samples characteristics (FAs or FTIR spectra absorbance) are based on predictors specific of a given cow and independent from its milk yield or body weight and could be much more useful to the dairy sector (Bittante and Cipolat-Gotet 2018).

Kandel et al. (2017) used calibration equations set up using SF₆ gas tracer on the milk FTIR spectra. Methane production estimated to have a heritability coefficient of 22–25%, and methane intensities of 17–18%, Shetty et al. (2017) also used FTIR milk spectra to predict the CH₄/CO₂ ratio, measured with IR sniffers, but the results were modest and the authors judged this method to be unfeasible for predicting EME traits.

Only van Engelen et al. (2015), like the present study, used measured milk fatty acid profiles to predict the methane yield (g/kg DMI) of 1,905 primiparous Holstein cows. These authors compared three calibration equations based on different FAs from three single feeding trials (50 observations) in respiration chambers (Dijkstra et al. 2011). Even though the predicting equations were characterised by relatively high accuracy (R² 0.63 to 0.73), the heritability estimates of these three predictions were very different from each other: 12%, 20% and 44%, respectively. It should be considered that those equations were based on only 50 observations of 3 experiments in one research centre. Our intra-herd heritability estimate for methane yield is 24% (Table 2). The across-herd heritability is obviously lower, especially for the traits characterised by the greater herd effect. Unlike van Engelen et al. (2015), we used prediction equations, also based on respiration chamber measurements, obtained from a meta-analysis (van Lingen et al. 2014), which combined data from 2 countries/research centres, 8 trials, 30 diets, and 146 observations. The results obtained should in any case be evaluated with prudence because the EME prediction equations used in this study (van Lingen et al. 2014) were characterised by moderate accuracy (R² 0.47 to 0.54), as often happens with meta-analysis on large and heterogeneous datasets. We preferred to use these latter equations because the variety of diets and feed supplements tested, together with the similar level of daily milk yield, appeared to warrant a better transferability to the conditions of our territory, as discussed in details in our previous work (Bittante et al. 2018). Unfortunately, we cannot compare the results obtained with the van Lingen’s et al. (2014) equations with those, also from a meta-analysis, proposed by van Gastelen et al. (2018) because of the lack in our database of one fatty acid per equation.

Effect of dairy system, parity and lactation stage on heritability of EME traits

The focal question regarding the representativeness of EME estimated using methods developed on a limited number of precise observations obtained in the respiratory chambers, i.e. an environment very different from commercial farms, cannot be solved directly. In the previous phenotypic analysis on this dataset (Bittante et al. 2018), we found that the results obtained were well in agreement with expectations from different dairy systems, parities and lactation stages based on animal’s physiology. Here we tested indirectly this question at a genetic level through a sensitivity analysis (Table 2). The Bayesian approach used allowed not only to obtain the central parameter (average) of the heritability estimates, but also their distribution. Using a relevance threshold of being at least 1.0 SD unit above or below the average heritability obtained on the entire dataset, very few estimates obtained on reduced datasets presented relevant differences respect to the reference value.

The exclusion from the dataset of ‘extreme’ farms, like those very traditional with tied cows or those modern with TMR, tended to increase the heritability values of all the EME traits, even though a statistical relevance was reached only for methane intensity-milk, whose heritability almost doubled respect to the reference on the entire dataset. If this could be due to a too large environmental variability (see HD values) or to some G × E interaction cannot be addressed and requires further research. In any case it seems that a more homogeneous dairy system could allow a greater accuracy of genetic evaluation of animals for EME traits and particularly CH₄/CM.

The effect of excluding very young or very old cows from the dataset was not relevant, whereas the stage of lactation affects heavily all the EME traits. In particular, in early lactation the emissions are much lower than later, as confirmed in the phenotypic study (Bittante et al. 2018). This is due to the fact that a relevant proportion of energy requirement of the fresh high yielding cow is covered by the mobilisation of body fatty reserves (negative energy balance), without any corresponding release of methane from the rumen. Also genetic parameters are affected by lactation stages, as shown by Vanrobays et al (2016) on infra-red predicted EME traits. Anyway, the exclusion of the cows in early lactation, and then of those possibly being in severe negative energy balance (and seldom included in respiration chamber experiments), seems not affects the heritability estimates of any EME traits. On the contrary, the exclusion from the dataset of the cows in late lactation exerted a relevant effect on the estimation of heritability of almost all EME traits (only methane production heritability was not affected). The heritability of methane yield was reduced, whereas that of the three methane intensities was strongly increased (Table 3). It is worth noting that in this phase of lactation generally milk yield is decreasing, like feed intake, whereas the cow is reconstituting the fatty body reserves. The entity and proportions of these phenomena depends heavily on the fact that the cow is pregnant or not and, in the first case, on the stage of pregnancy. The pregnancy in itself requires energy: the quantity of energy deposited in foetus and foetal annexes in lactating cows is not much high, but the efficiency of deposition of energy for pregnancy is particularly low (NRC 2001), meaning large energy losses, and how much they involve methane losses is not known and require further research.

The differences between the studies regarded also the breed of cows as we sampled Brown Swiss cows whereas the other studies were carried out on Holsteins. Even though it seems not probable that the breed of dairy cows with similar feeding regime and productive level affect the microbial activity in the rumen, it should be noted that Xue et al. (2011), comparing Jersey crossbreds with Holsteins purebreds, found many productive and metabolic differences between the two genotypes, but they observed that their energy loss from EME as a ratio of total energy intake remained very similar, also using diets of very different roughage: concentrate ratio.

Correlations between traits and selection opportunities for EME phenotypes

Knowledge of the genetic correlations among EME traits, and between them and their predictors is fundamental, especially where the aim is to design an optimal selection index, and only a few studies have reported these genetic parameters, almost always with large standard errors of estimates. What is still lacking is a reliable estimation of the genetic correlations between informative FAs and EME obtained with the golden standard method of respiration chambers that can go beyond the phenotypic correlations obtained by van Lingen et al. (2014). The knowledge of genetic correlations would allow the direct inclusion of informative fatty acids in a selection index aimed at improving ecological footprint of dairy sector. We found that milk yield is genetically correlated only with estimated daily methane production, although this is mainly because the latter is obtained from the former and because increased milk yield imply larger DMI and rumen loads. Milk quality, on the other hand, is negatively (favorably) correlated with methane yield, but not with methane intensity (with the obvious exception of the two %CY traits, and partly of milk protein content, with methane intensity per unit of cheese).

The possibility to use the heritabilities of estimated EME traits, the genetic correlations with other milk traits, and the availability of direct and indirect methods of prediction, other than those using respiration chambers, in large scale surveys, may open the way to genetic and/or genomic selection of dairy populations.

Current breeding programmes in many countries have an unfavourable effect on methane production (g/d), but a favourable effect on methane intensity (g/kg CM), as also demonstrated by Pryce and Bell (2017). In fact, selection for greater milk yield, and consequently greater feed intake and possibly live weight, is obviously leading to an increase in daily methane production per cow. The inclusion of methane production (g/d) with a negative weighing within a selection index would reduce the genetic progress in daily milk production because of the positive (unfavourable) genetic correlation between these two traits. The reduction of methane production per cow could be counterbalanced by an increase of the number of cows needed to fulfil the dairy-chain needs (Pryce and Bell. 2017). On the contrary, the increase in milk yield is said to increase cow’s efficiency and dilute methane production (g/d), thereby reducing methane intensity (g/kg CM), even though recent feeding experiments questioned this assumption (Flay et al. 2019). The main problem facing selection for improving cow’s efficiency and also for reducing daily methane production is that an inaccurate estimation of the variation of body reserves will lead to an unwanted selection for thin cows, increasing in this way fertility and longevity problems in the long time. For this reason, we included also body size and condition among the traits considered in this study. As shown in Table 4, the EME traits estimated from milk fatty acids were not relevantly correlated, from the genetic point of view, with the cow’s body size and condition, with the relevant exception of the negative correlations with methane production. This does not mean that increasing body size and condition will improve ecological footprint of dairy farms, but simply that high producers are less heavy and fat than low producers because of fatty body depot mobilisation. The genetic correlations between daily milk yield and hearth girth (−0.64), estimated live weight (−0.73) and BCS (−0.78), in fact, are all high and negative (data not shown).

If the market demand for milk is the leading driver of the dairy industry, methane intensity, and not methane production, should be the selection objective of the dairy chain. Moreover, where the dairy sector is oriented mainly to cheese production, like in Italy and many other European countries, greater account should be taken of methane intensity per unit of cheese in an attempt to optimise the relationships between economic and ecological sustainability in the dairy industry. It should also be considered that the new methane intensity traits (Table 1) calculated in terms of predicted EME per kg of fresh cheese (h² = 0.16) or per kg of cheese solids (h² = 0.20) exhibited greater heritabilities than the predicted EME per kg CM (h²= 0.13).

The economic value of EME traits and the breakeven prices of recording them were also modelled by Hansen Axelsson et al. (2015). These authors demonstrated that, when the cow’s entire career is taken into account, and not only its lactation periods, a viable and economic alternative to EME phenotyping of dairy cows is to use stayability after the 1st calving as an indirect indicator to dilute the EME of the replacement heifer on a longer productive lifespan.

Indirect EME predictions could be used for properly calibrating genomic selection. Genome-wide associations and genomic selection of EME traits have been more extensively studied in beef populations than in dairy cattle (Hayes et al. 2016; Manzanilla-Pech et al. 2016), mainly with an eye to the possibility of directly measuring these traits in candidate young bulls. In any case, Lassen et al. (2016) found very similar heritability values for methane production and methane intensity when they were estimated on the basis of pedigree information or genomic relationships. They also found some interesting correlations between fatty acid profile and EME traits, although they pointed out the need for further research on larger numbers of cows.

Conclusions

In the present study, we estimated the genetic parameters of EME traits obtained from equations based on a meta-analysis of their relationships with milk FA profiles. Heritability of estimated methane production (g/d) was similar to that of milk yield, whereas heritability of estimated methane intensity (g/kg CM) was similar to that of fat content. Estimated methane intensities per kg of fresh cheese and per kg of cheese solids, of interest to countries where a large proportion of the milk produced goes into cheese making, were more heritable than those per kg of corrected milk. Lastly, estimated methane yield (g/kg DMI) was the most heritable trait, together with the protein content and cheese yield of milk. Although results should be validated on larger population and different breeds, our estimates indicate the feasibility of selecting dairy cows for the improvement both milk production and the ecological footprint of dairy farming. However, further research is needed to study the representativeness of these indirect methods for estimating EME traits and especially on the effect of different dairy systems and feeding regimes and the complex interrelationships with feed efficiency, body reserve management and cow's fitness.

Ethics statement

The cows included in this study belonged to commercial private herds and were not subjected to any invasive procedures. Milk and blood samples were previously collected during the routine milk recording coordinated by technicians working at the Breeder Association of Trento Province (Italy) and therefore authorized by a local authority.

Acknowledgments

The authors thank the Superbrown Consortium of Bolzano and Trento (Trento, Italy) for their help with sampling and recording; the National Association of Brown Swiss Breeders (ANARB, Bussolengo, Verona, Italy) for their cooperation and for supplying the pedigree information; and Prof. M. Mele and G. Conte (Dipartimento di Scienze Agrarie, Alimentari, Agro-ambientali, Università di Pisa, via del Borghetto 80, 56124 Pisa, Italy) for their help in analysing and discussing the detailed fatty acid profiles.

Disclosure statement

No potential conflict of interest was reported by the authors.