Genomic Regions and Candidate Genes Affecting Response to Heat Stress with Newcastle Virus Infection in Commercial Layers Chicks Using Chicken 600K SNP Array

1. Introduction

Global warming is accelerating at a rapid pace. Climate change adversely affects the production, reproduction, and growth of livestock including poultry [1]. Heat waves around the world cause huge economic losses in the poultry industry [2]. With the increasing population, there is high demand for a more efficient and secure food production system. Globally, poultry production makes important contributions to nutrition and food security and provides affordable and sustainable animal-source protein and multi-nutrients, especially for communities living in resource-limited regions [3]. Poultry meat and eggs are rich in protein, vitamins, and minerals and low in saturated fat [4]. The market demand for poultry products is high and expected to rise as a result of the high feed efficiency of poultry compared to other livestock animals [5]. The US is the world’s largest poultry meat producer with 17% of global output [6]. Poultry eggs are also essential as a low cost source of high quality animal protein [4]. However, heat stress impacts both meat and egg production in the US, resulting in decreases in egg production by 0.5 to 7.2% annually, with lower feed intake (-16.4%), growth rate (-32.6%), higher feed-to-gain ratio (+25.6%) and mortality [4,7,8,9]. Heat stress is also considered a major welfare concern for both backyard and commercial producers [1].

In addition to the negative impacts on productivity, heat stress also depresses immunity in chickens [10]. A reduction in immune organ weights and lower levels of circulating antibodies have been documented in both broilers and layers experiencing heat stress [11,12]. Under heat stress, birds maintain thermal homeostasis by increasing reactive oxygen species (ROS) levels and in turn increasing heat shock protein (HSP) production to protect themselves from the adverse cellular effects of the ROS [13]. Furthermore, it has been shown that HSPs might be the primary drivers of the immune response with both heat stress and lipopolysaccharide treatments in a macrophage-like cell line [14]. In addition to the effects caused by heat stress, poultry are subject to infectious disease outbreaks. In some African countries, there is a higher incidence of Newcastle disease virus (NDV) outbreaks during periods of higher temperatures [15]. Therefore, studies to understand the combined effects of abiotic stress (heat stress) and biotic stress (NDV infection) are warranted.

The contribution of genetics to heat stress and NDV infection has been reported in different chicken lines and breeds, demonstrating the potential for genetic improvement of heat tolerance and disease resistance in the species. Quantitative trait loci (QTL) for body temperature, body weight, and breast muscle yield were identified using a Fayoumi advanced intercross population [8]. In the same intercross population, 61 QTLs were associated with significant increases in pH, base excess (BE), HCO3, TCO2, iCa2+, and sO2, and significant decreases in PCO2 and glucose after heat treatment [16]. In another GWAS utilizing a white egg layer line, genomic regions associated with mortality from heat stress, body temperature, cholesterol levels, glucose levels, and bacterial and viral infection were identified [17]. Few studies have attempted to identify QTLs associated with a response to heat stress combined with NDV infection in chickens. Our previous GWAS identified QTLs associated with viral load, antibody response, and body weight change in Hy-Line Brown (HLB) chicks subjected to heat stress and NDV infection [18]. Using a high-density (600K) single nucleotide polymorphism (SNP) array, this study identified several important candidate genes associated with viral replication and immune response [18]. Therefore, as the follow-up, the identification of QTLs and candidate genes associated with heat tolerance and NDV resistance are needed to provide potential selection markers for improving heat tolerance in chickens. The current study is a related work aiming to identify genomic regions and candidate genes for blood physiological parameters associated with heat stress and NDV infection responses using the Affymetrix SNP array in HLB laying chicks.

2. Results

2.1. Effects of Treatments on Blood Parameters

Physiological responses to heat stress between the treated and control groups were compared based on blood parameter measurements. Least square means and P-values are listed in Supplementary File 1 for all blood parameters at each time point. The acute heat (AH) treatment didn’t have any significant effects on any parameters. The chronic heat (CH) treatment significantly reduced the BE levels (P < 0.05) (Figure 1). Female birds had significantly higher Na⁺ levels than male birds with the CH treatment (P < 0.05). The glucose level was significantly increased with the chronic heat stress combined with NDV infection (CH&NDV) treatment (P<0.05) (Figure 2). Sex had significant effects on TCO₂, HCO₃, PO₂, and BE. The interaction between the CH&NDV treatment and sex had significant effects on the PO₂ levels.

2.2. Genetic Parameters

The heritability for all phenotypic traits estimated by ASReml is presented in Table 1. The HCO₃ level with the AH treatment had the highest estimated heritability of 0.26 ± 0.09. Carbon dioxide related blood gas parameters (PCO₂, TCO₂, HCO₃, and BE) had relatively higher estimated heritabilities (0.02 to 0.26) than other parameters in all treatments. In general, the level of electrolytes had very low heritabilities (0 to 0.11) for all treatments.

2.3. Genome-Wide Association Analysis

A GWAS analysis was conducted on blood parameter traits using the same numbers of SNPs from the same population as in our previously reported study [18]. Suggestive SNPs were identified based on 95% and 99.5% of the variation between SNP for 5% (P-value = 1.108E-06) and 10% (P-value = 2.035E-06) genome-wide significance. GWAS was performed for each trait at each treatment stage. Of all 11 traits at three treatment stages (33 traits), with 99.5% of the variation, 5 traits (AH_sO₂, CH_K⁺, CH&NDV_HCO₃, CH&NDV_Na⁺, and CH&NDV_TCO₂) had 9 suggestive SNPs at 5% and 11 suggestive SNPs at 10% significance (Table 2). When we relaxed the significance threshold to 95% of the variation, 7 traits (AH_sO₂, CH_pH, CH_BE, CH_K⁺, CH&NDV_Na⁺, CH&NDV_HCO₃, and CH&NDV_TCO₂) had 39 suggestive SNPs at 5% and 11 traits (AH_sO₂, CH_pH, CH_HCO₃, CH_BE, CH_ Glu, CH_Na⁺, CH_K⁺, CH_TCO₂, CH&NDV_HCO₃, CH&NDV_Na⁺, and CH&NDV_TCO₂) had 53 suggestive SNPs at 10% significance (Table 2). Based on the numbers of the identified SNPs and potential candidate genes they could be associated with, SNPs identified with 95% of variation and 5% significance were used for further downstream analysis. Manhattan plots for significant SNPs above the selected cutoff are presented in Figure 3 A-G. Manhattan plots for other traits at three stages are presented in Supplementary File 2. SNP rs312361512 (P-value=3.28E-07) on GGA4 was the only suggestive SNP identified for sO₂ with the AH treatment. The K⁺ at the CH stage had the largest number of identified SNPs (16) on GGA1, in which all 16 SNPs were clustered into 2 QTLs with 3 and 13 SNPs, respectively. The Na⁺ at the CH&NDV stage had 15 significant SNPs identified on GGA 12 and all these SNPs were clustered as one QTL (Chr12: 12.5-13.5 Mb). Three significant SNPs (rs316466702, rs31814125, and rs314444945) on GGA 3 were identified for HCO3 at the CH&DNV stage were also showing significance for TCO₂ at the same stage. One suggestive SNP (rs316524480) on the same chromosome was only identified for HCO₃ at the CH&NDV stage. All significant SNPs with GGA coordinates and significance are listed in Supplementary File 3.

A 1 Mb region of each SNP was examined to identify biologically relevant genes that may be responsible for the SNP associations for each phenotype [18]. For the 39 SNPs associated with 7 traits, there were 157 candidate genes within the boundary (1 Mb window centered by each SNP). Specific associated gene numbers are presented in Table 3. Detailed information on gene names and known functions is in Supplementary File 4. The SNP, rs312361512, was associated with sO₂ at the AH stage. There were 25 genes within a 0.5 MB region on either side of the SNP (rs312361512, associated with sO₂ at the AH stage) identified on GGA4, including hormone and iron transporter genes (solute carrier family 16 member 2 (SLC16A2) and magnesium transporter 1 (MAGT1)). Two QTLs associated with the K⁺ levels at the CH stage on GGA 1 were identified, with 63 genes located within the region, including candidate genes such as heat shock transcription factor 2 binding protein (HSF2BP), mitochondrial calcium uptake 2 (MICU2), intraflagellar transport 88 (IFT88), inducible T-cell costimulatory ligand (ICOSLG) and interleukin 17 D (IL17D). These genes are involved in metabolic and immune functions. Even though there were 15 SNPs identified as associated with the Na⁺ level for the CH&NDV treatment, only 14 candidate genes were in this region since these SNPs are clustered in a small region (0.096 Mb). Three SNPs associated with both CH&NDV_HCO₃ and CH&NDV_TCO₂ located on GGA3 with 13 genes including ELOVL fatty acid elongase 5 (ELOVL5) and other 6 glutathione transferase genes (glutathione transferase (GST), glutathione S-transferase-like (GSTL), glutathione S-transferase class-alpha (GSTA), GSTA2, GSTA3, and GSTA4. The SNP, rs316524480, was only identified by CH&NDV_HCO₃ and had three adjacent genes, 1-acylglycerol-3-phosphate O-acyltransferase 5 (AGPAT5), microcephalin 1 (MCPH1), and angiopoietin 2 (ANGPT2).

3. Discussion

Populations with different genetic backgrounds express distinct physiological responses to heat stress [19]. However, the underlying molecular mechanisms by which birds show differences in the heat stress response are not fully understood. Mapping genetic markers and identifying candidate genes associated with heat tolerance will provide novel insights to develop alternative approaches to mitigate climate change impacts. Combining information on positional candidate genes with that of functional candidate genes, which respond to the same treatements by differential gene expression, can give a more comprehensive understanding of the genetics associated with heat tolerance and NDV response. The use of different population structures can facilitate these studies. Populations with segregating genetic variation are a strong foundation for identifying structural variants associated with measured phenotypes in GWAS. However, complex population genetics may also yield more subtle treatment responses for complex phenotypes. In contrast, use of highly inbred research lines cannot be used for GWAS, but their very narrow genetic bases and distinct phenotypes can facilitate the identification of functional candidate genes. We integrated data of our current study with that of previous studies of inbred lines tested under the same treatment protocols, to determine overlap in positional and functional genes associated with heat stress and NDV. Based on identification of strong candidate markers or genes, breeding programs could be implemented to improve production traits in commercial lines, reducing morbidity and mortality and increasing productivity.

3.1. Blood Parameter Measurements

The combined heat stress and NDV infection treatment model developed by our group was previously applied to two genetically and phenotypically distinct inbred lines, Fayoumi (relatively resistant to heat stress and NDV infection) and Leghorn (relatively susceptible to heat stress and NDV infection) to characterize the different effects derived from the genetic background [20,21,22,23]. Physiological responses demonstrated changes in blood parameters between the two lines [24]. In contrast to this previous study utilizing inbred lines, the current study utilized a commercial brown egg production hybrid. The three treatments, AH, CH, and CH&NDV had only subtle effects on blood parameters in the HLB commercial cross birds (Table 4). Two parameters were significantly altered, which were a significant reduction of BE with the CH treatment and an increase in glucose levels with the CH&NDV treatment. The HLB birds responded to CH treatment with low levels of BE (similar to Fayoumis) and to CH&NDV treatment with high levels of glucose (similar to Leghorns) [24]. The BE was always negative in HLB birds and was significantly decreased with the CH treatment. Decreased BE indicates respiratory alkalosis and metabolic acidosis [25]. Fayoumi birds decreased levels of blood carbon dioxide related parameters to compensate for the respiratory alkalosis due to heat stress across all treatments [24]. HLB birds might have similar mechanisms to respond to heat stress by decreasing BE levels but without changing other carbon dioxide related parameters. Glucose levels in HLB birds were significantly increased with the CH&NDV treatment, which was similar to Leghorn birds. Increases in glucose levels are common responses to heat stress in birds [26,27], however, Fayoumi birds were able to keep stable glucose levels with treatments [24]. Therefore, based on these blood phenotypes, HLB birds are relatively heat resilient in general.

Heritabilities estimated from 11 traits during the three treatments, in general, were from low to moderate. Electrolytes and oxygen related traits (PO₂ and sO₂) had the lowest heritabilities (0-0.16). Carbon dioxide related parameters (PCO₂, TCO₂, HCO₃, and BE) had relatively higher heritabilities (0.08-0.26). These results were consistent with heritability estimations for these parameters reported in other studies [16,28]. Therefore, genetic selection to improve heat tolerance focusing on some of these blood parameters, in specific BE, may be feasible.

3.2. GWAS Analysis and Candidate Gene Identification

Genome-wide association studies on heat resistance have been widely studied [8,16,17,18,29,30,31,32,33,34]. SNPs in candidate genes associated with body temperature were identified by GWAS in a broiler population under acute heat stress. These significant SNPs were found to be related to apoptosis or responses to external stimuli and can be used as potential candidates for thermotolerance in chickens [32]. Heat stress associated QTLs were identified for blood components in an advanced intercross chicken line under heat stress in 2016 [16]. Heritabilities were estimated for blood parameters during heat stress, which was similar to the current study. Sixty-one QTLs and 999 candidate genes in the QTL regions were identified to be associated with heat stress [16]. It showed us the opportunities to find markers associated with multiple blood parameters that could genetically contribute to heat tolerance.

Our previous GWAS analyses focused on NDV infection and body weight associated traits in this HLB population [18]. The current GWAS aimed to identify genomic regions associated with heat tolerance, which revealed 39 significant SNPs for 7 traits located on GGA1, 3, 4, 6, 11, and 12. Among these significant SNPs, 1 SNP (rs15615489) was associated with CH_BE, which was significantly reduced by the CH treatment. Positional candidate genes were identified within 1 Mb of each significant SNP for the evaluated phenotypes. Expression level changes of these genes under heat stress in biologically relevant tissues may provide additional evidence that these genes could contribute to heat tolerance [35]. Our group has reported several transcriptome analyses on a variety of tissues of the two genetically distinct highly inbred lines, Fayoumi (resistant) and Leghorn (susceptible) with the same experimental treatments [20,21,22,23]. Therefore, differential gene expression data from the metabolic tissues of the liver, hypothalamus, and breast muscle, were integrated with the GWAS results to identify candidate genes associated with specific traits under heat stress and NDV infection.

The only suggestive SNP (rs312361512, Chr4: 12459540) showed an association with sO₂ during AH treatment of HLB. The blood sO₂ is the amount of oxyhemoglobin that indicates the capability of oxygen binding by the hemoglobin [36]. Therefore, candidate genes close to this SNP could be genetic markers associated with hemoglobin function during oxygen transportation. Within 1 Mb of this SNP, there were 25 genes identified. Twenty-one out of twenty-five genes are annotated genes, including some metabolism-related genes such as vascular endothelial growth factor receptor kdr-like (VGFR4), ras-like protein family member 11A-like, solute carrier family 16 member 2 (SLC16A2), and cytochrome oxidase subunit 7B (COX7B) genes. Of particular interest were two genes with differential expression in inbred lines under the same treatment: VGFR4 and COX7B. VGFR4 was upregulated in Leghorn livers with acute heat stress and COX7B was down-regulated in Fayoumi livers with chronic heat stress [37]. VGFR4 has tyrosine-protein kinase activity and is essential for new blood vessels [38]. COX7B is a component of the cytochrome c oxidase, which can catalyze the reduction of oxygen in the water [39]. Further investigation on these genes can help to elucidate their roles in regulating blood oxygen during heat stress.

With the CH&NDV treatment, three significant SNPs closely located on GGA3 were identified affecting both TCO₂ and HCO₃ levels in HLB. One specific SNP only associated with HCO₃ was identified a little further downstream of the same chromosome. TCO₂ and HCO₃ were all calculated from PCO₂ [40]. Both TCO₂ and HCO₃ are important for evaluating the acid-base and electrolyte imbalance. Therefore, these four suggestive SNPs could be candidate genetic markers for improving the homeostasis of blood acid-base and electrolyte balance. Within 1 Mb of the four significant SNPs, there were 16 genes identified in total. The major biological function enriched by these 16 genes was glutathione and xenobiotic metabolic processes through drug and glutathione metabolic pathways (Figure 4). Three of the 16 genes (AGPAT5, MCPH1, and ANGPT2) were specifically associated with the HCO₃-specific SNP (rs316524480). Of the remaining 13 genes, 6 of them are glutathione transferase genes (GSTs). Glutathione (GSH) is a powerful endogenous antioxidant [41]. In humans, GSH blood concentration can be used as a biomarker for oxidative damage protection [42]. GSH is involved in cellular defense mechanisms and the metabolism of xenobiotic compounds [43]. GSTs are a family of phase II detoxification enzymes that can eliminate toxic compounds by catalyzing the conjugation of GSH to electrophilic compounds [44], by which they can protect cells during oxidative stress [45]. GSTs are regulators of many stress-induced signaling pathways, such as the mitogen-activated protein kinase (MAPK) signaling pathway, the JNK signaling pathway, and the rapamycin (mTOR)-signaling pathway [44,45,46,47]. Modulation of blood GSH levels affects the HCO₃ levels in the plasma of rats [48]. Therefore, it is possible that blood HCO₃ and TCO₂ levels are highly correlated with GSH levels in chickens. More research is warranted to determine if blood concentrations of TCO₂ and HCO₃ could be used as biomarkers of oxidative damage from heat stress. Two out of the 16 genes, ANGPT2 and GCM1, showed differential expression patterns in the transcriptome analysis in the inbred lines [37]. ANGPT2 gene was down-regulated in Fayoumi breast muscle during CH&NDV treatment. This gene is always upregulated in multiple inflammatory diseases and is in the direct control of inflammation-related signaling pathways. It can even serve as a prognostic biomarker for acute respiratory distress syndrome [49]. Expression levels of ANGPT2 were decreased in treated Fayoumi birds compared to non-treated ones, which possibly indicated that Fayoumis reduced inflammation to protect themselves from the effects of heat stress and NDV infection.

With the highest number of significant SNPs (16 SNPs) associated with CH_K⁺, 63 genes were identified. The major biological function enriched by these genes was the cell-cell signaling by the GO analysis. Of these 63 genes, some of them are metabolism related such as DEAD-box helicase 10 (DDX10), heat shock transcription factor 2 binding protein (HSF2BP), solute carrier family 35 member F2 (SLC35F2). DDX10 was reported in a gene module regulating ubiquitin-mediated proteolysis during heat stress in chickens [50]. SLC35F2 was differentially expressed in chicken testis tissue in response to an acute heat stress treatment [51]. Some others are immune-related genes such as ICOSLG and IL17D. These results demonstrate that blood K⁺ levels may affect both metabolic and immune functions. When we overlaid the positional candidate gene list from the HLB study with the transcriptome data of the inbred line study, 10 genes showed differential expression in at least one contrast. Zinc finger DHHC-type containing 20 (ZDHHC20) gene was upregulated with the CH&NDV treatment in three contrasts: Fayoumi liver, Leghorn and Fayoumi breast muscle. ZDHHC20 was reported as induced by heat stress in chickens [52], which was consistent with our data. However, the relation between ZDHHC20 gene and the K⁺ level is unclear. Therefore, further investigation of the role of ZDHHC20 on heat stress is warranted.

The only significant SNP associated with the alteration of the BE levels due to the CH treatment in HLB, rs15615489, had 21 associated candidate genes. Within these candidate genes, zinc finger protein 507 (ZNF507) was significantly up-regulated in Fayoumi livers with chronic heat stress combined with NDV infection. ZNF507 is predicted to facilitate DNA binding and metal ion binding activities and can regulate transcription in mammals [53]. Whether higher expression levels of ZNF507 correlated with lower BE levels indicate heat resilience should be further investigated.

4. Materials and Methods

4.1. Experimental Population

Hy-Line Brown laying chicks (Hy-Line International) generated from 145 dams and pooled semen were utilized in this experiment, with three replicate trials (~180 birds per trial). Data for the current study were collected from the same animal trial we reported previously [18]. All animal protocols performed were approved by the Institutional Animal Care and Use Committee, University of California, Davis (IACUC #17853) [18]. In brief, chickens were housed in temperature and humidity-controlled chambers in a BSL-2 animal facility at the University of California, Davis. Twenty individuals per trial, randomly selected from 60 different dams, were housed in separate chambers and used as a control group and housed in one chamber. The remaining ~160 birds per trial were used as the treatment group and housed in another two chambers. Both control and treated groups were housed at the same temperature and humidity (32 °C for the first 7 days and then decreased to 28 °C) until day 14. From day 14 to the end of the experiment, the treatment groups received acute heat stress (35 °C and 50% humidity (AH), 4 hour post heat stress), chronic heat stress (32 °C and 50% humidity (CH), 6 days post heat stress), followed by chronic heat stress combined with NDV infection (32 °C, 50% humidity, and inoculated with 200 ul 10⁷ EID₅₀ of the La Sota strain of NDV (CH&NDV), 9 days post heat stress and 2 days post NDV inoculation) treatments. Six birds per group were euthanized for tissue collection (RNA-seq analysis in a separate study) on day 14 and day 23 from the treated and control groups [18].

4.2. Blood Parameter Measurements

Blood parameters were measured on all birds on day 14 of age (or AH for the treated group), day 20 of age (or CH for the treated group), and day 23 of age (or CH&NDV for the treated group). One ml of blood was collected into a heparinized syringe (2,000 unit/ml heparin in PBS, Sigma-Aldrich, St. Louis, MO) from each bird and analyzed immediately using an i-STAT Portable Blood Analyzer (Abbott Laboratories, San Diego, CA) [24]. Eleven blood parameters were measured by the i-STAT cartridge (CG8+), including four chemistry/electrolyte parameters (concentrations of sodium (Na⁺), potassium (K⁺), ionized calcium (Ca²⁺), and glucose (Glu)); and seven blood gas parameters (blood pH, carbon dioxide partial pressure (PCO₂), oxygen partial pressure (PO₂), total carbon dioxide (TCO₂), bicarbonate (HCO₃), base excess (BE), and oxygen saturation (sO₂)).

4.3. Genotyping

Genotyping was performed using the same methodology as previously described [18]. Genomic DNA was extracted and genotyped on the Axiom® 600K Genome-Wide Chicken Array Kit (Affymetrix CAT#902148) by Geneseek (Lincoln, NE). After the filtration of SNPs based on a call rate greater than 95% and a minor allele frequency of greater than 0.01, genotypes for 304,500 SNPs on 526 individuals were used in the downstream analysis. SNP positions were mapped to the genome based on the galGal 6 reference genome.

4.4. Data Analysis

The blood parameter data were analyzed individually for each time point using the following mixed linear model using the JMP 16 software (SAS Institute, 2022).

$$Y=\mu +T+S+Rep+T\times S+T\times Rep+Rep\left(C\right)+D+e$$

where Y is the dependent variable of various parameters, μ is the population mean of the measurements, T is the fixed effect associated with treatments, S is the fixed effect associated with sex, Rep is the random effect of the three replicates, T × Rep is the interaction effect of treatments and replicates, T × S is the interaction effect of treatments and sex, $Rep\left(C\right)$ is the random effect of chamber nested within replicates, D is the random effect of dams, and e is the random residual effect. All measurements in this study were represented as least square mean ± SE. P < 0.05 was considered statistically significant. All figures were generated by Prism 6 GraphPad software (La Jolla, CA).

Heritabilities of traits for each time point were estimated using ASReml4 using the same univariate animal model as used by [18]:

$$Y_{ijkl}=\mu +S_i+{Rep\left(C\right)}_{jk}+A_l+e_{ijkl}$$

Y is the dependent variable of the analyzed trait. Sex (S), and chamber, nested within replicate (Rep(C)), were fitted as fixed effects. Random effects included animal genetic effects (A) using a genomic relationship matrix [54] and residuals (e). Heritability was calculated as a ratio of animal variance to phenotypic variance.

Genome-wide association analyses were performed using the R package GenABEL described in our previous publication [18]. To establish significance thresholds, the number of independent tests was determined as the number of principal components that accounted for 95% of the variance in SNP genotypes, as described by Waide et al [55]. A total of 49,138 principal components were determined to account for 95% of the variance between SNPs. Based on the number of independent tests, 99.5 and 95% of the variance between SNPs were determined. The number of independent tests was used to determine the Bonferroni correction, the P-value for the 5 and 10% genome-wide significance levels were 1.018E-06 and 2.035E-06, respectively. The Affymetrix SNP IDs were converted into SNP ID’s using NCBI’s dbSNP database. A 1 Mb window centered around significant SNPs was used and SNPs were associated with genes using the Ensembl Variant Effect Predictor web tool [19]. Gene ontology (GO) analyses were carried out on identified candidate genes for each trait using DAVID [56,57]. GO terms in the biological process enriched by important SNP associated candidate genes with fold enrichment > 2 and FDR < 20% were considered significant overrepresentation. The linkage disequilibrium function in PLINK was used for the estimation of r² values for all pairs of SNP markers within the 1 Mb window size in R [58].

5. Conclusions

Heat stress and NDV infection generated only subtle changes in the blood parameters of HLB chicks. Relatively low to moderate heritabilities estimated for each trait demonstrate that it might be feasible, yet challenging, to genetically improve heat tolerance by using blood parameters alone as selection criteria. However, candidate genes within QTL regions identified in this study can provide more insight into identifying genetic markers, which may be used for genetic selection to improve heat tolerance in poultry. Specifically, a consensus candidate gene list was generated of genes meeting all three criteria: 1) located on or close to significant SNPs associated with significantly changed blood parameters; 2) has reported metabolic or immune functions; 3) differentially expressed in metabolic-related tissues with treatment. The 19 consensus genes identified in this study (Supplementary file 5) can form the basis of hypothesis-driven future investigations of molecular mechanisms underlying heat tolerance in chickens. Heat stress associated candidate SNPs and genomic regions identified can aid the development of targeted SNP panels to assist genomic selection to improve poultry production in regions under climate stress.

6. Patents

This section is not mandatory but may be added if there are patents resulting from the work reported in this manuscript.