Genome-Wide Association Study of Egg-Laying Traits and Egg Quality in LingKun Chickens

Egg production is the most important trait of laying hens. To identify molecular markers and candidate genes associated with egg production and quality, such as body weight at first oviposition (BWF), the number of eggs produced in 500 days (EN500), egg weight (EW), egg shell thickness (EST), egg shell strength (ESS), and Haugh unit (HU), a genome-wide analysis was performed in 266 LingKun Chickens. The results showed that thirty-seven single nucleotide polymorphisms (SNPs) were associated with all traits (p < 9.47 × 10−8, Bonferroni correction). These SNPs were located in close proximity to or within the sequence of the thirteen candidate genes, such as Galanin And GMAP Prepropeptide (GAL), Centromere Protein (CENPF), Glypican 2 (GPC2), Phosphatidylethanolamine N-Methyltransferase (PEMT), Transcription Factor AP-2 Delta (TFAP2D), and Carboxypeptidase Q (CPQ) gene related to egg-laying and Solute Carrier Family 5 Member 7 (SLC5A7), Neurocalcin Delta (NCALD), Proteasome 20S Subunit Beta 2 (PSMB2), Slit Guidance Ligand 3 (SLIT3), and Tubulin Tyrosine Ligase Like 7 (TTLL7) genes related to egg quality. Interestingly, one of the genes involved in bone formation (SLIT3) was identified as a candidate gene for ESS. Our candidate genes and SNPs associated with egg-laying traits were significant for molecular breeding of egg-laying traits and egg quality in LingKun chickens.


INTRODUCTION
It is well-known that egg-laying traits and egg quality are very important economic characteristics in the chicken industry and are greatly associated with both the breeding value and retail egg value. Currently, the release of the chicken quantitative trait locus (QTL) contains 16,271 QTLs from 367 publications, which represent 442 different traits (https://www.animalgenome.org/cgi-bin/QTLdb/ GG/index). However, despite a range of studies in this area, wide confidence intervals (CIs) for the positions of QTL remain that have rarely been replicated. A new research era was initiated with advances in sequencing technology and single nucleotide polymorphism (SNP) chip, and genomewide association study (GWAS) has become one of the most effective methods to detect genetic variation in livestock. In previous studies, two SNPs related to body weight at first oviposition (BWF) were located at 78.8 and 78. 9 Mb, respectively, on the GGA4 gene by GWAS of Jinghai Yellow chicken, and both were located in FAM184B gene (1). Liu et al. (2) performed a GWAS in 1,078 Rhode Island Red hens and found that the 117.87-118. 36 Mb region of the GGA1 was associated with the number of eggs produced (EN) in hens 37-72 weeks that include POLA1, PDK3, PRDX4, and APOO genes. Liu et al. (3) discovered that the GTF2A1 and CLSPN genes may be candidate genes associated with EN by influencing ovarian and uterine functions. For egg quality traits, Wolc et al. (4) reported that a 90 kb genomic region (169. 42-169.51 Mb) in GGA1 is associated with egg weight (EW) at 28-66 weeks of age, two promising genes (DLEU7 and MIR15A) can be mapped to this narrow significant region and may affect EW in a pleiotropic manner. Sun et al. (5) discovered a genomic region spanning from 57.3 to 71. 4 Mb in GGA1 that is significantly associated with eggshell quality and includes ITPR2, PIK3C2G, and NCAPG genes. In addition, previous studies have shown that a genomic region spanning from 8.95 to 9.31 Mb (≈0.36 Mb) on GGA13 was significantly associated with the Haugh unit (HU), Meanwhile, MSX2 and DRD1 genes were related to embryo development and egg production (6).
The LingKun chicken is a Chinese indigenous chicken breed, which has the characteristics of large body, early sexual maturity, good meat quality, disease resistance, especially excellent egg production performance, and good egg quality. However, the genetic mechanism of these dominant traits has not been wellanalyzed, so the market for Lingkun chicken has not been fully tapped, just partially reared in the Wenzhou area of Zhejiang Province. However, there are no reports on GWAS research on laying traits and egg quality of Lingkun chickens. In the present study, a GWAS was performed on 266 LingKun Chickens aged 72 weeks to identify molecular markers and candidate genes affecting reproductive traits by using wholegenome resequencing. Results could potentially benefit research into chicken reproductive traits. These promising loci and genes could be helpful to engineer practical breeding programs for the improvement of egg quality for old hens to meet the need of prolonging the laying cycle.

Egg Production and Quality of Lingkun Chickens
The basic statistics for egg-laying traits and egg quality in the Lingkun chicken for the GWAS are shown in Table 1.
There are six traits that include BWF, the number of eggs produced in 500 days (EN500), EW, egg shell thickness (EST), egg shell strength (ESS), and HU. Heritability analysis showed that eggshell strength exhibited the highest heritability, with a value of 0.51, while BWF showed the lowest heritability, with a value of 0.05. EW and EN500 showed moderate heritability, with values of 0.37 and 0.38, respectively. Correlation analysis between traits is shown in Table 2, EST and EW are significantly correlated (p = 0.192), and ESS and EST are significantly correlated (p = 0.555). Meanwhile, EST and ESS were negatively correlated with EN500. Meanwhile, Manhattan diagram (A) and QQ diagram (B) of correlation analysis results are shown in Figure 1B. It showed  that all the traits in this study had significant correlation SNPs, and the selection of analysis model was reasonable.

EST and ESS
Two significant SNPs were associated with EST, the most associated SNP C4415750T (p-value = 2.81E-08) was located on chromosome 23 at position 4,415,750 bp, which was found within the upstream of the Proteasome 20S Subunit Beta 2 (PSMB2) gene. A total of four SNPs (A5358928G, A5360528G, T5362672G, and G5362697C) were significantly associated with ESS, which were located within the 3,769 kb segment (between 5,358,928 and 5,362,697 bp) on chromosome 13, which were found within the Slit Guidance Ligand 3 (SLIT3) gene.

Haugh Unit
The associated SNP A17272906G (p-value = 5.09E-08) was located on chromosome 13 at position 17,272,906 bp, which was found within the downstream of Tubulin Tyrosine Ligase Like 7 (TTLL7) gene.

DISCUSSION
There has been an effective use of GWAS in egg-laying traits and egg quality traits in other species and narrow regions or SNPs associated with duck production trait have been revealed (7,8).
Here, we present a GWAS of egg-laying traits and egg quality in a LingKun chicken population from China. A total of 37 significant SNP loci were detected in this study. Although EW, ESS, ESS, and EST traits were significantly correlated, GWAS results did not find a common QTL.

Body Weight at First Oviposition
The BWF reflects the constitution of the reserved chicken, the greater the BWF, the greater the energy reserve of the chicken, the greater the EW, and the stronger the ability to continue laying eggs during the laying period. Hens begin to lay eggs after sexual maturity and are still in the development stage, and their weight can still increase by 30-40 g per week. After 20 weeks, about 40 weeks of age, the growth and hair breeding were basically stopped. Xie et al. (9) performed a GWAS with F2 chickens and found that the region 169-179 Mb on GGA1 was significantly related to 23 growth traits. Jin et al. (10) found that 18 SNPs reached 5% Bonferroni genome-wide significance with growth traits in Yancheng chickens, and these SNPs, which were located on four different chromosomes and in a region of 72.3-82.1 Mb on GGA4, had a significant effect on growth traits. In the present study, growth-related genes, GAL, CENPF, and GPC2, obesityrelated gene, PEMT, and CYP3A4, P450, TMEM200B genes were found to be associated with BWF. Growth and development are affected by the regulation level of hormones. The GAL gene encodes a neuroendocrine peptide, small neuropeptides that regulate a variety of physiological functions, the release of growth hormone and insulin, and adrenal gland secretion (11). CENPF has a potential role in regulating cell differentiation in skeletal myogenesis and embryogenesis (12)(13)(14). CENPFdepleted cells show delayed mitosis and reduced stability of centromeric microtubules (14,15). It is well-known that changes in body weight are related to the mitotic differentiation of cells. GPC2 belongs to the GPC (GPC1-6) family, which plays a variety of roles in growth factor signaling and cancer cell growth (16,17), regulates the growth and development process of the body. PEMT gene activity involves a number of physiological processes, such as lipid flow between the liver and plasma and the delivery of essential fatty acids to the blood and peripheral tissues through liver-derived lipoproteins. The PEMT gene encodes an enzyme that converts phosphatidylethanolamine to phosphatidylcholine through sequential methylation in the liver. PEMT deficiency reduces choline availability (18). In studies on childhood obesity genes, the PEMT gene was found to be associated with polyunsaturated fatty acids in human erythrocytes (19). In conclusion, de novo synthesis of choline through PEMT plays an important role in regulating systemic energy metabolism, and the PEMT gene may be a pharmacological target for the treatment of obesity and insulin resistance. We believe that PEMT has an indirect effect on BWF by influencing choline regulation. In addition, the CYP3A4 gene encodes a member of the cytochrome P450 superfamily of enzymes. The cytochrome P450 proteins are monooxygenases that catalyze many reactions involved in drug metabolism and synthesis of cholesterol, steroids, and other lipids. TMEM200B is a transmembrane protein overexpressed in the endometrium.

The Number of Eggs Produced in 500 Days
Regions associated with numbers of eggs produced or egg production rate have been detected on many chromosomes, while most of these regions were population specific, only the regions on GGA2, 4, and 5 were detected in multiple studies (20,21 (22) found a genome-wide significant locus (ss1985401199) located on the sex chromosome Z, which was associated with egg number in hens at 25-45 weeks. In the present study, four proximal genes (TFAP2D, DCDC2, NRSN1, and CPQ) were found, although most have not been previously reported in chicken. TFAP2D gene is the most diverse member of AP-2 family and plays a key role in development, such as the development of eyes, face, limbs, and neural tube, and is mainly expressed in the retina and brain during chicken embryo development (23). The 69.78 Mb (between 57.21 and 126.99) region on GGA2, which was considered as a QTL, was significantly correlated, which includes seven associated SNP loci, DCDC2 gene polymorphism or intron deletion may be involved in neural development and may affect signal transduction of primary cilia (24,25). NRSN1 gene is related to chicken brain synaptic development, which plays an important role in nerve organelle transport, nerve signal transduction, or nerve growth. Since the function of these genes in chickens is unclear, we infer they may regulate egg production by interacting with the central nervous system. In addition, the CPQ gene, which is also known as plasma glutamate carboxypeptidase (PGCP), is a secreted protein hydrolyzing the circulating peptides in the extracellular microenvironment (26) that liberates I-thyroxine (T4) from thyroglobulin (Tg) and is converted to 3,5,3 ′triiodothyronine (T3) by type I 5 ′ -deiodinase (DIO1) (27). The thyroid hormone T3 binds to thyroid hormone receptors (TRs), and after heterodimerization with retinoid, X receptor can act as a transcription factor mediating diverse physiological processes, such as embryonic development, cell differentiation, metabolism, and cell proliferation, by regulating the expression of the downstream target genes (28). The regulation of thyroid and gonad is realized through hypothalamus pituitary thyroid/ovarian axis. Therefore, it is suggested that the CPQ gene may affect egg production by regulating hormone levels.

Egg Weight
Egg weight displays a consecutive increase with the hen's age. It is known that the starting weight will affect the EW. Chickens that are fed a high-fat diet have large body weight, high liver fat content, fallopian tube weight, and EW. Liao et al. (22) discovered an SNP (ss1985401190) located on GGA4, which was significantly associated with EW. Liu et al. (4) found through GWAS that DLEU7 and MIR15A genes had a great influence on EW. In the present study, seven genes (SLC5A7, NCALD, UPP1, FHOD3, ENAH, CEP350, and RNF2) containing near SNPs are associated with EW. The biological process of the SLC5A7 gene includes affecting the development of uterine embryos. Its pathway includes the transport of glucose and other sugars, bile salts and organic acids, metal ions, and amine compounds mediated by transmembrane transport. These substances participate in metabolism and affect EW by affecting the formation of egg yolk. Meanwhile, the SLC5A7 gene responds to the depolarization of the cholinergic terminal by increasing the protein density in the synaptic constitution membrane, and choline is transported from the cell to the synaptic terminal to synthesize acetylcholine, which can mediate the synthesis of acetylcholine in cholinergic neurons. Sato et al. (29) showed that osteoblasts express specific acetylcholine receptors and cholinergic components and that acetylcholine plays a possible role in regulating the proliferation and differentiation of osteoblasts. In conclusion, the SLC5A7 gene can regulate the growth and development of hens and affect EW. NCALD is involved in the calcium signal pathway and G protein-coupled receptor signal pathway. Experiments have shown that NCALD is a potential hippocampal memory-related factor related to obesity (30). Total fat intake has an important impact on the low-density lipoprotein (LDL) peak particle diameter (LDL-PPD); LDL will lead to the increase of lipid levels in hens. When hyperlipidemia is high, it will cause insufficient blood supply to the placenta, thus affecting the laying weight. Rudkowska et al. (31) showed that NCALD with dietary fat intake influences the variation in the LDL-PPD. UPP1 gene is a dimeric enzyme and plays an essential role in pyrimidine salvage and regulation of uridine homeostasis. FHOD3 is a sarcomeric protein highly expressed in the cardiac tissue and required for the assembly of the contractile apparatus. Members of the ENAH gene family are involved in actin movement involved in a range of processes dependent on cytoskeleton remodeling and cell polarity. ENAH can activate mitogen-activated protein kinase (MAPK) and activate the p38mapk signaling pathway. The product of the CEP350 gene is a large protein with a cytoskeleton-associated protein-glycinerich domain (Cap-Gly), usually found in cytoskeleton-associated proteins. CEP350 gene was found to be expressed before oocyte maturation and disappeared after maturation. It is suspected to be related to the follicular formation (32). The protein encoded by the RNF2 gene is one of the Protein Crystal Growth (PcG) proteins, which is important for the transcriptional inhibition of various genes involved in the development and cell proliferation. Studies in mice have shown that the gene is involved in cell proliferation during early embryonic development. They regulate various developmental processes by promoting cell proliferation, inhibiting differentiation, and blocking cell fate regulation through many signaling pathways. It is speculated that these genes may affect EW during embryogenesis and organogenesis.

EST and ESS
Eggshell quality is an important index to reflect the damage resistance rate of eggs and is related to egg freshness and hatching rate. During the process of eggshell formation, uncalcified eggs are bathed in uterine fluid that plays regulatory role in eggshell calcification. The difference in protein abundance plays essential roles in influencing eggshell strength. Sun et al. (33) found that 15 proteins mainly related to eggshell matrix-specific proteins, calcium binding and transportation, protein folding and sorting, bone development or diseases, and thyroid hormone activity were considered to have a closer association with the formation of strong eggshell. A lots of QTL regions affecting eggshell quality were discovered by previous linkage studies and mainly distributed on GGA 1-9 and GGA11 (32). In the present study, the promising genomic region had no overlap with the previously reported regions; PSMB2 is a candidate gene for EST revealed by the present GWA analysis. PSMB2 is one of the fourteen 20 S proteasome subunits. Previous proteomic screening showed that there were a large number of binding proteins in eggshell matrix; the 20 S proteasome is the catalytic core of the 26 S proteasome complex that represents the major component of the adenosine triphosphate (ATP)/ubiquitin-dependent protein degradation pathway. This pathway regulates many processes in the cells, such as progress through the cell cycle, gene transcription, and metabolic pathway (34). PSMB2 has been found in the bone marrow of human and plays a role in regulating homologous recombination (HR) and DNA double-strand break (DSB) repair in human cells (35)(36)(37)(38). Studies in the pigs showed that proteasome inhibitors blocked the exit of maturing pig oocytes from the myocardial infarction (MI) stage (39). PSMB2 gene transcription had a positive effect on the meiosis of bovine oocytes (34). The above results show that the PSMB2 gene had an effect on follicular development and protein homeostasis, but its effect on eggshell quality is not clear. The results of the GWAS for the ESS trait showed that four SNP markers were located in a 3.77 kb region on GGA13. SLIT3 is one of the three SLIT ligands identified in the central nervous system, expressed and involved in muscle development and bone formation (40), angiogenesis (41), and stem cell migration (42), which are widely expressed in various cell types during embryogenesis and are involved in basic aspects of limb development. Knockdown of individual SLIT genes or Slit1 and Slit2 in combination resulted in a shortening of the length of the humerus (43). Kim et al. (44) found that β-catenin knockdown completely blocked Slit3-stimulated osteoblast migration and proliferation. The SLIT/ROBO pathway is associated with the development of regrading ovarian follicles in hens through autocrine or paracrine mode. Overexpression of SLIT3 can significantly reduce mRNA and protein expression levels of follicle-stimulating hormone (FSHR), growth and differentiation factor 9 (GDF9), osteogenesis acute regulatory protein (STAR), and cytochrome P450 11A1 (CYP11A1) in GC. The proliferation of primordial granulosa cells was significantly reduced. These results suggest that SLIT3 may inhibit the proliferation, differentiation, and follicular selection of granulosa cells (45). Meanwhile, the SLIT-ROBO pathway is also believed to play an important role in human angiogenesis and vascular system function. Reproductive tissues, such as ovaries, endometrium, and placenta, are rich in blood vessels and have dynamic and tightly regulated angiogenesis and remodeling. Eggshells are formed in the womb, and the rich blood vessels facilitate material exchange. Calcium for eggshell formation is provided by feed and bone tissue. The laying process of laying hens is a process of structural bone destruction. We hypothesized that SLIT3 regulates eggshell strength by affecting the proliferation, differentiation, absorption of bone cells, and changing the calcium content in bone tissue. All the above pieces of evidence suggested SLIT3 should be a crucial and promising candidate gene relating to eggshell quality.

Haugh Unit
Haugh unit is an important index to evaluate egg protein quality and reflect egg freshness. It is calculated according to EW and egg height. Liu

CONCLUSIONS
In this study, we report GWAS analysis for six egg-laying traits and egg quality in Lingkun chickens. A total of 37 significant SNPs were identified that include GAL, CENPF, GPC2, PEMT, TFAP2D, and CPQ genes related to egg-laying and SLC5A7, NCALD, PSMB2, SLIT3, and TTLL7 genes related to egg quality. In addition, the identified region on GGA13 for ESS was a promising QTL that could be potentially applied in marker assistant selection in a breeding program to improve egg quality, which needs further study. Our results found some potential genes, which might be used to improve egg production rate and egg quality.

Animals and Data Collection
The population selected for the study is an eighth generation pure line of Lingkun chicken from Wenzhou Golden Land Development Co. Ltd. The eighth generation was formed by 84 unrelated male chickens, each of which mated with 12 female chickens. There were 266 female offsprings that were chosen randomly from the same batch of chickens. The chickens were reared in stair-step cages under consistent conditions. We recorded their egg production every day. Eggs were collected in three consecutive days when hens were 72 weeks old. The average for 3 d was taken as the phenotypic value of each trait for every hen. Individual data for egg quality traits that included EW, EST, ESS, and HU were measured by instrument (Robotmation Co., Ltd., EMT-5200; Robotmation Co., Ltd., EFG-0503; Karl Deutsch, 1061). All phenotypic values of traits in genotyped individuals were tested for normality, and some abnormal values that were extremely deviating from normal distribution were deleted before conducting association tests.

Sample Preparation
Blood samples were collected from the wing vein of 266 female offspring using heparin sodium as anticoagulants.

Variant Discovery and Quality Control
A total of 1,509,826,877 SNPs were obtained. Quality control of the genotyping data was performed in PLINK (v1.9). During this process, markers were selected based on three conditions: The SNPs that exhibited minor allele frequencies (MAFs) lower than 5%, low minor allele frequency lower than 2%, and a Hardy-Weinberg equilibrium tests lower than 1.0 × 10 −8 were rejected. Finally, 266 individuals and 9,327,106 SNPs were identified in the present study. The experiment fast structure (https://github.com/ rajanil/fastStructure) was used to perform population structure analysis. It was based on the maximum likelihood method (MLM). In addition, a principal component analysis (PCA) for the independent 9,327,106 SNPs was performed in Genomewide Complex Trait Analysis (GCTA) (version 1.24). In order to reduce the influence of population stratification, the first and second components of the PCA were used in the model.

Statistical Analysis
The general linear regression model (GLM) in PLINK was used in this study. The model is as follows: where Y is the vector of observations; X is a design matrix of covariates containing all other fixed effects; b is a vector matrix of fixed effects that includes population structure, g is a vector of SNP effects with g ∼ N 0, Aσ 2 g , A is interpreted as the genetic relationship matrix (GRM), which was calculated using the 9,327,106 independent SNPs acquired by PLINK (v1.9) software, ε is a vector of errors with ε ∼ N 0, Iσ 2 ε . The analyses were performed using TASSEL v5.2 software (version 2.1). Significance thresholds were established from the estimated number of independent SNP markers and linkage disequilibrium (LD) blocks. LD block was defined as set D to 0.8 and the maximum interlocking area length to 500 KB. Using this approach, a total of 527,732 independent SNP markers and LD blocks were found. The threshold p-value of the 5% Bonferroni genome-wide significance was calculated based on the estimated number of independent SNPs and LD blocks for autosome SNPs. The threshold p-value of the 5% Bonferroni genome-wide significance was set at 9.47E-08 (0.05/527, 732), and the threshold p-value for significance was 1.89E-06 (1/527, 732).
The gene closest to the SNP with an extremely significant correlation is considered as a candidate gene. Sequence alignment was conducted through the National Center for Biotechnology Information (NCBI) and Ensembl databases to identify genes within 100 KB upstream and downstream of significant SNPs.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: National Genomics Data Center, China National Center for Bioinformation (accession: https://bigd.big.ac.cn/gsa/browse/CRA006206).

ETHICS STATEMENT
The experiment was approved by the Animal Care Committee of Anhui Agricultural University (No. SYDW-P20190600601).

AUTHOR CONTRIBUTIONS
LL, HC, SZ, YZ, SL, CW, YT, and TZ contributed to conception and design of the study. XL organized the database. WX performed the statistical analysis. JG wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.