Skip to main content
Download PDF

Weighted single-step GWAS and RNA sequencing reveals key candidate genes associated with physiological indicators of heat stress in Holstein cattle

Journal of Animal Science and Biotechnology volume 13, Article number: 108 (2022) Cite this article

Abstract

Background

The study of molecular processes regulating heat stress response in dairy cattle is paramount for developing mitigation strategies to improve heat tolerance and animal welfare. Therefore, we aimed to identify quantitative trait loci (QTL) regions associated with three physiological indicators of heat stress response in Holstein cattle, including rectal temperature (RT), respiration rate score (RS), and drooling score (DS). We estimated genetic parameters for all three traits. Subsequently, a weighted single-step genome-wide association study (WssGWAS) was performed based on 3200 genotypes, 151,486 phenotypic records, and 38,101 animals in the pedigree file. The candidate genes located within the identified QTL regions were further investigated through RNA sequencing (RNA-seq) analyses of blood samples for four cows collected in April (non-heat stress group) and four cows collected in July (heat stress group).

Results

The heritability estimates for RT, RS, and DS were 0.06, 0.04, and 0.03, respectively. Fourteen, 19, and 20 genomic regions explained 2.94%, 3.74%, and 4.01% of the total additive genetic variance of RT, RS, and DS, respectively. Most of these genomic regions are located in the Bos taurus autosome (BTA) BTA3, BTA6, BTA8, BTA12, BTA14, BTA21, and BTA24. No genomic regions overlapped between the three indicators of heat stress, indicating the polygenic nature of heat tolerance and the complementary mechanisms involved in heat stress response. For the RNA-seq analyses, 2627 genes were significantly upregulated and 369 downregulated in the heat stress group in comparison to the control group. When integrating the WssGWAS, RNA-seq results, and existing literature, the key candidate genes associated with physiological indicators of heat stress in Holstein cattle are: PMAIP1, SBK1, TMEM33, GATB, CHORDC1, RTN4IP1, and BTBD7.

Conclusions

Physiological indicators of heat stress are heritable and can be improved through direct selection. Fifty-three QTL regions associated with heat stress indicators confirm the polygenic nature and complex genetic determinism of heat tolerance in dairy cattle. The identified candidate genes will contribute for optimizing genomic evaluation models by assigning higher weights to genetic markers located in these regions as well as to the design of SNP panels containing polymorphisms located within these candidate genes.

Graphical Abstract

Introduction

A global warming trend has been observed over the last century [1], which has unfavorable effects on livestock production. Thus, heat stress is becoming a major challenge facing the global dairy industry [2]. Heat stress is the result of the sum of external forces acting on an animal that evokes a series of behavioral and physiological responses, including sweating, higher respiration rate, vasodilation with increased blood flow to skin surface, reduced metabolic rate, decreased dry matter intake, and altered water metabolism [3]. During periods of heat stress the hypothalamic-pituitary-adrenal (HPA) and the sympathetic-adrenal-medullary (SAM) axis are activated for maintaining homeostasis in response to stressful stimuli [4]. Under chronic stress, cortisol secretion associated with immune suppression [5, 6] leading to the animal becoming more susceptible to disease and immune challenges [7]. Heat stress is shown to severely alter the welfare, as well as productive and reproductive performance of dairy cows [8]. Due to the thermal hysteresis existing [9] for milk yield in dairy cattle, physiological performance traits are better indicators for monitoring heat load of individuals. We have established two scoring systems for respiration rate and salivation (drooling) in addition to measuring rectal temperature in Holstein population [10]. Respiration rate score (RR), rectal temperature (RT), and drooling score (DS) are lowly heritable traits controlled by numerous quantitative trait loci (QTL), each with a small effect, thereby making it more difficult to obtain fast genetic progress for these traits.

Genome-wide association studies (GWAS) are used to discover genomic regions associated with phenotypes of interest [11], and consequently, biological processes in which they are involved in. For instance, Luo et al. [12] identified candidate genes related to RT by using single-SNP regression GWAS. Subsequently, gene expression analyses were performed to validate the key functional genes identified [12]. A recent method known as single-step GWAS [13] is becoming the gold standard method for GWAS as it enables the simultaneous integration of genomic, phenotypic, and pedigree information in a single analysis. ssGWAS has been widely used for scanning QTLs associated with complex traits in dairy cattle such as milk yield and composition traits in Holstein cattle [14, 15]; fertility and reproductive disorders [16]; and heat tolerance [17]. Furthermore, the combination of GWAS and transcriptomics has been shown to significantly enhance the understanding of the genetic architecture of complex traits [18]. Understanding how genetic variation shapes the phenotypic variability of complex traits requires the identification of such genetic polymorphisms through GWAS combined with functional evaluation of the key candidate genes through RNA sequencing (RNA-seq). RNA-seq is a technique that evaluates the quantity and sequences of RNA in a sample based on next generation sequencing (NGS) tools. RNA-seq analyzes the full transcriptome, indicating which of the genes encoded in the DNA are turned on or off and to what extent [19]. To our best knowledge, no study has combined weighted single-step GWAS (WssGWAS) and RNA-seq for the determination of candidate genes and genetic variants associated with physiological indicators of heat stress in dairy cattle. Thus, the main objectives of this study were to: 1) identify genomic regions associated with RT, RR, and DS in Holstein cattle based on WssGWAS; 2) search for candidate genes in the QTL regions that explain greater proportions of the total additive genetic variance of each trait; and, 3) utilize RNA-seq and gene network analyses to investigate the biological processes shared by the candidate genes identified for the three indicators of heat tolerance.

Material and methods

Phenotype, pedigree and genotypic data

A total of 69,837 (RT), 40,760 (RR), and 40,889 (DS) phenotypic records were collected in 15,303 lactating Holstein cows from 2013 to 2020. The process of data collection, scoring protocol of RR and DS, and farm management have been described in details by Luo et al. [10]. The descriptive statistics for the phenotypic data and environmental index (temperature-humidity index, THI) are presented in Table 1. Each lactating cow was recorded twice a day for two consecutive days (07:00–11:00 and 14:00–18:00). A small proportion of cows were measured in more than one lactation (i.e., > 4 records/cow).

Table 1 Descriptive statistics of phenotypic records and environmental variables used in the analyses
Full size table

The pedigree file contained 38,101 animals, spanning over three generations. A total of 3200 animals (3119 cows and 81 bulls) were genotyped using the Illumina 150 K Bovine Bead chip (Illumina, Inc., San Diego, CA, USA), which contains 139,377 single nucleotide polymorphism (SNPs). The individuals with call rate higher than 0.9 were kept and the call rate of genotyped individuals is 0.977 before imputation. All of the genotyped animals were included in a genotype imputation analyses for imputing the missing SNPs. This was done using the Beagle5.1 software [20]. Genotype quality control kept SNPs with: minor allele frequency (MAF) greater than 0.05, no extreme departure from Hardy-Weinberg equilibrium (P-value greater than 10−6), known chromosome and genome position, and located in the autosomal chromosomes. After the quality control, 114,766 SNPs and 3200 animals were kept for further analyses.

Statistical analyses

A multiple-trait model was employed to estimate variance and covariance components, which can be described as follows:

$$\boldsymbol{y}=\boldsymbol{Xb}+\boldsymbol{Za}+\boldsymbol{Wpe}+\boldsymbol{e}$$

where y is the vector of phenotypic records (RT, RR, and DS); b is the vector of systematic effects including farm-year (for RT) or farm-year-scoring person (for TT and DS), parity (1, 2, or 3+), lactation stage (days in milk 1–50, 51–100, 101–150, 151–200, 201–250, 251–300, or > 300), milking status (data collected before milking, after milking, or unknown), and THI (as a continuous covariable in the model fitting a linear regression); a is the vector of random additive genetic effects; p is the vector of random permanent environmental effects; e is the vector of random residuals; and X, Z, and W are the incidence matrices linking phenotypic records to b, a, and pe, respectively. The model assumptions are:

$$\boldsymbol{a}\sim N\left(\mathbf{0},\boldsymbol{H}\otimes \left[\begin{array}{c}{\sigma}_{a(RT)}^2\ {\sigma}_{a\left( RT, RR\right)}\ {\sigma}_{a\left( RT, DS\right)}\\ {}{\sigma}_{a\left( RT, RR\right)}\ {\sigma}_{a(RR)}^2\ {\sigma}_{a\left( RR, DS\right)}\\ {}{\sigma}_{a\left( RT, DS\right)}\ {\sigma}_{a\left( RR, DS\right)}\ {\sigma}_{a(DS)}^2\ \end{array}\right]\right)$$
$$\boldsymbol{pe}\sim N\left(\mathbf{0},\boldsymbol{I}\otimes \left[\begin{array}{c}{\sigma}_{pe(RT)}^2\ {\sigma}_{pe\left( RT, RR\right)}\ {\sigma}_{pe\left( RT, DS\right)}\\ {}{\sigma}_{pe\left( RT, RR\right)}\ {\sigma}_{pe(RR)}^2\ {\sigma}_{pe\left( RR, DS\right)}\\ {}{\sigma}_{pe\left( RT, DS\right)}\ {\sigma}_{pe\left( RR, DS\right)}\ {\sigma}_{pe(DS)}^2\ \end{array}\right]\right)$$
$$\boldsymbol{e}\sim N\left(\mathbf{0},\boldsymbol{I}\otimes \left[\begin{array}{c}{\sigma}_{e(RT)}^2\ {\sigma}_{e\left( RT, RR\right)}\ {\sigma}_{e\left( RT, DS\right)}\\ {}{\sigma}_{e\left( RT, RR\right)}\ {\sigma}_{e(RR)}^2\ {\sigma}_{e\left( RR, DS\right)}\\ {}{\sigma}_{e\left( RT, DS\right)}\ {\sigma}_{e\left( RR, DS\right)}\ {\sigma}_{e(DS)}^2\ \end{array}\right]\right)$$

where \({\sigma}_{k\left(i,j\right)}^2\) is the variance of the effect k (a, pe, e) of trait i, σk(i, j) is the covariance between trait i and j for the effect k. H is the matrix that combines pedigree and genomic information [21], and I is an identity matrix. The inverse of H was calculated as [21]:

$${\boldsymbol{H}}^{-\mathbf{1}}={\boldsymbol{A}}^{-\mathbf{1}}+\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{G}}^{-\mathbf{1}}-{\boldsymbol{A}}_{\mathbf{22}}^{-\mathbf{1}}\end{array}\right]$$

where A is the numerator relationship matrix based on pedigree for all animals; A22 is the numerator relationship matrix for genotyped animals. The G matrix was computed as [22]:

$$\boldsymbol{G}=\frac{\boldsymbol{Z}\boldsymbol{D}{\boldsymbol{Z}}^{\prime }}{\sum_{i=1}^M2{p}_i\left(1-{p}_i\right)}$$

where Z is a matrix of gene content adjusted for allele frequencies (0, 1, or 2 for aa, Aa, and AA, respectively); D is a diagonal matrix of weights for SNP variances (initially D = I); M is the number of SNPs, and pi is the minor allele frequency of the ith SNP. Variance components were estimated using the average information-restricted maximum likelihood (AI-REML) procedure implemented in the AIREMLF90 package from the BLUPF90 family programs [23].

Weighted single-step genome-wide association study (WssGWAS)

The estimates of SNP effects and weights for the WssGWAS analyses (three iterations) were obtained according to Wang et al. [24]. The weight for each SNP was calculated as: \({d}_i={1.125}^{\frac{\left|{\hat{a}}_i\right|}{sd\left({\hat{a}}_i\right)}-2}\)[22], where i is the ith SNP. The percentage of the total additive genetic variance explained by the ith region was calculated as:

$$\frac{\mathrm{Var}\left(a_i\right)}{\upsigma_a^2}\times 100\%=\frac{\mathrm{Var}\left({\sum}_{j=1}^{10}{\mathrm{Z}}_j{\hat{\mathrm{u}}}_j\right)}{\upsigma_a^2}\times 100\%$$

Where ai is genetic value of the ith region that consists of contiguous 10 SNPs, \({\sigma}_a^2\) is the total additive genetic variance, Zj is a vector of gene content of the jth SNP for all individuals, and \({\hat{u}}_j\) is the marker effect of the jth SNP within the ith region.

Candidate genes detection and functional enrichment analyses

Genomic windows of 10 consecutive SNPs that explained 0.15% or more of the total additive genetic variance, based on the WssGWAS analyses, were considered to be associated with the studied traits. A Manhattan plot was created using the R software [25]. Genes were annotated on the basis of the starting and ending coordinates of each window (using the ARS-UCD1.2 assembly as the reference genome; GCA_002263795.2) by using the R package ‘Biomart’ of Ensembl [26]. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms were enriched via the “clusterProfiler” package [27].

RNA sequencing

Eight primiparous cows with DIM ranging from 135 to 144 d (mid-lactation and pregnant) were selected from 3200 genotyped individuals for RNA-seq. Blood samples were collected via the coccygeal venipuncture in the Spring (4 cows in April during the thermoneutral period with average daily temperature lower than 25 °C) and Summer (4 cows in July during heat stress period with average daily temperature higher than 25 °C) seasons, following the RT, TT, and DS recording on the same day. RT, TT, and DS were collected after blood samples collected to reduce artificial stimulation. The experimental design in current study is consisted with other previous research investigating heat stress by RNA-seq [28, 29]. The RNA was isolated from blood according to the manufacturer’s instructions of the TRIzol Reagent method [30]. The RNA concentration and quality were assessed using the Equalbit RNA BR Assay Kit (Invitrogen, CA, USA) and Nanodrop 2000 (Thermo, Massachusetts, USA). RNA integrity was detected by 1% agarose gel electrophoresis and used for library construction with 28S/18S > 1. For RNA-seq library, 2 μg of RNA was used for purification and fragment using NEBNext Poly(A) mRNA Magnetic Isolation Module (Cat No. E7490S, New England Biolabs Ltd., Hitchin, Herts, UK) then followed by cDNA library with NEBNext Ultra RNA Library Prep Kit for Illumina (Cat No. E7530S, New England Biolabs Ltd., Hitchin, Herts, UK). All libraries were quantified by Equalbit DNA BR Assay Kit (Invitrogen, CA, USA) and pooled equimolarly. They were finally submitted for sequencing using the NovaSeq 6000 System (Illumina, Inc., San Diego, CA, USA) which generated 150 bases paired-end reads.

Differential expression and functional analyses

The quality of the sequencing reads was evaluated using the FastQC software (v0.11.9) and global trimming using the Fastp [31]. All clean reads were mapped to the bovine genome of version ARS-UCD1.2 using the software STAR [32], and a Picard query [33] was carried out to eliminate duplicates. We investigated population structure through a principal component analysis (PCA) implemented in the PLINK software [34]. For differential expression gene screening, exact test based on quantile-adjusted conditional maximum likelihood (qCML) was performed using the edgeR [35] R package with criteria fold change ≥ 2 and 0.05 for the alpha of false discovery rate (4 samples vs. 4 samples). The heatmap was constructed using the pheatmap R package [36].

Results

Genetic parameter estimates

Estimates of variance components for RT, RR, and DS are shown in Table 2. The heritability estimates ranged from 0.03 (DS) to 0.06 (RT), with repeatability estimates ranging from 0.12 (DS) to 0.19 (RT), indicating the large environmental influence in these indicators of heat stress. Statistically significant genetic correlations were observed for RT with RR (0.22) and RR with DS (0.21).

Table 2 Genetic parameters, heritability, and repeatability for the evaluated physiological traits under heat stress
Full size table

QTL mapping and trait-related gene identification

Figure 1 presents the genetic variance of each genomic window after performing WssGWAS with three iterations. Fourteen, 19, and 20 genomic regions reached the pre-defined threshold (0.15%) and they explained 2.94%, 3.74%, and 4.01% of the total additive genetic variance for RT, RR, and DS, respectively. Most of these genomic regions are located in the Bos taurus autosome (BTA) BTA3, BTA6, BTA8, BTA12, BTA14, BTA21, and BTA24. However, we found no overlapping genomic regions between the three physiological traits in this study. The detailed information about these genomic regions related to the three physiological traits and candidate genes are presented in Table 3.

Fig. 1
figure 1

Proportion of the total additive genetic variance of 10-SNP genomic windows based on the weighted single-step genome association studies. A The Manhattan plot for rectal temperature. B The Manhattan plot for respiration rate score. C The Manhattan plot for drooling score. The red dashed lines represent the threshold 0.15% of the total additive genetic variance

Full size image
Table 3 Information of 10-SNP windows explaining more than 0.15% of the total additive genetic variance for rectal temperature, respiration rate score, and drooling score in Holstein cattle
Full size table

Many windows with small effects regulated RT, RR, and DS co-occurred, indicating that these three physiological indicators of heat stress are highly polygenic traits. A total of 54 protein-coding genes (27, 14, and 13 protein-coding genes located by the genomic region associated with RT, RS, and DS, respectively) were annotated in these genomic regions according to the Ensembl database. The GO enrichment and KEGG pathway analyses of the 54 protein-coding genes for the three physiological traits revealed 68 significant GO terms (Additional file 1: Fig. S1) including biological process as well as molecular functions and two significant KEGG pathway (purine metabolism and thiamine metabolism).

General features and genetic background of cows sampled for RNA sequencing analysis

The environmental THI and physiological performance (RT, RR, and DS) parameters during the two sample collection periods (April is the thermoneutral season in Beijing, named non-heat stress group, while July is the heat stress season, named heat stress group) were significantly different. The milk yield of the cows in the non-heat stress group (41.48 ± 3.52) was significantly higher than the cows in the heat stress group (35.25 ± 3.71). The description of environmental THI, physiological performance, and milk yield during periods of blood samples collecting are shown in Table 4. The genetic relationship of the dairy population evaluated is presented in Fig. 2, which shows a diverse genetic background of the cows (red and blue dots) sampled for the RNA-seq analysis.

Table 4 Summary statistics for environmental THI, physiological performance, and milk yield in April (non-heat stress) and July (heat stress)
Full size table
Fig. 2
figure 2

SNP-based principal component analysis of genotyped animals showing the diversity of animals selected for RNA sequencing samples. Red points represent the four heat stress season (HS) individuals, blue points represent the four non-heat stress season (NHS) individuals, and grey points represent the other individuals used in the weighted single-step genome-wide association study

Full size image

Identification of differentially expressed mRNAs by sequencing

There is a significant difference between the heat stress (HS) group and non-heat stress (NHS) group for RNA-seq counts based on principle component analysis and clustering structure (Fig. 3). The analysis of the differently expressed genes (DEGs) were detected based on the quantile-adjusted conditional maximum likelihood method. A total of 2627 significantly downregulated genes were found and 369 upregulated genes in the NHS group compared to HS group (Additional file 2: Table S1). The GO enrichment and KEGG pathway analysis of DEGs revealed 21 significant GO terms including Biological Process, Cellular Component, and Molecular Function and 86 significant KEGG pathways (Additional file 3: Fig. S2). There were 14 genes identified based on the WssGWAS that were also DEGs. Table 5 shows that two of the significant common genes are upregulated and the other ones are downregulated. Figure 4 shows the heatmap of the mRNA expression of these 14 genes using hierarchical cluster of log2 from the relative normalized expression in bovine blood transcriptome.

Fig. 3
figure 3

Cluster of heat stress and non-heat stress on the basis of read counts and volcano plot displaying differentially expressed genes. A, B Principle component analysis and dendrogram for samples collected under heat stress and non-heat stress on the basis of read counts. Red points in dot plot and red lines in dendrogram represent heat stress samples. Blue points in dot plot and blue lines in dendrogram represent heat stress samples. C Volcano plot showing significantly expressed genes. The red and green dots denote significantly downregulated and upregulated genes, respectively

Full size image
Table 5 List of overlapping genes between weighted single-step genome-wide association analysis and RNA sequencing
Full size table
Fig. 4
figure 4

Heat map showing the relative normalized expression of the mRNA of 14 genes in eight individuals between period of April-non-heat stress and July-heat stress using hierarchical cluster. Genes with increased expression are shown in red while those with decreased expression are shown in blue

Full size image

Discussion

The estimation of genetic parameters and detection of genomic regions for three physiological indicators of heat stress in Holstein cattle contributes to further understanding of the genetic architecture of heat tolerance in cattle. The heritability estimates for the three traits are low (0.03–0.06) in the population evaluated but statistically different than zero. The variance component results utilizing the H matrix are similar to earlier estimates for the same traits using the numerator relationship matrix based on pedigree (A matrix; [10], indicating that the size of the datasets used is large enough to accurately estimate genetic parameters for these three traits. The heritability estimated for RT in 1695 lactating Holstein cows during heat stress (0.17 ± 0.13; [37]) and in 3396 straight-bred and crossbred Romosinuano, Brahman, and Angus cattle (0.19 ± 0.03; [38] using bivariate models with coat score are higher than the estimates observed in our population. However, their accuracies of estimation are lower due to their smaller datasets.

For low-heritability traits, substantial genetic progress can be achieved by increasing intensity and accuracy of selection [39, 40]. Subsequently, genomic selection could significantly speed up the rates of genetic progress [41] by improving accuracy of selection and reducing generation interval. In our earlier study [12], single-SNP regression GWAS was performed for RT in a subset of the current population including 7598 Chinese Holstein cattle with 1114 genotyped cows. Ten SNPs (located on BTA3, BTA4, BTA8, BTA13, BTA14, and BTA29) were found to be significantly associated with RT and five positional candidate genes were identified (SPAG17, FAM107B, TSNARE1, RALYL, and PHRF1). In the present study we identified 53 trait-related genomic regions and 53 protein-coding genes through WssGWAS. Contrary to the single-SNP GWAS, more genomic regions were captured via WssGWAS (e.g., [42, 43], showing that WssGWAS can be more successful in detecting QTL. According to study of simulation [42], single-step GWAS (based on single-step genomic best linear unbiased prediction) could take care of correction for population structure in GWAS like other mixed linear models (e.g. EMMAX). Comparing to the method of WssGWAS, BayesB appears to overly shrink regions to zero, while overestimating the amount of genetic variation attributed to the remaining SNP effects in chicken population [24]. We compared two scenarios (weighting step including re-estimated GEBV and weighting step only with re-calculated SNP effect [24]) along with the inclusion of times of iteration for weighting (up to 5 iterations, results not shown). Although, the process of WssGWAS with only re-calculated SNP effects or more than three iterations resulted in too much noise (results not shown), which is consistent with prior studies based on simulated data in dairy cattle [44] and real data in broiler chickens [24]. Additionally, it is critical to select the number of SNPs included in the sliding windows when performing (W)ssGWAS as wider windows could capture more genomic regions [45, 46]. The SNP chip used has an average inter-marker spacing of 26.07 kb, and the average distance between adjacent SNPs with linkage disequilibrium (r2) higher than 0.4 was 200 kb in the studied population. Therefore, we selected 10-SNP windows in linkage disequilibrium when scanning QTL related to the three physiological traits.

RNA-seq analysis contributes to annotating new genes and splice variants, and provides cell- and context-specific quantification of gene expression. Research over the last few years has identified some of the physiological, metabolic, cellular, and molecular responses to heat stress in cattle [47, 48]. However, the level of gene expression depends on the physiological state of the organism and tissue analyzed. Cellular and transcriptomic adaptation of bovine granulosa cells were characterized to different heat stress intensities (39 °C, 40 °C, and 41 °C) in-vitro [49]. Several heat-responsive genes from different functional classes were identified and their associated pathways related to heat stress chaperons, cell death, apoptosis, hormonal synthesis, and oxidative stress. The expression of miRNAs in dairy cattle mammary gland under heat stress was investigated [50] (483 known bovine miRNAs and 139 novel miRNAs were identified), which detected the heat-dependent differential modulation of miRNAs. Decreased expression of heat response-associated genes in blood during HS was also observed very recently in lactating cows [28]. However, it was not considered that if metabolism of the milking stage from Spring to Summer directly affects the gene expression levels. In the current study of RNA-seq, we selected primiparous cows with similar DIM to avoid the effect of different lactation stage and the identified DEGs involved in the pathway of apoptosis, cellular senescence and autophagy, known to be affected by heat stress.

Validation of the results of GWAS is an essential component of the experimental program. It could increase the level of trust held by animal breeders in genetic improvement, and is important as confirmation of a scientific hypotheses [51]. There are many methods used to verify identified candidate genes/QTLs from GWAS results, e.g., SNP chip analysis [52], qRT-PCR [12], or association analysis in another population [53]. RNA-seq with the capabilities of high-throughput sequencing, could help us better understand the translation of genetic loci into biological mechanisms that underlie phenotypic expression of important traits [54]. Fourteen genes identified in WssGWAS were significantly and differently expressed between cows in heat stress period and thermoneutral conditions. PMAIP1 is a proapoptotic member of the BCL-2 protein family that acts as a proapoptotic sensitizer/de-repressor and regulates diverse cellular functions in autophagic cell death and metabolism [55, 56]. When comparing to heat shock of 38 °C and 43 °C, the expression of PMAIP1 was significantly and differentially upregulated in control (32 °C) mouse spermatogenic cells [57]. Another upregulated gene during heat stress based on the differential gene analysis is SBK1, which has been reported as part of widespread expression pattern involved in the protection of cells from apoptosis induced by the viral infection in ovarian cancer cells and promoting the cellular survival [58].

Five genes (TMEM33, GATB, CHORDC1, RTN4IP1, and BTBD7), out of the 12 significantly downregulated genes, were identified in prior studies involving heat stress/heat shock and cellular adaptive functions in the presence of stressors or disrupting molecular mechanisms. For instance, the overexpression of TMEM33 cause expression of endoplasmic reticulum stress-induced cell death signals increase, which leads to activation of the unfolded protein response signaling cascade and induction of an apoptotic cell death [59]. The gene GATB was identified to show significant transcriptional differences (downregulated gene) in response to heat shock where a temperature shift from 37 °C to 42 °C was accessed through the use of a microarray [60]. The GATB gene downregulation is shown to be responsible for the respiratory chain enzyme deficiencies [61], and thus disturbing the mitochondrial function and possesses complications for the much-needed energy intensive physiological mechanisms to dissipate incremental heat load in cattle. Furthermore, during the temperature increase phase, GATB composed an operon and was downregulated in Clostridium botulinum (temperature shift from 37 °C to 45 °C over a period of 15 min; [62]. GO annotations related to CHORDC1 include Hsp90 protein binding, and this gene was associated with the regulation of heat stress and immune response processes in rats [63]. Additionally, upregulation of CHORDC1 is shown to be primarily associated with dystrophic conditions [64]. Therefore, it can be assumed that the downregulation of CHORDC1 is a positive adaptive sign of mechanisms involving thermotolerance. RTN4IP1 (Reticulon 4 Interacting Protein 1) is included in the gene ontology of oxidoreductase activity and transferase activity. The proteins of RTN4IP1 were down-regulated at 26 °C compared to 18 °C in wild-type zebrafish [65]. The expression of RTN4IP1 tend to decrease in the presence of stressors and decreased HSPs response, which again points towards a sort of possible homeostatic cellular response mechanisms initiated by the heat stressed dairy cow. BTBD7 [BTB (POZ) domain-containing 7] is involved in a variety of biological functions [66]. BTBD7 has been shown to be induced by a matrix protein at sites of cleft progression and induce a transcription factor and suppress cell adhesion [67]. Other studies described BTBD7 as a cell growth suppressor protein (ZNF238 is expressed in postmitotic brain cells and inhibits brain tumor growth) and a promotor of angiogenesis [68]. This function points out to the possible conservation mechanisms of heat stressed cows showing high energetic mechanisms directed at attaining thermo-neutrality and mechanisms of vasodilation of peripheral blood distribution network [69]. Query of literature strongly support the involvement of the possible candidate genes proposed in this study, in the regulatory mechanisms directed at differential homoerotic and homeostatic mechanisms manifesting in various physiological modifications towards heat stress, as described in previous studies [70]. Given the relevant support from prior research studies about the possible involvement in heat stress/heat shock response, the identification of genomic regions and candidate genes through WssGWAS, and biological validation through RNA-seq analyses indicate seven genes (PMAIP1, SBK1, TMEM33, GATB, CHORDC1, RTN4IP1, BTBD7) as likely playing an important role in heat stress response. Further studies should investigate these candidate genes in more depth and in more controlled on-farm or in-vitro experiments in response to more divergent environmental heat loads.

Conclusions

Physiological indicators of heat tolerance are heritable and can be improved through direct selection. We identified 54 candidate genes associated with rectal temperature, respiration rate score, and drooling score in Chinese Holstein cattle using WssGWAS. A validation experiment based on RNA-seq of blood samples in Holstein cattle provided evidence of their possible role in the physiological responses to heat stress. The identified candidate genes (PMAIP1, SBK1, TMEM33, GATB, CHORDC1, RTN4IP1, BTBD7) may provide knowledge for developing genomic evaluation models by assigning higher weights to genetic markers located in these regions or the development of SNP panels contained polymorphisms in these genes.

Availability of data and materials

The data sets used during the current study are available from the corresponding author on reasonable request.

Abbreviations

RT:

Rectal temperature

RR:

Respiration rate score

DS:

Drooling score

QTL:

Quantitative trait loci

GWAS:

Genome-wide association studies

RNA-seq:

RNA sequencing

NGS:

Next generation sequencing

WssGWAS:

Weighted single-step GWAS

THI:

Temperature-humidity index

DIM:

Days in milk

AI-REML:

Average information-restricted maximum likelihood

KEGG:

Kyoto Encyclopedia of Genes and Genomes

GO:

Gene Ontology

DEGs:

Differently expressed genes

HS:

Heat stress

NHS:

Non-heat stress season

References

  1. Schär C, Vidale PL, Lüthi D, Frei C, Häberli C, Liniger MA, et al. The role of increasing temperature variability in European summer heatwaves. Nature. 2004;427(6972):332–6. https://doi.org/10.1038/nature02300.

    Article  CAS  PubMed  Google Scholar 

  2. Das R, Sailo L, Verma N, Bharti P, Saikia J, Imtiwati KR. Impact of heat stress on health and performance of dairy animals: a review. Vet World. 2016;9(3):260–8. https://doi.org/10.14202/vetworld.2016.260-268.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Polsky L, von Keyserlingk M. Invited review: effects of heat stress on dairy cattle welfare. J Dairy Sci. 2017;100(11):8645–57. https://doi.org/10.3168/jds.2017-12651.

    Article  CAS  PubMed  Google Scholar 

  4. Sejian V, Bhatta R, Gaughan JB, Dunshea FR, Lacetera N. Review: adaptation of animals to heat stress. Animal. 2018;12(s2):s431–44. https://doi.org/10.1017/S1751731118001945.

    Article  CAS  PubMed  Google Scholar 

  5. Ju XH, Xu HJ, Yong YH, An LL, Jiao PR, Liao M. Heat stress upregulation of toll-like receptors 2/4 and acute inflammatory cytokines in peripheral blood mononuclear cell (PBMC) of Bama miniature pigs: An in vivo and in vitro study. Animal. 2014;8(9):1462–8. https://doi.org/10.1017/S1751731114001268.

    Article  CAS  PubMed  Google Scholar 

  6. Jin Y, Hu Y, Han D, Wang M. Chronic heat stress weakened the innate immunity and increased the virulence of highly pathogenic avian influenza virus H5N1 in mice. J Biomed Biotechnol. 2011;2011:367846. https://doi.org/10.1155/2011/367846.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Bagath M, Krishnan G, Devaraj C, Rashamol VP, Pragna P, Lees AM, et al. The impact of heat stress on the immune system in dairy cattle: a review. Res Vet Sci. 2019;126:94–102. https://doi.org/10.1016/j.rvsc.2019.08.011.

    Article  CAS  PubMed  Google Scholar 

  8. Sammad A, Wang YJ, Umer S, Lirong H, Khan I, Khan A, et al. Nutritional physiology and biochemistry of dairy cattle under the influence of heat stress: consequences and opportunities. Animals. 2020;10(5):793. https://doi.org/10.3390/ani10050793.

    Article  Google Scholar 

  9. Parkhurst AM. Model for understanding thermal hysteresis during heat stress: a matter of direction. Int J Biometeorol. 2010;54(6):637–45. https://doi.org/10.1007/s00484-009-0299-z.

    Article  CAS  PubMed  Google Scholar 

  10. Luo H, Brito LF, Li X, Su G, Dou J, Xu W, et al. Genetic parameters for rectal temperature, respiration rate, and drooling score in Holstein cattle and their relationships with various fertility, production, body conformation, and health traits. J Dairy Sci. 2021;104(4):4390–403. https://doi.org/10.3168/jds.2020-19192.

    Article  CAS  PubMed  Google Scholar 

  11. Chen Z, Brito LF, Luo H, Shi R, Chang Y, Liu L, et al. Genetic and genomic analyses of service sire effect on female reproductive traits in Holstein cattle. Front Genet. 2021;12:713575. https://doi.org/10.3389/fgene.2021.713575.

  12. Luo H, Li X, Hu L, Xu W, Chu Q, Liu A, et al. Genomic analyses and biological validation of candidate genes for rectal temperature as an indicator of heat stress in Holstein cattle. J Dairy Sci. 2021;104(4):4441–51. https://doi.org/10.3168/jds.2020-18725.

    Article  CAS  PubMed  Google Scholar 

  13. Wang H, Misztal I, Aguilar I, Legarra A, Muir W. Genome-wide association mapping including phenotypes from relatives without genotypes. Genet Res. 2012;94:73–83. https://doi.org/10.1017/S0016672312000274.

    Article  CAS  Google Scholar 

  14. Freitas PHF, Oliveira HR, Silva FF, Fleming A, Miglior F, Schenkel FS, et al. Genomic analyses for predicted Milk fatty acid composition throughout lactation in north American Holstein cattle. J Dairy Sci. 2020;103(7):6318–31. https://doi.org/10.3168/jds.2019-17628.

    Article  CAS  PubMed  Google Scholar 

  15. Oliveira HR, Lourenco DAL, Masuda Y, Misztal I, Tsuruta S, Jamrozik J, et al. Single-step genome-wide association for longitudinal traits of Canadian Ayrshire, Holstein, and Jersey dairy cattle. J Dairy Sci. 2019;102(11):9995–10011. https://doi.org/10.3168/jds.2019-16821.

    Article  CAS  PubMed  Google Scholar 

  16. Guarini AR, Lourenco DAL, Brito LF, Sargolzaei M, Baes CF, Miglior F, et al. Genetics and genomics of reproductive disorders in Canadian Holstein cattle. J Dairy Sci. 2019;102(2):1341–53. https://doi.org/10.3168/jds.2018-15038.

    Article  CAS  PubMed  Google Scholar 

  17. Shi R, Brito LF, Liu A, Luo H, Chen Z, Liu L, et al. Genotype-by-environment interaction in Holstein heifer fertility traits using single-step genomic reaction norm models. BMC Genomics. 2021;22:193. https://doi.org/10.1186/s12864-021-07496-3.

    Article  CAS  Google Scholar 

  18. Zhao B, Luo H, Huang X, Wei C, Di J, Tian Y, et al. Integration of a single-step genome-wide association study with a multi-tissue transcriptome analysis provides novel insights into the genetic basis of wool and weight traits in sheep. Genet Sel Evol. 2021;53:56. https://doi.org/10.1186/s12711-021-00649-8.

  19. Ozsolak F, Milos PM. RNA sequencing: advances, challenges and opportunities. Nat Rev Genet. 2011;12(2):87–98. https://doi.org/10.1038/nrg2934.

    Article  CAS  PubMed  Google Scholar 

  20. Browning BL, Zhou Y, Browning SR. A one-penny imputed genome from next-generation reference panels. Am J Hum Genet. 2018;103(3):338–48. https://doi.org/10.1016/j.ajhg.2018.07.015.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Aguilar I, Misztal I, Johnson DL, Legarra A, Tsuruta S, Lawlor TJ. Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score. J Dairy Sci. 2010;93(2):743–52. https://doi.org/10.3168/jds.2009-2730.

    Article  CAS  PubMed  Google Scholar 

  22. VanRaden PM. Efficient methods to compute genomic predictions. J Dairy Sci. 2008;91(11):4414–23. https://doi.org/10.3168/jds.2007-0980.

    Article  CAS  PubMed  Google Scholar 

  23. Misztal I, Tsuruta S, Lourenco DAL, Masuda Y, Aguilar I, Legarra A, et al. Manual for BLUPF90 family programs: University of Georgia; 2018. http://nce.ads.uga.edu/wiki/doku.php?id=documentation.

    Google Scholar 

  24. Wang H, Misztal I, Aguilar I, Legarra A, Fernando RL, Vitezica Z, et al. Genome-wide association mapping including phenotypes from relatives without genotypes in a single-step (ssGWAS) for 6-week body weight in broiler chickens. Front Genet. 2014;5:134. https://doi.org/10.3389/fgene.2014.00134.

    PubMed  PubMed Central  Google Scholar 

  25. Team RC. R: A Language and environment for statistical computing. R Foundation for statistical computing, Vienna, Austria. 2021. URL https://www.r-project.org.

    Google Scholar 

  26. Durinck S, Spellman PT, Birney E, Huber W. Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nat Protoc. 2009;4(8):1184–91. https://doi.org/10.1038/nprot.2009.97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Yu G, Wang L, Han Y, He Q. ClusterProfiler: An R package for comparing biological themes among gene clusters. OMICS: J Integr Biol. 2012;16(5):284–7. https://doi.org/10.1089/omi.2011.0118.

    Article  CAS  Google Scholar 

  28. Liu G, Liao Y, Sun B, Guo Y, Deng M, Li Y, et al. Effects of chronic heat stress on mRNA and miRNA expressions in dairy cows. Gene. 2020;742:144550. https://doi.org/10.1016/j.gene.2020.144550.

    Article  CAS  PubMed  Google Scholar 

  29. Yue S, Wang Z, Wang L, Peng Q, Xue B. Transcriptome functional analysis of mammary gland of cows in heat stress and thermoneutral condition. Animals. 2020;10(6):1015.

  30. Rio DC, Ares MJ, Hannon GJ, Nilsen TW. Purification of RNA using TRIzol (TRI reagent). Cold Spring Harb Protoc. 2010;2010(6):t5439. https://doi.org/10.1101/pdb.prot5439.

    Article  Google Scholar 

  31. Chen S, Zhou Y, Chen Y, Gu J. Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90. https://doi.org/10.1093/bioinformatics/bty560.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21. https://doi.org/10.1093/bioinformatics/bts635.

    Article  CAS  PubMed  Google Scholar 

  33. Picard. URL http://broadinstitute.github.io/picard/. Accessed 9 Dec 2020.

  34. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81(3):559–75. https://doi.org/10.1086/519795.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012;40(10):4288–97. https://doi.org/10.1093/nar/gks042.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Raivo Kolde. Pheatmap: pretty heatmaps. 2019. URL https://CRAN.R-project.org/package=pheatmap.

    Google Scholar 

  37. Dikmen S, Cole JB, Null DJ, Hansen PJ. Heritability of rectal temperature and genetic correlations with production and reproduction traits in dairy cattle. J Dairy Sci. 2012;95(6):3401–5. https://doi.org/10.3168/jds.2011-4306.

    Article  CAS  PubMed  Google Scholar 

  38. Riley DG, Chase CC, Coleman SW, Olson TA. Genetic assessment of rectal temperature and coat score in Brahman, Angus, and Romosinuano crossbred and straightbred cows and calves under subtropical summer conditions. Livest Sci. 2012;148(1–2):109–18. https://doi.org/10.1016/j.livsci.2012.05.017.

    Article  Google Scholar 

  39. Muranty H, Troggio M, Sadok IB, Rifaï MA, Auwerkerken A, Banchi E, et al. Accuracy and responses of genomic selection on key traits in apple breeding. Hortic Res-England. 2015;2:15060. https://doi.org/10.1038/hortres.2015.60.

  40. Song H, Zhang J, Zhang Q, Ding X. Using different single-step strategies to improve the efficiency of genomic prediction on body measurement traits in pig. Front Genet. 2019;9:730. https://doi.org/10.3389/fgene.2018.00730.

  41. Wiggans GR, Cooper TA, VanRaden PM, Cole JB. Technical note: adjustment of traditional cow evaluations to improve accuracy of genomic predictions. J Dairy Sci. 2011;94(12):6188–93. https://doi.org/10.3168/jds.2011-4481.

    Article  CAS  PubMed  Google Scholar 

  42. Mancin E, Lourenco D, Bermann M, Mantovani R, Misztal I. Accounting for population structure and phenotypes from relatives in association mapping for farm animals: a simulation study. Front Genet. 2021;12:642065. https://doi.org/10.3389/fgene.2021.642065.

  43. Marques DBD, Bastiaansen JWM, Broekhuijse MLWJ, Lopes MS, Knol EF, Harlizius B, et al. Weighted single-step GWAS and gene network analysis reveal new candidate genes for semen traits in pigs. Genet Sel Evol. 2018;50:40. https://doi.org/10.1186/s12711-018-0412-z.

  44. Lourenco DAL, Fragomeni BO, Bradford HL, Menezes IR, Ferraz JBS, Aguilar I, et al. Implications of SNP weighting on single-step genomic predictions for different reference population sizes. J Anim Breed Genet. 2017;134(6):463–71. https://doi.org/10.1111/jbg.12288.

    Article  CAS  PubMed  Google Scholar 

  45. Cáceres P, Barría A, Christensen KA, Bassini LN, Correa K, Garcia B, et al. Genome-scale comparative analysis for host resistance against sea lice between Atlantic Salmon and Rainbow trout. Sci Rep-UK. 2021;11:13231. https://doi.org/10.1038/s41598-021-92425-3.

  46. Dikmen S, Cole JB, Null DJ, Hansen PJ. Genome-wide association mapping for identification of quantitative trait loci for rectal temperature during heat stress in Holstein cattle. PLoS One. 2013;8(7):e69202. https://doi.org/10.1371/journal.pone.0069202.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Srikanth K, Kwon A, Lee E, Chung H. Characterization of genes and pathways that respond to heat stress in Holstein calves through transcriptome analysis. Cell Stress Chaperones. 2017;22(1):29–42. https://doi.org/10.1007/s12192-016-0739-8.

    Article  CAS  PubMed  Google Scholar 

  48. Sengar GS, Deb R, Singh U, Junghare V, Hazra S, Raja TV, et al. Identification of differentially expressed microRNAs in Sahiwal (Bos Indicus) breed of cattle during thermal stress. Cell Stress Chaperones. 2018;23(5):1019–32. https://doi.org/10.1007/s12192-018-0911-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Khan A, Dou J, Wang Y, Jiang X, Khan MZ, Luo H, et al. Evaluation of heat stress effects on cellular and transcriptional adaptation of bovine granulosa cells. J Anim Sci Biotechno. 2020;11:25. https://doi.org/10.1186/s40104-019-0408-8.

  50. Li Q, Yang C, Du J, Zhang B, He Y, Hu Q, et al. Characterization of miRNA profiles in the mammary tissue of dairy cattle in response to heat stress. BMC Genomics. 2018;19:975. https://doi.org/10.1186/s12864-018-5298-1.

  51. Wen H, Luo H, Yang M, Augustino SMA, Wang D, Mi S, et al. Genetic parameters and weighted single-step genome-wide association study for supernumerary teats in Holstein cattle. J Dairy Sci. 2021;104(11):11867–77. https://doi.org/10.3168/jds.2020-19943.

    Article  CAS  PubMed  Google Scholar 

  52. Zhao BR, Fu XF, Tian KC, Huang XX, Jiang D, Yan B, et al. Identification of SNPs and expression patterns of FZD3 gene and its effect on wool traits in Chinese merino sheep (Xinjiang type). J Integr Agric. 2019;18(10):2351–60. https://doi.org/10.1016/S2095-3119(19)62735-8.

    Article  CAS  Google Scholar 

  53. Dikmen S, Wang XZ, Ortega MS, Cole JB, Null DJ, Hansen PJ. Single nucleotide polymorphisms associated with thermoregulation in lactating dairy cows exposed to heat stress. J Anim Breed Genet. 2015;132(6):409–19. https://doi.org/10.1111/jbg.12176.

    Article  CAS  PubMed  Google Scholar 

  54. Zhao B, Luo H, He J, Huang X, Chen S, Fu X, et al. Comprehensive transcriptome and Methylome analysis delineates the biological basis of hair follicle development and wool-related traits in merino sheep. BMC Biol. 2021;19:197. https://doi.org/10.1186/s12915-021-01127-9.

  55. Elgendy M, Sheridan C, Brumatti G, Martin SJ. Oncogenic ras-induced expression of Noxa and Beclin-1 promotes autophagic cell death and limits clonogenic survival. Mol Cell. 2011;42(1):23–35. https://doi.org/10.1016/j.molcel.2011.02.009.

    Article  CAS  PubMed  Google Scholar 

  56. Lowman XH, McDonnell MA, Kosloske A, Odumade OA, Jenness C, Karim CB, et al. The proapoptotic function of Noxa in human leukemia cells is regulated by the kinase Cdk5 and by glucose. Mol Cell. 2010;40(5):823–33. https://doi.org/10.1016/j.molcel.2010.11.035.

    Article  CAS  PubMed  Google Scholar 

  57. Janus P, Toma-Jonik A, Vydra N, Mrowiec K, Korfanty J, Chadalski M, et al. Pro-death signaling of cytoprotective heat shock factor 1: upregulation of NOXA leading to apoptosis in heat-sensitive cells. Cell Death Differ. 2020;27(7):2280–92. https://doi.org/10.1038/s41418-020-0501-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Wang P, Guo J, Wang F, Shi T, Ma D. Human SBK1 is dysregulated in multiple cancers and promotes survival of ovary cancer SK-OV-3 cells. Mol Biol Rep. 2011;38(5):3551–9. https://doi.org/10.1007/s11033-010-0465-8.

    Article  CAS  PubMed  Google Scholar 

  59. Sakabe I, Hu R, Jin L, Clarke R, Kasid UN. TMEM33: a new stress-inducible endoplasmic reticulum transmembrane protein and modulator of the unfolded protein response signaling. Breast Cancer Res TR. 2015;153(2):285–97. https://doi.org/10.1007/s10549-015-3536-7.

    Article  CAS  Google Scholar 

  60. Madsen ML, Nettleton D, Thacker EL, Edwards R, Minion FC. Transcriptional profiling of mycoplasma hyopneumoniae during heat shock using microarrays. Infect Immun. 2006;74(1):160–6. https://doi.org/10.1128/IAI.74.1.160-166.2006.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Friederich MW, Timal S, Powell CA, Dallabona C, Kurolap A, Palacios-Zambrano S, et al. Pathogenic variants in glutamyl-tRNAGln Amidotransferase subunits cause a lethal mitochondrial cardiomyopathy disorder. Nat Commun. 2018;9:4065. https://doi.org/10.1038/s41467-018-06250-w.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Liang W, Bi Y, Wang H, Dong S, Li K, Li J. Gene expression profiling of clostridium botulinum under heat shock stress. Biomed Res Int. 2013;2013:760904. https://doi.org/10.1155/2013/760904.

    Article  CAS  Google Scholar 

  63. Dou J, Khan A, Khan MZ, Mi S, Wang Y, Yu Y, et al. Heat stress impairs the physiological responses and regulates genes coding for extracellular exosomal proteins in rat. Genes-Basel. 2020;11(3):306. https://doi.org/10.3390/genes11030306.

    Article  CAS  PubMed Central  Google Scholar 

  64. Signorelli M, Ebrahimpoor M, Veth O, Hettne K, Verwey N, García‐Rodríguez R, et al. Peripheral blood transcriptome profiling enables monitoring disease progression in dystrophic mice and patients. EMBO Mol Med. 2021;13(4):e13328. https://doi.org/10.15252/emmm.202013328.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Toni M, Angiulli E, Miccoli G, Cioni C, Alleva E, Frabetti F, et al. Environmental temperature variation affects brain protein expression and cognitive abilities in adult zebrafish (danio Rerio): a proteomic and behavioural study. J Proteomics. 2019;204:103396. https://doi.org/10.1016/j.jprot.2019.103396.

    Article  CAS  PubMed  Google Scholar 

  66. Rouillard AD, Gundersen GW, Fernandez NF, Wang Z, Monteiro CD, McDermott MG, et al. The harmonizome: a collection of processed datasets gathered to serve and mine knowledge about genes and proteins. Database. 2016;2016:baw100. https://doi.org/10.1093/database/baw100.

  67. Sakai T. Development and regeneration of salivary gland toward for clinical application. Oral Sci Int. 2016;13(1):7–14. https://doi.org/10.1016/S1348-8643(15)00040-3.

    Article  Google Scholar 

  68. Tao YM, Huang JL, Zeng S, Zhang S, Fan XG, Wang ZM, et al. BTB/POZ domain-containing protein 7: epithelial-mesenchymal transition promoter and prognostic biomarker of hepatocellular carcinoma. Hepatology (Baltimore, Md). 2013;57(6):2326–37. https://doi.org/10.1002/hep.26268.

    Article  CAS  Google Scholar 

  69. Abbas Z, Sammad A, Hu L, Fang H, Xu Q, Wang Y. Glucose metabolism and dynamics of facilitative glucose transporters (GLUTs) under the influence of heat stress in dairy cattle. Metabolites. 2020;10(8):312. https://doi.org/10.3390/metabo10080312.

  70. Collier RJ, Baumgard LH, Zimbelman RB, Xiao Y. Heat stress: physiology of acclimation and adaptation. Anim Front. 2019;9:12–9.

    Article  Google Scholar 

Download references

Acknowledgements

The authors acknowledge the Beijing Dairy Cattle Center and Beijing Sunlon Livestock Development Company Limited for providing access to the datasets used and enabling data collection at their farms. The authors are also grateful to the members of the COWINFO group at the China Agriculture University for helping to measure rectal temperature in a large number of animals, especially Lu Cao, Zezhao Wang, Xinyi Yan, Hetian Huang, Rui Shi, and Ye Wang.

Funding

This research was funded by National Key R&D Program of China (2021YFD1200903, the earmarked fund for CARS36, the Program for Changjiang Scholar and Innovation Research Team in University (IRT_15R62).

Author information

Authors and Affiliations

  1. Laboratory of Animal Genetics, Breeding and Reproduction, Ministry of Agriculture of China, National Engineering Laboratory of Animal Breeding, College of Animal Science and Technology, China Agricultural University, No.2, Yuanmingyuanxi Road, Haidian District, Beijing, 100193, China

    Hanpeng Luo, Lirong Hu, Abdul Sammad, Yao Chang, Liwei Zhai & Yachun Wang

  2. Department of Animal Sciences, Purdue University, West Lafayette, IN, 47907, USA

    Luiz F. Brito

  3. College of Animal Science and Technology, Beijing University of Agriculture, Beijing, 100096, China

    Jinhuan Dou

  4. College of Animal Science, Xinjiang Agricultural University, Urumqi, 830052, China

    Longgang Ma

  5. Beijing Sunlon Livestock Development Company Limited, Beijing, 100029, China

    Gang Guo

  6. Beijing Dairy Cattle Center, Beijing, 100192, China

    Lin Liu

  7. College of Life Sciences and Bioengineering, Beijing Jiaotong University, Beijing, 100044, China

    Qing Xu

Authors
  1. Hanpeng Luo
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Lirong Hu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Luiz F. Brito
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Jinhuan Dou
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Abdul Sammad
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Yao Chang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Longgang Ma
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Gang Guo
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Lin Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Liwei Zhai
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Qing Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Yachun Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

HL, LB, QX, and YW conceived the study. HL, LH, JD, YC and LM analyzed the data. HL wrote the first draft of the manuscript. LB, YW, AS, and LZ edited the manuscript. GG and LL participated in the data collection. All the authors read and approved the final version of the manuscript.

Corresponding author

Correspondence to Yachun Wang.

Ethics declarations

Ethics approval and consent to participate

The study was approved by the ethical committee of the College of Animal Science and Technology, China Agricultural University, Beijing, China.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Supplementary Information

Additional file 1: Fig. S1

. The 68 significant GO terms revealed by 54 protein-coding genes associated with the three physiological traits.

Additional file 2: Table S1

. Differential gene expression between non-heat stress group and heat stress group.

Additional file 3: Fig. S2

. The significant GO terms and pathways were enriched by Differential gene expression between non-heat stress group and heat stress group.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Luo, H., Hu, L., Brito, L.F. et al. Weighted single-step GWAS and RNA sequencing reveals key candidate genes associated with physiological indicators of heat stress in Holstein cattle. J Animal Sci Biotechnol 13, 108 (2022). https://doi.org/10.1186/s40104-022-00748-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40104-022-00748-6

Keywords

Journal of Animal Science and Biotechnology

ISSN: 2049-1891