A method for implementing methane breeding values in Australian dairy cattle

Keywords: methane emission, sustainability, selection index, index weights.

It is recognised that within dairy systems, the largest source of greenhouse gas (GHG) emissions is from enteric fermentation and methane (CH₄) production (g/day). Over 30 years (between 1980 and 2010), Moate et al. (2016) estimated that production of enteric CH₄ has been almost static (185 000 t in 1980 versus 182 000 t in 2010). At the same time, milk production has increased, so that the intensity of CH₄ emissions (i.e. emissions per unit of product, in this case, milk) has declined considerably by 40%, from ~33.6 g CH₄/kg milk to 19.9 g CH₄/kg milk. Since milk production traits are a large part of the breeding goal, it follows that the reduction in GHG has been, in part, the result of genetic improvement in efficiency and dilution of emissions per litre of milk produced.

The Australian national breeding objective, known as the Balanced Performance Index (BPI), includes traits that contribute to cow profitability, farmer preferences and desired gains (Byrne et al. 2016). The BPI encompasses biological traits associated with milk production, longevity, fertility, feed efficiency and health. Feed-efficiency breeding values, known as feed saved, include genomic breeding values for residual feed intake and breeding values for maintenance requirements based on predicted cow bodyweight (Pryce et al. 2015). Thus, by selecting for production traits and feed saved simultaneously, it is anticipated that gross efficiency will improve in Australian dairy cows (Byrne et al. 2016).

There has been a lot of interest in recent years in developing estimated breeding values (EBVs) to reduce GHG from the livestock sector. However, a major limitation has been the quantity of high-quality CH₄ phenotypes measured on individual cows required to develop these EBVs. Currently, most published reliabilities are <10% (Manzanilla-Pech et al. 2021), which is insufficient for implementation of genomic prediction. However, it has been recognised that selecting for improved efficiency of milk production, longevity and fertility has had a beneficial effect on GHG emissions (Løvendahl et al. 2018). Consequently, it might be possible to develop an index that includes the contribution of these traits to GHG emissions, either separately, or in addition to CH₄ emissions. Breeders with an interest in selecting for reduced GHG emissions, without greatly affecting profitability, can select for high-BPI and low-GHG emitters.

An approach that could be suitable is to estimate the increment in CO₂-eq per unit change in EBVs that are current selection objectives, while accounting for interactions among traits, thereby calculating trait-specific environmental coefficients that are independent, such as those in Amer et al. (2018) and Richardson et al. (2021a) and that, for selection index purposes, are superior to non-independent coefficients, such as those calculated by Bell et al. (2013) and Pryce and Bell (2017).

The aim of this research was to (1) develop a GHG subindex using the gross CO₂-eq values derived by Richardson et al. (2021a) as subindex weights for the EBVs most strongly associated with GHG emissions (milk, fat and protein yield, longevity, and feed saved), (2) modify the GHG subindex to include residual CH₄ EBVs, as defined by Richardson et al. (2021b), from measured CH₄ data phenotypically corrected for energy-corrected milk and (3) validate the GHG subindex with CH₄ phenotypes corrected for fixed effects. Finally, we explored the effect of current selection objectives on CH₄ by regressing the derived subindexes for GHG on BPI and Health Weighted Index (HWI).

Previously calculated CH₄ coefficients (kg CO₂-eq/unit change in trait) were used as weights and applied to EBVs commonly used in selection and most strongly associated with emissions to derive two possible subindexes aimed to rank the environmental impact of individual animal on the basis of their genetic merit. The first index includes traits that are currently included in the national breeding objective and known to have a direct environmental impact, including production traits, survival and feed efficiency. The second includes these index traits as well as a direct CH₄ trait. The two subindexes were correlated with current national breeding indexes (BPI and HWI), as well as additional traits of interest.

Currently, there are two main indexes used for ranking dairy cattle in Australia. The BPI includes traits pertaining to production (milk, fat, protein), fertility, survival, health and feed efficiency. The HWI, which was developed using the same traits as in the BPI, places additional emphasis on traits related to health and efficiency, such as fertility and feed saved. Details of the index development are given in Byrne et al. (2016). DataGene (Melbourne, Vic., Australia) calculates the BPI and HWI routinely for all cows and bulls included in national genetic evaluations.

In total, 1712 individual cow CH₄ measurements were obtained from 464 cows measured over a 5-day period from 12 experimental batches across 5 years from 2013 to 2017. Cows were located at the Ellinbank SmartFarm (Melbourne, Vic., Australia) and milked twice per day. Measurements for CH₄ (g/day) were performed using the SF₆ tracer method previously described by Deighton et al. (2014). Records for CH₄ were averaged to obtain one observation per animal to account for day-to-day variation (Moate et al. 2016). Animals were between 68 and 187 days in milk at the start of measurement and across parities 1–9. Genotypes for the 464 cows used in this study were provided by DataGene Ltd. After editing by excluding genotype calls with a GenTrain score <0.6 using the methods described by Erbe et al. (2012), 47 162 single nucleotide polymorphism (SNP) markers were available for genomic analysis. EBVs for the 464 animals used in this study were estimated as part of the routine genetic evaluation service of DataGene Ltd and included milk, fat, protein, survival and feed saved, which are traits directly used in GHG subindex development.

Genotypes for the 3412 registered Holstein bulls used in this study were provided by DataGene Ltd, with processing and genotyping methods being consistent with the cow dataset. Bulls were born between 2010 and 2015. The EBVs used in this analysis included milk volume, milk fat, milk protein, survival and feed saved, as well as other traits of interest such as heat tolerance and liveweight and were accessed from the April 2020 official genetic evaluation run.

Richardson et al. (2021b) investigated nine definitions of residual CH₄ traits and concluded that CH₄ production corrected for energy-corrected milk is an appropriate selection candidate to reduce environmental impact without severely affecting other traits, such as production, health and fertility. Using this definition of residual CH₄ EBVs, we calculated genomic breeding values for bulls (or cows) that were not part of the reference dataset, by multiplying the vector of SNP effects with the genotype matrix of bulls. These SNP effects were derived using a process to back-calculate SNP effects from direct genomic values (DGVs) of cows in the reference population. In brief, using DGVs of CH₄ production corrected for energy-corrected milk of 464 Ellinbank cow reference population, 41 276 SNP effects were estimated as:

where Z is the genotype matrix (464 individuals × 41276 SNP), and a vector of descaled DGVs with a mean of 0 estimated using genomic best linear unbiased prediction (VanRaden 2008). The prediction equations of SNP effects were then used to predict residual CH₄ DGVs of the 3412 Holstein bulls.

Methane coefficients were previously calculated by Richardson et al. (2021a) on the basis of the approach used by Amer et al. (2018) and adapted to calculate the effect of a unit change in milk, fat, protein, feed saved, and survival traits on CO₂-eq emissions per cow and per kilogram of protein equivalents. Protein equivalents are a weighted aggregate of the product outputs from milk protein, milk fat, and milk volume weighted on the component value ratio relative to protein. Briefly, this method estimates the change in total emissions and product output caused by a 1 unit change in each index trait, resulting from either a direct emissions trait (CH₄ yield), changes in herd structure (fewer replacements), or the dilution effect of higher yields (milk production) and proliferation (more offspring/dam). The traits used in the current study were determined to have an independent effect on emissions, with the addition of fertility (Richardson et al. 2021a). However, as fertility is a primary reason for culling, the environmental impact of fertility is largely accounted for by the survival EBV, with minimal additional effects applying to extended lactations observed in seasonal calving systems (Workie et al. 2019; Richardson et al. 2021a). Therefore, only the survival CH₄ coefficient is considered in the index. The model was used in the current study to dynamically represent an Australian dairy herd and assess effects of changes in traits. The effects of a unit change of a trait on GHG emissions are shown in Table 1. For example, a unit change in milk protein is estimated to be associated with a GHG emission of 1.97 kg of CO₂-eq (3.70 kg DM/kg protein yield × 0.532 kg CO₂eq / kg DM = 1.97 kg CO₂eq/kg protein yield).

Table 1.  Responses in carbon dioxide equivalents (CO₂-eq) per unit change of key traits under selection
Values are means and standard deviations (σ) of breeding values in the bull population

Indexes for CO₂-eq emissions per cow per year were calculated by multiplying the EBVs of each trait in the breeding objective by the CH₄ coefficient (i.e. effect of a single unit change of the trait on CO₂-eq emissions). The sum of this is the total CO₂-eq index. EBVs selected for this study included milk, fat, protein, feed saved and survival (Amer et al. 2018; Richardson et al. 2021a). Gross CH₄ coefficients were applied to trait EBVs to calculate the GHG index as follows:

where w_i is the CH₄ coefficient (gross or intensity) of the ith trait and EBV_i is the estimated breeding value of the ith trait (milk, fat, protein, feed saved or fertility).

An additional index was developed that extends the GHG index described above through the addition of the EBV of residual CH₄ and was calculated as follows:

where EBV_RMP is the residual CH₄ EBV and GHG_index is described above.

The GHG indexes were correlated with the Australian national selection indexes, as well as with index traits currently estimated through the national genetic evaluation services (DataGene Ltd). Correlations were estimated using the cor.test() function in R.4.0.4 statistical programming (R Core Team 2013).

The GHG_index subindex was validated using a 12-fold cross-validation, where the correlation between GHG_index and CH₄ phenotypes for each experimental batch was independently estimated, that is, batch by batch. EBVs for index traits were provided by DataGene Ltd for milk, fat, protein, survival and feed saved on the 464 Ellinbank cows used in developing the residual CH₄ trait. For the 464 cows, first GHG_index was calculated and then the values were correlated with SF₆ CH₄ phenotypes. As described by Su et al. (2012), to calculate the predicted accuracies of the true breeding values, instead of the EBVs, the correlations between the GHG subindex and SF₆ CH₄ phenotypes was divided by the square root of the heritability of residual CH₄ production (h² = 0.21; Richardson et al. 2021b).

The size of the dataset was insufficient to conduct a validation analysis for the GHG⁺_index. This is because all data on 464 Ellinbank cows were required for development of the genomic prediction equation to estimate EBV_RMP, so an independent dataset was not available.

The GHG_index subindex was applied to 3412 Holstein bulls and had a mean value of 38.22 kg CO₂-eq, with a standard deviation of 76.39. When residual CH₄ was incorporated into the subindex, GHG⁺_index, the standard deviation was 76.62 with a mean of 37.65 kg CO₂-eq.

Correlations between the GHG subindexes, BPI and HWI are presented in Table 2. The two GHG-based indexes, with and without the RMP EBV, were highly correlated (0.99). The BPI was uncorrelated with the GHG_index subindex and GHG⁺_index subindex, with correlations of close to zero at –0.03 and –0.02 respectively. Favourable correlations were estimated between HWI and both GHG-based subindexes, with correlations of –0.35 with GHG_index and –0.36 with GHG⁺_index. These trends were consistent within the top 300 BPI bulls, where the average BPI of the top 30 BPI bulls was 333. Considering only bulls whose GHG index value was below zero, the average BPI for the top 30 BPI bulls decreased to 312, representing a low 21-point difference (Fig. 1).

Table 2.  Correlation between GHG indexes and BPI, HWI and other EBVs within the population of 3412 registered Holstein bulls

Fig. 1.  GHG_index subindex scores among the top 300 BPI bulls. The solid vertical line is the average BPI (333) in the top 30 bulls and the dotted vertical line in the average BPI (312) for the top 30 BPI bulls with low GHG index values.

Correlations between the GHG indexes and selected EBVs are presented in Table 2. The GHG_index and GHG⁺_index subindexes had strong favourable correlations to feed saved (–0.71 and –0.70 respectively), fertility (–0.52 and –0.51) and survival (–0.43 and –0.42), as well as to heat tolerance (–0.45 and –0.46). Low to moderate favourable correlations were observed between the indexes and mastitis resistance (–0.29 and –0.28) and RMP (0.01 and 0.07). The traits with the largest unfavourable correlation with both indexes were fat (0.66), protein (0.57) and liveweight (0.56). Milk, milking temperament, likability and overall type had low to moderately unfavourable correlations of 0.30, 0.17, 0.16 and 0.14 respectively.

When GHG index values were correlated with SF₆ CH₄ data (Fig. 2, Table 2), the average correlation of the 12 validation batches was 0.23 (ranging 0.07–0.41). This gives the GHG_index subindex an accuracy of ~0.50, assuming a moderate heritability of 0.21 (Richardson et al. 2021b).

Fig. 2.  Relationship between GHG_index and SF₆ methane phenotypes, corrected for days in milk, batch-year and parity, in 464 cows located at the Ellinbank SmartFarm (Melbourne, Vic., Australia).

We have developed two GHG-based subindexes that can be used to select for reduced emissions. The advantage of our approach is that it utilises the EBVs of traits currently estimated through the national genetic evaluation service and therefore does not require a direct CH₄ trait. Although we apply the method to Holsteins in the present paper, it can be adapted to derive weights for other dairy breeds, such as Jerseys and crossbreds, although the method should be validated in these breeds before implementation. This approach enables emission mitigation strategies through genetics to be implemented immediately, with the option of including a CH₄ trait when available. Our results showed that there is little to no advantage in including the current residual CH₄ trait, although this is expected to change with a higher-reliability residual CH₄ trait. It is expected that a residual CH₄ trait may be available in the future (Manzanilla-Pech et al. 2021; Richardson et al. 2021b). However, this is still a significant challenge, as CH₄ phenotypes are expensive and laborious to measure and combining international datasets is challenging, with multiple measuring apparatus and techniques being used in data CH₄ collection internationally, typically on a single breed (University of Guelph 2016).

The high correlation between the two calculated subindexes, GHG_index and GHG⁺_index, showed that the RMP trait has a small effect on the index. This was expected, as the RMP trait has a low accuracy and a standard deviation of 5.0, whereas, for example, the standard deviation of the survival component of the GHG subindex was 28.36. However, as the RMP EBVs become more accurate, it is expected that the effect of the RMP traits within the index will increase proportional to the standard deviation and that the two GHG subindexes will further differentiate. We assumed that because the GHG subindexes and residual CH₄ EBV are both measured in CO₂-eq and that the index components can be added together to calculate GHG_index⁺. However, further analysis is required to optimise the residual CH₄ EBVs within GHG_index⁺. Conceptually, the inclusion of RMP traits in the subindex is important (Knapp et al. 2014). The GHG subindex index will result in more environmentally efficient animals; however, selection on a direct CH₄ trait will capture the variation in CH₄ that exists between animals, resulting in a further reduction in CH₄ emissions. Further improvements to the GHG index and an increased availability of CH₄ phenotypes may help increase the accuracy of the index, and thereby the proportion of between-animal genetic variance captured by the index that can be utilised for selection. Although the effect of the RMP impact is minimal, these results have shown that we can use genetic selection to begin reducing emissions immediately. The GHG_index method offers an option to include a sustainability component in the national breeding objective, with minimal compromises to advances in BPI and profit. In dairy systems where it is not possible to obtain accurate CH₄ phenotypes, the GHG_index offers a valuable alternative that can be implemented simply and immediately, as it utilises EBVs currently available through the national genetic services.

The relationships between the GHG indexes (GHG_index and GHG⁺_index) and the current national selection indexes (BPI and HWI) are favourable. The strong favourable correlations between HWI and GHG subindexes support the hypothesis that selection for improved efficiency and survival is a viable mitigation option, as previous studies have identified (Wall et al. 2010; Pryce and Haile-Mariam 2020). Although the HWI and BPI are strongly correlated (0.93), some re-ranking does occur between the indexes. The HWI places greater selection emphasis on health and functional traits, such as 5% on feed saved and ~40% on fertility and survival. Comparatively, the relative index weight in BPI for feed saved, and fertility and survival are 2.5% and 30% respectively (DataGene). The difference in comparative weighting also explains the differences in correlation with the GHG subindexes. Both the GHG subindexes and HWI place additional weight on survival and efficiency traits. However, the GHG subindexes also penalise higher production (milk, fat and protein). Selection on the HWI would results in favourable gains in the GHG subindexes. For the BPI, the correlations with the GHG subindexes were negligible (–0.02). However, in the top 50 bulls, 10 ranked negative on the GHG subindexes and a considerable number of bulls within the top 300 rank negative or close to 0 (Fig. 1). When considering the HWI, the GHG subindexes have a favourable negative correlation (–0.36) meaning that the higher the HWI, the less GHGs are emitted. As expected, on the basis of the favourable correlation between HWI and GHG, all but one of the top HWI Holstein bulls (within the top 300 BPI bulls) had a negative GHG subindex value, with one bull that was positive having a value very close to zero. This presents an opportunity to select for bulls that are profitable and also low GHG emitters.

As expected, the GHG subindexes are strongly correlated with the traits used in the index development, as well as with traits independent of index development that are suspected to have mitigation properties. Fertility was not included in the GHG subindexes as the CH₄ coefficient previously calculated only accounted for a portion of the effects of fertility, with all effects related to a reduction in replacements being accounted for through the survival EBV (Richardson et al. 2021a). However, there is a strong favourable effect between the GHG subindexes and fertility. This large favourable effect is also seen with heat tolerance (Nguyen et al. 2017), a trait defined as a rate of decline in milk, fat and protein yields per unit increase in temperature–humidity index. Heat tolerance is not included in the BPI and has a strong negative correlation with production traits, similar to the GHG subindexes. This suggests that selection only for a lower GHG subindex would result in decreases in production. However, by combining the GHG subindex with HWI or BPI, we can ensure that both genetic gain in profitable traits and decreases in gross emissions are achieved. Historically, unfavourable correlations among traits have been overcome by simultaneous inclusion in selection indexes, such as with fat and protein yield, although genetic gain is seen at a slower rate (Miglior et al. 2005). Type and production traits are also unfavourably correlated with the GHG subindexes, whereas more functional traits such as calving ease and mastitis resistance are favourably correlated. This suggests that transitioning some selection emphasis from type and production traits is favourable for environmental impact, and as Miglior et al. (2017) suggested, the selection for robust animals is vital for the future sustainability of the dairy industry.

The validation methods of this research were limited by the small dataset size, as only the GHG_index subindex could be compared with the SF₆ CH₄ data. However, as the GHG_index and GHG⁺_index subindexes were highly related, it is expected that a similar accuracy may be obtained for the GHG⁺_index subindex. With an increase in dataset size to where a cross-validation may be conducted without jeopardising the accuracy of the RMP EBV or validity of the analysis, a validation may be also conducted on the GHG⁺_index subindex with a higher confidence. The reliability of the GHG_index is sufficient for including the subindex in selection strategies. However, the accuracy of the GHG_index to select for direct CH₄ is low. This is primarily due to the small dataset, as obtaining accurate CH₄ phenotypes on a large number of individual animals is challenging. In Australia, CH₄ phenotypes are currently measured on a reference population of 464 females, which is considerably lower than the suggested 5000-cow female reference population required for a moderately heritable trait to reach the same level of reliability (0.25) as the GHG subindex (González-Recio et al. 2014).

The GHG subindex may be used in practice to allow farmers to make environmentally conscious breeding decisions. The developed GHG_index subindex has many opportunities for implementation, with minimal compromise to profit (Fig. 1). Within the BPI, the GHG_index subindex could be included as an additional subindex and, although HWI is already favourably correlated with the GHG subindexes, inclusion of the GHG_index subindex in HWI would put additional weight on traits based on their environmental impact (i.e. favouring longevity and efficiency). Additionally, these weights may be used in culling decisions to differentiate between two animals of similar genetic merit, as the higher emitter may be removed from the herd or breeding stock. These GHG_index subindex values may be used on a magnitude or range basis. For example, the GHG_index subindex may be used to sort cows into high-, medium-, low- and negative-emitter ranges. While using only bulls or animals with negative GHG_index subindex values may not be possible at this time as this would result in increased inbreeding, farmers could use this range ranking system as a way of selecting for lower-emitting cows, leading towards a more environmentally conscious industry.

The traits within the GHG based subindexes are weighted based on their environmental impact. This weighting scheme is an alternative to the traditional economic index, which weights traits on the basis of profit. Alternative weighting options have been suggested, including social impact and farmer preference (Nielsen and Amer 2007). This has been implemented in the current Australian Breeding Objective through the HWI and previously through the (now discontinued) Type Weighted Index (Byrne et al. 2016). Using survey approaches described by Martin-Collado et al. (2015), weights were assigned to traits on the basis of the priority and preferences of farmers, which grouped farmers into three sections (profit based, health based, and type based). Multi-source weighting allows the weights within an index to be developed on the basis of a subset of conditions. Those conditions may be economic, social, environmental or preference based, with the final weight applied to a trait being an aggregate of these weights. This should be considered when including the GHG_index subindex in the national selection indexes (BPI and HWI) and not just as a supplementary selection decision tool. The GHG subindexes offer an additional option for weighting traits, on the basis of environmental impact. In theory, the GHG indexes are effectively weighting the traits twice. First, on the basis of economic analysis, and second, on the basis of environmental impact. However, in practice we are adjusting the index weight so that it is no longer purely economic based, but instead an aggregate weight of environmental and economic factors. This method may be especially useful to develop breeding strategies for traits with low economic values in breeds where expensive and laborious phenotypes, such as CH₄, are not currently collected. Ideally, further breed-specific validations should be conducted, but as the majority of data has been collected on Holsteins, validation in other breeds is not possible at this time. As interest to include traits with low direct economic value in national selection indexes grows, further investigation is required to determine optimal inclusion methods.

This research has proposed a GHG subindex using EBVs currently estimated through the national genetic service and previously derived index weights that describe the change in CO₂-eq associated with a unit change in each trait. The final set of EBVs selected for inclusion in the GHG subindex were milk, fat and protein yield, survival and feed saved as these traits had an independent effect on emissions. A further modification was to include a direct CH₄ trait in the GHG subindex, which was a RMP EBV. The two GHG subindexes were highly correlated, suggesting the RMP had a minimal impact on the GHG index, which was expected due to the low accuracy of the RMP EBV. The GHG subindexes had high reliabilities; however, the accuracy of the subindexes to predict CH₄ was reasonably low. The correlations between the GHG subindexes and current national selection indexes (HWI and BPI) suggest that selection strategies using HWI are expected to reduce emissions, whereas selection on BPI results in negligible changes in either GHG subindex. Type and production traits were unfavourably correlated with the GHG based subindexes, whereas traits associated with longevity, health, fertility and efficiency were favourable correlated. The GHG subindex may be included in the national selection index or used independently in practice by farmers to make environmentally conscious breeding and culling decisions, with minimal compromise to profit. Therefore, as the GHG subindex does not necessarily require a direct CH₄ trait to lower the environmental impact of dairy cattle, selection for more sustainable dairy cattle may be promptly implemented until sufficient data are collected on CH₄.

The authors declare no conflicts of interest.