Genomic Analysis Reveals Pleiotropic Alleles at EDN3 and BMP7 Involved in Chicken Comb Color and Egg Production

Artificial selection is often associated with numerous changes in seemingly unrelated phenotypic traits. The genetic mechanisms of correlated phenotypes probably involve pleiotropy or linkage of genes related to such phenotypes. Dongxiang blue-shelled chicken, an indigenous chicken breed of China, has segregated significantly for the dermal hyperpigmentation phenotype. Two lines of the chicken have been divergently selected with respect to comb color for over 20 generations. The red comb line chicken produces significantly higher number of eggs than the dark comb line chicken. The objective of this study was to explore potential mechanisms involved in the relationship between comb color and egg production among chickens. Based on the genome-wide association study results, we identified a genomic region on chromosome 20 involving EDN3 and BMP7, which is associated with hyperpigmentation of chicken comb. Further analyses by selection signatures in the two divergent lines revealed that several candidate genes, including EDN3, BMP7, BPIFB3, and PCK1, closely located on chromosome 20 are involved in the development of neural crest cell and reproductive system. The two genes EDN3 and BMP7 have known roles in regulating both ovarian function and melanogenesis, indicating the pleiotropic effect on hyperpigmentation and egg production in blue-shelled chickens. Association analysis for egg production confirmed the pleiotropic effect of selected loci identified by selection signatures. The study provides insights into phenotypic evolution due to genetic variation across the genome. The information might be useful in the current breeding efforts to develop improved breeds for egg production.


INTRODUCTION
Domestication is the strongest type of directional selection, resulting in the accumulation of numerous phenotypic variations in animals. During the last two decades, numerous genetic bases of traits accounting for the causative mutation of interesting phenotypes in domestic animals have been deciphered (Andersson and Georges, 2004). In livestock, it has been reported that vibrant appearance not only represents the physiognomy but also relates to development, health, and production traits. Such changes in unrelated phenotypes are potentially caused by pleiotropy or linkage between independent genetic architectures (Wright et al., 2010). Pleiotropy refers to a phenomenon in which a single gene or genetic variant affects multiple apparently unrelated phenotypic traits (Stearns, 2010). Several types of pleiotropy, including biological pleiotropy and spurious pleiotropy, have been defined (Hodgkin, 1998;Solovieff et al., 2013). Biological pleiotropy is described as a gene or genetic variant having a direct biological effect on more than one phenotypic trait, and the different causal variants are tagged by the same gene or genetic variant. Spurious pleiotropy is where the causal variants in different genes are in strong linkage disequilibrium (Hodgkin, 1998;Solovieff et al., 2013). Genome-wide association studies (GWASs) on human complex diseases (Solovieff et al., 2013;Pickrell et al., 2016) and livestock traits have revealed that pleiotropy is common in vertebrates (Rosengren Pielberg et al., 2008;Johnsson et al., 2012;Curik et al., 2013;Johnsson et al., 2014).
The comb of birds is a sexual ornament that is used as an indicator of sexual maturity and also related to mate choice in wild populations (Mukhtar and Khan, 2012). In female chicken, the size of comb is correlated to fecundity; for instance, domestic layers with large combs produce more eggs, which is potentially regulated by hormones (Wright et al., 2012), whereas, in male chicken, the comb size is inversely correlated to sperm viability (Navara et al., 2012). The type of comb is associated with sperm quality; for example, roosters with rose combs carrying homozygous allele have defective sperm motility (Crawford and Smyth, 1964;Imsland et al., 2012). The color of comb is also linked with sperm viability; for instance, it was reported that males with redder combs have higher percent of active sperms (Navara et al., 2012). In general, combs are colored red, yellow, or orange due to the blood engorged in the superficial dermis, which is richly vascularized, and deposition of carotenoid pigments in the epidermis (Stettenheim, 2000). However, several breeds, for instance, the well-known silky chicken, present dark or black combs owing to the heavy dermal melanization of the internal tissue, which is known as fibromelanosis (FM). Silky chicken is an extensively studied breed displaying FM. It has been reported that FM is caused by inverted duplication and combination of two genomic duplications separated by a >400-Kb region on chromosome 20. The endothelin 3 (EDN3) gene, which encodes a potent mitogen in melanocytes, is involved in the first duplication and is considered the causal gene of FM (Dorshorst et al., 2011;Shinomiya et al., 2012).
China has a wide range of chicken breeds with various phenotypes. Dongxiang blue-shelled chicken (Gallus gallus) is a specific indigenous breed that originated from Jiangxi Province, South China. Blue egg is the main characteristic of Dongxiang blue-shelled chicken, and the associated gene has been fixed, and the genetic mechanism of blue eggshell color has been deciphered . Besides the blue eggshell trait, another dermal hyperpigmentation trait has also been segregated in Dongxiang blue-shelled chicken population. There are two different comb colors in the chicken, namely, dark and red (Figure 1). Almost all the individuals with dark comb present a unique collection of dark phenotypes such as dark wattle, face, beak, skin, and muscle, whereas the red comb chickens have white or normal colored organs. The visible dark phenotypes in livestock are due to intense pigmentation of the dermal layer of the skin. Based on the comb color, two chicken lines were generated for over 20 generations since the early 2000, namely, dark comb line (DCL) and red comb line (RCL). Simultaneously, a mixed population was kept containing both red comb chickens and dark comb chickens and mated randomly during breeding, resulting in the random line (RAL). Except for the comb appearance, there are no obvious differences between the two chicken lines. As mentioned, the RCL has higher egg production than that of the DCL. Such a difference is found not only in just one breed, but in another Chinese local breed, Dehua black chicken also has two types of comb color. Egg production in red comb chickens is approximately 25% higher than that in dark comb chickens (Hongtao, 2006).
The study suggests that dark comb phenotype is associated with a high incidence of low egg production in chicken. Moreover, we speculate that a potential genetic pleiotropism might link the comb color and egg production. Therefore, in this study, we aimed to explore potential mechanisms involved in the relationship between comb color and egg production among chickens.

Ethical Statement
Regular quarantine inspection of the experimental farm in China Agricultural University was conducted. The blood samples were collected from brachial veins of chickens by standard venipuncture. The whole procedure was carried out strictly following the protocol approved by the Animal Care Committee of China Agricultural University.

Animals and Phenotype
The samples utilized in the study were from three lines of chicken in the experimental farm of China Agricultural University. A total of 178 female chickens were obtained from the following three lines: 86 chickens were from RCL, 30 chickens from DCL, and 62 chickens from RAL. In the RAL, 14 were red comb chickens and 48 were dark comb chickens ( Table 1). The number of eggs from weeks 17 to 60 was recorded for these chickens. The number of eggs was also recorded in another large population (N = ~500) from a conservation farm of Dongxiang blue-shelled chicken. Frozen section of comb tissues (N = 6 for DCL and RCL) were prepared and stained with hematoxylin and eosin following the protocol of a previous study (Han et al., 2015). The blood parameters were determined using an automated hematology analyzer and the ferroheme kit (Sino-uk Institute, Beijing, China).

Genotyping and Quality Control
The genomic DNA was extracted from blood samples using the avian blood DNA extraction kit (Tiangen Co., LTD., Beijing, China). Genotyping was carried out using the Affymetrix Axiom 600K Chicken Genotyping Array (Affymetrix, Inc., Santa Clara, CA, USA) with a total of 580,691 single nucleotide polymorphisms (SNPs) (Kranis et al., 2013); all genome coordinates are representative of the Gallus_gallus-5.0 reference assembly. The first-step quality control and genotype calling from the raw data in the form of CEL files were performed using Affymetrix Power Tools v1.16.0 software with the criterion of dish quality control of >0.82 and call rate of >97%. The second-step quality control was then implemented with plink v1.9 software (Purcell et al., 2007). The SNPs with unknown genomic positions and redundant coordinates were removed. The SNPs with a call rate of <90% and a minor allele frequency of <5% were excluded. The SNPs that deviated from Hardy-Weinberg equilibrium (P < 1E-6) were removed. Five individuals (one in DCL and four in RAL) with a call rate of <90% were removed. Finally, 173 chickens and 387,704 SNPs remained for the further analysis.

Analysis of Population Structure
Prior to the population structure analysis, all autosomal SNPs were pruned using the indep-pairwise option with a window size of 25 SNPs, step of five SNPs, and r 2 threshold of 0.2. After pruning, 36,853 independent SNPs remained. The principal component analysis (PCA) was carried out using GCTA software (Yang et al., 2011) to assess population structure using these independent SNPs.

Genome Wide Association Analysis
We conducted an association analysis of comb color phenotype with genotypic polymorphism using a univariate linear mixed model in Gemma software (Zhou and Stephens, 2012). Simultaneously, potential cryptic relatedness and hidden population stratification were corrected in the model. We estimated the centered relatedness matrix using independent SNPs implemented in Gemma software prior to the association analysis. The univariate linear mixed model for each SNP marker is as follows: where, y is an n-vector of binary comb color trait for all RAL individuals; W is an n × c matrix of covariates (fixed effects) including a column of 1s; α is a c-vector of the corresponding coefficients including the intercept; x is an n-vector of SNPs genotypes tested in association analysis, β is the effect size of the SNP; and u is an n-vector of random effects with a covariance structure as u ~ MVN n (0; λτ −1 K) (where τ −1 is the variance of the residual errors, λ is the ratio between the two variance components, and K is a known n × n relatedness matrix). MVN n denotes the n-dimensional multivariate normal distribution. ε is an n-vector of errors. We used Wald test statistic for each SNP testing.

Analysis of Signatures of Selection
Red comb chickens from RCL and dark comb chickens from DCL ( Table 1, Dataset2) were used to identify signatures of selection due to divergent selection. Two independent approaches, viz., Weir and Cockerham's F ST (Weir and Cockerham, 1984) and cross-population extended haplotype homozygosity (XP-EHH) (Sabeti et al., 2007), were implemented in the study. We scanned the whole genome to identify regions with increased genetic divergence (F ST ) between chickens with red and dark combs. We used VCFtools (v0.1.13) (Danecek et al., 2011) to calculate F ST and performed Z transformation. The genomic regions with a high ZF ST (1% quantiles) were analyzed for gene contents. The statistic of XP-EHH calculated for each SNP in the dataset was implemented in Selscan v-1.1.0b (Szpiech and Hernandez, 2014). The XP-EHH statistic was also transformed to Z distribution. Prior to the XP-EHH test, fastPHASE v1.4 software (Scheet and Stephens, 2006) was implemented to estimate missing genotypes and reconstructing haplotypes from unphased SNP genotype data. For DCL-RCL comparison, the XP-EHH results were standardized for each chromosome separately, and the 1% quantiles of the standardized XP-EHH distribution for all SNPs were used as threshold to identify outlier SNPs.

Gene Enrichment Analysis
The gene contents in the candidate regions were annotated according to the outliers tested by XP-EHH and F ST from the University of California-Santa Cruz database (galGal5). Firstly, we converted the chicken Ensembl transcript IDs to human orthologue gene IDs from Ensembl Gene 93 Database using the online tool BioMart (Kinsella et al., 2011). The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes pathway analyses were performed in the blue-shelled chickens with Metascape (http://metascape.org). Benjamini-Hochberg false discovery rate was used for correcting the P values. The candidate genes were classified into categories: molecular function, biological process (BP), and cellular component.

Association Analysis With Selected Single Nucleotide Polymorphism
The association analysis between the SNPs potentially under positive selection and the number of eggs from weeks 17 to 60 was carried out with the function -assoc in PLINK software (Purcell et al., 2007).

Verification of Duplications in Dark Comb Chickens
To test the boundaries of the two duplications causing the FM phenotype, a three primer diagnostic test for the presence of each of the duplications (Dorshorst et al., 2011) was used, following a touchdown thermal cycling protocol: 95°C for 5 min, 16 cycles of 95, 68 (−1°C/cycle), and 72°C for 30 s each, followed by 24 cycles of 95, 52, and 72°C for 30 s each.

Candidate Gene Expression
The total RNA was extracted from the ovary using Trizol reagent, followed by synthesis of cDNA from 1 μg of RNA using the reverse transcription kit (Tiangen Co., LTD., Beijing, China). Four pairs of primers designed for the four candidate genes (EDN3, BMP7, BPIFB3, and PCK1) and the housekeeping gene GAPDH are presented in Supplementary Table 1; all primers span one intron at least. The quantitative real-time polymerase chain reaction (qPCR) was performed for the four candidate genes from the DCL (N = 6) and RCL (N = 6) chickens. For the qPCR runs, the target genes and housekeeping gene were performed in three technical replicates, and the GAPDH gene was used as the internal control for normalizing the level of expression of the four target genes. The average relative quantitative values were calculated, and the differences in expression levels were analyzed by Student's t-test between DCL and RCL.

Difference in Egg Production Between Red Comb and Dark Comb Chickens
We calculated the number of eggs of different populations during the laying period in several generations. The results revealed that the number of eggs produced by RCL was significantly higher than that by DCL, almost throughout the laying period (P < 0.001) (Figure 2A). The number of eggs produced by the red comb chickens was significantly higher than that by the dark comb chickens in the RAL (P < 0.01) (Figure 2A). In another large population from the conservation farm of Dongxiang blue-shelled chicken, the number of eggs produced by RCL was significantly higher than that by DCL (P < 0.001) ( Figure 2B).

Dark Color of Chicken Comb Is Caused by Melanin Pigmentation
We determined the blood parameters, including the level of hemoglobin, heme, and iron (Fe) ion in the red comb and dark comb chickens. No significant difference was observed between the two chicken lines with respect to the blood parameters (Supplementary Figure 1). Furthermore, we prepared histologic section of the comb tissue obtained from the two chicken lines. We observed that melanin was deposited in the dark comb tissue and absent in the red comb tissue (Supplementary Figure 2). Our findings suggest that melanin instead of blood elements imparts the dark color of comb in chickens.

Population Structure Analysis
The PCA of 36,853 independent SNPs with r 2 < 0.2 using the first two principal components (Figure 3) showed that the three lines formed separate clusters and were well separated from each other. The samples from the RAL were closely clustered except for two outliers and one individual from DCL, probably due to the phenotype recording error, whereas the samples from RCL and DCL were slightly separated from each other, and both of them were far away from the RAL probably because of divergent selection for over 20 generations.

Identification of Genes Related to Comb Color by GWAS
To explore the genetic mechanism of comb color in Dongxiang blue-shelled chickens, we conducted a GWAS by the univariate method for comb color phenotype. As the PCA plot showed that there was a clear differentiation in the three lines ( Figure  3), we selected 58 chickens from the RAL ( Table 1, Dataset1), which were clustered together to perform the GWAS of comb color phenotype to reduce the false positive rate of association analysis induced by population stratification. A total of 138 genome-wide significant SNPs (Bonferroni adjusted) located on chromosome 20 were successfully identified for the trait (Figure 4, Table 2). The association SNPs located between positions 10,705,665 and 12,394,530 were entirely attributable to a chromosomal region (~1.7 Mb) harboring  24 genes including EDN3 and bone morphogenetic protein 7 (BMP7) (Supplementary Table 2). The most significant SNP AX-76199743 (Affymetrix Probe ID) is at position 10,971,168 on chromosome 20. This SNP and another seven SNPs with extremely low P values are close to EDN3 ( Table 2).

Detection of Signatures of Selection in Pairwise Populations of DCL and RCL
To determine whether a positive selection might have occurred recently in Dongxiang blue-shelled chicken, two independent methods F ST and XP-EHH were implemented to identify the signatures of selection between DCL and RCL. We first estimated the F ST to identify candidate genes that potentially regulate specific divergent traits in DCL and RCL. We standardized the F ST to Z scores and considered the upper tail (top 1%) of ZF ST distribution as the outlier SNPs and, thus, were considered to be potentially positively selected. There were 3,779 outlier SNPs with ZF ST of >3.5 between DCL and RCL (Figure 5A), and 34 outlier SNPs with extremely high ZF ST values (top 0.01%) were identified, which formed a potentially differentiated region on chromosome 20. A total of 813 genes were annotated according to the outliers based on the University of California-Santa Cruz database. Within the region with the highest ZF ST value, four genes, viz., PCK1, BPIFB3, EDN3, and BMP7 with important functions in physiological processes located on chromosome 20 were identified (Table 3).
We also employed the XP-EHH method to explore potential variants experiencing selection between DCL and RCL. Similar to the F ST method, we calculated and standardized XP-EHH data and considered the upper tail of Z XP-EHH distribution as mentioned previously. We identified 3,225 outlier SNPs with Z XP-EHH of >2.5 ( Figure 5B) and annotated 149 genes. The XP-EHH analysis also showed the extreme outliers with a high Z XP-EHH located on chromosome 20, comprising the four genes EDN3, BMP7, PCK1, and BPIFB3.
As expected, each approach detected regions displaying selection signatures according to the feature of genetic polymorphism data. We compared the overlap between outlier SNPs and genes in pairwise populations of DCL and RCL identified by the two methods. The region on chromosome 20 was mostly distinct from others in the genome. Of the 3,779 outliers detected by F ST sweep mapping, 305 overlapped the outliers that were detected by the XP-EHH method (Supplementary Figure 3B), and 253 of the overlapping outliers were detected by both methods are distributed on chromosome 20. Sixty-one genes were found to overlap, as detected by the two methods, which also included the four genes EDN3, BMP7, PCK1, and BPIFB3 on chromosome 20 (Supplementary Figure 3C). The distribution of potentially selected outlier SNPs in the genes was examined. Most of the outliers were located in introns, untranslated regions, and intergenic regions (Supplementary Figure 3A), suggesting that the potential target for selection was regulatory mutation.

Gene Enrichment Analysis
The potential selected regions contained 813 and 149 genes in the top 1% outliers detected by F ST and XP-EHH approaches, respectively. Analysis of gene enrichment with the clusters of genes presented significant categories potentially involved in melanin pigmentation and reproduction (Table 4). Twenty   Figure 4). The highly significant category was "regulation of protein kinase activity" (GO: 0045859) (P = 5.95E-8) containing the candidate genes PCK1 and EDN3; another significant category was "transmembrane receptor protein tyrosine kinase signaling pathway" (GO: 0007169) (P = 5.80E-9). Two significant categories "neural crest cell differentiation" (GO: 0014033) (P = 5.88E-4) and "neural crest cell migration" (GO: 0001755) (P = 1.29E-3) related to BPs of neural crest cell involved the two genes EDN3 and BMP7. We also found two significant categories related to skeletal system development, namely, "bone development" (GO: 0060348) (P = 2.32E-7) and "cartilage development" (GO: 0051216) (P = 1.46E-7), involving the genes BMP7 and bone morphogenetic protein receptor type 2 (BMPR2) ( Table 4). Among the candidate genes, we found that BMP7 and BMPR2 are also involved in the GO BPs category "reproductive system development" (GO: 0061458) (P = 6.89 E-4); BMP7 and PCK1 are involved in the "cellular response to steroid hormone stimulus" (GO: 0071383) (P = 1.92E-3) and "cellular response to hormone stimulus" (GO: 0032870) (P = 7.19E-3). The categories identified were closely related to reproduction.

Association Analysis Between Selected Loci and Egg Production
To test the potentially pleiotropic effects on egg production, we selected a cluster of SNPs underlying selection signatures identified by the independent selective sweep analysis to perform association analysis for egg production trait. Thirtyseven SNPs were identified to be significantly associated with the number of eggs. These SNPs are distributed closely on chromosome 20, which covered a 2.4-Mb region from 10,785,456 to 13,243,956 bp, including the EDN3 and BMP7 genes (Supplementary Table 3). This suggests that the genomic region has a strong effect on egg production.

Validation of Genomic Duplication in Associated Candidate Region
To verify whether dermal hyperpigmentation in the comb of Dongxiang blue-shelled chicken is caused by the inversion of genomic duplication, a multiple PCR diagnostic test was performed to detect the boundaries of the two duplications. We identified two boundaries of the two duplications with an inversion of the second duplication in DCL successfully, whereas no duplication was found in RCL (Supplementary Figure 5). We confirmed that the dark comb color in Dongxiang blue-shelled chicken is caused by the genomic rearrangement on chromosome 20.

Expression of Candidate Genes in the Ovary
To access the potential effect of the four candidate genes (EDN3, BMP7, PCK1, and BPIFB3) on egg production, the qPCR was performed using cDNA obtained from the ovary of dark comb and red comb layers. We found that EDN3 was overexpressed in the ovary of dark comb layers (P = 0.041, ~2.5-fold increase), suggesting that the genomic duplication close to the EDN3 gene act as an enhancer in regulating the function, whereas an approximate twofold decrease was observed in the expression of BMP7 in the ovary of dark comb layers, but the difference is marginally nonsignificant (P = 0.061). No significant differential expression was observed with respect to PCK1 and BPIFB3 in the ovary between DCL and RCL (Figure 6).

DISCUSSION
Long-term divergent selection has genome-wide effects on phenotypic variation, such as the body weight (Johansson et al., 2010;Lillie et al., 2018) and abdominal fat content (Zhang et al., 2012) in chickens. The chickens used in the present study have been separated for over 20 generations and subjected to divergent selection for comb color, resulting in a significant difference in the number of eggs and color of the comb, which can be an indicator of egg production. In the present study, through genomic analyses of three lines of Dongxiang blueshelled chicken, we elucidated the extent of genomic changes in the three lines resulting from selection and identified a genomic region associated with dermal hyperpigmentation in chicken, which is also involved in ovarian function. We demonstrated the pleiotropic alleles at EDN3 and BMP7 potentially explain the link between comb color and egg production. Dark appearance is widespread among livestock due to melanin migration and deposition in the dermal layer of tissues. Two loci are related to dermal hyperpigmentation, with the autosomal dominant nature of FM allele and the sex-linked inhibitor of dermal melanin (ID) locus on the Z chromosome (Dorshorst et al., 2011). A previous GWAS demonstrated that a locus of the EDN3 gene on chromosome 20 and a region on chromosome 21 are significantly associated with the comb color trait in two Swedish local chicken breeds (Johansson and Nelson, 2015). Furthermore, six SNPs on chromosome Z has been reported to segregate with the comb color phenotype in one of the breeds (Johansson and Nelson, 2015), which is inconsistent with the genomic region of chromosome Z associated with the ID locus in other breeds (Dorshorst et al., 2010). In the present study, we found that dermal hyperpigmentation in chicken comb is associated with the EDN3 gene, which is involved in genomic duplication. We also identified the BMP7 gene located on chromosome 20 with a potential role in melanin pigmentation. Previous studies have shown that EDN3 gene was associated with dermal hyperpigmentation (Dorshorst et al., 2011;Shinomiya et al., 2012), and the BMP7 gene was associated with hyperpigmentation of the visceral peritoneum (Luo et al., 2013). No other significant regions were found to be associated with comb color phenotype in the current study. This indicates that the locus of ID is homozygous in blue eggshell chickens. The different loci contributing to dark comb probably present a case of parallel selection at the molecular level under adaptive selection. As a matter of fact, the common phenotypes in different breeds are probably caused by various loci, such as frizzle feather (Guo et al., 2018), blue eggshell (Wragg et al., 2013), polydactyly (Zhang et al., 2016), and dwarfism (Wang et al., 2017) in different chicken breeds. Thus, considering the genetic architectures of different chicken populations is indispensable even though they have the same phenotypic trait.
As an important secondary sexual characteristic, the comb has a close relationship with the production and reproduction traits in chickens (Cornwallis and Birkhead, 2007;Mukhtar and Khan, 2012). In this study, we have shown that comb color was a significant predictor of egg production in Dongxiang blue-shelled chicken. The female chickens with red comb had higher number of eggs. We scanned genome-wide selective signatures in DCL and RCL and identified two genes, viz., EDN3 and BMP7, with pleiotropic effects on melanin pigmentation and reproduction, which potentially cause the link between comb color and egg production. Gene enrichment analysis showed significant categories involved in dermal hyperpigmentation and reproduction. Neural crest cell system and regulation of protein kinase activity are involved in the significant GO categories, which are associated with the process of melanocytes migration (Sommer, 2011;Shinomiya et al., 2012) and melanin synthesis (Hearing, 2011). GO categories that closely relate to bone development and cellular response to hormone stimulus are also involved. Bone metabolism has been shown to be linked with egg production, for instance, female chickens deposit calcium in the diaphysis, where it can be transferred into the eggshell (Johnsson et al., 2012;Wright et al., 2012). EDN3 and BMP7 genes are enriched in the categories, suggesting their important roles in melanin pigmentation and reproduction.
Endothelins (EDN) are 21-amino acid vasoactive peptides reported to play various roles in the reproduction system such as in steroidogenesis, folliculogenesis, and ovulation (Ervin et al., 2017). The EDN3 gene has been reported to influence the function of granulosa cells in several mammals. This can significantly suppress estrogen production induced by follicularstimulating hormone (FSH) (Calogero et al., 1998;Ervin et al., 2017). The expression of EDN3 can be regulated by gonadotropins in epithelial cells of the oviduct, particularly increased in the isthmus of this tissue (Jeoung et al., 2010). Moreover, EDN3 has been reported to play important roles in melanogenesis; it is required for the development of neural crest-derived pigment cells (like melanocytes) (Baynash et al., 1994), promoting neural crest cell proliferation and consequently increasing the number of melanocytes significantly (Lahav et al., 1996). In the present study, we found that the expression of EDN3 is increased in the ovary of dark comb layers (Figure 6). Based on the pleiotropic effects of EDN3, we inferred that the upregulated expression of EDN3 increased melanogenesis and inhibited folliculogenesis, resulting in the association between comb color and egg production. The BMP7 gene has known functions in folliculogenesis and ovulation in mammals (Lee et al., 2001;Shimasaki et al., 2004;Rossi et al., 2016) and chickens (Onagbesan et al., 2003). The BMP7 gene could increase the expression of follicle-stimulating hormone receptor (FSHR) gene in human granulosa cells and decrease the expression of luteinizing hormone receptor (LHR) gene (Shi et al., 2010). It also increases estradiol production by stimulating the activity of FSH (Shimasaki et al., 1999). In chickens, BMP7 can help follicular development by stimulating granulosa cell proliferation (Onagbesan et al., 2003). On the contrary, decreased BMP7 expression can disrupt melanocyte homeostasis in normal melanocytes and inhibit tumor growth in human uveal melanomas FIGURE 6 | Expression analysis of candidate genes. The quantitative realtime polymerase chain reaction (qPCR) analysis of four genes (EDN3, BMP7, PCK1, and BPIFB3) in the ovary tissue of dark comb and red comb chickens. The error bars indicate the standard error of the mean. *p < 0.05, ·p < 0.1. June 2019 | Volume 10 | Article 612 Frontiers in Genetics | www.frontiersin.org (Notting et al., 2007), whereas upregulation of BMP7 correlates with tumor progression (Hsu et al., 2005). In the present study, we found downregulation of BMP7 expression in the ovary of dark comb chickens, which might potentially be associated with egg production and melanin pigmentation.
During the past 10 years, GWAS have identified significant associations with many complex traits in humans and livestock; several variants affect multiple traits (Paaby and Rockman, 2013;Pickrell et al., 2016). Pleiotropy is a well-studied phenomenon with respect to a single mutation affecting multiple phenotypes. For instance, the comb mass has been shown to be linked with egg production and bone allocation in chicken. A locus containing two tightly linked genes produces pleiotropic effects on such related phenotypes (Johnsson et al., 2012). Rose-comb with homozygous loci in males has been associated with poor sperm motility. It has been reported that rose-comb is caused by a large genomic inversion, which induces ectopic expression of the MNR2 gene. Simultaneously, the genomic structure variation disrupts the CCDC108 gene located at the distal inversion breakpoint, affecting sperm motility (Imsland et al., 2012). Gray hair is associated with a high incidence of melanoma in horses. It has also been found that gray hair phenotype is caused by a 4.6-kb duplication in the STX17 gene constituting a cis-acting regulatory mutation, causing susceptibility to melanoma (Rosengren Pielberg et al., 2008). Several examples of structure variation contributing to single gene traits have previously been described in chicken (Wright et al., 2009;Wang et al., 2013;Dorshorst et al., 2015;Guo et al., 2016). This suggests that structural variation in the genome has range effects on the evolution of phenotypic diversity. Indeed, the structural genomic variations play important roles in adaptation and diversification both in animals and plants (Wellenreuther et al., 2019), for instance, a large structural variation (specifically an inversion) that harbors genes controlling color-pattern leads to an accentuated differentiation between ecotypes in stick insect (Lucek et al., 2019). As the structural variation along the genome with long fragment frequently spans more than one gene, it is plausible that the changes in genome can cause wide effects on multiple genes. We identified the genomic region encompassing two duplications, which has been reported to be the causal mutation of dermal hyperpigmentation (Dorshorst et al., 2011). We also identified two genes, namely EDN3 and BMP7, which are in or close to the region of duplication. We speculate that the structural variation probably influences the expression of genes harbored in the genomic region, although the effects of such genes must be validated through functional approaches.
In most avian species, the males with the most impressive sexual ornaments are favored by females, thereby siring more offspring than their competitors (Tobias et al., 2012). For instance, males with larger comb size were preferred by females because larger combs are a cue for females, indicating the sire quality of the male (Mukhtar and Khan, 2012). The relationship between comb color and egg production potentially provides a case of sexual selection. If females were choosing males with redder comb, then this would support phenotype-linked fertility hypothesis, which suggests a link between a male secondary sexual trait and the functional ability of its sperm (Sheldon, 1994), because comb color was positively related to sperm viability (Navara et al., 2012). Interestingly, the genetic basis of such sexual traits has been explicated in part by the mechanism of pleiotropy (Chenoweth and McGuigan, 2010). Although in the case of domesticated chicken, which is no longer subject to sexual selection, the preexisting genetic architecture of the ancestral red jungle fowl will be maintained by the current artificial selection (Johnsson et al., 2012) and could orchestrate trade-offs of sexual signals.

CONCLUSIONS
In conclusion, the results showed that there is a close relationship between comb color and egg production. The two loci on chromosome 20 (EDN3 and BMP7) potentially contributed to comb color phenotype. Further analyses of genomic data suggested that the GO terms related to melanogenesis and folliculogenesis are important targets of selection for hyperpigmentation and egg production in Dongxiang blue-shelled chickens. Both EDN3 and BMP7 genes play pleiotropic roles in hyperpigmentation and egg production in chicken. Our study provides novel insights on pleiotropy of chicken comb-related genes and necessitates further selection and breeding in egg production.

DATA AVAILABILITY STATEMENT
All the raw data of SNP genotyping for the study population were submitted to the National Center for Biotechnology Information Gene Expression Omnibus database with the Gene Expression Omnibus accession no. GSE124906.

ETHICS STATEMENT
Regular quarantine inspection of the experimental farm in China Agricultural University was conducted. The blood samples were collected from brachial veins of chickens by standard venipuncture. The whole procedure was carried out strictly following the protocol approved by the Animal Care Committee of China Agricultural University.