Estimates of the genetic contribution to methane emission in dairy cows: a meta-analysis

The present study aimed to perform a meta-analysis using the three-level model to integrate published estimates of genetic parameters for methane emission traits [methane yield (METY), methane intensity (METINT), and methane production (METP)] in dairy cows. Overall, 40 heritability estimates and 32 genetic correlations from 17 papers published between 2015 and 2021 were used in this study. The heritability estimates for METY, METINT, and METP were 0.244, 0.180, and 0.211, respectively. The genetic correlation estimates between METY and METINT with corrected milk yield for fat, protein, and or energy (CMY) were negative (− 0.433 and − 0.262, respectively). Also, genetic correlation estimates between METINT with milk fat and protein percentages were 0.254 and 0.334, respectively. Although the genetic correlation estimate of METP with daily milk yield was 0.172, its genetic correlation with CMY was 0.446. All genetic correlation estimates between METP with milk fat and protein yield or percentage ranged from 0.005 (between METP-milk protein yield) to 0.185 (between METP-milk protein percentage). The current meta-analysis confirmed the presence of additive genetic variation for methane emission traits in dairy cows that could be exploited in genetic selection plans.

estimates and genetic correlations. The logic for planning the present study was the requirement for integrating estimates from former studies to prepare summary genetic parameter estimates for methane emission traits to help establish breeding goals for dairy cattle. In addition, to prevent inappropriate correlated responses with other economically important traits, recognizing the genetic associations with production traits is essential before implementing methane into the breeding objectives of dairy cows. The present study aimed to perform a meta-analysis using the three-level model to integrate published estimates of genetic parameters for methane emission traits in dairy cows.

Materials and methods
Characterizing the scope of the meta-analysis study. The PRISMA (Preferred Reporting Items for Systematic Reviews and Meta-Analyses) guideline to address meta-analyses and systematic reviews was used in this study (Fig. 1) 15 . A systematic literature search using electronic databases of ISI Web of Knowledge (https:// apps. webof nowl edge. com), Google Scholar (https:// schol ar. google. com), and NCBI (https:// www. ncbi. nlm. nih. gov) was performed to recognize all citations reporting heritability estimates for methane emission traits and their genetic association with milk yield and composition of dairy cows. A comprehensive search was conducted with the following keywords and their synonyms or derivatives: "dairy cow", "methane emission", "enteric methane", "methane production", "methane intensity", "methane yield", "genetic parameters", "heritability", and "genetic correlation". Overall, 40 heritability estimates and 32 genetic correlations from 17 papers were used in this meta-analysis study. The included papers were published between 2015 and 2021. The characteristics of studies included in the database for conducting this meta-analysis are indicated in Table 1. The genetic parameter estimates were obtained from mixed animal models based on the Bayesian inference and restricted maximum likelihood (REML) estimation methodologies. Thus, parameter estimates from reduced models were  www.nature.com/scientificreports/ excluded. Only papers published in scientific index journals were used, and papers published in other sources were removed. The references in the above papers were also controlled. Methane emission traits considered in this study were methane yield per kg dry matter intake (METY, in g/kg), methane intensity per kg fat and protein corrected milk produced (METINT, in g/kg), and methane production as daily methane production per cow (METP, in g/day).
Data recorded. The data sets included information on the estimates of heritability for METY, METINT, and METP, and also genetic associations of methane emission traits with production traits [daily milk yield (DMY), corrected milk yield for fat, protein, and or energy (CMY), milk fat percentage (Fatp), milk protein percentage (Prop), milk fat yield (Faty), and milk protein yield (Proy)], and their standard errors. Other details registered were the year of publication, journal name, the number of observations or records, phenotypic mean and standard deviation, breed name, country of origin, parity number, years of data collection, univariate or multivariate analysis model, and the estimation method. When the same estimate was reported in different publications, based on the same database, only the most recent publication was included in the analysis. The analysis was performed exclusively for traits in which the parameter estimates were placed on not less than two distinct data sets. For articles in which the standard errors for the heritability or correlation estimates were not reported, approximated standard errors were derived by using the combined-variance method 30 , which is given by the following formula: where SE ij is the predicted standard error for the published parameter estimate for the ith trait in the jth article that has not reported the standard error, s ik is the published standard error for the parameter estimate for the ith trait in the kth article that has reported the standard error, n ik is the number of used records to predict the www.nature.com/scientificreports/ published parameter estimate for the ith trait in the kth article that has reported the standard error, and n´i j is the number of used records to predict the published parameter estimate for the ith trait in the jth article that has not reported the standard error.

Meta-analysis of heritability estimates and genetic correlations.
The database used in the current study has a hierarchical structure because some effect sizes were extracted from studies conducted by the same authors (Table 1). Indeed, the true underlying effects are expected to be more similar for such studies (i.e., effect sizes are likely correlated). A three-level meta-analysis model accounted for this dependency among studies 31 .
The random effects of the study and effect size were specified as a list of one-sided formulas in the random argument of the rma.mv function of the metafor package version 3.0-2 32 in R software. The REML method was used and the effect sizes with the same level within each grouping variable received the same random effect; otherwise, effect sizes were assumed to be independent. Forest plots were built to demonstrate the effect size for each study. Effect sizes in forest plots were the average estimates of heritability for methane emission traits or their genetic association with milk production traits with a 95% confidence interval.
Because of the non-normal nature of correlation estimates, almost all meta-analyses do not apply the published correlation estimates. Preferably, the reported correlation estimate is transformed to the Fisher's Z scale, and every analysis is conducted with the converted values 33 . Estimated parameters and their confidence intervals would then be transformed into correlations. The Fisher's Z transformation to obtain an approximate normal scale is represented as follows 33 : where r gij = the reported estimate of genetic correlation for the ith trait in the jth paper. The following formula was utilized to return to the original scale: where r * g ij = the re-transformed genetic correlation for the ith trait in the jth article and Z ij = the Fisher's Z transformation.
The 95% lower and upper limits for the estimated parameter would be calculated respectively for each trait as follows: where SE θ = the predicted standard error for the estimated parameter θ , given by: Heterogeneity. The I 2 statistic was utilized to assess heterogeneity as follows 34 : where Q = the χ 2 heterogeneity statistic and k is the number of studies. Q is the Q statistics given by the following formula: where w j = the parameter estimate weight [assumed as the inverse of published sampling variance for the param- ] in the jth article; θ j = the published estimate of the genetic parameter in the jth paper, θ = the weighted mean of the parameter in the population,, and k is the number of used papers. The I 2 statistic characterizes the variation among the studies because of heterogeneity. Negative values of I 2 are considered zero; therefore, I 2 values vary between 0 and 100% 35 . If the values of I 2 fall within the range of 0-40%, there is no concern about the heterogeneity. However, the I 2 values of 40-60% and 60-100% often represent moderate and substantial heterogeneities, respectively.
Publication bias. Egger's linear regression asymmetry was utilized to test the existence of publication bias.
Also, to demonstrate asymmetry, funnel plots were applied. When no publication bias exists, the funnel plot presents the symmetric distribution of effect sizes around the actual effect size.

Results
Descriptive statistics. The number of literature estimates, measurement units, the total number of records, weighted mean, standard deviation, and the coefficient of variation for methane emission traits and milk yield and composition of dairy cows are indicated in Table 2. The weighted coefficients of variation for methane emission traits were generally low to moderate and varied from 4.60 (for METINT) to 23.64% (for METP). In addition, the weighted coefficient of variation for milk yield and composition was low and varied from 7.14 (for FY) to 15.47% (for CMY).
Heritability estimates. Effect size and heterogeneity of the heritability estimates for methane emission traits obtained from the three-level model of the meta-analysis are presented in Table 3. The heritability estimates for METY, METINT, and METP were 0.244, 0.180, and 0.211, respectively. These estimates had low standard errors, and their 95% confidence intervals were small. Also, the heritability estimates for methane emission traits were significant (P < 0.05). The I 2 values showed minor heterogeneity for the heritability estimates of methane emission traits ( Table 3). The results of Egger's test for the occurrence of possible publication bias showed significant publication bias for METY (P = 0.088) and METP (P = 0.081), but non-significant publication bias (P = 0.125) was observed for METINT. The forest plots of individual studies for heritability estimates of METY, METINT, and METP in dairy cows are indicated in Fig. 2. Also, funnel plots of heritability estimates for METP and METINT are shown in Supplementary Figures S1 and S2, respectively.

Genetic correlation estimates. Effect size and heterogeneity of the genetic correlation estimates between
methane emission traits with milk yield and composition in dairy cows obtained from the three-level model of the meta-analysis are shown in Table 4. The genetic correlation estimates between METY-CMY and METINT-CMY were negative (− 0.433 and − 0.262, respectively). Genetic correlation estimates between METINT with Fatp and Prop were positive and moderate (Table 4). Although the genetic correlation estimate of METP with DMY was 0.172, its genetic correlation with CMY was 0.446. All genetic correlation estimates between METP with milk fat and protein yield or percentage ranged from 0.005 (between METP-milk protein yield) to 0.185 (between METP-milk protein percentage) ( Table 4). Except for genetic correlations between METP-CMY and METP-Faty, other genetic correlations were non-significant (P > 0.05). Hence, the 95% confidence interval of non-significant genetic correlation estimates included zero. The I 2 values indicated substantial heterogeneities for the genetic correlations between METINT-CMY and METP-DMY (Table 4). In addition, high heterogeneity was observed for the genetic correlation estimate between METINT-Prop. Also, moderate heterogeneity was observed for the genetic correlation between METP-CMY. Other genetic correlation estimates had negligible heterogeneities ( Table 4). The forest plots of individual studies for genetic correlation estimates between METP-CMY and METP-Faty are depicted in Fig. 3. Also, the forest plots of individual studies for other genetic correlation estimates are displayed in Supplementary Figures S3 to S12. The results of Egger's test showed significant  Table 3. Effect size and heterogeneity of the heritability estimates for methane emission traits of dairy cows obtained from three-level model of meta-analysis. *For traits, see Table 2. www.nature.com/scientificreports/ (P = 0.000) publication bias for the genetic correlation between METP and Fatp, but non-significant (P > 0.10) publication bias was observed for other genetic correlation estimates.

Discussion
This study is the first comprehensive meta-analysis of the heritabilities for methane emission traits and their genetic correlation with milk production traits in dairy cows. Only three methane emission traits of METY, METINT, and METP were considered in this study, and residual methane production was not considered because www.nature.com/scientificreports/ the number of genetic parameter estimates reported for this trait was lower than the required meta-analysis standards.
Knowledge of heritability and genetic correlations with other traits of economic importance is essential to include methane emission traits in breeding goals 36 . Present breeding goals for dairy cattle do not incorporate enteric methane traits. However, the genetic improvement of farm animals is especially a constructive methodology, providing cumulative and constant modifications in the traits of interest 36 . Therefore, determining the genetic variability of methane emission traits would be the first step in including these environment-related traits in the proposed selection indices 19 . In this context, three preconditions are required. First, methane emission traits must be adequately heritable to permit a somewhat quick and meaningful improvement. Second, enough genetic variation for these traits must be proven in the studied dairy cow population. Third, genetic associations of methane emission traits with other traits of interest must be understood. The genetic analysis must confirm these three preconditions 19 . Only the relationships with milk production traits were considered in this paper.
The generally low weighted coefficients of variation for the studied traits indicated the lower dispersion around the weighted phenotypic means. This result implied the weighted phenotypic means for these traits were accurate. The low weighted coefficients of variation for METINT and METY showed that the phenotypic variation Table 4. Effect size and heterogeneity of the genetic correlation estimates between methane emission traits with production traits in dairy cows obtained from the three-level model of meta-analysis. *For traits, see Table 2. r g : Genetic correlation.  www.nature.com/scientificreports/ for these traits is restricted biologically. However, the highest weighted coefficient of variation was observed for METP. This result indicates that higher phenotypic variation existed in this trait than in other traits. The current meta-analysis confirmed the presence of a genetic component for methane emission traits in dairy cows. The low heritability estimate for METINT indicated the minor influence of additive genetic effects on this trait. However, the average heritability estimates for METY and METP showed a medium effect of additive genes on these traits and a possibly suitable response to selection for them. Numerous variables may affect methane production, including measurement time, herd, diet composition, season, stage of lactation, and more. Therefore, it is necessary to adjust the relevant models of analysis for these effects when daily or lactational methane production is predicted. This adjustment will prevent bias in the phenotypic records, genetic predictions, and estimation of variance components 37 . Because meta-analysis combines published genetic parameter estimates reported by different studies, it is anticipated that the actual parameter could differ among the studies.
In general, most genetic correlation estimates had a wide 95% confidence interval which included zero. Therefore, these correlations must be interpreted with caution. The positive genetic correlation between METP and CMY likely reflects the association between methane production, energy intake, and milk yield 13 . This result would be expected because the increase in the genetic potential of animals to produce more milk increases methane emissions per animal because of an increase in feed consumption 38 . The positive genetic correlations between METP with CMY and Faty indicate that genetic and physiological mechanisms controlling these traits could be similar. It was reported that genes responsible for methane production also control lipid synthesis 39 . Milk composition and fatty acids are often used to predict methane production [40][41][42] . The positive genetic correlation suggests that breeding for increased milk production with higher fat content can increase methane production. The negative genetic correlations of METY and METINT with CMY indicated the opposite direction of changes when genetic selection is directly performed on milk production. Because of the negative genetic correlation between METINT-CMY, the positive genetic correlations of METINT with Fatp and Prop would be expected. The possible reasons for these positive genetic correlations would be the negative correlation between milk yield with Fatp and Protp and the dilution effect. It seems that selection for higher milk yield, and accordingly higher feed consumption and likely live weight, results in a rise in daily methane production per cow. If METP is included with a negative weight within a selection index, it would be possible to expect a decrease in the genetic gain in milk yield due to the unfavorable positive genetic correlation between these two traits.
Defining methane emission traits as ratios is a helpful metric for characterizing groups of animals, such as different treatment groups, herds, breeds, and species. However, ratio traits usually violate two statistical assumptions, which can affect the definition of the linear relationship (regression or correlation) between the two sets of traits, making them unsuitable for incorporation in genetic selection programs 43,44 . First, it is assumed that a ratio is independent, or uncorrelated, to its numerator or denominator. Second, it is assumed that the relationship between a ratio and its component traits is linear. Therefore, applying methane ratio traits in animal breeding and genetics is challenging because the complicated statistical characteristics of these ratio traits may cause unfavorable correlations (i.e., METINT may be unfavorably associated with energy-corrected milk, which is the denominator of the METINT ratio) 37 .
The meta-analysis of some genetic parameter estimates showed significant publication bias and or asymmetry in the funnel plots. Asymmetrical funnel plots may indicate publication bias or be due to exaggeration of treatment effects in small studies of low quality. Egger et al. 45 stated that the potential sources of asymmetry in funnel plots are generally grouped as follows: 1. Selection biases, 2. True heterogeneity (effect size differs according to study size), 3. Data irregularities, 4. Artefact (heterogeneity due to poor choice of effect measure), and 5. Chance.

Conclusion
The current meta-analysis confirmed the presence of a genetic component for methane emission traits in dairy cows that could be exploited in genetic selection plans. The positive genetic correlations between METP with CMY and Faty mean that a decrease in methane production should have adverse effects on milk and fat yields. Therefore, cows producing higher milk and milk fat are expected to emit more methane. The positive genetic correlations between METP with CMY and Faty indicate that genetic and physiological mechanisms controlling these traits could be similar.

Data availability
All data generated or analysed during this study are included in this published article (and its Supplementary  Information files).