Affiliation
Department of Animal Sciences, D.H. Barron Reproductive and Perinatal Biology Research Program, and Genetics Institute, University of Florida, Gainesville, Florida, United States of America

Figures

Abstract

Heat stress compromises production, fertility, and health of dairy cattle. One mitigation strategy is to select individuals that are genetically resistant to heat stress. Most of the negative effects of heat stress on animal performance are a consequence of either physiological adaptations to regulate body temperature or adverse consequences of failure to regulate body temperature. Thus, selection for regulation of body temperature during heat stress could increase thermotolerance. The objective was to perform a genome-wide association study (GWAS) for rectal temperature (RT) during heat stress in lactating Holstein cows and identify SNPs associated with genes that have large effects on RT. Records on afternoon RT where the temperature-humidity index was ≥78.2 were obtained from 4,447 cows sired by 220 bulls, resulting in 1,440 useable genotypes from the Illumina BovineSNP50 BeadChip with 39,759 SNP. For GWAS, 2, 3, 4, 5, and 10 adjacent SNP were averaged to identify consensus genomic regions associated with RT. The largest proportion of SNP variance (0.07 to 0.44%) was explained by markers flanking the region between 28,877,547 and 28,907,154 bp on Bos taurus autosome (BTA) 24. That region is flanked by U1 (28,822,883 to 28,823,043) and NCAD (28,992,666 to 29,241,119). In addition, the SNP at 58,500,249 bp on BTA 16 explained 0.08% and 0.11% of the SNP variance for 2- and 3-SNP analyses, respectively. That contig includes SNORA19, RFWD2 and SCARNA3. Other SNPs associated with RT were located on BTA 16 (close to CEP170 and PLD5), BTA 5 (near SLCO1C1 and PDE3A), BTA 4 (near KBTBD2 and LSM5), and BTA 26 (located in GOT1, a gene implicated in protection from cellular stress). In conclusion, there are QTL for RT in heat-stressed dairy cattle. These SNPs could prove useful in genetic selection and for identification of genes involved in physiological responses to heat stress.

This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.

Funding: Research was supported by Agriculture and Food Research Initiative Competitive Grants 2010-85122-20623 and 2013-68004-20365 from the USDA National Institute of Food and Agriculture, and a grant from Southeast Dairy Inc. Milk Checkoff Program. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Heat stress in dairy cattle can compromise a variety of physiological functions including milk yield [1], reproduction [2] and immune function [3]. One strategy for reducing the magnitude of heat stress is to select animals that are genetically resistant to heat stress. Thermotolerant strains or breeds have been developed in many species of livestock [4]. Most of the negative effects of heat stress on animal performance are a consequence of either the physiological adaptations that homeotherms undergo to regulate body temperature or the adverse consequences of failure to regulate body temperature [4]. Thus, selection for regulation of body temperature during heat stress could result in animals that are resistant to deleterious effects of heat stress. In fact, two specific genes associated with effects of heat stress on milk yield, the slick gene [5] and an allele of ATP1A1[6], are also associated with lower rectal temperature (RT).

Core body temperature during heat stress is a heritable trait in dairy cattle, with early estimates varying from 0.15 to 0.31 [7] and a more recent estimate indicating a value of 0.17 [8]. Reliability of genetic estimates for rectal temperature, such as for other genetically-controlled traits [9], [10], should be improved by genome wide association studies (GWAS) to identify SNPs associated with regulation of rectal temperature. Quantitative trait loci (QTL) can be identified for lowly-heritability traits and used to improve reliability of genetic estimates despite the gain in reliability being less than for more heritable traits [11], [12]. In addition, GWAS can be useful for understanding the underlying biology of a trait by identifying candidate genes in physical proximity to QTL [11], [13], [14].

The objective of the current study was to perform a GWAS for rectal temperature during heat stress in lactating Holstein cows and identify SNPs that serve as QTL for that trait.

Materials and Methods

Use of animals was approved by the University of Florida Institutional Animal Care and Use Committee (Approval No. 2009-03578).

Data Collection

The study was conducted with lactating Holstein cows located in Florida dairies at nine locations in Bell (29.75′° N, 82.86° W), Hague (29.77° N, 82.42° W), Okeechobee (27.24° N, 80.83° W) and Trenton (29.61° N, 82.82°W). The following dairies allowed access to cows and records: Alliance Dairy (Trenton, Florida), Hilltop Dairy (Trenton, Florida), Larson Dairy (Okeechobee, Florida; three locations), McArthur Dairy (Okeechobee, Florida; two locations), North Florida Holsteins (Bell, Florida) and the University of Florida Dairy Unit (Hague, Florida). All except the University of Florida Dairy Unit are private. For seven locations, cows were housed in free-stall barns equipped with fans and sprinklers. In the other two locations, some cows were housed in free-stall barns and other cows were housed in tunnel-ventilation barns. Cows were fed total mixed rations. Data on RT were collected from a total of 11,128 cows. A subset of these data have been analyzed previously [8], [12].

Rectal temperature was measured between 1500 and 1700 h during the months of June to September in 2007 and between 1400 and 1600 h in the months of June to September in 2010, 2011 and 2012. Each cow was measured once during the experiment. Rectal temperature was measured using a digital GLA M750 thermometer (GLA Agricultural Electronics, San Luis Obispo, CA, USA). Measurements were recorded while cows were under shade and, generally, while in head locks (for free-stall barns) or while resting (for tunnel-ventilation barns).

Measurements of dry bulb temperature (Tdb), relative humidity (RH), dew point temperature (Tdp) and black globe temperature (Tbg) inside the barn and Tbg outside the barn were measured at 1 min intervals during the period of data collection using a HOBO-U12 data logger (Tdb, RH, and Tdp) and HOBO Water Temp Pro V2 data loggers (Tbg) (Onset Company, Bourne, MA, USA). Measurements were taken at a position that was 2 m from the ground at a position in the center of the barn where cows were housed. Rectal temperature was matched with measurements of Tdb, RH, Tdp, and Tbg to the nearest minute at which environmental variables were recorded. Temperature humidity index (THI) [15] was calculated as follows:where Tdb = dry bulb temperature (°C), and RH = relative humidity (%).

Data collected from days in which the THI at the time of collection was ≥78.2 were retained for genomic analysis. This THI was chosen because it is an indication that heat stress was sufficient to cause hyperthermia [16]. After restricting records to those cows undergoing heat stress and requiring that cows have known sire and dam, records were available from 4,447 cows sired by 220 bulls. The associated pedigree files included 12,346 animals. Of these phenotyped animals, 1,451 animals (107 cows and 1,344 bulls) had available genotypes using the Illumina BovineSNP50 BeadChip (Illumina, Inc., San Diego, CA, USA). Animals and SNP with call rates <0.90, SNP with minor allele frequencies <0.05, monomorphic SNP, and animals with parent-progeny conflicts were dropped from the dataset as part of the quality control process, leaving genotypes for 39,759 SNP from 1,440 individuals. A pedigree file including 12,346 male and female ancestors traced back to 1960 was obtained from the National Dairy Database maintained by the Animal Improvement Programs Laboratory (Beltsville, MD, USA).

Rectal temperatures were matched with performance data from the test-day closest to the date on which RT measurements were taken. Summary statistics for daily milk yield on the days when RT was measured, parity, DIM, and environmental variables are shown in Table 1.

(Co)variance Components

(Co)variance components were estimated by restricted maximum likelihood [17] using the REMLF90 program [18]. Due to the limited number of observations available, genomic data were included for (co)variance components estimation. The estimates of Dikmen et al. [8] were used as starting values, and convergence was declared when the change between successive rounds of iteration was <10−6.

The single-trait model used for ssGBLUP analysis of RT waswhere y is a vector of RT phenotypes, b is a vector of fixed effects, a is a vector of additive genetic effects, pe is a vector of permanent environmental effects, and e is a vector of random residual effects. Fixed effects included parity, year of data collection, stage of lactation (separated into three stage of lactation classes: DIM <100; DIM between 100 and 200; and DIM >200), location of cows at measurement (resting or in headlocks), farm type (tunnel ventilation and regular free-stall barn), technician who measured RT, and farm. THI and milk yield at the time of measurement were included as covariates. Xb, Za, and Zpe are incidence matrices relating a record to fixed environmental effects in b and to random animal and permanent environmental effects in a and pe, respectively. Results were calculated using a Bayes A model as implemented in the BLUPF90 program [18] modified for genomic analyses [22]. Genome-wide association analyses for RT were conducted with individual SNP effects, as well as moving averages of 2, 3, 4, 5, and 10 consecutive SNP, using the POSTGSF90 package (I. Aguilar, INIA, Las Piedras, Canelones, Uruguay, 2012, personal communication).

Results and Discussion

Identification of QTL Associated with RT

Wang et al. [19] found in a simulation study that the correlation of QTL with the sum of adjacent SNPs increased up to 8 SNP, and then decreased as the number of summed SNP increased. This is because the closest SNP to a QTL is not always the best predictor of the QTL effect [25]. In the current study, 2, 3, 4, 5, and 10 adjacent SNP were averaged in order to identify consensus genomic regions associated with RT. A GWAS based on ssGBLUP was used in place of more traditional GWAS approaches because the latter do not directly use phenotypes of non-genotyped individuals. P-values were not used to declare regions as significant because such values are difficult to define and compare using classical hypothesis tests when shrunken estimators such as PTA are used [26]. Recent studies also have found that traditional GWAS methods often produce large numbers (e.g., hundreds) of significant effects, most of which cannot be validated in follow-up studies [27]. The principal objective of this study was to identify genomic regions associated with regulation of RT, and this is not dependent on tests of statistical significance that are more appropriate for candidate gene-based studies.

The proportion of SNP variance explained by each marker individually is shown in Figure S1 and File S1. The individual SNP results were very noisy, and only one region on BTA 26 showed a clear signal. The use of 2-SNP windows (Figure S2) resulted in smoother marker effects, but it was still difficult to differentiate between markers with large effects and those with small effects. The Manhattan plots from the 3- and 5-SNP windows (Figures 1 and 2) resulted in clearer, and generally consistent patterns. Results were similar for the 4-SNP windows (Figure S3).

The 10-SNP results (Figure S4) are similar to the results found using narrower windows. Nonetheless the 10-SNP analysis poses problems because the Illumina BovineSNP50 BeadChip has an average intermarker spacing of 49.4 kbp so that each point on the resulting Manhattan plot covers an average of 494 kbp. This is problematic because useful linkage dsequilibrium (LD) in the cow extends for less than 100 kb [28]. Genes or groups of genes that are not in LD with one another but which have large effects should be represented by separate peaks.

The 20 largest explanatory loci for RT for each of the analyses is listed in Tables 2 and 3 and S1 through S4. Most commonly, the largest proportion of variance was explained by markers flanking the region between 28,877,547 and 28,907,154 bp on Bos taurus autosome (BTA) 24. That region is flanked by a U1 spliceosomal RNA (U1) on the left (28,822,883 to 28,823,043 bp) and a cadherin-2 (NCAD) gene on the right (28,992,666 to 29,241,119 bp). The U1 small ribonucleoprotein is involved in postranscriptional modification and regulation of mRNA length [29], both of which could be related to changes in gene expression in cells exposed to elevated temperature [30], [31]. NCAD was more highly expressed in paratuberculosis-infected cattle than in noninfected animals [32], but it is not clear what relationship, if any, there is between NCAD and stress responses.

Two regions of interest were flagged on BTA 16. The 2- and 3-SNP analyses (Table S2 and 2) found that the SNP at 58,500,249 bp explained 0.08 and 0.11% of the SNP variance, respectively. That contig includes a small nucleolar RNA (SNORA19) from 58,520,021 to 58,520,149 bp; a ubiquitin-protein ligase (RFWD2) from 58,600,678 to 58,838,844 bp; and small Cajal body specific RNA 3 (SCARNA3) from 58,628,128 to 58,628,268 bp. Small nucleolar RNA like SNORA19 identify sites of pseudouridine formation [33], which may be involved in the initiation of translation. The SCARNA3 gene also encodes a small nucleolar RNA similar in structure to SNORA19. The RFWD12 gene encodes an E3 RING-type ubiquitin-protein ligase, which selects proteins for proteasomal degradation [34]. Many abnormal proteins are proteolyzed by the ubiquitin system as part of the stress response in eukaryotes [35]. Moreover, a total of 28 genes identified as being regulated by heat shock in the bovine embryo are associated with ubiquitin C [31].

The second region of interest on BTA 16 is located at about 35,272,426 bp. The closest genes to this area are centrosomal protein of 170 kDa (CEP170) and inactive phospholipase D5 (PLD5). Neither gene has an obvious relationship with physiological response to heat stress.

A region of BTA 5 at approximately 89,500,000 bp was consistently identified across analyses. Two genes flank the consensus region, solute carrier organic anion transporter family member 1C1 (SLCO1C1) and a phosphodiesterase (PDE3A). In humans, SLCO1C1 mediates the Na+-independent high-affinity transport of thyroxine and reverse triiodothyronine [36]. Perhaps this gene is involved in the mechanism which depresses plasma thyroxine concentrations in heat-stressed dairy cows [37]. Human heat shock protein 70 has been shown to increase the enzymatic activity of phosphodiesterase in heat-shocked cells [38] and there may be a similar association in cattle.

The 4- and 5-SNP analyses showed a peak near 64,400,000 bp on BTA 4 that explained 0.09–0.11% of the observed variance (Tables S3 and 3). This region includes the genes kelch repeat and BTB domain-containing protein (KBTBD2) and U6 snRNA-associated Sm-like protein LSM5 (LSM5). KBTBD2 encodes for a protein that participates in ubiquitination of proteins [39] while LSM5 is likely to participate in RNA processing and to form part of the stress granule seen in stressed-cells that contains mRNAs stalled in translation [40].

A region centered on the SNP at 20,290,497 bp on BTA 26 explained 0.10% to 0.22% of the observed variance in the 2-SNP to 5-SNP analyses (Tables 2 and 3 and S1–S3. The SNP is located in the glutamine-oxaloacetic transaminase, soluble (GOT1) gene. GOT1 is involved in synthesis of sulfur dioxide, which has been shown to reduce formation of reactive oxygen products and protect rat myocardial tissue in response to isoproterenol [41].

In an earlier study, Hayes et al. [42] identified a SNP on BTA 29 (BFGL-NGS-30169) that was associated with genetic variation in effects of heat stress on milk yield in Holsteins and Jerseys. A nearby SNP on BTA 29 (ARS-BFGL-NGS-107395 at 47,527,067 bp) was one of the top 30 markers in the 2- through 5-SNP analyses, ranking between the 19th and 35th largest SNP. It is located 355,843 bp upstream of the SNP identified by Hayes et al. [42] in a contig that includes only an annotated small nucleolar RNA (SNORD14).

There are small but significant genetic correlations of RT during heat stress with 305-d milk, fat, and protein yields, productive life, net merit, somatic cell score and daughter pregnancy rate [8]. Regardless of the analysis used to identify SNPs explanatory for RT, none of the 10 largest explanatory loci for RT were in common with 1,586 SNP markers related to 31 other traits in dairy cattle previously identified by Cole et al. [11]. Thus, it might be possible to use SNPs associated with RT to select for thermotolerance without inadvertent selection for other traits.

In conclusion, these data demonstrate that QTL exist to predict some of the genetic variation in RT. These QTL may prove important for genetic selection of thermotolerance. The challenge of performing GWAS for traits with low heritability when the number of phenotyped animals is low means that confirmation of the utility of these QTL is essential. In addition, several candidate genes for regulation of RT in dairy cattle were identified and one or more of these may play an important role in physiological adaptation to heat stress. Only one candidate gene, SLC01C1, which is involved in regulation of metabolic rate through transport of thyroxine, plays a known role in processes controlling body temperature. More commonly, candidate genes play roles that are important for stabilizing cellular function during stress. Among these are GOT1, which synthesizes the cytorotective compound sulphur dioxide, genes involved in protein ubiquitination (KBTBD2 and RFWD12), and genes involved in RNA metabolism (LSM5, SCARNA3, SNORA19, and U1).

The name, chromosome number, SNP location (BP), and proportion of SNP variance for rectal temperature explained by each SNP from a single-step GBLUP analysis in CSV format. Marker locations are based on the UMD 3.1 assembly of the bovine genome.