Genetic variations for the eggshell crystal structure revealed by genome-wide association study in chickens

Background Eggshell is a bio-ceramic material comprising columnar calcite (CaCO3) crystals and organic proteinaceous matrix. The size, shape and orientation of the CaCO3 crystals influence the microstructural properties of chicken eggshells. However, the genetic architecture underlying eggshell crystal polymorphism remains to be elucidated. Results The integral intensity of the nine major diffraction peaks, total integral intensity and degree of orientation of the crystals were measured followed by a genome-wide association study in 839 F2 hens. The results showed that the total integral intensity was positively correlated with the eggshell strength, eggshell thickness, eggshell weight, mammillary layer thickness and effective layer thickness. The SNP-based heritabilities of total integral intensity and degree of orientation were 0.23 and 0.06, respectively. The 621 SNPs located in the range from 55.6 to 69.1 Mb in GGA1 were significantly associated with TA. PLCZ1, ABCC9, ITPR2, KCNJ8, CACNA1C and IAPP, which are involved in the biological process of regulating cytosolic calcium ion concentration, can be suggested as key genes regulating the total integral intensity. Conclusions The findings greatly advance the understanding of the genetic basis underlying the crystal ultrastructure of eggshell quality and thus will have practical significance in breeding programs for improving eggshell quality. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-08103-1.


Background
In the global egg industry, the integrity of the eggshell is critical for the economic viability of poultry production. The shells must be strong enough to prevent failure during packing and transportation [1]. The eggshell also forms an embryonic chamber for the developing chick by providing mechanical protection, allowing gas exchange and preventing microbial contamination [2,3]. The process of eggshell formation is also a classic model for studying biomineralization [4]. The avian eggshell is composed of calcium carbonate, organic matrix and other minerals [5][6][7]. Several factors have been shown to affect eggshell structure [5,8,9] with most studies still based on optical observation [10]. Because of the difficulty of measurement, few studies have focused on the crystal ultrastructure of eggshells. X-ray crystallography has shown that calcium carbonate in avian eggshell is in the form of calcite crystals [11,12]. The microstructural characteristics influenced by the intricate interlacing of the crystal units have been reported to vary significantly from one eggshell to another, even within the same species [13][14][15]. The eggshell consists of ceramic material whose mechanical properties are determined by its microstructure [16]. Chicken eggshells consisting of highly-oriented crystals of abnormal sizes have been reported to be generally weaker than those consisting of smaller and less well-oriented crystals [14,17]. A comparison of eggshells at different ages showed that aged hens produced a greater variation in grain morphology and crystallographic texture than young hens [17]. The crystallographic texture of eggshells is also correlated with the thickness and strength of the shell. The crystal size and orientation has been found to be strongly positively correlated (r = 0.65) with the thickness of the eggshell [18]. The degree of orientation of crystals has been found to explain approximately 40 % of the variance in eggshell strength [17]. All these studies have indicated that the ultrastructure of eggshell was the key factor affecting eggshell quality.
Eggshell texture and ultrastructure are regulated by organic matrix proteins, which control the mineralization process and influence the biomechanical properties of the eggshell [19]. Association studies between polymorphisms of genes encoding shell proteins and shell characteristics have revealed that certain alleles were correlated with the mechanical properties of the eggshell [20,21]. Ovalbumin and ovo-transferrin (LTF) markers were found to be associated with crystal size, while ovocleidin-116 and ovocalyxin-32 (RARRES1) markers were associated with crystal orientation [18]. Because the crystal is the basic component of construction of the eggshell, the alleles of matrix proteins and genes known to be involved in eggshell formation have also been found to be associated with crystal size and orientation [7,18,21]. Therefore, genetic selection of hens with suitable eggshell ultrastructure properties is important for improving eggshell quality. However, the genetic mechanism of the eggshell crystal ultrastructure is still unclear.
In the present study, the crystal structure of chicken eggshells was determined using X-ray crystallography. To study the genetic basis of biomineralization, a genome-wide association study (GWAS) was performed in the designated F2 population by using a high-density 600k SNP chip for chickens. We aimed to explore the genetic architecture of eggshell crystal structure to identify candidate genes, which may be valuable for the genetic improvement of eggshell quality traits.

Phenotypic description and genetic parameters
The descriptive statistics for nine reflection peaks (A1-A9), total integral intensity (TA) and degree of orientation (OD) were shown in Table 1. The estimates of SNP-based heritability, the genetic and phenotypic correlations between eggshell crystal structure traits were shown in Table 2. The individual integral intensity exhibited higher phenotypic variation (11.71-21.23 %) than that of TA (9.40 %). Thus, the TA could provide a good estimate of the size of the eggshell crystals. The SNP-based heritability (h 2 snp ) of TA was 0.23. The phenotypic correlations between every individual integral intensity (A1-A9) were low (r=-0.15-0.35). The absolute values of the genetic correlations were all higher than those of the phenotypic correlations. The variation in OD was higher (23.72 %) than that of the other traits and the value of h 2 snp was low (0.06). In a previous study, we also measured the ultrastructure of these eggshells [22]. The phenotypic correlations between the crystal structure and eggshell ultrastructure are shown in Fig. 1. TA was positively correlated with eggshell strength (ESS), eggshell thickness (EST), eggshell weight (ESW), mammillary layer thickness (MT) and effective layer thickness (ET) (r = 0.14-0.50), but OD was not significantly correlated with these eggshell quality traits (p < 0.05). There was a high genetic correlation between the total peak area intensity and eggshell ultrastructure (r = 0.52), and a medium genetic correlation between the OD and the eggshell ultrastructure (r = 0.31) (Table 3). However, the genetic correlation of crystal structure with the eggshell ultrastructure all exhibited a large standard error (Table 3).

Genome-wide association study (GWAS)
The association analyses revealed a total of 621 genomewide significant SNPs associated with TA (Additional file 1: Table S2). No genome-wide significant locus was associated with OD. In addition, 1096 and 1 SNPs exhibited associations with TA and OD under suggestive levels of significance, respectively. The Manhattan and quantile-quantile (QQ) plots in Fig. 2 showed a global view of the P-values for all SNPs affecting the ultrastructure of the eggshell crystals. All significant SNPs associated with TA were in a 13.5-Mb region between 55.6 and 69.1 Mb on chromosome 1 (GGA1).  To identify independent SNPs, a further series stepwise conditional analyses were performed and the locus rs315989578, that was found to be significantly associated with TA was fitted into the model to examine these associations. The levels of significance for all loci around the most significant SNP decreased substantially (Fig. 3A) during the stepwise conditional Genome-wide Association (GWA) analysis, with no other independent loci being detected. The Fig. 3B illustrated the linkage disequilibrium between the regions associated with TA   Table S1). The result showed that no significant SNPs were associated with the A1, A2, A3, A4, A5, A7 and A8. 163 and 330 genome-wide significant SNPs were associated with the A6 and A9, respectively. These significant SNPs almost coincided with the significant SNPs for TA.

SNP annotation and promising genes associated with eggshell ultrastructure
The annotation of significant SNPs using the Variant Effect Predictor (VEP) tool supplied by Ensembl could help identify promising genes associated with eggshell crystal structure. Detailed information about the genes is summarized in Additional file 1: Table S2. In total, 74 potential candidate genes all with significant quantitative trait loci (QTLs) were identified. The results of enrichment analysis were listed in Additional file 1: Table S3. We found that six genes were significantly enriched in the blood circulation (GO:0008015): CACNA1C (voltage-dependent L-type calcium channel subunit alpha-1 C), IAPP (islet amyloid polypeptide), ITPR2 (inositol 1,4,5-trisphosphate receptor type 2), PLCZ1 (1-phosphatidylinositol 4,5-bisphosphate phosphodiesterase zeta-1), ABCC9 (ATP-binding cassette sub-family C member 9) and KCNJ8 (ATP-sensitive inward rectifier potassium channel 8). Among them, the CACNA1C, IAPP, ITPR2 and PLCZ1 also played considerable roles in the biological process of regulating cytosolic calcium ion concentration. We also found 26 genes that were  However, no information on the interaction network is available for these candidate genes.  copy of the effect allele at the rs314759160 site causing a 0.345 SD/allele decrease in TA. The phenotypic differences between the genotypes of those SNPs were shown in boxplots (Fig. 4), revealing that homozygotes of the effect allele and of the alternative alleles possessed the lowest and highest phenotypic values, respectively, whereas the values for the heterozygotes were intermediate. Data also showed that the most significant QTL associated with TA was the rs315989578 at the position 63.8 Mb on chromosome 1 (p-value 1.01e-11). However, this SNP located in the intergenic region, explained 6.45 % of the phenotypic variation (Additional file 1: Table S2).

Discussion
The eggshell of chicken is a bio-ceramic material comprising columnar calcite (CaCO 3 ) crystals and organic proteinaceous matrix [7]. Previous studies have found that the size and orientation of calcium carbonate crystals influenced the eggshell quality [17,18]. In the  present study, we used X-ray diffraction to measure the traits of calcium carbonate crystals in the eggshell. This method was based on the analysis of X-ray diffraction patterns formed by intact eggshells, with the intensity of the spots displayed in these patterns being directly related to the size of crystals in the eggshells [12]. Our results showed that the ESS, ESW, EST, ET and MT increased as the size of the crystals increased (Fig. 1), which was consistent with previous studies [18,21]. We speculated that in eggshells with larger crystals a faster crystal growth rate had occurred, causing the deposition of more calcites in the eggshell and thus increasing its strength, weight and thickness [23]. Moreover, the palisade layer was mainly composed of packed calcite crystals, being governed by the geometry of the nucleation sites, crystal growth rate and interactions between the growing crystals [24]. The thickness and strength of the eggshell would increase as the size of the calcite crystal increases [18]. Interestingly, our results have shown that the TA was highly negatively correlated with the OD (r = − 0.48, Table 2). This indicated that the crystal structure of the larger crystals was more randomly oriented, which was consistent with a previous study [18]. Table 3 showed that TA was highly genetically correlated with the microstructure traits ET, MT, ESS and EST, results consistent with those of Dunn's study [18]. We also noted that the genetic correlations were high between the A7, A8, ET and EST. However, this high correlation has not been reported in other studies. X-ray diffraction studies have showed that the ET was oriented to a few crystal planes [5]. So, we hypothesized that the arrangement of some specific crystal planes (such as A7 and A8) may have an important effect on the eggshell structure. However, further study is needed for the crystallographic structure mechanism. In the present study, OD was not correlated with the quality of the chicken eggshells. The estimated heritability of OD (h 2 = 0.06) was also low with a large estimation error, but was smaller than previous estimates [18], partly because of the higher variation coefficient for OD ( Table 2). The genetic correlations between OD and ET in the present study were high (r = 0.48) ( Table 3), and Dunn et al. also found a similar genetic correlation (r = 0.32) between OD and ET [18], but the estimation errors were large in both studies.
To determine whether this estimate is reliable, our experiments will need to be repeated with a larger sample. Therefore, further studies are needed to explore the mechanical aspects of crystallographic texture.
The eggshell is composed of four layers, with ET accounting for more than 70 % of the total thickness [25]. Therefore, changes in the effective layer of eggshells can significantly affect their properties [26]. In a previous study, the phenotypic correlation between the EST and ET of eggshells was 0.92, while the genetic correlation was 0.96 [22]. The association analysis of QTLs in the region between 57.3 and 71.4 Mb of GGA1 was significantly associated with ESS, ESW and EST [27]. We have also identified a 9-Mb QTL (GGA1: 59.4-68.5 Mb) in this region which was associated with the ET [22]. These results have indicated that this region of chromosome 1 had a great influence on the microstructure of eggshells, but the details are unclear. Interestingly, the present study has found that this region was also closely related to the TA. One reason may be that there was a strong correlation (r = 0.5) between the ET and TA (Fig. 1). However, no QTLs were identified as being associated with crystal growth orientation. We speculated that the size of the crystal was the dominant influence on the structure of the eggshell, which was independent of the orientation of crystal growth. Our results have also shown that GGA1 was closely related to the ultrastructure of the chicken eggshell.
The formation of eggshells is mainly achieved by biomineralization in the chicken uterus [16,28,29]. Two systems have been proposed as playing an important role in the biomineralization of eggshells: the ion transport system and the matrix protein system [30]. The ion transport system provides sufficient ions for eggshell formation and matrix proteins play a regulatory role in the process of eggshell mineralization [31,32]. In a previous study, ABCC9, ITPR2, KCNJ8 and WNK1 have been mapped as being associated with ET and EST [22], and all these genes are also involved in ion transport with their expression detected in the chicken uterus. Many studies have also reported that the eggshell matrix proteins were associated with EST [7,21]. Genetic associations for crystal measurements showed that OVAL and ovo-transferrin were associated with TA, while OC116 and OCX-32 were associated with OD [18]. In the present study, we have identified 76 candidate genes including ABCC9, ITPR2, KCNJ8 as being associated with crystal size and six genes with the process of blood circulation (Additional file 1: Table S3). ABCC9 and KCNJ8 encoded membrane-associated receptors that combined with the potassium channel have been reported to become smooth muscle ATP-sensitive potassium channels in the heart [33]. Inositol 1,4,5 trisphosphate receptors (ITPRs), a family of endoplasmic reticulum Ca 2+ channels, are essential for controlling intracellular Ca 2+ levels in virtually every type of mammalian cell [34]. ITPR2 has been found to be highly expressed in cardiomyocytes and to be associated with human heart failure [35,36]. As important ion transmembrane channels, these genes are reported to play a vital role in cardiac contraction [37]. It has been reported that the calcium required for eggshell calcification mainly depends on blood circulation, and more efficient calcium transport could promote the calcification of eggshells [38]. Our results thus suggested that the strength of cardiac contractility may affect the mineralization of eggshells.
Several genes related to intracellular calcium concentration were also identified. ITPR2 was also directly involved in the process of Ca 2+ signal transduction and released during bone growth [39]. Ca 2+ oscillation is mediated by the inositol 1,4,5-trisphosphate receptors 2 and 3 (ITPR2 and ITPR3) in the osteoclast differentiation [40]. Except for previously reported major candidate genes, we have also identified several new genes associated with crystal size, CACNA1C, IAPP and PLCZ1, which influence the function of regulating cytosolic calcium ion concentration. CACNA1C encoded the protein, voltage-dependent L-type calcium channel subunit alpha-1 C, which is essential for normal blood pressure regulation through the contraction of arterial smooth muscle cells [41]. CACNA1C is highly expressed in the brain, heart, jejunum, ovary, pancreatic beta-cells and vascular smooth muscle [41][42][43]. Ca 2+ influx through CACNA1C has also been reported to activate osteogenic transcriptional programs and promote mineralization [44]. We therefore propose that CACNA1C might play an important role in regulating the concentration of calcium in uterine fluid. Amylin is encoded by IAPP and is a member of the calcitonin family of hormones that is co-secreted with insulin by the pancreatic beta cells [45]. Amylin has been reported to use the calcitonin receptor to inhibit osteoclastogenesis in vivo [45]. There was no direct evidence that PLCZ1 participated in the release of calcium but it has been reported that it indirectly affected the transmembrane transport of calcium [46,47]. Interestingly, we also found that 22 candidate genes, including ITPR2, PLCZ1 and CACNA1C, were regulated by the transcription factor, Oct-1, which acted as downstream of notch signaling during radial glia formation. Oct-1 has been reported as forming a novel complex with the transcription factor (Runx2) to regulate the expression of the mammary gland-specific gene, betacasein. Runx2 is essential for expressing several bonespecific genes and is primarily considered as a master regulator of bone development [48]. This suggested that Oct-1 might also be related to the development of bone. The present study has shown that several genes related to bone growth were significantly associated with the strength of eggshell crystals. Except for eggshell formation, the biomineralization process in birds occurs in the bone calcification [49]. Some eggshell matrix proteins such as OC-116 have been found to be expressed in bone [50], which indicated that there may be similar biological pathways in the two mineralization processes. Our results further support this hypothesis. It is also worth noting that no known eggshell matrix protein was found in these candidate genes, because some matrix proteins such as ovocleidin- 17 have not yet been annotated in the reference genome. It will thus be necessary to obtain more information on the crystallization system for eggshells.

Conclusions
The integral intensity and degree of orientation of eggshell crystal were determined, then used for a GWA analysis with a 600 K high-density SNP array. The result showed that 621 genome-wide significant SNPs ranging from 55.6 to 69.1 Mb in GGA1 were associated with TA. Ultimately, six genes, PLCZ1, ABCC9, ITPR2, KCNJ8, CACNA1C and IAPP, were considered as promising genes that were implicated in the integral intensity of eggshell crystals and could improve eggshell quality. The present study has revealed the genetic basis of the crystallographic ultrastructure of eggshells and improved our understanding of its mechanical properties. These new findings will also provide a useful theoretical basis on biomineralization in chicken eggshells.

Experimental population
White Leghorn (WL) and Dongxiang (DX) chickens, representing a standard breed and a Chinese indigenous strain, respectively, were used to construct an F2 resource population for the present study. Six WL and six DX males were mated with 133 DX and 80 WL females, respectively. Then, 25 cocks and 407 hens from the WL/ DX cross and 24 cocks and 235 hens from DX/WL cross in the F1 generation were selected to produce the F2 generation. A total of 3749 chickens from the F2 population were reared in individual cages with ad libitum access to feed. Finally, 839 hens from 49 half-sib families and 365 full-sib families with sufficient phenotypic and pedigree information were chosen for SNP genotyping. All the experimental animals were hens and came from the same batch.

Phenotypic measurements
Eggs were collected from hens at 66 weeks of age on three successive days to ensure at least one egg per hen. The egg contents were removed then the eggshells were cleaned with tap water and dried naturally at room temperature. Pieces (approximately 3 × 5 mm) of eggshell were taken at the equator. The shells were measured using X-ray diffraction (R-Axis Spider XRD, Rigaku Corp., Tokyo, Japan). The measurement conditions were as follows: Mo Kα (λ = 0.71075 Å), 50 KV, 80 mA, a pin-hole collimator (0.8 mm in diameter) and an exposure time of 30 s per frame. The samples were mounted with the inner surface of the eggshell facing the incident X-ray beam. The beam passed through the sample and the transmission diffraction pattern was recorded on the detector. The software THCLXPD attached to the XRD was used to analyze the diffraction patterns by measuring the intensity of the reflection spots.
The integrated intensities were determined from each eggshell sample using nine different rings. These rings were oriented with a set of (hkl) crystallographic planes: 113, 104, 012, 110, 202, 018, 122, 208 and 300. The data from the nine reflection peaks (A1-A9) were added to give the TA to lower the data variability. The crystal size was estimated from the intensity of the reflection spots [51]. To determine the OD, a calcite powder standard (Joint Committee on Powder Diffraction File, PDF 05-0586) having a completely random orientation was used as a reference pattern for the crystals. The ratio between the integrated intensities of the calcite reflections from the eggshell and those of the standard sample (PDF 05-0586) provided a measure of the relative orientation of the eggshell crystals [9]. In the present study, the strongest integrated intensities of the standard sample (PDF 05-0586) were 113 and 104 (Additional file 1: Table S1), indicating that the crystals were preferentially oriented to these planes. The ratio (R std ) of the calcite powder standard sample (A 113 /A 104 ) was 0.18. Thus, the OD of eggshell crystals was calculated using Eq. (1): The value of OD represented the degree of the preferential orientation of the crystals. For a value of one, the sample consists of randomly oriented crystals and for a value > 1, the sample consists of highlyoriented crystals [17].
The data on eggshell quality traits (ESS, EST, ESW, MT and ET) from our previous publications was used in the correlation analysis [22,52]. The descriptive statistics of all the phenotypic data were calculated using R (version 3.3.1) (https://www.rproject.org/).

Genotyping and quality control (QC)
All blood samples were collected from brachial veins of chickens by standard venipuncture. Genomic DNA was extracted from whole blood samples using standard phenol-chloroform method and the 839 chickens were genotyped using the Affymetrix 600 K chicken SNP chip (Affymetrix Inc. Santa Clara, CA, USA). Array Power Tools (APT) (version 1.16.0; Affymetrix Inc) software was used for quality control and genotype calling. SNPs with a minor allele frequency of < 5 % and a Hardy-Weinberg equilibrium test with P < 1 × 10 − 6 were removed using the PLINK version 1.90 package. The BEAGLE version 4.0 procedure [53] was used to impute sporadic missing genotypes, and SNPs with an imputation quality score of R 2 > 0.5 were retained for the next analysis step. A total of 839 samples and 435,867 SNPs finally remained for the subsequent GWAS.

GWA analysis
The loci-based analysis was performed using the generalized linear mixed model implemented in GEMMA [54], where the kinship matrix is calculated using the standardization method. The mixed model was based mainly on the additive effect of sites: Where: Y is the vector phenotypes of eggshell; µ the overall mean; X the covariance matrix (containing the first ten PCA principal components obtained from analysis of the population substructure); b the estimator vector of fixed effects; Z is an incidence matrix associating u with y. u the additive polygenic effect (0, GV g ), with G the genomic kinship matrix, and the additive effect variance Vg; S the design matrix containing the corresponding SNP sites; α the substitution effect size corresponding to each site; and e the vector of random residual effects (N (0, IV e )), with I the identity matrix and the residual variance (V e ).
The genomic inflation factor (λ) was calculated using the R package, qqman (Version 0.1.2) [55]. Multiple test thresholds were calculated using the simpleM method [56]. Stepwise conditional analyses were conducted using the generalized linear mixed model of the R package with the strongest SNP identified being added at each step as a covariate. A total of 58,358 valid inspections was obtained. The genome-wide significance level was 8.57 × 10 − 7 (0.05/58,358), with a suggestive significance level of 1.71 × 10 − 5 (1.00/58,358). Principal component analysis was performed using the EIGENSOFT software (v.2.04) with the results showing that there was stratification in the experimental population. In order to correct the population stratification, we selected the top ten PCA [57] components as covariates in the GWAS analysis.
The SNP-based heritability of eggshell crystal traits, calculated using the residual maximum likelihood (REML) method in GCTA (version 1.24) [58], was estimated from the same mixed model as used in the GWAS analysis. The genetic correlations were estimated from a standard bivariate linear mixed model [59] using the Bivariate GREML analysis of GCTA [60]. We used Visscher's method [61] to calculate the power of genetic correlation. We calculated the phenotypic variance contribution (CPV) of those genome-wide significant SNPs based on the equation 2pqβ 2 =σ 2 , with allele frequencies p and q, β was corresponding effect size of SNP identified in association study, σ 2 was the phenotypic variance [62].

Linkage disequilibrium analysis and gene identification
We performed a series of linkage disequilibrium (LD) analyses to characterize the causative SNPs within strong LD regions by applying the solid spine algorithm in Haploview software version 4.2 [63]. The significant SNPs were functionally annotated, then the candidate genes related to the significant SNPs or genomic regions were identified using the Variant Effect Predictor [64]. The Biomart tools were supported by Ensembl based on the Galgal6 assembly.

Functional annotation
Gene Ontology annotation and KEGG pathway were completed using the Metascape website to assign genes with corresponding terms [65,66]. The identities of all the candidate genes were converted into those of the homologous genes in humans.