Edited by:
Reviewed by:
*Correspondence:
This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
The fatty acid profile of beef is a complex trait that can benefit from gene-interaction network analysis to understand relationships among loci that contribute to phenotypic variation. Phenotypic measures of fatty acid profile from triacylglycerol and phospholipid fractions of longissimus muscle, pedigree information, and Illumina 54 k bovine SNP genotypes were utilized to derive an annotated gene network associated with fatty acid composition in 1,833 Angus beef cattle. The Bayes-B statistical model was utilized to perform a genome wide association study to estimate associations between 54 k SNP genotypes and 39 individual fatty acid phenotypes within each fraction. Posterior means of the effects were estimated for each of the 54 k SNP and for the collective effects of all the SNP in every 1-Mb genomic window in terms of the proportion of genetic variance explained by the window. Windows that explained the largest proportions of genetic variance for individual lipids were found in the triacylglycerol fraction. There was almost no overlap in the genomic regions explaining variance between the triacylglycerol and phospholipid fractions. Partial correlations were used to identify correlated regions of the genome for the set of largest 1 Mb windows that explained up to 35% genetic variation in either fatty acid fraction. SNP were allocated to windows based on the bovine UMD3.1 assembly. Gene network clusters were generated utilizing a partial correlation and information theory algorithm. Results were used in conjunction with network scoring and visualization software to analyze correlated SNP across 39 fatty acid phenotypes to identify SNP of significance. Significant pathways implicated in fatty acid metabolism through GO term enrichment analysis included homeostasis of number of cells, homeostatic process, coenzyme/cofactor activity, and immunoglobulin. These results suggest different metabolic pathways regulate the development of different types of lipids found in bovine muscle tissues. Network analysis using partial correlations and annotation of significant SNPs can yield information about the genetic architecture of complex traits.
Beef is a nutritious source of protein, fat, vitamins, and minerals when appropriately included in the human diet. A large body of research indicates it is critical to maintain a balanced fatty acid intake to support a healthy blood lipid profile (
The fatty acid profile in beef cattle can be characterized by the abundance of 39 individual lipids of varying chain lengths and degrees of saturation (
Phenotypic correlations among fatty acids within and across lipid fractions have been published (
The Iowa State University and Oklahoma State University Institutional Review Boards approved the experimental protocols used in this study.
A total of 1,833 offspring of 155 Angus bulls represented by bulls (
Genomic DNA was extracted from the ground beef sample used for fatty acid composition and was genotyped with the Bovine SNP50 Infinium II BeadChip (Illumina, San Diego, CA, USA). Contemporary groups were defined based on cross-classifications of gender at harvest (bull, steer or heifer), finishing location (California, Colorado, Iowa, Texas), and harvest date, for a total of 33 groups. Contemporary groups were fit as fixed effects in genomic analyses. Effects of SNP on each trait were estimated using the Bayes-B option of GenSel accessed through the BIGSGUI Version 0.9.2 (
Significant pair-wise correlations among SNP effects were calculated across the 39 fatty acid phenotypes using an optimized implementation of the PCIT algorithm. Briefly, this algorithm begins by calculating third-order partial correlation coefficients among all possible trios of SNPs. Within each trio, a direct correlation between any pair of SNPs that does not exhibit a significant partial correlation to the third SNP is considered an independent association, and thus not significant for network co-association. This methodology identifies SNP co-associations that occur more frequently over the range of phenotypes given as input. The full mathematical model for identifying the most co-associated SNPs across multiple phenotypes is explained in detail by
All SNP were selected from 20 genomic windows of size 1 Mb with the largest posterior probability of association from the 16:0 fatty acid phenotype in order to build a matrix of SNP effects for the phospholipid fraction. A vector of posterior mean SNP effects for 418 SNP from 16:0 was augmented with the effects estimated for each of the other 38 fatty acids. PCIT network creation and visualization proceeded identically to the methods described for the triacylglycerol fraction.
Correlations among SNP were used to construct networks of SNP that exhibited common effects across multiple fatty acids. Correlation between SNP pairs with a non-zero partial correlation to another SNP were input into Cytoscape 3.0.2 (
Gene ontology (GO) enrichment was carried out using the DAVID v6.7 Functional Annotation Tool (
The posterior mean estimates for proportion of genetic variance explained by 1 Mb genomic windows as well as the posterior probability of association for selected lipids and lipid classes from the triacylglycerol fraction are shown in
Characterization of top 20 1 Mb genomic windows that account for variation in triacylglyceride lipids.
Trait | Map Pos | RS# Start | RS# End | # of SNP | σ2g (%) | PPA |
---|---|---|---|---|---|---|
14:0 | ||||||
19_51 | rs41923412 | rs109147235 | 25 | 34.55 | 1.00 | |
29_18 | rs42375315 | rs41589183 | 14 | 10.65 | 1.00 | |
10_19 | rs41647457 | rs110785951 | 24 | 3.21 | 0.98 | |
16:0 | ||||||
17_16 | rs109550465 | rs110684903 | 17 | 22.91 | 1.00 | |
7_56 | rs41614823 | rs42334377 | 17 | 16.26 | 1.00 | |
19_51 | rs41923412 | rs109147235 | 25 | 16.18 | 1.00 | |
1_80 | rs43245574 | rs110467946 | 21 | 10.44 | 1.00 | |
16_3 | rs41790571 | rs41633905 | 24 | 9.66 | 0.99 | |
29_18 | rs42375315 | rs41589183 | 14 | 5.87 | 1.00 | |
16:1 | ||||||
19_51 | rs41923412 | rs109147235 | 25 | 15.25 | 1.00 | |
29_18 | rs42375315 | rs41589183 | 14 | 4.72 | 1.00 | |
26_21 | rs109309604 | rs42086690 | 20 | 3.23 | 1.00 | |
18:0 | ||||||
29_18 | rs42375315 | rs41589183 | 14 | 5.95 | 1.00 | |
18:1 | ||||||
19_51 | rs41923412 | rs109147235 | 25 | 22.13 | 1.00 | |
8_103 | rs109285764 | rs41590918 | 18 | 6.43 | 0.97 | |
16_20 | rs110743197 | rs42542723 | 23 | 5.88 | 0.99 | |
29_18 | rs42375315 | rs41589183 | 14 | 4.96 | 1.00 | |
SFA | ||||||
19_51 | rs41923412 | rs109147235 | 25 | 15.58 | 1.00 | |
26_20 | rs42981135 | rs41623887 | 21 | 5.39 | 0.98 | |
MUFA | ||||||
19_51 | rs41923412 | rs109147235 | 25 | 16.63 | 1.00 | |
Other genomic regions of significance that appeared in both this data set that separately considers the two fractions and the total fatty acid analysis presented by
Several genomic regions of interest were found to explain relatively large proportions of genetic variance that were not detected in the data presented by
The posterior mean estimates for genetic variance explained by 1 Mb genomic windows as well as the posterior probabilities of association for selected lipids and lipid classes from the phospholipid fraction are shown in
Characterization of top 20 1 Mb genomic windows accounting for variation in phospholipid lipids.
Trait | Map Start | RS# Start | RS# End | # of SNP | σ2g (%) | PPA |
---|---|---|---|---|---|---|
14:0 | ||||||
16_4 | rs110257825 | rs109105804 | 26 | 1.79 | 0.68 | |
19_5 | rs41633989 | rs109106774 | 17 | 1.11 | 0.50 | |
16:0 | ||||||
21_36 | rs109143576 | rs42429437 | 22 | 2.54 | 0.72 | |
1_52 | rs41600017 | rs43711327 | 25 | 1.10 | 0.54 | |
16:1 | ||||||
4_95 | rs43412327 | rs42421263 | 20 | 1.39 | 0.68 | |
X_5 | rs109239523 | rs29023191 | 12 | 0.87 | 0.55 | |
18:0 | ||||||
24_29 | rs110012069 | rs42837712 | 24 | 4.27 | 0.78 | |
18:2, n-6 | ||||||
27_42 | rs42135519 | rs110741211 | 25 | 3.08 | 0.90 | |
18:3, n-3 | ||||||
9_12 | rs43731273 | rs110226869 | 21 | 3.52 | 0.99 | |
4_89 | rs109964815 | rs43586675 | 24 | 2.85 | 0.99 | |
18_31 | rs41574692 | rs41581150 | 11 | 2.93 | 0.98 | |
26_27 | rs109158324 | rs41636608 | 18 | 3.13 | 0.95 | |
1_24 | rs29011682 | rs29017639 | 27 | 1.81 | 0.88 | |
15_29 | rs29010888 | rs41582064 | 25 | 4.25 | 0.85 | |
17_70 | rs42794376 | rs109349100 | 21 | 2.02 | 0.79 | |
20:0 | ||||||
9_54 | rs109669651 | rs110216218 | 25 | 3.18 | 0.91 | |
25_6 | rs41592046 | rs110848072 | 25 | 2.04 | 0.74 | |
1_73 | rs109693922 | rs110962722 | 24 | 1.82 | 0.72 | |
20:1 | ||||||
9_104 | rs41636894 | rs41592131 | 22 | 3.42 | 0.93 | |
SFA | ||||||
3_114 | rs110764304 | rs109271147 | 30 | 4.3 | 0.63 | |
Several windows identified harbor candidate genes related to overall lipid and phospholipid metabolism. The genomic window on chromosome 16 at 4 Mb that accounted for 1.79% of the genetic variance in 14:0 harbors the candidate gene fructose-2,6-biphosphatase 2 (PFKFB2). This gene is known to be involved in synthesis and degradation of fructose-2,6-bisphosphate, which is a regulatory molecule involved in glycolysis in eukaryotes (
Some 394 of the 454 SNP entered into the PCIT for the triacylglycerol fraction were co-localized into 12 separate networks. Information detailing the top 5 network scoring results from the Cytoscape MCODE plugin for each triacylglycerol-derived network is in
MCODE results derived from network clustering with PCIT.
Fraction | Network | Score | Nodes | Edges |
---|---|---|---|---|
Triacylglycerol | ||||
1 | 53.26 | 56 | 1487 | |
2 | 37.24 | 92 | 1699 | |
3 | 23.56 | 58 | 559 | |
4 | 19.26 | 78 | 662 | |
5 | 18.18 | 38 | 278 | |
Phospholipid | ||||
1 | 33.91 | 52 | 987 | |
2 | 29.40 | 65 | 1034 | |
3 | 26.10 | 82 | 1064 | |
4 | 10.00 | 66 | 468 | |
5 | 7.40 | 51 | 365 | |
Candidate genes involved in fatty acid metabolism found within these networks include THRSP, Acyl-CoA synthetase-5 (ACSL5), glycerol-3-phosphate acyltransferase muscle-type (GPAM), and coiled coil domain-containing 57 (CCD57). The candidate genes THRSP and GPAM have been previously identified as playing a role in lipid metabolism in beef and dairy adipose via the PPAR pathway (
Annotated visualizations of the two highest scoring phospholipid networks with Cytoscape are in
Gene ontology term enrichment analysis was carried out for all genes located in the top 1 Mb regions shown in
DAVID Functional Annotation Clustering for significant 1 Mb windows identified through genome-wide association studies of the Triacylglycerol fraction.
Annotation Cluster 1 | Enrichment Score: 1.71 | ||||
---|---|---|---|---|---|
Category | Term | Count | % | FE | |
GOTERM_BP_FAT | GO:0048872∼homeostasis of number of cells | 6 | 4.32 | <0.01 | 11.52 |
GOTERM_BP_FAT | GO:0042592∼homeostatic process | 8 | 6.12 | <0.01 | 2.96 |
GOTERM_BP_FAT | GO:0030099∼myeloid cell differentiation | 3 | 2.54 | 0.04 | 8.49 |
GOTERM_BP_FAT | GO:0048534∼hemopoietic or lymph organ dev. | 6 | 2.65 | 0.05 | 6.02 |
GOTERM_BP_FAT | GO:0002520∼immune system development | 4 | 2.84 | 0.06 | 4.56 |
GOTERM_BP_FAT | GO:0030097∼hemopoiesis | 4 | 3.65 | 0.18 | 3.25 |
SP_PIR_KEYWORDS | nadp | 5 | 2.32 | <0.01 | 9.25 |
SP_PIR_KEYWORDS | nad | 6 | 4.17 | <0.01 | 6.54 |
INTERPRO | IPR016040:NAD(P)-binding domain | 6 | 3.21 | 0.02 | 5.66 |
UP_SEQ_FEATURE | active site: Proton acceptor | 6 | 5.84 | 0.02 | 4.57 |
SP_PIR_KEYWORDS | oxidoreductase | 8 | 5.65 | 0.02 | 1.99 |
GOTERM_BP_FAT | GO:0055114∼oxidation reduction | 3 | 6.98 | 0.03 | 2.26 |
The full list of GO terms is displayed in the Functional Annotation Chart results in
DAVID Functional Annotation Chart results for significant 1 Mb regions identified through genome-wide association studies of the Triacylglycerol fraction.
Category | Term | Count | % | FE | FDR | |
---|---|---|---|---|---|---|
GOTERM_BP_FAT | GO:0048872∼homeostasis of number of cells | 6 | 3.18 | <0.01 | 11.12 | 1.08 |
SP_PIR_KEYWORDS | nadp | 5 | 3.32 | <0.01 | 9.35 | 3.06 |
SP_PIR_KEYWORDS | nad | 6 | 4.56 | <0.01 | 6.44 | 3.50 |
GOTERM_BP_FAT | GO:0042592∼homeostatic process | 9 | 6.48 | <0.01 | 2.66 | 4.70 |
GOTERM_BP_FAT | GO:0016358∼dendrite development | 4 | 2.20 | <0.01 | 25.19 | 5.35 |
GOTERM_BP_FAT | GO:0000187∼activation of MAPK activity | 3 | 2.40 | 0.01 | 17.56 | 15.71 |
INTERPRO | IPR016040:NAD(P)-binding domain | 5 | 3.51 | 0.02 | 5.03 | 12.45 |
GOTERM_CC_FAT | GO:0016607∼nuclear speck | 3 | 2.95 | 0.02 | 12.48 | 18.97 |
UP_SEQ_FEATURE | active site:Proton acceptor | 5 | 4.65 | 0.02 | 4.15 | 20.14 |
GOTERM_BP_FAT | GO:0043406∼positive regulation of MAP kinase activity | 3 | 2.32 | 0.02 | 12.53 | 36.20 |
KEGG_PATHWAY | bta04130:SNARE interactions in vesicular transport | 5 | 2.65 | 0.03 | 9.26 | 45.25 |
GOTERM_BP_FAT | GO:0043113∼receptor clustering | 3 | 1.43 | 0.03 | 45.62 | 32.65 |
GOTERM_CC_FAT | GO:0043233∼organelle lumen | 6 | 5.85 | 0.03 | 2.35 | 32.54 |
GOTERM_BP_FAT | GO:0030099∼myeloid cell differentiation | 3 | 2.45 | 0.03 | 6.95 | 39.56 |
GOTERM_CC_FAT | GO:0005739∼mitochondrion | 3 | 6.84 | 0.03 | 3.62 | 38.32 |
GOTERM_BP_FAT | GO:0007172∼signal complex assembly | 2 | 2.10 | 0.03 | 59.63 | 47.65 |
GOTERM_CC_FAT | GO:0031974∼membrane-enclosed lumen | 8 | 6.56 | 0.03 | 3.26 | 39.25 |
SP_PIR_KEYWORDS | Chaperone | 4 | 2.78 | 0.04 | 5.64 | 42.33 |
GOTERM_BP_FAT | GO:0043405∼regulation of MAP kinase activity | 3 | 2.18 | 0.05 | 7.75 | 49.62 |
GOTERM_CC_FAT | GO:0016604∼nuclear body | 3 | 2.86 | 0.05 | 8.56 | 45.25 |
SMART | SM00241:ZP | 2 | 1.62 | 0.05 | 26.69 | 36.21 |
UP_SEQ_FEATURE | nucleotide phosphate-binding region:NAD | 3 | 2.32 | 0.05 | 7.56 | 35.26 |
Functional Annotation Clustering analysis for the top GWAS regions associated with phospholipid is in
DAVID Functional Annotation Clustering for significant 1 Mb regions identified through GWAS for the Phospholipid fraction.
Annotation Cluster 1 | Enrichment Score: 1.38 | ||||
---|---|---|---|---|---|
Category | Term | Count | % | FE | |
INTERPRO | IPR013783:Immunoglobulin-like fold | 6 | 7.82 | <0.01 | 6.32 |
INTERPRO | IPR007110:Immunoglobulin-like | 8 | 6.94 | <0.01 | 4.98 |
INTERPRO | IPR013106:Immunoglobulin V-set | 3 | 5.02 | 0.01 | 7.84 |
SP_PIR_KEYWORDS | Immunoglobulin domain | 5 | 5.19 | 0.01 | 7.89 |
INTERPRO | IPR003599:Immunoglobulin subtype | 4 | 3.56 | 0.13 | 5.62 |
SMART | SM00409:IG | 5 | 4.45 | 0.18 | 3.21 |
UP_SEQ_FEATURE | glycosylationsite: N-linked (GlcNAc…) | 6 | 8.95 | 0.20 | 2.12 |
SP_PIR_KEYWORDS | glycoprotein | 9 | 7.56 | 0.58 | 2.31 |
DAVID Functional Annotation Chart results for significant 1 Mb regions identified through GWAS for the Phospholipid fraction.
Category | Term | Count | % | FE | FDR | |
---|---|---|---|---|---|---|
INTERPRO | IPR013783:Immunoglobulin-like fold | 9 | 9.21 | <0.01 | 6.32 | 1.22 |
INTERPRO | IPR007110:Immunoglobulin-like | 5 | 8.45 | <0.01 | 6.56 | 4.56 |
INTERPRO | IPR013106:Immunoglobulin V-set | 6 | 6.54 | 0.01 | 9.01 | 10.10 |
SP_PIR_KEYWORDS | Immunoglobulin domain | 4 | 3.98 | 0.02 | 9.65 | 11.89 |
GOTERM_CC_FAT | GO:0048475∼coated membrane | 4 | 3.68 | 0.02 | 11.25 | 19.56 |
GOTERM_CC_FAT | GO:0030117∼membrane coat | 4 | 2.99 | 0.03 | 12.45 | 21.54 |
INTERPRO | IPR000920:MyelinP0 protein | 3 | 3.21 | 0.03 | 78.65 | 24.82 |
GOTERM_BP_FAT | GO:0046907∼intracellular transport | 6 | 6.32 | 0.03 | 5.05 | 35.18 |
GOTERM_MF_FAT | GO:0000166∼nucleotide binding | 11 | 11.25 | 0.03 | 2.03 | 36.32 |
GOTERM_CC_FAT | GO:0030663∼COPI coated vesicle membrane | 3 | 3.21 | 0.04 | 44.52 | 35.02 |
GOTERM_CC_FAT | GO:0030126∼COPI vesiclecoat | 2 | 2.65 | 0.04 | 42.17 | 35.26 |
GOTERM_CC_FAT | GO:0030137∼COPI-coated vesicle | 2 | 2.32 | 0.05 | 43.65 | 39.15 |
GOTERM_BP_FAT | GO:0006886∼intracellular protein transport | 4 | 4.84 | 0.05 | 5.89 | 54.55 |
Analysis of GWAS results for triacylglycerol and phospholipid fractions supports the conclusion that the triacylglycerol fraction is closely representative of the total fatty acid fraction. Significant genomic windows identified for triacylglycerol fraction overlapped with the results previously presented by (
An analysis of the genomic regions that affect the phospholipid fraction yielded few genes with a known biological association to lipid metabolism. Significant genomic regions identified explained lower percentages of genetic variance in comparison to the triacylglycerol. The low variation observed in the phospholipid fraction is likely due to the importance of the phospholipid membrane to biological function of the cell. Pathways prevalent in the phospholipid analysis appeared to be highly related to cell-to-cell adhesion, cellular trafficking, and coenzyme/cofactor activity. A larger dataset could possibly improve results when dealing with traits that exhibit low phenotypic variance. In conclusion, the combination of GWAS results with PCIT and network visualization represents a robust methodology for identifying candidate genes of interest for traits with multiple phenotypes and adequate phenotypic variance.
JB conducted the analysis and drafted the manuscript. JR conceived the analysis and assisted with the manuscript. DG assisted with the analysis and manuscript. QD assisted with data collection and analysis. DB assisted with he analysis and the manuscript. JK conceived the analysis and assisted with the manuscript. MS conceived the analysis and assisted with the manuscript. LK assisted with the methodology and and analysis. RM conceived the analysis and assisted with the manuscript.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
This project is supported by USDA-NIFA Award 2012-67015-19420 – Enhanced Bioinformatics to Implement Genomic Selection (e-BIGS). The authors also acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. URL:
The Supplementary Material for this article can be found online at: