2021 Volume 58 Issue 1 Pages 5-11
The aim of this study was to identify genes involved in comb development to provide insights into the molecular mechanism of chickens' comb formation. Fixation index (FST) and average number of base differences (π) of males with large and small combs were calculated based on whole-genome resequencing data. Chromosome regions with larger FST values and smaller π were considered candidate selection regions. Through further annotation of gene functions and pathways, we sought to screen possible selected genes associated with comb development. By screening whole genome resequencing data, FST and π were calculated using a 40 Kb sliding window strategy and eight regions were identified. Quantitative trait loci (QTL; FOX1 gene) related to comb length were found on chromosome 1. QTL (GLP1R, BTBD9, MIR6633, and MDGA1 genes) related to comb weight were found on chromosome 3. QTL (ALDH1A1, TMC1, and ANXA1 genes) associated with comb area were found on the Z chromosome. Nineteen genes, Wnt signaling pathway and neuroactive ligand-receptor interaction signaling pathway directly or indirectly related to comb growth and development were found through functional annotation and GO analysis. Among the selected genes LYN, GLP1R, FOX1, TBK1, STRAP, ST6GALNAC, and Wnt signaling pathways were related to immunity. MDGA1, BTBD9, MTSS1, SrGAPs, and neuroactive ligand receptor interaction signaling pathways related to neural function were screened. ALDH1A1, ANXAl, THBS, HIF-1α, and ACTN1 genes were related to heat dissipation. Among the selected genes FOX1, MDGAl, and ANXAl associated with immunity, neurological function, and heat dissipation function coincided with genes affecting the length, weight, and area of the comb. Comprehensive analysis suggested that comb development was due to multiple genes and signaling pathways.
Many avenues for selection in poultry breeding provide insights into the future selection of elite progeny. Among these, the comb is an important and reliable parameter for selection within a variety or breed. A longer redder comb in red jungle fowl males is preferentially chosen by females (Johnsen and Zuk, 1996). Comb color and size in grouse are indicators of T-cell-mediated immune responses (Mougeot, 2008) and the comb is generally an indicator of health status, disease resistance, and nutritional status in flocks of red jungle fowl (Zuk et al., 1990; Johnson et al., 1993; Mukhtar and Khan, 2012). Also, the comb is associated with reproductive performance and bone mass (Wright et al., 2008). Comb size and comb weight are associated with total egg production, egg number, and mean egg weight in chickens (Rydmel, 2010). In China, comb development is an important factor associated with consumer acceptance and price of live and frozen chickens. In breeding programs, it is common to eliminate the chickens with smaller combs and select those with larger combs. Therefore, there is a strong natural and artificial selection for comb size.
Selection leads to specific changes in the patterns of variation among selected loci and in the neutral loci linked to them. These genomic footprints of selection, known as selection signatures, and can be used to identify loci that have been subjected to selection. Based on genome-wide resequencing data, a large number of selected segments associated with domestication, fitness, and importance have been reported. Among these studies, research subjects include chicken (Elferink et al., 2012; Guo et al., 2016; Zhang et al., 2017; Lillie et al., 2018).
Wright's fixation index (FST) was initially developed by Lewontin and Krakauer (1973) to scan for signatures of selection. Negative selection usually decreases FST, while positive selection causes it to increase (Barreiro et al., 2008). Genomic regions under different selection or selection pressure eventually develop obvious genetic differences. Genomic regions with high levels of genetic differentiation can be identified by comparing FST values between distinct populations (Vitalis et al., 2001; Guo et al., 2016; Lillie et al., 2018).
Nucleotide variation, in the form of π (Tajima, 1983), can be used to quantify the degree of genetic differentiation among populations based on differences in allele frequencies. Both FST and π statistics are useful for the detection of selection signatures (Lillie et al., 2018).
This study was designed to compare male chicken with large and small combs. Based on whole-genome resequencing data and a public database, the selection signal for comb development was detected by comparing FST and π values between two groups. By comparing the genetic differentiation between two groups we expect to identify the genes related to comb development. This research could be useful in understanding the underlying mechanisms implicated in comb development and would contribute to more efficient breeding of chickens with larger combs.
All animal experiments were approved by the Animal Care and Use Committee at the Institute of Poultry Science, Chinese Academy of Agricultural Science (Approval ID: S20160605) and conducted at the institute. All of the experiments followed relevant guidelines and regulations set by the Ministry of Agriculture of the People's Republic of China.
SamplingAll males used in the present study were sourced from a Partridge Shank chicken breeding line with single comb in Poultry Science, Chinese Academy of Agricultural Science. All chickens were produced under the same conditions and raised using a standardized feeding method with free access to water. To obtain males with different comb sizes at 9 weeks of age for sequencing, comb areas and heights of 200 males were measured using the software for automatic measurement of comb development parameters (registration number: 2015SR071422). The means and standard deviations of comb areas and comb heights was 2291.90 ±481.30 mm2 and 23.06±6.61 mm, respectively.
Based on those results, 22 males with bigger combs (BC) and 19 males with smaller combs(SC) were selected for further analysis. The thresholds for the BC group were comb areas of more than 2000 mm2 and comb height of more than 34 mm, while those for the SC group were comb areas less than 1600 mm2 and comb height of less than 20 mm (Table 1).
| Sample | BC group | SC group | ||
|---|---|---|---|---|
| means | SD | means | SD | |
| Comb areas(mm2) | 2736.9 | 218.95 | 1846.89 | 166.22 |
| Comb heights(mm) | 40.34 | 3.63 | 19.65 | 1.75 |
Note: BC: Big comb; SC: Small comb
Blood samples were obtained from the wing vein. Genomic DNA was extracted from blood samples using the standard phenol/chloroform method. From total population (n=200), the bigger comb population (BC; n=22) and the smaller comb population (SC; n=19) were extracted before pooling the DNA of each population. Two pooled samples (BC and SC) were separately sequenced using a next generation sequencing technique on the Illumina HiSeq 2000 platform. A total of 140.4 Gb of pair-end reads of 100 bp were obtained. The following quality control criteria were implemented by removing 1) any reads with ≥ 10% unidentified nucleotides(N); 2) any reads with >10 nt aligned to the adapter sequence, allowing ≤ 10% mismatches; 3) any reads with > 50% bases having phred quality < 5; 4) putative PCR duplicates generated by PCR amplification in the library construction process. A total of 133.6 Gb high quality paired-end reads was mapped to the chicken genome (Gallus_gallus-5.0, http://hgdownload.soe.ucsc.edu/downloads.html#chicken) using BWA (Li and Durbin, 2009). Depth of reads of BC and SC populations were 55.41× and 53.2×, respectively, and > 90% coverage of the chicken genome was retained for single nucleotide polymorphism (SNP) calling via the Genome Analysis Toolkit (GATK; McKenna A et al., 2010). NGS data has been submitted to SRA in NCBI (SRR10303810 and SRR10303811). The BioProject accession number is PRJNA576318 and the Biosample accession number is SAMN12990250.
Analysis of Selection SignaturesTwo population genetic indices were estimated from the sequence data. Firstly, nucleotide diversity (π) was used to measure the degree of polymorphism within the population. Secondly, FST was determined as a measure of population differentiation due to genetic structure. In our analyses, π was calculated across 40,000 bp windows for each population pool using Popoolation (v1.2.2; Kofler et al., 2011a) with the following parameter set: window size of 10 K, step size of 5 K, minimum allele count of 2, minimum base quality of 20, minimum coverage of 10 and maximum coverage of 200, and the log value of π ratios was calculated. FST was calculated using PoPoolation2 (v1.201; Kofler et al., 2011b), with same parameter set as for π. Samtools (Li et al., 2009) was used to generate input files for PoPoolation 2. FST was calculated using the standard equation according to the principles of population genetics (Hartl and Clark, 2007). FST=Pi_total−Pi_within/ Pi_within, where Pi_within=(Pi_population1+Pi_population2)/2, and Pi=1−fA2−fT2−fC2−fG2 with fN being the frequency of nucleotide N (i.e. A, T, C or G), Pi_total is the total Pi for which allele frequencies in both populations are averages and Pi is calculated as above. The FST values were Z-transformed as follows : Z(FST)=(FST−µFST)/σFST, where µFST is the mean FST, and σFST is the standard deviation of FST (Karlsson et al., 2007). Putatively selected regions were located in fully overlapping windows with an extremely low π value (top 5% level) and extremely high Z(FST) values(top 5% level).
Candidate Gene Annotation for Selected Regions, GO Analysis, and Pathway AnalysisWe annotated genes in selected genomic regions using the chicken genome as a reference via the UCSC genome browser (http://genome.ucsc.edu/). The genes putatively under selection were submitted to David (http://david.ncifcrf.gov/) for enrichment analysis of the Gene Ontology (GO) and KEGG pathways. In addition, the chicken quantitative trait loci (QTL) database (https://www.animalgenome.org/cgi-bin/QTLdb/GG/index) was systematically tested, and analyzed to detect whether comb height, comb area, comb weight, comb length, early maturity, body weight growth, and other traits were enriched in selected regions.
The average Z(FST) was 2.202 with a standard deviation of 0.614. It can be seen from the Z(FST) that the populations with large and small sized combs exhibited genetic differentiations in some chromosomal regions, with the most significant interval located on chromosome 15 (see Fig. 1a). The genomic nucleotide diversity (π) of the chicken population was determined using the same size sliding window strategy. The average π was 1.533, and the standard deviation was 0.414. Whether these intervals were related to comb size was determined by comparing the overlapping areas between the smaller interval of π and the larger interval of Z (FST; Fig. 1b). In order to select the extreme candidate selected region, the overlapping region of the interval of Z (FST) ≥ 5 and the interval of 200 with the smallest negative value of π was selected, resulting in eight candidate regions, including 19 genes.
The distribution of Z (FST) and π on chromosome. Note: a. The positive end of the Z (FST) plotted on chromosomes. A dashed horizontal line indicates the cut-off Z (FST)=5. b. The negative end of π plotted on chromosomes. A dashed horizontal line indicates the cut-off for 200 smallest π.
To verify that the candidate regions were associated with comb length, comb thickness, comb area, and comb weight, the known QTL contained in the eight candidate selected regions were searched for in the Chicken QTL database (https://www.animalgenome.org). QTL (FOX1 gene) related to the length of the comb were found on chromosome 1, QTL (GLP1R, BTBD9, MIR6633, and MDGA1 genes) related to comb weight were found on chromosome 3, and QTL (ALDH1A1, TMC1, and ANXA1 genes) associated with the area of the comb were found on the Z chromosome (Table 2). These results indicated that the QTL associated with the comb development was enriched in the candidate region and also suggested that the method for detecting the selection signal was reliable.
| Chromosome | Gene | Region | π | Z (FST) | Comb QTL |
|---|---|---|---|---|---|
| 1 | FOXO1 | 170581722–170644997 | 1.1705 | 1.6818 | Comb length |
| 3 | GLP1R | 29377290–29456632 | 1.1655 | 3.6490 | Comb weight |
| 3 | BTBD9 | 29594123–29704099 | 1.3132 | 1.9223 | Comb weight |
| 3 | MIR6633 | 80493162–80493272 | 1.4758 | 2.2259 | Comb weight |
| 3 | MDGA1 | 30023079–30150365 | 1.1434 | 1.7914 | Comb weight |
| Z | ALDH1A1 | 35894470–35947479 | 1.8474 | 2.1396 | Comb area |
| Z | TMC1 | 35808065–35866858 | 1.9209 | 2.1256 | Comb area |
| Z | ANXA1 | 36026282–36040351 | 2.2242 | 2.0813 | Comb area |
A total of 26 genes were identified by genetic annotation of candidate regions using the online UCSC genome browser tool (Fig. 1). To clarify the function of these genes, PANTHER was used for functional annotation and GO analysis. Twenty genes were directly or indirectly related to comb growth and development (Table 3). Most of these genes are involved in cell growth and development or important basic signaling pathways. For example, LYN, SRGAP1, THBS1, WWC1, and MTSS1 are involved in cell migration and receptor binding. MDGA1, ST6GALNAC5, PLA2G4A, ADAM10, ZDHHC21, ACO1 are involved in Golgi action. The growth and development of the comb is inseparable from cell migration, receptor binding, and the role of the Golgi apparatus. The main function of the Golgi apparatus is to process, classify, and package proteins synthesized by the endoplasmic reticulum, which are then sent to specific parts of the cell or secreted outside the cell.
| Chromosome | Gene | Biological process |
|---|---|---|
| 1 | LYN, FOXO1 | cell migration (GO:0016477), receptor binding (GO:0005103), Golgi apparatus (GO:0005798) |
| 1 | SRGAP1 | cell migration (GO:0016478) |
| 1 | STRAP | negative regulation of epithelial cell migration (GO:0010633), receptor binding (GO:0005104) |
| 1 | TBK1 | phosphoprotein binding (GO:0051219), |
| 2 | WNT9A | Wnt signaling pathway (GO:0016056), irismorphogenesis (GO:0061073) |
| 2 | MTSS1 | receptor binding (GO:0005102) |
| 3 | MDGA1 | Golgi apparatus (GO:0005800) |
| 5 | THBS1 | cell migration (GO:0016477), peptide cross-linking (GO:0018150) |
| 5 | HIF1A | irismorphogenesis (GO:0061072) |
| 5 | ACTN1 | phosphoprotein binding (GO:0005794) |
| 8 | ST6GALNAC5 | Golgi apparatus (GO:0005794) |
| 8 | PLA2G4A | Golgi apparatus (GO:0005794) |
| 10 | ZDHHC21, ADAM10, ACO1 | Golgi apparatus (GO:0005794) |
| 11 | WWOX | Wnt signaling pathway (GO:0016057) |
| 20 | CSNK2A | Wnt signaling pathway (GO:0016055) |
| Z | ANXA1 | peptide cross-linking (GO:0018149) |
GO enrichment analysis showed that LYN, WNT9A, WWOX, and CSNK2A were involved in Wnt signaling pathway. The Wnt signaling pathway is a complex network of protein action that, while commonly functional in embryonic development and cancer, is also involved in the normal physiological processes of adult animals. The opening or closing of this pathway controls the expression of a large number of growth and metabolism-related genes (Xu et al., 2015).
The KEGG pathway analysis of these candidate genes was performed using David 6.7, and the neuroactive ligand-receptor interaction (gga04080) signaling pathway was enriched. Five genes, AGTR1, TACR3, TRPV1, PTH1R, and GLP1R, participate in this signal path.
A total of eight selection regions were identified in this study, including a number of QTL affecting comb length, area, and weight. QTL (FOX1 gene) related to comb length were found on chromosome 1, QTLs (GLP1R, BTBD9, MIR6633, MDGA1 genes) related to comb weight were found on chromosome 3, and QTLs associated with comb area (ALDH1A1, TMC1, ANXA1 genes) were found on the Z chromosome. Functional annotation and GO analysis were performed, and 19 genes were directly or indirectly related to comb development.
Among the selected genes LYN, GLP1R, FOX1, TBK1, STRAP, ST6GALNAC, and Wnt signaling pathways were related to immunity. MDGAl, BTBD9, MTSS1, SrGAPs, and neuroactive ligand receptor interaction signaling pathways related to neural function were screened. Genes related to heat dissipation function (ALDH1A1, ANXAl, THBS, HIF-1α, and ACTN1) were included in the selected regions. Among these 19 genes, some genes related to immunity and heat dissipation function, FOXO1, MDGAl, and ANXAl, coincide with genes affecting comb weight and area.
GO analysis showed that LYN, WNT9A, WWOX, and CSNK2A are involved in the Wnt signaling pathway. Its members have high homology and play important roles in regulating embryonic development, participating in cell proliferation and differentiation, as well as migration and polarity.
Among the genes involved in the Wnt signaling pathway, tyrosine protein kinase (LYN) may be involved in the body's immunity. The expression of this gene is regulated by the methylation of the promoter CpG site. The high methylation of this gene is important for maintaining normal cell growth and proliferation. LYN has an indispensable role in immunoglobulin-mediated signaling, particularly in establishing B cell tolerance (Hibbs et al., 1995).
As a transcriptional regulator, FOXO1 plays an important regulatory role in animal growth and development through transcriptional regulation and signal transcription. FOXO1, as a potential factor, has an important relationship with skeletal muscle differentiation and myotube formation in myoblast differentiation (Yang et al., 2008). Through genome-wide association analysis, four SNPs in the upstream and downstream flanking regions of the chicken FOXO1 gene were significantly correlated with growth traits such as body weight and average daily gain (Wang, 2012). In our study, the FOXO1 gene was identified to have a potential role in comb length QTL on chromosome 1.
The comb structure is divided into three layers: epidermis, dermis, and central layer. The epidermis is divided into stratum corneum and germinal layer. Merck's cells are contained in the germinal layer. Merkel cells establish prominent connections with nerve endings and play a key role in sensing external stimuli (Hitchcock et al., 2014). In the epidermis, collagen fiber bundles, adipose tissue, many larger arteries, and some nerves are arranged vertically. Our researches showed that MAM domain glycosylphosphatidylinositol anchor 1 (MDGA1) and neuroactive ligand receptor interaction signaling pathways were related to neural function. MDGA1 genes were identified to have a potential role in comb weight.
MDGA1, a unique cell surface glycoprotein, is similar to Ig-containing cell adhesion molecules that influence neuronal migration and process outgrowth. MDGA1 acts cell autonomously to control the migration of MDGA1-expressing superficial layer cortical neurons (Takeuchi and Dennis, 2006). The MDGA1 gene affects the development of the nervous system by changing the neuronal migration pathway and participates in the etiology and pathogenesis of schizophrenia (Gangwar et al., 2017).
The KEGG pathway analysis enriched the neuroactive ligand-receptor interaction signaling pathway, and the five genes (AGTR1, TACR3, TRPV1, PTH1R, and GLP1R) observed here were involved in this signaling pathway. The neuroactive ligand-receptor interaction signaling pathway is a collection of all receptors and ligands on the plasma membrane that are associated with intracellular and extracellular signaling pathways. Therefore, this signaling pathway is interesting because there are not only differentially expressed genes in this pathway, but are also physiologically the most closely related to neurological function.
The structure of the comb and the blood vessels distributed in the comb play an important role in the heat dissipation process of the chicken. The larger the comb area, the better the heat dissipation in a warmer climate, while small combs are beneficial in a cooler climate. Our research showed that the genes (ALDH1A1 and ANXA1) selected for the comb area are involved in heat dissipation. ALDH1A1 is increased in obese patients, negatively regulating their heat production. ALDH1A1 ping cells are called A1dh1a14 thermogenic fat cells. In vivo and in vitro experiments, all ALDH1A1 enzyme activities are coordinated with retinoic acid (RA) production (Reichert et al., 2011). Annexin A1 (also known as Annexin I, ANXA1) isa widely-functioning calcium and phospholipid binding protein, which participates in many cell activities such as anti-inflammatory reaction, cell differentiation and proliferation, cell death signal regulation, and apoptotic cell phagocytosis (Perretti and Solito, 2004; Sun et al., 2016), inhibiting cell migration to the site of inflammation (Damazo et al., 2011).
Based on genome-wide resequencing data, combined with population genetic differentiation index FST and population genomic nucleotide diversity (π) for system selection signal detection, a total of eight selection regions were identified. QTL (FOX1 gene) related to the length of the comb were found on chromosome 1. QTL (GLP1R, BTBD9, MIR6633, MDGA1 genes) related to comb weight were found on chromosome 3. QTL (ALDH1A1, TMC1, ANXA1 genes) associated with the area of the comb were found on the Z chromosome. Nineteen genes directly or indirectly related to comb growth and development were found through functional annotation and GO analysis. Comb development is affected by multiple genes and signaling pathways. Among the selected genes FOX1, MDGA1, and ANXA1 associated with immunity, neurological function, and heat dissipation function coincided with genes affecting the length, weight, and area of the comb.
We certify that there is no conflict of interest with any financial and personal relationships with other people or organization regarding the material discussed in the manuscript.
The current research was funded in part by the China Agriculture Research System (CARS-41-G02), National Natural Science Foundation of China (31402032), Jiangsu Province major agricultural new varieties creation project (PZCZ201728), Jiangsu Provincial Agricultural three new Project (JATS [2018]249), and the Agricultural Science and Technology Independent Innovation Funds (CX (18)3003).
