Advertisement

Genomic analyses for predicted milk fatty acid composition throughout lactation in North American Holstein cattle

Open ArchivePublished:May 14, 2020DOI:https://doi.org/10.3168/jds.2019-17628

      ABSTRACT

      Milk fat composition has important implications in the nutritional and processing properties of milk. Additionally, milk fat composition is associated with cow physiological and health status. The main objectives of this study were (1) to estimate genetic parameters for 5 milk fatty acid (FA) groups (i.e., short-chain, medium-chain, long-chain, saturated, and unsaturated) predicted from milk infrared spectra using a large data set; (2) to predict genomic breeding values using a longitudinal single-step genomic BLUP approach; and (3) to conduct a single-step GWAS aiming to identify genomic regions, candidate genes, and metabolic pathways associated with milk FA, and consequently, to understand the underlying biology of these traits. We used 629,769 test-day records of 201,465 first-parity Holstein cows from 6,105 herds. A total of 8,865 genotyped (Illumina BovineSNP50K BeadChip, Illumina, San Diego, CA) animals were considered for the genomic analyses. The average daily heritability ranged from 0.24 (unsaturated FA) to 0.47 (medium-chain and saturated FA). The reliability of the genomic breeding values ranged from 0.56 (long-chain fatty acid) to 0.74 (medium-chain fatty acid) when using the default τ and ω scaling parameters, whereas it ranged from 0.58 (long-chain fatty acid) to 0.73 (short-chain fatty acid) when using the optimal τ and ω values (i.e., τ = 1.5 and ω = 0.6), as defined in a previous study in the same population. Relevant chromosomal regions were identified in Bos taurus autosomes 5 and 14. The proportion of the variance explained by 20 adjacent single nucleotide polymorphisms ranged from 0.71% (saturated FA) to 15.12% (long-chain FA). Important candidate genes and pathways were also identified. In summary, our results contribute to a better understanding of the genetic architecture of predicted milk FA in dairy cattle and reinforce the relevance of using genomic information for genetic analyses of these traits.

      Key words

      INTRODUCTION

      The proportion of fat (i.e., lipids) has great economic and nutritional value in milk composition, and it is directly related to flavor and chemical-physical characteristics of the milk and milk products and to cheese-making properties (
      • Baer R.J.
      Alteration of the fatty acid content of milk fat.
      ;
      • Bergamaschi M.
      • Stocco G.
      • Valorz C.
      • Bazzoli I.
      • Sturaro E.
      • Ramanzin M.
      Cheesemaking in highland pastures: Milk technological properties, cream, cheese and ricotta yields, milk nutrients recovery, and products composition.
      ). In this regard, fatty acids (FA) are important components of milk fat, comprising 90% of its weight (
      • Samková E.
      • Špička J.
      • Pešek M.
      • Pelikánová T.
      • Hanuš O.
      Review: Animal factors affecting fatty acid composition of cow milk fat: A review.
      ). The FA vary in chain length (number of carbon atoms), saturation, and arrangement. Milk FA originating from de novo synthesis in the mammary gland are usually classified as short-chain FA (SCFA; from 4 to 8 or 10 carbons) and medium-chain FA (MCFA; 10 or 12 to 14 or 16 carbons) chain, whereas FA from the bloodstream and diet are usually classified as long-chain FA (LCFA; 16 to 18 carbons;
      • Cozma A.
      • Miere D.
      • Filip L.
      • Banc R.
      • Loghin F.
      A review of the metabolic origins of milk fatty acids.
      ;
      • Fleming A.
      • Schenkel F.S.
      • Chen J.
      • Malchiodi F.
      • Bonfatti V.
      • Ali R.A.
      • Mallard B.
      • Corredig M.
      • Miglior F.
      Prediction of milk fatty acid content with mid-infrared spectroscopy in Canadian dairy cattle using differently distributed model development sets.
      ). Fatty acids can also be classified as saturated (i.e., FA that have no double bonds between the individual carbon atoms) and unsaturated (i.e., FA that have at least one double bond). In brief, SFA account for 70 to 75% of total fat in cow milk and comprise FA containing from 4 to 18 carbons. The most common SFA are palmitic (C16:0), stearic (C18:0), and myristic (C14:0) acids (
      • MacGibbon A.K.H.
      • Taylor M.W.
      Composition and Structure of Bovine Milk Lipids. Advanced Dairy Chemistry. Volume 2 Lipids. 1–42.
      ). The UFA account for 25 to 30% of the total FA in milk and are characterized by the presence of a cis-double bond between carbons 9 and 10, in a carbon chain length that ranges from 10 to 18 (
      • Ntambi J.M.
      The regulation of stearoyl-CoA desaturase (SCD).
      ). The lower amount of UFA in milk compared with SFA is due, in particular, to biohydrogenation that occurs in the rumen (
      • Sauer F.D.
      • Fellner V.
      • Kinsman R.
      • Kramer J.K.G.
      • Jackson H.A.
      • Lee A.J.
      • Chen S.
      Methane output and lactation response in Holstein cattle with monensin or unsaturated fat added to the diet.
      ).
      Consumer choices today are based not only on the nutritional aspects of the food, but also on products known to promote better health or prevent disease. In this regards, the proportion of SFA in milk is a concern (
      • Bilal G.
      • Cue R.I.
      • Mustafa A.F.
      • Hayes J.F.
      Short communication: Genetic parameters of individual fatty acids in milk of Canadian Holsteins.
      ). This is because SFA are associated with increased risk of cardiovascular diseases, obesity, and weight gain, as they can increase levels of total and low-density lipoprotein cholesterol in the blood (
      • Temme E.H.
      • Mensink R.P.
      • Hornstra G.
      Comparison of the effects of diets enriched in lauric, palmitic, or oleic acids on serum lipids and lipoproteins in healthy women and men.
      ;
      • Haug A.
      • Høstmark A.T.
      • Harstad O.M.
      Bovine milk in human nutrition—A review.
      ;
      • Zong G.
      • Li Y.
      • Wanders A.J.
      • Alssema M.
      • Zock P.L.
      • Willett W.C.
      • Hu F.B.
      • Sun Q.
      Intake of individual saturated fatty acids and risk of coronary heart disease in US men and women: Two prospective longitudinal cohort studies.
      ;
      • Briggs M.A.
      • Petersen K.S.
      • Kris-Etherton P.M.
      Saturated fatty acids and cardiovascular disease: Replacements for saturated fat to reduce cardiovascular risk.
      ). In contrast, some FA have antiallergic, antimicrobial, anticarcinogenic and anti-inflammatory properties (
      • Williams C.M.
      Dietary fatty acids and human health.
      ;
      • Parodi P.
      Milk fat in human nutrition.
      ;
      • Haug A.
      • Høstmark A.T.
      • Harstad O.M.
      Bovine milk in human nutrition—A review.
      ;
      • Pepe G.
      • Tenore G.C.
      • Mastrocinque R.
      • Stusio P.
      • Campiglia P.
      Potential anticarcinogenic peptides from bovine milk.
      ;
      • Lordan R.
      • Zabetakis I.
      Invited review: The anti-inflammatory properties of dairy lipids.
      ). Therefore, increasing the proportion of favorable FA in milk can contribute to improved human health and disease resistance. In this context, selection for improved FA composition will enable the dairy industry to label products with a more favorable milk FA profile. This could have a positive effect on consumer preferences and, consequently, on industry profitability.
      In summary, milk fat composition can be altered in various ways, including dietary changes and genetic selection (
      • Palmquist D.L.
      • Beaulieu A.D.
      • Barbano D.M.
      Feed and animal factors influencing milk fat composition.
      ;
      • Kęsek M.
      • Szulc T.
      • Zielak-Steciwko A.
      Genetic, physiological and nutritive factors affecting the fatty acid profile in cows' milk—A review.
      ;
      • Narayana S.G.
      • Schenkel F.S.
      • Fleming A.
      • Koeck A.
      • Malchiodi F.
      • Jamrozik J.
      • Johnston J.
      • Sargolzaei M.
      • Miglior F.
      Genetic analysis of groups of mid-infrared predicted fatty acids in milk.
      ). Changes due to diet are faster to observe than those due to genetic selection; however, the latter results in permanent and cumulative changes at the population level. Thus, to select for milk FA composition, mid-infrared spectroscopy (MIR) technology has become an attractive method to accurately predict FA composition in milk (
      • Soyeurt H.
      • Dardenne P.
      • Dehareng F.
      • Lognay G.
      • Veselko D.
      • Marlier M.
      • Bertozzi C.
      • Mayeres P.
      • Gengler N.
      Estimating fatty acid content in cow milk using mid-infrared spectrometry.
      ;
      • Narayana S.G.
      • Schenkel F.S.
      • Fleming A.
      • Koeck A.
      • Malchiodi F.
      • Jamrozik J.
      • Johnston J.
      • Sargolzaei M.
      • Miglior F.
      Genetic analysis of groups of mid-infrared predicted fatty acids in milk.
      ). Thereafter, milk FA predicted based on MIR can be used to evaluate animal performance and to predict EBV of the animals along the lactation curve using random regression models (RRM; e.g.,
      • Narayana S.G.
      • Schenkel F.S.
      • Fleming A.
      • Koeck A.
      • Malchiodi F.
      • Jamrozik J.
      • Johnston J.
      • Sargolzaei M.
      • Miglior F.
      Genetic analysis of groups of mid-infrared predicted fatty acids in milk.
      ;
      • Freitas P.H.F.
      • Oliveira H.R.
      • Silva F.F.
      • Fleming A.
      • Schenkel F.S.
      • Miglior F.
      • Brito L.F.
      Short communication: Time-dependent genetic parameters and single-step genome-wide association analyses for predicted milk fatty acid composition in Ayrshire and Jersey dairy cattle.
      ).
      The calculation of EBV for each lactation day enables modeling of the shape of the lactation curve (
      • Oliveira H.R.
      • Brito L.F.
      • Lourenco D.A.L.
      • Silva F.F.
      • Jamrozik J.
      • Schaeffer L.R.
      • Schenkel F.S.
      Invited review: Advances and applications of random regression models: From quantitative genetics to genomics.
      ;
      • Schaeffer L.R.
      • Jamrozik J.
      • Kistemaker G.J.
      • Van Doormaal J.
      Experience with a test-day model.
      ), in addition to more precisely identifying external factors that may be affecting the animal production in different stages of lactation (
      • Jensen R.G.
      The composition of bovine milk lipids: January 1995 to December 2000.
      ). With the advancement of genomics and the availability of high-density SNP panels, it became possible to predict more accurate genomic EBV (GEBV) for young animals compared with the traditional parent average (
      • Meuwissen T.H.
      • Hayes B.J.
      • Goddard M.E.
      Prediction of total genetic value using genome-wide dense marker maps.
      ;
      • Hayes B.J.
      • Visscher P.M.
      • Goddard M.E.
      Increased accuracy of artificial selection by using the realized relationship matrix.
      ). Therefore, especially for sex-limited (e.g., milk production–related traits) and expensive-to-measure traits (e.g., direct individual FA), predicting more accurate GEBV at a young age using relatives' phenotypes predicted based on MIR is a very promising approach. In this context, we aimed to (1) estimate genetic parameters for 5 milk FA groups (i.e., short-chain, medium-chain, long-chain, SFA, and UFA) predicted based on milk spectra of first-parity Holstein cows; (2) perform genomic predictions over time for these 5 traits; and (3) conduct single-step GWAS to identify genomic regions and candidate genes associated with milk FA, and consequently, understand the underlying biology of these traits based on candidates genes and metabolic pathways of FA predicted by MIR.

      MATERIALS AND METHODS

      No animal care committee approval was necessary for the purposes of this study, as all information required was obtained from pre-existing databases.

      Phenotypic, Genotypic, and Pedigree Data

      The pedigree, phenotypic, and genomic data sets were provided by the Canadian Dairy Network (Guelph, ON, Canada). Milk test-day records used in this study were collected between January 2013 and July 2018, from Holstein cows in first lactation. Milk FA composition was predicted for specific days (i.e., test-days in which milk was recorded) from MIR spectra obtained from the routine milk recording systems, using the calibration equations developed by
      • Fleming A.
      • Schenkel F.S.
      • Chen J.
      • Malchiodi F.
      • Bonfatti V.
      • Ali R.A.
      • Mallard B.
      • Corredig M.
      • Miglior F.
      Prediction of milk fatty acid content with mid-infrared spectroscopy in Canadian dairy cattle using differently distributed model development sets.
      as well as the same instrument for spectra acquisition. The coefficients of determination of cross-validation (Rcv2) for the prediction equations were 0.93, 0.83, 0.73, 0.89, and 0.81 for SFA, UFA, SCFA, MCFA, and LCFA groups, respectively (
      • Fleming A.
      • Schenkel F.S.
      • Chen J.
      • Malchiodi F.
      • Bonfatti V.
      • Ali R.A.
      • Mallard B.
      • Corredig M.
      • Miglior F.
      Prediction of milk fatty acid content with mid-infrared spectroscopy in Canadian dairy cattle using differently distributed model development sets.
      ). In brief,
      • Fleming A.
      • Schenkel F.S.
      • Chen J.
      • Malchiodi F.
      • Bonfatti V.
      • Ali R.A.
      • Mallard B.
      • Corredig M.
      • Miglior F.
      Prediction of milk fatty acid content with mid-infrared spectroscopy in Canadian dairy cattle using differently distributed model development sets.
      analyzed the FA profile of individual milk samples using gas chromatography methodologies, and through a partial least squares regression method the authors developed equations to predict FA concentrations from MIR milk spectra. The FA were classified into 5 groups based on the length of the carbon chain and degree of saturation: (1) SCFA (4 to 10 carbons); (2) MCFA (12 to 16 carbons); (3) LCFA (17 to 22 carbons); (4) SFA (no double bonds); and (5) UFA (one or more double bonds).
      Only data from first-parity cows with records from 5 to 305 DIM and calving between 19 and 43 mo of age were kept in the data set for further analyses. The minimum number of records per cow ranged from 2 to 9. Milk FA predictions for samples with MIR spectral data above or below 3 standard deviations from the mean were considered outliers and removed from further analyses. In addition, test-day records assumed as an outlier in at least one FA group were removed from all other FA groups, because it indicates that the FA groups for that specific sample could have been poorly predicted. Herd-test-days with <4 cows, and records above or below 3 standard deviations from the mean, within each herd-test-day, were removed from the analyses. Animals were either genotyped with the Illumina BovineSNP50K BeadChip (Illumina Inc., San Diego, CA) or with a lower-density panel and accurately imputed to 50K using the FImpute software (
      • Sargolzaei M.
      • Chesnais J.P.
      • Schenkel F.S.
      A new approach for efficient genotype imputation using information from relatives.
      ). Genotypic quality control was performed using the preGSf90 software (
      • Aguilar I.
      • Misztal I.
      • Tsuruta S.
      • Legarra A.
      • Wang H.
      PREGSF90–POSTGSF90: Computational tools for the implementation of single-step genomic selection and genome-wide association with ungenotyped individuals in BLUPF90 programs.
      ;
      • Misztal I.
      • Tsuruta S.
      • Lourenco D.A.L.
      • Aguilar I.
      • Legarra A.
      • Vitezica Z.
      Manual for BLUPF90 family of programs.
      ), and SNP with Mendelian conflicts, duplicated or unknown position, call rate <0.95, and minor allele frequency <0.05 were removed. Individual genotypes with call rate <0.95 were also filtered out. After the phenotypic and genotypic quality controls, the data sets included 629,769 records from 201,465 animals in 6,105 herds, and 8,865 animals genotyped with 44,368 SNP. The pedigree file included 10 generations back for all phenotyped animals and contained 7,587,436 individuals.

      Variance Component Estimation

      Variance components for all traits were estimated through Bayesian inference of single-trait RRM under a Markov chain Monte Carlo framework using the GIBBS2F90 software (
      • Misztal I.
      • Tsuruta S.
      • Strabel T.
      • Auvray B.
      • Druet T.
      • Lee D.
      BLUPF90 and related programs (BGF90).
      ). The following RRM was used:
      y = + Zhh + Zaa + Zpp + e,


      where y is the vector of longitudinal observations (FA content; g/100 g of milk); β is the vector of systematic effects, which included herd-test-day and the systematic regression coefficients for age-season of calving effect; h is the vector of random regression polynomial coefficients for herd-year of calving effect; a is the vector of random regression coefficients for animal additive genetic effect; p is the vector of regression coefficients for permanent environmental effect; and e is the residual vector; X, Zh, Za, and Zp are the corresponding incidence matrices for the mentioned effects. Both systematic and random effects were modeled through third-order Legendre orthogonal polynomials (
      • Kirkpatrick M.
      • Lofsvold D.
      • Bulmer M.
      Analysis of the inheritance, selection and evolution of growth trajectories.
      ), which were suggested by
      • Narayana S.G.
      • Schenkel F.S.
      • Fleming A.
      • Koeck A.
      • Malchiodi F.
      • Jamrozik J.
      • Johnston J.
      • Sargolzaei M.
      • Miglior F.
      Genetic analysis of groups of mid-infrared predicted fatty acids in milk.
      using a subset of the same population.
      The data distribution was given by
      y|β,h,a,p,D,T,P,σe2~N(Xβ+Zhh+Zaa+Zpp,Iσe2),


      where D, T, and P are the covariance matrices associated with random regression polynomial coefficients for herd-year of calving, for animal additive genetic, and for permanent environmental effects; I is the identity matrix, and σe2 is the residual variance, assumed as a scaled inverted chi-squared distribution σe2|ve,Se~χ-2(ve,Se), where v is the number of degrees of freedom and S is the scale parameter. The prior distribution for the vector of systematic effects was β~N(0,βI) Σβ being a known diagonal matrix with values 1e+10 (large variances) to represent vague prior knowledge. For the random effects, the following prior distributions were assumed: h~N(0,DI), a~N(0,TA), and p~N(0,PI), I being the identity matrix and A the additive pedigree-based relationship matrix. Furthermore, we assumed that D, T, and P follow the inverted Wishart (IW) distributions, IW(vh, Vh), IW(va, Va), and IW(vp, Vp), respectively, where v is the number of degrees of freedom and V is the scale matrix.
      The length of the Gibbs sampler chain was 270,000 cycles, considering a burn-in period of 30,000 cycles, and a sampling interval (thin) of 15 cycles for all traits. The convergence criterion was based on graphical analysis and Raftery and Lewis criterion (
      • Raftery A.L.
      • Lewis S.
      Comment: One long run with diagnostics: Implementation strategies for Markov chain Monte Carlo.
      ) available in the Bayesian Output Analysis (
      • Smith B.J.
      boa: An R Package for MCMC output convergence assessment and posterior inference.
      ) package of the R software (
      • R Core Team
      R: A language and environment for statistical computing.
      ).
      The posterior marginal distribution samples for heritabilities over DIM were obtained at each Gibbs sampler iteration (q) as follows:
      h2(q)=σg2(q)σa2(q)+σh2(q)+σp2(q)+σe2(q),


      where σg2 is the additive genetic variance, σh2 is the herd-year variance, σp2 is the permanent environmental variance, and σe2 is the residual variance. The h2 estimate for each trait was obtained as the average of daily h2 estimates.

      Genomic Prediction of Breeding Values

      Single-Step Genomic BLUP Approach.

      The same RRM used for the estimation of variance components were used in the single-step genomic BLUP (ssGBLUP) analyses, except that A was replaced by the H matrix. Because the direct estimation of H is computationally demanding, H−1 was calculated as in
      • Aguilar I.
      • Misztal I.
      • Johnson D.
      • Legarra A.
      • Tsuruta S.
      • Lawlor T.
      Hot topic: A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score.
      :
      H-1=A-1+(000τ(0.95G+0.05A22)-1-wA22-1),


      where A−1 is the inverse of A, G−1 is the inverse of the genomic relationship matrix [calculated using the first method proposed by
      • VanRaden P.M.
      Efficient methods to compute genomic predictions.
      and the observed allele frequencies], and A22-1 is the inverse of the portion of A related to genotyped animals. As usual, G−1 and A22-1 are not in the same scale, so scaling factors (τ and ω) were used to make the matrices more compatible. The values reported by
      • Oliveira H.R.
      • Lourenco D.A.L.
      • Masuda Y.
      • Misztal I.
      • Tsuruta S.
      • Jamrozik J.
      • Brito L.F.
      • Silva F.F.
      • Schenkel F.S.
      Application of single-step genomic evaluation using multiple-trait random regression test-day models in dairy cattle.
      for milk fat yield, in the same Holstein population, were used in this study and considered, to differentiate from the default values, to be the optimal values (i.e., τ = 1.5 and ω = 0.6). The ssGBLUP analyses were performed using the BLUPF90 software (
      • Misztal I.
      • Tsuruta S.
      • Lourenco D.A.L.
      • Aguilar I.
      • Legarra A.
      • Vitezica Z.
      Manual for BLUPF90 family of programs.
      ). The variance components previously estimated were assumed as the true variance components for the analyzed population. Animals born between 1957 and 2009 (n = 8,489 genotyped animals and 158,765 animals with phenotypes) and between 2010 and 2012 (n = 376 genotyped animals) were grouped as the training and validation populations, respectively. Thus, the phenotypes of animals born between 2010 and 2012 were excluded from the analyses, to create a reduced data set to be used to predict the GEBV.

      Reliability and Bias of Genomic Predictions.

      The GEBV of the validation animals were used to evaluate the reliability (r2) and bias of genomic predictions for each trait. The reliability of genomic predictions was estimated as the squared Pearson correlation coefficient (r) between daily GEBV, predicted using the reduced data set, and daily EBV predicted using the complete data set (i.e., without excluding phenotypes for the validation animals or their descendants). Using the full-data set EBV instead of pseudo-phenotypes to validate the performance of the genomic prediction relies on the fact that animals in the validation population have reliable breeding values (0.84 ± 0.03). Full-data set EBV have also been used to validate the performance of genomic predictions using RRM in other studies (e.g.,
      • Oliveira H.R.
      • Lourenco D.A.L.
      • Masuda Y.
      • Misztal I.
      • Tsuruta S.
      • Jamrozik J.
      • Brito L.F.
      • Silva F.F.
      • Schenkel F.S.
      Application of single-step genomic evaluation using multiple-trait random regression test-day models in dairy cattle.
      ,
      • Oliveira H.R.
      • Lourenco D.A.L.
      • Masuda Y.
      • Misztal I.
      • Tsuruta S.
      • Jamrozik J.
      • Brito L.F.
      • Silva F.F.
      • Cant J.P.
      • Schenkel F.S.
      Single-step genome-wide association for longitudinal traits of Canadian Ayrshire, Holstein, and Jersey dairy cattle.
      ). The genomic prediction bias was calculated by obtaining the regression coefficient (b1) estimated using a linear regression of EBV (from the complete data sets) on GEBV (from the reduced data sets); that is, EBV = b0 + b1 × GEBV + e. Only animals from the validation population were used to calculate r2 and b1. To analyze the advantages of including genomic information in the genetic evaluation of milk FA, the parent average (PA) was predicted for the validation animals; PA was used to calculate r2 and b1 using daily PA (predicted using the reduced data set) and daily EBV (predicted using the complete data set) for the animals in the validation population.

      GWAS and Functional Analyses.

      Genome-wide association analyses were performed using the postGSf90 software (
      • Aguilar I.
      • Misztal I.
      • Tsuruta S.
      • Legarra A.
      • Wang H.
      PREGSF90–POSTGSF90: Computational tools for the implementation of single-step genomic selection and genome-wide association with ungenotyped individuals in BLUPF90 programs.
      ), considering the variance explained by moving windows of 20 adjacent SNP. The postGSf90 software back-solves the additive genomic random regression coefficients (i.e., the GEBV for the additive random regression coefficients) to SNP effects and it can be described as follows (
      • Wang H.
      • Misztal I.
      • Aguilar I.
      • Legarra A.
      • Fernando R.L.
      • Vitezica Z.
      • Okimoto R.
      • Wing T.
      • Hawken R.
      • Muir W.M.
      Genome-wide association mapping including phenotypes from relatives without genotypes.
      ):
      uˆc=Z[ZZ]-1aˆgc,


      where uˆc is the vector of estimated SNP solutions for the cth random regression coefficient; Z is the matrix containing the centered genotypes; and aˆgc is the vector of GEBV for the cth random regression coefficient, estimated in the ssGBLUP analyses, which contains the cth random regression coefficients for all the genotyped animals. Further, to calculate the effect of each SNP along the lactation, the SNP-specific solutions for all random regression coefficients (c = 1, 2, 3, 4) were combined into a vector (uˆk=[uˆk1uˆk2uˆk3uˆk4]) and used to estimate the SNP effects for all DIM (from 5 to 305), as suggested by
      • Oliveira H.R.
      • Lourenco D.A.L.
      • Masuda Y.
      • Misztal I.
      • Tsuruta S.
      • Jamrozik J.
      • Brito L.F.
      • Silva F.F.
      • Schenkel F.S.
      Application of single-step genomic evaluation using multiple-trait random regression test-day models in dairy cattle.
      :
      SNPˆk=Tuˆk,


      where SNPˆk is the vector that contains the effect of the kth SNP on the trait estimated for each DIM, T is a matrix of covariates for each DIM, associated with the third-order Legendre orthogonal polynomials, and uˆk is the vector of SNP solutions for all random regression coefficients related to the kth SNP. The percentage of genetic variance explained by the ith region was calculated as described by
      • Wang H.
      • Misztal I.
      • Aguilar I.
      • Legarra A.
      • Muir W.M.
      Genome-wide association mapping including phenotypes from relatives without genotypes in a single-step (ssGWAS) for 6-week body weight in broiler chickens.
      :
      Var(ai)σa2×100=Var(j=120Zjuˆj)σa2×100,


      where ai is the genetic value of the ith region that consisted of 20 consecutive SNP, σa2 is the total genetic variance, Zj is the matrix that contains the centered genotypes for the jth SNP, and Ûj is the marker effect of the jth SNP within the ith region.
      A threshold of 1% of the total genetic variance explained by each genomic window was used to define the important genomic regions associated with the traits included in this study. Positional candidate genes were mapped using the Biomart tool (
      • Kinsella R.J.
      • Kahari A.
      • Haider S.
      • Zamora J.
      • Proctor G.
      • Spudich G.
      • Almeida-King J.
      • Staines D.
      • Derwent P.
      • Kerhornou A.
      • Kersey P.
      • Flicek P.
      Ensembl BioMarts: A hub for data retrieval across taxonomic space.
      ) embedded in the Ensembl Genes database version 96 (http://useast.ensembl.org/index.html). Based on the start and end chromosomal positions, important genomic regions were further investigated to understand the biological processes related to the studied traits and to define the most likely functional candidate genes. Complete gene functions were obtained from the National Center for Biotechnology Information database (www.ncbi.nlm.nih.gov/gene/) based on the Bos taurus reference genome ARS-UCD1.2 (https://www.ncbi.nlm.nih.gov/assembly/GCF_002263795.1/). The biological functions and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (
      • Kanehisa M.
      KEGG: Kyoto Encyclopedia of Genes and Genomes.
      ;
      • Kanehisa M.
      • Furumichi M.
      • Tanabe M.
      • Sato Y.
      • Morishima K.
      KEGG: New perspectives on genomes, pathways, diseases and drugs.
      ,
      • Kanehisa M.
      • Sato Y.
      • Furumichi M.
      • Morishima K.
      • Tanabe M.
      New approach for understanding genome variations in KEGG.
      ) in which these genes are involved were assessed using the Database for Annotation, Visualization and Integrated Discovery v6.8. (DAVID, https://david.ncifcrf.gov/home.jsp;
      • Huang D.W.
      • Sherman B.T.
      • Lempicki R.A.
      Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.
      ,
      • Huang D.W.
      • Sherman B.T.
      • Lempicki R.A.
      Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of large gene lists.
      ).

      RESULTS

      Descriptive Statistics

      Table 1 shows the descriptive statistics for the 5 groups of milk FA. Higher proportions of MCFA, followed by LCFA and SCFA (lowest) were observed. On average, the proportions of LCFA, MCFA, SCFA, SFA, and UFA were 33.85, 55.5, 10.60, 73.2, and 26.8%, respectively.
      Table 1Descriptive statistics for the different groups of milk fatty acids in the Holstein breed
      Fatty acid group
      Fatty acid groups were defined as long-chain fatty acids (LCFA; 17 to 22 carbons), medium-chain fatty acids (MCFA; 12 to 16 carbons), short-chain fatty acids (SCFA; 4 to 10 carbons), SFA (no double bond), and UFA (one or more double bonds).
      (g/100 g of milk)
      MeanSDMinimumMaximum
      LCFA1.3310.390.1434.480
      MCFA2.1830.440.5134.497
      SCFA0.4170.090.0630.848
      SFA2.8960.530.7025.689
      UFA1.0650.260.1713.596
      1 Fatty acid groups were defined as long-chain fatty acids (LCFA; 17 to 22 carbons), medium-chain fatty acids (MCFA; 12 to 16 carbons), short-chain fatty acids (SCFA; 4 to 10 carbons), SFA (no double bond), and UFA (one or more double bonds).

      Genetic Parameters

      The average daily h2 estimates for each FA are shown in Table 2, and the complete pattern of daily h2 across the lactation is shown in Figure 1. The average daily h2 ranged from 0.24 (UFA) to 0.47 (SFA and MCFA). In general, h2 observed at the beginning of the lactation was similar for all 5 FA groups (~0.20). For LCFA, the h2 estimates increased until the first third of the lactation, stabilized until the end of the lactation, and finally slightly decreased. For MCFA and SCFA, the h2 pattern was similar, presenting an increase from the middle of the lactation to the end. Following the same pattern of MCFA and SCFA, the h2 curve for SFA increased at the beginning of lactation and continued to increase until reaching the maximum value at the end of lactation. For UFA and LCFA, daily h2 values had a slight increase at the beginning of the lactation and after reaching its peak, decreased until the end of the curve.
      Table 2Mean and SD of daily heritability estimates for 5 groups of milk fatty acids in the Holstein breed
      Fatty acid group
      Fatty acid groups were defined as long-chain fatty acids (LCFA; 17 to 22 carbons), medium-chain fatty acids (MCFA; 12 to 16 carbons), short-chain fatty acids (SCFA; 4 to 10 carbons), SFA (no double bond), and UFA (one or more double bonds).
      MeanSD
      LCFA0.250.07
      MCFA0.470.10
      SCFA0.310.09
      SFA0.470.09
      UFA0.240.07
      1 Fatty acid groups were defined as long-chain fatty acids (LCFA; 17 to 22 carbons), medium-chain fatty acids (MCFA; 12 to 16 carbons), short-chain fatty acids (SCFA; 4 to 10 carbons), SFA (no double bond), and UFA (one or more double bonds).
      Figure thumbnail gr1
      Figure 1Daily heritability estimates of long-chain, medium-chain, short-chain, saturated, and unsaturated milk fatty acid groups from 5 to 305 DIM of first-parity Holstein cows.

      Reliability and Bias of Genomic Predictions

      Average validation r2 and their standard deviations considering both analyses, with and without (i.e., considering the default values) optimal τ and ω, are shown in Table 3. When considering τ and ω equal to 1 (default), the average r2 ranged from 0.564 (LCFA) to 0.737 (SCFA); when using the optimal τ and ω values, the GEBV r2 ranged from 0.583 (LCFA) to 0.732 (SCFA). The inclusion of optimal values of τ and ω slightly changed (in both directions) the average r2, as observed for LCFA and UFA where there was an increase in r2; however, a slight decrease was observed for the other 3 traits. For all traits, the inclusion of genomic information increased the average r2 in both scenarios (with and without the optimal values for τ and ω).
      Table 3Mean validation reliability (and SD) estimated for parent average (r
      SCFA = short-chain fatty acids; MCFA = medium-chain fatty acids; LCFA = long-chain fatty acids.
      PA) and genomic breeding values predicted using (r
      SCFA = short-chain fatty acids; MCFA = medium-chain fatty acids; LCFA = long-chain fatty acids.
      GEBVτω) or not (r
      SCFA = short-chain fatty acids; MCFA = medium-chain fatty acids; LCFA = long-chain fatty acids.
      GEBV) the optimal values
      The optimal values assumed were τ = 1.5 and ω = 0.6.
      for the scaling factors τ and ω, for 5 groups of milk fatty acids
      Trait
      SCFA = short-chain fatty acids; MCFA = medium-chain fatty acids; LCFA = long-chain fatty acids.
      r
      SCFA = short-chain fatty acids; MCFA = medium-chain fatty acids; LCFA = long-chain fatty acids.
      PA
      r
      SCFA = short-chain fatty acids; MCFA = medium-chain fatty acids; LCFA = long-chain fatty acids.
      GEBV
      r
      SCFA = short-chain fatty acids; MCFA = medium-chain fatty acids; LCFA = long-chain fatty acids.
      GEBVτω
      SCFA0.59 (0.02)0.74 (0.02)0.73 (0.01)
      SFA0.55 (0.02)0.70 (0.02)0.68 (0.01)
      MCFA0.55 (0.02)0.70 (0.02)0.68 (0.01)
      UFA0.44 (0.05)0.59 (0.04)0.62 (0.03)
      LCFA0.41 (0.03)0.56 (0.03)0.58 (0.03)
      1 The optimal values assumed were τ = 1.5 and ω = 0.6.
      2 SCFA = short-chain fatty acids; MCFA = medium-chain fatty acids; LCFA = long-chain fatty acids.
      The pattern of GEBV r2 curves over time using default and optimal values for τ and ω are shown in Figure 2. The GEBV r2 were approximately constant over DIM for all groups of FA. Slight differences in the pattern of the GEBV r2 curves were observed when including, or not, τ and ω in the analyses. The average b1 (indicator of bias) and its respective standard deviations are shown in Table 4. The average b1 for PA ranged from 0.690 (UFA) to 1.062 (SCFA). When no optimal values of τ and ω were used (i.e., the default values were applied), the b1 coefficients ranged from 0.756 (LCFA) to 0.905 (SCFA). When the optimal values of τ and ω were used, b1 ranged from 1.075 (LCFA) to 1.172 (SCFA). In summary, an improvement in the bias (i.e., b1 closer to 1) was observed for all traits when optimal τ and ω were used. The pattern of bias over time using the default and the optimal values for τ and ω are shown in Figure 3. Similar pattern of bias over time were found when using, or not, the optimal values of τ and ω. In general, GEBV predicted over time without the optimal values were deflated, and GEBV predicted over time using the optimal values were inflated.
      Figure thumbnail gr2
      Figure 2Daily reliabilities using scaling factors τ and ω equal to 1 (default) and optimal values for long-chain, medium-chain, short-chain, saturated, and unsaturated milk fatty acid groups from 5 to 305 DIM of first-parity Holstein cows
      Table 4Regression coefficients (and SD) estimated for parent average (bPA), and genomic EBV (GEBV) predicted using (bGEBVτω) or not (bGEBV) the optimal values
      The optimal values assumed were τ = 1.5 and ω = 0.6.
      for the scaling factors τ and ω, for 5 groups of milk fatty acids
      Trait
      SCFA = short-chain fatty acids; MCFA = medium-chain fatty acids; LCFA = long-chain fatty acids.
      bPAbGEBVbGEBVτω
      SCFA1.062 (0.025)0.905 (0.046)1.172 (0.028)
      SFA1.030 (0.018)0.883 (0.035)1.140 (0.021)
      MCFA1.054 (0.018)0.759 (0.032)1.164 (0.020)
      UFA0.690 (0.052)0.800 (0.052)1.120 (0.033)
      LCFA0.965 (0.021)0.756 (0.045)1.075 (0.022)
      1 The optimal values assumed were τ = 1.5 and ω = 0.6.
      2 SCFA = short-chain fatty acids; MCFA = medium-chain fatty acids; LCFA = long-chain fatty acids.
      Figure thumbnail gr3
      Figure 3Daily bias (b1) using scaling factors τ and ω equal to 1 (default) and optimal values for long-chain, medium-chain, short-chain, saturated, and unsaturated fatty acid groups from 5 to 305 DIM of first-parity Holstein cows.

      GWAS and Functional Analyses

      A total of 39 genomic windows were associated with the 5 FA groups. These genomic regions were located in chromosomes BTA5, BTA13, and BTA14, and the maximum proportion of the additive genetic variance explained was 15.12% (LCFA) on BTA14. The largest variances explained were found in BTA14, indicating that this chromosome contains important genes associated with milk FA groups. All common candidate genes mapped through Ensembl considering 20-SNP windows for the SNP selected in each trait, using default and optimal τ and ω values, are shown in Table 5. Table 6 shows the different candidate genes discovered between the 2 scenarios (default and optimal) for τ and ω. The biological pathways in which these genes are involved are shown in Table 7. The main biological pathways are carbon metabolism, lipid metabolism, cellular lipid metabolic process, and lipid catabolic process. The Manhattan plots of the proportion of additive genetic variance explained for each third-order Legendre coefficients are shown in Supplemental Figures S1 to S10 (https://doi.org/10.3168/jds.2019-17628). The Manhattan plots show that the peaks of the markers were similar when scaling factors with optimal values or default values (equal to 1) of τ and ω were used.
      Table 5Chromosome information and common candidate gene symbols mapped through Ensembl for SNP considering the 20-SNP windows selected in each trait between default and optimal values of the scaling factors τ and ω
      The default values were τ and ω equal to 1, and the optimal values assumed were τ = 1.5 and ω = 0.6.
      Group
      LCFA = long-chain fatty acids; MCFA = medium-chain fatty acids; SCFA = short-chain fatty acids.
      Common genes
      LCFABTA14: LY6D,
      Genes previously reported to be related to milk fat traits.
      LYNX1, LYPD2, SLURP1, TSNARE1, SLC45A4, PTK2,
      Genes previously reported to be related to milk fat traits.
      AGO2,
      Genes previously reported to be related to milk fat traits.
      ARC, JRK, TRAPPC9,
      Genes previously reported to be related to milk fat traits.
      FAM135B
      Genes previously reported to be related to milk fat traits.
      MCFABTA5: EPS8, RERG, ACO2,
      Genes previously reported to be related to milk fat traits.
      PLBD1 BTA13: ITCH BTA14: LYNX1, SLURP1, TSNARE1, PTK2,
      Genes previously reported to be related to milk fat traits.
      AGO2
      Genes previously reported to be related to milk fat traits.
      SCFABTA5: ERP27, ART4, PLBD1 BTA14: LY6D,
      Genes previously reported to be related to milk fat traits.
      LYNX1, LYPD2, SLURP1, TSNARE1, SLC45A4, PTK2,
      Genes previously reported to be related to milk fat traits.
      AGO2,
      Genes previously reported to be related to milk fat traits.
      ARC, JRK, TRAPPC9,
      Genes previously reported to be related to milk fat traits.
      FAM135B
      Genes previously reported to be related to milk fat traits.
      SFABTA5: RERGL, EPS8, RERG, ARHGDIB, GUCY2C BTA14: LY6D,
      Genes previously reported to be related to milk fat traits.
      LYNX1, LYPD2, SLURP1, TSNARE1, ARC, JRK, SLC45A4, PTK2,
      Genes previously reported to be related to milk fat traits.
      AGO2,
      Genes previously reported to be related to milk fat traits.
      TRAPPC9,
      Genes previously reported to be related to milk fat traits.
      KCNK9
      Genes previously reported to be related to milk fat traits.
      UFABTA14: LY6D,
      Genes previously reported to be related to milk fat traits.
      LYNX1, LYPD2, SLURP1, TSNARE1, SLC45A4, PTK2,
      Genes previously reported to be related to milk fat traits.
      AGO2,
      Genes previously reported to be related to milk fat traits.
      ARC, JRK, TRAPPC9,
      Genes previously reported to be related to milk fat traits.
      KCNK9
      Genes previously reported to be related to milk fat traits.
      1 The default values were τ and ω equal to 1, and the optimal values assumed were τ = 1.5 and ω = 0.6.
      2 LCFA = long-chain fatty acids; MCFA = medium-chain fatty acids; SCFA = short-chain fatty acids.
      3 Genes previously reported to be related to milk fat traits.
      Table 6Chromosome information and different candidate gene symbols mapped through Ensembl for SNP considering the 20-SNP windows selected in each trait between default and optimal values of the scaling factors τ and ω
      The default values were τ and ω equal to 1, and the optimal values assumed were τ = 1.5 and ω = 0.6.
      Group
      LCFA = long-chain fatty acids; MCFA = medium-chain fatty acids; SCFA = short-chain fatty acids.
      Default τ and ωOptimal τ and ω
      LCFABTA14: ARC, JRK, TRAPPC9,
      Genes previously reported to be related to milk fat traits.
      FAM135B
      Genes previously reported to be related to milk fat traits.
      MCFABTA5: RERGL, ARHGDIB, ERP27, ART4, SMCO3 BTA13: ITCH BTA14: LY6D,
      Genes previously reported to be related to milk fat traits.
      LYPD2, SLC45A4, TRAPPC9,
      Genes previously reported to be related to milk fat traits.
      FAM135B
      Genes previously reported to be related to milk fat traits.
      BTA5: GUCY2C BTA14: ARC, JRK, GPR20
      SCFABTA5: RERGL, EPS8, RERG, ARHGDIB, ERP27; BTA13: DIP2C, PRNP, PRNDBTA5: ERP27, WBP11, GUCY2C BTA14: FAM135B
      Genes previously reported to be related to milk fat traits.
      SFABTA5: RERGL, EP27, ART4, SMCO3, PLBD1BTA14: ST3GAL1, NDRG1, SLA
      UFABTA14: FAM135B
      Genes previously reported to be related to milk fat traits.
      1 The default values were τ and ω equal to 1, and the optimal values assumed were τ = 1.5 and ω = 0.6.
      2 LCFA = long-chain fatty acids; MCFA = medium-chain fatty acids; SCFA = short-chain fatty acids.
      3 Genes previously reported to be related to milk fat traits.
      Table 7Biological functions and pathways (Kyoto Encyclopedia of Genes and Genomes database) of selected genes from long-chain, medium-chain, short-chain, saturated, and unsaturated milk fatty acids
      Biological pathways related to milk fat traits.
      Gene symbolGene nameBiological pathway
      ACO2Aconitase 2Citrate cycle (tricarboxylic acid cycle), Glyoxylate and dicarboxylate metabolism, Carbon metabolism,
      Biological pathways related to milk fat traits.
      2-Oxocarboxylic acid metabolism, biosynthesis of amino acids
      ARCActivity regulated cytoskeleton associated proteinCytoskeleton organization, Anterior/posterior pattern specification, Cell migration, Regulation of cell morphogenesis
      ARHGDIBRho GDP dissociation inhibitor betaNeurotrophin signaling pathway, Vasopressin-regulated water reabsorption
      DIP2CDisco interacting protein 2 homolog CLipid metabolism,
      Biological pathways related to milk fat traits.
      Secondary metabolites biosynthesis, transport, and catabolism
      FAM135BFamily with sequence similarity 135 member BCellular lipid metabolic process
      Biological pathways related to milk fat traits.
      ITCHItchy E3 ubiquitin protein ligaseUbiquitin mediated proteolysis, Endocytosis, tumor necrosis factor signaling pathway, Non-alcoholic fatty liver disease (NAFLD)
      KCNK9Potassium two pore domain channel subfamily k member 9Aldosterone synthesis and secretion
      PLBD1Phospholipase B domain containing 1Lipid catabolic process,
      Biological pathways related to milk fat traits.
      Phospholipase B-like
      PRNPPrion proteinPrion diseases
      PTK2Protein tyrosine kinase 2ErbB signaling pathway, Chemokine signaling pathway, PI3K-Akt signaling pathway, Axon guidance, Vascular endothelial growth factor (VEGF) signaling pathway, Focal adhesion, Leukocyte transendothelial migration, Regulation of actin cytoskeleton, Bacterial invasion of epithelial cells
      1 Biological pathways related to milk fat traits.
      A total of 444 SNP were selected as the most relevant (top 1%) for each trait. When considering τ and ω equal to 1 (default), the proportion of variance explained by the 20-SNP window ranged from 0.72% (UFA) to 15.12% (LCFA); when using the optimal τ and ω values, the proportion of variance explained ranged from 0.71% (SFA) to 15.10% (LCFA). In relation to SNP effects over DIM, when considering default values for τ and ω, it ranged (in module) from 0.57% (SCFA) to 14.96% (SFA). When using the optimal τ and ω values, the SNP effect over DIM, it ranged from 0.59% (SCFA) to 14.97% (SFA). Both the largest SNP effects and the genomic windows explaining the highest proportion of the additive genetic variance were found for the LCFA group. The trajectory of SNP effects over DIM, estimated for the top 10 SNP of each trait considering τ and ω equal to 1 and optimal values, are shown in Supplemental Figure S11 (https://doi.org/10.3168/jds.2019-17628) and Figure 4, respectively. The chromosome and position of the top 10 SNP are shown in Supplemental Table S1 (https://doi.org/10.3168/jds.2019-17628). The pattern of the SNP effects over lactation for LCFA, SCFA, and UFA tended to remain constant, with small changes in magnitude. Higher deviations in the pattern of SNP effects were observed for MCFA and SFA, where the effects of SNP increased as lactation progressed.
      Figure thumbnail gr4
      Figure 4Single nucleotide polymorphism effect patterns of the top 10 SNP over DIM for long-chain (LCFA), medium-chain (MCFA), short-chain (SCFA), saturated (SFA), and unsaturated (UFA) milk fatty acids using optimal values of scaling factors τ and ω. SNPSol = effect of the SNP.

      DISCUSSION

      Descriptive Statistics and Genetic Parameters

      Current recommendations are for consumers to reduce SFA intake and consume more UFA (
      • FAO (Food and Agriculture Organization of the United Nations)
      Fats and fatty acids in human nutrition. Report of an expert consultation. 91.
      ). Coupled with an appropriate feeding system, a reduction in the proportion of SFA in milk and an increase in the proportion of FA that promote good health, may therefore be of great interest in dairy cattle (
      • Fleming A.
      Phenotypic and genetic variation of milk fat components incorporating mid-infrared technology.
      ). Previous authors have indicated that higher prediction model accuracies were obtained when defining the FA measurement unit on a per-milk basis versus a per-fat basis (
      • Soyeurt H.
      • Dardenne P.
      • Dehareng F.
      • Lognay G.
      • Veselko D.
      • Marlier M.
      • Bertozzi C.
      • Mayeres P.
      • Gengler N.
      Estimating fatty acid content in cow milk using mid-infrared spectrometry.
      ;
      • Rutten M.J.M.
      • Bovenhuis H.
      • Hettinga K.A.
      • van Valenberg H.J.F.
      • van Arendonk J.A.M.
      Predicting bovine milk fat composition using infrared spectroscopy based on milk samples collected in winter and summer.
      ;
      • Fleming A.
      • Schenkel F.S.
      • Chen J.
      • Malchiodi F.
      • Bonfatti V.
      • Ali R.A.
      • Mallard B.
      • Corredig M.
      • Miglior F.
      Prediction of milk fatty acid content with mid-infrared spectroscopy in Canadian dairy cattle using differently distributed model development sets.
      ). This might be due to the fact that milk samples contain different amounts of total fat, and samples with the same relative concentrations of FA can contain very different total quantities of the FA (
      • Fleming A.
      Phenotypic and genetic variation of milk fat components incorporating mid-infrared technology.
      ). Therefore, the unit used in this study was grams per 100 g of milk.
      • Fleming A.
      • Schenkel F.S.
      • Chen J.
      • Malchiodi F.
      • Bonfatti V.
      • Ali R.A.
      • Mallard B.
      • Corredig M.
      • Miglior F.
      Prediction of milk fatty acid content with mid-infrared spectroscopy in Canadian dairy cattle using differently distributed model development sets.
      reported mean values of FA based on GC analysis of 2.81, 1.10, 0.464, 2.09, and 1.40 g/dL of milk for SFA, UFA, SCFA, MCFA, and LCFA, respectively. These estimates are in a similar range to the averages (based on MIR predictions) observed in this study. Milk FA are highly correlated with milk fat percentage and fat yield (e.g.,
      • Fleming A.
      • Schenkel F.S.
      • Malchiodi F.
      • Ali R.A.
      • Mallard B.
      • Sargolzaei M.
      • Jamrozik J.
      • Johnston J.
      • Miglior F.
      Genetic correlations of mid-infrared-predicted milk fatty acid groups with milk production traits.
      ). If the breeding goal is to alter a specific FA profile, then expressing the traits on a grams per 100 g of fat basis could be a good approach. However, we might be interested in increasing the amount of certain FA in the milk (independently of fat yield). Because predictions based on grams per 100 g of fat are less accurate (e.g.,
      • Soyeurt H.
      • Dardenne P.
      • Dehareng F.
      • Lognay G.
      • Veselko D.
      • Marlier M.
      • Bertozzi C.
      • Mayeres P.
      • Gengler N.
      Estimating fatty acid content in cow milk using mid-infrared spectrometry.
      ;
      • Rutten M.J.M.
      • Bovenhuis H.
      • Hettinga K.A.
      • van Valenberg H.J.F.
      • van Arendonk J.A.M.
      Predicting bovine milk fat composition using infrared spectroscopy based on milk samples collected in winter and summer.
      ;
      • Fleming A.
      • Schenkel F.S.
      • Chen J.
      • Malchiodi F.
      • Bonfatti V.
      • Ali R.A.
      • Mallard B.
      • Corredig M.
      • Miglior F.
      Prediction of milk fatty acid content with mid-infrared spectroscopy in Canadian dairy cattle using differently distributed model development sets.
      ), this could greatly limit their use in breeding programs. Therefore, we used grams per 100 g of milk in this study.
      The heritability estimates obtained for all 5 FA groups indicated that they are under moderate genetic control; therefore, it is genetically possible to manipulate milk FA proportions. The higher h2 estimates observed for de novo synthesized SCFA compared with dietary- and adipose-derived LCFA was expected. Our results were mostly in agreement with the literature (
      • Bastin C.
      • Gengler N.
      • Soyeurt H.
      Phenotypic and genetic variability of production traits and milk fatty acid contents across days in milk for Walloon Holstein first-parity cows.
      ;
      • Narayana S.G.
      • Schenkel F.S.
      • Fleming A.
      • Koeck A.
      • Malchiodi F.
      • Jamrozik J.
      • Johnston J.
      • Sargolzaei M.
      • Miglior F.
      Genetic analysis of groups of mid-infrared predicted fatty acids in milk.
      ;
      • Hein L.
      • Sørensen L.P.
      • Kargo M.
      • Buitenhuis A.J.
      Genetic analysis of predicted fatty acid profiles of milk from Danish Holstein and Danish Jersey cattle populations.
      ). The h2 estimates were similar to those found by
      • Narayana S.G.
      • Schenkel F.S.
      • Fleming A.
      • Koeck A.
      • Malchiodi F.
      • Jamrozik J.
      • Johnston J.
      • Sargolzaei M.
      • Miglior F.
      Genetic analysis of groups of mid-infrared predicted fatty acids in milk.
      , using data from the same population of first-parity Canadian Holsteins (data subset) and FA MIR prediction equations for LCFA and UFA (0.23 and 0.21, respectively), but were slightly higher for MCFA (0.32), SCFA (0.24) and SFA (0.33). In contrast, the values for these FA groups were similar to those reported by
      • Bastin C.
      • Gengler N.
      • Soyeurt H.
      Phenotypic and genetic variability of production traits and milk fatty acid contents across days in milk for Walloon Holstein first-parity cows.
      in Holstein cattle from the Walloon Region of Belgium, who found h2 of 0.44 and 0.43 for MCFA and SFA, respectively. These differences may be because the present study had more observations than the aforementioned and the distinct populations used in each study. In addition, the accuracy of the calibration equations can affect the results, as reported by
      • Fleming A.
      • Schenkel F.S.
      • Chen J.
      • Malchiodi F.
      • Bonfatti V.
      • Ali R.A.
      • Mallard B.
      • Corredig M.
      • Miglior F.
      Prediction of milk fatty acid content with mid-infrared spectroscopy in Canadian dairy cattle using differently distributed model development sets.
      . In the current data set, a limited number of test-days were available for some animals because MIR spectra was recorded for only a proportion of all milk samples during DHI milk recording. Consequently, there were fewer records per cow and the traits may not have been as well modeled along one lactation compared with
      • Bastin C.
      • Gengler N.
      • Soyeurt H.
      Phenotypic and genetic variability of production traits and milk fatty acid contents across days in milk for Walloon Holstein first-parity cows.
      .

      Reliability and Bias of Genomic Predictions

      The τ and ω parameters are used to make G−1 more compatible with A22-1 and A−1 (
      • Misztal I.
      • Bradford H.L.
      • Lourenco D.A.L.
      • Tsuruta S.
      • Masuda Y.
      • Legarra A.
      Studies on inflation of GEBV in single-step GBLUP for type.
      ;
      • Oliveira H.R.
      • Lourenco D.A.L.
      • Masuda Y.
      • Misztal I.
      • Tsuruta S.
      • Jamrozik J.
      • Brito L.F.
      • Silva F.F.
      • Schenkel F.S.
      Application of single-step genomic evaluation using multiple-trait random regression test-day models in dairy cattle.
      ). In this context,
      • Misztal I.
      • Bradford H.L.
      • Lourenco D.A.L.
      • Tsuruta S.
      • Masuda Y.
      • Legarra A.
      Studies on inflation of GEBV in single-step GBLUP for type.
      showed that τ is used to account for the reduced genetic variance, whereas ω is used to account for different depths of pedigree. For all traits, the GEBV r2 were higher than the GEBV r2 calculated for PA, even when no scaling factors were used. These results are in agreement with other studies, such as
      • Oliveira H.R.
      • Lourenco D.A.L.
      • Masuda Y.
      • Misztal I.
      • Tsuruta S.
      • Jamrozik J.
      • Brito L.F.
      • Silva F.F.
      • Schenkel F.S.
      Application of single-step genomic evaluation using multiple-trait random regression test-day models in dairy cattle.
      and
      • Baba T.
      • Gotoh Y.
      • Yamaguchi S.
      • Nakagawa S.
      • Abe H.
      • Masuda Y.
      • Kawahara T.
      Application of single-step genomic best linear unbiased prediction with a multiple-lactation random regression test-day model for Japanese Holsteins.
      , which showed that the use of genomic information through ssGBLUP increases the r2 of breeding values. It is also worth noting that the animals used for these analyses have highly accurate EBV for all traits (and consequently highly accurate PA).
      With regards to b1, the use of optimal τ and ω scaling factors yielded less biased GEBV (i.e., b1 closer to 1) compared with analyses using the default values (τ and ω = 1). In general, GEBV predicted over time with and without the optimal values were inflated and deflated, respectively. However, GEBV estimated for the majority of traits in this study would be accepted for use in the international comparisons performed by Interbull (
      • Interbull
      Genetic evaluation.
      ), if estimated with or without the scaling factors. Compared with PA, GEBV were more biased for the majority of traits. In this context,
      • Lourenco D.A.L.
      • Fragomeni B.O.
      • Bradford H.L.
      • Menezes I.R.
      • Ferraz J.B.S.
      • Aguilar I.
      • Tsuruta S.
      • Misztal I.
      Implications of SNP weighting on single-step genomic predictions for different reference population sizes.
      showed that the use of a genomic relationship matrix weighting markers according to the analyzed trait can generate more reliable and less biased estimates. Thus, posterior studies using the weighted ssGBLUP are needed.
      Among the advantages of the method proposed in this study to evaluate FA is the fact that GEBV and PA are predicted for all DIM. Thus, different selection criteria can be defined in the breeding schemes. In this regard, it becomes possible to change the pattern of the genetic (and phenotypic) FA lactation curves. Increase in GEBV r2 can be achieved by increasing the numbers of both genotyped and phenotyped animals and the use of more accurate FA prediction equations (developed using a larger number of animals with gold standard FA measurements). Different τ and ω values may affect the GEBV obtained and, therefore, may likely change the SNP effects. Here, optimal τ and ω values were used that maximized the observed accuracy of the GEBVs for milk fat yield in the same Holstein population.

      GWAS and Functional Analyses

      This study supports the polygenic nature of the studied traits, with many regions across the genome making small contributions to total genetic variation. The coefficients of the Legendre polynomial, despite presenting differences in the patterns of markers identified, have no biological explanation, because they are used to model the lactation curve. Thus, some coefficients allow identification of different markers at different points in the curve. Important contributions were found on BTA14. This chromosome is well known to have a substantial number of QTL for milk and fat production traits (
      • Viitala S.M.
      • Schulman N.F.
      • de Koning D.J.
      • Elo K.
      • Kinos R.
      • Virta A.
      • Virta J.
      • Mäki-Tanila A.
      • Vilkki J.H.
      Quantitative trait loci affecting milk production traits in Finnish Ayrshire dairy cattle.
      ;
      • Wang H.
      • Jiang L.
      • Liu X.
      • Yang J.
      • Wei J.
      • Xu J.
      • Zhang Q.
      • Liu J.-F.
      A post-GWAS replication study confirming the PTK2 gene associated with milk production traits in Chinese Holstein.
      ), including DGAT1, a gene involved in the last step of the synthesis of triacylglycerol and that has a major effect on milk fat content (
      • Grisart B.
      Positional candidate cloning of a QTL in dairy cattle: Identification of a missense mutation in the bovine DGAT1 gene with major effect on milk yield and composition.
      ;
      • Cruz V.A.R.
      • Oliveira H.R.
      • Brito L.F.
      • Fleming A.
      • Larmer S.
      • Miglior F.
      • Schenkel F.S.
      Genome-wide association study for milk fatty acids in Holstein cattle accounting for the DGAT1 gene effect.
      ). Despite its importance, DGAT1 was not found in this study because of the absence of markers located in its proximity. In addition, DGAT1 has changed position in the new cattle reference genome (ARS-UCD1.2): from BTA14:1,795,425–1,804,838 Mb (Bos_taurus_UMD_3.1.1) to BTA14:603,981–614,153 Mb (ARS-UCD1.2). However, the peak observed at the beginning of BTA14 is likely in high linkage disequilibrium with markers located in DGAT1. In fact, conducting a GWAS for milk FA in the same population,
      • Cruz V.A.R.
      • Oliveira H.R.
      • Brito L.F.
      • Fleming A.
      • Larmer S.
      • Miglior F.
      • Schenkel F.S.
      Genome-wide association study for milk fatty acids in Holstein cattle accounting for the DGAT1 gene effect.
      identified DGAT1 when using the previous reference genome.
      Candidate genes were identified and shown to have important biological functions associated with the traits under study, including lipid metabolism and cholesterol synthesis. For instance, PTK2, associated with all FA, is related to several signal transduction pathways, such as cell motility and microtubule stability (
      • Harte M.T.
      • Hildebrand J.D.
      • Burnham M.R.
      • Bouton A.H.
      • Parsons J.T.
      p130Cas, a substrate associated with v-Src and v-Crk, localizes to focal adhesions and binds to focal adhesion kinase.
      ;
      • Ezratty E.J.
      • Partridge M.A.
      • Gundersen G.G.
      Microtubule-induced focal adhesion disassembly is mediated by dynamin and focal adhesion kinase.
      ), and an important trait statistic related to fat percentage in Chinese Holstein cattle (
      • Wang H.
      • Jiang L.
      • Liu X.
      • Yang J.
      • Wei J.
      • Xu J.
      • Zhang Q.
      • Liu J.-F.
      A post-GWAS replication study confirming the PTK2 gene associated with milk production traits in Chinese Holstein.
      ). Similarly to PTK2, TRAPPC9 is associated with all the FA traits.
      • Wang X.
      • Ma P.
      • Liu J.
      • Zhang Q.
      • Zhang Y.
      • Ding X.
      • Jiang L.
      • Wang Y.
      • Zhang Y.
      • Sun D.
      • Zhang S.
      • Su G.
      • Yu Y.
      Genome-wide association study in Chinese Holstein cows reveal two candidate genes for somatic cell score as an indicator for mastitis susceptibility.
      reported that TRAPPC9 plays a role in regulation of cell differentiation and is associated with mastitis susceptibility in Chinese Holstein cattle. Moreover, ACO2 is related to the second step of the tricarboxylic acid cycle (
      • Palombo V.
      • Milanesi M.
      • Sgorlon S.
      • Capomaccio S.
      • Mele M.
      • Nicolazzi E.
      • Ajmone-Marsan P.
      • Pilla F.
      • Stefanon B.
      • D'Andrea M.
      Genome-wide association study of milk fatty acid composition in Italian Simmental and Italian Holstein cows using single nucleotide polymorphism arrays.
      ). This gene encodes a protein that acts in the mitochondrion to catalyze the interconversion of citrate to isocitrate via cis-aconitate.
      Some genes were found to be related to all FA groups. For instance, SLURP1 was reported by
      • Zhao A.
      • Jiang F.
      • Yang G.
      • Liu H.
      • Li L.
      Sfrp5 interacts with Slurp1 to regulate the accumulation of triglycerides in hepatocyte steatosis model.
      , who suggested that this gene may play an important role in the regulation of triglycerides through secreted frizzled-related protein 5, an anti-inflammatory adipocytokine secreted by adipocytes in mices.
      • Kim B.S.
      • Jung J.S.
      • Jang J.H.
      • Kang K.S.
      • Kang S.K.
      Nuclear Argonaute 2 regulates adipose tissue-derived stem cell survival through direct control of miR10b and selenoprotein N1 expression.
      reported that AGO2 has a leading function in microRNA-induced RNA silencing, which is critical for stem cell self-renewal, development, and other functions in humans. Also related to all traits, TSNARE1 is reported to be associated with inflammatory processes in humans (
      • Ligthart S.
      • Marzi C.
      • Aslibekyan S.
      • Mendelson M.M.
      • Conneely K.N.
      • Tanaka T.
      • Colicino E.
      • Waite L.L.
      • Joehanes R.
      • Guan W.
      • Brody J.A.
      • Elks C.
      • Marioni R.
      • Jhun M.A.
      • Agha G.
      • Bressler J.
      • Ward-Caviness C.K.
      • Chen B.H.
      • Huan T.
      • Bakulski K.
      • Salfati E.L.
      • Fiorito G.
      • Wahl S.
      • Schramm K.
      • Sha J.
      • Hernandez D.G.
      • Just A.C.
      • Smith J.A.
      • Sotoodehnia N.
      • Pilling L.C.
      • Pankow J.S.
      • Tsao P.S.
      • Liu C.
      • Zhao W.
      • Guarrera S.
      • Michopoulos V.J.
      • Smith A.K.
      • Peters M.J.
      • Melzer D.
      • Vokonas P.
      • Fornage M.
      • Prokisch H.
      • Bis J.C.
      • Chu A.Y.
      • Herder C.
      • Grallert H.
      • Yao C.
      • Shah S.
      • McRae A.F.
      • Lin H.
      • Horvath S.
      • Fallin D.
      • Hofman A.
      • Wareham N.J.
      • Wiggins K.L.
      • Feinberg A.P.
      • Starr J.M.
      • Visscher P.M.
      • Murabito J.M.
      • Kardia S.L.R.
      • Absher D.M.
      • Binder E.B.
      • Singleton A.B.
      • Bandinelli S.
      • Peters A.
      • Waldenberger M.
      • Matullo G.
      • Schwartz J.D.
      • Demerath E.W.
      • Uitterlinden A.G.
      • van Meurs J.B.J.
      • Franco O.H.
      • Chen Y.-D. I.
      • Levy D.
      • Turner S.T.
      • Deary I.J.
      • Ressler K.J.
      • Dupuis J.
      • Ferrucci L.
      • Ong K.K.
      • Assimes T.L.
      • Boerwinkle E.
      • Koenig W.
      • Arnett D.K.
      • Baccarelli A.A.
      • Benjamin E.J.
      • Dehghan A.
      DNA methylation signatures of chronic low-grade inflammation are associated with complex diseases.
      ;
      • Arpón A.
      • Milagro F.I.
      • Ramos-Lopez O.
      • Mansego M.L.
      • Riezu-Boj J.-I.
      • Martínez J.A.
      Methylome-wide association study in peripheral white blood cells focusing on central obesity and inflammation.
      ).
      Biological pathways related to lipid catabolism and metabolism have been identified and present important processes in the animal body. The pathway “Lipid catabolism” has FA as a direct product. The FA are catabolized in the β-oxidation process, sequentially removing 2-carbon acetyl groups from the fatty acid chain ends, reducing nicotinamide adenine dinucleotide (NAD+) and flavin adenine dinucleotide (FAD) to produce NADH and FADH2, respectively, and these electrons can be used in energy production by the oxidative phosphorylation process (
      • Stryer L.
      Fatty acid metabolism.
      ). Another important biological pathway identified is related to “Lipid metabolism.”
      • Keller H.
      • Dreyer C.
      • Medin J.
      • Mahfoudi A.
      • Ozato K.
      • Wahli W.
      Fatty acids and retinoids control lipid metabolism through activation of peroxisome proliferator-activated receptor-retinoid X receptor heterodimers.
      showed that FA, especially polyunsaturated FA, have great potential in activating the peroxisome proliferator-activated receptors, the receptors responsible for regulating β-oxidation of FA, thus preventing them from accumulating at lethal levels in blood.

      CONCLUSIONS

      All 5 groups of MIR-predicted milk fatty acids (SCFA, MCFA, LCFA, SFA, and UFA) are under moderate genetic control in the North American Holstein breed, with heritability estimates varying throughout lactation. Genetic improvement can be made in all of these traits through selection. Moderately to highly accurate genomic breeding values can be obtained for daily points throughout lactation using ssGBLUP. Inclusion of scaling factors improved the accuracy for MCFA, UFA, and LCFA. Various genomic regions are associated with the traits, but many regions across the genome have small individual effects on the total genetic variation of each trait. Candidate genes and pathways were identified, which helps in our understanding of the genetic background of milk fatty acid composition. These results contribute to a better understanding of the genetic architecture of predicted milk fatty acid in dairy cattle and will be of great value for the implementation of genomic selection for milk fatty acid composition.

      ACKNOWLEDGMENTS

      The first author acknowledges the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES; Brasília, DF, Brazil). Funding for the collection of milk spectra was provided by Agriculture and Agri-Food Canada and by additional contributions from Dairy Farmers of Canada, the Canadian Dairy Network (Ottawa, ON, Canada), and the Canadian Dairy Commission under the Agri-Science Clusters Initiative (Dairy Farmers of Canada, Agriculture and Agri-Food Canada, the Canadian Dairy Network, and the Canadian Dairy Commission). The authors have not stated any conflicts of interest.

      REFERENCES

        • Aguilar I.
        • Misztal I.
        • Johnson D.
        • Legarra A.
        • Tsuruta S.
        • Lawlor T.
        Hot topic: A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score.
        J. Dairy Sci. 2010; 93 (20105546): 743-752
        • Aguilar I.
        • Misztal I.
        • Tsuruta S.
        • Legarra A.
        • Wang H.
        PREGSF90–POSTGSF90: Computational tools for the implementation of single-step genomic selection and genome-wide association with ungenotyped individuals in BLUPF90 programs.
        in: Proc. 10th World Congr. Genet. Appl. Livest. Prod., Vancouver, BC, Canada. American Society of Animal Science, Champaign, IL2014: 1-3
        • Arpón A.
        • Milagro F.I.
        • Ramos-Lopez O.
        • Mansego M.L.
        • Riezu-Boj J.-I.
        • Martínez J.A.
        Methylome-wide association study in peripheral white blood cells focusing on central obesity and inflammation.
        Genes (Basel). 2019; 10 (31212707): 444
        • Baba T.
        • Gotoh Y.
        • Yamaguchi S.
        • Nakagawa S.
        • Abe H.
        • Masuda Y.
        • Kawahara T.
        Application of single-step genomic best linear unbiased prediction with a multiple-lactation random regression test-day model for Japanese Holsteins.
        Anim. Sci. J. 2017; 88 (27925408): 1226-1231
        • Baer R.J.
        Alteration of the fatty acid content of milk fat.
        J. Food Prot. 1991; 54: 383-386
        • Bastin C.
        • Gengler N.
        • Soyeurt H.
        Phenotypic and genetic variability of production traits and milk fatty acid contents across days in milk for Walloon Holstein first-parity cows.
        J. Dairy Sci. 2011; 94 (21787950): 4152-4163
        • Bergamaschi M.
        • Stocco G.
        • Valorz C.
        • Bazzoli I.
        • Sturaro E.
        • Ramanzin M.
        Cheesemaking in highland pastures: Milk technological properties, cream, cheese and ricotta yields, milk nutrients recovery, and products composition.
        J. Dairy Sci. 2016; 99: 9631-9646
        • Bilal G.
        • Cue R.I.
        • Mustafa A.F.
        • Hayes J.F.
        Short communication: Genetic parameters of individual fatty acids in milk of Canadian Holsteins.
        J. Dairy Sci. 2014; 97 (24290819): 1150-1156
        • Briggs M.A.
        • Petersen K.S.
        • Kris-Etherton P.M.
        Saturated fatty acids and cardiovascular disease: Replacements for saturated fat to reduce cardiovascular risk.
        Healthcare (Basel). 2017; 5: 29
        • Cozma A.
        • Miere D.
        • Filip L.
        • Banc R.
        • Loghin F.
        A review of the metabolic origins of milk fatty acids.
        Not. Sci. Biol. 2013; 5: 270-274
        • Cruz V.A.R.
        • Oliveira H.R.
        • Brito L.F.
        • Fleming A.
        • Larmer S.
        • Miglior F.
        • Schenkel F.S.
        Genome-wide association study for milk fatty acids in Holstein cattle accounting for the DGAT1 gene effect.
        Animals (Basel). 2019; 9 (31752271): 997
        • Ezratty E.J.
        • Partridge M.A.
        • Gundersen G.G.
        Microtubule-induced focal adhesion disassembly is mediated by dynamin and focal adhesion kinase.
        Nat. Cell Biol. 2005; 7 (15895076): 581-590
        • Fleming A.
        Phenotypic and genetic variation of milk fat components incorporating mid-infrared technology.
        (PhD Thesis) Department of Animal Biosciences, University of Guelph, Guelph, ON, Canada2016
        • Fleming A.
        • Schenkel F.S.
        • Chen J.
        • Malchiodi F.
        • Bonfatti V.
        • Ali R.A.
        • Mallard B.
        • Corredig M.
        • Miglior F.
        Prediction of milk fatty acid content with mid-infrared spectroscopy in Canadian dairy cattle using differently distributed model development sets.
        J. Dairy Sci. 2017; 100 (28434722): 5073-5081
        • Fleming A.
        • Schenkel F.S.
        • Malchiodi F.
        • Ali R.A.
        • Mallard B.
        • Sargolzaei M.
        • Jamrozik J.
        • Johnston J.
        • Miglior F.
        Genetic correlations of mid-infrared-predicted milk fatty acid groups with milk production traits.
        J. Dairy Sci. 2018; 101 (29477537): 4295-4306
        • FAO (Food and Agriculture Organization of the United Nations)
        Fats and fatty acids in human nutrition. Report of an expert consultation. 91.
        FAO, Rome, Italy2010
        • Freitas P.H.F.
        • Oliveira H.R.
        • Silva F.F.
        • Fleming A.
        • Schenkel F.S.
        • Miglior F.
        • Brito L.F.
        Short communication: Time-dependent genetic parameters and single-step genome-wide association analyses for predicted milk fatty acid composition in Ayrshire and Jersey dairy cattle.
        J. Dairy Sci. 2020; 103: 5263-5269
        • Grisart B.
        Positional candidate cloning of a QTL in dairy cattle: Identification of a missense mutation in the bovine DGAT1 gene with major effect on milk yield and composition.
        Genome Res. 2002; 12 (11827942): 222-231
        • Harte M.T.
        • Hildebrand J.D.
        • Burnham M.R.
        • Bouton A.H.
        • Parsons J.T.
        p130Cas, a substrate associated with v-Src and v-Crk, localizes to focal adhesions and binds to focal adhesion kinase.
        J. Biol. Chem. 1996; 271 (8662921): 13649-13655
        • Haug A.
        • Høstmark A.T.
        • Harstad O.M.
        Bovine milk in human nutrition—A review.
        Lipids Health Dis. 2007; 6 (17894873): 25
        • Hayes B.J.
        • Visscher P.M.
        • Goddard M.E.
        Increased accuracy of artificial selection by using the realized relationship matrix.
        Genet. Res. (Camb.). 2009; 91 (19220931): 47-60
        • Hein L.
        • Sørensen L.P.
        • Kargo M.
        • Buitenhuis A.J.
        Genetic analysis of predicted fatty acid profiles of milk from Danish Holstein and Danish Jersey cattle populations.
        J. Dairy Sci. 2018; 101 (29248226): 2148-2157
        • Huang D.W.
        • Sherman B.T.
        • Lempicki R.A.
        Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.
        Nat. Protoc. 2009; 4 (19131956): 44-57
        • Huang D.W.
        • Sherman B.T.
        • Lempicki R.A.
        Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of large gene lists.
        Nucleic Acids Res. 2009; 37 (19033363): 1-13
        • Interbull
        Genetic evaluation.
        • Jensen R.G.
        The composition of bovine milk lipids: January 1995 to December 2000.
        J. Dairy Sci. 2002; 85 (11913692): 295-350
        • Kanehisa M.
        KEGG: Kyoto Encyclopedia of Genes and Genomes.
        Nucleic Acids Res. 2000; 28: 27-30
        • Kanehisa M.
        • Furumichi M.
        • Tanabe M.
        • Sato Y.
        • Morishima K.
        KEGG: New perspectives on genomes, pathways, diseases and drugs.
        Nucleic Acids Res. 2017; 45 (27899662): D353-D361
        • Kanehisa M.
        • Sato Y.
        • Furumichi M.
        • Morishima K.
        • Tanabe M.
        New approach for understanding genome variations in KEGG.
        Nucleic Acids Res. 2019; 47 (30321428): D590-D595
        • Keller H.
        • Dreyer C.
        • Medin J.
        • Mahfoudi A.
        • Ozato K.
        • Wahli W.
        Fatty acids and retinoids control lipid metabolism through activation of peroxisome proliferator-activated receptor-retinoid X receptor heterodimers.
        Proc. Natl. Acad. Sci. USA. 1993; 90 (8384714): 2160-2164
        • Kęsek M.
        • Szulc T.
        • Zielak-Steciwko A.
        Genetic, physiological and nutritive factors affecting the fatty acid profile in cows' milk—A review.
        Anim. Sci. Pap. Rep. 2014; 32: 95-105
        • Kim B.S.
        • Jung J.S.
        • Jang J.H.
        • Kang K.S.
        • Kang S.K.
        Nuclear Argonaute 2 regulates adipose tissue-derived stem cell survival through direct control of miR10b and selenoprotein N1 expression.
        Aging Cell. 2011; 10 (21241449): 277-291
        • Kinsella R.J.
        • Kahari A.
        • Haider S.
        • Zamora J.
        • Proctor G.
        • Spudich G.
        • Almeida-King J.
        • Staines D.
        • Derwent P.
        • Kerhornou A.
        • Kersey P.
        • Flicek P.
        Ensembl BioMarts: A hub for data retrieval across taxonomic space.
        Database (Oxford). 2011; 2011 (21785142)bar030
        • Kirkpatrick M.
        • Lofsvold D.
        • Bulmer M.
        Analysis of the inheritance, selection and evolution of growth trajectories.
        Genetics. 1990; 124 (2323560): 979-993
        • Ligthart S.
        • Marzi C.
        • Aslibekyan S.
        • Mendelson M.M.
        • Conneely K.N.
        • Tanaka T.
        • Colicino E.
        • Waite L.L.
        • Joehanes R.
        • Guan W.
        • Brody J.A.
        • Elks C.
        • Marioni R.
        • Jhun M.A.
        • Agha G.
        • Bressler J.
        • Ward-Caviness C.K.
        • Chen B.H.
        • Huan T.
        • Bakulski K.
        • Salfati E.L.
        • Fiorito G.
        • Wahl S.
        • Schramm K.
        • Sha J.
        • Hernandez D.G.
        • Just A.C.
        • Smith J.A.
        • Sotoodehnia N.
        • Pilling L.C.
        • Pankow J.S.
        • Tsao P.S.
        • Liu C.
        • Zhao W.
        • Guarrera S.
        • Michopoulos V.J.
        • Smith A.K.
        • Peters M.J.
        • Melzer D.
        • Vokonas P.
        • Fornage M.
        • Prokisch H.
        • Bis J.C.
        • Chu A.Y.
        • Herder C.
        • Grallert H.
        • Yao C.
        • Shah S.
        • McRae A.F.
        • Lin H.
        • Horvath S.
        • Fallin D.
        • Hofman A.
        • Wareham N.J.
        • Wiggins K.L.
        • Feinberg A.P.
        • Starr J.M.
        • Visscher P.M.
        • Murabito J.M.
        • Kardia S.L.R.
        • Absher D.M.
        • Binder E.B.
        • Singleton A.B.
        • Bandinelli S.
        • Peters A.
        • Waldenberger M.
        • Matullo G.
        • Schwartz J.D.
        • Demerath E.W.
        • Uitterlinden A.G.
        • van Meurs J.B.J.
        • Franco O.H.
        • Chen Y.-D. I.
        • Levy D.
        • Turner S.T.
        • Deary I.J.
        • Ressler K.J.
        • Dupuis J.
        • Ferrucci L.
        • Ong K.K.
        • Assimes T.L.
        • Boerwinkle E.
        • Koenig W.
        • Arnett D.K.
        • Baccarelli A.A.
        • Benjamin E.J.
        • Dehghan A.
        DNA methylation signatures of chronic low-grade inflammation are associated with complex diseases.
        Genome Biol. 2016; 17 (27955697): 255
        • Lordan R.
        • Zabetakis I.
        Invited review: The anti-inflammatory properties of dairy lipids.
        J. Dairy Sci. 2017; 100 (28342603): 4197-4212
        • Lourenco D.A.L.
        • Fragomeni B.O.
        • Bradford H.L.
        • Menezes I.R.
        • Ferraz J.B.S.
        • Aguilar I.
        • Tsuruta S.
        • Misztal I.
        Implications of SNP weighting on single-step genomic predictions for different reference population sizes.
        J. Anim. Breed. Genet. 2017; 134 (28833593): 463-471
        • MacGibbon A.K.H.
        • Taylor M.W.
        Composition and Structure of Bovine Milk Lipids. Advanced Dairy Chemistry. Volume 2 Lipids. 1–42.
        • Meuwissen T.H.
        • Hayes B.J.
        • Goddard M.E.
        Prediction of total genetic value using genome-wide dense marker maps.
        Genetics. 2001; 157 (11290733): 1819-1829
        • Misztal I.
        • Bradford H.L.
        • Lourenco D.A.L.
        • Tsuruta S.
        • Masuda Y.
        • Legarra A.
        Studies on inflation of GEBV in single-step GBLUP for type.
        Interbull Bull. 2017; XX: 38-42
        • Misztal I.
        • Tsuruta S.
        • Lourenco D.A.L.
        • Aguilar I.
        • Legarra A.
        • Vitezica Z.
        Manual for BLUPF90 family of programs.
        University of Georgia, Athens2014
        • Misztal I.
        • Tsuruta S.
        • Strabel T.
        • Auvray B.
        • Druet T.
        • Lee D.
        BLUPF90 and related programs (BGF90).
        in: Proc. 7th World Congr. Genet. Appl. Livest. Prod. Montpellier, France. Editions Quae, Montpellier, France2002: 21-22
        • Narayana S.G.
        • Schenkel F.S.
        • Fleming A.
        • Koeck A.
        • Malchiodi F.
        • Jamrozik J.
        • Johnston J.
        • Sargolzaei M.
        • Miglior F.
        Genetic analysis of groups of mid-infrared predicted fatty acids in milk.
        J. Dairy Sci. 2017; 100 (28342614): 4731-4744
        • Ntambi J.M.
        The regulation of stearoyl-CoA desaturase (SCD).
        Prog. Lipid Res. 1995; 34 (7480063): 139-150
        • Oliveira H.R.
        • Brito L.F.
        • Lourenco D.A.L.
        • Silva F.F.
        • Jamrozik J.
        • Schaeffer L.R.
        • Schenkel F.S.
        Invited review: Advances and applications of random regression models: From quantitative genetics to genomics.
        J. Dairy Sci. 2019; 102 (31255270): 7664-7683
        • Oliveira H.R.
        • Lourenco D.A.L.
        • Masuda Y.
        • Misztal I.
        • Tsuruta S.
        • Jamrozik J.
        • Brito L.F.
        • Silva F.F.
        • Cant J.P.
        • Schenkel F.S.
        Single-step genome-wide association for longitudinal traits of Canadian Ayrshire, Holstein, and Jersey dairy cattle.
        J. Dairy Sci. 2019; 102 (31477296): 9995-10011
        • Oliveira H.R.
        • Lourenco D.A.L.
        • Masuda Y.
        • Misztal I.
        • Tsuruta S.
        • Jamrozik J.
        • Brito L.F.
        • Silva F.F.
        • Schenkel F.S.
        Application of single-step genomic evaluation using multiple-trait random regression test-day models in dairy cattle.
        J. Dairy Sci. 2019; 102 (30638992): 2365-2377
        • Palmquist D.L.
        • Beaulieu A.D.
        • Barbano D.M.
        Feed and animal factors influencing milk fat composition.
        J. Dairy Sci. 1993; 76 (8326036): 1753-1771
        • Palombo V.
        • Milanesi M.
        • Sgorlon S.
        • Capomaccio S.
        • Mele M.
        • Nicolazzi E.
        • Ajmone-Marsan P.
        • Pilla F.
        • Stefanon B.
        • D'Andrea M.
        Genome-wide association study of milk fatty acid composition in Italian Simmental and Italian Holstein cows using single nucleotide polymorphism arrays.
        J. Dairy Sci. 2018; 101 (30243637): 11004-11019
        • Parodi P.
        Milk fat in human nutrition.
        Aust. J. Dairy Technol. 2004; 59: 3-59
        • Pepe G.
        • Tenore G.C.
        • Mastrocinque R.
        • Stusio P.
        • Campiglia P.
        Potential anticarcinogenic peptides from bovine milk.
        J. Amino Acids. 2013; 2013 (23533710)939804
        • R Core Team
        R: A language and environment for statistical computing.
        R Foundation for Statistical Computing, Vienna, Austria2018
        • Raftery A.L.
        • Lewis S.
        Comment: One long run with diagnostics: Implementation strategies for Markov chain Monte Carlo.
        Stat. Sci. 1992; 7: 493-497
        • Rutten M.J.M.
        • Bovenhuis H.
        • Hettinga K.A.
        • van Valenberg H.J.F.
        • van Arendonk J.A.M.
        Predicting bovine milk fat composition using infrared spectroscopy based on milk samples collected in winter and summer.
        J. Dairy Sci. 2009; 92 (19923625): 6202-6209
        • Samková E.
        • Špička J.
        • Pešek M.
        • Pelikánová T.
        • Hanuš O.
        Review: Animal factors affecting fatty acid composition of cow milk fat: A review.
        S. Afr. J. Anim. Sci. 2012; 42: 83-100
        • Sargolzaei M.
        • Chesnais J.P.
        • Schenkel F.S.
        A new approach for efficient genotype imputation using information from relatives.
        BMC Genomics. 2014; 15: 478
        • Sauer F.D.
        • Fellner V.
        • Kinsman R.
        • Kramer J.K.G.
        • Jackson H.A.
        • Lee A.J.
        • Chen S.
        Methane output and lactation response in Holstein cattle with monensin or unsaturated fat added to the diet.
        J. Anim. Sci. 1998; 76: 906-914
        • Schaeffer L.R.
        • Jamrozik J.
        • Kistemaker G.J.
        • Van Doormaal J.
        Experience with a test-day model.
        J. Dairy Sci. 2000; 83 (10821590): 1135-1144
        • Smith B.J.
        boa: An R Package for MCMC output convergence assessment and posterior inference.
        J. Stat. Softw. 2007; 21: 1-37
        • Soyeurt H.
        • Dardenne P.
        • Dehareng F.
        • Lognay G.
        • Veselko D.
        • Marlier M.
        • Bertozzi C.
        • Mayeres P.
        • Gengler N.
        Estimating fatty acid content in cow milk using mid-infrared spectrometry.
        (J. Dairy Sci. 89: 3690-3695.) (16899705)
        • Stryer L.
        Fatty acid metabolism.
        in: Biochemistry. 4th ed. W. H. Freeman and Company, New York, NY1995: 603-628
        • Temme E.H.
        • Mensink R.P.
        • Hornstra G.
        Comparison of the effects of diets enriched in lauric, palmitic, or oleic acids on serum lipids and lipoproteins in healthy women and men.
        Am. J. Clin. Nutr. 1996; 63 (8644684): 897-903
        • VanRaden P.M.
        Efficient methods to compute genomic predictions.
        J. Dairy Sci. 2008; 91 (18946147): 4414-4423
        • Viitala S.M.
        • Schulman N.F.
        • de Koning D.J.
        • Elo K.
        • Kinos R.
        • Virta A.
        • Virta J.
        • Mäki-Tanila A.
        • Vilkki J.H.
        Quantitative trait loci affecting milk production traits in Finnish Ayrshire dairy cattle.
        J. Dairy Sci. 2003; 86 (12778594): 1828-1836
        • Wang H.
        • Jiang L.
        • Liu X.
        • Yang J.
        • Wei J.
        • Xu J.
        • Zhang Q.
        • Liu J.-F.
        A post-GWAS replication study confirming the PTK2 gene associated with milk production traits in Chinese Holstein.
        PLoS One. 2013; 8 (24386238)e83625
        • Wang H.
        • Misztal I.
        • Aguilar I.
        • Legarra A.
        • Muir W.M.
        Genome-wide association mapping including phenotypes from relatives without genotypes in a single-step (ssGWAS) for 6-week body weight in broiler chickens.
        Front. Genet. 2014; 20: 134
        • Wang H.
        • Misztal I.
        • Aguilar I.
        • Legarra A.
        • Fernando R.L.
        • Vitezica Z.
        • Okimoto R.
        • Wing T.
        • Hawken R.
        • Muir W.M.
        Genome-wide association mapping including phenotypes from relatives without genotypes.
        Genet. Res. (Camb.). 2012; 94 (22624567): 73-83
        • Wang X.
        • Ma P.
        • Liu J.
        • Zhang Q.
        • Zhang Y.
        • Ding X.
        • Jiang L.
        • Wang Y.
        • Zhang Y.
        • Sun D.
        • Zhang S.
        • Su G.
        • Yu Y.
        Genome-wide association study in Chinese Holstein cows reveal two candidate genes for somatic cell score as an indicator for mastitis susceptibility.
        BMC Genet. 2015; 16 (26370837): 111
        • Williams C.M.
        Dietary fatty acids and human health.
        Ann. Zootech. 2000; 49: 165-180
        • Zhao A.
        • Jiang F.
        • Yang G.
        • Liu H.
        • Li L.
        Sfrp5 interacts with Slurp1 to regulate the accumulation of triglycerides in hepatocyte steatosis model.
        Biochem. Biophys. Res. Commun. 2019; 512 (30879770): 256-262
        • Zong G.
        • Li Y.
        • Wanders A.J.
        • Alssema M.
        • Zock P.L.
        • Willett W.C.
        • Hu F.B.
        • Sun Q.
        Intake of individual saturated fatty acids and risk of coronary heart disease in US men and women: Two prospective longitudinal cohort studies.
        BMJ. 2016; 355 (27881409)i5796