In this study, we downloaded publicly available whole-genome resequencing data of 484 Chinese duck individuals from the National Center for Biotechnology Information Sequence Read Archive database prior to November 26, 2023. We adopted a stringent and uniform pipeline to analyze raw genome data and filter single nucleotide polymorphism (SNP)s. Briefly, we first trimmed adaptors and removed low-quality reads via Trimmomatic (v0.39) with the following parameters: LEADING: 3 TRAILING: 3 SLIDINGWINDOW: 4:15 MINLEN: 36. We then mapped the clean reads to the Pekin duck (Anas platyrhynchos) reference genome (ZJU1.0) via the Burrows-Wheeler Aligner (v0.7.17) [10] with the parameters mem -M–R. We used SAMtools (1.13) [11] to sort bam files and the MarkDuplicates module in GATK (4.2.3.0) [12] to remove duplicates. We used mosdepth (0.3.3) [13] to calculate the sequencing depth of each sample and selected samples with sequencing depths greater than 3×. Finally, we called variants through the modules HaplotypeCaller, CombineGVCF, and GenotypeGVCF in GATK (4.2.3.0) and filtered SNPs via the module SelectVariants with the following parameters: QD<2.0 || MQ<40.0 || FS>60.0 || SOR>3.0 || MQRankSum<−12.5 || ReadPosRankSum<−8.0. We utilized BCFtools (v1.9) [11] with the parameters --min-af 0.01:minor -e ‘F_MISSING>0.9’ -m 2 -M 2 -v snps to obtain biallelic SNPs, and these SNPs were then more rigorously filtered via PLINK (v1.9) [14] with the parameters minor allele frequency<0.05 and Hardy-Weinberg equilibrium <0.000001.

After removing 68 samples with sequencing depths below 3×, we filtered 416 high-quality individual ducks with an average depth of 7× coverage for the following analysis. These samples were derived from Bioproject IDs: PRJNA599025, PRJNA658213, PRJNA888455, PRJNA409830, PRJNA419832, PRJNA450892, and PRJNA549423. We categorized these 416 ducks into wild populations (72), local populations (159), and commercial populations (185) (Table 1). The wild breed included 72 Mallard ducks, whereas the commercial populations included Pekin ducks, Mapleleaf ducks and Cherry Valley ducks. The local populations included 18 breeds, i.e., JinDing duck, Liancheng White duck, MaWang duck, Putian Black duck, SanSui duck, ShanMa duck, ShaoXing duck, Taiwan duck, Yulin-Wu duck, DongLan duck, YouXian duck, GaoYou duck, Ji’an Red duck, Longsheng-Cui duck, WenQiao duck, Rongshui-Xiang duck, Xilin-Ma duck, and Yulin-Ma duck.

After quality control, we obtained 5,140,908 high-quality SNPs, and these SNPs were annotated via the ANNOVAR package [15] on the basis of the gene annotations of the Pekin duck reference genome (ZJU1.0). We divided these SNPs into five classes: exon regions, intron regions, splicing sites, upstream and downstream regions and intergenic regions.

We performed principal component analysis (PCA) for 416 ducks by using PLINK (v1.9), with 5,140,908 SNPs shared among the three populations. We constructed an unrooted neighbor-joining (NJ) phylogenetic tree via iTOL [16] on the basis of the identity-by-state genetic distance matrix among 416 individuals calculated via PLINK (v1.9). Population genetic structure was inferred via ADMIXTURE software (v-1.3.0) [17], which employs a maximum likelihood-based method to estimate individual ancestries from multilocus SNP genotype datasets. We then performed visualization by using TBtools-II [18] while setting a predefined range of genetic clusters (K) from 2 to 12 to encompass the maximum number of lineages.

We estimated the population genetic diversity of the wild, local and commercial populations. The genetic diversity included the effective population size (Ne), inbreeding coefficient (F), expected heterozygosity (He), observed heterozygosity (Ho), nucleotide diversity (Pi), and regions of homozygosity (ROH). The effective population size (Ne) of different populations was estimated via SNeP (v1.1) [19] with independent SNPs calculated via PLINK (v1.9) with the parameters --indep-pairwise 50 5 0.2, and 72 samples from each population were used to avoid the effect of the number of samples. We selected 72 individuals from commercial and local populations on the basis of the corresponding proportion of each breed included in each population. We calculated F (--het), He (--hardy), Ho (--hardy) and ROH (--homozyg --homozyg-window-snp 50 --homozyg-window-het 1 --homozyg-snp 50 --homozyg-window-missing 5 --homozyg-window-threshold 0.05 --homozyg-density 500 --homozyg-gap 1000) for each duck of the wild population, local population, and commercial population, respectively. We calculated F, He, Ho, and ROH for each duck from the wild population, local population, and commercial population, respectively. We used VCFtools (v1.17) [12] to compute Pi values by scanning all genomic regions within 50-kb sliding windows with a step size of 25 kb. We compared the patterns of linkage disequilibrium (LD) among wild ducks, local ducks, and commercial ducks through pairwise LD estimates, which were measured as r² values calculated from PopLDdecay (v3.42) [20]. The samples used for the LD estimates were the same as those used for the Ne calculation.

First, we calculated genetic differentiation coefficient (Fst) values for the pairwise populations among the wild population, local population and commercial population using VCFtools (v1.17) with 50-kb sliding windows and a step size of 25 kb. To identify selective genomic signals in commercial or local populations compared with those in the other two populations, we then derived the locus-specific branch length (LSBL) [21] from the pairwise Fst distances calculated as follows: LSBL_CM–(LC+WT) = (Fst_CM–LC + Fst_CM–WT − Fst_LC–WT)/2 and LSBL_LC–(CM+WT) = (Fst_CM–LC + Fst_LC–WT − Fst_CM–WT)/2, where Fst_CM–LC refers to the Fst value between the commercial population and the local population, Fst_CM–WT refers to the Fst value between the commercial population and the wild population, and Fst_LC–WT refers to the Fst value between the local population and the wild population. We selected the top 1% of LSBL values of the genomic windows as candidate regions. Moreover, for the local population, we also analyzed selective signatures between egg and dual-purpose breeds by calculating Fst and Pi values from VCFtools (v1.17) with 50-kb sliding windows and a step size of 25 kb. We annotated the genes related to these important genomic regions, and we utilized the clusterProfiler R package [22] for the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of these candidate genes. The GO and KEGG pathways were selected from the top 10 enriched pathways in terms of Benjamini-adjusted p values. For some importantly selected genomic regions, we employed LSBL and Pi values with a 10-kb sliding window and a 5-kb step size across wild, local, and commercial populations to conduct comparisons and to identify candidate genes within these selected genomic regions.

We obtained 416 ducks with a wide and representative geographic distribution across various regions in China, yielding 5,140,908 high-quality SNPs distributed across 29 chromosomes (Figures 1A, 1B). As shown in Figure 1C, among these SNPs, we identified 4,321,515 SNPs in the wild population, 4,491,942 in the local population, and 3,746,643 in the commercial population. There were 2,832,649 SNPs shared among these three populations, accounting for approximately 55.1% of the total SNPs. The commercial population had 195,701 unique SNPs, the local population had 222,579 unique SNPs, and the wild population had 135,613 unique SNPs.

In the functional annotation of SNPs, the percentages of SNPs located in introns, intergenic regions, introns and exons of noncoding RNAs, untranslated regions3 regions, downstream, and exonic regions were very similar among the wild, local, and commercial populations. For example, the majority of SNPs are within introns, accounting for approximately 54.2%, 54.06%, and 54% of the total number of SNPs in the wild, local, and commercial populations, respectively. In contrast, the proportions of SNPs located in exons were relatively low, i.e., approximately 1.4% for these three populations (Figure 1D).

PCA clearly revealed that individuals from the wild, local, and commercial populations were clustered, with PC1 and PC2 accounting for 30.61% and 11.57% of the total variation, respectively (Figure 2A). The NJ tree also confirmed that the ducks were divided into three groups: wild, local, and commercial populations (Figure 2B). Moreover, ADMIXTURE analysis indicated that the wild, local, and commercial populations were clearly distinguishable at K = 3, with a log-likelihood value of −1,637,271,699.91 and a cross-validation error of 0.49702 (Figure 2C).

As shown in Figure 3A, in terms of Ho, He, and Pi, commercial duck breeds presented the lowest values (He = 0.22±0.15, Ho = 0.19±0.13, and Pi = 0.0011±0.0009), suggesting a relatively high degree of inbreeding. This was supported by the inbreeding coefficient, where commercial breeds also had the highest value (F = 0.26±0.15). The commercial duck breeds also presented the longest average ROH length (1905.33±1371.00 kb) among the three populations. The commercial population had the smallest Ne at generation 999, except that the Ne of the wild population was lower than that of the commercial population before generation 150, while the local population consistently had the largest Ne across all generations (Figure 3B). Compared with the other two populations, the commercial population presented a greater extent of LD (r²) (Figure 3C). These results were consistent with our general understanding of the wild, local, and commercial populations.

As shown in Figure 4A, we identified 152 genes by using the top 1% LSBL_CM–(LC+WT) genome regions of the commercial duck breeds. We subsequently conducted KEGG and GO enrichment analyses for these genes (Figures 4D, 4E). For the KEGG pathways, the mTOR signaling pathway and ErbB signaling pathway were significantly enriched (p<0.05). The GO enrichment analysis revealed that the enriched pathways were highly related to muscle development, such as cartilage development, connective tissue development, myofibrils and contractile fibers. Furthermore, we analyzed two important genomic regions by calculating LSBL and Pi with a 10-kb sliding window and 5-kb step and identified the candidate genes LRP6, BORCS5, MANSCI, BCL2L14, WNT7B, and PRR5 in the region Chr1: 70.25 to 74.00 Mb (Figure 4B) and BRD9 and EDN1 in the region Chr2: 97.10 to 99.76 Mb (Figure 4C).

Similarly, we detected 277 genes involved by using the top 1% of LSBL_LC–(CM+WT) values (Figure 5A). Through KEGG enrichment analysis, other types of O-glycan biosynthesis and mucin type O-glycan biosynthesis were significantly enriched (p<0.05), which may be linked to the immune system and disease resistance of the local population (Figure 5E). The significantly enriched GO pathways were regulation of axonogenesis, regulation of neuron projection development, protein glycosylation, and macromolecule glycosylation (Figure 5D). These pathways have implications for neurodevelopment and protein glycosylation processes. By using LSBL and Pi values with a 10-kb sliding window and 5-kb step for two important genomic regions, we identified the candidate gene NCAM2 in the region Chr1: 103.60 to 107.58 Mb (Figure 5B) and MPHOSPH6 in the region Chr12: 3.53 to 6.95 Mb (Figure 5C).

We identified 415 genes in the overlap of the top 5% Fst and Pi ratios between egg-laying and dual-purpose ducks in the local population. For these genes, KEGG pathway analysis revealed that the ovarian steroidogenesis pathway and protein digestion and absorption pathway were significantly enriched (p<0.05) (Figure 6C). GO analysis revealed that the regulation of monoatomic anion transport were the most significantly enriched terms (Figure 6B). We identified two candidate genes, i.e., PTGS2 and PLA2G4A, which are involved in the ovarian steroidogenesis pathway, that are positively selected in the egg laying population, whereas KCNK16, KCNK5 and KCNK17 are positively selected in dual-purpose ducks and are involved in the fundamental processes of protein digestion and absorption and monoatomic anion transport.