Genetic patterns and genome-wide association analysis of eggshell quality traits of egg-type chicken across an extended laying period

INTRODUCTION

Eggs are rich in vitamins, minerals, protein, and lipids, being an affordable and easily digestible food for humans. The eggshell can protect the internal contents during the transportation and storage of eggs and can prevent bacterial invasion. Hence, enhancing eggshell quality has always been a key focus in poultry breeding.

Due to the improvements of the production performance of egg-type chickens, the laying period has been extended from the initial 72 wk to 80 wk of age (Liu et al., 2018). Besides, the breeding goal of achieving a production of 500 qualified eggs by 100 wk of age represents a novel trend in poultry breeding, which can reduce the necessity for parental stocks, thus curbing the number of reared pullets in commercial poultry farms (Bain et al., 2016). However, a key challenge is the decline in eggshell quality during the late laying period.

Eggshell quality is a comprehensive concept, relating to the traits being mostly quantitative, such as eggshell strength, eggshell thickness, and eggshell color. Conventional breeding approaches have become inadequate to meet the breeding demands. Advances in single nucleotide polymorphism (SNP) arrays and sequencing technologies have rendered Genome-Wide Association Analysis (GWAS) a highly effective approach for detecting the effective genetic variants in farm animals. Liu (2011) conducted the first GWAS with the Illumina 60K SNP array, revealing the genetic associations with egg quality traits in chickens. Sun et al.(2015) discovered numerous variants and genes that associated with egg yolk weight and ovary weight in chickens. However, most research on extending the laying period has only focused on the nutritional needs of hens during the late laying period (Lv et al., 2019). Hence, there is a lack of genetic research on eggshell quality traits throughout the entire laying period, especially in the late laying period. Furthermore, previous research on eggshell quality generally focuses on the peak laying period, indicating a research gap in the study of age-dependent eggshell quality.

Here, we employed a total of 3,261 Rhode Island Red chickens to investigate the genetic basis underlying the age-dependent eggshell quality by GWAS, SNP enrichment analysis and functional annotation of cross-tissue regulatory elements, laying a valued scientific foundation for chicken molecular breeding.

MATERIALS AND METHODS

).

Estimation of Genetic Parameters and Contribution to Phenotypic Variance (CPV)

In this study, SNP-based genetic parameters were estimated by the software of GCTA v1.93.2. Restricted maximum likelihood (REML) was employed for variance component estimation. The variance components model for explaining the genetic variation using whole-genome SNP markers was defined as follows:(1)

where y represents the phenotype values, β denotes fixed effects, X is the correlation matrix associated with β, g is individual genetic effects, and ε is random error. V represents the covariance matrix, A is the genomic relationship matrix (GRM) computed by GCTA to represent the genetic relatedness between individuals, σ2 represents the additive genetic variance of the SNPs, I is the identity matrix, and σ2 ε is the residual variance. Additionally, we calculated the phenotypic variance contribution of SNPs that reached the genome-wide significance in GWAS analysis based on the genetic matrix.

Genotyping and Quality Control

Genomic DNA was extracted from whole blood samples using standard phenol/chloroform methods and the 3262 qualified hens were genotyped with a custom Illumina 50K BeadChip, named PhenoixChip-I (Liu et al., 2021). Out of the detected 44,387 SNPs, only autosomal SNPs were retained for further analysis. Quality control was performed using PLINK v1.90 software, with the criteria of individual missing rate (mind) < 0.05, locus missing rate (geno) < 0.05, minor allele frequency (MAF) > 0.01, and Hardy Weinberg equilibrium (HWE) < 1e–6. The missing genotypes were imputed using Beagle v4.0 software. Additionally, LD pruning was conducted on the entire dataset using PLINK v1.90 software with the command “-indep-pairwise 25 5 0.2″. Finally, a total of 3,259 individuals and 41,838 SNPs distributed across 33 autosomes were retained, meeting the requirements for subsequent genome-wide analysis. SNP annotation was performed with the Gallus-gallus 6.0 reference genome using the SnpEff software.

Genome-Wide Association Analysis

The kinship matrix was built through the independent SNP markers. The principal components were calculated from the linear combination of markers by the eigenvectors of the kinship matrix, which were treated as covariates and included in the model of GWAS as fixed effects to control the population structure effects.

The univariate GWAS was first implemented with a linear mixed model to account for the associations between each trait and the effective SNPs, which was supplied by GEMMA software. The statistical model was as follows:(2)

where y is the phenotypic values of n individuals in the population of interest; X is a matrix of covariates (fixed effects: top 5 principal components and a column of 1 s) controlling for population structure; α is a vector of corresponding effects that compose the intercept; Z is the marker genotypes; β is the corresponding marker’s effect; Wμ is the Random Effect; and e is the vector of random residuals.

Additionally, the P value from the likelihood ratio test was chosen as a benchmark to evaluate the significance of the association between SNPs and traits. The threshold for genome-wide significance was established using a modified Bonferroni correction implemented through the R package simpleM. A total of 11,154 valid SNPs were carried out, and the threshold was defined as 4.48E–06 (0.05/11,154) and 8.97E–05 (1.00/11,154), respectively. Finally, the results were visualized through the R package “qqman”.

QTL Region Definition and Enrichment Analysis

We used GWAS results to define quantitative trait loci (QTL) as a chromosomal region where the distance between adjacent pairs of significant variants was less than 1 Mb (175.55 ∼ 176.88 Mb; 22.76 ∼ 24.97 Mb). Within each locus, we identified the most significant variant as the lead variant. We only considered a maximum distance of 0.5 Mb on either side of the lead variant Enrichment analysis and it was conducted by an R package, LOLA. Located chromosomal regions were searched directly on a website created by Pan et al. (2023) and functionally annotated for cross-tissue regulatory elements. The website (http://genome.ucsc.edu/s/zhypan/galGal6_FAANG_V1

RESULTS

Phenotype Statistics and Estimation of Genetic Parameters

Phenotype data regarding 4 eggshell quality traits, namely EST, ESS, ESC, and ESG, were collected at 7 age points spanning from 36 to 90 wk of age. The descriptive statistics suggested the decreased phenotype values and increased phenotypic variations of all traits along with the aging process, especially after 55 wk of age (Table 1). Notably, ESG exhibits the highest coefficient of variation, ranging from 19.59 to 25.89%, in comparison to the other traits which range from 4.84 to 23.63%.

Table 1. Descriptive statistics of eggshell quality traits along with aging process

Traits	N	Mean	SD	CV (%)	Min	Max	h2(SE)
EST36	2,919	0.2834	0.02	7.04%	0.2130	0.3800	0.26 (0.03)
EST46	843	0.2964	0.02	8.34%	0.2130	0.3870	0.26 (0.06)
EST55	2,800	0.2821	0.02	7.16%	0.2300	0.3850	0.23 (0.03)
EST66	2,667	0.2837	0.02	8.03%	0.2127	0.3970	0.21 (0.03)
EST72	2,597	0.2832	0.02	8.55%	0.2270	0.4030	0.17 (0.02)
EST80	2,478	0.2838	0.02	8.13%	0.2130	0.3830	0.08 (0.02)
EST90	2,287	0.2873	0.02	7.47%	0.2370	0.4000	0.11 (0.02)
ESS36	2,907	3.530	0.66	18.73%	1.078	5.616	0.32 (0.03)
ESS46	843	3.455	0.46	13.20%	1.516	5.033	0.18 (0.05)
ESS55	2,786	3.413	0.70	20.43%	1.013	6.826	0.22 (0.03)
ESS72	2,574	3.269	0.68	20.66%	1.009	5.593	0.22 (0.03)
ESS80	2,432	2.982	0.65	21.94%	1.013	5.276	0.22 (0.03)
ESS90	2,243	2.861	0.69	23.63%	1.004	5.779	0.15 (0.03)
ESC36	2,919	58.81	3.39	5.77%	45.5	78.21	0.37 (0.03)
ESC46	843	57.58	2.79	4.84%	45.34	74.81	0.30 (0.06)
ESC55	2,802	60.13	3.56	5.93%	49.36	80.28	0.44 (0.03)
ESC66	2,674	59.64	3.67	6.16%	47.51	81.24	0.38 (0.03)
ESC72	2,611	60.81	3.80	6.25%	46.83	77.40	0.44 (0.03)
ESC80	2,482	60.24	3.90	6.47%	48.79	81.67	0.40 (0.03)
ESC90	2,295	58.45	4.22	7.22%	45.65	80.47	0.34 (0.03)
ESG36	2,918	3.224	0.65	20.26%	1.230	5.820	0.30 (0.03)
ESG46	843	3.115	0.61	19.59%	1.570	5.380	0.33 (0.06)
ESG55	2,801	2.392	0.55	22.82%	1.260	5.380	0.43 (0.03)
ESG66	2,671	2.494	0.66	25.89%	1.260	5.650	0.31 (0.03)
ESG72	2,611	2.606	0.59	22.76%	1.360	5.580	0.39 (0.03)
ESG80	2,479	2.484	0.59	23.93%	1.310	7.010	0.30 (0.03)
ESG90	2,294	2.550	0.53	20.92%	1.320	4.710	0.08 (0.02)

Abbreviations: EST, eggshell thickness; ESS, eggshell strength; ESC, eggshell color; ESG, eggshell gloss. Numbers represent the age in wk. N, number of samples. Mean, average value. SD, standard deviation. CV (%), coefficient of variation. Min, minimum value. Max, maximum value. h²(SE), SNP-based heritability (standard error).

The SNP-based estimates of heritability for each trait, as well as the genetic and phenotypic correlations, were calculated (Additional file 1: Table S1). The analysis of heritability reveals that ESC exhibits relatively high heritability (0.30-0.4), while ESS, EST, and ESG display moderate levels of heritability (0.15-0.32; 0.08-0.26; 0.08-0.43) (Table 1). An overview of these heritability measures suggests that the traits at the late laying period exhibit lower heritability compared to the early and peak laying periods, with ESG notably declining from 0.30 at 36 wk of age to 0.08 at 90 wk of age. Correlation analysis reveals a significant positive correlation between EST and ESS, both genetically (0.31-0.41) and phenotypically (0.72-0.75). Additionally, EST and ESG display a notably high negative genetic correlation ranging from -0.26 to -0.42. Conversely, ESC and ESS exhibit a moderately high positive genetic correlation (0.17-0.23).

In traditional breeding schedule, the phenotype value of eggshell quality traits at 36 wk of age is used to select the top breeding chickens, while the genetic correlation between 36 and 90 wk of age varied greatly among different traits (Figure 1). The genetic correlation of ESC between 36 and 90 wk of age reached to 0.97, indicating that early selection can still have an excellent genetic improvement on ESC at 90 wk of age. In contrast, the genetic correlations for the other traits between 36 to 90 wk of age declined to 0.7-0.8 (ESS: 0.77; ESG: 0.86; EST: 0.74), suggesting that the selection time for these traits should be postponed to attain the same genetic effect.

1. Download : Download high-res image (1MB)
2. Download : Download full-size image

Figure 1. Estimates of SNP-based heritability (on the diagonal) and genetic (above the diagonal) and phenotypic correlations between eggshell quality traits. Diagonal, SNP-based heritability (standard error); Upper triangle, SNP-based genetic correlation; Lower triangle, phenotypic correlation; Top left corner, age.

Genome-Wide Association Study

The GWA tests were conducted for all eggshell quality traits using a univariate linear mixed model and identified a total of 200 unique candidate SNPs that were genome-wide significantly or suggestively associated with these 4 eggshell quality traits (Additional file 2: Table S2). Among them, 4 SNPs located on GGA1 spanning from 175.95-176.05 Mb were significantly associated with ESS36 (P value < 4.48E-06). In addition, for ESS46, one SNP on GGA1 (1_23,249,882_G_A) also reached the genome-wide significant level (P value < 4.48E-06). The manhattan and QQ plots for ESS36 and ESS46 are shown in Figure 2. For other traits, the genome-wide suggestive variants were identified (P < 8.97E-05). Specifically, a total of 35 SNPs were screened for EST, 92 SNPs for ESS, 42 SNPs for ESC, and 30 SNPs for ESG (Table S2 and Figure S1). With the aging process, there is an alteration in the positions of significant loci (Table 2). For example, considering ESS, potentially significant SNPs were identified on GGA1 in the peak laying period (36 and 46 wk of age), GGA3 and GGA8 in the mid-laying period (55 and 66 wk of age), and located on GGA10 during the late-laying period (72, 80, and 90 wk of age).

1. Download : Download high-res image (845KB)
2. Download : Download full-size image

Figure 2. Manhattan and Q-Q plots for GWASs of ESS36 and ESS46. Each dot on this figure corresponds to a SNP within the dataset, while the horizontal red and blue lines denote the genome-wide significant (4.48E-06) and suggestive significant thresholds (8.97E-05), respectively. The Manhattan plot contains -log10 observed P-values for genome-wide SNPs (y-axis) plotted against their corresponding position on each chromosome (x-axis), while the Q-Q plot contains expected -log10-transformed P-values plotted against observed -log10-transformed P-values. GIF denotes the genomic inflation factor indicating the degree of population stratification.

Table 2. The distribution of significant loci of eggshell quality traits with aging process

Trait	Laying period	Chromosome	No. of SNP
ESS	Pre	1	46
ESS	Mid	3;8	20
ESS	Late	10	26
EST	Pre	1;5	3
EST	Mid	3;6;13;17	18
EST	Late	3;5;7;8	14
ESC	Pre	2;5;13;16	8
ESC	Mid	12	12
ESC	Late	1;4	22
ESG	Pre	2;3;4	20
ESG	Mid	1; 9	4
ESG	Late	2;13;18	6

Abbreviations: Pre, pre-laying period, 36 and 46 wk of age; Mid, mid-laying period, 55 and 66 wk of age; Late, late-laying period,72, 80 and 90 wk of age. Chromosome, the chromosome where the main significant loci located. No. of SNP, the number of the SNPs that reached the genome-wide suggestive threshold (P < 8.97E-05).

SNP Enrichment Analysis

By combining all 5 epigenetic marks across 23 tissues via ChromHMM, Pan et al. predicted 15 distinct chromatin states. These chromatin states mainly represented promoters (TssA, TssAHet, and TssBiv), TSS-proximal transcribed regions (TxFlnk, TxFlnkWk, and TxFlnkHet), enhancers (EnhA, EnhAMe, EnhAWk, EnhAHet, and EnhPoi), accessible islands (ATAC islan), repressed regions (Repr and ReprWk), and quiescent regions (Qui). Of all the identified chromatin states, enhancer activities were the most dynamic, whereas promoter activities remained relatively consistent among tissues(Pan et al., 2023).

To determine whether the genome-wide significant SNPs exhibit enrichment in specific chromatin states, an analysis of SNP enrichment was performed. The analysis revealed a higher enrichment of active regulatory elements in comparison to the 15Qui region. Specifically, the regions E1(9.3%; 9.5%), E2(7.8%; 7.2%), E3(8.2%; 8.2%), E4(10.3%; 9.7%), and E5(8.9%; 8.7%) exhibited notable enrichment (Figure 3). Notably, E1 and E2 correspond to promoters, while E3, E4, and E5 represent TSS-proximal transcribed regions. Previous studies have consistently observed the conservation of promoters in the chromatin state, with E1(TssA) displaying the highest level of enrichment in DNA sequence level among different tissues, which aligns with the findings of this study. Besides, the SNP enrichment analysis indicated that E1(TssA), as the most prominent region, primarily contains genes involved in a variety of disorders related to brain development, immune response, and intestinal functions. Additionally, this study provided evidence of its association with ESS.

1. Download : Download high-res image (356KB)
2. Download : Download full-size image

Figure 3. Single nucleotide polymorphism enrichment results for ESS36 and ESS46. The violin diagram, horizontal axis, E represents the classification of regulatory element functional levels, where numbers 1 to 15 indicate the strength of regulatory ability; the vertical axis is the P-value.

Functional annotation of regulatory elements cross-tissues can provide valuable insight into the multifaceted functional information to elucidate the results of GWAS and facilitate the discoverer of genetic molecular mechanisms underlying economic traits. To explore whether tissue-specific regulatory elements exerted greater regulatory effects on GWAS SNPs compared to other tissues, an SNP enrichment analysis was conducted using tissue-specific EnhAs for traits ESS 36 and ESS46 (Figure 4).

1. Download : Download high-res image (350KB)
2. Download : Download full-size image

Figure 4. The regulatory elements involved in ESS36. The vertical axis represents different tissues, while the horizontal axis represents the SNP region that reaches the GWAS significant level.

The results demonstrated a significant enrichment of enhancer elements specific in shell gland and kidney tissues for ESS36. These findings emphasize the critical involvement of genes undergoing active regulation in these tissues to influence eggshell formation, and ultimately its strength. This observation aligns with prior research highlighting the significant role of the hypothalamic-pituitary-gonadal axis in egg production among laying hens(Tan et al., 2021).

DISCUSSION

The eggshell plays a crucial role in the storage and incubation of eggs, with its synthesis process affected by various factors, including stress, age, medications, diseases, feed, and environmental conditions (Solomon, 2010). In recent years, the focus and challenge have revolved around maintaining eggshell quality during the late stages of egg laying. During this phase, the eggshell quality notably declines compared to the peak laying period, resulting in substantial economic losses for the poultry industry. Consequently, the extension of the laying period heavily depends on breeding strategies that prioritize enhancing eggshell quality.

The results of descriptive statistics (Table 1) indicate a consistent decline in various eggshell quality traits as the aging process of chicken. This downward trend can be attributed to the decline in physiological functions and reproductive capabilities of the laying hens (Rattanawut et al., 2018; Gan et al., 2020). Previous studies have shown an increase in egg weight and volume as laying hens aging, but without significant changes in eggshell weight (Molnar et al., 2016). However, the decline in eggshell quality, specifically EST, and ESS, leads to an increase in egg damage, which negatively impacts egg storage, transportation, and incubation processes (Tumova and Gous, 2012).

Eggshell color (ESC), although not directly impacting the internal quality or nutritional value of eggs, remains a significant factor considered by consumers when evaluating egg quality. The brown coloration of chicken eggshells, influenced by the protoporphyrin-IX pigment, is regulated by multiple genetic factors. Hence, there can be considerable variations in eggshell color even within the same breed (Zheng et al., 2014). As the aging of hens, their ability to secrete protoporphyrin-IX pigments decreases to varying degrees, resulting in changes in laying performance and eggshell color, leading to an increase in the coefficient of variation of ESC.

The ESG is closely associated with embryo development, as a smooth eggshell surface plays a crucial role in preventing moisture from blocking eggshell pores (Li et al., 2020). This water resistance is vital in reducing surface tension and facilitating gas exchange for the developing embryo (Igic et al., 2015). Furthermore, a high reflectance on the eggshell surface serves as a protective measure against potential solar radiation damage posed to the embryo (Samiullah et al., 2016). The present study discovered that ESG maintains comparatively high values from 72 to 90 wk of age. This is believed to be associated with the decline in egg production rate and the lengthening of egg-laying intervals during the later laying stage. This suggests that more substances are progressively deposited on the outermost layer of the eggshell within the hen’s oviduct over time, contributing to increased brightness.

Estimated genetic parameters suggest that ESC (Zhang et al., 2005) and ESG are highly heritable traits, whereas ESS and EST (Liu et al., 2018) exhibit moderate heritability, which aligns with previous research (Liu et al., 2018). Compared to the early laying period, there is a significant decrease in the heritability of traits during the late laying phase. This decrease could be attributed to a wide range of organ aging in hens during this particular stage, resulting in an increase in phenotypic variations of eggshell-related traits (Liu et al., 2019). Meanwhile, this indicates that with the aging process, genetic factors have a relatively weaker impact on eggshell quality, so we can improve eggshell quality by controlling external environmental conditions (Cheng and Ning, 2023). Previous studies had found that dietary mulberry-leaf flavonoids (MF) could ameliorate the eggshell quality of aged hens by improving antioxidative capability and Ca deposition in the shell gland of uterus (Huang et al., 2022). In addition, eggshell quality can be improved by organic trace elements that also improve eggshell ultrastructure, mineral deposition in the eggshell, antioxidant capacity, and immune function in the late laying period of laying hens(Qiu et al., 2020; Zhang et al., 2021).

To date, many studies have focused on unraveling the genetic determinants of eggshell quality traits, with the majority of investigations centered on F2 hybrid populations or specific egg laying stages. This study represents the pioneering use of a univariate mixed linear model for GWAS analysis throughout the entire laying period (up to 90 wk of age) in the Rhode Island Red chicken population. GWAS and SNP annotations indicate FRY as the key candidate gene influencing ESS36. Previous studies have found that FRY encodes a microtubule-binding protein that plays a role in diverse biological processes such as cell skeleton assembly and cellular movement (Dadousis et al., 2021). Besides, as a key component of chromosome separation in metaphase, FRY inhibits the growth of cancer cells in a Hippo pathway-dependent manner in some malignancies (Liu et al., 2019; Irie et al., 2020; Mai et al., 2022). Dadousis et al., (2021) reported the involvement of FRY in the regulation of body weight in broiler chickens. In this study, FRY was found to be related to eggshell strength. Our finding suggests that FRY functions as a pleiotropic gene, influencing not just body weight accumulation, but also the eggshell strength during certain developmental stages in chickens. There is currently a lack of in-depth functional investigations concerning these genes in chickens. Based on respective studies in humans, we postulate that these genes may modulate eggshell strength through the regulation of cellular skeleton development. Future research should consider investigating the functional mechanisms of the FRY gene during the eggshell formation process in chickens, as well as its specific correlation with eggshell strength.

The SNP enrichment analysis illustrates a higher enrichment level of regulatory elements in active regions compared to inactive regions. This is particularly evident in regions E1 to E5 (promoters and proximal transcribed regions to Transcription Start Sites [TSS]). Among these, E1(TssA) as the promoter, exhibits the most significant enrichment and shows a correlation with ESS. The functional annotation of cross-tissue regulatory elements reveals enrichment of enhancer elements for ESS36, specifically in shell gland and kidney tissues. Previous studies have shown that they are associated with calcium transport. The kidneys and shell gland are rich in calcium-binding protein (CaBP), and studies have shown that when calcium is deficient, eggshell quality decreases and CaBP increases in the kidneys and uterine glands (Bar and Hurwitz, 1984). The observed enrichment in the shell gland could be largely attributed to the process of eggshell formation, which occurs in the shell gland over a span of approximately 18 to 20 h (Wolfenson et al., 1982; Poyatos et al., 2020). Furthermore, enrichment in the kidney is mainly attributed to its role as a major organ for calcium transport, which indirectly has a major impact on eggshell formation (Xin et al., 2021). Overall, the thorough annotation of regulatory elements in the chicken genome will facilitate the analysis of the molecular mechanisms underlying complex traits and provide a valuable resource for improving economically important traits in poultry.

CONCLUSIONS

This study found that eggshell quality traits have lower heritability during the late laying period compared to the early and peak laying periods. Through GWAS analysis, FRY and PCNX2 have been identified as significant candidate genes for eggshell strength. SNPs reached the genome-wide suggestive threshold (P < 8.97E-05) explaining 9.59% and 0.48% of the phenotypic variations to ESS46 and ESS36, respectively. SNP enrichment analysis and cross-tissue regulatory elements functional annotation implied an enrichment of enhancer elements specific to shell gland and kidney tissues. These findings provide valuable insights for the development of breeding strategies aimed at extending the laying period.

Source: Science Direct