Advertisement

Derivation and genome-wide association study of a principal component-based measure of heat tolerance in dairy cattle

Open ArchivePublished:March 29, 2017DOI:https://doi.org/10.3168/jds.2016-12249

      ABSTRACT

      Heat stress represents a key factor that negatively affects the productive and reproductive performance of farm animals. In the present work, a new measure of tolerance to heat stress for dairy cattle was developed using principal component analysis. Data were from 590,174 test-day records for milk yield, fat and protein percentages, and somatic cell score of 39,261 Italian Holstein cows. Test-day records adjusted for main systematic factors were grouped into 11 temperature-humidity index (THI) classes. Daughter trait deviations (DTD) were calculated for 1,540 bulls as means of the adjusted test-day records for each THI class. Principal component analysis was performed on the DTD for each bull. The first 2 principal components (PC) explained 42 to 51% of the total variance of the system across the 4 traits. The first PC, a measure of the level at which the curve is located, was interpreted as a measure of the level at which the DTD curve was located. The second PC, which shows the slope of increasing or decreases DTD curves, synthesized the behavior of the DTD pattern. Heritability of the 2 component scores was moderate to high for level across all traits (range = 0.23–0.82) and low to moderate for slope (range = 0.16–0.28). For each trait, phenotypic and genetic correlations between level and slope were equal to zero. A genome-wide association analysis was carried out on a subsample of 423 bulls genotyped with the Illumina 50K bovine bead chip (Illumina, San Diego, CA). Two single nucleotide polymorphisms were significantly associated with slope for milk yield, 4 with level for fat percentage, and 2 with level and slope of protein percentage, respectively. The gene discovery was carried out considering windows of 0.5 Mb surrounding the significant markers and highlighted some interesting candidate genes. Some of them have been already associated with the mechanism of heat tolerance as the heat shock transcription factor (HSF1) and the malonyl-CoA-acyl carrier protein transacylase (MCAT). The 2 PC were able to describe the overall level and the slope of response of milk production traits across increasing levels of THI index. Moreover, they exhibited genetic variability and were genetically uncorrelated. These features suggest their use as measures of thermotolerance in dairy cattle breeding schemes.

      Key words

      INTRODUCTION

      The improvement of an animal's ability to cope with adverse environmental conditions is one of the great challenges of animal breeding for the future (
      • Bernabucci U.
      • Lacetera L.H.
      • Baumgard R.P.
      • Rhoads B.
      • Ronchi B.
      • Nardone A.
      Metabolic and hormonal acclimation to heat stress in domesticated ruminants.
      ). Among the traits that contribute to define animal adaptability to environmental variation, tolerance to heat stress plays a major role. Heat stress can be defined as the condition where the animal is not able to adequately dissipate the excess of endogenous or exogenous heat to maintain body thermal balance (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effect of heat stress in Italian Holstein dairy cattle.
      ). In dairy cattle, it is known that heat stress results in relevant economic losses due to reduced milk production and reproduction performance (
      • Aguilar I.
      • Misztal I.
      • Tsuruta S.
      Short communication: Genetic trends of milk yield under heat stress for US Holsteins.
      ;
      • Nardone A.
      • Ronchi B.
      • Lacetera N.
      • Ranieri M.S.
      • Bernabucci U.
      Effects of climate changes on animal production and sustainability of livestock systems.
      ;
      • Biffani S.
      • Bernabucci U.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      Short communication: Effect of heat stress on nonreturn rate of Italian Holstein cows.
      ). Increasing concern about tolerance to heat stress for dairy animals in temperate areas is a consequence of both climate change and higher metabolic heat production by high-yielding animals (
      • Kadzere C.T.
      • Murphy M.R.
      • Silanikove N.
      • Maltz E.
      Heat stress in lactating dairy cows: a review.
      ;
      • Hansen P.J.
      Exploitation of genetic and physiological determinants of embryonic resistance to elevated temperature to improve embryonic survival in dairy cattle during heat stress.
      ;
      • Segnalini M.
      • Nardone A.
      • Bernabucci U.
      • Vitali A.
      • Ronchi B.
      • Lacetera N.
      Dynamics of the temperature-humidity index in the Mediterranean basin.
      ).
      If tolerance to heat stress is a quite straightforward concept, its systematic measure remains problematic. On the other hand, a quantification of this trait is fundamental if it is to be considered a potential selection goal in breeding programs.
      Some physiological traits are related to the ability of the animal to cope with heat stress. For example, rectal temperature and respiration rate increase when animals are exposed to warm environment (
      • Dikmen S.
      • Cole J.B.
      • Null D.J.
      • Hansen P.J.
      Heritability of rectal temperature and genetic correlations with production and reproduction traits in dairy cattle.
      ;
      • Perano K.M.
      • Usack J.G.
      • Angenent L.T.
      • Gebremedhin K.G.
      Production and physiological responses of heat-stressed lactating dairy cattle to conductive cooling.
      ;
      • Garner J.B.
      • Douglas M.L.
      • Williams S.R.O.
      • Wales W.J.
      • Marett L.C.
      • Nguyen T.T.T.
      • Reich C.M.
      • Hayes B.J.
      Genomic selection improves heat tolerance in dairy cattle.
      ). These traits exhibit a genetic component; for example, a moderate heritability and associations with SNP and candidate genes have been reported for rectal temperature (
      • Dikmen S.
      • Cole J.B.
      • Null D.J.
      • Hansen P.J.
      Genome-wide association mapping for identification of quantitative trait loci for rectal temperature during heat stress in Holstein cattle.
      ,
      • Dikmen S.
      • Wang X.Z.
      • Ortega M.S.
      • Cole J.B.
      • Null D.J.
      • Hansen P.J.
      Single nucleotide polymorphisms associated with thermoregulation in lactating dairy cows exposed to heat stress.
      ). However, the inclusion of these heat tolerance indicator traits in large-scale phenotype recording systems for selecting thermotolerant animals appears rather problematic in terms of logistics and costs. An alternative is to evaluate heat tolerance by measuring changes of milk production traits under warm environmental conditions (
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      ;
      • Carabaño M.J.
      • Logar B.
      • Bormann J.
      • Minet J.
      • Vanrobays M.L.
      • Diaz C.
      • Tychon B.
      • Gengler N.
      • Hammami H.
      Modeling heat stress under different environmental conditions.
      ;
      • Nguyen T.T.T.
      • Bowman P.J.
      • Haile-Mariam M.
      • Pryce J.E.
      • Hayes B.J.
      Genomic selection for tolerance to heat stress in Australian dairy cattle.
      ). In dairy cattle populations involved in selection programs, milk production data could be easily retrieved from dairy recording systems and associated with climate data provided by weather stations. The variable most frequently used to evaluate heat stress conditions is the temperature-humidity index (THI). The approach commonly used to evaluate heat tolerance relies on the so-called broken line model (
      • Ravagnolo O.
      • Misztal I.
      Genetic component of heat stress in dairy cattle, parameter estimation.
      ). It assumes the existence of a comfort zone limited by an upper threshold value (TH0), beyond which the production linearly decreases as THI increases (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effect of heat stress in Italian Holstein dairy cattle.
      ;
      • Carabaño M.J.
      • Bachaga K.
      • Ramon M.
      • Diaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: selection tools.
      ).
      In statistical models, tolerance to heat stress might be fitted according to a reaction norm model (
      • Kolmodin R.
      • Bijma P.
      Response to mass selection when the genotype by environment interaction is modeled as a linear reaction norm.
      ), where the phenotype is expressed as a linear function of an environmental variable (for example THI or temperature). Very often, the environmental variable effect is a dummy variable, set to zero when THI < TH0 and to THI − TH0 when THI > TH0 (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effect of heat stress in Italian Holstein dairy cattle.
      ). Some studies adopted the fixed value of 72 for TH0 (
      • Ravagnolo O.
      • Misztal I.
      Genetic component of heat stress in dairy cattle, parameter estimation.
      ;
      • Bohmanova J.
      • Misztal I.
      • Tsuruta S.
      • Norman H.D.
      • Lawlor T.J.
      Short communication: Genotype by Environment Interaction Due to Heat Stress.
      ;
      • Aguilar I.
      • Misztal I.
      • Tsuruta S.
      Genetic components of heat stress for dairy cattle with multiple lactations.
      ), but recently different TH0 have been estimated across traits, parities and geographical regions (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effect of heat stress in Italian Holstein dairy cattle.
      ;
      • Biffani S.
      • Bernabucci U.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      Short communication: Effect of heat stress on nonreturn rate of Italian Holstein cows.
      ).
      Some studies on tolerance of heat stress have used individual production curves along different THI levels corrected for fixed factors as a measure of heat tolerance (
      • Hayes B.J.
      • Bowman P.J.
      • Chamberlain A.J.
      • Savin K.
      • van Tassell C.P.
      • Sostengard T.S.
      • Goddard M.E.
      A validated genome wide association study to breed cattle adapted to an environment altered by climate changes.
      ;
      • Carabaño M.J.
      • Logar B.
      • Bormann J.
      • Minet J.
      • Vanrobays M.L.
      • Diaz C.
      • Tychon B.
      • Gengler N.
      • Hammami H.
      Modeling heat stress under different environmental conditions.
      ). Average curves of bull progeny for milk production traits across different THI levels, named as daughter trait deviations (DTD), have been recently used as phenotypes in a genomic selection study on tolerance to heat stress (
      • Nguyen T.T.T.
      • Bowman P.J.
      • Haile-Mariam M.
      • Pryce J.E.
      • Hayes B.J.
      Genomic selection for tolerance to heat stress in Australian dairy cattle.
      ).
      For genetic purposes, individual effects for heat tolerance are usually fitted with an intercept and a slope, representing the overall level of production and the response of the animal to heat stress, respectively. Main concerns about these approaches are on the use of a common threshold across all animals and the assumption of linearity for the production decay after TH0 (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effect of heat stress in Italian Holstein dairy cattle.
      ;
      • Carabaño M.J.
      • Bachaga K.
      • Ramon M.
      • Diaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: selection tools.
      ). On the other hand, estimation of individual thresholds (
      • Sánchez J.P.
      • Misztal I.
      • Aguilar I.
      • Zumbach B.
      • Rekaya R.
      Genetic determination of the onset of heat stress on daily milk production in the US Holstein cattle.
      ) is more realistic though it is more computationally demanding. Individual change points of production patterns for increasing THI levels have been fitted also with Legendre polynomials in random regression models (
      • Brügemann K.
      • Gernand E.
      • von Borstel U.U.
      • König S.
      Genetic analyses of protein yield in dairy cows applying random regression models with time-dependent and temperature x humidity-dependent covariates.
      ;
      • Carabaño M.J.
      • Bachaga K.
      • Ramon M.
      • Diaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: selection tools.
      ,
      • Carabaño M.J.
      • Logar B.
      • Bormann J.
      • Minet J.
      • Vanrobays M.L.
      • Diaz C.
      • Tychon B.
      • Gengler N.
      • Hammami H.
      Modeling heat stress under different environmental conditions.
      )
      Several papers that evaluated the effect of heat stress on milk reported an unfavorable genetic relationship between production and heat tolerance (
      • Sánchez J.P.
      • Misztal I.
      • Aguilar I.
      • Zumbach B.
      • Rekaya R.
      Genetic determination of the onset of heat stress on daily milk production in the US Holstein cattle.
      ;
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effect of heat stress in Italian Holstein dairy cattle.
      ;
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      ). These results were confirmed also by the strong negative correlations (−0.85 and −0.75) between genomic breeding value for milk DTD-derived heat tolerance and EBV for milk yield in Australian Holsteins and Jerseys respectively (
      • Nguyen T.T.T.
      • Bowman P.J.
      • Haile-Mariam M.
      • Pryce J.E.
      • Hayes B.J.
      Genomic selection for tolerance to heat stress in Australian dairy cattle.
      ). Such correlations are the result of the increased metabolic heat production that occurs in high-producing cows and that exacerbate the effects of the external heat. This represents a severe constraint to an efficient selection for improving heat tolerance without negative consequences on production. The aggregation of the 2 traits into a selection index may help selection, even though the definition of optimal economic weights could remain a theoretical issue and the negative correlation undoubtedly will reduce the selection response on each individual trait. An alternative could be the use of a measure of tolerance to heat stress that is not correlated with production levels. The use of a model-free approach, able to disentangle main features of DTD without imposing specific constraints, is an appealing option for assessing proper variables to study tolerance to heat stress. Principal component analysis (PCA) is a multivariate statistical technique able to synthesize complex patterns as the lactation curves for dairy traits in 2 variables with a clear technical meaning (
      • Macciotta N.P.P.
      • Vicario D.
      • Cappio-Borlino A.
      Use of multivariate analysis to extract latent variables related to level of production and lactation persistency in dairy cattle.
      ,
      • Macciotta N.P.P.
      • Gaspa G.
      • Bomba L.
      • Vicario D.
      • Dimauro C.
      • Cellesi M.
      • Ajmone-Marsan P.
      Genome-wide association analysis in Italian Simmental cows for lactation curve traits using a low-density (7K) SNP panel.
      ). Principal component analysis can, therefore, be conveniently used to analyze DTD curves for extracting new variables able to synthesize the pattern.
      In the present work, a PCA approach was tested to derive indicator variables of tolerance to heat stress from milk production data in dairy cattle. Moreover, a genome-wide association study (GWAS) using a medium-density (50K) SNP panel was used for investigating the genetic determinism of these new variables.

      MATERIALS AND METHODS

      Data

      Data were 590,174 test-day (TD) records for milk yield (MY), fat (FP) and protein (PP) percentages, and SCS (
      • Ali A.K.A.
      • Shook G.E.
      An optimum transformation for somatic cell concentration in milk.
      ) of 39,261 Italian Holstein cows (first, second, and third parity) from 484 farms, collected from 2001 to 2007. Data were recorded by the Italian Breeders Association according to International Committee for Animal Recording standards (http://www.icar.org/index.php/publications-technical-materials/recording-guidelines/). Age at calving classes were established for each parity according to the following thresholds: 20 to 36 (17 classes), 31 to 50 (20 classes), and 42 to 65 (24 classes) months of age for first-, second-, and third-parity cows, respectively. All cows had first-lactation data and a minimum of 8 TD records per lactation (from 5–305 DIM). A minimum of 24 records per herd-year of calving were required. Cows were sired by 4,184 AI bulls.
      Daily weather information were collected from 35 meteorological stations located no more than 5 km from the considered herd. The THI index (
      • Kelly C.F.
      • Bond T.E.
      Bioclimatic factors and their measurements.
      ) was then calculated as:
      THI=(1.8×AT+32)(0.550.55×RH)×[(1.8×AT+32)58],


      where AT is the maximum daily temperature, expressed in degrees Celsius, and RH is the minimum relative humidity, expressed as a percentage.

      Statistical Analysis

      Test-day records were first analyzed with the following mixed linear model
      y=month(year)+age+DIM×parity+herd(year)+e,
      [1]


      where y is the record for MY, FP, PP, or SCS; month(year) is the fixed effect of the month of calving (12 mo) nested within the year of calving (7 yr, 2001–2007); age is the fixed effect of age class in months (61 classes, from 20 to 65 mo); DIM × parity is the interaction between the fixed effect of the DIM class (10 intervals of 30 d each) and the fixed effect of parity (3 parities, 1–3); herd(year) is the random effect of the herd (458) nested within calving year; and e his the random residual.
      Residuals of model [1] are therefore production data adjusted for main systematic factors, except from additive genetic and THI effects. On the basis of THI values, records were grouped into 11 THI classes (1 = 50–52, 2 = 53–54, …, 11 = >79). Distribution of records across THI classes is reported in Figure 1. Means of residuals were calculated for each bull and THI class for obtaining DTD (
      • Nguyen T.T.T.
      • Bowman P.J.
      • Haile-Mariam M.
      • Pryce J.E.
      • Hayes B.J.
      Genomic selection for tolerance to heat stress in Australian dairy cattle.
      ). The DTD plotted against the THI class express the sensitivity of bull's daughter production performance for increasing THI levels.
      Figure thumbnail gr1
      Figure 1Distribution of test-day records across different classes of temperature-humidity index (THI).
      The following step was to derive a measure able to summarize the shape of these curves that could be used as dependent variable in a GWAS. A PCA was carried out on the 11 points of the DTD curves, considered as different variables; for example, the first bull for milk yield had 11 records (i.e., DTD1_MY, DTD2_MY, …, DTD11_MY). Only bulls that had the complete set of 11 DTD values for each trait were considered (1,540 for MY, 1,513 for FP, 1,536 for PP, and 1,535 for SCS, respectively). The PCA was carried out using the SAS PRINCOMP procedure (
      • SAS Institute
      ). The number of principal components to be retained was based on their eigenvalue, and on their relationships with the original variables.
      Principal component (PC) scores were then calculated for each bull and treated as new phenotypes for performing either genetic parameter estimation or GWAS. Variance components for PC scores were estimated with the following multi-trait animal model:
      y=μ+animal+e,


      where y is a vector of PC scores for MY, FP, PP, and SCS, respectively; μ is the overall mean; animal is the random additive genetic effect; and e is the residual term. The following (co)variance structure was assumed for the random effects:
      var|ae|=|AG000IR0|


      where A is the numerator relationship matrix, G0 is the matrix of (co)variances for additive effects, R0 is a diagonal matrix of residual variances corresponding to each trait, and I is an identity matrix. The pedigree file had 21,685 animals, including the 1,540 sires with DTD in the data set.
      The model was solved using the program AIREML90 (
      • Misztal I.
      • Tsuruta S.
      • Strabel T.
      • Auvray B.
      • Druet T.
      • Lee D.H.
      BLUPF90 and related programs (BGF90).
      ). Considering that PC scores were calculated starting from average yields per bull, the ratio was
      h2=σA2σA2+σE2,


      where σA2 and
      σE2


      are the additive genetic and the residual variances, respectively; h2 represents an approximation of the true heritability because averaging affects the variability of the response (different number of TD records per bull). Thus, obtained values have been properly called as pseudo-heritability.
      Of the 1,540 bulls considered for the PC score calculation, 423 were genotyped with the Illumina 50K bovine bead chip (Illumina, San Diego, CA). Monomorphic SNP (7,140) and SNP with a call rate <95% (1,045) were discarded. In total, 45,546 SNP were retained for the analysis. Genome-wide scan was performed on PCA scores with the GenABEL R package (
      • Aulchenko Y.S.
      • Ripke S.
      • Isaacs A.
      • van Duijn C.M.
      GenABEL: An R package for genome-wide association analysis.
      ), using the GRAMMAR procedure. First, an additive polygenic model was fitted to obtain individual residuals using the genomic relationship matrix. Then, SNP association was tested using a linear model on residuals of the first step. The SNP statistical significance was corrected for the stratification of the population using the genomic control (GC) option (
      • Amin N.
      • van Duijn C.M.
      • Aulchenko Y.S.
      A genomic background based method for association analysis in related individuals.
      ). The GC-corrected P-values (GC_Pi) were further corrected for multiple testing using either (1) the Bonferroni correction, obtained as GC_Pi × m (where m is the number of performed tests); (2) and calculating the false discovery rate (FDR), as (GC_Pi × m)/m0, where m0 is the number of tests having the GC P-values lower or equal to GC_Pi. A marker was declared significantly associated with a trait when the FDR was <0.10.
      Gene discovery analysis was carried out considering windows of 0.5 Mb surrounding the significant marker (0.25 Mb up- and downstream, respectively). Genes were derived from UCSC Genome Browser Gateway (http://genome.ucsc.edu/). Both SNP and gene positions were obtained from the UMD3.1 bovine genome assembly (
      • Zimin A.V.
      • Delcher A.L.
      • Florea L.
      • Kelley D.R.
      • Schatz M.C.
      • Puiu D.
      • Hanrahan F.
      • Pertea G.
      • Van Tassell C.P.
      • Sonstegard T.S.
      • Marçais G.
      • Roberts M.
      • Subramanian P.
      • Yorke J.A.
      • Salzberg S.L.
      A whole-genome assembly of the domestic cow.
      ).

      RESULTS

      PCA

      Eigenvectors and eigenvalues of the first 2 PC extracted from all the 4 considered phenotypes are reported in Table 1. The first 2 PC explained from 42 to 51% of the total variance of the system across the 4 traits. The choice of retaining only the first 2 PC was motivated by the magnitude of single eigenvalues, even though the amount of explained variance was not particularly relevant. A common criterion used for retaining PC is that the eigenvalue should be greater than 1. In the present work, for all 4 traits only the first eigenvalue fulfilled this requirement, and the second eigenvalue was very close to 1 (Table 1).
      Table 1Eigenvectors and eigenvalues of the first 2 principal components (level and slope) extracted from the correlation matrix of daughter trait deviations for milk yield (DTD_MY), fat (DTD_FP) and protein (DTD_PP) percentages, and SCS (DTD_SCS)
      Temperature-humidity index intervalDTD_MYDTD_FPDTD_PPDTD_SCS
      LevelSlopeLevelSlopeLevelSlopeLevelSlope
      50–520.250.450.240.760.26−0.480.220.49
      53–550.280.410.290.370.29−0.370.250.55
      56–580.320.220.300.170.31−0.200.330.07
      59–610.310.140.31−0.250.32−0.070.300.12
      62–640.320.170.310.030.32−0.180.340.09
      65–670.320.000.31−0.020.33−0.060.300.01
      68–700.32−0.040.32−0.140.310.100.34−0.11
      71–730.32−0.130.32−0.280.310.220.31−0.06
      74–760.31−0.310.31−0.020.310.110.33−0.26
      77–790.30−0.340.31−0.300.290.360.32−0.14
      >790.25−0.550.28−0.110.250.600.27−0.57
      Eigenvalue4.430.964.250.914.790.913.600.98
      Eigenvalue %409398438339
      The first principal component (PC1) showed positive and moderate eigenvector coefficients or loadings (ranging from 0.25 to 0.35) with all the original variables (Table 1). Thus, PC1 can be considered as a measure of the level at which the curve is located and it was named level. Bulls with large scores for this component have their DTD pattern located on a higher level. The second principal component (PC2) exhibited larger loadings (up to 0.76 for DTD_FP) with both positive and negative signs. In particular, PC2 showed positive values for the first part of THI interval and negative for the second part for MY, SCS, and, even if less definite, for FP; on the contrary, PC2 for PP showed the opposite trend. For its structure, PC2 defined the shape of the DTD curve. Therefore, larger values (or smaller in the case of protein content) of PC2 scores characterize DTD curves with a decreasing pattern, whereas smaller (or larger) PC2 scores indicate increasing patterns. This component was named slope. The interpretation of level and slope meaning may be better inferred from the average DTD curves for different PC1 and PC2 classes (Figures 2). Only MY data were reported for brevity, but the other traits showed the same pattern.
      Figure thumbnail gr2
      Figure 2Average curves of daughter trait deviations (DTD) for milk yield of groups of bulls of different level (a) and slope (b) score classes (♦ = <−2; ▪ = −2 to −1; Δ = −1 to 0; • = 0 to 1; continuous line = 1 to 2; dotted line = >2). Points are plotted for the average DIM on each test day. THI = temperature-humidity index.
      To simplify the comparison between slope for different traits, scores of slope for PP were multiplied by −1. Pearson correlations among PC scores (Table 2) confirm the expected orthogonality between level and slope within each trait. Sign and magnitude of correlations between PC scores of different traits confirm the meaning of the new variables extracted; for example, the negative correlations between level for milk yield and for fat and protein percentage and the positive correlation between the last 2 traits.
      Table 2Pearson correlations between the scores of the first principal components (level and slope) extracted from the correlation matrix of daughter trait deviations for milk yield (MY), fat (FP) and protein (PP) percentage, and SCS
      ItemLevel MYSlope MYLevel FPSlope FPLevel PPSlope PPLevel SCSSlope SCS
      Level MY1.000.00−0.37−0.03−0.340.00−0.07−0.07
      Slope MY1.000.01−0.12−0.01−0.20−0.04−0.16
      Level FP1.000.000.520.02−0.050.04
      Slope FP1.000.030.23−0.020.05
      Level PP1.000.000.050.02
      Slope PP1.000.02−0.15
      Level SCS1.000.00
      Slope SCS1.00

      Genetic Parameter Estimation

      All new variables exhibited genetic variability (Table 3). In particular, pseudo-heritability had a moderate to high level across all traits (values ranged between 0.32 for SCS and 0.82 for FP, respectively) and low to moderate slope (ranging from 0.16 for SCS to 0.28 for PP, respectively). Moreover, values for slope of FP and PP were consistent with the proportion of variability explained by one of the canonical variables obtained for the eigendecomposition of the additive (co)variance matrix of random regression coefficients for these traits (
      • Carabaño M.J.
      • Bachaga K.
      • Ramon M.
      • Diaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: selection tools.
      ). Pseudo-genetic correlations confirm at a genetic level the substantial orthogonality between the level and slope components within each trait. Values of pseudo genetic correlation rg between level values across the different traits were large (absolute values >0.65), with the exception of comparisons involving SCS. Large pseudo-genetic correlations were also observed between slope values for all traits, with the exception of the correlation between FP and SCS (Table 3).
      Table 3Pseudo-heritability (on the diagonal) and pseudo-genetic correlations (out of the diagonal) between the scores of the first principal components extracted from the correlation matrix of daughter trait deviations for milk yield (MY), fat (FP) and protein (PP) percentage, and SCS (SE of pseudo-heritabilities in parentheses)
      ItemLevel MYSlope MYLevel FPSlope FPLevel PPSlope PPLevel SCSSlope SCS
      Level MY0.52 (0.07)0.02−0.67−0.05−0.670.01−0.05−0.19
      Slope MY0.24 (0.04)0.04−0.56−0.05−0.83−0.17−0.85
      Level FP0.82 (0.11)0.010.750.04−0.100.14
      Slope FP0.23 (0.01)0.060.88−0.100.36
      Level PP0.74 (0.11)0.000.110.10
      Slope PP0.28 (0.01)−0.070.69
      Level SCS0.32 (0.07)0.04
      Slope SCS0.16 (0.05)

      GWAS

      Eight SNP were significantly associated with the considered traits (FDR <0.10; Bonferroni P-value <0.08; Table 4). Two SNP were associated with the slope of MY. Four, including the 3 top significant SNP, were associated with the level of FP and 2 for PP (1 to level and 1 to slope). No significant association was found for principal components extracted from DTD for SCS.
      Table 4Markers significantly associated [false discovery rate (FDR) <0.10] with scores of principal components extracted from daughter trait deviations for milk yield (MY) and fat (FP) and protein (PP) contents
      SNPBTAPosition (bp)P-GC_P Bonf
      Level of significance of the test adjusted with Bonferroni correction.
      FDRTrait
      ARS-BFGL-NGS-4939141,801,1160.0000010.000001Level FP
      ARS-BFGL-NGS-57820141,651,3110.0000210.000011Level FP
      ARS-BFGL-NGS-107379142,054,4570.0003350.000112Level FP
      ARS-BFGL-NGS-296782622,383,6450.0009730.000973Slope MY
      Hapmap30383-BTC-005848141,489,4960.0492440.012311Level FP
      ARS-BFGL-NGS-192755114,818,2060.0538300.053830Slope PP
      Hapmap32110-BTA-153952635,555,2470.0779470.038974Slope MY
      Hapmap32435-BTC-0121881456,075,4350.0811550.081155Level PP
      1 Level of significance of the test adjusted with Bonferroni correction.

      MY

      The 2 markers significantly associated with DTD_MY were both related to slope (Table 4 and Figure 3a), the PC that expresses the shape of individual curves for increasing levels of THI. The first SNP was located on BTA26 at approximately 22.3 Mb. A possible candidate gene located within the interval defined by this marker is the β-transducin repeat containing E3 ubiquitin protein ligase (BTRC; Table 5), reported to be associated with milk production (
      • Raven L.A.
      • Cocks B.G.
      • Kemper K.E.
      • Chamberlain A.J.
      • Goddard M.E.
      • Hayes B.J.
      Targeted imputation of sequence variants and gene expression profiling identifies twelve candidate genes associated with lactation volume, composition and calving interval in dairy cattle.
      ) and leg morphology (
      • van den Berg I.
      • Fritz S.
      • Rodriguez S.
      • Rocha D.
      • Boussaha M.
      • Lund M.S.
      • Boichard D.
      Concordance analysis for QTL detection in dairy cattle: A case study of leg morphology.
      ) in cattle, and with growth rate in chicken (
      • Zhang G.X.
      • Fan Q.C.
      • Zhang T.
      • Wang J.Y.
      • Wang W.H.
      • Xue Q.
      • Wang Y.J.
      Genome-wide association study of growth traits in the Jinghai Yellow chicken.
      ). The 0.5-Mb window includes also genes involved in folliculogenesis in cattle, such as fibroblast growth factor 8 (FGF8), found in a selection sweep study in dairy cattle (
      • Kemper K.E.
      • Saxton S.J.
      • Bolormaa S.
      • Hayes B.J.
      • Goddard M.E.
      Document selection for complex traits leaves little or no classic signatures of selection.
      ), meningioma expressed antigen 5 (hyaluronidase; MGEA5), and taurus Kv channel interacting protein 2 (KCNIP2;
      • Hatzirodos N.
      • Irving-Rodgers H.F.
      • Hummitzsch K.
      • Harland M.L.
      • Morris S.E.
      • Rodgers R.J.
      Transcriptome profiling of granulosa cells of bovine ovarian follicles during growth from small to large antral sizes.
      ). Of interest is also Hermansky-Pudlak syndrome 6 (HPS6), a gene related to pigmentation in humans (
      • Sturm R.A.
      • Duffy D.L.
      Human pigmentation genes under environmental selection.
      ). A marker located on BTA26 and significantly associated with sweating rate was reported by
      • Dikmen S.
      • Cole J.B.
      • Null D.J.
      • Hansen P.J.
      Genome-wide association mapping for identification of quantitative trait loci for rectal temperature during heat stress in Holstein cattle.
      for Holsteins, although the map position is about 2.0 Mb from the SNP identified in our study.
      Figure thumbnail gr3
      Figure 3Genome-wide association study of the scores of the second principal component (slope) for milk yield (MY; a), fat percentage (FP; b), and protein percentage (PP; c). On the vertical axis is the negative logarithm of the P-value corrected for the stratification of the population. On the horizontal axis are the SNP ordered by their position and by chromosome. The line corresponds to a false discovery rate of 0.10. GC = genomic control.
      Table 5Putative candidate genes located in the 0.5-Mb interval surrounding the significant SNP
      BTAPosition (Mb)Gene name and symbolTrait
      Level and slope of principle components for milk yield (MY), fat percentage (FP), and protein percentage (PP).
      FunctionAssociations with phenotypic traits
      2622.09–22.17β-Transducin repeat containing E3 ubiquitin protein ligase, BTRCSlope MYIntracellular protein degradationMilk production, leg morphology, growth rate
      2622.37–22.38Fibroblast growth factor 8, FGF8Slope MYMitogenic and cell survivalFolliculogenesis
      2622.39–22.42Meningioma expressed antigen 5 (hyaluronidase), MGEA5Slope MYProtein glycosilationFolliculogenesis
      2622.42–22.43TaurusKv channel interacting protein 2, KCNIP2Slope MYIon transportationFolliculogenesis
      2622.62–22.63Hermansky-Pudlak syndrome 6, HPS6Slope MYMelanosome developmentPigmentation
      635.10–35.95Coiled-coil serine rich protein 1, CCSER1Slope MYCell divisionBeef traits
      141.79–1.81Diacylglycerol O-acyltransferase 1, DGAT1Level FPFat metabolismMilk fat content
      141.81–1.83Heat shock transcription factor 1, HSF1Level FPResponse to cell stressHeat stress tolerance, Mycobacterium bovis susceptibility
      141.56–1.60Rho GTPase activating protein 39, ARHGAP39Level FPNervous system developmentMilk fat content
      141.50–1.51Ribosomal protein L8, RPL8Level FPHomeostasisHeat stress (fish)
      142.23–2.24Mitogen-activated protein kinase 15, MAPK15Level FPCell proliferationSCS
      141.49–1.50Zinc finger protein 34, ZNF34Level FPTranscription regulationMilk fat
      511.45–11.46Malonyl-CoA-acyl carrier protein transacylase, MCATSlope PPFatty acid metabolismHeat stress (chicken)
      511.49–11.50Sorting and assembly machinery component, SAMM50Slope PPMitochondrial proteinTriglyceride levels (humans)
      511.45–11.46Translocator protein, TSPOSlope PPCholesterol transportationFolliculogenesis
      1 Level and slope of principle components for milk yield (MY), fat percentage (FP), and protein percentage (PP).
      The second marker associated with slope for MY was located on BTA6 at approximately 35.5 Mb (Table 4). In this region maps the coiled-coil serine rich protein 1 (CCSER1; Table 5), a gene involved in the mechanism of cell division found to be associated with birth weight in Brangus cattle (
      • Saatchi M.
      • Schnabel R.D.
      • Taylor J.F.
      • Garrick D.J.
      Large-effect pleiotropic or closely linked QTL segregate within and across ten US cattle breeds.
      ) and with Na concentration in muscle of Nelore cattle (
      • Tizioto P.C.
      • Taylor J.F.
      • Decker J.E.
      • Gromboni C.F.
      • Mudadu M.A.
      • Schnabel R.D.
      • Coutinho L.L.
      • Mourão G.B.
      • Oliveira P.S.N.
      • Souza M.M.
      • Reecy J.M.
      • Nassu R.T.
      • Bressani F.A.
      • Tholon P.
      • Sonstegard T.S.
      • Alencar M.M.
      • Tullio R.R.
      • Nogueira A.R.A.
      • Regitano L.C.A.
      Detection of quantitative trait loci for mineral content of Nelore longissimus dorsi muscle.
      ).
      • Dikmen S.
      • Cole J.B.
      • Null D.J.
      • Hansen P.J.
      Genome-wide association mapping for identification of quantitative trait loci for rectal temperature during heat stress in Holstein cattle.
      ,
      • Dikmen S.
      • Wang X.Z.
      • Ortega M.S.
      • Cole J.B.
      • Null D.J.
      • Hansen P.J.
      Single nucleotide polymorphisms associated with thermoregulation in lactating dairy cows exposed to heat stress.
      found markers significantly associated with rectal temperature and respiration rate on BTA6, but at a position about 10 Mb from the marker flagged in the present study.

      FP

      The 4 significant SNP detected for level of FP were located on BTA14 in an interval of approximately 0.6 Mb (Table 4). The top significant SNP (ARS-BFGL-NGS-4939) was reported to be significantly associated with milk fat in Italian, US, and German Holsteins (
      • Cole J.B.
      • Wiggans G.R.
      • Ma L.
      • Sonstegard T.S.
      • Lawlor Jr., T.J.
      • Crooker B.A.
      • Van Tassell C.P.
      • Yang J.
      • Wang S.
      • Matukumalli L.K.
      • Da Y.
      Genome-wide association analysis of thirty one production, health, reproduction and body conformation traits in contemporary U.S. Holstein cows.
      ;
      • Wang X.
      • Wurmser C.
      • Pausch H.
      • Jung S.
      • Reinhardt F.
      • Tetens J.
      • Thaller G.
      • Fries R.
      Identification and dissection of four major QTL affecting milk fat content in the German Holstein-Friesian population.
      ;
      • Capomaccio S.
      • Milanesi M.
      • Bomba L.
      • Cappelli K.
      • Nicolazzi E.L.
      • Williams J.L.
      • Ajmone-Marsan P.
      • Stefanon B.
      Searching new signals for production traits through gene-based association analysis in three Italian cattle breeds.
      ), in Italian Simmental cattle (
      • Macciotta N.P.P.
      • Gaspa G.
      • Bomba L.
      • Vicario D.
      • Dimauro C.
      • Cellesi M.
      • Ajmone-Marsan P.
      Genome-wide association analysis in Italian Simmental cows for lactation curve traits using a low-density (7K) SNP panel.
      ), and in a multibreed population (
      • Raven L.A.
      • Cocks B.G.
      • Hayes B.J.
      Multibreed genome wide association can improve precision of mapping causative variants underlying milk production in dairy cattle.
      ). The 0.5-Mb window surrounding this top marker flagged a zone characterized by a high density of annotated genes. In particular, in that region maps one of the most important genes affecting milk fat content and yield in cattle (Table 5), the diacylglycerol o-acyltransferase 1 (DGAT1) (
      • Grisart B.
      • Coppieters W.
      • Farnir F.
      • Karim L.
      • Ford C.
      • Berzi P.
      • Cambisano N.
      • Mni M.
      • Reid S.
      • Simon P.
      • Spelman R.
      • Georges M.
      • Snell R.
      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.
      ). However, on BTA14 at approximately 1.81 to 1.83 Mb maps the heat shock transcription factor 1 (HSF1), a protein that is involved in the mechanism of response to heat stress (
      • Guettouche T.
      • Boellmann F.
      • Lane W.S.
      • Voellmy R.
      Analysis of phosphorylation of human heat shock factor 1 in cells experiencing a stress.
      ). Among the genes that map in the interval defined by the second marker (Table 4), of interest is the Rho GTPase activating protein 39 (ARHGAP39), involved in the development of the central nervous system (
      • Ma S.
      • Nowak F.V.
      The RhoGAP domain-containing protein, Porf-2, inhibits proliferation and enhances apoptosis in neural stem cells.
      ) (Table 5). This gene has been reported to be associated with milk fat composition in Danish (
      • Buitenhuis B.
      • Janss L.L.G.
      • Poulsen N.A.
      • Larsen L.B.
      • Larsen M.K.
      • Sørensen P.
      Genome-wide association and biological pathway analysis for milk-fat composition in DanishHolstein and Danish Jersey cattle.
      ) and North American (
      • Nayeri S.
      • Sargolzaei M.
      • Abo-Ismail M.K.
      • May N.
      • Miller S.P.
      • Schenkel F.
      • Moore S.S.
      • Stothard P.
      Genome-wide association for milk production and female fertility traits in Canadian dairy Holstein cattle.
      ) Holstein, and to SCS (
      • 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.
      ) in Chinese Holstein. Another interesting gene located in this region is the ribosomal protein L8 (RPL8), involved in the cellular mechanisms of homeostasis (
      • Katz M.J.
      • Gandara L.
      • De Lella Ezcurra A.L.
      • Wappner P.
      Hydroxylation and translational adaptation to stress: Some answers lie beyond the STOP codon.
      ), associated with response to acute heat stress in the fish Lates calcarifer (
      • Newton J.R.
      • De Santis C.
      • Jerry D.R.
      The gene expression response of the catadromous perciform barramundi Latescalcarifer to an acute heat stress.
      ).
      The third marker located on BTA14 displayed the mitogen-activated protein kinase 15 (MAPK15) (Table 5), suggested as a candidate gene for SCS in a study on Chinese Holstein (
      • 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.
      ). The last significant SNP for PC1 of DTD_FP was also located on BTA14, very close to the second significant SNP. Close to this marker maps the zinc finger protein 34 (ZNF34), found to be associated with milk fat percentage in Chinese Holstein (
      • Jiang L.
      • Liu X.
      • Yang J.
      • Wang H.
      • Jiang J.
      • Liu L.
      • He S.
      • Ding X.
      • Liu J.
      • Zhang Q.
      Targeted resequencing of GWAS loci reveals novel genetic variants for milk production traits.
      ) and Italian Simmental (
      • Macciotta N.P.P.
      • Gaspa G.
      • Bomba L.
      • Vicario D.
      • Dimauro C.
      • Cellesi M.
      • Ajmone-Marsan P.
      Genome-wide association analysis in Italian Simmental cows for lactation curve traits using a low-density (7K) SNP panel.
      ). No significant SNP were detected for the slope component for fat percentage (Figure 3b)

      PP

      A SNP significantly associated with slope of PP was found on BTA5, at approximately 114.8 Mb (Table 4 and Figure 3c). In the 0.5-Mb interval flanking this marker is an interesting gene, malonyl-CoA-acyl carrier protein transacylase (MCAT; Table 5), expressed in the mitochondrion and involved in fatty acid metabolism. This gene has been found to be differentially expressed in chicken embryos exposed to heat challenge (
      • Loyau T.
      • Hennequet-Antier C.
      • Coustham V.
      • Berri C.
      • Leduc M.
      • Crochet S.
      • Sannier M.
      • Duclos M.J.
      • Mignon-Grasteau S.
      • Tesseraud S.
      • Brionne A.
      • Métayer-Coustard S.
      • Moroldo M.
      • Lecardonnel J.
      • Martin P.
      • Lagarrigue S.
      • Yahav S.
      • Collin A.
      Thermal manipulation of the chicken embryo triggers differential gene expression in response to a later heat challenge.
      ). Other genes of interest located in this interval are the sorting and assembly machinery component (SAMM50), a mitochondrial protein found associated with serum triglyceride levels in humans (
      • Kitamoto T.
      • Kitamoto A.
      • Yoneda M.
      • Hyogo H.
      • Ochi H.
      • Nakamura T.
      • Teranishi H.
      • Mizusawa S.
      • Ueno T.
      • Chayama K.
      • Nakajima A.
      • Nakao K.
      • Sekine A.
      • Hotta K.
      Genome-wide scan revealed that polymorphisms in the PNPLA3, SAMM50, and PARVB genes are associated with development and progression of nonalcoholic fatty liver disease in Japan.
      ), and the translocator protein (TSPO), a gene upregulated in atretic bovine follicles (
      • Hatzirodos N.
      • Irving-Rodgers H.F.
      • Hummitzsch K.
      • Harland M.L.
      • Morris S.E.
      • Rodgers R.J.
      Transcriptome profiling of granulosa cells of bovine ovarian follicles during growth from small to large antral sizes.
      ).
      The other significant marker associated with PP was related to level. It was located on BTA14, in a region where no annotated genes have been retrieved in the UCSC genome database.

      DISCUSSION

      Traits able to describe efficiently the response of animals to heat stress are rather problematic to be routinely measured in dairy cattle populations involved in breeding programs.
      • Dikmen S.
      • Cole J.B.
      • Null D.J.
      • Hansen P.J.
      Heritability of rectal temperature and genetic correlations with production and reproduction traits in dairy cattle.
      estimated that 13 to 17% of the variation in rectal temperature in cows during heat stress is due to genetic differences. However, it will not be practical to select cows for heat tolerance based on rectal temperature directly because this trait is not recorded on dairy farms. On the contrary, patterns of milk production traits can be conveniently modeled to estimate tolerance to heat stress (
      • Carabaño M.J.
      • Bachaga K.
      • Ramon M.
      • Diaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: selection tools.
      ;
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effect of heat stress in Italian Holstein dairy cattle.
      ;
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      ). In particular, DTD have been proposed as a proxy of individual response of bulls across increasing levels of heat load (
      • Hayes B.J.
      • Bowman P.J.
      • Chamberlain A.J.
      • Savin K.
      • van Tassell C.P.
      • Sostengard T.S.
      • Goddard M.E.
      A validated genome wide association study to breed cattle adapted to an environment altered by climate changes.
      ). The reliability of heat tolerance proxies based on production traits was underlined in a recent study on dairy cattle subjected to heat challenge under a controlled environment, namely climate chambers (
      • Garner J.B.
      • Douglas M.L.
      • Williams S.R.O.
      • Wales W.J.
      • Marett L.C.
      • Nguyen T.T.T.
      • Reich C.M.
      • Hayes B.J.
      Genomic selection improves heat tolerance in dairy cattle.
      ). Heat-tolerant cows, ranked according to a genomic breeding value of heat tolerance based on milk yield, showed lower values of physiological indicators (core body temperature, rectal and intravaginal temperature) than heat-sensitive cows.
      In the present work, DTD calculated for 4 different milk production traits were analyzed with PCA. This approach was able to extract 2 new variables, which explained approximately 50% of the original variance. Interestingly, they were related to the level of the DTD curve and to its slope, respectively. Our interpretation of PC meaning is in agreement with the outcome of eigendecomposition of coefficient matrix of random regression models used to estimate heat tolerance (
      • Carabaño M.J.
      • Bachaga K.
      • Ramon M.
      • Diaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: selection tools.
      ) and with previous reports on PCA carried out on milk production traits in cattle (
      • Macciotta N.P.P.
      • Vicario D.
      • Cappio-Borlino A.
      Use of multivariate analysis to extract latent variables related to level of production and lactation persistency in dairy cattle.
      ,
      • Macciotta N.P.P.
      • Gaspa G.
      • Bomba L.
      • Vicario D.
      • Dimauro C.
      • Cellesi M.
      • Ajmone-Marsan P.
      Genome-wide association analysis in Italian Simmental cows for lactation curve traits using a low-density (7K) SNP panel.
      ).
      Most of the genetic models used to study the heat stress effect on dairy traits use the reaction norm model approach (
      • Kolmodin R.
      • Bijma P.
      Response to mass selection when the genotype by environment interaction is modeled as a linear reaction norm.
      ;
      • Shariati M.M.
      • Su G.
      • Madsen P.
      • Sorensen D.
      Analysis of milk production traits in early lactation using a reaction norm model with unknown covariates.
      ). In particular, for each animal the THI effect is fitted as a general intercept plus a slope (
      • Hayes B.J.
      • Bowman P.J.
      • Chamberlain A.J.
      • Savin K.
      • van Tassell C.P.
      • Sostengard T.S.
      • Goddard M.E.
      A validated genome wide association study to breed cattle adapted to an environment altered by climate changes.
      ;
      • Sánchez J.P.
      • Misztal I.
      • Aguilar I.
      • Zumbach B.
      • Rekaya R.
      Genetic determination of the onset of heat stress on daily milk production in the US Holstein cattle.
      ;
      • Carabaño M.J.
      • Bachaga K.
      • Ramon M.
      • Diaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: selection tools.
      ). The results of the model-free PCA approach used in the present study basically confirm theoretical assumptions of the reaction norm model. Also, the predominance of level over slope in terms of variance explained is in agreement with previous reports. In particular, slope eigenvalue (about 10%) is not far from values reported for the second eigenfunction for fat and protein yield obtained from the decomposition of (co)variance matrix of random regression models by
      • Carabaño M.J.
      • Bachaga K.
      • Ramon M.
      • Diaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: selection tools.
      . These suggest that the behavior of milk production traits across increasing THI levels can be partitioned into 2 main components, one basically related to the overall production genetic potential of the animal and the other to the individual specific response.
      The slope component could be proposed as a measure of individual tolerance to heat stress; however, the interpretation of its values deserves further discussion. From Figures 3b it can be noted that animals with negative slope scores exhibited, on average, lower DTD in the first part of the THI scale (i.e., in the comfort zone) and higher in the second part. Thus, selection in favor of this kind of response pattern would result in less productive animals under comfort conditions (i.e., the most frequent at least in temperate areas). However, it should be remembered that graphs depicted in Figures 2 represent deviations from the average response pattern (
      • Carabaño M.J.
      • Logar B.
      • Bormann J.
      • Minet J.
      • Vanrobays M.L.
      • Diaz C.
      • Tychon B.
      • Gengler N.
      • Hammami H.
      Modeling heat stress under different environmental conditions.
      ). Thus, animals exhibiting a curve as the one reported in Figure 2b are not expected to respond with a decrease of production to lower THI levels, but to have a quite constant average level of production in this point of the curve. Moreover, provided the orthogonality between the 2 traits, the crossed distribution of individual patterns among the classes of the 2 PC (Table 6) shows values in all the cells. Although intermediate classes (i.e., those having values between −1 and 1 for both the components) are the most abundant, 10 bulls belong to the class having the most positive and negative values for level and slope, respectively. The DTD_MY patterns of some of these bulls are reported in Figure 4; they exhibit the typical ascending pattern of this slope class, with the points mostly having values higher than zero or showing great variability.
      Table 6Absolute frequencies of individual patterns of daughter trait deviations for milk yield across different classes of level and slope principal component scores
      LevelSlope
      ≤−2−2 to −1−1 to 00 to 11 to 2≥2
      ≤−210247281216
      −2 to −1
      Lower limit of the class is included.
      4158280254
      −1 to 0532138137314
      0 to 1427128137283
      1 to 27219084168
      ≥2103270622715
      1 Lower limit of the class is included.
      Figure thumbnail gr4
      Figure 4Individual patterns of daughter trait deviations (DTD) for milk yield, which have scores >2 and <−2 for the level and slope components, respectively. THI = temperature-humidity index.
      Another point of agreement between the results of the present study and those of other studies on heat tolerance based on the broken line model can be found in the structure of PC eigenvectors. It is worth noting that the inversion of eigenvector coefficient sign of slope (Table 1) occurs in the THI class 68 to 70 for DTD of MY, PP, and SCS, and in the class 65 to 67 for FP. These values basically agreed with estimates of TH0 threshold reported by some authors.
      • Sánchez J.P.
      • Misztal I.
      • Aguilar I.
      • Zumbach B.
      • Rekaya R.
      Genetic determination of the onset of heat stress on daily milk production in the US Holstein cattle.
      estimated a THI threshold of 72 for milk yield in US Holsteins using hierarchical models.
      • Carabaño M.J.
      • Logar B.
      • Bormann J.
      • Minet J.
      • Vanrobays M.L.
      • Diaz C.
      • Tychon B.
      • Gengler N.
      • Hammami H.
      Modeling heat stress under different environmental conditions.
      estimated TH0 thresholds of 72 to 73 for milk yield in Holstein populations of different European countries. The lower threshold of fat percentage compared with the other traits reported in the present study also agrees with previous reports. In our previous study, THI thresholds of 76, 73, and 74 were estimated in first-, second-, and third-parity Holstein cows, respectively, for MY (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effect of heat stress in Italian Holstein dairy cattle.
      ).
      In view of a possible implementation of PCA-derived measures of tolerance to heat stress in breeding programs, the 2-stage approach suggested in the current study has the limitation that only bulls having complete records (i.e., the 11 points of the DTD curve) could be considered. Such a requirement strongly reduced the number of animals for which PC could be computed (for example, from 4,184 to 1,540 for MY in the present study). To overcome this problem, the correlation matrix between different points of the DTD curve could be reconstructed by using a random regression model (RRM) in which records at different THI classes are treated as repeated measures for each sire. Thus, DTD_MY were fitted with a mixed model having the same structure of [1] but with the sire effect fitted as random. The covariance within animal was accounted for by an 11 × 11 unstructured matrix of between sire effects for each THI level. To avoid convergence problems, bulls having 7 or more records were considered. A total of 35,992 records belonging to 3,697 bulls were used. From Table 7, it can be seen that observed Pearson correlations among DTD_MY and the correlation matrix estimated by the RRM are very similar. Finally, PCA carried out on the RRM estimated correlation matrix yielded basically the same results (Table 8) of the 2-step approach, both in terms of variance explained by the first 2 PC and of their eigenvector structure (i.e., the first can be considered a level PC, the second a slope). Such a concordance of results opens interesting perspectives on the possible use of this heat tolerance indices on large scale.
      Table 7Pearson correlations among daughter trait deviations for milk yield (above the diagonal) and estimated unstructured correlation matrix estimated with random regression model (under the diagonal)
      THI = temperature-humidity index.
      ItemTHI1THI2THI3THI4THI5THI6THI7THI8THI9THI110THI11
      THI10.410.370.350.360.370.330.290.310.300.23
      THI20.430.420.400.400.400.370.360.350.330.25
      THI30.380.440.430.460.440.430.410.430.360.26
      THI40.350.410.430.410.440.420.420.400.340.33
      THI50.350.400.440.410.440.440.410.430.370.30
      THI60.340.390.420.420.420.440.440.420.430.35
      THI70.320.370.410.410.430.400.430.450.430.34
      THI80.280.350.380.400.370.410.390.420.410.43
      THI90.270.300.350.360.390.390.440.420.450.37
      THI100.230.260.330.270.310.370.430.420.470.41
      THI110.180.220.220.260.260.340.340.410.410.43
      1 THI = temperature-humidity index.
      Table 8Eigenvectors and eigenvalues of the first 2 principal components (PC) extracted from the correlation matrix estimated by fitting a random regression model to daughter trait deviations for milk yield (DTD_MY)
      Temperature humidity index intervalPC1PC2
      50–520.26−0.40
      53–550.30−0.36
      56–580.31−0.27
      59–610.31−0.21
      62–640.31−0.18
      65–670.32−0.03
      68–700.320.07
      71–730.310.19
      74–760.310.30
      77–790.290.42
      >790.270.51
      Eigenvalue4.661.15
      Eigenvalue %4210
      Previous studies on the genetic basis of tolerance to heat stress indicated that this trait has a moderate genetic variability and that it could be split into an intercept and a slope component that are genetically related. The PC-based measure of heat tolerance proposed in the present work does show genetic variability. In particular, the pseudo-heritability of the slope component is similar to heritability values reported in other studies for the slope parameter of the reaction norm model. Thus, the 2 PC extracted from DTD of 4 milk traits could be considered as possible breeding criteria when selecting for improved heat tolerance in cattle. However, compared with previous measurements of tolerance to heat stress based on milk production data, a distinguishing feature of the PC is their phenotypic and genetic orthogonality. The absence of any genetic relationship between level and slope (Table 2) suggests that an independent selection of the 2 main aspects of DTD patterns (i.e., level of production and heat tolerance) should be feasible. A simultaneous selection for improving both heat tolerance and dairy traits could be achieved also by implementing a selection index in which suitable economic weight have to be determined (
      • Nguyen T.T.T.
      • Bowman P.J.
      • Haile-Mariam M.
      • Pryce J.E.
      • Hayes B.J.
      Genomic selection for tolerance to heat stress in Australian dairy cattle.
      ). However, provided the unfavorable genetic correlation between heat tolerance and production, a smaller selection response is expected for each single trait in comparison with the use of the 2 uncorrelated level and slope components.
      The amount of variance accounted for by the slope parameter in reaction norm models is used as an indicator of the genotype × environment interaction (
      • Kolmodin R.
      • Bijma P.
      Response to mass selection when the genotype by environment interaction is modeled as a linear reaction norm.
      ;
      • Shariati M.M.
      • Su G.
      • Madsen P.
      • Sorensen D.
      Analysis of milk production traits in early lactation using a reaction norm model with unknown covariates.
      ). In the present study, slope and level additive variance ratios were 0.10 for MY, 0.05 for FP, 0.06 for PP, and 0.13 for SCS. These values confirm results of
      • Santana M.L.
      • Bignardi A.B.
      • Pereira R.J.
      • Stefani G.
      • El Faro L.
      Genetics of heat tolerance for milk yield and quality in Holsteins.
      , who concluded that a genotype × environment interaction due to heat stress is more relevant for milk yield and SCS than for fat and protein percentage.
      Association analysis highlighted a limited number of significant markers. This was not an unexpected outcome due to the severe correction needed to account for multiple testing when high-throughput platforms are used, to the complex biology underlying the physiological response to heat stress, and to the limited size of the sample of animals genotyped.
      Half of the significant SNP were associated with the level variable for FP and located on BTA14. This result is quite common in genomic studies on Holstein cattle, mainly due to the genetic architecture of the trait, which is largely influenced by a single segregating gene, DGAT1 (
      • Grisart B.
      • Coppieters W.
      • Farnir F.
      • Karim L.
      • Ford C.
      • Berzi P.
      • Cambisano N.
      • Mni M.
      • Reid S.
      • Simon P.
      • Spelman R.
      • Georges M.
      • Snell R.
      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.
      ). The relevant influence of this gene on FP was previously observed in the Italian Holstein population (
      • Fontanesi L.
      • Calò D.G.
      • Galimberti G.
      • Negrini R.
      • Marino R.
      • Nardone A.
      • Ajmone-Marsan P.
      • Russo V.
      A candidate gene association study for nine production traits in Italian Holstein sires.
      ).
      However, it is worth mentioning that a gene encoding for a heat shock transcription factor (HSF1) maps on BTA14 very close to the top significant marker detected in the present study. In humans, HSF1 mediates the expression of the heat shock and stress proteins in response to physical and chemical stresses (Guettoche et al., 2005). This gene has been recognized to have a central role in coordinating thermal tolerance in cattle (
      • Collier R.J.
      • Collier J.L.
      • Rhoads R.P.
      • Baumgard L.H.
      Invited review: Genes involved in the bovine heat stress response.
      ). Differences in the expression of HSF1 in the liver were found between cows calving in spring and summer, respectively (
      • Shahzad K.
      • Akbar H.
      • Vailati-Riboni M.
      • Basiricò L.
      • Morera P.
      • Rodriguez-Zas S.L.
      • Nardone A.
      • Bernabucci U.
      • Loor J.J.
      The effect of calving in the summer on the hepatic transcriptome of Holstein cows during the peripartal period.
      ). Moreover, associations between this gene and thermotolerance have been detected in Chinese Holstein (
      • Li Q.L.
      • Ju Z.H.
      • Huang J.M.
      • Li J.B.
      • Li R.L.
      • Hou M.H.
      • Wang C.F.
      • Zhong J.F.
      Two novel SNPs in HSF1 gene are associated with thermal tolerance traits in Chinese Holstein cattle.
      ). An overexpression of HSF1 has been found in buffalo during summer under tropical environment (
      • Kumar A.
      • Ashrafa S.
      • Sridhar Gouda T.
      • Grewalc A.
      • Singha S.V.
      • Yadavb B.R.
      • Upadhyaya R.C.
      Expression profiling of major heat shock protein genes during different seasons in cattle (Bos indicus) and buffalo (Bubalus bubalis) under tropical climatic condition.
      ) and is reported to be associated with genetic susceptibility to Mycobacterium bovis infection in dairy cattle (
      • Richardson I.W.
      • Berry D.P.
      • Wiencko H.L.
      • Higgins I.M.
      • More S.J.
      • McClure J.
      • Lynn D.J.
      • Bradley D.G.
      A genome-wide association study for genetic susceptibility to Mycobacterium bovis infection in dairy cattle identifies a susceptibility QTL on chromosome 23.
      ). Differential expression of heat shock proteins has been related to in vitro fertilization rate and blastocyst rate of bovine embryos (
      • Zhang B.
      • Peñagaricano F.
      • Driver A.
      • Chen H.
      • Khatib H.
      Differential expression of heat shock protein genes and their splice variants in bovine preimplantation embryos.
      ). The HSP70.1 polymorphism has been associated with cellular thermos tolerance in Holstein lactating cows (
      • Basiricò L.
      • Morera P.
      • Primi V.
      • Lacetera N.
      • Nardone A.
      • Bernabucci U.
      Cellular thermotolerance is associated with heat shock protein 70.1 genetic polymorphisms in Holstein lactating cows.
      )
      The relationships between tolerance to heat stress and fertility has been confirmed by associations found in the present study. Candidate genes found in the 0.5-Mb intervals defined by the significant markers for the slope component are involved in cellular regulation mechanism, fertility, and weight at birth.
      • Nguyen T.T.T.
      • Bowman P.J.
      • Haile-Mariam M.
      • Pryce J.E.
      • Hayes B.J.
      Genomic selection for tolerance to heat stress in Australian dairy cattle.
      , using the DTD deviation as indicators of heat stress in Australian Holsteins, found favorable correlation between DTD and fertility.
      Among putative candidate genes that have been detected in the present study for the SLOPE trait, one (MCAT) is involved in fatty acid metabolism. This result is in agreement with findings by
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      , who highlighted a relationship between milk fatty acid content and tolerance to heat stress in cattle. It should pointed out that milk fatty acid profile is an index of animal energy balance (
      • 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.
      ) and is strongly related to the diet, which may be affected by climatic conditions.
      • Nardone A.
      • Lacetera N.
      • Bernabucci U.
      • Ronchi B.
      Composition of colostrum from dairy heifers exposed to high air temperatures during late pregnancy and the early postpartum period.
      found greater proportions of long-chain fatty acids in colostrum produced by heifers under heat stress conditions. Those authors demonstrated that the higher proportion of long-chain fatty acids was due to the reduced synthesis of short- and medium-chain fatty acids in the mammary gland cells. Thus, the role of milk FA as potent biomarkers for evaluating individual thermotolerance could be hypothesized.
      Previous GWAS studies carried out on DTD and rectal temperature have highlighted genomic regions associated with heat stress tolerance. These results were not confirmed in the present study, even though significant markers were seen for slope for milk yield on BTA 5, 6, and 26 [i.e., on the same chromosomes where
      • Dikmen S.
      • Cole J.B.
      • Null D.J.
      • Hansen P.J.
      Genome-wide association mapping for identification of quantitative trait loci for rectal temperature during heat stress in Holstein cattle.
      found significant associations for rectal temperature]. The detection of a limited number of significant markers and a poor repeatability of results across studies and populations is a major issue in GWAS studies carried out on livestock species. Sample size, genetic differences among populations (i.e., level of linkage disequilibrium, allelic SNP frequencies) were mentioned as main reasons for such a lack of concordance between experiments. Moreover, the severe corrections of significance levels due to huge number of repeated tests further reduce the number of detected markers. Finally, a more general issue that has been raised in GWAS carried out in humans is that not all the genetic variation of a trait is captured by available markers (i.e., the so-called problem of missing heritability;
      • Gusev A.
      • Bhatia G.
      • Zaitien N.
      • Vilhjalmsson B.J.
      • Diogo D.
      • Stahl E.A.
      • Gregersen P.K.
      • Worthington J.
      • Klareskog L.
      • Rayachaudhri S.
      • Plenge R.M.
      • Pasanius B.
      • Price A.L.
      Quantifying missing heritability at known GWAS loci.
      ). On the other hand, in the case of GWAS for heat tolerance traits, the role of the phenotypes should be carefully considered.
      • Dikmen S.
      • Wang X.Z.
      • Ortega M.S.
      • Cole J.B.
      • Null D.J.
      • Hansen P.J.
      Single nucleotide polymorphisms associated with thermoregulation in lactating dairy cows exposed to heat stress.
      found that 1 out of 4 significant SNP previously detected using a phenotype derived from milk yield was also associated with a physiological indicator of heat stress (i.e., sweating rate). However, SNP validation is problematic also within the same trait.
      • Hayes B.J.
      • Bowman P.J.
      • Chamberlain A.J.
      • Savin K.
      • van Tassell C.P.
      • Sostengard T.S.
      • Goddard M.E.
      A validated genome wide association study to breed cattle adapted to an environment altered by climate changes.
      , on a second independent data set, confirmed only 2 out of 42 SNP significantly associated with the slope component of a reaction norm model fitted to milk yield in a Holstein population. A further issue in genetic studies of tolerance to heat stress is represented by the methodologies used to obtain the environmental variable. The THI calculations, for example, could differ in the kind of variable used (i.e., daily maximum or average temperature, minimum or average relative humidity) and in the time lag with the day of the test. Such heterogeneity of measures could be one of the reasons for the differences between studies in the estimates of the THI upper threshold for the comfort zone. The existence of all these sources of variation that may possibly affect results of studies on the genetic dissection of heat tolerance traits should lead scientists to make efforts to increase the power of their experiments, validate results in independent populations, and in harmonize methodologies for calculating environmental variables. Furthermore, alternatives to traditional methods for measuring physiological indicators of heat tolerance should be found. Recent achievements in precision farming, such as the use of microequipment directly installed on the animal or the use of indirect predictors of physiological traits (as for example midinfrared milk spectra), could provide interesting tools.

      CONCLUSIONS

      The PCA was able to derive 2 new variables able to describe the overall level and the slope of response of milk production traits across increasing levels of THI index. These 2 new phenotypes are uncorrelated and may therefore provide an option for overcoming the problem of the negative correlation between heat tolerance and production level that has been found within the context of milk and weather record analysis. The genetic background of the 2 PC was investigated and some putative candidate genes were proposed. A useful numerical property of the 2 extracted variables is their orthogonality. This feature makes the use of the second PC as a measure of thermotolerance for breeding purposes, which is particularly appealing because many authors have stressed the need for using measures of tolerance to heat stress uncorrelated from the production level. Moreover, this variable could be derived from data that are routinely recorded in breeding programs and, therefore, used on large scale could be proposed.

      ACKNOWLEDGMENTS

      The research was funded by the Italian Ministry of Agriculture (SELMOL grant). Authors thank the Italian Holstein Association (ANAFI, Rome, Italy) for providing milk production and pedigree data.

      REFERENCES

        • Aguilar I.
        • Misztal I.
        • Tsuruta S.
        Genetic components of heat stress for dairy cattle with multiple lactations.
        J. Dairy Sci. 2009; 92 (19841230): 5702-5711
        • Aguilar I.
        • Misztal I.
        • Tsuruta S.
        Short communication: Genetic trends of milk yield under heat stress for US Holsteins.
        J. Dairy Sci. 2010; 93 (20338455): 1754-1758
        • Ali A.K.A.
        • Shook G.E.
        An optimum transformation for somatic cell concentration in milk.
        J. Dairy Sci. 1980; 63: 487-490
        • Amin N.
        • van Duijn C.M.
        • Aulchenko Y.S.
        A genomic background based method for association analysis in related individuals.
        PLoS One. 2007; 2 (18060068): e1274
        • Aulchenko Y.S.
        • Ripke S.
        • Isaacs A.
        • van Duijn C.M.
        GenABEL: An R package for genome-wide association analysis.
        Bioinformatics. 2007; 23: 1294-1296
        • Basiricò L.
        • Morera P.
        • Primi V.
        • Lacetera N.
        • Nardone A.
        • Bernabucci U.
        Cellular thermotolerance is associated with heat shock protein 70.1 genetic polymorphisms in Holstein lactating cows.
        Cell Stress Chaperones. 2011; 16 (21274669): 441-448
        • 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
        • Bernabucci U.
        • Biffani S.
        • Buggiotti L.
        • Vitali A.
        • Lacetera N.
        • Nardone A.
        The effect of heat stress in Italian Holstein dairy cattle.
        J. Dairy Sci. 2014; 97 (24210494): 471-486
        • Bernabucci U.
        • Lacetera L.H.
        • Baumgard R.P.
        • Rhoads B.
        • Ronchi B.
        • Nardone A.
        Metabolic and hormonal acclimation to heat stress in domesticated ruminants.
        Animal. 2010; 4 (22444615): 1167-1183
        • Biffani S.
        • Bernabucci U.
        • Vitali A.
        • Lacetera N.
        • Nardone A.
        Short communication: Effect of heat stress on nonreturn rate of Italian Holstein cows.
        J. Dairy Sci. 2016; 99 (27108174): 5837-5843
        • Bohmanova J.
        • Misztal I.
        • Tsuruta S.
        • Norman H.D.
        • Lawlor T.J.
        Short communication: Genotype by Environment Interaction Due to Heat Stress.
        J. Dairy Sci. 2008; 91 (18218772): 840-846
        • Brügemann K.
        • Gernand E.
        • von Borstel U.U.
        • König S.
        Genetic analyses of protein yield in dairy cows applying random regression models with time-dependent and temperature x humidity-dependent covariates.
        J. Dairy Sci. 2011; 94 (21787948): 4129-4139
        • Buitenhuis B.
        • Janss L.L.G.
        • Poulsen N.A.
        • Larsen L.B.
        • Larsen M.K.
        • Sørensen P.
        Genome-wide association and biological pathway analysis for milk-fat composition in DanishHolstein and Danish Jersey cattle.
        BMC Genomics. 2014; 15 (25511820): 1112
        • Capomaccio S.
        • Milanesi M.
        • Bomba L.
        • Cappelli K.
        • Nicolazzi E.L.
        • Williams J.L.
        • Ajmone-Marsan P.
        • Stefanon B.
        Searching new signals for production traits through gene-based association analysis in three Italian cattle breeds.
        Anim. Genet. 2015; 46 (25997511): 361-370
        • Carabaño M.J.
        • Bachaga K.
        • Ramon M.
        • Diaz C.
        Modeling heat stress effect on Holstein cows under hot and dry conditions: selection tools.
        J. Dairy Sci. 2014; 97 (25262182): 7889-7904
        • Carabaño M.J.
        • Logar B.
        • Bormann J.
        • Minet J.
        • Vanrobays M.L.
        • Diaz C.
        • Tychon B.
        • Gengler N.
        • Hammami H.
        Modeling heat stress under different environmental conditions.
        J. Dairy Sci. 2016; 99 (26923054): 3798-3814
        • Cole J.B.
        • Wiggans G.R.
        • Ma L.
        • Sonstegard T.S.
        • Lawlor Jr., T.J.
        • Crooker B.A.
        • Van Tassell C.P.
        • Yang J.
        • Wang S.
        • Matukumalli L.K.
        • Da Y.
        Genome-wide association analysis of thirty one production, health, reproduction and body conformation traits in contemporary U.S. Holstein cows.
        BMC Genomics. 2011; 12 (21831322): 408
        • Collier R.J.
        • Collier J.L.
        • Rhoads R.P.
        • Baumgard L.H.
        Invited review: Genes involved in the bovine heat stress response.
        J. Dairy Sci. 2008; 91 (18218730): 445-454
        • Dikmen S.
        • Cole J.B.
        • Null D.J.
        • Hansen P.J.
        Heritability of rectal temperature and genetic correlations with production and reproduction traits in dairy cattle.
        J. Dairy Sci. 2012; 95 (22612974): 3401-3405
        • Dikmen S.
        • Cole J.B.
        • Null D.J.
        • Hansen P.J.
        Genome-wide association mapping for identification of quantitative trait loci for rectal temperature during heat stress in Holstein cattle.
        PLoS One. 2013; 8 (23935954): e69202
        • Dikmen S.
        • Wang X.Z.
        • Ortega M.S.
        • Cole J.B.
        • Null D.J.
        • Hansen P.J.
        Single nucleotide polymorphisms associated with thermoregulation in lactating dairy cows exposed to heat stress.
        J. Anim. Breed. Genet. 2015; 132 (26198991): 409-419
        • Fontanesi L.
        • Calò D.G.
        • Galimberti G.
        • Negrini R.
        • Marino R.
        • Nardone A.
        • Ajmone-Marsan P.
        • Russo V.
        A candidate gene association study for nine production traits in Italian Holstein sires.
        Anim. Genet. 2014; 45 (24796806): 576-580
        • Garner J.B.
        • Douglas M.L.
        • Williams S.R.O.
        • Wales W.J.
        • Marett L.C.
        • Nguyen T.T.T.
        • Reich C.M.
        • Hayes B.J.
        Genomic selection improves heat tolerance in dairy cattle.
        Sci. Rep. 2016;
        • Grisart B.
        • Coppieters W.
        • Farnir F.
        • Karim L.
        • Ford C.
        • Berzi P.
        • Cambisano N.
        • Mni M.
        • Reid S.
        • Simon P.
        • Spelman R.
        • Georges M.
        • Snell R.
        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
        • Guettouche T.
        • Boellmann F.
        • Lane W.S.
        • Voellmy R.
        Analysis of phosphorylation of human heat shock factor 1 in cells experiencing a stress.
        BMC Biochem. 2005; 6 (15760475): 4
        • Gusev A.
        • Bhatia G.
        • Zaitien N.
        • Vilhjalmsson B.J.
        • Diogo D.
        • Stahl E.A.
        • Gregersen P.K.
        • Worthington J.
        • Klareskog L.
        • Rayachaudhri S.
        • Plenge R.M.
        • Pasanius B.
        • Price A.L.
        Quantifying missing heritability at known GWAS loci.
        PLoS Genet. 2013; 9 (24385918): e1003993
        • Hammami H.
        • Vandenplas J.
        • Vanrobays M.-L.
        • Rekik B.
        • Bastin C.
        • Gengler N.
        Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
        J. Dairy Sci. 2015; 98 (25958288): 4956-4968
        • Hansen P.J.
        Exploitation of genetic and physiological determinants of embryonic resistance to elevated temperature to improve embryonic survival in dairy cattle during heat stress.
        Theriogenology. 2007; 68 (17482669): S242-S249
        • Hatzirodos N.
        • Irving-Rodgers H.F.
        • Hummitzsch K.
        • Harland M.L.
        • Morris S.E.
        • Rodgers R.J.
        Transcriptome profiling of granulosa cells of bovine ovarian follicles during growth from small to large antral sizes.
        BMC Genomics. 2014; 15 (24422759): 24
        • Hayes B.J.
        • Bowman P.J.
        • Chamberlain A.J.
        • Savin K.
        • van Tassell C.P.
        • Sostengard T.S.
        • Goddard M.E.
        A validated genome wide association study to breed cattle adapted to an environment altered by climate changes.
        PLoS One. 2009; 4 (19688089): e6676
        • Jiang L.
        • Liu X.
        • Yang J.
        • Wang H.
        • Jiang J.
        • Liu L.
        • He S.
        • Ding X.
        • Liu J.
        • Zhang Q.
        Targeted resequencing of GWAS loci reveals novel genetic variants for milk production traits.
        BMC Genomics. 2014; 15 (25510969): 1105
        • Kadzere C.T.
        • Murphy M.R.
        • Silanikove N.
        • Maltz E.
        Heat stress in lactating dairy cows: a review.
        Livest. Prod. Sci. 2002; 77: 59-91
        • Katz M.J.
        • Gandara L.
        • De Lella Ezcurra A.L.
        • Wappner P.
        Hydroxylation and translational adaptation to stress: Some answers lie beyond the STOP codon.
        Cell. Mol. Life Sci. 2016; 73 (26874685): 1881-1893
        • Kelly C.F.
        • Bond T.E.
        Bioclimatic factors and their measurements.
        in: A Guide to Environmental Research in Animals. Natl. Acad. Sci., Washington, DC1971: 7
        • Kemper K.E.
        • Saxton S.J.
        • Bolormaa S.
        • Hayes B.J.
        • Goddard M.E.
        Document selection for complex traits leaves little or no classic signatures of selection.
        BMC Genomics. 2014; 15 (24678841): 246
        • Kitamoto T.
        • Kitamoto A.
        • Yoneda M.
        • Hyogo H.
        • Ochi H.
        • Nakamura T.
        • Teranishi H.
        • Mizusawa S.
        • Ueno T.
        • Chayama K.
        • Nakajima A.
        • Nakao K.
        • Sekine A.
        • Hotta K.
        Genome-wide scan revealed that polymorphisms in the PNPLA3, SAMM50, and PARVB genes are associated with development and progression of nonalcoholic fatty liver disease in Japan.
        Hum. Genet. 2013; 132 (23535911): 783-792
        • Kolmodin R.
        • Bijma P.
        Response to mass selection when the genotype by environment interaction is modeled as a linear reaction norm.
        Genet. Sel. Evol. 2004; 36 (15231233): 435-454
        • Kumar A.
        • Ashrafa S.
        • Sridhar Gouda T.
        • Grewalc A.
        • Singha S.V.
        • Yadavb B.R.
        • Upadhyaya R.C.
        Expression profiling of major heat shock protein genes during different seasons in cattle (Bos indicus) and buffalo (Bubalus bubalis) under tropical climatic condition.
        J. Therm. Biol. 2015; 51 (25965018): 55-64
        • Li Q.L.
        • Ju Z.H.
        • Huang J.M.
        • Li J.B.
        • Li R.L.
        • Hou M.H.
        • Wang C.F.
        • Zhong J.F.
        Two novel SNPs in HSF1 gene are associated with thermal tolerance traits in Chinese Holstein cattle.
        DNA Cell Biol. 2011; 30 (21189066): 247-254
        • Loyau T.
        • Hennequet-Antier C.
        • Coustham V.
        • Berri C.
        • Leduc M.
        • Crochet S.
        • Sannier M.
        • Duclos M.J.
        • Mignon-Grasteau S.
        • Tesseraud S.
        • Brionne A.
        • Métayer-Coustard S.
        • Moroldo M.
        • Lecardonnel J.
        • Martin P.
        • Lagarrigue S.
        • Yahav S.
        • Collin A.
        Thermal manipulation of the chicken embryo triggers differential gene expression in response to a later heat challenge.
        BMC Genomics. 2016; 17 (27142519): 329
        • Ma S.
        • Nowak F.V.
        The RhoGAP domain-containing protein, Porf-2, inhibits proliferation and enhances apoptosis in neural stem cells.
        Mol. Cell. Neurosci. 2011; 46 (21185940): 573-582
        • Macciotta N.P.P.
        • Gaspa G.
        • Bomba L.
        • Vicario D.
        • Dimauro C.
        • Cellesi M.
        • Ajmone-Marsan P.
        Genome-wide association analysis in Italian Simmental cows for lactation curve traits using a low-density (7K) SNP panel.
        J. Dairy Sci. 2015; 98 (26387014): 8175-8185
        • Macciotta N.P.P.
        • Vicario D.
        • Cappio-Borlino A.
        Use of multivariate analysis to extract latent variables related to level of production and lactation persistency in dairy cattle.
        J. Dairy Sci. 2006; 89 (16840636): 3188-3194
        • Misztal I.
        • Tsuruta S.
        • Strabel T.
        • Auvray B.
        • Druet T.
        • Lee D.H.
        BLUPF90 and related programs (BGF90).
        in: Proceedings of the 7th World Congress on Genetics Applied to Livestock Production, Montpellier, France. INRA, France2002: 743-744
        • Nardone A.
        • Lacetera N.
        • Bernabucci U.
        • Ronchi B.
        Composition of colostrum from dairy heifers exposed to high air temperatures during late pregnancy and the early postpartum period.
        J. Dairy Sci. 1997; 80 (9178123): 838-844
        • Nardone A.
        • Ronchi B.
        • Lacetera N.
        • Ranieri M.S.
        • Bernabucci U.
        Effects of climate changes on animal production and sustainability of livestock systems.
        Livest. Sci. 2010; 130: 57-69
        • Nayeri S.
        • Sargolzaei M.
        • Abo-Ismail M.K.
        • May N.
        • Miller S.P.
        • Schenkel F.
        • Moore S.S.
        • Stothard P.
        Genome-wide association for milk production and female fertility traits in Canadian dairy Holstein cattle.
        BMC Genet. 2016; 17 (27287773): 75
        • Newton J.R.
        • De Santis C.
        • Jerry D.R.
        The gene expression response of the catadromous perciform barramundi Latescalcarifer to an acute heat stress.
        J. Fish Biol. 2012; 81 (22747805): 81-93
        • Nguyen T.T.T.
        • Bowman P.J.
        • Haile-Mariam M.
        • Pryce J.E.
        • Hayes B.J.
        Genomic selection for tolerance to heat stress in Australian dairy cattle.
        J. Dairy Sci. 2016; 99 (27037467): 2849-2862
        • Perano K.M.
        • Usack J.G.
        • Angenent L.T.
        • Gebremedhin K.G.
        Production and physiological responses of heat-stressed lactating dairy cattle to conductive cooling.
        J. Dairy Sci. 2015; 98 (26074243): 5252-5261
        • Ravagnolo O.
        • Misztal I.
        Genetic component of heat stress in dairy cattle, parameter estimation.
        J. Dairy Sci. 2000; 83: 2126-2130
        • Raven L.A.
        • Cocks B.G.
        • Hayes B.J.
        Multibreed genome wide association can improve precision of mapping causative variants underlying milk production in dairy cattle.
        BMC Genomics. 2014; 15: 62
        • Raven L.A.
        • Cocks B.G.
        • Kemper K.E.
        • Chamberlain A.J.
        • Goddard M.E.
        • Hayes B.J.
        Targeted imputation of sequence variants and gene expression profiling identifies twelve candidate genes associated with lactation volume, composition and calving interval in dairy cattle.
        Mamm. Genome. 2016; 27: 81-97
        • Richardson I.W.
        • Berry D.P.
        • Wiencko H.L.
        • Higgins I.M.
        • More S.J.
        • McClure J.
        • Lynn D.J.
        • Bradley D.G.
        A genome-wide association study for genetic susceptibility to Mycobacterium bovis infection in dairy cattle identifies a susceptibility QTL on chromosome 23.
        Genet. Sel. Evol. 2016; 48 (26960806): 19
        • Saatchi M.
        • Schnabel R.D.
        • Taylor J.F.
        • Garrick D.J.
        Large-effect pleiotropic or closely linked QTL segregate within and across ten US cattle breeds.
        BMC Genomics. 2014; 15 (24906442): 442-458
        • Sánchez J.P.
        • Misztal I.
        • Aguilar I.
        • Zumbach B.
        • Rekaya R.
        Genetic determination of the onset of heat stress on daily milk production in the US Holstein cattle.
        J. Dairy Sci. 2009; 92 (19620687): 4035-4045
        • Santana M.L.
        • Bignardi A.B.
        • Pereira R.J.
        • Stefani G.
        • El Faro L.
        Genetics of heat tolerance for milk yield and quality in Holsteins.
        Animal. 2017; 11 (27532229): 4-14
        • SAS Institute
        SAS/STAT version 9.2. SAS Inst. Inc., Cary, NC2008
        • Segnalini M.
        • Nardone A.
        • Bernabucci U.
        • Vitali A.
        • Ronchi B.
        • Lacetera N.
        Dynamics of the temperature-humidity index in the Mediterranean basin.
        Int. J. Biometeorol. 2011; 55 (20524014): 253-263
        • Shahzad K.
        • Akbar H.
        • Vailati-Riboni M.
        • Basiricò L.
        • Morera P.
        • Rodriguez-Zas S.L.
        • Nardone A.
        • Bernabucci U.
        • Loor J.J.
        The effect of calving in the summer on the hepatic transcriptome of Holstein cows during the peripartal period.
        J. Dairy Sci. 2015; 98 (26074246): 5401-5413
        • Shariati M.M.
        • Su G.
        • Madsen P.
        • Sorensen D.
        Analysis of milk production traits in early lactation using a reaction norm model with unknown covariates.
        J. Dairy Sci. 2007; 90 (18024770): 5759-5766
        • Sturm R.A.
        • Duffy D.L.
        Human pigmentation genes under environmental selection.
        Genome Biol. 2012; 13 (23110848): 248
        • Tizioto P.C.
        • Taylor J.F.
        • Decker J.E.
        • Gromboni C.F.
        • Mudadu M.A.
        • Schnabel R.D.
        • Coutinho L.L.
        • Mourão G.B.
        • Oliveira P.S.N.
        • Souza M.M.
        • Reecy J.M.
        • Nassu R.T.
        • Bressani F.A.
        • Tholon P.
        • Sonstegard T.S.
        • Alencar M.M.
        • Tullio R.R.
        • Nogueira A.R.A.
        • Regitano L.C.A.
        Detection of quantitative trait loci for mineral content of Nelore longissimus dorsi muscle.
        Genet. Sel. Evol. 2015; 47 (25880074): 15
        • van den Berg I.
        • Fritz S.
        • Rodriguez S.
        • Rocha D.
        • Boussaha M.
        • Lund M.S.
        • Boichard D.
        Concordance analysis for QTL detection in dairy cattle: A case study of leg morphology.
        Genet. Sel. Evol. 2014; 46 (24884971): 31
        • 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
        • Wang X.
        • Wurmser C.
        • Pausch H.
        • Jung S.
        • Reinhardt F.
        • Tetens J.
        • Thaller G.
        • Fries R.
        Identification and dissection of four major QTL affecting milk fat content in the German Holstein-Friesian population.
        PLoS One. 2012; 7 (22792397): e40711
        • Zhang B.
        • Peñagaricano F.
        • Driver A.
        • Chen H.
        • Khatib H.
        Differential expression of heat shock protein genes and their splice variants in bovine preimplantation embryos.
        J. Dairy Sci. 2011; 94 (21787952): 4174-4182
        • Zhang G.X.
        • Fan Q.C.
        • Zhang T.
        • Wang J.Y.
        • Wang W.H.
        • Xue Q.
        • Wang Y.J.
        Genome-wide association study of growth traits in the Jinghai Yellow chicken.
        Genet. Mol. Res. 2015; 14 (26634498): 15331-15338
        • Zimin A.V.
        • Delcher A.L.
        • Florea L.
        • Kelley D.R.
        • Schatz M.C.
        • Puiu D.
        • Hanrahan F.
        • Pertea G.
        • Van Tassell C.P.
        • Sonstegard T.S.
        • Marçais G.
        • Roberts M.
        • Subramanian P.
        • Yorke J.A.
        • Salzberg S.L.
        A whole-genome assembly of the domestic cow.
        Genom. Biol. 2009; 10 (19393038): R42