The descriptive statistics for studied traits is presented in Table 1. The results of variance analysis showed that the contemporary group effect was only statistically significant for PME, butyric acid and isovaleric acid (p < 0.05), while the effect of age was only statistically significant for PME.
Genome-wide association studies for volatile fatty acids
After removing the genotypes with low imputation accuracy (R2 < 0.30) and other criteria (e.g. call rate < 95%, MAF < 0.05, and Chi-square < 10-6 of Hardy-Weinberg equilibrium), a total of 6,583,595 SNPs on 29 Bos taurus autosomes (BTAs) were remained. The number of quality-control passed imputed variants used for association analyses ranged from 125,756 on BTA25 to 443,774 on BTA1.
The results demonstrated that 76 (on BTA3, 6, and 17), 3 (on BTA12), 5 (on BTA3 and 17), 274 (on BTA2, 4, 5, 9, 10, 11, 12, 22, 24, 25, 26 and 27) and 285 (BTA2, 3, 4, 9, 10, 11, 12 and 28) SNPs passed the suggestive-association threshold (p < 10−5) for acetic acid (%), butyric acid (%), propionic acid (%), isovaleric acid (%) and valeric acid (%) traits, respectively (Figure 1). Alemu et al. [11] reported a strong correlation (-0.85) between acetic and propionic acid. Our study also showed that there were some regions (e.g. BTA3) associated with both acetic and propionic acid traits. Further, using the chromosome-wide Bonferroni correction threshold, we found two SNPs significantly associated with isovaleric acid (BTA9 and 28) and 29 SNPs with valeric acid (BTA5, 11, and 25). However, we did not discover any significant association for other VFA traits based on chromosome-wide Bonferroni correction. Moreover, using the recommended significant threshold of 5 × 10−8, 18 SNPs (two on BTA5 and 16 on BTA 25) were associated with valeric acid variation, two of which on BTA5 were found to be significant at the genome-wide Bonferroni correction threshold of 7.59 × 10-9 (Table 2).
The isoacids (e.g. isobutyric, 2-methylbutyric, isovaleric acid, and straight-chain valeric acid) produced in the digestion process of ruminants are mainly generated from the degradation products of the amino acids valine, isoleucine, leucine, and proline [12]. Isoacids are either required by, or stimulate the growth of ruminal cellulolytic bacteria. Therefore, they play a critical role in microbial fermentation [12]. According to the previous studies, milk production in dairy cattle could increase by 5 to 10% with the addition of isoacids to the diet [13]. Therefore, the existence of significant association between host genome and both iso-valeric acid, and valeric acid traits can help the dairy farmers to improve these acids availability in rumen for the ruminal cellulolytic bacteria using genetic selection.
Genome-wide association studies for predicted methane emission
The GWAS results for PME, PME per kg milk, and PME per kg fat are presented in Manhattan plots in Figure 2. According to suggestive regions (p < 10-5), we found 53 SNPs (BTA8, 10, 15, 18, 23, and 26) for PME; 41 SNPs (BTA15) for PME per kg milk; and 308 SNPs (BTA1, 2, 4, 5, 8, 10, 13, 14, 16, 19, 20, 21, 26, and 28) for PME per kg fat. Although, the association studies did not show any significance SNPs for PME trait, the chromosome-wide Bonferroni correction identified five significant SNPs on BTA15 for PME per kg milk, 20 significant SNPs on BTA4, 13, 19, and 21 for PME per kg fat. PME per animal that were not adjusted per unit of product showed a deceased variation compared to the adjusted PME per kg of milk or fat i.e. two animals with an equal rate of methane emission do not necessarily have an equal efficiency for methane production (methane emission per unit of product). As a result, the power of the association studies to capture quantitative trait loci (QTL) without considering the adjustment of PME per unit of product can be diminished. Using the recommended Bonferroni correction threshold (p < 5 × 10−8), we identified two SNPs on BTA15 for PME per kg milk, 14 SNPs on BTA4, 13, 19, and 28 for PME per kg fat (Table 2); and two SNPs on BTA4 and 19 passed the genome-wide Bonferroni correction for PME per kg fat (Table 2).
Given that animals in this study were sampled based on highest and lowest percentiles of milk production EBV with an equal number of animals, this strategy can enhance the ability of QTL detection especially in population with small sample size. However, it has been reported that Bonferroni correction can increase the number of type II errors, particularly when the sample size is small [14]. Thus, chromosome-wide Bonferroni correction can be considered more accurate than genome-wide Bonferroni correction [15]. Since, the sample size in our GWAS was relatively small, the power of finding SNPs that are associated with PME traits was limited, and thus a few number of SNPs passed the genome-wide Bonferroni correction stringent p-value threshold.
Cattle have a considerable ability to produce high-quality proteins from non-edible plant cell wall components for human consumption [16]. However, it is well established that gastrointestinal microbiota contributes to feed digestion and enteric methane emission in ruminants [17-20], as approximately 17% of global methane emissions generated through ruminants [8]. Methane emission is a difficult and expensive trait to measure and is not routinely measured in dairy cattle. Furthermore, this trait has not been directly considered in the selection index of dairy cattle. However, the detrimental environmental effects of cattle production can be reduced significantly by improving the efficiency of cattle in converting feed to milk and meat, rather than wasting energy as enteric methane. To improve this efficiency in dairy cattle, it is crucial to understand the sources of genetic variation in methane production among individual animals, and the genetic architecture of methane production. van Engelen et al. [21] reported that the heritability of predicted methane yields ranges from 0.12 to 0.44 in Dutch Holstein-Friesian cows, and suggested that PME based on milk fatty acids can be reduced using genetic selection programs. de Haas et al. [7] estimated the heritability of PME and PME per fat- and protein-corrected milk yield to be 0.35 and 0.58, respectively. However, they reported that seven SNPs associated with PME in Holstein cattle could explain only 0.2% of the total genetic variance [7]. In another GWAS, 3,304 significant SNPs were identified for PME (p < 0.005) in which 630 of them were also associated with weight-at-test and dry-matter intake [8]. It has been estimated that 19% to 23% of methane emissions per kg of milk can be reduced and converted to production when cows are selected based on milk fat plus protein (19%), or based on the average genetic merit of milk fat plus protein yield (23%) [22]. According to Yan et al. [23], an effective method for mitigating methane emissions in dairy cows is to select the animals according to increased milk-production efficiency and increased energy-utilization performance.
Understanding the genetic variation and underlying genes associated with PME allows reduction of methane production in cattle through marker-assisted or genomic selection. It has been reported that genomic selection can reduce PME up to approximately 5% over 10 years [24]. The results of our study showed a large variation of PME traits that were used for mitigation of methane emissions in dairy cattle might. Thus, by including the QTL associated with the PME traits in breeding programs, the rate of methane emission per animal is expected to be reduced. Moreover, reduction in methane emissions can also improve feed-conversion efficiency in dairy cattle. In addition, unlike other strategies developed to reduce methane emissions in dairy cattle such as change in dietary formulations, feed additives, and anti-methanogen vaccines, genetic improvements have the benefits of being heritable, cumulative, and permanent.
Post-GWAS bioinformatics analysis
The candidate genes were identified within 1 Mb of most significant SNPs based on p < 5 × 10−8. For valeric acid, two and 16 candidate genes were discovered around 5:16795260 and 25:37967076 SNPs, respectively (Table 3). For PME per kg milk, four candidate genes around 15:25797132 SNP, and for PME per kg fat 11, 3, 32 and 1 candidate genes were discovered around 4:115131249, 13:81673732, 19:24494923 and 28:21771233 SNPs, respectively (Table 4). Ten genes of 17 candidate genes for valeric acid remained for Gene Ontology (GO) term analysis as known and non-ambiguous genes. Based on the GO terms, 6 promising candidate genes were significantly clustered in organelle organization (GO:0004984, p = 3.9 × 10-2) including BAIAP2L1, SMURF1, ARPC1A, ARPC1B, BHLHA15 and TRRAP. Two of the most promising candidate genes (ARPC1A and TRRAP) were related to residual feed intake (RFI) [25, 26]. It has been observed that low-RFI pigs had a tendency toward lower valeric acid (p = 0.07) and isovaleric acid (p = 0.09) concentrations [27]. Therefore, available valeric acid concentration might be influenced by these genes in cattle. Twenty-five genes of 40 candidate genes for methane emission trait remained for GO term analysis as known and non-ambiguous genes. Based on the GO terms, 17 candidate genes were significantly clustered in olfactory receptors activity (GO:0004984, p = 4 × 10-10) including LOC538966, LOC522582, LOC540082, LOC618593, LOC508980, LOC509525, LOC509526, LOC617122, LOC532238, OR1E1, LOC618124, OR1G1, LOC511509, LOC526294, LOC615901, LOC618112 and LOC101902679. Olfactory receptors genes influenced feeding behavior such as food preferences and feed intake [28]. Furthermore, they are associated with RFI and modulated olfactory transduction pathways in pig [29]. Since there was a positive genetic correlation (ranging from 0.18 to 0.84) between RFI and PME in Holstein-Friesian cows [7], olfactory receptors genes may affect methane mission per animal. In addition, five candidate genes including CYP51A1 on BTA 4, PPP1R16B on BTA 13, and NTHL1, TSC2, and PKD1 on BTA 25 suggested by Pszczola et al. [30], involved in a number of metabolic processes that might be related to methane emission. PKD1 gene was associated with development of the digestive tract [30].
The results of gene networks analyses for valeric acid and methane emission traits were shown in Figures 3 and 4, respectively. There were some contributions among candidate genes by co-expression, pathway, physical interactions and shared protein domains for valeric acid and methane emission traits. Most of these contributions were related to physical interactions between genes that was indicative of protein-protein interaction and if two genes showed the same protein-protein interaction, their products were linked together. Therefore, identified candidate genes in our study had a significant protein-protein interaction either together or with associated genes.
The summary of significant SNPs (p < 5 × 10−8) associated with PME and valeric acid traits that are in close distance to reported QTLs is presented in Table 5. Given that methane emission is a complex trait influenced by many genes, and is not measured routinely in dairy cattle (unlike the traits of milk production and body weight), there are very limited number of QTLs associated with methane emission in the QTLdb website. Most of the QTLs that are in close distance to SNPs associated with PME traits have been found to influence milk production and its components, body weight, and RFI, rather than methane production. However, all SNPs associated with valeric acid in our study were close to the reported QTLs for milk production and its components, and RFI. The correlation coefficients between EBV for methane intensity with milk yield, fat yield, protein yield, and somatic cell score were estimated as -0.68, -0.13, -0.47, and 0.07, respectively [31]. Further, a positive genetic correlation between RFI and PME (ranging from 0.18 to 0.84) was reported in Holstein-Friesian cows [7]. Herd et al. [32] also reported a high correlation between daily methane production and yearling weight (0.85), and dry-matter intake (0.79) in beef cattle. These reports indicated that some common genes and QTLs influenced traits related to methane-emission, feed-efficiency and milk-production traits that was confirmed by our study. However, further studies are needed to validate the results of our GWAS and investigate the biological effects of validated SNPs.