Genome-wide association study for intramuscular fat deposition and composition in Nellore cattle

Abstract

Background

Meat from Bos taurus and Bos indicus breeds are an important source of nutrients for humans and intramuscular fat (IMF) influences its flavor, nutritional value and impacts human health. Human consumption of fat that contains high levels of monounsaturated fatty acids (MUFA) can reduce the concentration of undesirable cholesterol (LDL) in circulating blood. Different feeding practices and genetic variation within and between breeds influences the amount of IMF and fatty acid (FA) composition in meat. However, it is difficult and costly to determine fatty acid composition, which has precluded beef cattle breeding programs from selecting for a healthier fatty acid profile. In this study, we employed a high-density single nucleotide polymorphism (SNP) chip to genotype 386 Nellore steers, a Bos indicus breed and, a Bayesian approach to identify genomic regions and putative candidate genes that could be involved with deposition and composition of IMF.

Results

Twenty-three genomic regions (1-Mb SNP windows) associated with IMF deposition and FA composition that each explain ≥ 1% of the genetic variance were identified on chromosomes 2, 3, 6, 7, 8, 9, 10, 11, 12, 17, 26 and 27. Many of these regions were not previously detected in other breeds. The genes present in these regions were identified and some can help explain the genetic basis of deposition and composition of fat in cattle.

Conclusions

The genomic regions and genes identified contribute to a better understanding of the genetic control of fatty acid deposition and can lead to DNA-based selection strategies to improve meat quality for human consumption.

Background

Many consumers associate consumption of fat from beef with coronary heart disease, diabetes and obesity, due to the presence of cholesterol, high concentration of saturated fatty acids (SFA), and low concentration of polyunsaturated fatty acids (PUFA). However, consumption of fatty acids is necessary for human nutrition [1]. Beef has high nutritional value from children to seniors, is a rich source of protein (essential amino acids), iron, zinc, B vitamins and essential polyunsaturated fatty acids such as linoleic and linolenic acid [2]. Beef fat also has a high concentration of monounsaturated fatty acids (MUFA), whose melting point is low and can reduce the concentration of bad cholesterol (LDL) in blood circulation [3].

The amount of fatty acid and its composition in beef varies by breed, nutrition, sex, age and carcass finishing level [4]. The difficulties associated with determining intramuscular fat (IMF) deposition and composition as well as the limited knowledge on the genetic mechanisms that control these traits has limited genetic progress in the production of healthier beef.

The development of high-density bovine genotyping [5] and their use in genome-wide association studies (GWAS) have allowed identification of genomic regions associated with phenotypes of interest. The technique of GWAS exploits differences in allele frequencies of thousands of polymorphic markers available in unrelated individuals who possess different phenotypes (for example, deposition and composition of intramuscular fat), and leads to the identification of markers associated with a given phenotype [6]. Bayesian approaches have been applied to GWAS to detect significant quantitative trait loci (QTL) for traits of economic importance. One such approach uses multiple regression (evaluating marker effects simultaneously), treating marker effects as random to reduce overestimation bias of significant QTL effects, generating the actual posterior distribution of QTL effects given the data which can provide richer inference than can be obtained by simply constructing p-values as well as providing an alternative to the use of p-values to avoid false positives [7–9].

Brazilian beef is exported and consumed in more than 100 countries [10]. Purebred and crossbred Nellore cattle, which are of Bos indicus descent, are the predominant source of beef in Brazil. Previous research has documented that muscle and fat tissues from Bos indicus cattle develop in a different manner than in Bos taurus breeds [11–14]. However, studies documenting the genetics of fatty acid deposition and composition in Bos indicus breeds are limited.

In this study we performed a GWAS using high-density single nucleotide polymorphism (SNP) chips (770 k) and Bayesian methods (Bayes B) to identify genomic regions associated with fat deposition and fatty acid composition (FA) in Nellore beef.

Results and discussion

Intramuscular fat deposition and composition

Modern consumers are concerned with their overall health and often desire to reduce their caloric intake. This has increased the demand for lean meat production with a healthier fatty acid composition, which would comprise a lower proportion of SFA and greater proportion of MUFA. The amount of fat deposition in meat represented as IMF (mean = 2.77%) reported in this feedlot-finished study was greater than for pasture-finished counterparts, as expected, but lower than normally observed in Continental and English breeds [15–17]. Despite of that, IMF observed in this work was within limits of reasonable amount of fat to assure acceptable quality levels for consumers according to Nuernberg et al. [18].

The fatty acid composition observed for the most abundant FAs were: C14:0 at 3.54%, C16:0 at 26.69%, C18:0 at 14.98%, C16:1 cis-9 at 3.31%, C18:0 at 14.98%, C18:1 cis-9 at 37.46%, C18:2 cis-9 cis 12 at 1.60%, SFA at 47.23%, MUFA at 48.34%, and PUFA at 2.87% (Table 1). The proportion of MUFA was higher than SFA, and oleic acid (C18:1 cis-9) was the most abundant single fatty acid (37.46%).The FA composition presented in this work is similar to those reported in the literature for Nellore or other Bos indicus breeds [19–21]. This population also presented, in relation to Bos taurus breeds, average composition of fatty acids which is in agreement with reports from the USDA [22]. However, in this study a lower quantity of PUFA was observed in general, but not for C18:2 cis-9 trans-11, consequently a lower ratio of PUFA/SFA (6.08%). Similar MUFA and PUFA results have been reported in previous studies that utilized Bos indicus steers [20, 23].

Table 1 Descriptive statistics, variance components and heritability by GBLUP for IMF deposition and composition in Nellore

The omega 6: omega 3 (n-6:n-3) ratio was high (4.84), which has a beneficial effect on patients with chronic diseases. A ratio of 4:1 has been associated with a 70% decrease in total mortality in the secondary prevention of cardiovascular disease [24]. Omega 6 and 3 fatty acids are essential cell membrane and tissue constituents that are important for many biological functions. These fatty acids are also essential for synthesis of prostaglandins, thromboxane, leukotriene, hydroxyl fatty acids, and lipoxins that are involved with inflammatory response. Humans and other mammals can convert n-6 to n-3 using desaturation enzymes but this conversion is slow and there is competition between n-6 and n-3 fatty acids for the desaturation enzymes [25].

Heritability

Descriptive statistics and heritabilities estimated using a genomic relationship G matrix are in Table 1. Heritabilities estimated in this study varied from low (<0.10 for C12:0, C16:0, C18:1 cis-11, C18:1 cis-12, C18:2 cis-9 trans-11, C20:1, C20:3 n-6, C20:4 n-6, C20:5 n-3 and AI, respectively) to moderate (up to 0.29 for IMF, C14:0, C14:1 cis-9, C16:1 cis-9, C17:0, C17:1, C18:0, C18:1 cis-9, C18:1 trans-6, 7, 8, C18:1 trans-10, 11, 12, C18:2 cis 9 cis-12 n-6, C18:2 trans-11 cis-15, C18:3 n-6, C18:3 n-3, C22:5 n-3, C22:6 n-3, SFA, MUFA, PUFA, Sn-3, Sn-6 and n-6:n-3). For C15:0, C18:1 cis-13, C18:1 cis-15, C18:1 trans-16, C20:2 the heritability estimates were zero. In Angus [26] and Japanese Black cattle [27] estimates of heritability for IMF fat deposition and composition traits were higher than in this study. The lower values of heritability reported for this population could be explained by the reduced sample size [28] or lower amount of genetic variation in the population [29].

The estimates of genomic heritability from Bayes B approach (Table 2) were also low to moderate (from 0.09 to 0.46) and generally similar to values obtained in ASReml using a genomic relationship matrix.

Table 2 Posterior means of variance components for IMF deposition and composition in Nellore by Bayes B

Genome wide association studies and genomic regions identified

GWAS and genomic regions identified were reported for those traits that showed genomic heritability ≥ 0.10: IMF, 33 different FA and ratios, and one FA index. The 1-Mb SNP window regions that explained more than 1% of the genetic variance were used to search for putative candidate genes (PCG) and are presented on Table 3. Therefore, 35 traits were used for GWAS and these represented 23 different 1-Mb genomic regions (Table 3). These regions were distributed over 12 different chromosomes: 2, 3, 6, 7, 8, 9, 10, 11, 12, 17, 26 and 27 and the corresponding PCG in these regions are reported in Table 3. Intramuscular fat was one of the traits that presented moderate heritability (0.25), however no region that explained more than 1% of the genetic variance for this trait was identified.

Table 3 QTL regions associated with fatty acid composition in Nellore steers by Bayes B

Saturated fatty acids

Eight genomic regions (1 Mb windows) explained more than 1% of genotypic variation for C12:0, C14:0, C16:0, and C18:0 (Table 3). These regions overlap with QTL previously reported for marbling score [28], backfat thickness [29], carcass weight and body weight in Angus cattle [30]. Additional file 1 shows in more detail the Manhattan plot of the proportion of genetic variance explained by window across the 29 autosomes for the important saturated fatty acids for beef palatability: C12:0 (lauric acid), palmitic acid (C16:0), stearic acid (C18:0).

BTA12 at 60 Mb QTL region was associated with C12:0 fatty acid and harbors the PCG SLIT and NTRK-like family, member 6 (SLITRK6). This gene is related to the lipopolysaccharide receptor complex (GO:0016021). Lipopolysaccharide (endotoxin) is an important component of gram-negative cell walls, which can start the inflammatory process by binding to the SLITRK6 complex. This complex is located at the surface of innate immune cell, which presents an important role in obesity associated with insulin resistance [31].

BTA3 at 6 Mb and BTA9 at 36 Mb were associated with C14:0 fatty acid. In the first region, no PCG was identified. In the second region, the following PCGs were identified: guanine nucleotide binding protein (G protein), gamma 11 (GNG11) and regulator of G-protein signaling 5 (RGS5). These genes are associated with G proteins that have been implicated in the regulation of body weight and metabolic function, hyperinsulinemia, impaired glucose tolerance and resistance to insulin in mice [32].

BTA3 at 26 Mb and BTA12 at 13 Mb were associated with C16:0 fatty acid. In the first QTL region the tryptophanyl tRNA synthetase 2, mitochondrial (WARS2) gene was identified. In BTA12 at 13 Mb no PCG was identified. WARS2 gene encodes an essential enzyme that catalyzes aminoacylation of tRNA with tryptophan. The WARS2 protein contains a signal peptide for mitochondrial import, OMIM: 604733 [33], and is involved with regulation of fat distribution in human visceral and subcutaneous fat [34].

BTA3 at 25 Mb, BTA10 at 74 Mb and BTA11 at 105 Mb QTL regions were associated with C18:0 fatty acid. On BTA3 at 25 Mb four PCG were identified: 3-hydroxy-3-methylglutaryl-Coenzyme A synthase 2, mitochondrial (HMGCS2), phosphoglycerate dehydrogenase (PHGDH), hydroxy-delta-5-steroid dehydrogenase, 3 beta- and steroid delta-isomerase 1 (HSD3B1) and hydroxyacid oxidase 2, long chain (HAO2). HMGCS2 was associated with fatty acid oxidation and ketogenesis in HepG2 cells [35] and is directly regulated by PPAR-α gene that is an important key regulator of β-oxidation. On the other hand the PPAR-γ expression induces lipogenesis by PXR activation in mice liver [36]. PHGDH is related to the regulation of gene expression according to gene ontology (GO) terms (GO:0010468), HSD3B1 gene is associated with steroid biosynthesis (GO:0006694) and metabolic process (GO: 0008202) and HAO2 is associated with fatty acid oxidation (GO:0019395).

On BTA10 at 74 Mb QTL region, dehydrogenases/reductases (SDRs) (DHRS7) was identified. DHRS7 catalyzes the oxidation/reduction of a wide range of substrates, including retinoid and steroids [37] and has high expression level in adipocyte and skeletal muscle [38]. In addition, this gene is responsible for the final step in cholesterol production, the conversion of 7-dehydrocholesterol to cholesterol [39].

On BTA11 at 105 Mb region, NUP214 was identified. This gene is involved with intermembrane transport. The nucleoporins are characterized by phenylalanine-glycine rich (FG) repeat sequences and the FG domains have an unfolded structure and are responsible for interaction with importin-cargo complexes that move through the pore. One of these importin-cargo complexes is sterol regulatory element binding proteins-sterol-sensing accessory factor (SREBP-SCAP). This factor enters the nucleus, then binds to sterol regulatory elements (SRE) in the promoter regions of genes, whose products mediate the synthesis of cholesterol and fatty acids [40].

Monounsaturated fatty acids

Ten genomic regions (1 Mb region) explained more than 1% of genotypic variation for monounsaturated fatty acids, which relates C14:1 cis-9, C16:1 cis-9, C18:1 cis-9, and C18:1 trans-16 (Table 3). These regions overlap with QTL reported for carcass weight, marbling score in Angus [30], docosahexaenoic acid content in Charolais x Holstein crossbred cattle [41] and palmitoleic acid in dairy cattle [42]. Manhattan plots of the proportion of genetic variance explained by each 1-Mb window (2,527 windows) across the 29 autosomes for the most important fatty acid for beef quality and human health: myristoleic acid (C14:1 cis-9), palmitoleic acid (C16:1 cis-9), and oleic acid (C18:1 cis-9) and are in Figures 1, 2 and 3, respectively.

Figure 1

Manhattan plot of the genome-wide association study for C14:1 cis 9 (myristoleic acid) in Nellore. The X-axis represents the chromosomes, and the Y-axis shows the proportion of genetic variance explained by SNP window from Bayes B analysis.

Figure 2

Manhattan plot of the genome-wide association study result for C16:1 cis9 (palmitoleic acid) in Nellore. The X-axis represents the chromosomes, and the Y-axis shows the proportion of genetic variance explained by SNP window from Bayes B analysis.

Figure 3

Manhattan plot of the genome-wide association study result for C18:1 cis9 (oleic acid) in Nellore. The X-axis represents the chromosomes, and the Y-axis shows the proportion of genetic variance explained by SNP window from Bayes B analysis.

BTA2 at 26 Mb, BTA11 at 27 Mb , and BTA12 at 66 Mb were associated with C14:1 cis-9 fatty acid. On BTA2 at 26 Mb region, glutamate decarboxylase 1 (GAD1) and specificity protein 5-transcription factor (Sp5) were identified. The GAD1 gene is involved with food behavior and insulin secretion. It was shown to be associated with morbid obesity [43] in humans, and body weight and daily gain in cattle [44]. Sp5 transcription factors are involved in the regulation of pyruvate kinase, lactate dehydrogenase and fatty acid synthase in cancer cells [45, 46]. On BTA11 at 27 Mb region, ATP-binding cassete, sub-family G (WHITE), member 5 (ABCG5) was identified. The ABCG family members were associated with cellular lipid-trafficking in macrophages and hepatocytes. This gene also presents an important role in the PPARγ and LXR pathways [47]. On BTA12 at 66 Mb region, glypican 6 (GPC6) was identified. Glypicans family are involved with the control of cell growth and cell division. Other glypican gene, GPC4 was associated with body fat distribution and obesity in humans [48]. A previous GWAS using chicken population reported an association between GPC6 and body weigh, which has a stronger correlation with body fatness and obesity [49].

BTA3 at 25 Mb, BTA7 at 80 Mb and BTA10 at 74 Mb were associated with C16:1 cis-9 fatty acid. Both BTA3 at 25 Mb and BTA10 at 74 Mb regions were also associated with saturated fatty acids, where HMGCS2, PHGDH, HSD3B1, HAO2, and DHRS7 were described above. On BTA7 at 80 Mb region, no PCG associated with this trait was identified.

BTA3 at 25 Mb, BTA3 at 26 Mb, and BTA12 at 13 Mb were associated with C18:1 cis-9. The first two regions were also associated with saturated fatty acids (C18:0 and C16:0, respectively) and C16:1 cis-9. On BTA12 at 13 Mb region, there were no annotated genes [50].

BTA6 at 50 Mb and BTA17 at 41 Mb were associated with C18:1 trans-16. On BTA 6 at 50 Mb region, no PCG associated with this trait was identified, and on BTA17 at 41 Mb region one PCG was found: Rap guanine nucleotide exchange factor (RAPGEF2), that selectively and non-covalently interacts with diacylglycerol, a diester of glycerol and two fatty acids (GO:0019992).

Four genomic regions were associated with saturated and monounsaturated fatty acids. BTA3 at 26 Mb and BTA12 at 13 Mb were associated with C16:0 and C18:1 cis-9. BTA3 at 25 Mb was associated with C18:0 and C18:1 cis-9. BTA3 at 25 and BTA10 at 74 Mb were associated with C18:0 and C16:1 cis-9. We also observed a negative phenotypic correlation between C16:0 and C18:1 (r = -0.65) and between C18:0 and C18:1(r = -0.72). Palmitic acid (C16:0) carbon chain skeleton is a predominant source for fatty acid elongation and desaturation, which might generate palmitoleic (C16:1 cis-9), stearic (C18:0), and oleic (C18:1 cis-9) fatty acids [51]. Examination of the 30 markers with the larger effect in regions of BTA3 at 26 Mb and BTA12 at 13 Mb (Table 4) and BTA3 at 25 Mb (Additional file 2) revealed that alleles positively associated with saturated fatty acid were negatively associated with unsaturated fatty acids. This indicates that QTL in these regions are directly or indirectly involved with elongation and/or desaturation.

Table 4 Top 30 markers effect in BTA3 at 26 Mb and BTA12 at 13 Mb QTL associated with C16:0 and C18:1 cis-9 in Nellore steers

Polyunsaturated fatty acids

Nine genomic regions (1 Mb window) explained > 1% of genotypic variation for group of polyunsaturated fatty acids, which relates C18:2 cis-9 cis12 n-6, C18:2 trans-11 cis-15, C18:3 n-3, C18:3 n-6, C20:3 n-6, C20:5 n-3, and C22:5 n-3 (Table 3). These overlapped with QTL reported for marbling score, body weight in Angus cattle [30], intramuscular fat, saturated fatty acid content, stearic acid in Fleckvieh bulls [52] and RFI in half-sib families from Angus and Charolais [53]. Additional file 3 shows in more detail the Manhattan plot of the proportion of genetic variance explained by window across the 29 autosomes for the important fatty acid to human health: γ-linolenic acid (C18:3 n-3), α-linolenic acid (C18:3 n-6) and total of n-3 (n-3).

BTA3 at 72 Mb, BTA 8 at 68 Mb, and BTA26 at 20 Mb were associated with C18:2 trans-11 cis-15 fatty acid. On BTA3 at 72 Mb QTL region, two PCG related to lipid metabolism were observed: aquaporin 7 (AQP7) and lysil oxidase-like2 (LOXL2). AQP7 is involved with the PPAR signaling pathway [54], while LOXL2 is associated with lean body mass in mouse [55]. On BTA 8 at 68 Mb and BTA26 at 20 Mb no PCG was associated with this trait.

BTA10 at 49 Mb region was associated with C20:5 n-3 fatty acid, which harbors the RAR-related orphan receptor A (RORA) gene that is related to steroid hormone receptor activity. When combined with a steroid hormone, it produces the signal within the cell to initiate a change in cell activity or function (GO:0003707).

BTA3 at 27 Mb and BTA10 at 29 Mb QTL regions explained more than 1% of genetic variance for C22:3 n-3. The first region (BTA3 at 27 Mb) harbors two PCG: SPAG17 and WDR3. These genes are involved with nucleus membrane (cellular components), the lipid bilayers that surround the nucleus and that form the nuclear envelope excluding the intermembrane space (GO:0031965). In the second region no PCG involved with this trait was identified. In other regions associated with C18:3 n-3 (BTA7 at 34 Mb), C18:3 n-6 (BTA17 at 24 Mb), and C20:3 n-6 (BTA27 at 26 Mb) no annotated genes were found [56].

The total of omega-3 (n-3) and omega-6 (n-6) GWAS result is consistent with individual polyunsaturated n-3 and n-6 fatty acids, where the same QTL regions were associated with these traits (BTA3 at 27 Mb and BTA 27 at 26 Mb).

GWAS have been used to investigate complex traits in many species including livestock [6]. Deposition and composition of fat in mammalians are complex traits and are influenced by many loci throughout the genome [57]. GWAS provides one approach to understand the genetic variation in complex traits, by identifying regions that can be fine-mapped to identify individual loci responsible for variation [58]. In the present GWAS two interesting QTL regions BTA3 at 25 Mb and BTA3 at 26 Mb were associated with saturated, monounsaturated, and polyunsaturated fatty acids in Bos indicus cattle. These regions harbor interesting PCG, which are involved with lipid metabolism in different species.

Previous studies reported CCAAT/enhancer binding proteins (C/EBPβ), peroxisome proliferator-activated receptor gamma (PPARγ or PPARG), carnitine palmitoyltransferase–1 beta (CPT–1β), stearoylcoenzyme A desaturase (SCD), AMP-activated protein kinase alpha (AMPKα), and G-coupled protein receptor 43 (GPR43) genes to be related to fat deposition and fatty acid composition in Angus [59]. In this GWAS and in a recent GWAS publication using Angus cattle, the regions that harbor these genes were not associated with fat deposition and fatty acid composition [60]. However, the positional candidate genes HMGCS2, PHGDH, HSD3B1 and HAO2 identified in the present study are involved in the PPAR signaling pathway in human (PATH: hsa03320), carbon metabolism (PATH:bta01200), lipid metabolism and steroid hormone biosynthesis (PATH: bta00140), carbohydrate metabolism and glycosylate and decarboxylase metabolism (PATH: bta00630) according to Kegg: Kyoto Encyclopedia of Genes and Genomes [61], respectively.

This is the first GWAS for intramuscular fat deposition and composition in Nellore. GWAS results in Bos taurus (Angus and Japanese Black cattle) reported different genomic regions associated with fat deposition and composition than those reported herein [62, 63]. Previously, it has been reported that a QTL on BTA19 was associated with fatty acid composition in Bos taurus breeds. This region on BTA19 harbors fatty acid synthase gene (FASN), which is an enzyme involved with de novo synthesis of long-chain fatty acid in mammalian, lipogenesis. FASN has been suggested as a candidate gene for fat traits in beef cattle [60, 62, 63].

Differences in SNP allele frequencies and linkage disequilibrium profile (LD between SNPs and causal variants) may explain the different marker effects between Bos indicus and Bos taurus cattle [64, 65]. This explanation has been confirmed by Bolormaa and collaborators [66], who compared marker effects using GWAS for Bos taurus and Bos indicus animals, which demonstrated that a SNP effect depends on the origin of alleles and the QTL segregation. The QTL segregation could result from mutation lost or fixation of alleles in one of the breeds, and also that the mutations occurred after divergence of these breeds [67]. Furthermore, it is possible that the differences in physiological and metabolism factors could contribute to the observed differences between different breeds [12].

Conclusion

The present study using BovineHD BeadChip (770 k) identified several 1-Mb SNP regions and genes within these regions that were associated with IMF and FA. The values of genomic heritabilities described in this study have not been reported before for Nellore cattle, and this information is important for breeding programs interesting in improving these traits. IMF deposition and composition are considered complex traits (polygenic) and are moderately heritable. In this study it is apparent that IMF composition are affected by many loci with small effects. Identification of several genomic regions and putative positional genes associated with lipid metabolism reported here should contribute to the knowledge of the genetic basis of IMF and FA deposition and composition in Nellore cattle (Bos indicus breed) and lead to selection for those traits to improve human nutrition and health.

Methods

Animals were handled and managed according to Institutional Animal Care and Use Committee Guidelines from Brazilian Agricultural Research Corporation – EMBRAPA approved by the president, Dr. Rui Machado.

Animals and phenotypes

Nellore steers (386) bred in the Brazilian Agricultural Research Corporation (EMBRAPA/Brazil) experimental breeding herd between 2009 and 2011 were available for this study. Steers were sired by 34 unrelated sires, and were selected to represent the main genealogies used in Brazil according to the National Summary of Nellore produced by the Brazilian Association of Zebu Breeders (ABCZ) and National Research Center for Beef Cattle. Animals were raised in feedlots under identical nutrition and handling conditions until slaughter at an average age of 25 months [68]. Steaks (2.54 cm thick) from the Longissimus dorsi muscle between the 12th and 13th ribs were collected 24 hours after slaughter.

Muscle samples (~100 g) were lyophilized and ground for IMF and FA analysis. The IMF was obtained using an Ankom XT20 extractor as described [69]. FA analysis was conducted as described by Hara and Radin [70], except the hexane to propanol ratio was increased to 3:2. Approximately 4 g of LD muscle was lyophilized, ground in liquid nitrogen, mixed with 28 mL of hexane/propanol (3:2 vol/vol) and homogenized for 1 min. Samples were vacuum filtered and 12 ml sodium sulfate (67 mg mL^{− 1}) solution was added and agitated for 30 s. The supernatant was transferred to a tube with 2 g of sodium sulfate and insufflated with N₂, after which the tube was sealed and incubated at room temperature for 30 min. Subsequently, the liquid was transferred to 10 mL test tube, insufflated with N₂, sealed and kept at − 20°C until dry with N₂ for methylation. The extracted lipids were hydrolyzed and methylated as described by Christie [71], except that hexane and methyl acetate were used instead of hexane:diethyl ether:formic acid (90:10:1). Around 40 mg of lipids were transferred to a tube containing 2 mL of hexane. Subsequently, 40 μL of methyl acetate were added, the sample agitated, and 40 μL of methylation solution (1.75 mL of methanol/0.4 mL of 5.4 mol/L of sodium metoxyde) were added. This mixture was agitated for 2 min and incubated for 10 min at room temperature. Then 60 μL of finishing solution (1 g of oxalic acid/30 mL of diethyl ether) were added and the mixture was agitated for 30 s, after which 200 mg of calcium chloride was added. The sample was then mixed and incubated at room temperature for 1 h. Samples were centrifuged at 3200 rpm, for 5 min at 5°C. The supernatant was collected for determination of fatty acids. Fatty acid methyl esters were quantified with a gas chromatograph (ThermoFinnigan, Termo Electron Corp., MA, USA) equipped with a flame ionization detector and a 100 m Supelco SP-2560 (Supelco Inc., PA, USA) fused silica capillary column (100 m, 0.25 mm and 0.2 μm film thickness). The column oven temperature was held at 70°C for 4 min, then increased to 170°C at a rate of 13°C min⁻¹, and subsequently increased to 250°C at a rate of 35°C min⁻¹, and held at 250°C for 5 min. The gas fluxes were 1.8 mL min⁻¹ for carrier gas (He), 45 mL min⁻¹ for make-up gas (N₂), 40 mL min⁻¹ for hydrogen, and 450 mL min⁻¹ for synthetic flame gas. One μL sample was analyzed. Injector and detector temperatures were 250 and 300°C, respectively. Fatty acids were identified by comparison of retention time of methyl esters of the samples with standards of fatty acids butter reference BCR-CRM 164, Anhydrous Milk Fat-Producer (BCR Institute for Materials and Reference Measurements) and also with commercial standard for 37 fatty acids Supelco TM Component FAME Mix (cat 18919, Supelco, Bellefonte, PA). The nomenclature of fatty acids follow IUPAC Compendium [72]. Fatty acids were quantified by normalizing the area under the curve of methyl esters using Chromquest 4.1 software (Thermo Electron, Italy). Fatty acids were expressed as a weight percentage (mg/mg). These analyses were performed at the Animal Nutrition and Growth Laboratory at ESALQ, Piracicaba, São Paulo, Brazil.

DNA extraction and genotypic data

DNA was isolated from blood as described by Tizioto et al. [68]. Genotyping was performed at the Bovine Functional Genomics Laboratory ARS/USDA and Genomics Center at ESALQ, Piracicaba, São Paulo, Brazil using BovineHD 770 k BeadChip (Infinium BeadChip, Illumina, San Diego, CA) according to manufacturer’s protocol. Genotypes were obtained in Illumina A/B allele format and used to represent a covariate value at each locus coded as 0, 1, or 2, representing the number of B alleles. Missing genotypes, represented < 0.2% of genotypes and were replaced with the average covariate value at that locus. Initial visualization and data analysis was performed by GenomeStudio Data Analysis Software [73]. The SNPs with call rate ≤ 95%, minor allele frequency (MAF) ≤ 5%, those located on sex chromosomes and those not mapped in the Bos taurus UMD 3.1 assembly were removed. After filtering, a total of 449,363 SNP were utilized in GWAS.

Descriptive statistics and heritability

Descriptive statistics for IMF and FA were estimated using PROC MEANS and normality tests were performed using PROC UNIVARIATE in SAS (Ver. 9.3; SAS Inst. Cary, NC). SAS PROC MIXED was used to test independent sources for significance. Fixed effects included contemporary group classes (animals with the same origin, birth year and slaughter date) and hot carcass weight as a covariate. Animal and residuals were fitted as random effects. Restricted maximum likelihood was used to estimate animal and residual variance components, heritability and standard error (SE) using ASREML software [74]. The model used in single-trait analyses of all traits was, y = Xb + Zu + e, where y is the vector of observations representing the trait of interest (dependent variable), X and Z are the design or incidence matrices for the vectors of fixed and random effects in b and u, respectively, and e was the vector of random residuals. The variance of vector u was Gσ²_m for the genomic analyses where G is the genomic relationship matrix derived from SNP markers using allele frequencies as suggested by VanRaden [75], with σ²_m being the marker-based additive genetic variance.

Genome wide association study

Associations between SNP and phenotypes (IMF and FA) were obtained using Bayes B, which analyzed all SNP data simultaneously and assumed a different genetic variance for each SNP locus [76, 77]. The prior genetic and residual variances were estimated using Bayes C [78], with π being 0.9997. The model equation was:

\mathit{y}=\mathit{Xb}+{\displaystyle \sum _{\mathit{j}=1}^{\mathit{k}}}{\mathit{a}}_{\mathit{j}}{\mathit{\beta }}_{\mathit{j}}{\mathit{\delta }}_{\mathit{j}}+\mathit{e},

where y was the vector of the phenotypic values, X was the incidence matrix for fixed effects, b was the vector of fixed effects defined above, k was the number of SNP loci (449,363), a _j was the column vector representing the SNP covariate at locus j coded as the number of B alleles, β _j was the random substitution effect for locus j, which conditional on σ² _β was assumed to be normally distributed N (0, σ² _β ) when δ _j = 1 but β _j = 0 when δ _j = 0, with δ _j being a random 0/1 variable indicating the absence (with probability π) or presence (with probability 1-π) of locus j in the model, and e was the vector of the random residual effects assumed normally distributed N (0, σ² _e ). The variance σ² _β (or σ² _e ) was a priori assumed to follow a scaled inverse Chi-square with v _β = 4 (or v _e = 10) degrees of freedom and scale parameter S² _β (or S² _e ). The scale parameter for markers was derived as a function of the assumed known genetic variance of the population, based on the average SNP allele frequency and number of SNP assumed to have nonzero effects based on parameter π being 0.9997. This procedure used GenSel software [8] to obtain the posterior distributions of SNP effects using Markov chain Monte Carlo (MCMC). This comprised a burn-in period of 1,000 iterations from which results were discarded, followed by 40,000 iterations from which results were accumulated to obtain the posterior mean effect of each SNP. In the Bayesian variable selection multiple-regression models with π = 0.9997 about 100-150 SNP markers were fitted simultaneously in each MCMC iteration. Inference of associations in these multiple-regression models was based on 1-Mb genomic windows rather than on single markers [8, 79]. Genomic windows were constructed from the chromosome and base-pair positions denoted in the marker map file [8] based on UMD3.1 bovine assembly.

The SNP effects from every 40^th post burn in iteration were used to obtain samples from the posterior distribution of the proportion of variance accounted for by each window from 1,000 MCMC samples of genomic merit for each animal following Onteru et al. [79] and Peters et al. [7]. In the present study there were 2,527 1 Mb SNP windows across the 29 autosomes. The proportion of genetic variance explained by each window in any particular iteration was obtained by dividing the variance of window BV by the variance of whole genome BV in that iteration. The window BV was computed by multiplying the number of alleles that represent the SNP covariates for each consecutive SNP in a window by their sampled substitution effects in that iteration.

All traits were used for GWAS (Additional file 4) and genomic heritability estimate, but only the ones with genomic heritability ≥ 0.10 were reported. Fatty acids were indexed as groups of saturated, monounsaturated, polyunsaturated fatty acid, total of saturated fatty acid (SFA), total monounsaturated (MUFA), total of polyunsaturated (PUFA), total of omega 3 (n-3) and total of omega 6 (n-6). Genome windows with the highest posterior mean proportion of genetic variance ≥1% were considered the most important regions associated with the traits, and were declared the most promising QTL regions.

Positional candidate genes were investigated for 1 Mb windows using the Cattle Genome Browser [50] and UCSC Genome Browser [80], which allowed visualization of SNP based on Bos taurus genome assembly UMD 3.1. Animal QTL database (Animal QTLdb) was used to search for published QTL and trait mapping data [56]. Gene annotations were retrieved from Ensembl Genes 71 Database using Biomart software [54, 81]. The functional classification of genes was done using DAVID [54] and BioGPS [38] online annotation databases. Those genes reported to be involved in fatty acid and lipid metabolism were selected as positional candidate genes.

Availability of supporting data

The data sets supporting the results of this article are included within the article and its additional files.