The introduction of single nucleotide polymorphism (SNP) arrays have significantly advanced the genetic analysis of economically important traits in chickens, offering unprecedented insights into their complex genetic architecture. In genome-wide association studies (GWAS), SNP arrays serve as powerful tools for identifying novel loci associated with phenotypic variation across different chicken lines, including egg production and quality [1, 2], growth [3, 4], meat quality [5], and disease resistance [6]. These advancements have identified quantitative trait loci (QTLs) and candidate genes linked to key traits of laying chickens, such as egg production and quality, and body weight, that can drive genetic improvements in poultry breeding [7–13].

Despite the achievements, challenges remain. Although SNP arrays are a cost-effective and widely used genotyping tool, their medium density limits the resolution for fine-mapping causal variants. One way to overcome this limitation is through imputation, where lower density SNP arrays are imputed to whole-genome sequence (WGS) level using reference panels of related sequenced animals. In chickens, several studies have incorporated imputation in GWAS to enhance detection power [14–16]. This approach provides denser genomic information, and has been shown to provide a substantial boost in GWAS power from a simulation study, enabling the performance of commercial SNP arrays close to an ideal “complete” chip [17]. Other researchers suggested that WGS-based GWAS enables the discovery of novel variants that were missed by array-based GWAS, but the gains in detection power are modest since they often reinforce signals from known QTLs [18–20]. Nevertheless, there is broad consensus that increasing sample size remains one of the most effective strategies for improving detection power in GWAS [17, 20–22].

Another fundamental limitation for GWAS resolution is that a detected QTL contains numerous genes rather than pinpointing candidate genes with precision, due to the extent of linkage disequilibrium (LD) and generally long haplotypes in structured populations. However, studies have demonstrated that repeated intercrossing of outbred lines can enhance QTL mapping resolution for growth traits such as body weight [23–26], which suggests that introducing recombination events help reduce QTL intervals. GWAS methodologies can also be adapted to improve resolution and interpretation. Including covariates in GWAS is a common approach to control for confounding factors such as population structure and environmental influences that could otherwise bias the identification of the true genetic associations. Covariates such as principal components, age, sex, batch effect, leading SNPs, and correlated traits were generally considered to avoid stratification [16, 26–29]. Beyond just addressing confounding, multi-trait GWAS approaches have been increasingly explored, allowing to detect pleiotropic loci that may not reach the significance thresholds in single-trait GWAS analyses [20, 30, 31].

While traditional GWAS is effective for detecting loci with strong effects, it faces limitations in detection of polygenic traits where many SNPs contribute small individual effects, or when LD complicates the identification of causal variants. Bayesian methods offer a solution by incorporating prior probabilities and modelling sparse SNP effects, thereby refining the large QTLs identified through GWAS. A sequential approach using GWAS to identify QTLs, followed by a Bayes-based approach for fine-mapping and integration of functional enrichment, has successfully revealed candidate genes for important traits in cattle [32–34].

The objective of this study was to dissect the genetic architecture underlying production traits through integration of different GWAS methodologies and Bayesian fine-mapping tool, focused on a backcrossing and intercrossing scheme between Araucana chickens and high-performing White Layer (WL) hens.

Animals and Phenotypes

The animals used in this study are owned by the breeding company Lohmann Breeders GmbH. This study focused on a backcrossing scheme between distinct breeds, developed as part of the EU Horizon 2020 project Innovative Management of Animal Genetic Resources (IMAGE; https://www.imageh2020.eu/). This backcrossing scheme was designed to introgress blue eggshell colour from Araucana chickens into a high-performing White Layer line. The blue eggshell colour is determined by the autosomal inheritance of the SLCO1B3 gene (1:65167287–65185448) [35], “BB”, “Bw” and “ww” refer to the genotypes of the blue eggshell insertion where “B” refers to the color allele. Six Araucana males were crossed with ten high-performing White Layer (WL) females from a White Leghorn line to produce F1, followed by two back-crossing generations with WL hens to produce BC1 and BC2, and intercrossing within BC2 to produce the IC generation (Fig. 1). Only males of F1 and BC1 were selected to mate with white layer hens for producing BC1 and BC2. Parents from F1, BC1 and BC2 were selected as the carrier of the color allele (Bw).

Phenotypic traits were recorded in hens from BC1, BC2, IC and the WL founders of the BC cohorts (FounderBC1, FounderBC2), data is available in Supplementary Data 1. Only eggs with intact shell were selected for measuring egg quality traits. Egg number (EN) was recorded in thirteen 28-day periods, starting from 20 weeks (w) of age as the first period (EN1; 20w to 24w), measured as the total number of eggs during the corresponding period. There was no record for EN9 where the selection of individuals took place. Egg weight (EW) was measured at 30, 40, 50 and 70 weeks of age as the average weight of a 3-day clutch, although “FounderBC1” and “BC2” did not have data for EW50. Body weight (BW) was a onetime measurement at 32 weeks of age (Table 1).

Genotyping and Array-based GWAS

Araucana founders (6 males) were sequenced and WL founders (10 females) were genotyped using Affymetrix Axiom 600k array [36]. Remaining animals were genotyped for 52k SNPs using a custom Affymetrix array (ThermoFisher, Santa Clara, CA, USA). Overlapping SNPs were combined and genotypes were mapped to the chicken reference genome GRCg6a. Genotype data were filtered using PLINK 1.9 [37], by removing samples and SNPs with missing call rate greater than 15 %. The iltered genotype data contained 46,605 SNPs and 1,639 individuals with 1,024 females and 615 males.

Genome-wide association studies (GWAS) for EN, EW and BW measured throughout different periods were performed based on the filtered genotype data and different sample groups, including all the animals with non-missing phenotypic records (WL founders for BC1 and BC2, BC1, BC2, IC) denoted as group “all”, and group “offspring” included BC1, BC2 and IC, and by each generation. BW32 was included as a quantitative covariate in the analyses of EW traits, and the corresponding EW measurement was included as a covariate in the analyses of BW32.

Array-based GWAS was conducted for each trait including the generation as a fixed effect, using single-marker regression with polygenic correction, and the genetic relationship matrix (GRM) for each GWAS was computed between all the animals (1,639) using array genotypes, both by GCTA (version v1.94.1) [38–40]. Additional minor allele frequency (MAF) filtration was applied on each GWAS to remove SNPs with MAF lower than 0.01. The filtered genotyped data was fitted to the mixed linear model:

where is a vector of phenotypes of animals; is a vector of fixed effects; is the design matrix for fixed effects; is a vector of random genetic effects that captures the total genetic variation; is a vector of residuals. We applied a strict genome-wide Bonferroni-corrected threshold, calculated as divided by the total number of tested SNPs for each analysis [41].

Imputation and accuracy

The filtered genotype data were phased and imputed using Beagle 5.4 [42, 43] with a reference panel of nine Araucana (of which six were the founders for the IMAGE backcross) and 161 White Leghorn samples and commercial White Layer lines, retrieved from a global reference set of around 5k chickens, prepared by the Chicken Genomic Diversity Consortium [44]. Since the Genotype data was mapped to GRCg6a, the reference sequence was lifted from GRCg7b to GRCg6a by bcftools-1.20 [45] using a custom chain file prepared according to [46]. The snakemake pipeline for the liftover is available on https://doi.org/10.5281/zenodo.15343279.

The imputation produced nearly 11 million biallelic variants. Imputation accuracy of the filtered data was estimated in 5-fold cross validations with two repetitions, by masking variants partially from the array-genotyped data in each fold and compared with the imputed SNPs correspondingly. The mean concordance and Pearson’s correlation were 0.95 and 0.92, respectively. The Imputation Quality Score (IQS), which adjusts for the minor allele frequency [47], was 0.88 on average. For the following analyses, we considered 4,868,808 (44.4 %) ariants having good predicted imputation quality (dosage r-square; DR2 > = 0.8).

Imputed-based GWAS

GWAS based on the 4,868,808 imputed variants was conducted for each trait including the generation as a fixed effect, using a fast mixed linear model based tool in GCTA [48]. First, a sparse GRM was generated between all the animals (1639) using the imputed genotypes with a cutoff value of 0.05, then the GWAS was performed using the sparse GRM previously generated for each trait and different sample groups with non-missing phenotypic records. Additional minor allele frequency (MAF) filtration was applied on both analyses to remove SNPs with MAF lower than 0.01. The genome-wide Bonferroni-corrected threshold became stricter (0.05 divided by the total number of tested SNPs for each analysis) as the number of SNPs significantly increased. Additionally, BW32 was included as a quantitative covariate in the analyses of EW traits, and vice versa.

Multi-trait metaGWAS

A multi-trait metaGWAS was conducted on all EN measurements (EN1 – EN13, excluding EN9 due to missing data) using the summary statistics from imputation-based GWAS. The formula suggested by Bolormaa et al. [30], where each SNP was assigned a t-value calculated as the allele effect divided by the standard error obtained from the imputed-based GWAS , while is an vector of t-values of for the traits, is a transpose of the , is an inverse of the correlation matrix between all the traits calculated based on the signed t-values of the SNPs. The genome-wide Bonferroni-corrected threshold was applied to identify significant variants.

Phenotypic correlations and genetic parameters

Phenotypic correlations between traits after adjusted by generation means, were calculated as Pearson’s correlations by R (version 4.3.2) [49] using function cor with pairwise.complete.obs, the corresponding p-values for the phenotypic correlations were calculated using cor.mtest from R package “corrplot” (version 0.92) [50].

Genetic parameters and correlations between traits were estimated from variance captured by the imputed variants including generation as a fixed effect, with the same model as described in the Array GWAS section, using restricted maximum likelihood analysis. We used GREML and Bivariate GREML analysis implemented in GCTA [39, 51, 52] that utilised the GRM computed between all animals (1639) using the imputed genotypes with MAF filtration of 0.01 [38]. GCTA only considered individuals with non-missing phenotypes, therefore the number of samples included in the analysis varied depending on the corresponding phenotypic records.

Fine-mapping analysis

We defined the QTL for fine-mapping based on the distance between the first and the last significant SNP that passed the Bonferroni threshold identified by imputed-based GWAS in group “all”, with a 1 Mb extension both upstream and downstream. All QTL regions were limited to a maximum size of 11 Mb. QTLs of 17 traits for fine-mapping can be found in the Supplementary Data 2.

We used two Bayesian-based fine-mapping methods, the BFMAP (Version 0.65) [32] and FINEMAP (Version 1.4.2) [53]. In our analysis, we observed that FINEMAP exhibited sensitivity to minor alterations in input data that has led to inconsistency in the results, therefore FINEMAP results are not discussed further here. Meanwhile, BFMAP demonstrated a high degree of stability in credible set (CS) identification, providing robust and interpretative results for further discussion.

The phenotypes and imputed genotypes of the individuals for which both records exist, including generation as a covariate, were used as input for the following BFMAP analyses. The QTL heritability was estimated from the proportion of variance explained by SNPs in the QTL region by , based on the regional GRM computed by the QTL SNPs between all animals according to Yang’s formula [38]. The overall GRM was computed by BFMAP based on the imputed genotypes with MAF filtration of 0.01 between all animals (1639), and were included during fine-mapping together with the estimated QTL heritability. The Bayesian model developed for fine-mapping is:

where is a vector of a phenotype from n animals; is a vector of covariate effects (generation) and is the corresponding design matrix; is a vector of variant effects and is the corresponding genotype coding matrix; is a vector of polygenic effect considering the population structure; is the residual. Probability of the null model P(D|M₀) and the probability of model with variant inclusion given the data P(D|M) can be computed to obtain the scaled Bayes factor and the corresponding p-value for a SNP. BFMAP approximates flat priors for SNP effects [32].

For fine-mapping, BFMAP first performs an iterative forward selection process to identify independent signals, by sequentially adding lead SNPs that were evaluated conditional on all previously selected lead SNPs to the model. The process continues until no variants meet the software’s predefined threshold (2logBF + 1 < 2logm_eff), with each selected lead SNP represents the proxy of an independent signal. Next, each signal is expanded to include variants in high LD (default r² > 0.3) with the signal proxy. The final lead SNP for each signal is then re-evaluated and determined through repositioning within each signal set. Finally, a posterior probability of causality is assigned to each variant within a signal conditioning on the lead SNPs from all other signals while excluding the lead SNP of that signal. PPC values are then sorted in descending order, starting from top, a list of variants with the cumulative sum of their PPC over 95% consists a credible variant set (CS) for each signal [32]. The summed PPC for a gene can also be generated based on the , variants with 1 kb extension both upstream and downstream of the gene were also included [54]. All genes with summed PPC over 0.01 were documented in Supplementary Data 3. Afterwards we combined all signal-specific CS into a final CS for a QTL.

Functional enrichment of BFMAP fine-mapping results

Functional enrichment of the fine-mapping results from BFMAP was conducted using the R package “gemrich” (version 0.1.3) [54], which is designed to integrate functional annotations with the forward selection approach of BFMAP. Variant functional annotations were obtained with the Ensembl Variant Effect Predictor (VEP, version 112) [55], applied to variants included in the CS determined by BFMAP. The categories analyzed included Consequences, Impact, and Biotype. The number of variants analyzed in each category could vary slightly dependent on the availability of annotation.

Signal filtration was applied using a p-value threshold that targeted the lead variant of the signal, thereby a signal was either retained or filtered out as a whole. We applied a threshold of , as recommended by the author. All variants within retained signals were included in the following analyses.

The estimated probability of causal variants being in a category C was generated using maximum likelihood based on the produced . The enrichment factor for each category quantified how enriched the candidate causal variants in that category compared to its background frequency : . To ensure robust enrichment estimates, we excluded runs where fewer than 10 variants remained after signal filtration, or where fewer than two category types were present; these cases were flagged as “Failed” in the “mle_status” column, with the reason noted in the “reason” column. In addition, categories representing less than 1% of variants were grouped into a single “remaining” category. Parameters for enrichment estimation were documented in Supplementary Data 4.

For a total of variants remaining, we then renormalized for each variant to account for the corresponding category enrichment , weighted by the total PPC of variants:

Variants within filtered signals, along with their and , are provided in Supplementary Data 5. When computing the summed PPC for a gene based on the filtered variants’ , we included variants within the gene plus 1 kb upstream and downstream. Genes with summed PPC greater than 0.01 were reported in Supplementary Data 6.

Workflow

The pipeline of this work was implemented with Snakemake 7.32.4 [56]. The scripts, snakefiles, working conda environment definitions in yaml format, and the dependency analytics graph (DAG) and the rulegraph of the pipeline, are available on https://doi.org/10.5281/zenodo.17473280.

Phenotypical characteristics

Descriptive statistics for phenotypes in all generations are shown in Table 1. A smaller mean number of eggs laid (11.9) and a larger coefficient of variance (CV % at 76.1%) were observed during early period of EN measurement (EN1; 20 w – 23 w), while during later periods (EN2 – EN13) the mean number of eggs laid ranged from 21.7 to 25.7 and CV % ranged from 17.2% to 35%. During all EN measurements there were always some individuals (from 0.8% to 19.6% of the population) with none egg laid. From 1.8% to 39.7% of the population laid 28 eggs across the periods. The mean weight and CV % for BW at 32 weeks were 1786.2 grams and 11.6%, respectively. For EW measured at different weeks, their CV % were similar average at 7.7%, the mean egg weight was increasing from 57.7 to 64.0 grams as the hens aged.

Visualization of the phenotypic correlations (, lower triangle) and genetic correlations (, based on imputed SNPs, upper triangle) among all traits, after adjusting for generation effects, together with their estimated heritabilities, is presented in Fig. 2.

Overall, showed moderate to strong positive correlations among EN measurements, although these correlations weakened as the interval between measurements increased. In contrast, many across EN failed to converge after adjusting for generation effects, specifically from EN4 to EN6. This period coincided with peak laying performance, during which reduced variance limited reliable estimation. Interestingly, EN1 exhibited strong positive with the start and end of the laying period, and moderate positive with the mid-to-late laying phase. A similar though weaker pattern was also observed for EN2.

Generally, EN showed negative correlations with both EW and BW in our analysis, with showing stronger correlations than . For EW70, correlations with EN gradually strengthened over the laying period, while for BW32, EW30 and EW40, both and were strongest at the beginning and end of the laying period, and weakest during the mid-to-late phase.

EW measurements at different weeks of age were consistently and positively correlated, with ranging from 0.64 to 0.78 and from 0.88 to 0.96. However, estimates involving EW50 did not converge. BW32 also showed positive correlations with EW traits, with the strongest observed with EW40 (0.43) and the strongest with EW 30 (0.57).

Estimated genetic heritability was exceptionally high for EN1 that reached 0.4, and dropped sharply for EN2, EN10 and EN11, ranging from 0.12 to 0.15, and was close to zero for other EN periods. In contrast, BW32 had the highest heritability at 0.8, while EW30 to EW70 showed moderately high values, ranging from 0.57 to 0.69.

All array-based and imputed GWAS results across traits and groups are presented as Manhattan plots in Supplementary Figs. 1.

Both array and imputed GWAS indicated a polygenic architecture for EN measurements across different sample groups, with many loci each contributing small effects to EN. In general, imputed GWAS identified more significant variants than array GWAS, particularly in the “FounderBC2” and “FounderBC1” groups. During the mid-to-late laying phase (EN7 – EN11), significant variants were detected across more sample groups, most notably on chr 2 and 12 (Fig. 3).

For EW30, both array and imputed GWAS revealed significant peaks on chr 1 and 4 in the "all" and "offspring" groups, with an additional peak on chr 2 detected in the “offspring” group which just below the significance threshold in the “all” group. Interestingly, in the “IC” group, a significant signal on chr 4 was detected in the imputed GWAS, which applied a stricter significance threshold, while the same signal in the array GWAS fell just below significance. Conversely, array GWAS identified significant signals on chr 4 in the "BC1" and "BC2" groups that did not surpass the more stringent significance threshold in the imputed analysis. Emerging signals on chr 1 and 2 were also observed in the “BC2” group in both GWAS approaches.

Across EW40, EW50, and EW70, both array and imputed GWAS consistently identified a significant peak on chr 4 that persisted across all EW measurements in the "all" and "offspring" groups, as well as in the "IC" group for EW40 and EW50, and partially in the “BC1” group. The significant peak on chr 1 detected for EW30 was also observed for EW40 in the "all" and "offspring" groups in both GWAS analyses, and in the "BC2" group based on imputed GWAS. Notably, for EW70 in the “BC1” and “BC2” groups, additional suggestive peaks were observed on chr 1 at a distinct position and on chr 7, respectively, both falling just below the genome-wide significance threshold in the two GWAS analyses.

For BW32, three major significant peaks on chr 1, 4 and 27 were detected by both array and imputed GWAS in the "all," "offspring," and "IC" groups, while peaks on chr 1 and 4 were identified in the "BC2" group. In the "BC1" group, array GWAS identified significant signals on chr 1 and 4, but only the chr 4 signal remained significant in the imputed GWAS. No significant variants were detected in the WL founders, although the chr 4 signal nearly reached significance in the "FounderBC2" group, and a spurious signal on chr 1 emerged in the “Founder_BC1” group, both from the imputed GWAS. Notably, the major peaks on chr 1 and 4 identified for BW32 overlapped with those associated with multiple EW measurements (Fig. 4).

Multi-trait metaGWAS for EN

Given the moderate to strong correlations among EN measurements (Fig. 2), and the polygenic background observed for EN in both array-based and imputed GWAS, we conducted a multi-trait metaGWAS to increase statistical power and identify pleiotropic loci influencing EN across the laying period. The metaGWAS identified several significant variants on chr 2 associated with EN traits from EN1 to EN13 (Fig. 5).

Including correlated traits as covariates in GWAS

Given the high correlations between BW32 and EW measurements, we performed both array and imputed GWAS for EW traits while including BW32 as a quantitative covariate, and vice versa. This approach provided additional insights into the genetic architecture of these traits.

When BW32 was included as a covariate in EW GWAS, signals that previously overlapped between BW32 and EW40, EW50, and EW70 on chr 1 and 4 were suppressed and no longer reached significance across any sample group in either GWAS approach. For EW30, covariate adjustment helped clarified the presence of spurious associations without markedly altering most p-values. For example, a signal on chr 2 became significant in the “all” group and remained significant in the “offspring” group in array GWAS. Meanwhile, the overlapping signal on chr 1 was strongly suppressed and disappeared across sample groups, whereas the signal on chr 4 showed greatly increased p-values and remained close to the significance threshold only in the “all” and “offspring” groups. (Fig. 4). Interestingly, for EW70, the previously suggestive peak on chr 1 at a distinct position in the “BC1” group was retained with comparable p-values, and additional suggestive peak was observed on chr 21 in the imputed GWAS, although both remained just below the genome-wide significance threshold.

Conversely, including EW measurements as covariates in the GWAS for BW32 generally increased p-values for the major signals on chr 1 and 4 (Fig. 4). The locus on chr 27 also showed increased p-values when EW50 and EW70 were included. Such reduction in statistical power sometimes preventing these signals from reaching genome-wide significance. Nevertheless, the three major peaks on chr 1, 4, and 27 were largely retained in the "all" and "offspring" groups, except for the chr 1 peak, which dropped below the significance threshold when EW50 was used as the covariate. Across different EW covariates, the chr 1 signal was suppressed in the “IC” group but retained in “BC2”, while the chr 4 signal was sometimes suppressed in the “IC” and “BC2” groups and completely lost in the “BC1” group.

All array and imputed GWAS results including covariate among groups between BW32 and EW30, EW40, EW50, EW70 are available as Manhattan plots in Supplementary Figs. 2.

Fine-mapping GWAS-identified QTLs and enrichment analysis

Fine-mapping analysis was conducted to refine the QTL regions identified through imputation-based GWAS, with the aim of pinpointing potential causal variants and candidate genes. EN traits exhibited polygenic signals, with only a few variants reaching significance and no clear peaks, making these results unreliable for QTL identification. Given the moderate to strong correlation among EN measurements, a multi-trait metaGWAS was performed to detect pleiotropic loci associated with EN1 – EN13. Fine-mapping was subsequently applied to the significant loci identified in this meta-analysis. The fine-mapping results—including the QTLs analysed for 17 traits, their estimated heritabilities, enrichment factors, and the probability of including causal variants within functional annotated categories, as well as the filtered variants and candidate genes—are summarised in Supplementary Data 2–6.

Body Weights at Week 32

The estimated for the three QTLs associated with BW32 on chr 1, 4, and 27 were comparable, at 0.39, 0.43, and 0.38, respectively. Before functional enrichment analysis, the CS for these QTLs contained 3,616 (QTL-chr1), 87 (QTL-chr4), and 47 (QTL-chr27) variants.

For QTL-chr1, two signals were identified: one spanning the entire QTL region and another confined to a narrower region. After applying signal filtration using a p-value threshold of , only the refined signal was retained, consisting of 160 variants. Among these, a variant located at 170,911,032 bp exhibited the highest PPC, increasing from 43% pre-enrichment to 92.3% and 43.9% post-enrichment under the “Consequence” and “Impact” categories, respectively. This variant was annotated as an intergenic variant and modifier, but without “Biotype” information. The signal highlighted three genes — FAM124A, ENSGALG00000052822, and ENSGALG00000053256 — with PPC values ranging from 3.8% to 43% before enrichment. After renormalization by “Impact”, PPC values for these genes remained largely unchanged, however, only ENSGALG00000053256 persisted in the “Consequence” category, where its PPC increased markedly to 92.3%. Conversely, only ENSGALG00000052822 was retained in the “Biotype” category, with a modest PPC increase from 37.2% to 39.1%.

For QTL-chr4, a single signal localized near the NCAPG-LCORL region and its upstream area was identified and retained after filtration, comprising 87 variants. Enrichment estimation was not feasible for the “Consequence” and “Impact” categories due to a lack of annotation variability. When enriched by the “Biotype” category, PPC values for variants annotated as lncRNA decreased from 0.3% to near zero, while those annotated as protein coding remained stable (about 1.3%). At the gene level, NCAPG and LCORL displayed high PPC values of 90.8% and 88.3%, which further increased to 95.1% and 92.5% after enrichment by “Biotype” category.

On QTL-chr27, one refined signal was detected and retained after filtration, comprising 47 variants. Due to uniform “Biotype” annotation (all variants as protein coding), enrichment estimation was not applicable for this category. Variants PPC values ranged from 0.8% to 3.8% before, and nearly zero to 4.3% after enrichment, consistently highlighted IGF2BP1 with the highest PPC values of 61.2% pre-enrichment, 69.4% (“Consequence”) and 64.4% (“Impact”) post-enrichment. Additional genes— GIP, SNF8, UBE2Z, and ENSGALG00000001525 were also identified, with PPC values ranging from 2.2% to 16.3% across categories before and after enrichment.

Egg Weight (EW)

Two QTLs were identified for EW30 on chr 1 and 4, with estimated of 0.15 and 0.73 and CS sizes of 800 and 455 variants, respectively. On QTL-chr1, two initial signals spanning the region highlighted several genes, namely SETDB2, ENSGALG00000050514, FNDC3A, and ENSGALG00000031392, with PPC values ranging from 10.7% to 44.2%. After signal filtration, only one signal retained, comprising 655 variants, with the top PPC increasing from 1.3% to 6.2%. ENSGALG00000053103 became the top gene in the “Consequence” category, with PPC rising from 8.7% to 41.5%. ENSGALG00000050514 was retained and ranked highest across the “Impact” and “Biotype” categories, with PPC increasing from 11.8% to 11.9% and 19.1%, respectively. Notably, ENSGALG00000031392 was pinpointed by both initial signals with different PPC (44.2% and to 3.4%), but showed lower PPC in the retained signal. For QTL-chr4, two signals were detected and retained after filtration: one localized near SLIT2 and LDB2, and another in an intergenic region. The combined signal comprised 455 variants, with the top PPC of 9.9%. Enrichment analysis did not converge for the “Consequence” category and could not be evaluated for “Impact” due to limited annotation variability. The retained signal highlighted LDB2, SLIT2, ENSGALG00000054285, and ENSGALG00000032384, with PPC values unchanged after enrichment in the “Biotype” category (1.3% to 65.4%).

Two QTLs were identified for EW40 on chr 1 and 4, with estimated of 0.06 and 0.20 and CS sizes of 3,355 and 206 variants, respectively. On QTL-chr1, two initial signals were detected, with THSD1 showing the highest PPC (21.2%). After signal filtration, one signal remained (1,245 variants), and the top PPC increased from 1.0% to 7.4%. The gene set also shifted, with ENSGALG00000050514 emerging as the top gene across the “Consequence”, “Impact” and “Biotype” categories (PPC increasing from 7.9% to 26%, 7.9%, and 12.0%, respectively). For QTL-chr4, a single retained signal localized near the NCAPG-LCORL region and its upstream area (206 variants, top PPC of 2.1%). These variants highlighted NCAPG (30%) and LCORL (25.2%) pre-enrichment. Post-enrichment PPC values dropped to 5.7% for NCAPG, and < 1% for LCORL across the “Consequence” and “Biotype” categories. The “Impact” category could not be evaluated due to lack of annotation variability.

For EW50, one QTL was identified on chr 4 (estimated = 0.24, CS = 186 variants). A single signal near the NCAPG-LCORL region and its upstream area, as well as on LDB2, was retained after filtration. Pre-enrichment PPC values highlighted NCAPG (7%), LCORL (2.6%), and LDB2 (11.8%). Post-enrichment, only NCAPG remained, with PPC of 5.2% across the “Consequence” and “Biotype” categories (“Impact” category not evaluable).

For EW70, one QTL was identified on chr 4, overlapping the EW50 region (estimated = 0.21, CS = 1,277 variants). Two signals were initially detected: one overlapping the EW50 locus and another spanning a downstream region. These variants highlighted several genes, notably TAPT1 (35.2%) and PROM1 (29.4%) pre-enrichment. After filtration, only the first signal was retained (176 variants), continuing to highlight NCAPG (3.7%), LCORL (1.7%), and LDB2 (4.9%) pre-enrichment. Similar to EW50, only NCAPG remained after enrichment, with PPC of 2.1% across the “Consequence” and “Biotype” categories (“Impact” category not evaluable).

Overall, three genes located on chr 1 that were previously pinpointed by BW32 were also observed for EW30 and EW40 both pre- and post-enrichment, although their QTL-chr1 regions slightly differed. Among other genes pinpointed on chr 1 for EW30 and EW40, about half overlapped pre-enrichment, and most overlapped post-enrichment. On chr 4, NCAPG and LCORL were consistently identified across BW32 and EW traits (except EW30), both pre- and post-enrichment, while LDB2 was repeatedly highlighted in EW traits (except EW40) before and after enrichment.

Egg Number (EN)

The multi-trait metaGWAS for EN1 to EN13 identified one QTL on chr 2. The estimated for this QTL across different EN measurements was generally low, with a maximum of 0.1, and CS sizes ranging from 31 to 1,807 variants. These variants consistently highlighted ENSGALG00000049955 (PPC = 7 − 27.8% for EN3 – EN13, declining to 1.8% in EN2), and GARS (PPC = 32.5–58.3% for EN3 – EN12, declining to 2% in EN13). Several genes overlapped between the beginning and end of the laying period, including XIRP1, OXSR1, DLEC1, ENSGALG00000005884, and CTDSPL (their PPC = 1.1–18.1%), while ENSGALG00000050406 and ENSGALG00000049302 were mainly detected in the late laying period (PPC = 1.0–3.3%). Additional genes were detected specific to a period, e.g., CRHR2 for EN2, PLCD1 for EN3, and MYD88 for both EN2 and EN3.

After signal filtration, the CS sizes were reduced to 35–38 variants, with no variants retained for EN2 and EN10 – EN13. Following enrichment, GARS and ENSGALG00000049955 continued to be highlighted across EN3 – EN8. GARS showed the highest PPC under the “Impact” category (57.7–96.1%), followed by “Biotype” (1.8–48.7%), but could not be evaluated under “Consequence”. In contrast, ENSGALG00000049955 had the highest PPC under “Biotype” (39.3–96.1%), followed by “Consequence” (13.3–13.6%) and “Impact” (10.5–17.7%). Notably, three genes – ENSGALG00000005934, ENSGALG00000006018 and OXSR1 – were uniquely identified in EN1 after enrichement in the “Impact” category (PPC = 2.5–5.6%).

Enhancing GWAS through Backcrossing Generations and Imputation

A backcrossing scheme involves the introduction of multiple recombination events that can break down long-range LD blocks, which allows for a higher resolution of QTL mapping and refinement of QTLs [23]. In this study, our results showed that GWAS identified more and stronger signals in “IC” than other generations, while GWAS including individuals from all generations thereby resulting to a larger sample size, have identified a greater number of QTLs compared to single-generation analyses. Furthermore, imputation of array-genotyped data significantly increases the density of variants across the genome, which strengthens the statistical power of QTLs that were identified in array GWAS. These findings highlight the benefits of combining breeding strategies with advanced genomics tools to enhance the resolution and detection power of GWAS.

Including correlated traits as covariates in GWAS reveals pleiotropic and trait-specific loci

Covariates play a critical role in GWAS by accounting for confounding factors that might obscure true genetic signals. Specifically, when two traits are correlated, conditioning GWAS on one trait by including it as a covariate can help to disentangle their shared and distinct genetic architectures. In this study, BW32 and EW traits exhibited strong phenotypic and genetic correlations across different time points. To clarify the genetic basis of these traits, we performed GWAS for EW traits while including BW32 as a covariate, and vice versa.

When BW32 was included as a covariate in the GWAS for EW at four time points, the previously identified QTLs on chr 1 disappeared and diminished on chr 4, suggesting that these associations were primarily driven by the effects of BW32 rather than by direct effects on EW traits. In contrast, when EW traits were included as covariates in the GWAS for BW32, the coinciding signals on chr 1 and 4 generally remained significant, although with reduced statistical power (i.e., larger p-values). This attenuation also extended to the signal on chr 27 when EW50 and EW70 were included. Together, these results indicate that the QTLs on chr 1 and 4 are pleiotropic — influencing both BW32 and EW traits — with stronger direct effects on BW32 and secondary, indirect effects on EW traits, prominently on chr 1. Meanwhile, the locus on chr 27 appears to specifically affect BW32, while also showing correlated effects with EW at later developmental stages.

Notably, signals on other chromosomes remained with comparable p-values, suggesting genomic regions that contribute specifically to the trait of interest. In some cases, covariate adjustment even helped clarify spurious associations by reducing or stabilizing p-values, as observed in the GWAS analyses of EW30 and EW70 when BW32 was included as a covariate. Further investigation into the biological mechanisms underlying both pleiotropic and trait-specific QTLs — for example, through functional annotation or gene expression analyses — may provide valuable insights into how these loci influence the traits.

These findings underscore the importance of considering correlated traits as covariate in GWAS to disentangle complex genetic relationships, refine association signals and enhance the accuracy and interpretability of GWAS results.

High variation of heritability estimation

The large variation in estimates observed for both phenotypes and QTL regions in this backcrossing scheme reflects the complex genetic architecture underlying the traits studied, as well as their evaluation in a heterogeneous population. Overall, tended to be higher for BW32, EW at different time points, and EN during the early stage of production (highest for EN1, moderate for EN2). In contrast, EN during the mid-to-late laying stages exhibited substantially lower estimates (low for EN3 — EN13, moderate for EN10 and EN11), reaching zero during the peak laying period (EN4 — EN6). This decline in with advancing age has been previously reported in chicken and is likely due to the increasing influence of environmental factors on EN over time [57, 58]. Our findings partially support this interpretation, as we observed a moderate increase in during EN10 and EN11 compared with the surrounding low periods. Furthermore, the absence of a clear trend of changes of environmental variation in EW with age [59] aligns with the observation of relatively stable estimates for EW across various ages. Interestingly, GWAS results revealed more suggestive or significant variants during periods of moderate-to-low (EN3, EN4, EN7 — EN11), while a highly polygenic architecture was observed during the period with the highest (EN1, = 0.40). For EN, the QTL on chr 2 derived from the multi-trait metaGWAS explained most variance for EN4, EN7, and EN8 ( = 0.03 — 0.09), which were also the stages with the suggestive or significant peak on chr 2 in GWAS, while the genome-wide for these traits were low (0.00 — 0.06). Conversely, EN1 exhibited the highest (0.40) but a smaller contribution from this QTL (0.02). This suggests that the genetic control of EN may shift across the laying cycle, with early production influenced by many small-effect loci, and later stages affected by a few loci with small-effect together with environmental factors.

We occasionally observed regional (QTL-based) estimates for a trait exceeding the genome-wide trait (e.g., for EW30, genome-wide = 0.57, QTL-chr4 = 0.73). This likely reflects local haplotype and linkage disequilibrium in the backcross population, which can lead to inflated variance estimates when the QTL GRM captures not only the genetic effects, but also part of the background covariance due to ancestry. Nevertheless, QTL-based and trait-level were strongly correlated (Pearson’s r = 0.74), indicating that the regional estimates broadly reflect trait genetic architecture. For EW, the QTL on chr 1 — previously reported to disappear after adjusting for BW32 — explained a smaller proportion of variance (EW30 = 0.15; EW40 = 0.06). In contrast, the QTL on chr 4, which showed attenuation in GWAS with BW32 covariate, accounted for substantial variation at the early stage (EW30 = 0.73) and remained moderately influential through later stages (EW40 — EW70 0.2). For BW32, QTLs on chr 1 and 27 explained comparable proportions of variance ( around 0.38), while the QTL on chr 4 explained slightly more ( = 0.43). This pattern further supports the pleiotropic nature of QTLs on chr 1 and 4 for BW and EW traits, with the effect on chr 4 being more prominent.

Fine-mapping pinpoints promising genes for multiple traits

In this study, fine-mapping analyses using a Bayesian approach identified several promising candidate genes with pleiotropic effects across multiple traits. Notably, NCAPG and LCORL were highlighted as influencing BW32 and all EW measurements except EW30, supported by comparatively high PPC values. NCAPG encodes a subunit of the condensin complex essential for chromosomal regulation during mitosis and meiosis, while LCORL encodes a ligand-dependent nuclear receptor corepressor-like protein involved in transcriptional regulation via RNA polymerase II. Both genes have been extensively linked to growth traits across species. In cattle, the NCAPG–LCORL locus is robustly associated with growth, feed intake, and carcass traits, as supported by genomic association studies, gene expression analyses, and haplotype mapping [60–62]. In horses, this locus influences body size and conformation traits with strong GWAS signal [63, 64]. Selective sweep analysis has shown strong selection signature of LCORL associated with stature in pigs [65]. Recently, Bai et al. confirmed causal mutations in the NCAPG–LCORL region affecting body weights, through combined analyses of selective sweep across domesticated species and gene-editing studies in mice [66].

Despite strong evidence for direct effects on body weight, direct biological links between NCAPG–LCORL and egg weight remain limited. In our study, both pre- and post-enrichment PPC values for NCAPG–LCORL were consistently lower for EW than for BW, showing a sharp decline with increasing age. For example, while PPC values for EW40 reached approximately 30%, they dropped to below 7% in EW50 and EW70 before enrichment, and further declined to below 6% after signal filtration and enrichment from EW40 to EW70, highlighting a substantial reduction in their contribution. Combined with the strong phenotypic and genetic correlations observed between BW and EW, and the loss of these significant signals in GWAS when BW32 was included as a covariate, these findings suggests that NCAPG–LCORL may influence EW primarily through genetic effects mediating body weight, rather than a direct impact on egg weight itself.

Similarly, Fine-mapping of BW32 highlighted IGF2BP1 (insulin-like growth factor 2 mRNA-binding protein 1). IGF2BP1 is a well-known candidate gene associated with growth and body size across various livestock species [67–71]. In chickens, causal promoter variants in IGF2BP1 regulating body size were identified via pan-genome analysis, particularly between 3 to 12 weeks of age [72]. Additionally, IGF2BP1 has also been shown to promote fat and muscle development in chickens, especially at 8 weeks and 120 days, respectively (73, 74). It’s expression level has been reported to decline with age in poultry [69, 74, 75], which may partly explain its reduced effect on body weight at later developmental stages. Interestingly, IGF2BP1 alone has been reported to inhibit cell proliferation unless coupled with post-translational modification, suggesting that its biological activity may be context-dependent [75]. In our study, BW was measured at 32 weeks of age, which is later than in previous studies, yet IGF2BP1 was still consistently pinpointed with relatively high PPC values that even increased slightly after functional enrichment. This finding underscores the continued importance of IGF2BP1 and its potential long-term impact in influencing body weight.

Apart from IGF2BP1, Fine-mapping for BW32 also identified several nearby protein-coding genes, including GIP (involved in glucose response), SNF8 (protein transport), UBE2Z (apoptosis), and ENSGALG00000001525 (autophagosome regulation; NCBI gene name for CALCOCO2). These genes have been previously reported to be associated with breast muscle weights in chicken [76] and abdominal fat percentage [77], suggesting they may contribute to body weight variation through their roles in energy metabolism.

Fine-mapping identified several lncRNA genes on chr 1 associated with BW32, EW30 and EW40, including ENSGALG00000053256, ENSGALG00000052822, and FAM124A. These genes have been reported to be associated with body weight in chickens [78, 79], as well as with stature [80] and skeletal muscle production [81]. This suggests that their influence is primarily on body weight, which aligns with the substantially higher PPC values computed for BW32 compared to EW40 in our fine-mapping results.

On chr 4, fine-mapping pinpointed several genes associated with EW across different periods, including LDB2 and ENSGALG00000032384, as well as genes specific to EW30 (SLIT2, ENSGALG00000054285) and EW70 (TAPT1, PROM1, FBXL5, CD38). Among these, LDB2, SLIT2 and TAPT1 have been widely reported in the literature [3, 12, 24, 71, 82–88], while others such as PROM1, CD38 and ENSGALG00000054285 have been less frequently reported in chickens [82, 83]. Many of these genes are primarily linked to growth traits or body weight. Given their moderate to low PPC values, the loss of corresponding GWAS signals for EW after accounting for BW, and their biological functions — including transcription regulation (LDB2), regulation of cellular protein (FBXL5), cell signalling and regulation (SLIT2, CD38, PROM1, TAPT1) — these loci likely contribute to egg weight indirectly through body growth pathways.

For EN, fine-mapping identified several genes associated across and within specific laying periods. After signal filtration, two genes were consistently linked from EN3 to EN8, while others were unique to EN1. However, these genes generally did not overlap with findings reported in prior studies, in line with the observation that many QTLs identified for egg production in chickens re rarely replicated [1, 2, 57]. Despite this, their molecular and biological functions provide insights into potential mechanisms influencing reproductive traits.

At the early stages, functional annotations indicated involvement in oxidative stress response (OXSR1), transmembrane transport (ENSGALG00000005934, ENSGALG00000006018), and actin cytoskeleton organization (XIRP1) for EN1; cell surface receptor signalling (CRHR2) for EN2; lipid metabolism (PLCD1) for EN3; and Immune signalling (MYD88) for both EN2 and EN3.

At later stages, apart from OXSR1 and XIRP1, additional processes include defence response to tumour cells (DLEC1), protein phosphorylation (ENSGALG00000005884), and cell cycle regulation and differentiation via tumour suppressor activity (CTDSPL). Genes consistently identified across multiple periods, both before and after enrichment, were involved in protein biosynthesis (GARS) and included non-coding elements (ENSGALG00000049955).

Together, the biological processes associated with these genes suggest that reproductive performance may be influenced through diverse pathways, particularly those related to cell development, metabolism and immune regulation.

Overall, fine-mapping of GWAS-identified QTLs has successfully uncovered promising genes for multiple traits. However, it is important to keep in mind that the coverage and quality of imputed genotypes can limit the resolution of fine-mapping and downstream analyses such as signal filtration by a p-value threshold, and functional enrichment. In this study, imputation substantially increase the number of variants from ~ 46.6k to ~ 4.9M, yet the resulting dataset is still incomparable to the completeness of WGS data. Furthermore, the concordance between imputed and array genotypes did not exceed 96%, which introduces uncertainty. This limitation may lead to missing causal variants, particularly structural or rare variants, despite the efforts to extend QTL regions by 2Mb upstream and downstream to capture as many informative markers as possible. These factors should be carefully considered when interpreting fine-mapping results, especially in the context of imputed data.

In conclusion, our study demonstrated that GWAS achieves greater detection power and statistical strength when all individuals from a backcrossing scheme are combined and genotype imputation is applied to increase marker density. Moreover, incorporating correlated traits as covariates in GWAS offers valuable insights, helping to distinguish pleiotropic loci from trait-specific signals.

This study is the first to apply Bayesian fine-mapping to GWAS-identified QTLs in chickens, leveraging a backcrossing scheme to successfully pinpoint key genes underlying economically important traits. Our findings not only confirmed previously reported loci but also provided new insights into the complexity of trait inheritance. Nevertheless, the coverage and accuracy of genotype imputation may limit the fine-mapping resolution. Overall, this study highlights the power of Bayesian fine-mapping approaches to refine GWAS signals and deepen our understanding of the genetic architecture of complex traits, offering promising targets for future validation and breeding applications.