Advertisement

Genome-wide association studies for heat stress response in Bos taurus × Bos indicus crossbred cattle

Open ArchivePublished:July 03, 2019DOI:https://doi.org/10.3168/jds.2018-15305

      ABSTRACT

      Heat stress is an important issue in the global dairy industry. In tropical areas, an alternative to overcome heat stress is the use of crossbred animals or synthetic breeds, such as the Girolando. In this study, we performed a genome-wide association study (GWAS) and post-GWAS analyses for heat stress in an experimental Gir × Holstein F2 population. Rectal temperature (RT) was measured in heat-stressed F2 animals, and the variation between 2 consecutive RT measurements (ΔRT) was used as the dependent variable. Illumina BovineSNP50v1 BeadChip (Illumina Inc., San Diego, CA) and single-SNP approach were used for GWAS. Post-GWAS analyses were performed by gene ontology terms enrichment and gene-transcription factor (TF) networks, generated from enriched TF. The breed origin of marker alleles in the F2 population was assigned using the breed of origin of alleles (BOA) approach. Heritability and repeatability estimates (± standard error) for ΔRT were 0.13 ± 0.08 and 0.29 ± 0.06, respectively. Association analysis revealed 6 SNP significantly associated with ΔRT. Genes involved with biological processes in response to heat stress effects (LIF, OSM, TXNRD2, and DGCR8) were identified as putative candidate genes. After performing the BOA approach, the 10% of F2 animals with the lowest breeding values for ΔRT were classified as low-ΔRT, and the 10% with the highest breeding values for ΔRT were classified as high-ΔRT. On average, 49.4% of low-ΔRT animals had 2 alleles from the Holstein breed (HH), and 39% had both alleles from the Gir breed (GG). In high-ΔRT animals, the average proportion of animals for HH and GG were 1.4 and 50.2%, respectively. This study allowed the identification of candidate genes for ΔRT in Gir × Holstein crossbred animals. According to the BOA approach, Holstein breed alleles could be associated with better response to heat stress effects, which could be explained by the fact that Holstein animals are more affected by heat stress than Gir animals and thus require a genetic architecture to defend the body from the deleterious effects of heat stress. Future studies can provide further knowledge to uncover the genetic architecture underlying heat stress in crossbred cattle.

      Key words

      INTRODUCTION

      Heat stress is an important factor in the global dairy industry, because it is responsible for large losses in milk production, growth, and reproductive performance, and negatively affects animal health and welfare (
      • Nardone A.
      • Ronchi B.
      • Lacetera N.
      • Ranieri M.S.
      • Bernabucci U.
      Effects of climate changes on animal production and sustainability of livestock systems.
      ). Climate change may induce increases in average temperature and decreases in rainfall, leading to worsened environmental conditions for livestock production in tropical areas. Harsh environmental conditions are common in Brazil, where two-thirds of the territory is located in tropical regions, with prevalent high temperatures due to solar radiation and grazing systems dependent on the rainy season (

      da Cruz, L. V., D. de S. R. Angrimani, B. R. Rui, and M. A. da Silva. 2011. Efeitos do estresse térmico na produção leiteira: revisão de literatura. Rev. Científica Eletrônica Med. Veterinária.

      ).
      Briefly, when animals cannot dissipate sufficient body heat to prevent an increase in their body temperature, heat stress occurs (
      • West J.W.
      Effects of heat-stress on production in dairy cattle.
      ). Modern high-producing dairy cows are more susceptible to the deleterious effects of heat stress (
      • Renaudeau D.
      • Collin A.
      • Yahav S.
      • de Basilio V.
      • Gourdine J.L.
      • Collier R.J.
      Adaptation to hot climate and strategies to alleviate heat stress in livestock production.
      ) due to their very high metabolic rate associated with milk production. Bos indicus animals have a greater ability to regulate their body temperature, because they have undergone natural selection for thousands of years under high temperatures in India. Part of this ability is the result of lower metabolic rates and higher capacity for heat loss through thermoregulatory mechanisms, which helps in heat stress prevention (
      • Hansen P.J.
      Physiological and cellular adaptations of zebu cattle to thermal stress.
      ). In tropical areas, an alternative to overcome heat stress in dairy production systems is the use of crossing or synthetic breeds containing some portion of Bos indicus genetics, such as the Girolando breed, which is derived from crossing Gir and Holstein animals.
      Respiration rate (RR) acts as a powerful thermoregulatory mechanism that helps maintain body temperature through evaporative cooling. This physiological parameter is important in the prediction of heat stress in dairy cattle (
      • McDowell R.E.
      • Hooven N.W.
      • Camoens J.K.
      Effect of climate on performance of Holsteins in first lactation.
      ). Rectal temperature (RT) is an indicator of thermal equilibrium; an increase in this parameter expresses failure or exhaustion of thermoregulatory mechanisms, resulting in heat stress (
      • de Andrade Ferrazza R.
      • Mogollón Garcia H.D.
      • Vallejo Aristizábal V.H.
      • de Souza Nogueira C.
      • Veríssimo C.J.
      • Sartori J.R.
      • Sartori R.
      • Pinhiero Ferreira J.C.
      Thermoregulatory responses of Holstein cows exposed to experimentally induced heat stress.
      ).
      Molecular markers can help in understanding heat stress adaptive parameters of dairy cattle that are relevant to production systems in tropical regions. This approach allows the identification of genomic regions and genes associated with heat stress response, becoming an additional attribute in the process of genetic selection of heat-tolerant animals (
      • 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.
      ).
      • 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 QTL for RT close to genes involved in cell protection against heat stress in Holstein and concluded that some SNP may be useful in genetic selection and identification of genes controlling heat tolerance. In a subsequent study,
      • 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.
      identified specific genetic markers responsible for genetic variation of thermoregulation, showing the importance of genome-wide association studies (GWAS) in finding markers and genes involved in physiological response to heat stress. However, studies about genes underlying RT, via gene-transcription factor (TF) networks and gene ontology (GO) enrichment analyses, in Bos indicus × Bos taurus animals have not been reported yet.
      Recently,
      • Vandenplas J.
      • Calus M.P.L.
      • Sevillano C.A.
      • Windig J.J.
      • Bastiaansen J.W.M.
      Assigning breed origin to alleles in crossbred animals.
      developed an approach that enables the assignment of breed of origin of alleles (BOA) in crossbred animals. Results from this approach can be used to estimate SNP effects depending on the BOA and to improve genomic prediction or GWAS in crossbred animals. According to
      • Sevillano C.A.
      • Napel J.
      • Guimarães S.E.F.
      • Silva F.F.
      • Calus M.P.L.
      , estimated effects and explained variance of SNP strongly associated with crossbred performance are different depending upon which parental breed donated them. Therefore, including breed-specific SNP effect in the genomic evaluation model allowed, in some cases, a better prediction for crossbreed performance (
      • Sevillano C.A.
      • Vandenplas J.
      • Bastiaansen J.W.M.
      • Bergsma R.
      • Calus M.P.L.
      Genomic evaluation for a three-way crossbreeding system considering breed-of-origin of alleles.
      ).
      • Otto P.I.
      • Guimarães S.E.F.
      • Verardo L.L.
      • Azevedo A.L.S.
      • Vandenplas J.
      • Soares A.C.C.
      • Sevillano C.A.
      • Veroneze R.
      • de Fatima M.
      • Pires A.
      • de Freitas C.
      • Prata M.C.A.
      • Furlong J.
      • Verneque R.S.
      • Martins M.F.
      • Panetto J.C.C.
      • Carvalho W.A.
      • Gobo D.O.R.
      • da Silva M.V.G.B.
      • Machado M.A.
      Genome-wide association studies for tick resistance in Bos taurus × Bos indicus crossbred cattle: A deeper look into this intricate mechanism.
      used the BOA approach to evaluate the origin of marker alleles of the candidate genes identified for tick resistance among the same population used in the current study. Those authors observed that most animals classified as resistant for tick infestation showed 2 alleles from the Gir breed, whereas susceptible animals presented 2 alleles from the Holstein breed.
      In this study, we performed a GWAS in a Gir × Holstein F2 experimental population to identify SNP associated with heat stress response and used gene networks to identify the most likely candidate genes. The origins of alleles of selected genes were assigned to help understand the gene mechanisms involved in heat stress response.

      MATERIALS AND METHODS

      Data

      The data were obtained from previous Embrapa research, in which an experimental F2 population was produced by crossing 4 Holstein bulls with 27 Gir cows to generate 150 F1 (1/2 Gir:1/2 Holstein) animals. From this population, 65 F1 females were mated to 4 F1 males to generate 376 F2 animals; the pedigree file contains a total of 476 animals. All F2 animals were raised together in the Embrapa Dairy Cattle experimental station, located in the southeast of Brazil (
      • Machado M.A.
      • Azevedo A.L.S.
      • Teodoro R.L.
      • Pires M.A.
      • Peixoto M.G.C.D.
      • de Freitas C.
      • Prata M.C.A.
      • Furlong J.
      • da Silva M.V.G.B.
      • Guimarães S.E.F.
      • Regitano L.C.A.
      • Coutinho L.L.
      • Gasparin G.
      • Verneque R.S.
      Genome wide scan for quantitative trait loci affecting tick resistance in cattle (Bos taurus × Bos indicus).
      ).
      A total of 341 Gir × Holstein F2 animals were subjected to a heat chamber at 42°C and 60% relative humidity (RH), after an adaptation period of 12 h at 22°C and 50% RH. Each evaluation included 6 animals, and during the whole period that the animals were housed in the heat chamber, water and food were not offered, aiming to reduce their effects on the animals' responses. Parental Gir, Holstein, and F1 animals were not subjected to heat stress.
      The RT and RR of the F2 animals were evaluated in 2 replicates. The first evaluation was performed after the adaptation period, and the second evaluation 6 h after the heat chamber reached 42°C. To measure RT, a digital clinical thermometer was inserted approximately 7.5 cm into each animal's rectum. For RR measurement, flank movements were counted for 30 s. The difference between rectal temperature measurements was calculated, to estimate the variation of rectal temperature (ΔRT) in each animal as follows:
      ΔRT = RTARTB,


      where RTA is the rectal temperature after heat stress and RTB is the rectal temperature before heat stress (after the adaptation period). The variation of respiration rate (ΔRR) was defined similarly to ΔRT. Each animal was subjected twice to the heat chamber to perform measurement of ΔRT and ΔRR.
      Animals were evaluated in contemporary groups aged from 10 and 14 mo. Additional traits that might affect heat stress evaluations were assessed during this period: coat color, coat thickness and length, and density and type of hair (
      • Machado M.A.
      • Azevedo A.L.S.
      • Teodoro R.L.
      • Pires M.A.
      • Peixoto M.G.C.D.
      • de Freitas C.
      • Prata M.C.A.
      • Furlong J.
      • da Silva M.V.G.B.
      • Guimarães S.E.F.
      • Regitano L.C.A.
      • Coutinho L.L.
      • Gasparin G.
      • Verneque R.S.
      Genome wide scan for quantitative trait loci affecting tick resistance in cattle (Bos taurus × Bos indicus).
      ). Ancestry was estimated using Admixture software (
      • Alexander D.H.
      • Novembre J.
      • Lange K.
      Fast model-based estimation of ancestry in unrelated individuals.
      ) with a reduced panel (7,425 SNP), pruned by linkage disequilibrium (LD) between subsequent markers.
      All animals of the experimental population were genotyped with the Illumina BovineSNP50 v1 BeadChip (Illumina, San Diego, CA), including parental Gir, Holstein, F1, and F2 populations. Genotype quality control was implemented using the SNPStats package in R software (
      • Solé X.
      • Guinó E.
      • Valls J.
      • Iniesta R.
      • Moreno V.
      SNPStats: A web tool for the analysis of association studies.
      ). Samples with a call rate <0.90 and heterozygosity of 3.0 standard deviations above or below the observed mean were removed. For quality control of mapped SNP, only autosomal SNP with call rate >0.90 and minor allele frequency (MAF) >0.03 were considered. After quality control, the genotype file contained 40,283 SNP markers.

      Genome-Wide Association Analyses

      A single-SNP GWAS was performed with ASReml software (
      • Gilmour A.R.
      • Gogel B.
      • Cullis B.
      • Thompson R.
      ASReml user guide release 3.0.
      ), applying the following model:
      yijklmnop=μ+CGi+CCj+HTk+SNPl+Sm+β1(CLijklmnopCL¯)+β2(CTijklmnopCT¯)+β3(HDijklmnopHD¯)+β4(VRRijklmnopVRR¯)+β5(BCijklmnopBC¯)+an+peo+eijklmnop,


      where yijklmnop is the variation of rectal temperature (ΔRT); μ is the general mean; CGi is the effect accounting for contemporary groups; CCj is the effect accounting for coat color (j ranges from 1 to 4); HTk is the effect accounting for hair type (k ranges from 1 to 4); SNPl is the fixed effect of SNP genotype, where l was coded as 0, 1, or 2 copies of one of the alleles; Sm is the effect accounting for season (dry or rainy season); CLijklmnop is a covariate for the effect of coat length; CTijklmnop is a covariate for the effect of coat thickness; HDijklmnop is a covariate for the effect of hair density; VRRijklmnop is a covariate for the effect of ΔRR, calculated by the difference between both measurements; BCijklmnop is a covariate for the estimated indicine percent of the animal; β1, β2, β3, β4, and β5 are the regression coefficients an is a random additive genetic effect of animal n, assuming a~N(0,Aσa2), with additive pedigree relationship matrix A and the additive genetic variance σa2; peo is the random permanent environmental effects, assuming pe~N(0,Iσpe2) with identity matrix I and permanent variance σpe2; and eijklmnop is the random residual effect, assuming e~N(0,Iσe2), where σe2 is the residual variance. The variance components were estimated based on same model without the SNP effect.
      Although the additional evaluated traits show correlation with RT control, in this work we selected ΔRT as the main trait. For this reason, we included these additional traits only as covariates in the model to look for significant SNP that could be linked only to heat stress response, disregarding the effects of these additional traits on ΔRT. The descriptive statistics of phenotype and covariates are shown in Table 1.
      Table 1Descriptive statistics of phenotype and covariates included in the model to evaluate the variation of rectal temperature (ΔRT) in an experimental Gir × Holstein F2 population evaluated for heat stress response
      VariableSeason
      Season: 1 = rainy; 2 = dry.
      Parameters
      NumberMeanSDCV (%)MinimumMaximum
      ΔRT13142.650.6624.781.004.45
      23392.010.5728.220.553.40
      Coat length13201.040.2221.490.511.89
      23341.810.3821.070.793.28
      Coat thickness13223.630.7320.071.486.00
      23414.681.1524.581.739.89
      Hair density1320184.9155.7030.1282.00412.00
      2332178.3556.6431.7667.67375.00
      ΔRR
      Variation of respiration rate.
      1316104.8819.9319.0158.00154.00
      234388.3428.4832.2312.00160.00
      BC
      Breed composition.
      3690.480.0510.260.310.66
      1 Season: 1 = rainy; 2 = dry.
      2 Variation of respiration rate.
      3 Breed composition.
      After the association analysis, the genome-wide false discovery rate (FDR) was applied to avoid false positives due to multiple testing. The R package qvalue (
      • Dabney A.
      • Storey J.D.
      • Warnes G.R.
      qvalue: Q-value estimation for false discovery rate control. R package version 1.36.0.
      ) was used to provide the P-values corrected for FDR (q-values) for the SNP association tests. Associations with a q-value ≤0.05 were considered significant.
      A single-SNP GWAS for the ΔRR trait was performed with the same model used to evaluate the ΔRT trait, in which the ΔRT trait was included as covariate. Nevertheless, no significant SNP was found to be associated with ΔRR.

      QTL Regions

      We defined QTL regions based on the location of the significant SNP and the average LD (measured as the squared correlation coefficient, r2) in the evaluated population. Haploview software (
      • Barrett J.C.
      • Fry B.
      • Maller J.
      • Daly M.J.
      Haploview: Analysis and visualization of LD and haplotype maps.
      ) was used to calculate the pairwise LD between markers. Because the significant SNP identified in the GWAS analysis were located on chromosome 17, we used the LD of this chromosome to fix the QTL size. In addition, the LD between significant SNP identified in GWAS (Supplemental Figures S1 and S2; https://doi.org/10.3168/jds.2018-15305) were also used to define QTL regions. An average r2 ≥ 0.15 was found for SNP within a distance of 500 kb. Thus, all significant SNP located within a region of 500 kb were considered as belonging to the same QTL region.

      Gene Search and Generation of Gene Networks

      Putative candidate genes within QTL regions were identified based on the UMD 3.1 assembly of the bovine genome (
      • 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, Bos taurus.
      ), using the National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov).
      The gene set was used for GO enrichment analysis via the ClueGO application in Cytoscape (
      • Bindea G.
      • Mlecnik B.
      • Hackl H.
      • Charoentong P.
      • Tosolini M.
      • Kirilovsky A.
      • Fridman W.H.
      • Pagès F.
      • Trajanoski Z.
      • Galon J.
      ClueGO: A Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks.
      ) to generate the GO. This analysis was performed based on a unilateral hypergeometric test with Bonferroni correction and aims to analyze the gene set on the search for a functional GO term or pathway that might show correlation with selected genes.
      The same gene set was used to perform the analysis of regulatory sequences. For that, promoter sequences (FASTA format) were obtained by considering 3,000 bp upstream and 300 bp downstream from genes' transcription start site in the bovine genome (
      • Soares A.C.C.
      • Guimarães S.E.F.
      • Kelly M.J.
      • Fortes M.R.S.
      • Silva F.F.E.
      • Verardo L.L.
      • Mota R.
      • Moore S.
      Multiple-trait genomewide mapping and gene network analysis for scrotal circumference growth curves in Brahman cattle.
      ). These data were screened for the identification of enriched TF using the TFM-explorer software (http://bioinfo.lifl.fr/TFM/TFME/), which uses weighting matrices from the JASPAR database (http://jaspar.binf.ku.dk/;
      • Sandelin A.
      • Alkema W.
      • Engström P.
      • Wasserman W.W.
      • Lenhard B.
      JASPAR: An open-access database for eukaryotic transcription factor binding profiles.
      ) to detect potential transcription factor binding sites (TFBS) of a set of gene sequences and search for over-represented TFBS. In addition, this software extracts significant groups—TFBS regions of the selected gene sequences associated with a factor—by calculating a score function threshold that was chosen to generate a P-value ≤ 10−3 for each position in each sequence, as described in
      • Touzet H.
      • Varré J.-S.
      Efficient and accurate P-value computation for Position Weight Matrices.
      .
      The TF list was then analyzed in Cytoscape software (
      • Shannon P.
      • Markiel A.
      • Ozier O.
      • Baliga N.S.
      • Wang J.T.
      • Ramage D.
      • Amin N.
      • Schwikowski B.
      • Ideker T.
      Cytoscape: A software environment for integrated models of biomolecular interaction networks.
      ), using the Biological Networks Gene Ontology tool (BiNGO;
      • Maere S.
      • Heymans K.
      • Kuiper M.
      BiNGO: A Cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks.
      ) to determine significantly overrepresented functional GO terms. Based on the over-represented biological processes in the BiNGO tool and evidence from literature data related to the investigated traits, we were able to identify the TF most related to heat stress, allowing the generation of gene-TF interaction networks.
      Gene-TF network analysis was performed in Cytoscape software (
      • Shannon P.
      • Markiel A.
      • Ozier O.
      • Baliga N.S.
      • Wang J.T.
      • Ramage D.
      • Amin N.
      • Schwikowski B.
      • Ideker T.
      Cytoscape: A software environment for integrated models of biomolecular interaction networks.
      ) based on the number of TFBS (connections with the selected TF) in each gene to determine the associations with ΔRT, which aided in the identification of the most likely candidate genes.

      Assignment of Breed of Origin of Alleles

      To assign the breed origin of marker alleles in the F2 population, the BOA approach (
      • Vandenplas J.
      • Calus M.P.L.
      • Sevillano C.A.
      • Windig J.J.
      • Bastiaansen J.W.M.
      Assigning breed origin to alleles in crossbred animals.
      ) was used with the parameter settings recommended by
      • Sevillano C.A.
      • Vandenplas J.
      • Bastiaansen J.W.M.
      • Calus M.P.L.
      Empirical determination of breed-of-origin of alleles in three-breed cross pigs.
      . The BOA approach consists of 3 steps. (1) The purebred and crossbred animals' genotypes were phased with AlphaPhase 1.1 software (
      • Hickey J.M.
      • Kinghorn B.P.
      • Tier B.
      • Wilson J.F.
      • Dunstan N.
      • van der Werf J.H.J.
      A combined long-range phasing and long haplotype imputation method to impute phase for SNP genotypes.
      ). Haplotype phasing was performed using pedigree and 9 combinations of core and tail lengths (350:50, 250:100, 300:100, 350:100, 150:200, 200:200, 250:200, 300:200, and 350:200), in which each combination was run in both “Offset” and “NotOffset” modes, allowing each allele to be considered 18 times through different haplotypes of variable length. “Offset” analysis shifts the beginning of each core to halfway along the first core, creating 50% overlaps between cores (
      • Vandenplas J.
      • Calus M.P.L.
      • Sevillano C.A.
      • Windig J.J.
      • Bastiaansen J.W.M.
      Assigning breed origin to alleles in crossbred animals.
      ). (2) The unique haplotypes for each pure breed were determined. In this step, it is necessary to be able to observe 80% of a haplotype's copies in a specific breed, to assign its breed of origin. (3) Breed origin was assigned to each allele observed in F2 animals based on the determined unique haplotypes.
      After performing the BOA approach, 10% of F2 animals with the lowest breeding values for ΔRT were designated as low-ΔRT, and the 10% with the highest breeding values for ΔRT were designated as high-ΔRT. Considering the LD between SNP and genes, the physical positions of candidate genes were annotated and used as coordinates for selection of adjacent SNP, which were classified as HH (both alleles originated from the Holstein breed), GG (both alleles originated from the Gir breed), or HG/GH (one allele originated from the Holstein and one from the Gir breed), based on the results of the BOA approach. Finally, the alleles of each selected and classified SNP were quantified to evaluate the gametic segregation and check for prevalence of alleles from a certain breed when animals are classified as low ΔRT or high ΔRT. For that, a chi-squared test was used to evaluate the expected gametic segregation of 25% HH, 50% HG/GH, and 25% GG in both groups. Those SNP with a difference in prevalence of alleles from a certain breed, which did not show the expected gametic segregation, were considered significant (P ≤ 0.05).

      RESULTS AND DISCUSSION

      Genetic Parameters

      We found heritability estimates of 0.13 ± 0.08 and repeatability estimates of 0.29 ± 0.06 for the ΔRT trait. During the heat stress evaluation period, the RT was evaluated at 42°C and 60% RH, in which the temperature-humidity index (THI) was 97.
      • 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 heritability of 0.17 ± 0.13 for RT, recorded between 1500 and 1700 h at different THI, in Holstein cattle during summer in north Florida. In beef cattle, heritability estimates of 0.22 for RT, recorded only during summer months when the ambient temperature was >30°C, were found in Brahman animals, and estimates of 0.14 were found in crossbred animals showing genetic backgrounds of different tropical breeds (
      • Porto-Neto L.R.
      • Reverter A.
      • Prayaga K.C.
      • Chan E.K.F.
      • Johnston D.J.
      • Hawken R.J.
      • Fordyce G.
      • Garcia J.F.
      • Sonstegard T.S.
      • Bolormaa S.
      • Goddard M.E.
      • Burrow H.M.
      • Henshall J.M.
      • Lehnert S.A.
      • Barendse W.
      The genetic architecture of climatic adaptation of tropical cattle.
      ).
      Most studies use THI to estimate the presence or absence of the heat stress phenotype in a given environment (
      • Biffani S.
      • Bernabucci U.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      Short communication : Effect of heat stress on nonreturn rate of Italian Holstein cows.
      ;
      • Hagiya K.
      • Hayasaka K.
      • Yamazaki T.
      • Shirai T.
      • Osawa T.
      • Terawaki Y.
      • Nagamine Y.
      • Masuda Y.
      • Suzuki M.
      Effects of heat stress on production, somatic cell score and conception rate in Holsteins.
      ). However, evaluation of the response of each animal individually during the heat stress is relevant, as animals respond in different ways to environmental stimuli, which is accentuated in the F2 crossbred animals used in the current study because each animal possesses a unique genetic composition.

      QTL Regions

      In the current study, the GWAS for ΔRT identified 6 significant SNP (q ≤ 0.05) in the F2 population (Figure 1; Table 2). These SNP were distributed over 3 QTL regions located on Bos taurus chromosome 17 (BTA17). These QTL regions were defined based on the average LD between markers at the BTA17 and the LD between significant SNP.
      Figure thumbnail gr1
      Figure 1Manhattan plot from the association analysis of variation of rectal temperature (ΔRT) in the experimental Gir × Holstein F2 population. Significant SNP (q ≤ 0.05) are shown as blue triangles, and the gray dots denote SNP that were not significant in genome-wide association studies.
      Table 2Significant SNP located on chromosome 17 associated with variation of rectal temperature (ΔRT) in an experimental Gir × Holstein F2 population evaluated for heat stress response
      SNPQTL
      QTL region where the significant SNP is located.
      Pos
      Position on the chromosome (in Mbp).
      Effect
      Allele substitution effect.
      −log10(P-value)q-value
      False discovery rate-based q-value.
      MAF
      Minor allele frequency.
      BTA.41874.no.rs171,165,3420.155.450.050.33
      ARS.BFGL.NGS.78729171,231,1590.155.430.050.34
      ARS.BFGL.NGS.10248171,260,8510.155.430.050.34
      ARS.BFGL.NGS.58770274,292,3190.185.750.010.24
      ARS.BFGL.NGS.24012374,948,9210.175.440.030.27
      ARS.BFGL.NGS.18349374,998,3490.175.440.030.27
      1 QTL region where the significant SNP is located.
      2 Position on the chromosome (in Mbp).
      3 Allele substitution effect.
      4 False discovery rate-based q-value.
      5 Minor allele frequency.
      Several QTL for RT single records have also been reported on other chromosomes. In purebred cattle (Angus, Simmental, and Piedmontese), QTL were identified on BTA1, 8, 10, 11, 12, 20, 22, 23, 25, 26, and 27 (
      • Howard J.T.
      • Kachman S.D.
      • Snelling W.M.
      • Pollak E.J.
      • Ciobanu D.C.
      • Kuehn L.A.
      • Spangler M.L.
      Beef cattle body temperature during climatic stress: A genome-wide association study.
      ). However, the same authors annotated genes related to heat stress on BTA10, 12, 22, 23, and 25 only.
      • Porto-Neto L.R.
      • Reverter A.
      • Prayaga K.C.
      • Chan E.K.F.
      • Johnston D.J.
      • Hawken R.J.
      • Fordyce G.
      • Garcia J.F.
      • Sonstegard T.S.
      • Bolormaa S.
      • Goddard M.E.
      • Burrow H.M.
      • Henshall J.M.
      • Lehnert S.A.
      • Barendse W.
      The genetic architecture of climatic adaptation of tropical cattle.
      identified significant SNP associated with RT on BTA6 and 17 in a Brahman and tropical composite population. In a GWAS for RT single records evaluated at different THI in heat-stressed Holstein animals,
      • 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.
      identified QTL on BTA4, 5, 16, 24, and 26, using the single-step genomic BLUP approach. In a subsequent study, based on previously selected SNP related to thermotolerance,
      • 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.
      identified significant SNP associated with RT single records on BTA4, 6, and 24. Based on selected SNP located in coding regions of genes previously associated with reproduction, production or health traits, the same authors identified significant SNP associated with RT single records on BTA6, 7, 9, 10, 11, 15, 17, 18, 22 and 25. Because the regions on BTA17 were shown to be associated with RT in the current study as well as in previous studies, and also because these regions contain candidate genes that have been shown to be involved in heat stress responses (
      • Porto-Neto L.R.
      • Reverter A.
      • Prayaga K.C.
      • Chan E.K.F.
      • Johnston D.J.
      • Hawken R.J.
      • Fordyce G.
      • Garcia J.F.
      • Sonstegard T.S.
      • Bolormaa S.
      • Goddard M.E.
      • Burrow H.M.
      • Henshall J.M.
      • Lehnert S.A.
      • Barendse W.
      The genetic architecture of climatic adaptation of tropical 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.
      ), it is possible that these BTA17 regions may have a role in cattle heat stress.

      Candidate Genes

      In the current study, 54 genes encoding proteins with different functions were annotated for ΔRT. The gene sets were grouped according to their GO terms enrichment in a GO network to better understand their functions (Figure 2). In this network, we identified genes involved in Janus kinase/signal transducers and activators of transcription (JAK/STAT) pathway, the interleukin 6 family cytokine (LIF), oncostatin M (OSM), the DGCR8 microprocessor complex subunit (DGCR8), and a gene linked to the response to oxygen radical biological process, thioredoxin reductase 2 (TXNRD2). The OSM and LIF genes are cytokines that belong to the interleukin 6 (IL6) family, and they are closely related to the structures and functions of this cytokine family (
      • Tanaka M.
      • Miyahima A.
      Onconstatin M, a multifunctional cytokine.
      ). Upon binding to cell surface receptors, IL6-type cytokines activate members of the JAK tyrosine kinase family, resulting in the activation of members of the STAT TF family, which are involved in the JAK/STAT pathway (
      • Hong D.S.
      • Banerji U.
      • Tavana B.
      • George G.C.
      • Aaron J.
      • Kurzrock R.
      Targeting the molecular chaperone heat shock protein 90 (HSP90): Lessons learned and future directions.
      ). During heat shock, this signaling pathway can be quickly activated by OSM and LIF cytokines, aiding in reduction of the effects stimulated by heat shock as well as in the regulation of the expression of heat shock protein (HSP) genes Hsp70 and Hsp90 (
      • Chatterjee M.
      • Jain S.
      • Stühmer T.
      • Andrulis M.
      • Ungethüm U.
      • Kuban R.J.
      • Lorentz H.
      • Bommert K.
      • Topp M.
      • Krämer D.
      • Müller-Hermelink H.K.
      • Einsele H.
      • Greiner A.
      • Bargou R.C.
      STAT3 and MAPK signaling maintain overexpression of heat shock proteins 90α and β in multiple myeloma cells, which critically contribute to tumor-cell survival.
      ;
      • Allegra A.
      • Sant'Antonio E.
      • Penna G.
      • Alonci A.
      • D'Angelo A.
      • Russo S.
      • Cannavò A.
      • Gerace D.
      • Musolino C.
      Novel therapeutic strategies in multiple myeloma: Role of the heat shock protein inhibitors.
      ;
      • Stephanou A.
      • Latchman D.S.
      Transcriptional modulation of heat-shock protein gene expression.
      ). In particular, STAT III shows important biological functions associated with apoptosis control, cell survival, and activation and inhibition of the immune and inflammatory responses (
      • Loor J.J.
      Genomics of metabolic adaptations in the peripartal cow.
      ).
      Figure thumbnail gr2
      Figure 2Functional networks showing gene interactions (triangles nodes) related to variation of rectal temperature (ΔRT) and the relationship across genes and their sub-networks related to positive regulation of tyrosine phosphorylation of signal transducer and activator of transcription family 3 (STAT 3) protein in the experimental Gir × Holstein F2 population. JAK = Janus kinase.
      The DGCR8 gene, together with nuclear RNase III enzyme Drosha, comprise the microprocessor complex, which is essential in the processing of microRNA in animals (
      • Triboulet R.
      • Chang H.M.
      • Lapierre R.J.
      • Gregory R.I.
      Post-transcriptional control of DGCR8 expression by the Microprocessor.
      ;
      • Faller M.
      • Toso D.
      • Matsunaga M.
      • Atanasov I.
      • Senturia R.
      • Chen Y.
      • Zhou Z.H.
      • Feng G.
      DGCR8 recognizes primary transcripts of microRNAs through highly cooperative binding and formation of higher-order structures.
      ). According to
      • Knuckles P.
      • Carl S.H.
      • Musheev M.
      • Niehrs C.
      • Wenger A.
      • Bühler M.
      RNA fate determination through cotranscriptional adenosine methylation and microprocessor binding.
      , DGCR8 and the N6-adenosine-methyltransferase Mettl3 are reallocated to heat-shock genes under acute temperature stress. They work in consonance to mark Hsp70 mRNA, allowing subsequent RNA degradation to control the timing and magnitude of the heat-shock response. Previous studies described a significant increase in the protein expression levels of HSP in cattle subjected to heat stress (
      • Belhadj Slimen I.
      • Najar T.
      • Ghram A.
      • Abdrrabba M.
      Heat stress effects on livestock: Molecular, cellular and metabolic aspects, A review.
      ;
      • Hu H.
      • Zhang Y.
      • Zheng N.
      • Cheng J.
      • Wang J.
      The effect of heat stress on gene expression and synthesis of heat-shock and milk proteins in bovine mammary epithelial cells.
      ). Clearing of HSP after heat stress is very important, to avoid accumulation of the proteins and presumably worse recovery after stress (
      • Knuckles P.
      • Carl S.H.
      • Musheev M.
      • Niehrs C.
      • Wenger A.
      • Bühler M.
      RNA fate determination through cotranscriptional adenosine methylation and microprocessor binding.
      ).
      The TXNRD2 gene, together with the thioredoxin 2 (TRX2) and peroxiredoxin 3 (PRX3) genes, comprise the mitochondrial thioredoxin system, which provides a primary line of redox regulation through high control over reactive oxygen species (ROS) emission from mitochondria (
      • Stanley B.A.
      • Sivakumaran V.
      • Shi S.
      • McDonald I.
      • Lloyd D.
      • Watson W.H.
      • Aon M.A.
      • Paolocci N.
      Thioredoxin reductase-2 is essential for keeping low levels of H2O2 emission from isolated heart mitochondria.
      ;
      • Yoshioka J.
      Thioredoxin reductase 2 (Txnrd2) regulates mitochondrial integrity in the progression of age-related heart failure.
      ), which is a central source of H2O2. Animals exposed to heat stress show an increase in RR, and consequently H2O2 production exceeds its scavenging, leading to oxidative stress (
      • Bernabucci U.
      • Ronchi B.
      • Lacetera N.
      • Nardone A.
      Markers of oxidative status in plasma and erythrocytes of transition dairy cows during hot season.
      ;
      • Srikandakumar A.
      • Johnson E.H.
      Effect of heat stress on milk production, rectal temperature, respiratory rate and blood chemistry in Holstein, Jersey and Australian Milking Zebu cows.
      ). Thioredoxin reductase 2 (TRXR2) supplies nicotinamide adenine dinucleotide phosphate (NADPH) electrons to PRX3, which controls H2O2 levels and establishes mitochondrial redox homeostasis (
      • Stanley B.A.
      • Sivakumaran V.
      • Shi S.
      • McDonald I.
      • Lloyd D.
      • Watson W.H.
      • Aon M.A.
      • Paolocci N.
      Thioredoxin reductase-2 is essential for keeping low levels of H2O2 emission from isolated heart mitochondria.
      ;
      • Aon M.A.
      • Stanley B.A.
      • Sivakumaran V.
      • Kembro J.M.
      • O'Rourke B.
      • Paolocci N.
      • Cortassa S.
      Glutathione/thioredoxin systems modulate mitochondrial H2O2 emission: An experimental-computational study.
      ). Inhibition of TRXR2 leads to impaired redox homeostasis and results in increased levels of H2O2, endangering cell functions (
      • Prasad R.
      • Chan L.F.
      • Hughes C.R.
      • Kaski J.P.
      • Kowalczyk J.C.
      • Savage M.O.
      • Peters C.J.
      • Nathwani N.
      • Clark A.J.L.
      • Storr H.L.
      • Metherell L.A.
      Thioredoxin reductase 2 (TXNRD2) mutation associated with familial glucocorticoid deficiency (FGD).
      ).
      The same gene set was used to perform the analysis of regulatory sequences, through which, based on the biological process and literature review, we were able to identify 3 heat stress–associated TF: nuclear factor kappa B (NF-κB), aryl hydrocarbon receptor nuclear translocator (ARNT), and STAT III (Table 3). A gene-TF network was generated with the 3 selected TF and genes that show potential binding sites for these TF (Figure 3).
      Table 3Enriched transcription factors (TF) associated with genes identified for cattle variation of rectal temperature (ΔRT), based on gene ontology biological process and literature review; data obtained from an experimental Gir × Holstein F2 population evaluated for heat stress response
      TF
      ARNT = aryl hydrocarbon receptor nuclear translocator; NF-κB = nuclear factor kappa B; STAT III = signal transducer and activator of transcription family 3.
      Biological processLiterature review
      Cited studies are just a sample of the vast available literature.
      ARNTRegulation of homeostatic processInvolved in negative regulation of necrosis; can be regulated by heat shock proteins before and after hypoxia (
      • Belenichev I.F.
      • Kolesnik Y.M.
      • Pavlov S.V.
      • Sokolik E.P.
      • Bukhtiyarova N.V.
      Disturbance of HSP70 chaperone activity is a possible mechanism of mitochondrial dysfunction.
      )
      NF-κBAnti-apoptosisHeat stress can inhibit activity of NF-kB and induce massive apoptosis (
      • Belardo G.
      • Piva R.
      • Santoro M.G.
      Heat stress triggers apoptosis by impairing NF-kappaB survival signaling in malignant B cells.
      )
      STAT IIIHomeostatic processCell survival regulator (
      • Terui K.
      • Enosawa S.
      • Haga S.
      • Zhang H.Q.
      • Kuroda H.
      • Kouchi K.
      • Matsunaga T.
      • Yoshida H.
      • Engelhardt J.F.
      • Irani K.
      • Ohnuma N.
      • Ozaki M.
      Stat3 confers resistance against hypoxia/reoxygenation-induced oxidative injury in hepatocytes through upregulation of Mn-SOD.
      ) that participates in intracellular ROS homeostasis (
      • He G.
      • Yu G.Y.
      • Temkin V.
      • Ogata H.
      • Kuntzen C.
      • Sakurai T.
      • Sieghart W.
      • Peck-Radosavljevic M.
      • Leffert H.L.
      • Karin M.
      Hepatocyte IKKb/NF-kB inhibits tumor promotion and progression by preventing oxidative stress-driven STAT3 activation.
      ;
      • Li L.
      • Cheung S.-h.
      • Evans E.L.
      • Shaw P.E.
      Modulation of gene expression and tumor cell growth by redox modification of STAT3.
      )
      1 ARNT = aryl hydrocarbon receptor nuclear translocator; NF-κB = nuclear factor kappa B; STAT III = signal transducer and activator of transcription family 3.
      2 Cited studies are just a sample of the vast available literature.
      Figure thumbnail gr3
      Figure 3Gene-transcription factor (TF) network: genes located in the QTL regions associated with variation of rectal temperature (ΔRT; green circle nodes) and their associated TF (yellow octagon nodes). Node size corresponds to the network analyses (Cytoscape;
      • Shannon P.
      • Markiel A.
      • Ozier O.
      • Baliga N.S.
      • Wang J.T.
      • Ramage D.
      • Amin N.
      • Schwikowski B.
      • Ideker T.
      Cytoscape: A software environment for integrated models of biomolecular interaction networks.
      ) in which larger nodes denote a higher edge density associated with the number of TF binding sites. Blue square nodes show the gene ontology biological processes related to TF. Data were obtained from an experimental Gir × Holstein F2 population evaluated for heat stress response. ARNT = aryl hydrocarbon receptor nuclear translocator; NF-κB = nuclear factor kappa B; STAT III = signal transducer and activator of transcription family 3.
      The most enriched TF was NF-κB, linked to the highest number of genes. The response to heat shock has been implicated in the negative regulation of the NF-κB signaling pathway (
      • Yenari M.A.
      • Liu J.
      • Zheng Z.
      • Vexler Z.S.
      • Lee J.E.
      • Giffard R.G.
      Antiapoptotic and anti-inflammatory mechanisms of heat-shock protein protection.
      ). According to
      • Belardo G.
      • Piva R.
      • Santoro M.G.
      Heat stress triggers apoptosis by impairing NF-kappaB survival signaling in malignant B cells.
      , the effects of hyperthermic treatment are related to increase in temperature above physiological conditions and duration of exposure. These authors observed that heat stress induced an inhibition of constitutive NF-κB activity in HS-Sultan cells, a B-cell lymphoma directly related to temperature increase, in addition to inhibiting the activity of NF-κB and inducing massive apoptosis in other types of aggressive B-cell neoplasms
      • Belardo G.
      • Piva R.
      • Santoro M.G.
      Heat stress triggers apoptosis by impairing NF-kappaB survival signaling in malignant B cells.
      . In a study regarding crosstalk among cytokines and induced hyperthermia using genomic approaches,
      • Janus P.
      • Stokowy T.
      • Jaksik R.
      • Szoltysek K.
      • Handschuh L.
      • Podkowinski J.
      • Widlak W.
      • Kimmel M.
      • Widlak P.
      Cross talk between cytokine and hyperthermia-induced pathways: identification of different subsets of NF-κB-dependent genes regulated by TNFα and heat shock.
      found a substantial overlap among the set of genes potentially regulated by heat shock protein or heat shock transcription factor 1 (HSFI) and the set of genes potentially regulated by tumor necrosis factor alpha (TNFα) or NF-κB. This allowed the identification of predictable inhibitory effects of heat shock on the expression of classical NF-κB target genes and new patterns of activation (or co-activation) related to responses to different types of stress and hormonal stimuli. In addition, it helped to identify repression (or co-repression) in subsets of genes associated with processes mediated by NF-κB signaling and humoral immune responses, as well as chemokine and cytokine responses. Thus, pleiotropic effects of heat stress on the regulation of NF-κB–dependent genes should be expected.
      The ARNT TF, also known as hypoxia inducible factor 1 subunit beta (HIFIB), is a subunit of transcription factor hypoxia-inducible factor 1 (hif-1;
      • Ahn J.E.
      • Zhou X.
      • Dowd S.E.
      • Chapkin R.S.
      • Zhu-Salzman K.
      Insight into hypoxia tolerance in cowpea bruchid: Metabolic repression and heat shock protein regulation via hypoxia-inducible factor 1.
      ), a key regulator of the developmental and physiological networks required for the maintenance of O2 homeostasis. Past studies have found that hif-1 plays an important role in preventing overproduction of mitochondrial ROS in hypoxia conditions (
      • Semenza G.L.
      Hypoxia-inducible factor 1: Regulator of mitochondrial metabolism and mediator of ischemic preconditioning.
      ). After preconditioning at intermediate temperatures, HSFI promotes survival under extreme thermal stress (
      • Kourtis N.
      • Nikoletopoulou V.
      • Tavernarakis N.
      Small heat-shock proteins protect from heat-stroke-associated neurodegeneration.
      ), and according to
      • Horowitz M.
      • Assadi H.
      Heat acclimation-mediated cross-tolerance in cardioprotection: Do HSP70 and HIF-1alpha play a role?.
      its expression increased after heat acclimation-mediated cross-tolerance. Meanwhile, ARNT is involved in the negative regulation of necrosis and can be stabilized and have its expression regulated as well as its lifetime increased by heat shock proteins before and after hypoxia (
      • Belenichev I.F.
      • Kolesnik Y.M.
      • Pavlov S.V.
      • Sokolik E.P.
      • Bukhtiyarova N.V.
      Disturbance of HSP70 chaperone activity is a possible mechanism of mitochondrial dysfunction.
      ).
      The less enriched TF in the gene-TF network (Figure 3), STAT III regulates the expression of Hsp70 and Hsp90 genes, as previously described; it is also an important regulator of cellular survival after apoptotic stimuli (
      • Terui K.
      • Enosawa S.
      • Haga S.
      • Zhang H.Q.
      • Kuroda H.
      • Kouchi K.
      • Matsunaga T.
      • Yoshida H.
      • Engelhardt J.F.
      • Irani K.
      • Ohnuma N.
      • Ozaki M.
      Stat3 confers resistance against hypoxia/reoxygenation-induced oxidative injury in hepatocytes through upregulation of Mn-SOD.
      ) and directly protects cells from oxidative stress (
      • Li L.
      • Wei W.
      • Zhang Y.
      • Tu G.
      • Zhang Y.
      • Yang J.
      • Xing Y.
      SirT1 and STAT3 protect retinal pigmented epithelium cells against oxidative stress.
      ). This TF is sensitive to intracellular oxidants and can be activated in response to ROS accumulation. Furthermore, it participates in intracellular ROS homeostasis and shows an inverse relationship with NF-κB (
      • He G.
      • Yu G.Y.
      • Temkin V.
      • Ogata H.
      • Kuntzen C.
      • Sakurai T.
      • Sieghart W.
      • Peck-Radosavljevic M.
      • Leffert H.L.
      • Karin M.
      Hepatocyte IKKb/NF-kB inhibits tumor promotion and progression by preventing oxidative stress-driven STAT3 activation.
      ;
      • Li L.
      • Cheung S.-h.
      • Evans E.L.
      • Shaw P.E.
      Modulation of gene expression and tumor cell growth by redox modification of STAT3.
      ).
      From the gene-TF network, candidate genes for ΔRT were highlighted based on TFBS. The LIF, DGCR8, and TXNRD2 genes, which also share ontology in the GO enrichment analysis, are part of the most enriched gene group in this network (Figure 3).

      Breed of Origin of Alleles for Candidate Genes

      Based on the GO and gene-TF networks, we were able to select 4 candidate genes to assess the origin of alleles associated with ΔRT in cattle. Considering the LD between SNP and genes, the physical positions of LIF, OSM, TXNRD2, and DGCR8 genes were annotated and used as coordinates for selection of adjacent SNP. We selected a total of 7 SNP: 2 associated with the LIF and OSM genes, 3 associated with the TXNRD2 gene, and 2 associated with the DGCR8 gene. The animals were grouped as low ΔRT or high ΔRT based on their genomic breeding values for ΔRT, and the alleles of the selected SNP were classified and quantified according to their breed of origin (HH, GG, or HG/GH).
      Based on chi-squared test results obtained from low- and high-ΔRT F2 animals, we observed that all selected SNP for allele origin were significant (P < 0.05) in both evaluated groups (Table 4) and showed a deviation of the expected gametic segregation for each genotypic class: 1HH to 2HG/GH to 1GG. These results suggest a prevalence of alleles from one breed in the low- and high-ΔRT groups. On average, 53% of low-ΔRT animals had 2 alleles from the Holstein breed, whereas 6.4% had 2 alleles from the Gir breed. On the other hand, in high-ΔRT animals, the average proportion of animals for HH and GG were 2.9% and 49.9%, respectively. Based on these results, we can observe that Holstein alleles are more frequent in low-ΔRT animals, and Gir alleles are more abundant in high-ΔRT animals.
      Table 4Percentage of individuals classified by breed of origin of each SNP allele located in the candidate genes (LIF, OSM, TXNRD2, and DGCR8) associated with variation of rectal temperature (ΔRT) in low- and high-ΔRT animals from an experimental Gir × Holstein F2 population evaluated for heat stress response
      All SNP are significant for allele origin in the chi-squared test (P < 0.05).
      AlleleLIF/OSMTXNRD2DGCR8
      SNP1SNP2SNP3SNP4SNP5SNP6SNP7
      Low-ΔRT animals
       HH
      Both alleles originated from Holstein breed.
      55.355.351.451.451.454.551.4
       GG
      Both alleles originated from Gir breed.
      7.97.95.75.75.76.15.7
       HG/GH
      One allele originated from Holstein and one from Gir breed.
      36.836.842.942.942.939.442.9
      High-ΔRT animals
       HH10.010.00.00.00.00.00.0
       GG52.552.548.148.150.048.150.0
       HG/GH37.537.551.951.950.051.950.0
      1 All SNP are significant for allele origin in the chi-squared test (P < 0.05).
      2 Both alleles originated from Holstein breed.
      3 Both alleles originated from Gir breed.
      4 One allele originated from Holstein and one from Gir breed.
      It is known that Bos indicus animals have a greater ability to regulate their body temperature, which helps in heat stress prevention (
      • Hansen P.J.
      Physiological and cellular adaptations of zebu cattle to thermal stress.
      ;
      • Cattelam J.
      • Martinez M.
      Estresse térmico em bovinos.
      ) compared with Bos taurus animals, which are affected by heat stress (
      • Renaudeau D.
      • Collin A.
      • Yahav S.
      • de Basilio V.
      • Gourdine J.L.
      • Collier R.J.
      Adaptation to hot climate and strategies to alleviate heat stress in livestock production.
      ). In the heat chamber, F2 animals were subjected to high temperatures and RH, which may have prejudiced the evaporation mechanisms, responsible for 80% of body heat loss to the environment (
      • Shearer J.K.
      • Beede D.K.
      Thermoregulation and physiological responses of dairy cattle in hot weather.
      ;
      • West J.W.
      Effects of heat-stress on production in dairy cattle.
      ), and favored the initiation of heat stress. Given that the genes identified in our study play important functions in controlling the effects of current heat stress, and not heat stress prevention per se, it is reasonable to assume that the significant SNP alleles were derived from the Holstein breed, which requires a more efficient genetic architecture than Gir animals to defend the body from deleterious effects of heat stress.

      CONCLUSIONS

      In summary, the ΔRT trait is a relevant phenotype for individual evaluation of heat stress response. A GWAS allowed the identification of significant SNP associated with ΔRT in the experimental F2 Gir × Holstein population. Gene-TF networks and GO enrichment analyses allowed identification of candidate genes (LIF, OSM, TXNRD2, and DGCR8) that are biologically related to heat stress response. According to BOA analysis, Holstein breed alleles could be associated with a more complex response to heat stress effects, which can be explained by the fact that Holstein animals are more affected by heat stress than Gir animals, requiring an intricate genetic architecture to defend the body from the deleterious effects of heat stress. Future studies are needed to further our understanding of the genetic architecture underlying the heat stress response in crossbred cattle.

      ACKNOWLEDGMENTS

      We thank EMBRAPA Dairy Cattle Research Center, Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), and the Ministério da Ciência Tecnologia e Inovação (MCTI)/Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq)/Instituto Nacional de Ciência e Tecnologia—Ciência Animal (INCT-CA) for their financial support, Embrapa Dairy Cattle Research Center for providing the data needed to carry out this study, and Mario P. L. Calus (Wageningen University and Research Animal Breeding and Genomics, 6700 AH Wageningen, the Netherlands) for his assistance with BOA methodology. The authors thank and acknowledge all the team effort and the visionary work done by the late Mário Luiz Martinez, former leader of this project at Embrapa Dairy Cattle Research Center.

      Supplementary Material

      REFERENCES

        • Ahn J.E.
        • Zhou X.
        • Dowd S.E.
        • Chapkin R.S.
        • Zhu-Salzman K.
        Insight into hypoxia tolerance in cowpea bruchid: Metabolic repression and heat shock protein regulation via hypoxia-inducible factor 1.
        PLoS One. 2013; 8 (23593115)
        • Alexander D.H.
        • Novembre J.
        • Lange K.
        Fast model-based estimation of ancestry in unrelated individuals.
        Genome Res. 2009; 19 (19648217): 1655-1664
        • Allegra A.
        • Sant'Antonio E.
        • Penna G.
        • Alonci A.
        • D'Angelo A.
        • Russo S.
        • Cannavò A.
        • Gerace D.
        • Musolino C.
        Novel therapeutic strategies in multiple myeloma: Role of the heat shock protein inhibitors.
        Eur. J. Haematol. 2011; 86 (21114539): 93-110
        • Aon M.A.
        • Stanley B.A.
        • Sivakumaran V.
        • Kembro J.M.
        • O'Rourke B.
        • Paolocci N.
        • Cortassa S.
        Glutathione/thioredoxin systems modulate mitochondrial H2O2 emission: An experimental-computational study.
        J. Gen. Physiol. 2012; 139 (22585969): 479-491
        • Barrett J.C.
        • Fry B.
        • Maller J.
        • Daly M.J.
        Haploview: Analysis and visualization of LD and haplotype maps.
        Bioinformatics. 2005; 21 (15297300): 263-265
        • Belardo G.
        • Piva R.
        • Santoro M.G.
        Heat stress triggers apoptosis by impairing NF-kappaB survival signaling in malignant B cells.
        Leukemia. 2010; 24 (19924145): 187-196
        • Belenichev I.F.
        • Kolesnik Y.M.
        • Pavlov S.V.
        • Sokolik E.P.
        • Bukhtiyarova N.V.
        Disturbance of HSP70 chaperone activity is a possible mechanism of mitochondrial dysfunction.
        Neurochem. J. 2011; 5: 251-256
        • Belhadj Slimen I.
        • Najar T.
        • Ghram A.
        • Abdrrabba M.
        Heat stress effects on livestock: Molecular, cellular and metabolic aspects, A review.
        J. Anim. Physiol. Anim. Nutr. (Berl.). 2016; 100 (26250521): 401-412
        • Bernabucci U.
        • Ronchi B.
        • Lacetera N.
        • Nardone A.
        Markers of oxidative status in plasma and erythrocytes of transition dairy cows during hot season.
        J. Dairy Sci. 2002; 85 (12362449): 2173-2179
        • 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; (27108174): 5837-5843
        • Bindea G.
        • Mlecnik B.
        • Hackl H.
        • Charoentong P.
        • Tosolini M.
        • Kirilovsky A.
        • Fridman W.H.
        • Pagès F.
        • Trajanoski Z.
        • Galon J.
        ClueGO: A Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks.
        Bioinformatics. 2009; 25 (19237447): 1091-1093
        • Cattelam J.
        • Martinez M.
        Estresse térmico em bovinos.
        Rev. Port. Ciências Veterinárias. 2013; 108: 96-102
        • Chatterjee M.
        • Jain S.
        • Stühmer T.
        • Andrulis M.
        • Ungethüm U.
        • Kuban R.J.
        • Lorentz H.
        • Bommert K.
        • Topp M.
        • Krämer D.
        • Müller-Hermelink H.K.
        • Einsele H.
        • Greiner A.
        • Bargou R.C.
        STAT3 and MAPK signaling maintain overexpression of heat shock proteins 90α and β in multiple myeloma cells, which critically contribute to tumor-cell survival.
        Blood. 2007; 109 (17003370): 720-728
      1. da Cruz, L. V., D. de S. R. Angrimani, B. R. Rui, and M. A. da Silva. 2011. Efeitos do estresse térmico na produção leiteira: revisão de literatura. Rev. Científica Eletrônica Med. Veterinária.

        • Dabney A.
        • Storey J.D.
        • Warnes G.R.
        qvalue: Q-value estimation for false discovery rate control. R package version 1.36.0.
        • 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: 409-419
        • Faller M.
        • Toso D.
        • Matsunaga M.
        • Atanasov I.
        • Senturia R.
        • Chen Y.
        • Zhou Z.H.
        • Feng G.
        DGCR8 recognizes primary transcripts of microRNAs through highly cooperative binding and formation of higher-order structures.
        RNA. 2010; 16 (20558544): 1570-1583
        • de Andrade Ferrazza R.
        • Mogollón Garcia H.D.
        • Vallejo Aristizábal V.H.
        • de Souza Nogueira C.
        • Veríssimo C.J.
        • Sartori J.R.
        • Sartori R.
        • Pinhiero Ferreira J.C.
        Thermoregulatory responses of Holstein cows exposed to experimentally induced heat stress.
        J. Therm. Biol. 2017; 66 (28477912): 68-80
        • Gilmour A.R.
        • Gogel B.
        • Cullis B.
        • Thompson R.
        ASReml user guide release 3.0.
        VSN International Ltd., Hemel Hempstead, UK2009
        • Hagiya K.
        • Hayasaka K.
        • Yamazaki T.
        • Shirai T.
        • Osawa T.
        • Terawaki Y.
        • Nagamine Y.
        • Masuda Y.
        • Suzuki M.
        Effects of heat stress on production, somatic cell score and conception rate in Holsteins.
        Anim. Sci. J. 2017; 88 (27113198): 3-10
        • Hansen P.J.
        Physiological and cellular adaptations of zebu cattle to thermal stress.
        Anim. Reprod. Sci. 2004; 82–83 (15271465): 349-360
        • He G.
        • Yu G.Y.
        • Temkin V.
        • Ogata H.
        • Kuntzen C.
        • Sakurai T.
        • Sieghart W.
        • Peck-Radosavljevic M.
        • Leffert H.L.
        • Karin M.
        Hepatocyte IKKb/NF-kB inhibits tumor promotion and progression by preventing oxidative stress-driven STAT3 activation.
        Cancer Cell. 2010; 17 (20227042): 286-297
        • Hickey J.M.
        • Kinghorn B.P.
        • Tier B.
        • Wilson J.F.
        • Dunstan N.
        • van der Werf J.H.J.
        A combined long-range phasing and long haplotype imputation method to impute phase for SNP genotypes.
        Genet. Sel. Evol. 2011; 43 (21388557): 12
        • Hong D.S.
        • Banerji U.
        • Tavana B.
        • George G.C.
        • Aaron J.
        • Kurzrock R.
        Targeting the molecular chaperone heat shock protein 90 (HSP90): Lessons learned and future directions.
        Cancer Treat. Rev. 2013; 39 (23199899): 375-387
        • Horowitz M.
        • Assadi H.
        Heat acclimation-mediated cross-tolerance in cardioprotection: Do HSP70 and HIF-1alpha play a role?.
        Ann. N. Y. Acad. Sci. 2010; 1188 (20201904): 199-206
        • Howard J.T.
        • Kachman S.D.
        • Snelling W.M.
        • Pollak E.J.
        • Ciobanu D.C.
        • Kuehn L.A.
        • Spangler M.L.
        Beef cattle body temperature during climatic stress: A genome-wide association study.
        Int. J. Biometeorol. 2014; 58 (24362770): 1665-1672
        • Hu H.
        • Zhang Y.
        • Zheng N.
        • Cheng J.
        • Wang J.
        The effect of heat stress on gene expression and synthesis of heat-shock and milk proteins in bovine mammary epithelial cells.
        Anim. Sci. J. 2016; 87 (26467738): 84-91
        • Janus P.
        • Stokowy T.
        • Jaksik R.
        • Szoltysek K.
        • Handschuh L.
        • Podkowinski J.
        • Widlak W.
        • Kimmel M.
        • Widlak P.
        Cross talk between cytokine and hyperthermia-induced pathways: identification of different subsets of NF-κB-dependent genes regulated by TNFα and heat shock.
        Mol. Genet. Genomics. 2015; 290 (25944781): 1979-1990
        • Knuckles P.
        • Carl S.H.
        • Musheev M.
        • Niehrs C.
        • Wenger A.
        • Bühler M.
        RNA fate determination through cotranscriptional adenosine methylation and microprocessor binding.
        Nat. Struct. Mol. Biol. 2017; 24 (28581511): 561-569
        • Kourtis N.
        • Nikoletopoulou V.
        • Tavernarakis N.
        Small heat-shock proteins protect from heat-stroke-associated neurodegeneration.
        Nature. 2012; 490 (22972192): 213-218
        • Li L.
        • Cheung S.-h.
        • Evans E.L.
        • Shaw P.E.
        Modulation of gene expression and tumor cell growth by redox modification of STAT3.
        Cancer Res. 2010; 70 (20807804): 8222-8232
        • Li L.
        • Wei W.
        • Zhang Y.
        • Tu G.
        • Zhang Y.
        • Yang J.
        • Xing Y.
        SirT1 and STAT3 protect retinal pigmented epithelium cells against oxidative stress.
        Mol. Med. Rep. 2015; 12 (25847123): 2231-2238
        • Loor J.J.
        Genomics of metabolic adaptations in the peripartal cow.
        Animal. 2010; 4 (22444613): 1110-1139
        • Machado M.A.
        • Azevedo A.L.S.
        • Teodoro R.L.
        • Pires M.A.
        • Peixoto M.G.C.D.
        • de Freitas C.
        • Prata M.C.A.
        • Furlong J.
        • da Silva M.V.G.B.
        • Guimarães S.E.F.
        • Regitano L.C.A.
        • Coutinho L.L.
        • Gasparin G.
        • Verneque R.S.
        Genome wide scan for quantitative trait loci affecting tick resistance in cattle (Bos taurus × Bos indicus).
        BMC Genomics. 2010; 11 (20433753): 280
        • Maere S.
        • Heymans K.
        • Kuiper M.
        BiNGO: A Cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks.
        Bioinformatics. 2005; 21 (15972284): 3448-3449
        • McDowell R.E.
        • Hooven N.W.
        • Camoens J.K.
        Effect of climate on performance of Holsteins in first lactation.
        J. Dairy Sci. 1975; 59: 965-971
        • 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
        • Otto P.I.
        • Guimarães S.E.F.
        • Verardo L.L.
        • Azevedo A.L.S.
        • Vandenplas J.
        • Soares A.C.C.
        • Sevillano C.A.
        • Veroneze R.
        • de Fatima M.
        • Pires A.
        • de Freitas C.
        • Prata M.C.A.
        • Furlong J.
        • Verneque R.S.
        • Martins M.F.
        • Panetto J.C.C.
        • Carvalho W.A.
        • Gobo D.O.R.
        • da Silva M.V.G.B.
        • Machado M.A.
        Genome-wide association studies for tick resistance in Bos taurus × Bos indicus crossbred cattle: A deeper look into this intricate mechanism.
        J. Dairy Sci. 2018; 101 (30243625): 11020-11032
        • Porto-Neto L.R.
        • Reverter A.
        • Prayaga K.C.
        • Chan E.K.F.
        • Johnston D.J.
        • Hawken R.J.
        • Fordyce G.
        • Garcia J.F.
        • Sonstegard T.S.
        • Bolormaa S.
        • Goddard M.E.
        • Burrow H.M.
        • Henshall J.M.
        • Lehnert S.A.
        • Barendse W.
        The genetic architecture of climatic adaptation of tropical cattle.
        PLoS One. 2014; 9 (25419663)e113284
        • Prasad R.
        • Chan L.F.
        • Hughes C.R.
        • Kaski J.P.
        • Kowalczyk J.C.
        • Savage M.O.
        • Peters C.J.
        • Nathwani N.
        • Clark A.J.L.
        • Storr H.L.
        • Metherell L.A.
        Thioredoxin reductase 2 (TXNRD2) mutation associated with familial glucocorticoid deficiency (FGD).
        J. Clin. Endocrinol. Metab. 2014; 99 (24601690): E1556-E1563
        • Renaudeau D.
        • Collin A.
        • Yahav S.
        • de Basilio V.
        • Gourdine J.L.
        • Collier R.J.
        Adaptation to hot climate and strategies to alleviate heat stress in livestock production.
        Animal. 2012; 6 (22558920): 707-728
        • Sandelin A.
        • Alkema W.
        • Engström P.
        • Wasserman W.W.
        • Lenhard B.
        JASPAR: An open-access database for eukaryotic transcription factor binding profiles.
        Nucleic Acids Res. 2004; 32 (14681366): D91-D94
        • Semenza G.L.
        Hypoxia-inducible factor 1: Regulator of mitochondrial metabolism and mediator of ischemic preconditioning.
        Biochim. Biophys. Acta. 2011; 1813 (20732359): 1263-1268
        • Sevillano C.A.
        • Napel J.
        • Guimarães S.E.F.
        • Silva F.F.
        • Calus M.P.L.
        Effects of alleles in crossbred pigs estimated for genomic prediction depend on their breed-of-origin. 2018; 19: 740
        • Sevillano C.A.
        • Vandenplas J.
        • Bastiaansen J.W.M.
        • Bergsma R.
        • Calus M.P.L.
        Genomic evaluation for a three-way crossbreeding system considering breed-of-origin of alleles.
        Genet. Sel. Evol. 2017; 49 (29061123): 75
        • Sevillano C.A.
        • Vandenplas J.
        • Bastiaansen J.W.M.
        • Calus M.P.L.
        Empirical determination of breed-of-origin of alleles in three-breed cross pigs.
        Genet. Sel. Evol. 2016; 48 (27491547): 55
        • Shannon P.
        • Markiel A.
        • Ozier O.
        • Baliga N.S.
        • Wang J.T.
        • Ramage D.
        • Amin N.
        • Schwikowski B.
        • Ideker T.
        Cytoscape: A software environment for integrated models of biomolecular interaction networks.
        Genome Res. 2003; 13 (14597658): 2498-2504
        • Shearer J.K.
        • Beede D.K.
        Thermoregulation and physiological responses of dairy cattle in hot weather.
        Agri-Practice. 1990; 11: 5-17
        • Soares A.C.C.
        • Guimarães S.E.F.
        • Kelly M.J.
        • Fortes M.R.S.
        • Silva F.F.E.
        • Verardo L.L.
        • Mota R.
        • Moore S.
        Multiple-trait genomewide mapping and gene network analysis for scrotal circumference growth curves in Brahman cattle.
        J. Anim. Sci. 2017; 95 (28805926): 3331-3345
        • Solé X.
        • Guinó E.
        • Valls J.
        • Iniesta R.
        • Moreno V.
        SNPStats: A web tool for the analysis of association studies.
        Bioinformatics. 2006; 22 (16720584): 1928-1929
        • Srikandakumar A.
        • Johnson E.H.
        Effect of heat stress on milk production, rectal temperature, respiratory rate and blood chemistry in Holstein, Jersey and Australian Milking Zebu cows.
        Trop. Anim. Health Prod. 2004; 36 (15563029): 685-692
        • Stanley B.A.
        • Sivakumaran V.
        • Shi S.
        • McDonald I.
        • Lloyd D.
        • Watson W.H.
        • Aon M.A.
        • Paolocci N.
        Thioredoxin reductase-2 is essential for keeping low levels of H2O2 emission from isolated heart mitochondria.
        J. Biol. Chem. 2011; 286 (21832082): 33669-33677
        • Stephanou A.
        • Latchman D.S.
        Transcriptional modulation of heat-shock protein gene expression.
        Biochem. Res. Int. 2011; 2011 (21152185)
        • Tanaka M.
        • Miyahima A.
        Onconstatin M, a multifunctional cytokine.
        Rev. Physiol. Biochem. Pharmacol. 2003; 149 (12811586): 39-52
        • Terui K.
        • Enosawa S.
        • Haga S.
        • Zhang H.Q.
        • Kuroda H.
        • Kouchi K.
        • Matsunaga T.
        • Yoshida H.
        • Engelhardt J.F.
        • Irani K.
        • Ohnuma N.
        • Ozaki M.
        Stat3 confers resistance against hypoxia/reoxygenation-induced oxidative injury in hepatocytes through upregulation of Mn-SOD.
        J. Hepatol. 2004; 41 (15582129): 957-965
        • Touzet H.
        • Varré J.-S.
        Efficient and accurate P-value computation for Position Weight Matrices.
        Algorithms Mol. Biol. 2007; 2 (18072973): 15
        • Triboulet R.
        • Chang H.M.
        • Lapierre R.J.
        • Gregory R.I.
        Post-transcriptional control of DGCR8 expression by the Microprocessor.
        RNA. 2009; 15 (19383765): 1005-1011
        • Vandenplas J.
        • Calus M.P.L.
        • Sevillano C.A.
        • Windig J.J.
        • Bastiaansen J.W.M.
        Assigning breed origin to alleles in crossbred animals.
        Genet. Sel. Evol. 2016; 48 (27549177): 61
        • West J.W.
        Effects of heat-stress on production in dairy cattle.
        J. Dairy Sci. 2003; 86 (12836950): 2131-2144
        • Yenari M.A.
        • Liu J.
        • Zheng Z.
        • Vexler Z.S.
        • Lee J.E.
        • Giffard R.G.
        Antiapoptotic and anti-inflammatory mechanisms of heat-shock protein protection.
        Ann. N. Y. Acad. Sci. 2005; 1053 (16179510): 74-83
        • Yoshioka J.
        Thioredoxin reductase 2 (Txnrd2) regulates mitochondrial integrity in the progression of age-related heart failure.
        J. Am. Heart Assoc. 2015; 4 (26199229): 1-3
        • 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, Bos taurus.
        Genome Biol. 2009; 10 (19393038)