Lactation in mammals is an energy-demanding process and requires an increase in feed intake. When feed is utilised to produce milk in cows, heat is generated as a by-product of the metabolic processes, thus an increase in milk production will result in increased heat generation. High milk yielding cows are therefore more likely to suffer from heat stress than low yielding cows in hotter climates, and heat stress is a factor that negatively affects milk production, especially for cows of high genetic merit (Thatcher & Collier, Reference Thatcher, Collier and Morrow1986; Kadzere et al. Reference Kadzere, Murphy, Silanikove and Maltz2002; Bryant et al. Reference Bryant, Lόpez-Villalobos, Pryce and Holmes2007).
One measure of heat load for animals is the Temperature Humidity Index (THI) which provides a simple metric that combines temperature and humidity measures (Dikmen & Hansen, Reference Dikmen and Hansen2009), with the suggestion that THI values above 72 consistently cause heat stress in dairy cows (Ravagnolo & Misztal, Reference Ravagnolo and Misztal2000). This figure is now widely cited (Kendall et al. Reference Kendall, Verkerk, Webster and Tucker2007). Use of such a threshold value suggests there is a ‘tipping-point’ for thermal index after which a decline of production occurs, and that all cows have the same threshold and rate of decline in every environment. However, these assumptions may be wrong, with Sánchez et al. (Reference Sánchez, Misztal, Aguilar, Zumbach and Rekaya2009) proposing that the threshold and production decline rate are both different and variable for individual cows. Brügemann et al. (Reference Brügemann, Gernand, von Borstel and König2011) and Hammami et al. (Reference Hammami, Bormann, M'hamdi, Montaldo and Gengler2013) reported lower thresholds for a temperature-humidity index (THI) for production traits in temperate environments, compared with those reported under tropical and subtropical conditions (Bohmanova et al. Reference Bohmanova, Misztal, Tsuruta, Norman and Lawlor2008; Aguilar et al. Reference Aguilar, Misztal and Tsuruta2009; Boonkum et al. Reference Boonkum, Misztal, Duangjinda, Pattarajinda, Tumwasorn and Sanpote2011), which suggests that we may need to have greater consideration of genotype by environment interactions in understanding thermal stress and its effect on milk production.
It has been suggested that the performance of lactating animals is limited by their capacity to dissipate excessive body heat and thus maintain homeothermy (Król & Speakman, Reference Król and Speakman2003; Speakman & Król, Reference Speakman and Król2005). Accordingly, increases in milk production would require decreases in heat generation from other metabolic processes, such as thermogenesis in brown adipose tissue (BAT). In rodents, the thermogenic activity of BAT is suppressed during lactation (Trayhurn et al. Reference Trayhurn, Douglas and McGuckin1982; Nizielski et al. Reference Nizielski, Billington and Levine1993) and Król et al. (Reference Król, Martin, Huhtaniemi, Douglas and Speakman2011) reported that thermogenic genes in mice are down-regulated during lactation, with the greatest level of down-regulation (a two-thirds decrease) being observed for the uncoupling protein-1 (UCP1) gene (UCP1).
UCP1 is a member of the mitochondrial carrier-protein family. It is predominantly expressed in BAT and plays important roles in maintaining body temperature, regulating metabolic rate and controlling energy expenditure via both non-shivering thermogenesis (NST) and diet-induced thermogenesis (Cannon & Nedergaard, Reference Cannon and Nedergaard2004; Kozak et al. Reference Kozak, Koza and Anunciado-Koza2010). It does this by uncoupling the respiratory chain, thus allowing for faster substrate oxidation with reduced ATP synthesis.
UCP1-ablated mouse studies have produced contrasting results. In one study the UCP1-ablated mice had increased sensitivity to cold exposure, but failed to become obese on either a standard or a high-fat diet, when housed at room temperature (Enerback et al. Reference Enerback, Jacobsson, Simpson, Guerra, Yamashita, Harper and Kozak1997). However, in a contrasting study (Feldmann et al. Reference Feldmann, Golozoubova, Cannon and Nedergaard2009), the authors suggested that UCP1-ablation does induce obesity even in mice fed on a normal diet, and vastly increased obesity on high-fat diets, when the mice were housed at thermoneutrality. These authors reconciled their findings with the earlier model by suggesting that under the conditions that the mice were housed in the first study they were effectively too cold, and thus to preserve their body temperature they had to increase their metabolism and food intake. At thermoneutrality (~30 °C for mice) UCP1-ablation therefore becomes obesogenic. Given its role in thermogenesis and fat deposition, and the change in its expression level during lactation, UCP1 may be hypothesised to have a role in milk production and variation in UCP1 may, therefore, affect milk traits in cows.
UCP1 structure is conserved across mammals and the gene contains six small coding exons separated by five introns. It spans over 6 kb in total. Variation in UCP1 has been extensively investigated in humans and a number of single nucleotide polymorphisms (SNPs) in the promoter region and exons have been associated with obesity risk and/or type 2 diabetes mellitus risk (Jia et al. Reference Jia, Tian, Cao, Tao, Zhang, Gao, Ge, Lin and Jois2010; Brondani et al. Reference Brondani, Assmann, Duarte, Gross, Canani and Crispim2012). In sheep, variation in UCP1 has been associated with subcutaneous carcass fat depth and lean meat yield (Yang et al. Reference Yang, Forrest, Zhou, Hodge and Hickford2014).
A cDNA sequence for bovine UCP1 has been known for over two decades (Casteilla et al. Reference Casteilla, Bouillaud, Forest and Ricquier1989), but there is only one study reporting variation in the bovine UCP1 intron 1 (Sonstegard & Kappes, Reference Sonstegard and Kappes1999). No investigation has been undertaken into variation in other gene regions, or the effect of this genetic variation, should it exist, on production traits.
In this study, we use PCR-single strand conformational polymorphism (PCR-SSCP) analysis to investigate variation in four regions of bovine UCP1, and then use two of these variable regions to construct haplotypes spanning the promoter, 5′-UTR and exon 2. The impact of the common UCP1 haplotypes on milk production traits was then investigated in a herd of 612 Holstein-Friesian × Jersey (HF × J) cross dairy cows in a pasture-based production system.
Materials and methods
Cows and DNA samples
Six hundred and twelve HF × J cross (KiwiCross™) dairy cows that were three to nine years old (i.e. 2nd to 8th parity) and had commenced lactation within ten weeks of the study were investigated. Of these, 415 cows had measurements from one lactation and 197 cows had measurements from two lactations. This breed cross is not a simple 50 : 50 cross of parental genetics, and may involve 50 : 50 cross bulls being mated to cows of either pure parental breed, or to another HF × J cross cow. All the cows were maintained on the Lincoln University Dairy Farm (LUDF; Canterbury, New Zealand) and grazed over their lactation on perennial ryegrass-white clover pasture. A blood sample from each cow was collected onto an FTA card and allowed to air dry. Genomic DNA was purified from a 1·2 mm punch of the dried blood spot using a two-step washing procedure as described by Zhou et al. (Reference Zhou, Hickford and Fang2006).
Milk sampling and phenotype measurement
All cows were milked twice daily and their daily milk yield (MY) in kilograms wet weight was recorded using Tru-test milk meters (Tru-test Ltd; New Zealand). Samples of milk (from both the morning and afternoon milkings) were collected once a month, over the 6 months from September to February (Spring to Summer). These samples were analysed for fat percentage (FP) and protein percentage (PP) using Fourier-Transform Infra-Red Spectroscopy (MilkoScan FT 120 Foss, Hillerød, Denmark). The averages of the FP and PP measurements were used for association analyses. A summary of the data is given in Table 1.
A summary of cow and milk trait data
Screening for variation in bovine UCP1
Variation in bovine UCP1 was screened for using PCR-SSCP. Primers (Table 1) were designed to amplify four regions of bovine UCP1 and based on the Bos taurus UMD 3·1 assembly (NCBI assembly accession AC_000174). These were Region-1 (upstream region), Region-2 (upstream region including the 5′-UTR), Region-3 (exon 1 and its flanking sequences), and Region-4 (exon 2 and its flanking sequences). These primers were synthesised by Integrated DNA Technologies (Coralville, IA, USA).
PCR amplification was performed in a 15-μl reaction containing the genomic DNA on one 1·2-mm punch of FTA paper, 0·25 µm of each primer, 150 µm of each dNTP (Bioline, London, UK), 2·5 mm of Mg2+, 0·5 U of Taq DNA polymerase (Qiagen, Hilden, Germany) and 1 × reaction buffer supplied. The thermal profile consisted of 2 min at 94·0 °C, followed by 35 cycles of 30 s at 94·0 °C, 30 s at the annealing temperature shown in Table 1 and 30 s at 72·0 °C, with a final extension of 5 min at 72·0 °C. Amplification was carried out in S1000 thermal cyclers (Bio-Rad, Hercules, CA, USA).
Amplicons were visualised by electrophoresis in 1% agarose gels (Quantum Scientific, Queensland, Australia), using 1 × TBE buffer (89 mm Tris, 89 mm boric acid, 2 mm Na2EDTA) containing 200 ng/ml of ethidium bromide.
A 0·7-μl aliquot of each amplicon was mixed with 7 µl of loading dye (98% formamide, 10 mm EDTA, 0·025% bromophenol blue, 0·025% xylene-cyanol). After denaturation at 95·0 °C for 5 min, the samples were rapidly cooled on wet ice and then loaded onto 16 cm × 18 cm, 14% acrylamide:bisacrylamide (37·5 : 1) (Bio-Rad) gels. Electrophoresis was performed in 0·5 × TBE buffer for 18 h using Protean II xi cells (Bio-Rad) under conditions specifically optimised for the individual gene regions being targeted (see Table 2). Gels were silver-stained according to the method of Byun et al. (Reference Byun, Fang, Zhou and Hickford2009).
PCR primers and PCR-SSCP conditions for the four regions of bovine UCP1
† Nucleotide coordinates refer to bovine sequence AC_0001740.
Sequencing and sequence analysis
For each region, PCR amplicons representing different SSCP banding patterns from cows that appeared to be homozygous were sequenced in both directions at Lincoln University, New Zealand. Sequence alignments, translations and comparisons were carried out using DNAMAN (version 5.2·10, Lynnon BioSoft, Vaudreuil, Canada).
Haplotype determination
Haplotypes spanning the two most variable regions (Region-2 and Region-4) were determined using the haplotyping approach described by Zhou et al. (Reference Zhou, Li, Liu, Wang, Luo and Hickford2014). This required the design of an additional four haplotype-specific primers (see Table 3), based on the sequence information revealed with the initial primers.
In this approach, each of the two regions was genotyped by PCR-SSCP. For those cows that were homozygous at both regions, or had only one heterozygous region, their haplotypes could be inferred directly. For those cows that were heterozygous in both regions, their haplotypes were deduced by haplotype-specific PCR using a three-primer (one forward and two reverse primers) approach. PCR amplifications were carried out using the conditions described above, but with the annealing temperatures shown in Table 3 and an extension time of 90 s.
Haplotype-specific PCR primers† and annealing temperatures used for resolving haplotypes of bovine UCP1
F, forward; R, reverse.
† Nucleotide mismatches are introduced to improve PCR specificity in the haplotyping procedure and these are underlined. A nucleotide introduced to reduce primer complementarity is shown in lower case.
Statistical analyses
All statistical analyses were performed using Minitab version 16 (Minitab Inc., PA, USA).
Two series of General Linear Models (GLMs) were used to assess the effect of the presence or absence of the common haplotypes (i.e. those with frequencies over 10%) on MY, FP and PP. The first series of GLMs were performed using single-haplotype presence/absence models to determine which haplotypes would be included in the second series of GLMs where all the haplotypes that had an association with P < 0·200 (and which could therefore potentially impact on the trait being tested), were included in the models.
For those haplotypes that were found to be associated with milk traits, pairwise comparisons between genotypes (once again only if the genotype frequency was over 10%) were then performed using a Tukey test.
The variables ‘cow age’ (counted in multiples of one year from birth) and ‘weeks in milk’ (counted in multiples of seven from the day of calving) affected many of the measured milk traits, so these were included in all of the GLMs as explanatory factors.
We were unable to factor ‘sire’ into the GLMs, as the majority of cows were of unknown paternity as sire genetics (semen straws) were purchased from a commercial semen producer on the basis of quantitative genetic evaluation of performance for key dairy traits. These straws often contain mixed semen, so individual sire identity is impossible to ascertain without access to all the seed stock. However, this weakness in the data can be partially off-set in two ways. Firstly, the cows studied ranged in age between 3 and 9 years and as it is unlikely that the same sires were used in the production of semen straws over this period, the cows are unlikely to share paternity. Secondly, if we concede that cows of the same age may share paternity, the inclusion of ‘cow age’ in the GLMs as an explanatory variable, would have at least partially accommodated this effect.
Results
Variation detected in bovine UCP1
For each region analysed, an amplicon of the expected size was obtained. The amplicons exhibited polymorphism upon SSCP analysis, with two PCR-SSCP patterns (named M and N) observed for Region-1 (upstream region), three patterns (named a, b and c) observed for Region-2 (upstream region), two patterns (named α and β) observed for Region-3 (exon 1 and its flanking sequences), and three patterns (named A, B and C) observed for Region-4 (exon 2 and its flanking sequences) (Fig. 1).
Variation in bovine UCP1. Four regions (R1, R2, R3 and R4) of the gene were screened for sequence variation (a). Two, three, two and three PCR-SSCP patterns (b) representing two, three, two and three different nucleotide sequences respectively (c) were detected in R1, R2, R3 and R4. Shaded boxes represent the exons (coding sequence) while the 5′UTR is shown as an unshaded box. Vertical lines represent the SNPs with the non-synonymous SNPs being indicated by *. The nucleotide numbering follows the HGVS nomenclature (http://www.hgvs.org/mutnomen). The coordinates for the four regions investigated refer to Table 2.
DNA sequencing and sequence analysis of PCR amplicons representative of the different SSCP patterns confirmed that these amplicons were derived from the expected gene regions and that the different SSCP patterns represented different nucleotide sequences. Three single nucleotide polymorphisms (SNPs) were detected in Region-1 (upstream region), four SNPs were detected in Region-2 (upstream region including the 5′-UTR), three SNPs were detected in Region-3 (exon 1 coding sequence and its flanking sequences) and three SNPs were detected in Region-4 (exon 2 and its flanking sequences).
There were two non-synonymous SNPs detected, one located in exon 1 (which would lead to the amino acid substitution p.Ala6Phe) and the other located in exon 2 (which would lead to the amino acid substitution p.Ile49Thr). In total, thirteen SNPs were identified in the regions of bovine UCP1 that were analysed.
Haplotypes identified across Region-2 to Region-4
Of the four regions investigated, Region-2 and Region-4 were the most variable gene regions, with each having three different sequences. These two regions were subsequently chosen to construct extended haplotypes. In total, six haplotypes were identified. The most common haplotype was a-B, followed by c-A and b-B. Haplotype c-C was rare and only occurred at a frequency of 0·4% in the cows investigated (Table 4).
Genotypes and haplotypes identified across Region-2 and Region-4 of UCP1 in the cows investigated
Effect of variation in UCP1 on milk production traits
Association analyses were carried out for haplotypes a-B, b-B and c-A which occurred at frequencies over 10% in the cows investigated.
The presence of haplotypes b-B and c-A was found to be associated with increased MY, and tended to be associated with decreased FP (Table 5). No associations were observed for haplotype a-B, or with other traits.
Association of presence or absence of UCP1 haplotypes with variation in milk production traits (mean ± SE)*
MY, milk yield; FP, fat percentage; PP, protein percentage.
* Predicted means and standard error from GLM including ‘cow age’ and ‘weeks in milk’ as factors. P < 0·05 in bold and 0·05 ≤ P < 0·1 in italics.
Of the four common genotypes (a-B/a-B, a-B/b-B, a-B/c-A and c-A/c-B), that had a frequency of over 10%; cows with genotype a-B/b-B produced more milk, but with a lower percentage of fat and protein, than cows with genotype a-B/a-B (Table 6). Cows of genotype a-B/c-A tended to produce more milk, but the milk produced by these cows had a lower fat content than cows of genotype a-B/a-B (Table 6).
Pairwise comparisons between genotypes of UCP1 on milk production traits (Mean ± se)*
MY, milk yield; FP, fat percentage; PP, protein percentage.
* Predicted means and se from GLM including ‘cow age’ and ‘weeks in milk’ as factors. P < 0·05 in bold and 0·05 ≤ P < 0·1 in italics.
Discussion
This study reports genetic variation in four regions of bovine UCP1. While SNPs in intron 1 of bovine UCP1 have been previously reported (Sonstegard & Kappes, Reference Sonstegard and Kappes1999), all the SNPs identified in this study were new. The presence of SNPs in this gene appears to be consistent with SNPs being reported for human UCP1 and sheep UCP1, and they are found in similar regions of the gene (Brondani et al. Reference Brondani, Assmann, Duarte, Gross, Canani and Crispim2012; Yuan et al. Reference Yuan, Liu, Liu, Qiao and Wu2012; Yang et al. Reference Yang, Forrest, Zhou, Hodge and Hickford2014). Overall the number of sequence variants identified for each region was low, when compared to the number of SNPs found; this being a consequence of the linkage of the SNPs. In total only six different haplotypes were identified across two of the variable regions. It is however likely that more SNPs and haplotypes may be identified, when more regions of the gene are analysed, and when more cows from different breeds are investigated.
Variation in bovine UCP1 was found to mainly affect MY, but it also had a minor effect on FP and PP. This is consistent with previous reports of moderate negative phenotypic correlations between MY, and FP and PP (Ahlborn & Dempfle, Reference Ahlborn and Dempfle1992). It suggests that selection for either MY, or FP and PP, may lead to deterioration in the other trait(s), and thus that higher milk yields will be typically associated with more diluted milk.
Haplotype b-B had an effect on milk traits, while haplotype a-B did not. As these haplotypes are identical in Region-4 (exon 2 and part of its intron flanking sequences), and are only different in Region-2, this suggests that the different effect of these two haplotypes is due to the variation in Region-2 of the gene. Given that the SNPs identified in Region-2 were located either within the 5′-UTR, or in sequences close to the putative TATA box, they may either affect RNA polymerase II binding and hence the expression of the gene, or the binding of the 30S subunit of the ribosomal RNA. Alternatively the SNPs identified in Region-2 may be in linkage disequilibrium with variation in other regions of the gene of functional or structural significance.
It is interesting to note that New Zealand (NZ) dairy cows investigated in this study had relatively low frequencies of haplotypes b-B and c-A, these occurring at frequencies of 11·1 and 21·7%, respectively. These haplotypes were found to favour increased MY, but they were unfavourable to FP and PP. The low frequency of these haplotypes in these cows would be consistent with the breeding objectives of the DairyNZ Breeding Worth (BW) genetic evaluation system (www.dairynz.co.nz/animal/animal-evaluation/about-nzael/), the system that is used in purchasing semen straws and in breeding replacement cows for the LUDF herd. This system defines a breeding index that is comprised of seven key traits. It places positive economic value on milk fat and protein percentage (i.e., FP and PP), but negative economic value on milk yield. The negative rating for milk yield is only lightly weighted in the index.
Under this breeding system, these haplotypes that favour milk yield may have been selected against. This is supported by the fact that the average milk yield of NZ dairy cows is low and ranked 37th in the world in 2011 (http://www.milkproduction.com/Library/Editorial-articles/World-Milk-Production/). This appears to be in agreement with the haplotype frequencies of UCP1 found in these NZ dairy cows, but it also strongly reflects the NZ dairy industry being primarily perennial ryegrass-white clover pasture-based, with animals being grazed outdoor year round, at low overall feeding allowance (Macdonald et al. Reference Macdonald, Verkerk, Thorrold, Pryce, Penno, McNaughton, Burton, Lancaster, Williamson and Holmes2008).
The low frequency of haplotypes favouring milk yield in NZ dairy cows may also be because UCP1 affects other traits that are of importance in the NZ dairy industry. Ablation of UCP1 induced obesity in mice (Feldmann et al. Reference Feldmann, Golozoubova, Cannon and Nedergaard2009) and in humans variation in UCP1 has been reported to be associated with body fat accumulation and body weight gain (Jia et al. Reference Jia, Tian, Cao, Tao, Zhang, Gao, Ge, Lin and Jois2010). In sheep variation in UCP1 has been associated with carcass fat grading score (Yang et al. Reference Yang, Forrest, Zhou, Hodge and Hickford2014), so taken together the findings of our research and the literature suggest there may be some genetic trade-off between body condition, or fat traits and milk yield.
It has been reported previously that dairy cow body condition score (BCS) has a negative genetic correlation with milk yield and this correlation becomes stronger as lactation progresses (Loker et al. Reference Loker, Bastin, Miglior, Sewalem, Schaeffer, Jamrozik, Ali and Osborne2012; Battagin et al. Reference Battagin, Sartori, Biffani, Penasa and Cassandro2013). This suggests that UCP1 may affect the ability of cows to store fat and thus cows that store less fat, produce more milk, and vice versa. As the ablation of UCP1 induces obesity in mice (Feldmann et al. Reference Feldmann, Golozoubova, Cannon and Nedergaard2009), it might be concluded that b-B and c-A haplotypes in cattle are less obesogenic than the other haplotypes, and thus more energy is used to increase milk yield and not stored as fat.
The thermogenic effect of UCP1 is another factor that may need to be considered. Although NZ does not have a hot climate, cows in NZ may be exposed to mild heat stress conditions at certain times in summer (Bluett et al. Reference Bluett, Fisher and Waugh2000). While the b-B and c-A haplotypes were associated with increased milk yield, they may make the cows more prone to heat stress and possibly also less likely to store fat if excess energy is available. Future studies should therefore not only investigate milk traits, but also BCS or some other measure of body fat storage (e.g. via ultrasound measurement of back-fat) and variation in cow temperature over lactation, and/or in warmer temperatures. Renna et al. (Reference Renna, Lussiana, Malfatto, Mimosi, Mattaglini, Darnhoefer and Grotzer2010) previously reported a decrease in Saturated Fatty Acid (SFA) content and an increase in unsaturated FA, monounsaturated FA, and polyunsaturated FA contents at high THI values, and more recently Hammami et al. (Reference Hammami, Vandenplas, Vanrobays, Rekik, Bastin and Gengler2015) reported that while most yield and fat traits had phenotypic declines as THI increased, C18 : 0, C18 : 1 cis-9, and four fatty acid groups (unsaturated FA, monounsaturated FA, polyunsaturated FA, and long-chain FA) increased with THI. C18 : 1 cis 9 levels are known to reflect body reserve mobilisation, suggesting a next step may be to investigate the effect of UCP1 variation on variation in milk FA constituents.
Financial support from the Ministry of Science and Innovation (DRCX 0802; Dairy Systems for Environmental Protection; NZ) and Lincoln University Gene-Marker Laboratory are acknowledged. We acknowledge the support of the AGMARDT Postdoctoral Fellowship to Long Cheng. We thank Guo Yang, Innocent Rugoho, Helen Hague, Thomas Vazhathara and Freeman Fang from Lincoln University (NZ) for technical support and useful discussion; and Peter Hancox from Lincoln University Dairy Farm for provision of cows and assistance with sampling.
