Genomic Regions Associated with Early/Late Feathering in Daheng Broilers Identified Using Genome-Wide Association Studies

1Farm Animal Genetic Resources Exploration and Innovation Key Laboratory of Sichuan Province, Sichuan Agricultural University, Chengdu, Sichuan 611130, China 2Animal Breeding and Genetics Key Laboratory of Sichuan Province, Sichuan Animal Science Academy, Chengdu, Sichuan 610066, China Article Information Received 02 March 2023 Revised 26 March 2023 Accepted 09 April 2023 Available online 28 July 2023 (early access)


INTRODUCTION
E F and LF are sex-linked phenotypes, and one day of age is considered the best period for automatic sex determination (Zhao et al., 2016). The alleles (K-k+) are sex genes, located on the short arm of the chicken's Z chromosome (Pal and Singh, 1997). The early and late feathers are controlled by a pair of alleles (K-k), and the K allele are incomplete dominant to k+ allele (Elferink et al., 2008). Among them, the K allele was related to LF, which retard the emergence of primary and secondary flight feathers, while the k+ allele resulting in the earliest emergence of feathers. Based on the length of primary feathers and primary-convert feathers at hatching, the EF phenotype is roughly divided into two types: more than 5mm and the difference between 2-5mm, while the classification and distribution of the LF phenotype are different. Some studies have divided it into five types, namely C (longer and isochron), D (shorter and isochron), E (large gap and isochron), F (small gap and isochron > G), but the distribution ratio has not been reported (Song et al., 2003). The establishment of early and late feathered strains by phenotypic selection or molecular identification is of great significance for chicken breeding and scientific management, which is reflected in the rapid operation, low cost and low chick stress response.
Although the LF/EF phenotype is beneficial for chick sexing, the relationship between the feather rate phenotype and performance is also being explored. Since plumage have different heat preservation ability at different temperatures, different feather types indirectly affect feed conversion rate and fat deposition, and thus affect carcass weight of poultry (Fotsa et al., 2001;Merkley and Lowe,

O n l i n e F i r s t A r t i c l e
1988). LF/EF phenotype affects egg production and quality.
Researches have shown that the egg production of brown layers at sexual maturity in LF group was lower than that in EF group (Collins et al., 1970). Compared with yellowfeathered broilers with LF phenotype, laying rate and hatching performance of EF phenotype were significantly increased (Goger et al., 2017). In addition, the lesion score indicated that the resistance of chicks to E. coli inoculation was associated with the K allele (Dunnington et al., 1986). Therefore, it is of great significance to study the mechanism of feather difference for sex identification and economic character improvement of chickens.
In recent years, sequencing technology has profoundly changed the way we understand genomes, enriching the genetic and biological inquiry into complex traits (Fuentes-Pardo and Ruzzante, 2017;Visscher et al., 2012). Based on the low cost of genotyping, high-throughput single nucleotide polymorphism (SNP) genotyping assays are increasingly used in genome-wide association studies in chickens and as a breeding tool (genome prediction and selection) (Groenen et al., 2011). GWAS is a genome-wide trait mapping method to locate genetic loci associated with traits by testing the significance of association between marker polymorphisms and phenotypic diversity of target traits (Sharma et al., 2015). In this context, Onteru et al. (2013) performed GWAS using a dense SNP panel to identify SNPs associated with FE-related traits and incorporated them into the prediction model (Onteru et al., 2013). Previous studies have identified two candidate genes associated with the chicken feather K allele: prolactin receptor (PRLR) (Bu et al., 2013) and sperm flagellum 2 (SPEF2) (Elferink et al., 2008). Both have been identified as biomarkers of chicken feathering. However, more candidate genes and more comprehensive identification methods for feathering traits need to be developed. In the actual production process, although the gender of chickens is mainly identified according to the development of wing feathers, but the success rate of identification is still worth improving.
In this study, genotyping-by-sequencing (GBS) technology was used to develop genome-wide SNP markers for early and late feathering traits in Daheng broilers (55 samples in group K and 54 samples in group M). On this basis, GWAS analysis was carried out to find out the molecular genetic markers significantly related to this trait, so as to provide theoretical support for later molecular breeding technology.

Animal samples
All the chickens used in the experiments were provided by Sichuan Daheng Poultry Breeding Co. LTD, which were raised in a traditional environment, with the same temperature and humidity and with free access to water and feed. Phenotypes of the feather rate in one-day-old chicks were distinguished according to Xie et al. (1985). Briefly, 1-day-old chicks with a primary feather more than 2 mm longer than the primary-covert feathers is an EF phenotype, otherwise defined as a LF phenotype, which can be divided into inverted, equal and slightly long types (Fig. 1). Then, blood samples were collected from a total of 99 1-day-old Daheng broilers which provided by Sichuan Daheng Poultry Breeding Co. LTD, including 55 EF phenotype chickens and 54 LF phenotype chickens, according to the method of Luo et al. (2012). In brief, chickens were sterilized with alcohol, followed by 1 ml venous blood drawn using a 2 ml medical syringe and stored immediately in a -20℃ refrigerator. Chicken early-feathering phenotype. (B) Chicken latefeathering phenotype. One-day-old chicks with primary feathers longer than primary-covert feathers by more than 2 mm were identified with the EF, otherwise, they would be identified with the LF.

DNA library construction and quality control
The DNA concentration was accurately quantified by Qubit and concentrated to 50 ng/μL (Wang et al., 2017). The entire library was prepared using a Covaris ultrasonic crushing machine with random interruptions according to Illumina's official procedures. Then, Qubit2.0 was used for initial quantification, the library was diluted to 1.5ng/ μl, and the insert size of the library was detected using Agilent 2100. Illumina NovaSeq6000 sequencing platform was used for PE150 mode sequencing; the obtained Raw data were transformed into Raw sequencing sequences by base recognition, called Raw Reads, which were filtered for quality using Fastp software (Chen et al., 2018). The following filtering parameters were adopted: Remove the reads that (i) contained adaptor and (ii) more than 3 unknown bases (N bases) and (iii) those, in which one end O n l i n e

F i r s t A r t i c l e
had more than 20% low quality base reads (sequencing quality value< 5).

SNP genotyping
The BWA software ) was used to compare Clean Reads to the chicken reference genome (GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b_ genomic.fa), and then SAM tools software  was used to convert the initial alignment results into BAM format and sort them.

Variant detection and annotation
The haplotype caller model of GATK software ) was used to detect SNP and InDel in a single sample. Finally, the genotype GVCF files of individual samples were merged using the combine variants function of GATK software, and the genotype GVCFs function was used for population SNP detection. SNP quality filtering was performed using Vcf tools (Danecek et al., 2011) and Plink software (Purcell et al., 2007). The set was then filtered using the following criteria: (i) maximum alleles number <2; (ii) mapping quality >20; (iii) depth of the variant position >3.

Genome-wide association study (GWAS)
All GWAS analyses were performed using Plink (Purcell et al., 2007) software. The SNPS of autosomes and sex chromosomes were analyzed separately using the method of previous studies (O'Sullivan et al., 1991). Since EF/LF trait was a categorical variable, the case/ control model was used to analyze the differences in the frequencies of alleles at SNP loci between groups K and M. The Bonferroni method was used to correct for multiple comparisons to obtain statistically significant differences (P < 1.32e-07).

SNPs loci detection
We designed primers to verify the three significant SNPs loci on Z chromosome, aiming to find a reference method for the molecular identification of sex. All primers sets were synthesized by Sangon Biotech Primer Design Center (Shanghai, China) (Supplementary Table SI). First generation sequencing is used for sequence determination.

Genome sequencing data analysis
A total of about 331 Gb of raw data was obtained, and the Quality score of raw reads was concentrate mainly on 36, and most of them >30 ( Fig. 2A and Supplementary  Table SII). Base composition of reads: base composition deviation occurred in the first 1-10bp; <10bp, the frequencies of the four bases were similar, and there was little difference in position (Fig. 2B). Reads without unique read alignments were discarded; the difference between A/T ratio and G/C ratio at any location was less than 20%.

SNP detection and type distribution
The rare SNPs, loci and SNPs with high sample missing rate in the population will affect the accuracy of subsequent population genetic analysis and GWAS analysis. SNPs were detected, and their quality was filtered using the Vcf tool and Plink software. In total, 11626 high quality SNPs were obtained, which were distributed in 33 autosomes and Z chromosomes. The number of SNPs per chromosome is shown in Supplementary Table SIII. The position of SNP on chromosome is shown in Figure 3, in which exons are 776, accounting for 3.35%, and introns are the most, accounting for 55.32%. Subsequently, 14 major autosomes and Z chromosomes were selected, and the density distribution of SNP in these chromosomes was analyzed in detail (Fig. 4).

Genome-wide association study
We used the Plink software (Version 1.90b3.34; Purcell et al., 2007) for the genome-wide association study. 10 significant SNPs were found in the whole genome, among which the three points with the strongest signal were all on the Z chromosome and located in the genomic regions known to be related to this trait (Fig. 5). For these 10 significant SNP loci, we studied their specific position O n l i n e

F i r s t A r t i c l e
on the chromosome, allele frequencies among different groups and other information (Supplementary Table SIV).

Filtering of differentially expressed genes and functional enrichment analysis
According to the 10 significant SNP sites selected, all candidate genes within 1M region around the SNP sites were found, and a total of 87 candidate genes were identified, of which 67 were protein-coding RNAs and 20 were long non-coding RNAs (lncRNA) ( Supplementary  Table SV). Surprisingly, no known candidate genes were found for the three SNPs on the Z chromosome.
All the candidate genes identified in this study may play a role in the phenotype of chicken feather percentage. In order to understand the function and regulatory mechanism of these genes, GO analysis (Fig. 6, Supplementary Table SVI) and KEGG pathway annotation (Table I) were performed. From GO analysis, we can see that these genes are enriched in the following processes: embryonic forelimb morphogenesis; proteasome ubiquitin independent protein catabolic process; peripheral nervous system neuron development, etc. Furthermore, Candidate genes were enriched not only in common pathways such as MAPK signaling pathway, mTOR signaling pathway, and p53 signaling pathway, but also in C-type lectin receptor signaling pathway and Neuroactive ligand-receptor interaction. Fig. 6. The GO analysis with candidate genes between early-feathering and late-feathering Daheng broilers.

Bayesian model for predicting chicken feather pattern
Firstly, each candidate SNP and population frequency information were determined by referring to the sequencing results (Supplementary Table SVII). Then, the probability of being LF under a certain observed genotype was estimated by Bayesian model: Here, p(G) is the probability of a certain genotype in the whole population, p(S) is the proportion of LF in the population (p (S)=0.5), and p (G│S) is the probability of a certain genotype in the LF. The statistical results of different loci are shown in Table II. Multiple loci can be combined for joint estimation, which provides reference for identification of early and late feather molecules.

DISCUSSION
As an important genetic trait in modern poultry breeding, EF/LF phenotype is widely used by poultry breeders all over the world. Not only is it commonly used for sex determination, but it is also associated with growth rate, egg production and disease resistance (O'Sullivan et al., 1991). In addition, the integrity and maturity of feathers is an important factor affecting the carcass beauty of local chickens, especially broilers. At present, it is generally believed that this trait is a pair genetic quality trait, which is controlled by four complex alleles Kn, Ks, K and k on Z chromosome (Dunnington et al., 1986). The late feather (K) is 1dominant to the early feather (k), and the recessive sequence of genes is Kn >Ks >K >k (Somes, 1969). The duplicated gene dPRLR-SPEF2 is localized in the short arm of Z chromosome, which can be used as the molecular basis of EF phenotype (Elferink et al., 2008). Studies have shown that there is no close linkage between it and ev21 endogenous virus insertion (Guo et al., 2014). Therefore, the interlacing situation and degree between late feathering line and early feathering line can be used as the discrimination standard to allow auto-sexing of chicks. PRLR is involved in the regulation of the development, metabolism and circulation of skin tissues and their appendages (mammary gland, sebaceous gland, hair, etc.) (Craven et al., 2001); and its regulatory role in the growth cycle of animal coat and seasonal hair moult has been widely studied. Similarly, SPEF2 has recently attracted more and more attention as a candidate gene for EF/LF phenotype due to its partial sequence duplication with PRLR (Elferink et al., 2008). However, feather growth is not only associated with the Z chromosome alone, but also affected by related genes on autosomes. In 1946, Jones identified a gene on an autochromosome that controls wing feather growth in chicks and named it tarded, denoted by t. While, many other genes on autosomes that regulate feather growth rate need to be identified and explored. By comparing the genomes of late-feathered chickens and early-feathered chickens, we identified 11626 high-quality SNPs distributed in 33 autosomes and Z chromosomes, which provides a new complement for the study of candidate SNPs related to chicken feather growth.
In this study, we used the Case/Control model to analyze whether there was a significant difference in gene frequencies between the EF and LF phenotype groups at SNP loci. Ten significant SNP loci were found in the whole genome, among which the three points with the strongest signal were all on the Z chromosome and located in the genomic regions known to be related to this trait. The chicken genetic linkage map shows that the K locus is located on the short arm of the Z chromosome and thus can be used for automatic sex identification at hatching, even in the embryonic stage (Hamoen et al., 2001). However, a search of the 1M region around the three SNPS on the Z chromosome yielded no known candidate genes. This may be due to differences in sequencing depth or filtering methods. Interestingly, we found more emerging candidate genes on autosomes that may influence feather development. Among them, Davoodi et al. (2022) identified TBC1D22A located on chromosome 1 as a candidate gene for feather color through 60 k high-density SNP array analysis (Davoodi et al., 2022). Feather color is not only considered a social signal for chickens, but also a breeding recognition tool for breeders. In addition, it has been shown that PRSS23, as a key candidate gene for two different feather follicles, may be involved in the regulation of feather morphology (Xu et al., 2022). Miraculously, among our candidate genes, HOXD family genes were clustered in large numbers, almost including HOXD3-HOXD13. Vertebrate Hox genes encode transcription factors that function in the development of a variety of organs and structures (Schep et al., 2016). Primarily, the HOX gene regulates the development of skin derivatives, such as mammary glands in mammals (Awgulewitsch, 2003;Chuong et al., 1990) and hair follicles and feathers in poultry (Carroll and Capecchi, 2015). Early studies have suggested that the expression of CHOXC-8 and CHOXD-13 in chick skin and skin appendage specifications may be involved in dorsal feather and foot scale formation in chick embryos (Kanzler et al., 1997). Alasdair I Reid found that Hoxb-4, Hoxa-7 and Hoxc-8 were expressed temporally and spatially in a collinear manner in the skin (Reid and Gaunt, 2002). Thus, the HOX gene provides regionally restricted location cues during embryonic skin modeling and also conveys general signals required for feather morphology, which may have a significant impact on transcriptional activity in the developing chick skin.
In this study, we detected multiple pathways that have been confirmed to be related to feather or hair production and development, such as Wnt signaling pathway (Xie et al., 2020), mTOR signaling pathway (Gao et al., 2019), p53 signaling pathway (Qiu et al., 2022), and MAPK signalling pathway (Fang et al., 2018). These pathways not only affect feather development, but also regulate embryonic development, epithelial appendage formation, and hair follicle repair and regeneration. Additionally, we also identified some pathways that have not yet been identified but may also be involved in avian feather development. Metabolic pathways were determined to be associated with the anaerobic degradation of feather keratin and coloration of chicken feathers; Sphingolipid metabolism also acts on quantitative changes in feather color (Wu et al., 2022;Davoodi et al., 2022). Overall, we identified a number of genes and signaling pathways that are associated with or potentially play a role in feather pattern in poultry.
Molecular marker-assisted selection (MAS) is mainly used for breeding new varieties. However, the early identification of sex is equally important for accelerating the speed of breed improvement and commercialization. Bayesian prediction allow to estimate the effect of molecular markers through differential contraction, which can describe the relationship between complex molecular markers and target traits more flexibly (Gianola et al., 2009). Studies have used Bayesian models to identify the relationship between heart-type fatty acid binding proteins and fatty acid metabolism in chickens (Li et al., 2010). In addition, Bayesian models have practical applications O n l i n e

F i r s t A r t i c l e
Genomic Regions Associated with Early/Late Feathering in Daheng Broilers 7 in poultry epigenetics and stress (Videla Rodriguez et al., 2022). Here we sequenced three significant SNPs on the Z chromosome using next-generation sequencing; identified candidate SNPs and population frequency information. Roughly, an individual with GC genotype at chrZ: 11263583 locus was found in the population, with 74.19% certainty that the chicken was the LF; 100% of CC genotype was EF trait.

CONCLUSION
In conclusion, a total of 11626 high-quality SNPs were identified. On this basis, GWAS analysis was carried out to identify molecular genetic markers significantly related to this trait, screen out candidate genes that may be related to feather type, and search for biological pathways to play a role. Therefore, our results may provide a genomic and molecular evidence that autosomes, in addition to sex chromosomes, also have an impact on emergence rate. The combination of population gene frequency information and Bayesian probability model with apparent traits is beneficial to improve the accuracy of early poultry sex discrimination.

IRB approval
All animal experiments were approved by the Institutional Animal Care and Use Committee of Sichuan Agricultural University (Certification No. SYXK2019-187).

Ethical statement
The animals and experimental methods used in the experiments were studied in accordance with the guidelines of the Ministry of Agriculture of China. The procedures followed in animal experiments in this study were approved by the Animal Care and Use Committee of Sichuan Academy of Animal Sciences. Strive to minimize the suffering of animals.
All animal experiments were approved by the Institutional Animal Care and Use Committee of Sichuan Agricultural University (Certification No. SYXK2019-187).

Supplementary material
There is supplementary material associated with this article. Access the material online at: https://dx.doi. org/10.17582/journal.pjz/20230302070341

Statement of conflict of interest
The authors have declared no conflict of interest.