Implications
Excessive abdominal fat in Holstein cows with similar body condition scores (BCS) is not caused by a single Mendelian locus, and the trait appears to be at least moderately heritable. Selection to reduce the prevalence of excessive abdominal fat is potentially possible, which may allow the prevention of metabolic diseases related to fat deposition dynamic in dairy cattle.
Introduction
The pregnant dairy cow undergoes a series of complex metabolic and physiological changes as parturition approaches. As a result, periparturient disorders may occur affecting the future performance of dairy cows (Goff and Horst, Reference Goff and Horst1997). Dry matter intake begins to decrease a few weeks before parturition with the lowest level occurring at calving, whereas simultaneously nutrient demands increase due to the exponential growth of the fetus and the initiation of lactation (Ingvartsen and Andersen, Reference Ingvartsen and Andersen2000). Consequently, the cow typically experiences a state of negative energy balance and a characteristic mobilization of fat from adipose tissue stored in different areas of the body. Extreme lipid mobilization leads to an increased release from adipose tissue and uptake of non-esterified fatty acids by the liver with a pathologic accumulation of hepatic triglycerides (Drackley, Reference Drackley1999; De Koster and Opsomer, Reference De Koster and Opsomer2013). In addition to an energy source, adipose tissue also has immune, endocrine, regenerative, mechanical and thermal functions. The fuel and non-fuel functions of adipose tissue depots vary with depot size and body-fat distribution (Tchkonia et al., Reference Tchkonia, Thomou, Zhu, Karagiannides, Pothoulakis, Jensen and Kirkland2013).
Fat is deposited under the epidermis and around vital organs, where it may play immunologically defensive and mechanically protective roles (Duffaut et al., Reference Duffaut, Zakaroff-Girard, Bourlier, Decaunes, Maumus, Chiotasso, Sengenès, Lafontan, Galitzky and Bouloumié2009). Once inflamed, adipose tissue shifts from storing to releasing fatty acids, potentially driven in part through local proinflammatory cytokine release. Regional differences in pre-adipocytes replication, differentiation, abundance, susceptibility to apoptosis and gene expression may contribute to regional variation in fat-tissue function (Tchkonia et al., Reference Tchkonia, Thomou, Zhu, Karagiannides, Pothoulakis, Jensen and Kirkland2013). Variation in gene expression between different fat depots has been demonstrated in dairy cattle (Saremi et al., Reference Saremi, Sauerwein, Dänicke and Mielenz2012).
Based upon on these findings in humans and lab animals, it is noteworthy that genome-wide association studies have suggested that body-fat distribution is associated with variation in genes involved in pattern formation during embryonic development, angiogenesis, pre-adipocyte signaling and adipocyte development (Heid et al., Reference Heid, Jackson, Randall, Winkler, Qi, Steinthorsdottir, Thorleifsson, Zillikens, Speliotes and Mägi2010). These associations with genetic variants in orthologous genes may also occur in dairy cattle. Indeed, a New Zealand study found that Holstein strains (New Zealand v. US cows) had different fatty acid profiles in their adipose tissues and milk, which may differentially affect the metabolic status of the adipose depots (Meier et al., Reference Meier, Verkerk, Kay, Macdonald and Roche2013). Other studies have demonstrated that desaturase and lipase activities differ between abdominal and subcutaneous fat, supporting the hypothesis of a preferential mobilization of abdominal fat in dairy cows, especially when they develop left displacement of the abomasum (Hostens et al., Reference Hostens, Fievez, Leroy, Van Ranst, Vlaeminck and Opsomer2012). This may imply that some of the genes responsible for variation in abdominal fat accumulation are also related to variation in the risk of development of this digestive condition in Holstein cattle. Indeed, some bovine practitioners have observed that it is becoming very common to find cows with a displacement of the abomasum and excessive omental fat, but with a normal BCS during the surgical correction of this digestive disorder (Melendez P and Poock S, unpublished data). Therefore, identifying risk factors, such as genetic variants predisposing cows to this type of periparturient diseases, would be beneficial to the dairy industry. This might enable the identification of genetic markers that are predictive of variation in a visceral fat deposition, which in turn could permit the establishment of proper management and selection strategies to prevent and control fat-tissue-related disorders. Consequently, the objective of this study was to assess the extent to which genetic variation was responsible for differences in the degree of visceral adiposity in US Holstein cows. We further sought to determine if there were variants with large phenotypic effects responsible for abnormally high visceral fat deposition in cows with a normal subcutaneous fat deposition.
Material and methods
Facilities and animal care and use committee
This study was approved by the Animal Care and Use Committee of the University of Missouri, protocol number 8443.
The study population consisted of Holstein cows processed in a beef harvesting facility located in Green Bay, WI, USA. This plant was selected because it processes more than 2000 dairy cows per day and provided access to a large number of animals with similar characteristics for which their carcass omental, visceral and subcutaneous fat could be straightforwardly appraised. In addition, all the cows came from dairies of the state of Wisconsin that have similar nutritional management and were confined in free-stall systems. A preliminary visit to the abattoir was conducted by one of the investigators to establish a phenotyping protocol to identify cows with different propensities for fat deposition (visceral v. subcutaneous). The sampling sessions were subsequently performed over two visits (2 days each) during September 2016.
Experimental design
Adult Holstein cows were selected from the processing line based upon the presence of permanent unworn incisor teeth before the neck emerging (4 to 6 years old). Therefore, there was minor variation in age. After the hide had been removed and the carcass eviscerated, the animals were evaluated for BCS (Ferguson et al., Reference Ferguson, Galligan and Thomsen1994) and an assessment of the amount of omental fat at the level of the insertion of the lesser omentum over the pylorus area was conducted. A visual score after measuring a fold of the lesser omentum with a caliper was established. Score 0 was defined as an omental fold <5 mm in thickness and minimum fat deposition observed throughout the entire omentum. Score 1 was defined as an omental fold ⩾20 mm in thickness and marked fat deposition observed throughout the entire omentum. Animals with an omentum thickness between 5 and 20 mm were not eligible for selection and were excluded from the study. The omentum was evaluated using a visual score, as the thickness of the omentum could not be measured due to the required viewing location and the accelerated pace of the slaughter line. The evaluator was not permitted to stand in the area of evisceration.
Further, only animals with normal BCS between 2.75 and 3.25 were included.
A group of 100 Holstein cows with an omental fat score of 0 and a group of 100 cows with an omental fat score of 1 were selected. The purpose of this sampling protocol was to obtain two groups of animals with extreme differences in omental fat deposition but for which animals within both groups had a moderate subcutaneous fat deposition (BCS 2.75 to 3.25). Although Holstein steers are known to have increased amounts of omental fat relative to beef breed steers at the same live weight (Schaefer, Reference Schaefer2005), no studies appear to have reported the heritability or genetic variance associated with variation in the amount of omental fat within Holstein cattle. Consequently, it was not possible for us to design this experiment based on power calculations concerning the probabilities of detection of quantitative trait loci (QTL) of certain magnitudes of effect. Although the study is clearly underpowered for the detection of QTL underlying variation in omental fat, we hypothesized that the phenotype of excessive omental fat might represent a metabolic disorder that could be caused by the segregation of a single locus of large effect. It is well known that sample sizes required to association map Mendelian loci can be much smaller than those required to map QTL, and n=200 is more than sufficient for this purpose when there is a large allele frequency between the case and control groups and a large effect due to the causal locus (Altshuler et al., Reference Altshuler, Daly and Lander2008).
After the animals had been selected, a small piece of muscle from the neck area was collected into a sterile container for DNA extraction.
Genetic analysis
Muscle samples were submitted to a commercial laboratory (GeneSeek a Neogen Company, Lincoln, NE, USA) for interrogation of genome-wide genomic variation through the use of the Illumina BovineHD Beadchip containing over 777 K single nucleotide polymorphism (SNP) probes that span the bovine genome. The visceral/subcutaneous fat deposition phenotype was treated as a binary trait based on the classification of cases (Score 1) and controls (Score 0) in the study population. Quality control was performed on the genotypic data and genotypes for SNPs and samples were retained for analysis if the following criteria were met: (1) the minor allele frequency (MAF) >5%; (2) the percentage of missing genotypes across all samples <10%; (3) there was no evidence for a strong deviation from Hardy–Weinberg equilibrium (P<1E-6) (Wigginton et al., Reference Wigginton, Cutler and Abecasis2005); and (4) sample SNP call rate was >90%. The quality control of the data was deemed to be adequate as using an individual genotyping call rate ≤85% allows for high-quality genotype calls and minimal genotyping errors regardless of MAF (Purfield et al., Reference Purfield, McClure and Berry2016). The sex check option available in PLINK 1.9 (Purcell et al., Reference Purcell, Neale, Todd-Brown, Thomas, Ferreira, Bender, Maller, Sklar, de Bakker, Daly and Sham2007) was used to impute sex with a maximum F statistic of 0.5 for females and a minimum value of 0.75 for males, resulting in five individuals being identified as males. After filtering, 589 732 SNPs and 189 cattle were retained.
Genome-wide Efficient Mixed Model Association (GEMMA) software was used to calculate a genomic relationship matrix, and the centered relatedness matrix option was used as suggested by the authors (Zhou and Stephens, Reference Zhou and Stephens2012). Genome-wide association analyses were conducted to test potential associations between fat deposition and genomic variation. A univariate mixed linear model analysis was performed in GEMMA to identify SNPs significantly associated with variation in visceral fat deposition. Models were fit with sex as a covariate.
As denoted by Zhou and Stephens, the software fits the following univariate mixed model:
$$\eqalign {{\rm y}{\equals}&{\bf W\balpha }{\plus}{\bf x}\beta {\plus}{\bf u}{\plus}{\bf {\bvarepsilon}} \hbox{;} \,\,{\bf u}\,\sim\,{\rm MVN}_{n} \left( {0,\,\lambda \tau ^{{{\minus}1}} {\bf K}} \right),\, \bvarepsilon \, \cr & \hskip-3pt \sim\, {\rm MVN}_{n} \left( {0,\,\tau ^{{{\minus}1}} {\bf I}_{n} } \right),$$
where y is an n-vector of quantitative traits (or binary trait labels) for the n individuals; W=(w1, …, wc) is an n × c matrix of covariates including a column of 1 s; α is a c-vector of the corresponding coefficients including the intercept; x is an n-vector of marker genotypes; β is the effect size of the marker; u is an n-vector of random effects; ε is an n-vector of errors; τ −1 is the variance of the residual errors; λ is the ratio between the two variance components; K is a known n × n relatedness matrix; and I n is an n × n identity matrix. MVN n denotes the n-dimensional multivariate normal distribution (Zhou and Stephens, Reference Zhou and Stephens2012).
GEMMA tests the alternative hypothesis H1: β≠0 against the null hypothesis; H0: β=0 for each SNP in turn, using one of the three commonly used test statistics (Wald, likelihood ratio or score). GEMMA obtains either the maximum likelihood estimate or the restricted maximum likelihood estimate of λ and β and outputs the corresponding P-value. In addition, GEMMA estimates the PVE by typed genotypes or ‘chip heritability’ (Zhou and Stephens, Reference Zhou and Stephens2012).
In addition, pathway analyses were used to identify any metabolic networks that might be enriched for visceral fat deposition-associated candidate genes. To facilitate this, a list containing all genes located within a 250 kb window on either side of the significantly associated variants were supplied to DAVID Bioinformatics Resources 6.8 (Huang et al., Reference Huang, Sherman and Lempicki2009a and Reference Huang, Sherman and Lempicki2009b).
Results
Q–Q plots produced for the models revealed that, as expected, the study was underpowered but possessed sufficient power to detect the largest effect QTLs (Figure 1). Figure 1 also shows that including the effect of sex in the model slightly improved the power of the association analysis. The chip heritability was 0.686 for the model that excluded the effect of animal sex and 0.692 for the model that included sex, and for this model, the estimated additive genetic variance was 0.427, and the residual (environmental and non-additive genetic) variance estimate was 0.074.
(Color online) Q–Q plots for models (a) including and (b) excluding the effect of animal sex in dairy cattle. The observed −log10Wald P-value is shown on the y-axis and the expected −log10Wald P-value is shown on the x-axis.
In total, 11 SNPs defining four QTL regions were found to be significantly associated with variation in a visceral fat deposition at a Wald’s P<10−5 (Table 1). Two of the QTL were detected with four and five significantly associated SNPs, respectively; whereas, the QTLs detected on BTA12 and BTA19 (at 38 Mb) were each detected with only one significantly associated SNP. The associated SNPs in the QTLs with more than one marker were found to be in linkage disequilibrium (r 2=1 and D'=1) (Purcell et al., Reference Purcell, Neale, Todd-Brown, Thomas, Ferreira, Bender, Maller, Sklar, de Bakker, Daly and Sham2007). Figure 2 shows the genome-wide association studies Manhattan plots for the models including and excluding the effect of animal sex. Genes annotated in the UMD3.1 reference assembly within the four QTL regions are reported in Table 2 and their biological functions are reported in Table 3. Other than the cluster of HOXB genes on BTA19, no enriched gene ontology terms were found within the gene networks harboring these genes when supplied to DAVID using either the Bos taurus or human gene ontology databases.
(Color online) Manhattan plots representing two genome-wide association analyses for visceral adiposity (a) including an effect for animal sex or (b) excluding the effect of sex in dairy cattle. Both plots show −log10Wald P-values on the y-axis, and the genome-wide significance threshold of P<10−5 is shown in blue.
Single nucleotide polymorphisms (SNPs) significantly associated with the degree of visceral adiposity in US Holstein cows
Mb=mega base; rs=reference SNP identification number; BTA=Bos taurus autosomal chromosome.
List of bovine genes supplied to DAVID Bioinformatics Resources 6.8
Mb=mega base; SNPs=single nucleotide polymorphisms; FGF14=fibroblast growth factor 14; BTA=Bos taurus autosomal chromosome; rs=reference SNP identification number; GDPD1=glycerophosphodiester phosphodiesterase domain containing 1; YPEL2=Yippee-like 2; DHX40=DEAH (Asp-Glu-Ala-His) box polypeptide 40; CLTC=clathrin, heavy chain (Hc); PTRH2=peptidyl-tRNA hydrolase 2; VMP1=vacuole membrane protein 1; HOXB13=homeobox B13; PRAC2=HOXB cluster antisense RNA 5; HOXB9=homeobox B9; HOXB8=homeobox B8; HOXB7=homeobox B7; HOXB6=homeobox B6; HOXB5=homeobox B5; HOXB4=homeobox B4; HOXB3=homeobox B3; HOXB2=homeobox B2; HOXB1=homeobox B1; SKAP1=Src kinase associated phosphoprotein 1.
Functions of annotated bovine genes within the four identified quantitative trait locus (QTL) regions
CLTC=clathrin, heavy chain (Hc); DHX40=DEAH (Asp-Glu-Ala-His) box polypeptide 40; FGF14=fibroblast growth factor 14; GDPD1=glycerophosphodiester phosphodiesterase domain containing 1; HOXB=homeobox B Family; PRAC2=HOXB cluster antisense RNA 5; PTRH2=peptidyl-tRNA hydrolase 2; SKAP1=Src kinase associated phosphoprotein 1; VMP1=vacuole membrane protein 1; YPEL2=Yippee-like 2.
Discussion
Our study used 777 k genotypes to test potential associations between genetic variation and the degree of visceral adiposity in Holstein cows. The chip heritability represents the proportion of phenotypic variation explained by the markers included on the genotyping assay and is thus dependent on the number of markers and their allele frequency spectrum, which generally represents common SNPs within taurine cattle. The chip heritability and variance components obtained in this study are biased estimates for omental fat deposition in Holstein cattle because the sampling scheme required equal numbers of animals to represent the tails of the phenotypic distribution defined as low and excessive visceral fat deposition classes. However, our sampling scheme was deliberately chosen in order to maximize the likelihood of detecting trait-associated SNPs using a relatively small sample size. In this regard, we were clearly successful and produced large estimates for the chip heritability and additive genetic variation, which indicates that the trait is heritable. However, we did not find evidence for a single Mendelian locus of large effect underlying the study phenotype, which was based on our hypothesis that the accumulation of excessive omental fat represented a Mendelian defect in Holstein cattle. We found evidence for quantitative inheritance of the phenotype, but rather surprisingly given the lack of power in our design due to small sample size; we found 11 significantly associated SNPs at P<10−5 which predicted four QTLs underlying this phenotype (Figure 2). These loci may well include type 1 errors as two of the QTLs were each identified by only one significantly associated SNP and the network analysis for the genes located in the QTL regions failed to reveal any enriched gene ontology terms within the networks in which these genes participate. Notably, one of the QTL regions identified by a single associated SNP contained a cluster of HOXB genes and HOXC family members have previously been shown to be associated with variation in fat distribution in humans (Heid et al., Reference Heid, Jackson, Randall, Winkler, Qi, Steinthorsdottir, Thorleifsson, Zillikens, Speliotes and Mägi2010). The HOXB gene family, in particular, has been found to be associated with the regulation of tissue differentiation, cell cycle control and embryonic patterning (Bene and Wittbrodt, Reference Bene and Wittbrodt2005).
Based upon on our results, this study supports the hypothesis that body-fat distribution is associated with variation in genes, perhaps most likely involved in pattern formation during embryonic development, angiogenesis, pre-adipocyte signaling and adipocyte development, all pointing toward heritable differences in adipocyte characteristics, as occurs in humans (Heid et al., Reference Heid, Jackson, Randall, Winkler, Qi, Steinthorsdottir, Thorleifsson, Zillikens, Speliotes and Mägi2010). Indeed, it has been shown that there are diverse genetic Holstein lines that had different fatty acid profiles in their adipose tissues, as well as in milk, which may differentially affect the metabolic status of the adipose depots (Meier et al., Reference Meier, Verkerk, Kay, Macdonald and Roche2013). In this sense, depending upon the type and distribution of adipose tissue (visceral v. subcutaneous), it has been demonstrated that the metabolism of these two types of fat is different. Visceral fat is much more metabolically active than subcutaneous fat. Indeed, the content of adiponectin varies within these two types of fat depots. This adipokine has an anti-lipolytic effect and stimulates insulin sensitivity. The synthesis of adiponectin is lower in retroperitoneal fat than subcutaneous fat, suggesting a higher rate of lipolysis for retroperitoneal fat (Singh et al., Reference Singh, Häussler, Heinz, Akter, Saremi, Müller and Rehage2014). In addition, the hormone-sensitive lipase protein expression and extent of its phosphorylation were significantly higher in retroperitoneal fat compared with subcutaneous fat at 21 days postpartum in Holstein cows (Locher et al., Reference Locher, Meyer, Weber, Rehage, Meyer, Dänicke and Huber2011). Other studies have also demonstrated that desaturase and lipase activities are higher in abdominal than subcutaneous fat in dairy cows, especially when they develop left displacement of the abomasum (LDA) (Hostens et al., Reference Hostens, Fievez, Leroy, Van Ranst, Vlaeminck and Opsomer2012). As a result of the aforementioned studies in dairy cattle, the following hypothesis is established: because abdominal fat is more active, an excessive accumulation of omental deposits due to genetic factors may predispose cows to a larger release of non-esterified fatty acids from abdominal triglycerides during the postpartum period. In addition, lower adipokine expression may decrease insulin sensitivity. The cow will respond with a higher uptake of fatty acids by the liver. Hepatic non-esterified fatty acids upregulate expression of gluconeogenic, beta-oxidative, lipogenic and ketogenic genes, promoting hyperglycemia, hyperlipidemia and ketosis. In turn, β-OH butyrate has been shown to decrease peripheral insulin-stimulated glucose uptake, and may contribute to the development of insulin resistance in individuals with excessive fat accumulation in the liver (Geisler and Renquist, Reference Geisler and Renquist2017). The excess of ketone bodies may also predispose the cow to a lower dry matter intake, creating a vicious cycle of excessive fat mobilization and ketogenesis. In addition, high levels of non-esterified fatty acids are also related to hypocalcemia. There will be rumen and abomasum atony due to low blood calcium concentrations. Considering a lower rumen content and a prominent space in front of the rumen after the calf is delivered, the flaccid abomasum will fill with gas and will lay down under the rumen and migrate from the right to the left side of the abdominal cavity, causing the typical left displacement of the abomasum (Melendez and Risco, Reference Melendez and Risco2016). This digestive condition may be a consequence of excessive abdominal fat accumulation or may imply that some of the genes responsible for variation in abdominal fat accumulation, which is a more metabolic active tissue, are also related to variation in the risk of development of LDA in Holstein dairy cattle (Mömke et al., Reference Mömke, Sickinger, Lichtner, Doll, Rehage and Distl2013). Consequently, diseases such as ketosis and LDA that are intimately related to fat metabolism might be approached through genetic management. These periparturient disorders result in significant economic losses to the dairy industry through reductions in reproductive performance and milk yield during the subsequent lactation, cost of treatments and increased risk of culling. As an example, a cost of an LDA in a multiparous cow can reach a value of US$ 630, a value between a half and a third of the price of the animal (Liang et al., Reference Liang, Arnold, Stowe, Harmon and Bewley2017). Therefore, the genetic differences between subcutaneous and visceral fat deposition, demonstrated in the present study, might suggest that the identification of potentially genetic markers in the future, included in new traits, predicting a higher degree of visceral fat deposition, could enable a genomic selection strategy to prevent and control adipose tissue related disorders. However, the most important challenge for future genomic selection will be the incorporation of these new traits such as health characteristics, and feed efficiency into selection indexes (Weller et al., Reference Weller, Ezra and Ron2017).
One of the shortcomings of this study was that the sampling of animals was conducted without knowing data about their physiological state, such as being lactating or dry, healthy or sick, etc. In addition, during the sampling process, five males were mistakenly selected due to the dynamic and rapid movement through the slaughter line. However, the use of this slaughterhouse permitted the evaluation of a significant number of animals for their abdominal and subcutaneous fat content in a fast way. Despite the sampling error in relation to gender (five males) and the lack of homogeneity in the physiological characteristics of the animals, it was possible to demonstrate genomic differences in relation to the distribution of abdominal fat in Holstein cattle. The males were retained in the evaluated population and effect of sex was added into the model to prevent further reduction of the sample size. This was a pilot study that will allow to establish new research hypothesis in more controlled studies, especially related to the gene expression of key enzymes related to the metabolism of the visceral fat in the Holstein breed.
It is concluded that excessive omental fat in Holstein cows with similar BCS is not caused by a single Mendelian locus and the trait appears to be at least moderately heritable. Selection to reduce the prevalence of excessive omental fat is potentially possible but would require the generation of predicted transmitting abilities from extensively collected phenotypes or the development of molecular estimates of predicted transmitting abilities produced from a much larger and random sample of Holsteins.
Acknowledgments
The authors thank the slaughterhouse for allowing us to conduct this study. This work was supported by the USDA National Institute of Food and Agriculture (AHFF), project 1006957.
Declaration of interest
The authors declare no potential conflict of interest for this study.
Ethics statement
This study was approved by the Animal Care and Use Committee of the University of Missouri, protocol number 8443.
Software and data repository resources
Genome-wide Efficient Mixed Model Association (GEMMA) software was used to calculate a genomic relationship matrix (Zhou and Stephens, Reference Zhou and Stephens2012). DAVID Bioinformatics Resources 6.8 were also used (Huang et al., Reference Huang, Sherman and Lempicki2009a and Reference Huang, Sherman and Lempicki2009b). Data are available from the corresponding author for non-commercial uses.
