Genome-wide scan for selective footprints and genes related to cold tolerance in Chantecler chickens.

The Chantecler chicken, a unique Canadian indigenous breed, is well adapted to extremely cold environments. However, its genetic characteristics have not been well studied. Here, we analyzed the whole genomes of 10 Chantecler chickens and 121 worldwide chickens, which indicated that Chantecler chickens were derived from commercial chickens and exhibit a high level of inbreeding. Based on a genome-wide scan, we identified two vital candidate regions containing ME3 and ZNF536, which are related to fat metabolism and nervous system in cold adaptation, respectively. We also found that the genetic mechanism of cold adaptation in Chantecler chickens differed from that of chickens from other cold regions, such as northern China. Our study indicated that specialized commercial chickens in the early 20th century contained sufficient genetic diversity to adapt to extreme cold environments over a very short time. These findings enrich our understanding of the adaptive potential of commercial species.

INTRODUCTION

Elucidating the genetic mechanisms underlying organismal adaptation to extreme environments, such as cold climate, is of great interest. Livestock and poultry, such as cattle, goats, sheep, chickens, and ducks, which are distributed worldwide, have adapted to various environmental conditions and production systems through both natural and artificial selection. Climate-mediated selection pressure has shaped phenotypic diversity and left a genetic “imprint” on the genomes of animals raised in different geographical regions. Understanding the genetic basis of adaptation has become a central focus of evolutionary biology. Rapid progress in high-throughput technologies, computational methods, and analytical techniques has made it possible to study candidate genes and causal variations that contribute to the genetic architecture of complex phenotypes in different organisms (Ai et al., 2015; Rubin et al., 2010; Yang et al., 2016).

To understand the genetic basis of animal adaptations to the natural environment, studies have focused on various species, including humans, living in different climates, e.g., chickens (Wang et al., 2015), yaks (Qiu et al., 2012), and canids (Li et al., 2014; Miao et al., 2017; Signore et al., 2019) at high elevation (Tibetan Plateau of China), camels (Wu et al., 2014) in desert regions, and Inuits (Fumagalli et al., 2015), Siberians (Hallmark et al., 2019), and Yakutian horses (Librado et al., 2015) in cold environments. Low temperature is a major environmental factor that limits animal growth and can threaten animal survival. Animals have evolved various physiological and biochemical mechanisms to adapt to cold environments. To date, studies on humans (Fumagalli et al., 2015) and different animal species, such as sled dogs (Sinding et al., 2020), sheep (Lv et al., 2014), and goats (Cai et al., 2020), have identified several candidate genes with major effects on cold adaptation. For example, mitochondrial uncoupling protein genes UCP1 and UCP3 are related to non-shivering thermogenesis (NST) in brown adipose tissue (BAT) and skeletal muscle (Cannon & Nedergaard, 2004). Several potential candidate genes in the human genome, i.e., FADS1, FADS2, FADS3, TBX15, and WARS2, are associated with adaptation to a polyunsaturated fatty acid (PUFA)-rich diet and cold climates (Fumagalli et al., 2015). Some specific haplotype genes, including TRPC4,TRPV2,CACNA1A,SLC25A40, andAPOO, are potentially related to arctic adaptation in sled dogs (Sinding et al., 2020). Comprehensive analysis of sheep suggests that strong selection of the TBC1D12 gene could act as a key point for studying the biological processes and potential mechanisms of genetic adaptation to the environment (Lv et al., 2014). Modern goat populations from southern and northern China show two major genes, FGF5 and EDA2R, associated with hair follicle development, which may be related to local adaptive evolution (Cai et al., 2020). Although chickens are widely distributed and adapted to the harsh environmental stresses prevalent in many rural areas, there are scant published studies on their adaptation to cold climates.

The Chantecler (CA) chicken is a unique dual-purpose breed from Canada and is well adapted to the cold climate of Quebec, where the average temperature in winter ranges from about −10 to −25 ℃. They exhibit several unique cold-adapted characteristics such as dense feathers and a thick down layer against the body. Furthermore, they have a distinctive small cushion comb and wattles, which may help them withstand harsh climate conditions (Ekarius, 2007). However, the excellent adaptability of CA chickens to extreme cold environments remains poorly explored. Thus, research on the mechanisms underlying their adaptation to such environments will provide new insight into cold adaptation in poultry.

In the present study, we resequenced the whole genome of 10 CA chickens and downloaded the resequencing data of 121 different chickens distributed worldwide (Li et al., 2017; Wang et al., 2015, 2020; Yi et al., 2014) to analyze the mechanisms underlying cold adaptation in CA chickens and to explore their genetic characteristics. These results should provide a theoretical basis for improvement in chicken breeds.

MATERIALS AND METHODS

Sample collection and genome sequencing

We sampled a total of 10 CA chickens from a local breeder in Quebec, Canada. DNA was extracted from the blood of each individual using the standard phenol-chloroform protocol. DNA samples were sequenced with Illumina HiSeq technology at the McGill University and Genome Quebec Innovation Centre. Moreover, whole-genome resequencing data of 121 chickens from all over the world were collected, including Europe and America (Commercial, C (n=10)), Middle East (Iran, ME (n=30)), East Asia (China, EA (n=25)), Tibet, China (TB (n=10)), Yunnan, China, and Southeast Asia (Yunnan and Indonesia, YN.SEA (n=35)), South Asia (India, SA (n=10)), and Indonesia (green jungle fowl, GJF (n=1)) (Supplementary Table S1).

Sequence alignment and variant calling

Raw sequencing reads were trimmed using Trimmomatic v0.38 (Bolger et al., 2014) to filter low-quality bases and sequences using default parameters. All clean reads were then mapped to the chicken reference genome (GRCg6a) using the Burrows-Wheeler Aligner “BWA-MEM” algorithm with default parameters (Li & Durbin, 2009). On average, 99.17% of the reads sequenced in this study were mapped and the final average sequencing coverage was 20.37× (ranging from 15.36× to 25.26×) per individual. Using Picard Tools, reads were sorted and merged, and all duplicate reads were removed successively. Variant calling of sequence data was performed using the Genome Analysis Toolkit (GATK, v3.6-0-g89b7209) (McKenna et al., 2010). The applied criteria to all single nucleotide polymorphisms (SNPs) were: (1) <1/3× mean sequencing depth (over all included individuals)<3×; (2) mapping quality (MQ)>40.0; (3) Quality by Depth (QD)>2.0; (4) Fisher Strand (FS)>60.0; (5) MQRankSum>−12.5; (6) ReadPosRankSum>−8. Haplotype phasing was implemented using BEAGLE v4.1 ( Browning & Browning, 2007). Identified SNPs were further classified by ANNOVAR (Wang et al., 2010) based on the gene annotation of GRCg6a.

Population genetics analysis

An individual-based neighbor-joining (NJ) tree was constructed from the matrix of pairwise genetic distances of autosomal SNP data from 81 chickens (10 CA chickens, 10 samples from each group of C, ME, TB, and SA, 15 samples from each group of EA and YN.SEA, and one sample of GJF) using PLINK software (Purcell et al., 2007). The NJ tree was visualized by iTOL (Letunic & Bork, 2019). Using the smartpca program in EIGENSOFT v5.052 (Patterson et al., 2006), principal component analysis (PCA) was carried out based on pruned SNP data and eigenvector significance was detected by the Tracy-Widom test. Genetic structure clustering was deduced with the model-based assignment program ADMIXTURE v1.3.0 (Alexander et al., 2009) by estimating the ancestry of each individual using the genome-wide unlinked SNP dataset. The number of assumed genetic clusters (K) ranged from 2 to 7, with 200 bootstrap replicates. Using the qp3Pop program in the AdmixTools package (Patterson et al., 2012), we calculated outgroup f3 statistics, using GJF as the outgroup.

Run of homozygosity (ROH), genetic diversity, and linkage disequilibrium (LD) detection

The ROH for each group was identified using the --homozyg option implemented in PLINK. The following criteria were applied for ROH identification: (1) required minimum density (--homozyg-density) 50; (2) number of missing calls of a window (--homozyg-window-missing) 5; (3) number of heterozygotes of a window (--homozyg-window-het) 3. Length and average number of ROHs of each group were estimated, with length divided into four categories, i.e., 0.5–1 Mb, 1–2Mb, 2–4 Mb, and >4 Mb. Nucleotide diversity of each group was also calculated with a window size of 40 kb and step size of 20 kb using VCFtools ( Danecek et al., 2011). The decay of LD was calculated via the squared correlation coefficient (r2) using PopLDdecay (Zhang et al., 2019). Genomic inbreeding coefficients based on ROH (FROH) were estimated for each bird, according to: FROH=LROH/Ltotal (McQuillan et al., 2008).

Whole-mitochondrial genome phylogeny

The BAM files of the mitochondrial genomes of 10 CA chickens were transformed to FASTQ files and then assembled using Mapping Iterative Assembler v1.0 (Briggs et al., 2009). Using the maximum-likelihood (ML) model within IQ-TREE (Nguyen et al., 2015), we analyzed the mitochondrial phylogenetic relationships between CA chickens and 58 previously published wild fowls and domestic chickens (Miao et al., 2013) (Supplementary Table S2). TIM2+F+R2 was chosen as the best-fit model for the complete mitochondrial genome sequences using ModelFinder (Kalyaanamoorthy et al., 2017). Bootstrap support values of ML analysis were generated with 1 000 replicates (-bb 1000). The final tree topology was visualized using iTOL (Letunic & Bork, 2019).

Genome-wide selective sweep analysis

To identify genomic regions of selective sweeps associated with cold adaptation, we used several methods to investigate selection signatures, including population differentiation index (Fst) (Weir & Cockerham, 1984), nucleotide diversity (π), and cross-population extended haplotype homozygosity (XPEHH) (Szpiech & Hernandez, 2014). These three methods are often used in positive natural selection analyses. Statistics were calculated using a sliding window approach with a window size of 40 kb and step size of 20 kb. We calculated the average Fst, π-Ratio, and XPEHH values of SNPs in each window and used the outlier method to obtain windows with the top 1% ofFst, π-Ratio, and XPEHH values. Protein-coding genes in these outlier windows were annotated using ANNOVAR software (Wang et al., 2010) and were treated as candidate positively selected genes (PSGs). Using the KOBAS v3.0 tool (Xie et al., 2011), we performed Gene Ontology (GO) pathway analysis of the potential candidate genes to better understand their biological functions. Pathways with a false discovery rate (FDR)-corrected P-value<0.05 were considered significantly enriched.

RESULTS

Genome resequencing and variation

A total of 10 CA chickens were collected for genome resequencing (average depth of 20.37×). To facilitate comparisons with other chicken breeds, we combined our data with available whole-genome resequencing data of 121 individuals from seven geographical regions according to their locations, including C, ME, EA, TB, YN.SEA, SA, and GJF (phylogenetic outgroups) (Figure 1A; Supplementary Table S1). Functional annotation of the polymorphic sites identified a total of 4 681 038 SNPs in the 10 CA genomes, which included most intergenic regions (50.3%) and intronic regions (33.8%). Exons contained 2.0% of the total SNPs, with 23 911 nonsynonymous SNPs and 67 932 synonymous SNPs.

Figure 1

Population genetic analyses of CA chickens

A: Geographic variation in monthly mean surface temperature for location of all chickens used in the study. B: NJ tree of CA chickens sequenced in this study and publicly available whole-genome sequences of representative chickens worldwide. C: PCA with first (PC1) and second (PC2) principal components. D: Population structure analysis. E: Outgroup f3 statistics (GJF; CA, X), with higher f3 values suggesting more ancestral alleles shared by CA and X, and thus their closer relationship. CA, Chantecler chickens; C, commercial chickens; ME, domestic chickens from Middle East; EA, domestic chickens from East Asia; TB, Tibetan chickens; YN.SEA, domestic chickens from Yunnan, China, and Southeast Asia; SA, domestic chickens from South Asia; GJF, green jungle fowl from Indonesia.

Population genetic structure of CA chickens

To determine the population structure and genetic relationships of the different chicken groups, we conducted a series of analyses, including phylogenetic reconstructions, PCA, and Bayesian clustering, using whole-genome resequencing data. In total, 81 chickens were used to perform genetic relationships (see Supplementary Table S1 for detailed sample information). The NJ tree, rooted with GJF, showed that the CA group was closest to the C group, followed by the ME group (Figure 1B). Similarly, the genetic relationships between CA and the other chickens were also confirmed by PCA (Figure 1C). The PCA results showed that the CA group deviated from the other groups, suggesting that the CA group has unique characteristics. The admixture results provided additional corroborating evidence for the genetic relationships of CA and other chickens. In clustering analysis, when K=3, the CA group shared similar ancestries as chickens from the C and ME groups; when K=4, the CA group was separated independently (Figure 1D). We also used outgroup f3 (Figure 1E) to investigate the genetic relationships between CA and other chicken groups. Results from the outgroup f3 statistics unequivocally indicated that CA chickens have a higher genetic affinity to chickens from the C group, followed by the ME chickens. These results indicate that CA chickens are genetically closest to the C group, followed by the ME group.

To examine variations in mitochondrial DNA (mtDNA), we constructed a ML phylogenetic tree based on the complete mitochondrial genome sequences of 10 CA chickens and 58 downloaded chickens (Supplementary Table S2) that cover almost all known mitochondrial haplotypes. Results showed that CA chickens belonged to the mtDNA haplotype E1 (Supplementary Figure S1), which is the most widely distributed matrilineal lineage (Miao et al., 2013).

Patterns of genomic variations and LD

To explore the patterns of genomic variation, LD, and FROH in CA and other chickens, we analyzed ROH size, average number of ROHs, π, LD, and FROH. The SNPs (6 Mb) of each population were used for ROH analysis. The genetic variation parameters showed a basically consistent pattern in the different populations. Specifically, CA, C, and TB chickens exhibited larger ROHs and higher ROH numbers, whereas the opposite genomic variations were observed in the YN.SEA and SA chickens (Figure 2A, B). In addition, CA chickens not only exhibited lower π values but also a slow decay rate and visibly high level of LD (Figure 2C, D). The genomic diversity and LD patterns strongly suggest that CA chickens are likely affected by inbreeding or population founder effects. To further explore the reasons for the above observations, we calculated FROH. Of note, the higher FROH found for the CA chickens indicated a higher level of inbreeding in this group (Figure 2E).

Figure 2

Analysis of genomic variation, LD, and FROH for CA chickens

Length (A) and average number (B) of ROHs in CA chickens and other representative chickens of worldwide geographical distribution. C: Distributions of π in different groups. D: Decay of LD in different groups, with one line per breed. E: FROH in different chicken groups.

Genome-wide selective sweep analyses of CA chickens

To identify candidate population stratification genes potentially related to cold adaptation in the CA group, we performed genome-wide selective sweep analyses based on the genomic variation data in this study. Typically, those genomic regions related to traits of interest that have undergone selection pressure in one group should show strong differences in comparison to other groups, including high Fst, significantly reduced π, and long-range extended haplotype homozygosity (EHH) (Sabeti et al., 2006). Therefore, we used the abovementioned statistical methods to identify potential candidate genes that may play a significant role in the regulation of stress response to cold adaptation in CA chickens. We used 10 and 30 chickens from the CA and ME groups, respectively, for genome-wide selective sweep analyses. In the analyses, the top 1% was used as a threshold to identify potential candidate regions, and genes annotated in these regions were considered as potential candidate genes. In total, 310 genes with significantly higher FstCA-to-ME values (top 1%) were identified as candidate PSGs in the CA chickens (Figure 3A; Supplementary Table S3). As high population differentiation is affected by both positive selection and demographic history (Sabeti et al., 2006), the π-RatioME/CA and XPEHHCA-to-ME statistical methods were also adopted in this study. These two methods are less sensitive to population demographic events than other approaches (Innan & Kim, 2008). In total, 247 and 212 PSGs in the top 1% were identified using the π-Ratio and XPEHH methods, respectively (Figure 3B, C; Supplementary Tables S4, S5). Overall, 622 unique PSGs were found using the three methods, with only 36 shared by all three (Figure 3D; Supplementary Table S6). Unsurprisingly, the overlap among detected PSGs using different statistical methods was underwhelming. Different statistical methods based on different principles can yield contradictory results from the same dataset, which may lead to this rare overlap phenomenon (Sabeti et al., 2006). Importantly, all 36 PSGs were found by at least one of the above methods when CA chickens were regarded as the target population compared to domestic chickens from Indonesia (ID, Supplementary Table S6). These results indicate that all 36 PSGs are reliable PSGs in the CA chickens.

Figure 3

Genome-wide selection scan for PSGs in CA chickens using sliding window analysis (40 kb window size, 20 kb step size, kb increment, 99th percentile cutoff)

A: Selection signatures in CA chickens for FstCA-to-ME. B: Selection signatures in CA chickens for π-RatioME/CA. C: Selection signatures in CA chickens for XPEHHCA-to-ME. Threshold (top 1%) of FstCA-to-ME, π-RatioME/CA, and XPEHHCA-to-ME is marked with a horizontal black line. D: Number of PSGs identified in CA chickens by the three methods listed in each Venn diagram component.

To further elucidate the genetic mechanisms related to several candidate genes, functional enrichment analyses were applied to the three different PSG lists. GO analysis of the PSGs showed several significantly enriched categories related to various metabolic processes (corrected P<0.041), biological processes (corrected P<0.022), cellular processes (corrected P<0.0009), and biosynthetic processes (corrected P<0.047) (Supplementary Figures S2, S3 and Tables S7, S8), which may contribute to thermostatic maintenance in cold environments. Based on gene annotation, the identified PSGs were also related to different physiological functions and metabolic processes, such as neuronal development or migration, lipid metabolism, hair bundle morphogenesis, vasoconstriction and vasodilation, developmental processes, kidney development, skin pigmentation, and circadian rhythm. These processes are reported to play important roles in cold-induced thermoregulation (Adolph & Molnar, 1946; Barnett, 1959; Boulant & Dean, 1986; Lynch et al., 2015; Masoro, 1966; Smith, 1962; Yang et al., 2018).

ME3 and ZNF536 genes in CA chickens

To determine which of the 36 PSGs shared by the three methods in the CA group were most likely to be related to cold adaptation, we conducted haplotype and nonsynonymous mutation analyses and carefully searched the literature. Interestingly, we found a strong selective sweep spanning a 250 kb region (190 090 000–190 340 000 bp) on chromosome 1, which exhibited high FstCA-to-ME, π-RatioME/CA, and XPEHHCA-to-ME values (Figure 4A). This region contained five PSGs, including PRSS23, ME3,FAM181B,PRCP, andDDIAS, which may be related to cold adaptability in CA chickens. Further analyses revealed that the CA chickens had a core haplotype in this region, which was distinct from that of other chickens (Figure 4B). In particular, the genotype frequency of the nonsynonymous mutation in the ME3 coding region (c. A287G), which results in an amino acid substitution (p.H96R), was 90% in the CA chickens but was less than 36% in the other chickens (Figure 4C). In addition, the genotype frequency of a nonsynonymous mutation (p.A9S) in the DDIAS gene was 90% in the CA chickens, but less than 32% in the other chickens (Figure 4D). ME3 is an isoform of the malic enzyme, which is an oxidative decarboxylase that catalyzes the conversion of malate to pyruvate, simultaneously generating NADPH (Jiang et al., 2013; Pongratz et al., 2007). Malic enzymes participate in pyruvate cycling and insulin secretion (Pongratz et al., 2007). Compared to other tissues, malic enzymes are more abundant in lipogenic, white and brown adipose, and liver tissues, and may play a role in lipogenesis by providing NADPH (Chung et al., 1999; Wise & Ball, 1964). Our results indicated that the ME3 gene, which is related to fat production, likely plays an important role in cold adaptation in CA chickens. Although DDIAS functions as an anti-apoptotic protein in response to DNA repair and PRSS23 expression is regulated by estrogens, the other functions of these two genes remain to be elucidated (Huan et al., 2014; Im et al., 2016).PRCP andFAM181B may participate in cold adaptation by regulating angiogenesis and nervous system development (Hagedorn, 2011; Marks et al., 2016). In addition, ZNF536, a gene on chromosome 11, showed strong positive selection in the CA chickens, as exhibited by the high FstCA-to-ME, π-RatioME/CA, and XPEHHCA-to-ME values (Figure 5A), while shared haplotypes were rarely observed in the CA and other chickens (Figure 5B). A nonsynonymous mutation (c.G3685A) was also found in the ZNF536 coding region, generating an amino acid substitution (p.A1229T) in a conserved domain (Supplementary Figure S4). The derived genotype frequency in the CA chickens was 80% but was less than 8% in the other chickens (Figure 5C). ZNF536 is a highly conserved zinc finger protein that negatively regulates neuronal differentiation (Qin et al., 2009). It plays an important role in the development of forebrain neurons involved in social behavior and stress (Thyme et al., 2019). Thus, ZNF536 likely participates in cold adaptation in CA chickens.

Figure 4

Selection signals in 250 kb region of chromosome 1 in CA chickens

Figure 5

Selective signature on ZNF536 in CA chickens

A: Zoom of selective analyses in 250 kb region of chromosome 1. FstCA-to-ME, π-RatioME/CA, and XPEHHCA-to-ME were plotted as lines based on sliding window analysis (4 kb window and 2 kb increment). B: Degree of haplotype sharing between CA and other chickens in 250 kb region of chromosome 1, which contained five genes. Structure diagram of ME3 (C) and DDIAS (D) genes, with red bars representing exons, and allele frequency of mutant ME3 and DDIAS loci. Blue, gray, and orange represent homozygous mutant, heterozygous mutant, and homozygous wild-type, respectively.

A: Zoom of selective analyses in ZNF536 region. FstCA-to-ME, π-RatioME/CA, and XPEHHCA-to-ME were plotted as lines based on sliding window analysis (4 kb window and 2 kb increment). B: Structure diagram ofZNF536 gene, with red bars representing exons, and degree of haplotype sharing between CA and other chickens in ZNF536 gene. C: Allele frequency of mutant ZNF536 loci. Blue, gray, and orange represent homozygous mutant, heterozygous mutant, and homozygous wild-type, respectively.

Parallelism in cold-tolerant chickens

To explore whether the genetic mechanisms underlying cold tolerance in chickens in different regions are the same, we conducted genome-wide selective sweep analyses on representative cold-resistant chickens (i.e., Bian and Dagu chickens) in northern China (NCN). We used eight Bian chickens, eight Dagu chickens, and 30 ME chickens for analysis of selection signals. Based on selection signals analysis of the chickens, we detected 208, 1 023, and 228 PSGs in the top 1% using FstNCN-to-ME, π-RatioME/NCN, and XPEHHNCN-to-ME, respectively (Supplementary Figure S5 and Tables S9–S11). GO analyses of the three PSG lists showed that the significantly enriched categories not only included cell-related processes (corrected P<0.043), biosynthetic processes (corrected P<0.049), and metabolism-related processes (corrected P<0.049), but also regulation of neuronal differentiation, genesis, and development (corrected P<0.036) (Supplementary Figures S6, S7 and Tables S12, S13). Although the cold-tolerant NCN and CA chickens shared 68 PSGs, none of the 36 PSGs detected using the three methods in the CA chickens were found in the PSG list of NCN chickens (Supplementary Figure S8 and Table S6). These results suggest that parallel genetic mechanisms are likely responsible for cold adaptation in these two chicken breeds.

DISCUSSION

The primary breeding objective of animal husbandry has been to increase growth and production traits (e.g., meat, eggs, and milk) to meet the rapidly increasing demands of human consumption. Although climate has a profound impact on the survival of organisms, limited previous research has focused on animal adaptations to extreme climate conditions, such as cold adaptation. Recently, with increasing societal needs and technology development, as well as changes in the environment, animal adaptability has attracted increasing attention from both a breeding and biological perspective. Several studies have been conducted on cold adaptability in humans and animals (Fumagalli et al., 2015; Hallmark et al., 2019; Librado et al., 2015; Liu et al., 2014; Lynch et al., 2015). However, very few studies have focused on cold adaptability in poultry. In this study, 10 cold-tolerant CA chickens were used for whole-genome resequencing to understand the genetic mechanisms underlying cold adaptability in chickens. Our results represent a valuable resource and provide an excellent opportunity for exploring cold adaptability in chickens.

Although a pre-Columbian Polynesian introduction of chickens to the Americas has been suggested (Storey et al., 2007), the viewpoint that chickens were introduced to the Americas by European (Portuguese or Spanish) colonizers after 1500 AD has also been recognized (Herrera et al., 2020). The CA chicken was developed in the early 20th century for cold weather with high egg and meat production. Here, population genetic analysis indicated that CA chickens shared more ancestral components with C chickens (Figure 1) than the other chickens studied. Overall, 41 of the 622 PSGs identified in the CA chickens appear to be related to cold adaptation in humans and other animals. For example, TBX15, which plays a role in adipocytes (Gburcik et al., 2012), is reportedly related to cold adaptation in Inuits (Fumagalli et al., 2015), and was identified as a PSG in the CA chickens. TRPC4 is involved in the temperature sensitivity pathway (Hofmann et al., 2002) and has been identified in sled dogs (Sinding et al., 2020) and CA chickens. OCA2 andTYR, two key genes for skin lightening in humans (Shriver et al., 2003), were identified in CA chickens. We found limited PSG overlap between our study and previous reports, which may be due to parallel genetic mechanisms underlying cold adaptation in CA chickens and other animals. Our study and other results (Ai et al., 2015; Librado et al., 2015; Lv et al., 2014) show that cold adaptation is a complex process involving the nervous system, fat metabolism, hair formation, kidney development, vasoconstriction and vasodilation, pigmentation, and circadian rhythm. In selective sweep analysis, we identified 24 PSGs in the CA chickens involved in fat metabolism, including ME3, TBX15, LRP2, PPARG, PNPLA2, and CERS6, consistent with previous studies stating that fat metabolism plays an important role in cold adaptation in Greenlandic Inuits, Siberians, Yakutian horses, polar bears, and woolly mammoths (Fumagalli et al., 2015; Hallmark et al., 2019; Librado et al., 2015; Liu et al., 2014; Lynch et al., 2015). Of note, ME3, which was identified by several methods in our study (Figure 4), likely plays an important role in the cold adaptation of CA chickens. In addition, 58 of the identified PSGs in CA chickens, including ZNF536, CTNNA2, TIAM2, DISP3, KCNH7, and FBXO2, participate in neural processes, such as regulation of neuronal development, migration, differentiation, and neurogenesis, consistent with previous findings that the neural system plays an important role in environmental adaptation in Chinese domestic pigs (Ai et al., 2015). In the current study, ZNF536 was preliminarily shown to be related to cold adaptation in CA chickens (Figure 5). DDIAS, an anti-apoptotic protein involved in DNA repair (Im et al., 2016), may participate in cold adaptation in CA chickens, although other functions related to this gene need to be explored. The distinctive physical attributes of CA chickens, e.g., thicker and fluffier plumage, make them perfectly adapted to cold climates (Ekarius, 2007). Here, FZD3, which is involved in hair follicle development (Hung et al., 2001), was identified by selective sweep analysis in CA chickens. However, further research on this vital PSG is needed. We believe that ME3 andZNF536 are more reliable PSGs and are more likely to participate in cold adaptation, although further functional experiments are needed to understand their role in CA chicken cold adaptation.

In this study, we performed selective sweep analysis between two cold-tolerant chicken breeds from Canada (CA chickens) and China (Bian and Dagu chickens) with chickens from other parts of the world. Several loci on different chromosomes were identified as potentially involved in cold adaptation traits. Our results suggest that parallel genetic mechanisms may underlie cold adaptation in these chicken breeds, although additional physiological and functional experiments are needed for verification. Cold adaptation in chickens appears to have evolved in parallel through different pathways involving different complements of genes, similar to high-altitude adaptation in Tibetan chickens (Wang et al., 2015). Overall, our study revealed the population genetic structure and some PSGs in CA chickens that may be important in cold adaptation, contributing to our understanding of cold adaptation in poultry and providing a molecular basis for poultry breeding.

DATA AVAILABILITY

The raw FASTQ sequences were deposited in the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/) database under the BioProject accession No. PRJNA720223, and also deposited in the Genome Sequence Archive (https://ngdc.cncb.ac.cn/gsa) and Science Data Bank (https://www.scidb.cn/en) with accession No. CRA005002 and 10.11922/sciencedb.01142, respectively.

Supplementary data to this article can be found online.

Click here for additional data file.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRIBUTIONS

X.Z. and Y.J. conceived and supervised the study. N.Y.X. performed most analyses with contributions from M.L., M.G., and P.P.B. W.S. and J.L. prepared the samples. N.Y.X. wrote the manuscript with input from the other authors. X.Z., Y.J., and H.A.N. revised the manuscript. All authors read and approved the final version of the manuscript.