Growth resilience to weather variation in commercial free-ranging chickens in Ethiopia

TL;DR

This study explores how commercial chickens in Ethiopia adapt to weather changes and identifies genetic factors that could help breed more resilient chickens.

Contribution

Novel resilience phenotypes were developed to describe chicken growth responses to weather fluctuations, revealing genetic markers linked to environmental adaptation.

Findings

Variations in temperature, humidity, and precipitation strongly affect chicken growth resilience.

Genomic markers for resilience traits are linked to genes involved in lipid metabolism and adipocyte homeostasis.

Significant genomic variation was found for resilience to temperature and humidity changes.

Abstract

The poultry industry in sub-Saharan Africa is a rapidly developing sector mostly based on smallholder farming. Increased demand for poultry-derived products, driven by the growing economy and population, has intensified importations of highly productive exotic breeds and crossbreeding with local ecotypes. However, commercial chickens with exotic genes often struggle to adapt to the local climate under smallholder farmers management. Understanding the chicken response to weather changes is crucial for developing selection schemes that ensure proper adaptation. In the present study, we derived individual phenotypes for growth resilience of commercial free-ranging chickens to changing weather conditions in Ethiopia. In addition, we performed genomic association analyses to assess the genetic background of these phenotypes and identify potential candidate genes of interest.

Novel resilience phenotypes describing changes in chicken growth profiles in response to weather fluctuation were developed. Variations in daily air temperature, relative humidity and amount of precipitation had the strongest impact on growth. Significant genomic variance was detected for growth resilience to changes in air temperature measurements and a temperature-humidity index. Genomic markers correlated with these resilience traits were mostly located within or near candidate genes associated with lipid metabolism and adipocyte homeostasis. Some of these genes have been previously linked to animal responses to environmental stressors in other species.

The phenotypes of growth resilience of chickens to changing weather conditions exhibited significant genomic variation. The outcomes of this study may facilitate the genomic selection of commercial chickens that are not only highly productive, but also capable of maintaining their production levels under varying weather conditions.

The online version contains supplementary material available at 10.1186/s12864-025-11561-6.

Figures

Population reaction norms at the phenotypic level for each weather variable. Air temperatures are in degrees centigrade (^o^C), relative humidity in percentage, and amount of precipitation in millimetres (mm); growth is in grams per day

Full text

Background

The agricultural sector in sub-Saharan Africa employs about 64% of the labour force and contributes to a substantial proportion of the Gross Domestic Product [1]. Within this sector, the poultry sub-sector is characterised by a diverse array of production systems, ranging from small family farms to large-scale commercial enterprises. The sector plays a vital role, not only as a food protein source, but also as a driver of social development and income. Approximately 80% of the total poultry production is associated with smallholder farming [2] and around 70% of the total chicken population are of indigenous breeds and ecotypes, although these proportions are highly variable across different regions [3].

In Ethiopia, the current national poultry production sector comprises around 56 million birds producing 77,000 tons of meat annually, of which 88% emanate from smallholder farms where chickens range freely [4]. Although the majority of these chickens are of indigenous ecotypes, the use of commercial crosses and imported exotic breeds is rapidly increasing due to their higher productivity [5], currently accounting for about 18% of the reared population [6]. This increase is fuelled by an expanding economy, a rapidly growing consumer population, and urbanisation. Indeed, demand for chicken meat in Ethiopia is estimated to increase by 268% in the next 25 years [7], prompting considerable pressure on the existing production systems. However, imported exotic breeds may be maladapted to the new local conditions and environmental challenges [8]. Climate change is a key challenge whose effect in Ethiopia is expected to be critically evident in the near future, with average air temperatures potentially increasing by 2.9 degrees Celsius by 2050 [9], accompanied by intense and unpredictable fluctuation in weather events [10]. These changes constitute a serious environmental stressor on chicken growth and production capacity, and commercial crosses based on exotic genes may be particularly vulnerable.

Consequently, breeding goals and methods must evolve to ensure selected individuals are capable of coping with the new challenges and avoid losses in chicken numbers and productivity [11]. A recent study on Ethiopian indigenous chickens [12] highlighted their genetic resilience to diverse environmental factors including air temperature, rainfall and land cover, and linked candidate genes to genetic adaptation to local climatic challenges. Thus, there is scope to develop novel chicken resilience phenotypes to inform existing and future breeding [13] and conservation [14] programmes aiming to mitigate the effect of climate change.

The aim of the present study was to investigate the impact of changing weather conditions on chicken growth in a popular commercial cross reared in Ethiopia by (i) deriving and examining novel chicken growth resilience phenotypes to weather volatility and (ii) assessing the genomic architecture of these new phenotypic traits of chickens. We deployed random regression methods on longitudinal data to derive resilience phenotypes and used genotype-by-sequencing data and genomic association analyses to estimate genomic parameters and identify potential molecular markers and candidate genes associated with chicken resilience to weather.

Material and methods

Animal data

Individual data on 1,590 chickens were available for the present study. All chickens were T451A Sasso, which is a slow-growing commercial dual-purpose type suitable for smallholder village conditions [15], resulting from a cross between T44 males and S51A females [16–18].

Data collection had taken place under a previous study [16] at the International Livestock Research Institute (ILRI) poultry research facility in Addis Ababa, Ethiopia, between 2019 and 2021 (Table 1). In that previous study, chickens grew in four batches outdoors in specially designed paddocks emulating the free-ranging semi-scavenging conditions that prevail in smallholder farms [13]. This phase started at the chicken age of 56 days, which represents the average age when chickens are sold to village smallholders for rearing in free range and lasted eight weeks. Live body weight was measured in grams and recorded weekly on each individual chicken together with the corresponding tag identifier, date of measurement, batch number and sex of the bird. Table 1. Summary of data used in the studyBatch numberNumber of chickensRecording start dateRecording end dateLocal season143618/12/201912/02/2020Bega-Belg242913/07/202007/09/2020Kiremt340614/10/202009/12/2020Bega431524/12/202018/02/2021Bega-BelgLocal seasons pertain to the calendar months of weekly chicken body weight measurement and recording: Bega (October to January), Belg (February to May) and Kiremt (June to September)

For the purposes of the present study, average daily growth rate was calculated for each chicken and week based on the live body weight difference between two consecutive weeks divided by seven. Thus, seven weekly measures of daily growth rate were available for each chicken.

Chicken genotypes from low-pass whole-genome sequencing and imputation had been produced and quality assured in the previous study [16]. Briefly, individual blood samples of 100μl were collected from each chicken, preserved in absolute ethanol, transferred onto QIAcard FTA Elute Micro cards, and transported to Neogen Genomics (USA) for low-pass whole-genome sequencing (0.5X coverage). SNP imputation was conducted by Gencove (USA) using a reference panel of 583 high-coverage (30X) chicken genomes, including the studied Sasso type, mapped to the Gallus gallus 6 assembly (GRCg6a), with a panel density of approximately 29 million variants. After assessing the accuracy of imputation, which was greater than 90% [16], quality control removed samples with call rates less than 90%, Single Nucleotide Polymorphism (SNP) markers with minor allele frequencies less than 0.02, and SNPs located on the sex chromosomes. Linkage Disequilibrium pruning (r^2^ > 0.8) was implemented. Duplicate samples with an identity-by-state estimate greater than 98% were also removed. After quality control, the final dataset for the present study comprised 1,586 genotyped chickens and 2,940,003 SNPs across 33 autosomal chromosomes.

Weather data

Daily records of maximum, minimum and average air temperature, relative humidity, and amount of precipitation for the period of study (2019–2021) were extracted from three main weather stations in Addis Ababa (Supplementary File 1) as well as satellites (NASA and Copernicus EU program) and radar [19].

Average daily air temperature and relative humidity records were combined to develop a temperature-humidity index (THI) using the Finocchiaro method [20]. THI was added to the list of weather variables.

The profile of each weather variable was explored by calendar year and the following three local seasons: Bega (from October to January), Belg (from February to May) and Kiremt (from June to September). Bega is often characterised by mild average and maximum daily air temperatures, relatively low minimum temperatures at night, and low humidity. The Belg season is characterised by variable humidity, cold nights, and mild daily average temperatures, but certain days tend to reach relatively high maximum temperatures. Kiremt is the wet season, with frequent heavy rain, high humidity, and relatively mild air temperatures.

Weekly means and variances of the daily values of each weather variable were calculated and matched to the corresponding weekly measurements of chicken daily growth for the ensuing data analyses.

Population reaction norms at the phenotypic level

A regression model (1) was used to describe average daily growth change in response to changing weather, manifested separately in each weather variable, at the population level [21]:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${y}_{ij}=X+f\left(\beta ,{X}_{j}\right)+{{a}_{i}+ e}_{ij}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${y}_{ij}$$\end{document} is the weekly performance record (average daily growth rate) of individual chicken i in relation to the corresponding weekly mean or variance of weather variable* j*; X corresponds to a set of fixed effects including the sex of chicken and a combined batch number and local season effect (6 levels); \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f\left(\beta ,{X}_{j}\right)$$\end{document} is a covariate that corresponds to the population norm function describing the relationship between the average daily growth of all chickens and the corresponding mean or variance of weather variable j in the same week; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${a}_{i}$$\end{document} corresponds to the deviation of individual chicken \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i$$\end{document} from the population average growth; and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${e}_{ij}$$\end{document} corresponds to the random residual.

The function linking chicken growth to the weather variables in model (1) was a second-degree Legendre polynomial, implying a quadratic non-linear relationship between chicken growth profile and weather change. No genotypic information was included in this step and, therefore, population norms were expressed at the phenotypic level.

Model (1) equations were solved using the EM-REML algorithm implemented in the BLUPF90 suite of software [22] and produced estimates of the three coefficients (intercept, linear and quadratic) of the population norm.

Individual deviation norms at the phenotypic level

A random regression model (2) was used to assess changes in individual chicken growth profiles in response to changing weather [21]:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${y}_{ij}=X+f\left(\beta ,{X}_{j}\right)+{f}_{i}\left({a}_{i},{X}_{j}\right)+{e}_{ij}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${f}_{i}\left({a}_{i},{X}_{j}\right)$$\end{document} corresponds to a second-degree Legendre polynomial function describing the relationship between the daily growth of chicken i and weather variable j, expressed as a deviation from the population norm \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f\left(\beta ,{X}_{j}\right)$$\end{document} ; all other effects are as in model (1). Again, as no genotypic information was included, individual deviation norms were expressed at the phenotypic level.

Solutions of the model (2) equations for the individual deviation coefficients (intercept, linear and quadratic) were derived using the EM-REML algorithm implemented in BLUPF90 [22]. Pearson correlations of the intercept with the linear and quadratic coefficients were estimated. Approximate standard errors of these correlations were also calculated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(1-{r}^{2})\sqrt{(n-3)}$$\end{document} , where r is the correlation estimate and n is the sample size [23].

Derivation of growth resilience phenotypes

Individual reaction norms were calculated by adding the individual deviations from model (2) to the corresponding population norm from model (1). One individual reaction norm was calculated for each chicken and the weekly mean of each weather variable, and another for each chicken and the weekly variance of the weather variable.

The points across a reaction norm curve correspond to chicken growth at the respective values of the weather variable. The slope of the tangent to the reaction norm at a given point of the curve represents the change in the chicken’s growth profile in response to changes in the specific weather value. This slope constitutes a potential growth resilience phenotype. Theoretically, there can be an infinite number of such phenotypes on each individual. For the purposes of the present study, we derived two resilience phenotypes for each chicken and individual reaction norm by calculating the derivatives at two points of the quadratic norm, each on either side of the absolute maximum value.

Genomic analyses of resilience phenotypes

Variance component and heritability estimates of each growth resilience phenotype were derived with mixed model (3), using the GCTA software [24]:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{y}}={\varvec{W}}\boldsymbol{\alpha }+{\varvec{Z}}{\varvec{u}}+{\varvec{\varepsilon}}$$\end{document}

where y represents the vector of growth resilience phenotypic values, W is the incidence matrix for vector α of fixed effects, Z is the incidence matrix for vector u of random polygenic effects (distributed as a multivariate normal distribution MVN(0,VgG) with G being the chicken genomic relatedness matrix and Vg the genomic variance of the trait), and ε represents the vector of residual errors (distributed as MVN(0,VeI) with I being the identity matrix and Ve the residual variance of the trait).

As the growth resilience phenotypes had already been corrected for the effects of sex of chicken, batch number and local season during the reaction norm calculation with models (1) and (2), there was no need to consider these effects again. However, we fitted the first three principal components (PC) derived from a PC analysis of the genotypes with software GCTA [24] as fixed effects in model (3) to account for any remaining population structure after including the genomic relatedness matrix G.

The statistical significance of the variance component and heritability estimates derived with model (3) was assessed with the likelihood ratio test [25].

Subsequently, genome-wide association studies were performed on growth resilience phenotypes that had exhibited statistically significant (P < 0.05) genomic variance in the analysis with model (3). The following linear mixed model and GEMMA software [26] were used:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{y}}={\varvec{W}}\boldsymbol{\alpha }+{\varvec{x}}{\varvec{\beta}}+{\varvec{Z}}{\varvec{u}}+{\varvec{\varepsilon}}$$\end{document}

where x represents the vector of genotypes at a SNP locus (coded as 0/1/2 according to the number of copies of the minor allele), β is the regression coefficient of the phenotype on the genotypes, and all other effects are as described in model (3).

A Bonferroni correction was applied for multiple testing to determine a genome-wide (P < 0.05) and suggestive (one false positive per genome scan) significance threshold, resulting in final threshold values of 1.7E-8 and 3.4E-7, respectively. The proportion of phenotypic variance explained by each of the identified significant SNPs was calculated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[{2\beta }^{2}\times maf\times (1-maf)] / [{2\beta }^{2}\times maf\times (1-maf)+{\left(\text{se}\left(\beta \right)\right)}^{2}\times 2N\times maf\times (1-maf)]$$\end{document} , where maf was the SNP minor allele frequency, β and se(β) the SNP effect and standard error, and N the sample size [27]. Furthermore, annotated genes neighbouring the identified SNPs were investigated using the Ensembl database.

Results

Weather data exploration

Supplementary files 2 and 3 illustrate results from the exploratory analysis of seasonal variation of weather variables within and across years for the study period (2019–2021).

Air temperatures ranged between 4 and 33°C. Average daily temperatures were generally within the 15–20°C range but could dip lower during the cold local season (Bega) and rise higher during the hot local season (Belg). Minimum daily temperatures were relatively low, averaging below 10°C during the Bega local season. Maximum daily temperatures were usually mild, with an average below 25°C during the Bega and Kiremt local seasons but reaching higher values during the Belg season. Daily relative humidity was mostly in the range of 30%-60% during the Bega and Belg seasons but increased considerably during the wet local season (Kiremt), with values reaching over 75%. In terms of average daily precipitation, high values were relatively frequent during the Kiremt season, reflecting strong rainfalls (over 30 mm).

Table 2 summarises weekly means and variances of the studied weather variables corresponding to the weekly chicken growth measurements. Table 2. Summary statistics for weekly means and variances of weather variables during the study period by local seasonMaximum air temperature (Weekly mean of daily values)Maximum air temperature (Weekly variance of daily values)Local SeasonMinimumMedianMeanMaximumMinimumMedianMeanMaximum****Bega22.2023.5723.4925.020.180.981.063.02Belg23.3725.4725.5427.250.691.201.945.60Kiremt19.7720.9821.1022.470.301.551.623.19Minimum air temperature (Weekly mean of daily values)Minimum air temperature (Weekly variance of daily values)Local SeasonMinimumMedianMeanMaximumMinimumMedianMeanMaximum****Bega7.339.279.4912.730.251.281.464.55Belg7.7511.0210.9212.630.371.171.081.70Kiremt12.2213.1212.9513.600.170.540.561.14Average air temperature (Weekly mean of daily values)Average air temperature (Weekly variance of daily values)Local SeasonMinimumMedianMeanMaximumMinimumMedianMeanMaximum****Bega14.7016.0016.1117.280.050.380.431.41Belg15.8717.7217.3218.270.250.290.662.47Kiremt15.4516.2516.2017.220.120.370.420.75Average relative humidity (Weekly mean of daily values)Average relative humidity (Weekly variance of daily values)Local SeasonMinimumMedianMeanMaximumMinimumMedianMeanMaximum****Bega47.1561.7059.1270.923.8621.3033.64147.44Belg41.6354.1254.9267.5013.9647.0856.13117.55Kiremt79.8383.2083.3186.973.276.789.9123.94Average THI (Weekly mean of daily values)Average THI (Weekly variance of daily values)Local SeasonMinimumMedianMeanMaximumMinimumMedianMeanMaximum****Bega14.6315.6915.7316.660.030.220.260.88Belg15.3917.0116.5717.250.070.110.331.41Kiremt15.3716.0416.0216.920.090.280.370.60Average precipitation (Weekly mean of daily values)Average precipitation (Weekly variance of daily values)Local SeasonMinimumMedianMeanMaximumMinimumMedianMeanMaximum****Bega0.000.050.373.450.000.013.5358.15Belg0.000.000.241.380.000.000.452.59Kiremt3.008.979.6316.974.91114.63133.47353.67Temperatures are in degrees centigrade (^o^C), relative humidity in percentage and amount of precipitation in millimetresTHI Temperature-humidity index

Population reaction norms

Estimates of the intercept, linear and quadratic coefficients of the population norm are presented in Table 3. Population norms for each weather variable are illustrated in Fig. 1. Table 3. Estimates and standard errors of the population norm coefficients by weather variablea_0_ (SE)a_1_ (SE)a_2_ (SE)Weekly meanMaximum daily temperature10.25 (0.44)7.46 (0.41)-0.27 (0.28)Minimum daily temperature13.75 (0.41)1.60 (0.28)-0.61 (0.21)Average daily temperature12.51 (0.45)1.84 (0.25)1.05 (0.30)Average daily humidity12.10 (0.43)-2.92 (0.40)-5.10 (0.35)Average daily THI13.27 (0.46)1.78 (0.24)-0.31 (0.34)Average daily precipitation8.24 (1.04)-5.92 (0.61)0.51 (0.43)Weekly varianceMaximum daily temperature14.01 (0.41)-1.17 (0.28)-2.87 (0.22)Minimum daily temperature14.45 (0.40)-1.19 (0.25)3.39 (0.21)Average daily temperature14.58 (0.42)-0.91 (0.29)-1.38 (0.25)Average daily humidity14.17 (0.38)-1.69 (0.25)1.01 (0.21)Average daily THI14.74 (0.43)-1.05 (0.29)-1.44 (0.26)Average daily precipitation13.26 (0.59)-2.78 (0.29)-1.39 (0.29)a_o_, a_1_, a_2_ refer to the intercept, linear, and quadratic terms, respectivelyFig. 1Population reaction norms at the phenotypic level for each weather variable. Air temperatures are in degrees centigrade (^o^C), relative humidity in percentage, and amount of precipitation in millimetres (mm); growth is in grams per day

All linear coefficients and most quadratic coefficients were significantly different from zero, thereby reflecting a significant effect of the corresponding weather variable on chicken growth. Maximum air temperature, average precipitation and average relative humidity had the strongest effect on the overall chicken growth profile change.

Individual deviation reaction norms

Table 4 presents summary statistics of the intercept, linear, and quadratic coefficients of the individual deviation norms. The arithmetic means of these coefficients are not included as they were practically zero, since they represented the individual chicken deviation from the population norm. Table 4. Estimates and summary statistics of the coefficients for the individual deviation norms by weather variablea0a1a2MSEMσMSEMσMSEMσWeekly meanMaximum daily temperature0.0580.0391.564-0.0010.0060.230-0.0240.0110.423Minimum daily temperature0.0610.0391.541-0.0150.0080.3220.0040.0040.159Average daily temperature0.0500.0381.526-0.0020.0010.042-0.0600.0160.648Average daily humidity0.0450.0341.366-0.0640.0210.847-0.0030.0070.273Average daily THI0.0430.0401.5820.0040.0010.055-0.0380.0090.377Average daily precipitation0.0540.0381.511-0.0130.0110.456-0.0130.0080.338Weekly variance****Maximum daily temperature0.0610.0431.7000.0050.0030.114-0.0020.0020.079Minimum daily temperature0.0570.0411.637-0.0020.0020.0750.0170.0090.348Average daily temperature0.0620.0431.7050.0040.0030.100-0.0070.0050.205Average daily humidity0.0480.0401.594-0.001< 0.0010.017< 0.001< 0.001< 0.001Average daily THI0.0620.0411.6520.0020.0010.044-0.0060.0040.170Average daily precipitation0.0400.0361.443-0.0070.0050.214< 0.0010.0010.021a_o_, a_1_, a_2_ are intercept, linear, and quadratic terms, respectively; summary statistics are median (M), standard error of the mean (SEM) and standard deviation (σ) calculated across all individuals

Correlations of the intercept with the linear and quadratic coefficients are shown in Table 5. These estimates were significantly different from zero (P < 0.05). Table 5. Product-moment correlations and their standard errors between the coefficients of the individual deviation norms by weather variablera0_a1** (SE)r_a0_a2_ (SE)Weekly meanMaximum daily temperature**0.79 (0.01)-0.96 (< 0.02)Minimum daily temperature-0.98 (< 0.01)0.98 (< 0.01)Average daily temperature-0.25 (0.02)0.06 (0.03)Average daily humidity-0.51 (0.02)-0.87 (0.01)Average daily THI-0.72 (0.01)0.62 (0.02)Average daily precipitation-1.00 (< 0.01)-1.00 (< 0.01)**Weekly varianceMaximum daily temperature**1.00 (< 0.01)-1.00 (< 0.01)Minimum daily temperature0.52 (0.02)0.97 (0.01)Average daily temperature1.00 (< 0.01)-1.00 (< 0.01)Average daily humidity-0.16 (0.03)-0.84 (0.01)Average daily THI0.98 (< 0.01)-1.00 (< 0.01)Average daily precipitation-1.00 (< 0.01)-0.86 (0.01)r_a0_a1_ and r_a0_a2_ correspond to the correlations between the intercept and the linear term and between the intercept and the quadratic term, respectively

Resilience phenotypes

Table 6 presents summary statistics for the chicken growth resilience phenotypes to each weather variable developed in the present study. Table 6. Summary statistics for the growth resilience phenotypes by weather variableP****˄****˅Mµ (SE)σCV%Weekly meanMaximum daily temperature22 (^o^C)1.522.822.212.20 (0.01)0.209.226 (^o^C)1.232.791.991.99 (0.01)0.2311.4Minimum daily temperature8 (^o^C)0.331.670.950.95 (0.01)0.2121.911 (^o^C)0.120.590.340.34 (< 0.01)0.0618.3Average daily temperature15 (^o^C)-1.973.550.300.26 (0.01)0.47182.717 (^o^C)-2.103.831.401.45 (0.01)0.5034.8Average daily humidity50 (%)0.160.480.310.31 (0.001)0.0414.270 (%)-0.43-0.14-0.29-0.29 (0.001)0.0413.5Average daily THI15.5-0.546.222.502.44 (0.02)0.6225.416.5-2.073.501.051.10 (0.01)0.5146.6Average daily precipitation3 (mm)-0.89-0.67-0.79-0.79 (0.01)0.033.313 (mm)-0.96-0.23-0.62-0.61 (0.01)0.1219.0Weekly variance****Maximum daily temperature1.18 (^o^C)1.261.931.631.62 (< 0.01)0.106.34.18 (^o^C)-2.59-2.47-2.53-2.53 (< 0.01)0.010.5Minimum daily temperature1.16 (^o^C)-3.87-2.40-3.19-3.18 (< 0.01)0.247.43.16 (^o^C)0.842.961.971.96 (< 0.01)0.2914.8Average daily temperature0.53 (^o^C)-0.062.431.271.26 (0.01)0.4031.61.53 (^o^C)-1.89-1.49-1.71-1.70 (< 0.01)0.063.6Average daily humidity35 (%)-0.06-0.06-0.06-0.06 (< 0.01) < 0.010.695 (%)0.020.030.020.02 (< 0.01) < 0.013.1Average daily THI0.4-0.062.401.241.23 (0.01)0.4032.21.2-7.13-4.19-5.75-5.76 (0.01)0.468.1Average daily precipitation100 (mm)-0.01-0.00-0.01-0.01 (< 0.01) < 0.0120.3250 (mm)-0.03-0.02-0.03-0.03 (< 0.01) < 0.015.3Resilience phenotypes pertain to specific points (P) on the population norm curve of the corresponding weather variable; summary statistics include minimum (˄), maximum (˅), Median (M), mean (μ) and standard error (SE), standard deviation (σ), and coefficient of variation (CV, expressed in percentage)

The resilience phenotypes associated with the weekly means of average daily temperature and THI change exhibited the highest variability, manifested in the corresponding coefficient of variation estimates. The lowest variation was mainly observed for resilience phenotypes associated with weekly variances of the weather variables.

Genomic analyses of resilience phenotypes

The variance component and heritability estimates of the growth resilience phenotypes are shown in Table 7. Table 7. Variance component and heritability estimates and standard errors (SE) for the growth resilience phenotypes by weather variableP****VGSE(VG)****VESE(VE)h^2^SE(h^2^)P-valueWeekly meanMaximum daily temperature22 (^o^C)0.0030.0030.0380.0030.070.070.1626 (^o^C)< 0.0010.0030.0510.0040.010.070.45Minimum daily temperature8 (^o^C)0.0040.0030.0400.0030.090.070.1111 (^o^C)< 0.001< 0.0010.004< 0.0010.110.070.05Average daily temperature15 (^o^C)0.0320.0150.1810.0140.150.070.0117 (^o^C)0.0370.0170.2100.0170.150.070.01Average daily humidity50 (%)< 0.001< 0.0010.002< 0.001< 0.010.060.5070 (%)< 0.001< 0.0010.002< 0.001< 0.010.070.50Average daily THI15.50.0770.0280.3000.0260.210.07< 0.0116.50.0550.0190.2030.0180.210.07< 0.01Average daily precipitation3 (mm)< 0.001< 0.0010.001< 0.0010.050.070.2113 (mm)0.001< 0.0010.0130.0010.080.070.11**Weekly variance*Maximum daily temperature1.18 (^o^C)0.001< 0.0010.0100.0010.100.070.094.18 (^o^C)< 0.001< 0.001< 0.001< 0.001< 0.010.070.50Minimum daily temperature1.16 (^o^C)0.0040.0040.0520.0040.070.070.163.16 (^o^C)0.0060.0060.0780.0030.070.070.16Average daily Temperature0.53 (^o^C)0.0150.0110.1440.0110.0970.070.091.53 (^o^C)< 0.001< 0.0010.004< 0.0010.0670.070.17Average daily humidity35 (%)< 0.001< 0.001< 0.001< 0.0010.0570.070.2295 (%)< 0.001< 0.001< 0.001< 0.0010.0070.070.50Average daily THI0.40.0140.0110.1430.0110.0970.070.091.20.0170.0150.1990.0150.080.070.12Average daily precipitation100 (mm)< 0.001< 0.001< 0.001< 0.0010.0970.070.09250 (mm)< 0.001< 0.001< 0.001< 0.0010.0870.070.12Resilience phenotypes pertain to specific points (P) on the population norm curve of the corresponding weather variable; parameters include the genomic variance (V_G_), residual variance (V_E_), heritability (h^2^) and corresponding standard errors (SE); P-values refer to heritability estimates. Significant P-values (P < 0.05) from the Likelihood Ratio Test are indicated with *

Statistically significant genomic variance estimates (P < 0.05) were derived for five growth resilience phenotypes pertaining to the weekly means of average daily air temperature, minimum daily temperature, and THI. Corresponding heritability estimates for these traits were moderately low, ranging from 0.11 to 0.21. The other resilience phenotypes did not exhibit significant genomic variation (P > 0.05).

No inflation was detected in the ensuing genome-wide association study, with the lambda factor values ranging from 0.998 to 1.003. No genome-wide significant SNPs were detected for any of the studied phenotypes. However, three genome-wide suggestive SNPs were found for the growth resilience phenotypes related to weekly mean THI and another two suggestive SNPs were associated with response to the weekly mean of minimum daily temperature (Table 8). In all cases, multiple trailing SNPs supported the corresponding hits relatively well (Supplementary Files 4, 5 and 6). Two of the former SNPs were common to both THI resilience phenotypes at the two points of the reaction norm. The proportion of variance explained by each SNP is included in Table 8. Table 8. Suggestive single nucleotide polymorphisms detected for growth resilience phenotypesPChrPosition (in bp)Beta effect and standard errorMAFP-valuePVE (%)Weekly meanMinimum temperature11 (^o^C)423,738,2160.013 ± 0.0020.4681.2E-71.7511 (^o^C)3107,013,0670.039 ± 0.0070.0251.5E-71.72Average THI15.5343,321,960-0.132 ± 0.0250.3841.5E-71.7315.514,752,2060.393 ± 0.0750.0221.9E-71.7016.5343,321,9600.109 ± 0.0210.3841.6E-71.6716.514,752,206-0.320 ± 0.0620.0222.7E-71.6516.5239,779,7280.290 ± 0.0560.0282.9E-71.66Resilience phenotypes pertain to specific points (P) on the population norm curve of the corresponding weather variable. Chr, MAF and PVE correspond to the chromosome number, minor allele frequency and the percentage of phenotypic variance explained by the marker, respectively

Discussion

The present study aimed to derive and examine individual growth resilience phenotypes to weather changes of commercial chickens reared in smallholder farm conditions in Ethiopia and to investigate the genomic architecture of these new phenotypic traits.

Weather variable profiles

The venue of the study was the ILRI poultry research facility in Addis Ababa, Ethiopia, which is located in an elevated region at 2,382 m above sea level with a weather profile reflective of the Ethiopian highlands. The study took place across three distinct local seasons: Belg, Bega and Kiremt. The Belg local season (February to May) was characterised by the highest maximum daily temperatures, often exceeding 25°C but rarely reaching 30°C. This season also presented the lowest average daily humidity, with values commonly in the range of 40% to 60%. The Bega local season (October to January) featured the lowest minimum daily temperatures, with values usually under 10°C, and a relative humidity only slightly higher than the one observed in Belg. Both Belg and Bega were generally dry seasons with very little precipitation in the form of rain. On the contrary, the Kiremt local season (June to September) was characterised by high daily relative humidity and considerable precipitation levels. However, air temperatures were milder than in the other two seasons.

While the prevailing air temperatures in the region can be considered moderate, the relatively short range of the chicken thermo-neutrality zone (18–24°C) [28, 29] implies the possibility of observing effects of both heat stress and, perhaps more likely, cold stress in the chickens of study. The impact of relative humidity is not well known, with most studies indicating that values ranging from 25 to 75% play a minor role, if any, unless accompanied by high temperatures [30, 31]. In the present study, the latter was captured in the calculated THI. Similarly, strong precipitation can be related to an increase in the transmission of parasitic diseases such as coccidiosis, thereby affecting chicken growth [32].

Reaction norms and growth resilience phenotypes to weather change

Intercepts and linear coefficients of the population reaction norms were always significantly different from zero (P < 0.05). The signs of the linear coefficients meant that chicken growth generally benefitted from increasing air temperature and THI and decreasing humidity and precipitation. The negative signs of the linear coefficients associated with the weekly variance of the weather variables suggested that growth was challenged by increased weather instability. The quadratic coefficients were also statistically significant (P < 0.05) except for the weekly mean of maximum air temperature, THI, and amount of precipitation. When significant, this coefficient revealed a non-linear reaction of average chicken growth to changing weather events. In case of the latter, the direct numerical interpretation of some of these results can be challenging. A linear coefficient represents the expected change in growth per unit change in a weather variable, and in the absence of a quadratic term, this rate of change remains constant across the entire range of weather values. However, a significant quadratic coefficient indicates the convex or concave nature of the relationship, suggesting that the rate of change is not uniform across the range of values. For example, a 1°C increase in average daily temperature was associated with an increase in average daily growth of 1.84 g within the studied range of 14–18°C. The significant quadratic coefficient of 1.05, however, indicated a convex relationship, where growth only increased when temperatures exceeded 15.5°C. In contrast, the interpretation of the effects of THI and precipitation were more straightforward, as the quadratic terms were not significant. For these variables, a unit increase within their mean value ranges (15–17 and 0–17 mm, respectively) resulted in a daily growth increase of 1.78 g for THI and a decrease of 5.92 g for precipitation.

For the individual deviation norms, the linear coefficient on average relative humidity and the quadratic coefficient on average air temperature exhibited the highest variance, suggesting that individual chickens may respond differently to these weather effects. The variation of the linear and quadratic coefficients on the weekly variances of weather events was low, implying that most chickens may respond similarly to this challenge.

Exploring the correlations of the intercept with the linear and quadratic coefficients of the individual deviation norms provided an insight into the possible phenotypic relationship between the actual level of chicken growth (intercept) and changes in growth profile in response to weather fluctuations (linear and quadratic coefficients). In the present study, all these phenotypic correlations were significantly different from zero (P < 0.05). For the weekly means of weather variables, negative linear and positive quadratic correlations were found for minimum and average air temperature and THI, suggesting that fast growing chickens may be relatively resilient at low values of these variables but quickly become adversely affected as values increase. This result complements the positive linear and negative quadratic correlation observed for maximum air temperature, implying that fast growing chickens may benefit from low values but the benefit diminishes quickly as air temperature increases. Conversely, slow growing chickens seem to be more resilient to high temperatures, but also more vulnerable to the cold. These results are in line with the expectation that small chickens are more susceptible to stress due to low air temperature, as their high surface area compared to volume may lead to faster rates of heat loss [33]. Similarly, fast growing chickens are expected to be more susceptible to heat stress [34], but the range of maximum air temperatures observed in the present study was moderate and only slightly over the thermo-neutral zone of chickens. Thus, although we observe a diminishing benefit of high temperatures for fast growing chickens, we still do not see the full adverse effect of heat stress on growth.

Negative linear and negative quadratic correlations for relative humidity and amount of precipitation indicate that fast growing chickens are likely to be more affected by increasing values of these weather variables, implying higher humidity and heavy rain, than slowly growing chickens.

The correlation estimates pertaining to weekly variances of weather variables collectively suggested that fast growing chickens might generally be less susceptible to variation in daily air temperature, but more susceptible to variation in humidity and precipitation. We speculate that, while larger chickens may tolerate weekly weather variation, particularly during short-term cold stress events, they face potential difficulties when looking for shelter under sporadic heavy precipitation. Also, reduced feed intake while sheltering and not scavenging may have a greater impact on larger birds with higher metabolic demand.

Notably, the associations between chicken growth and resilience manifested in all these correlations are at the phenotypic level, and future research should investigate their genetic component.

Novel phenotypes of chicken growth resilience were derived here based on the reaction norm coefficients, reflecting changes in individual chicken growth profiles in response to changing weather conditions. Large variations in resilience would suggest that individual chickens react differently to this environmental challenge. We observed the largest phenotypic variation in growth resilience in response to average air temperature and THI changes. Interestingly, these weather variables exhibited relatively small fluctuations across the timeline of the study, manifested in the low coefficient values in the respective population norm. As such, it is important to clarify that a significant population norm coefficient may indicate an overall strong effect of the weather variable on growth, but it does not necessarily imply variation in the individual resilience phenotypes. Similarly, there is a distinction between variation in the coefficients of the individual deviation norms and variation in the respective resilience phenotypes, as the former pertains to the entire curve of the norm while the latter is estimated at specific points of the curve.

Genomic analyses of growth resilience phenotypes

Significant genomic variance and heritability estimates were derived for growth resilience to changes in average air temperature, minimum air temperature and THI. These estimates ranged from 0.11 to 0.21 and are considered moderately low, in line with previous studies on fitness-related traits in chickens [35] and other species [36, 37]. Despite the low heritability, the presence of genomic variation among individual chickens supports the possibility of selective breeding to enhance growth resilience within future breeding programmes [38].

Genome-wide association analyses revealed five SNP markers that were suggestively associated with the growth resilience phenotypes that exhibited significant genomic variance.

The SNP identified on chromosome 4 in association with growth resilience to changes in minimum air temperature was located within the region of gene ENSGALG00000047398, a long non-coding RNA (lncRNA). In recent years, several studies have identified lncRNAs as major players in the regulation of gene expression in biological networks associated with diverse processes [39], including metabolic response to heat and cold stress in both plants [40, 41] and animals [42, 43]. These processes also include the regulation of adipose tissue development and activity [44].

The SNP on chromosome 3, also associated with growth resilience to changes in minimum air temperature, was not found within any coding region but was in linkage disequilibrium (within 200 Kb) with the PINX1 gene, which codifies for a telomerase inhibitor. Telomerase inhibitors prevent telomere elongation and help maintain telomere homeostasis [45]. Air temperature stress is known to affect telomeric integrity [46] and previous studies have associated leucocyte telomere attrition rates with weather conditions and productivity in other farm animal species [47]. It is plausible that telomerase inhibitors play a significant role in the resilience of animal production to weather change through an effect on telomere integrity and cell senescence.

The SNP on chromosome 3 associated with THI change is located within the coding region of the phosphodiesterase 10A gene (PDE10A). Phosphodiesterases are involved in the regulation of cyclic AMP and cyclic GMP in multiple processes [48]. Previous studies linked phosphodiesterase 10A with osteogenic differentiation [49], lung inflammation [50] and vascular remodelling [51]. Furthermore, it has been proposed that the inhibition of the enzyme may protect the individual from diet-induced obesity and insulin resistance [52]. Interestingly, studies in mice [53, 54] implicated PDE10A in thermoregulation and obesity: inhibition of phosphodiesterase 10A appeared to stimulate thermogenic gene expression and reduce hedonic feeding. Thus, PDE10A could possibly be related to energy homeostasis in response to weather change.

The SNP on chromosome 1 associated with THI change was found within the coding region for the inter-alpha-trypsin inhibitor heavy chain (ITIH5). This protein belongs to a family of plasma inhibitors that have an important role in inflammation and carcinogenesis [55]. ITIH5 reportedly inhibits the growth and metastasis of melanoma in humans [56]. A recent study [57] linked ITIH5 with the regulation of adipose tissue homeostasis, which may suggest a role in the coping mechanism to environmental stress, particularly considering that the non-shivering thermogenesis of the brown adipose tissue may be activated by cold stimulation [58]. Furthermore, genes involved in lipid metabolism are reportedly highly enriched under heat stress in the avian liver [59], where most of the lipogenesis occur [60]. These genes may also contribute to chicken response to varying weather conditions. Finally, ITIH5 is a p-53 responsive gene [56], and heat stress can induce apoptosis by activating p-53 mediated mitochondrial pathways [61]. This indicates that ITIH5 may have a role in cell response to weather stressors.

The SNP on chromosome 2, associated with growth resilience to changes in THI, is located within the coding region for TGFβR2 (transforming growth factor beta receptor 2). A differential expression of this growth factor was previously identified in a study related to heat stress of Holstein cows and apoptosis processes in fibroblasts [62]. Furthermore, as previously mentioned, heat stress can induce apoptosis by activating p-53 mediated pathways, and the inactivation of p-53 is also known to distort TGFβ signalling and promote tumour growth and malignancy [63].

Noticeably, the SNPs on chromosomes 1 and 3 affected resilience to THI change at both sides of the reaction norm and the effect had the opposite sign in the two cases. This is not necessarily surprising, as response to stress is highly dependent on gene expression, with different patterns activating and deactivating metabolic networks [64].

The proportion of variance explained by each SNP was between 1.65% and 1.73%, within the range of previous reports on other chicken traits [65, 66], although our estimates might be slightly inflated due to the Beavis effect [67]. This result, together with the overall illustration of all individual SNP effects (Supplementary Files 4, 5 and 6) suggest a largely polygenic architecture for the studied resilience phenotypes. This is typical for fitness-related traits [68], where multiple genes, each contributing a relatively modest effect, collectively control the genetic background and influence the resilience phenotypes under investigation. Selective breeding based on estimated genomic breeding values is the best approach to genetic enhancement of chicken performance for such traits [69]. Nevertheless, results on significant genetic variants may be incorporated into the genomic evaluation process to increase the accuracy of the ensuing selection [70].

Finally, research is warranted on the biological mechanisms through which weather fluctuations influence growth, including heat stress, immunosuppression, and reduced feed intake, as well as the potential interactions among them. Given their complexity, chicken resilience to weather change becomes a truly multifactorial trait, beyond the polygenic architecture reported in the present study. Thus, a deeper understanding of the underlying biological mechanisms is crucial for improving the accuracy of predictive models and further refining the discovery of relevant genomic variants. Other environmental factors, particularly those experienced during incubation and early chick development, may also play a role in shaping future thermal resilience. While batch effects in the present study accounted for these conditions to some extent, there may still be unexplained variation affecting the physiological response to environmental stress. Furthermore, seasonal variations in disease prevalence add yet another layer of complexity, affecting chicken performance independently of weather fluctuations. Recognising these multifaceted influences is essential for a holistic analytical approach towards understanding the specific mechanisms of environmental resilience.

Conclusions

To our knowledge, this is the first phenotypic and genetic study of growth resilience in commercial chickens in Ethiopia. We have shown that random regression models can be used to derive resilience phenotypes describing changes in the chicken growth profiles in response to varying weather conditions. The significant genomic variance and candidate genes reported here may be used to inform breeding programmes aimed at addressing climate change and sustainably improving chicken production under ever-evolving environmental challenges. Future studies should investigate the genetic correlation of the new growth resilience phenotypes with other chicken traits in the breeding goals. Furthermore, additional studies should focus on other chicken ecotypes in sub-Saharan Africa, including local breeds, to reveal possible distinctive adaptation patterns.

Supplementary Information

Supplementary Material 1. Supplementary Material 2. Supplementary Material 3. Supplementary Material 4. Supplementary Material 5. Supplementary Material 6.