A total of 113 samples of Eurasian wigeon representing five populations and 92 samples of American wigeon from three flyways were included in analysis of mtDNA variability (Fig. 1; Table S1). In addition to published sequences (Peters et al. 2005; Peters et al. 2014; Kulikova et al. 2019), we sampled 15 Eurasian wigeons and 41 American wigeons. DNA was extracted using a DNeasy Blood and Tissue kit following the manufacturer’s protocols (Qiagen, Valencia, CA, USA). A total of 659–661 base pairs (bp) of mitochondrial DNA control region were sequenced (domains I and II) using primers L78 and H774 (Sorenson et al. 1999) following standard protocols (McCracken et al. 2001). All sequences have been deposited in GenBank (accession numbers PV787184- PV787239; see also Supplementary Material Table S1 for all samples numbers). Relationships among mtDNA haplotypes sequenced for American and Eurasian wigeons were reconstructed and visualized with a median-joining haplotype network in the program NETWORK v.10.2.0.0 (Bandelt et al. 1999). Pairwise nucleotide diversity (π), absolute divergence (d_XY), and the fixation index (Ф_ST) were calculated in R package POPGENOME v.2.7.5 (Pfeifer et al. 2014). ARLEQUIN v.3.5.2.2 (Excoffier and Lischer 2010) was used to calculate analysis of molecular variance (AMOVA) and pairwise mismatch distributions with Rogers’ (1995) model of sudden population expansion.

For the analysis of nuclear DNA, we obtained tissue samples from 70 Eurasian wigeons and 40 American wigeons collected across Asia from Western Siberia to Western Beringia and across North America (Fig. 1; Table S2). DNA was extracted as described above. Double-digest restriction associated DNA sequencing (ddRAD-seq) libraries were prepared following the protocol of DaCosta and Sorenson (2014). In brief, genomic DNA was digested using SbfI and EcoRI restriction enzymes, and Illumina TruSeq compatible barcodes were ligated for future demultiplexing. Ligated DNA fragments 300 to 450 bp in length were extracted from 2% low-melt agarose gels and purified using a MinElute gel extraction kit (Qiagen) and amplified using standard PCR. Magnetic AMPure XP beads (Beckman Coulter, Inc.) were applied to purify PCR products that were quantified using real-time PCR and an Illumina library quantification kit (KAPA Biosystems). Equimolar concentrations of each individual library were pooled and sequenced (150 bp reads) on an Illumina HiSeq 2500 at TUCF Genomics, Tufts University (Medford, MA, USA).

The raw Illumina reads were processed using the DaCosta and Sorenson’s (2014) computational pipeline (Python scripts available at https://github.com/BU-RAD-seq/ddRAD-seq-Pipeline; also see Lavretsky et al. 2015). For each individual, identical reads were combined into a single read while retaining read counts and the highest quality score for each position. Reads with an average Phred score of <20 were removed. The retained identical reads were concatenated and clustered into putative loci using USEARCH v.5 (Edgar 2010) with an identity threshold of 0.85. Loci were mapped into chromosomes by using BLASTn v.2 (Altschul et al. 1990) and mapped to a mallard (Anas platyrhynchos) reference genome (accession numbers SS263068950–SS263191362; Huang et al. 2013; Kraus et al. 2011).

The aligned sequences were then genotyped using python scripts available with the DaСosta and Sorenson (2014) pipeline. Homozygotes were defined when ≥93% of the reads were identical, whereas heterozygotes were scored if a second haplotype was represented by at least 29% of sequence reads or if as few as 20% of reads were consistent with a second allele and that haplotype was represented in other individuals. Genotypes were flagged if none of these criteria were met, or more than two haplotypes met the criteria, or if they were represented by ≤5 reads; for those genotypes, we retained only the allele represented by the majority of reads and scored the second allele as missing data. We retained all loci that contained ≤10% missing genotypes and ≤5% flagged genotypes.

A representative sequence from each of the final alignments was aligned to the reference mallard genome (assembly v.ZJU1.0, accession no. GCA_015476345.1). This was done to localize each locus and separate autosomal and Z-linked loci for downstream analyses. Only uniquely aligned reads were selected for further analysis.

Analysis of population structure was done using three methods. First, we used principal coordinates analysis (PCoA) based on the Euclidian distances between individual genotypes and implemented by dudi.pco in the R software package ADEGENET v.2.1.3 (Jombart 2008). The two-dimensional PCoA plots were drawn using the GGPLOT2 package version 3.5.1 (Wickham 2016). Then, maximum likelihood-based individual assignment probabilities were obtained in ADMIXTURE v.1.3.0 (Alexander et al. 2009). Data formatting for ADMIXTURE was done with PLINK v.1.07 (Purcell et al. 2007). We analyzed autosomal and Z chromosome linked biallelic SNPs separately with 10-fold cross-validation performed in each ADMIXTURE analysis and with a quasi-Newton algorithm employed to accelerate convergence (Zhou et al. 2011). We ran ADMIXTURE for K populations of 1–10 and applied a block relaxation algorithm for the point estimation that results in termination of analyses once the change in the log-likelihood of the point estimations increased by <0.0001. The optimum K was based on the lowest average of CV-errors across the analyses per K value. ADMIXTURE outputs were processed and visualized with CLUMPAK v.1.1. ADMIXTURE plots were produced with GGPLOT2 v.3.5.1. Next, we examined population subdivision by using the Bayesian clustering method implemented in STRUCTURE v.2.3.4 (Falush et al. 2003), which was run using an admixture model and correlated allele frequencies among populations without prior information regarding species or population designation. Twenty replicates for each value of K in the range of 1–10 were run with 200,000 steps of the MCMC after a 50,000-step burn-in for each run. We applied CLUMPAK v.1.1 (Kopelman et al. 2015) to process STRUCTURE output files and to determine the optimum K based on calculation of ΔK (Evanno et al. 2005). Although both STRUCTURE and ADMIXTURE use a similar underlying model to assign individuals to K clusters, they differ in their computational approach. STRUCTURE uses a Bayesian framework with Markov Chain Monte Carlo (MCMC) sampling, while ADMIXTURE uses a fast optimization algorithm based on maximum-likelihood estimation. By comparing the results from both methods, we ensured that the identified ancestry components were not an artifact of a specific algorithm’s assumptions or convergence behavior.

Finally, we utilized the R package POPGENOME v.2.7.7 (Pfeifer et al. 2014) to calculate composite species pairwise and population pairwise estimates of relative divergence (Φst), absolute divergence (d_XY), and nucleotide diversity (π) for concatenated FASTA files of autosomal and Z-chromosome linked ddRAD-seq loci. Φst was selected over traditional Fst because it incorporates the number of mutational differences between haplotypes, providing a more biologically realistic measure of genetic distance for our sequence-based ddRADseq data (Excoffier et al. 1992). We concatenated RAD loci to facilitate analysis in POPGENOME, which requires continuous sequence alignments. While this approach artificially links physically unlinked loci, it provides valid estimates for genome-wide summary statistics like π and d_XY, which are calculated as averages across sites. POPGENOME estimates the relative genetic distance between populations (Φst) using the statistic by Excoffier et al. (1992). For comparison, we also calculated the statistic of Hudson et al. (1992). Analysis of molecular variance (AMOVA) was run in ARLEQUIN v.3.5.2.2 to examine genetic differentiation within and among eight populations segregated from the two species. We also estimated the ratio of effective population sizes of Z chromosome and autosomes (Irwin 2018) calculating the ratio of adjusted Z diversity (Z diversity divided by an estimate of 1.1 for substitution rate ratio of Z vs. autosomes) to autosome diversity (π_Z/π_A).

We calculated the pairwise per locus values of Ф_ST, d_XY, and π in the r package POPGENOME as described above. Z-linked and some autosomal pairwise Ф_ST values were plotted by chromosomal position in EXCEL (i.e., Manhattan plots). To identify putative loci under selection we used 16,548 SNPs that is quite sufficient number for the task (Lotterhos and Whitlock 2015). We employed two complementary genome-scan approaches: a principal component-based method and a differentiation-based method. First, we used PCADAPT v.4.4.1 (Privé et al. 2020), which operates without a priori population assignments by identifying loci that are outliers in the multivariate genetic space defined by principal components. This method is particularly effective at detecting selection in the presence of complex population structure. We made two separate analyses, one with autosomal loci and another with Z-sex chromosome linked loci. Each analysis was performed using K=10 principal components, retaining loci with a false discovery rate < 5%. Second, we utilized the Bayesian approach implemented in BAYESCAN v.2.1 (Foll and Gaggiotti 2008), which explicitly models the interplay between population-specific effects and locus-specific effects of selection. This method calculates posterior probabilities for each locus by comparing a model including selection to a neutral model. BAYESCAN analyses included 20 pilot runs of 5000 steps each, followed by 100,000 burn-in and 200,000 sampling steps with a thinning interval of 10. The probability of false discovery rate (qval) was set at 0.01 and 0.05. The consensus set of outlier loci identified by both methods was considered to represent high-confidence candidates for being under divergent selection, thereby reducing the rate of false positives that can arise from the assumptions of any single method.

American wigeon had almost two-fold higher nucleotide diversity for mtDNA control region than Eurasian wigeon (Table 1); haplotype diversity showed a similar trend but to a lesser extent: 0.923 and 0.779, respectively. Mismatch distributions for Eurasian wigeon and American wigeon haplotypes were distributed normally and did not differ from Rogers’ (1995) model of sudden population expansion (Ps > 0.1). Five Eurasian wigeon populations had similar values of nucleotide as well as of haplotype diversity; the same was true for three American wigeon populations studied (Table S3). Genetic diversity values in populations of American wigeon were higher than in populations of Eurasian wigeon.

We recovered two expected deeply divergent species-specific mtDNA haplogroups (Peters et al. 2014) supported by high Ф_ST value of 0.88 (Fig. 2; Table 1). Two American wigeons: one from Alaska, USA, and the other from Saskatchewan, Canada, shared haplotypes with Eurasian wigeons in the Eurasian clade, and one Eurasian wigeon from Western Beringia (Anadyr) had a common haplotype of the American clade (Fig. 2). There was no within species grouping according to population designation as haplotypes were broadly shared among ducks from different populations within each haplogroup. AMOVA supported overall weak population structure in both species: 0.3% of observed genetic variation was partitioned among populations, 11.7% within populations, and 88% between species. As expected, based on the haplotype network (Fig. 2) and AMOVA results, we recovered very low relative differentiation among populations of American wigeon (Ф_ST = 0.0018–0.0030; Table S4c). However, in Eurasian wigeon the samples from Siberia and the North American Atlantic Flyway were differentiated from samples collected in Western Beringia, Russian Far East and the North American Pacific Flyway (Ф_ST = 0.14–0.39), but there was no differentiation among samples within these groups (Ф_ST = 0–0.0007; Table S4c). AMOVA confirmed these results and thus supported the presence of some barrier to mtDNA gene flow between Eurasian wigeons from Siberia and the Atlantic Flyway on the one hand and Western Beringia, Russian Far East, and Pacific Flyway on the other hand: 82% of the genetic variation was within populations, 0% was partitioned between populations within groups, and 18% was observed between groups.

After quality filtering, we recovered 3281 ddRAD-seq loci, with a mean depth of 199.1 reads per locus per individual. Among these loci, 3092 loci (388,590 aligned base pairs; 37,182 SNPs) were assigned to autosomes and 189 loci (23,513 aligned base pairs; 1567 SNPs) were assigned to the Z chromosome. The loci were evenly distributed across chromosomes (Table S5) with the number of loci per chromosome being proportional to the chromosome size (R² = 0.974).

Eurasian and American wigeon had similar autosomal and Z-chromosome nucleotide diversities (Table 1). For autosomes absolute divergence (d_XY) closely approximated the nucleotide diversities while for Z chromosome d_XY was slightly higher than nucleotide diversities. The estimated π_Z/π_A ratios were 0.442 and 0.443 for Eurasian and American wigeon, respectively. Genetic differentiation was higher for the Z chromosome (overall Ф_STZ = 0.192) in comparison to autosomes (overall Ф_STA = 0.0386) with the overall Ф_STZ/ Ф_STA ratio of 4.97. We note that another Ф_ST metric by Hudson et al. (1992), yielded nearly identical results (0.035 for autosomes and 0.196 for the Z chromosome), confirming the robustness of our finding of greater divergence on the Z chromosome. In contrast to Ф_ST, absolute divergence was higher for autosomes (d_XY = 0.0082) than for the Z chromosome (d_XY = 0.0048). On a locus-by-locus basis, 1.36 % and 0.48% of autosomal loci exhibited high (0.15 < Ф_ST < 0.25) and very high (Ф_ST > 0.25) divergence, respectively, whereas for Z chromosome the percent of high and very high divergent loci were 7.41% and 3.17%, respectively.

Plotting the first two principal coordinates from the PCoA clearly separated Eurasian and American wigeons with the first coordinate axis playing the main part in species separation (Fig. 3). However, PCoA failed to differentiate any groups inside these species. There was one Eurasian wigeon from Alaska, USA (North American (NA) Pacific) that occupied an intermediate position between Eurasian and American wigeons when using autosomal and Z-chromosome markers (Fig. 3a,b). PCoA based on autosomal markers identified one more Eurasian wigeon from California, USA (NA Pacific), that clustered between the main Eurasian cluster and the intermediate Eurasian wigeon from Alaska (Fig. 3a).

ADMIXTURE results were based on a total of 15,991 biallelic autosomal SNPs and a total of 553 biallelic Z-chromosome SNPs. We also made ADMIXTURE analysis with a single biallelic SNP randomly chosen from each locus with a total of 2664 biallelic autosomal SNPs and 153 Z-chromosome SNPs (Fig. 4). The results of full SNPs and single SNPs datasets analyses were similar. The optimal number of populations (K) was two for both autosomal and Z-chromosome loci (Fig. S1). At K = 2, ADMIXTURE results were consistent with PCoA, clearly distinguishing Eurasian from American wigeons. Again, analysis of autosomal SNPs revealed two Eurasian wigeons from the NA Pacific (California and Alaska) with mixed ancestry. One of them, a wigeon from Alaska, also received a mixed assignment in the Z-loci analysis (Fig. 4). Increasing K values up to 10 for both types of markers did not provide any additional resolution of population structure in the two species. STRUCTURE results were concordant with ADMIXTURE analyses. The two-population model (K = 2; Fig. S2) was best supported by delta K calculation for both autosomal and Z chromosome markers (Fig. S3). STRUCTURE also failed to resolve additional population structure at K > 2 (Fig. S2). Two Eurasian wigeons from Alaska and California showed evidence of admixture from American wigeon at autosomal loci with Q of 0.493 and 0.189, respectively, and the Alaskan individual had assignment to the sister species at Z loci with Q of 0.468. These two admixed Eurasian wigeons occupied intermediate position between Eurasian and American wigeon clusters in the PCoA plot (Fig. 3). Both putative hybrids shared mtDNA haplotypes with Eurasian wigeon.

A low level of population genetic differentiation across Z-chromosome loci was observed in both Eurasian and American wigeon. Thus, relative divergence ranged from –0.007 to 0.004 between Eurasian and from –0.001 to 0.005 between American populations, thus effectively zero in both species. In contrast, genetic differentiation between species was high and varied from 0.166 to 0.208 (Table S4a). Values of relative divergence across autosomal loci were higher than across Z-linked loci in both Eurasian and American wigeon and ranged from 0.004 to 0.021 (Table S4b). Lower levels of relative divergence were found in pairwise comparisons of NA Pacific and North Asian (Western Beringia and Far East) populations of Eurasian wigeons (Ф_ST = 0.004–0.006) and between NA Pacific and Central populations of American wigeons, whereas the most differentiated Eurasian wigeon population was the NA Atlantic (Ф_ST = 0.018–0.021). Values of autosomal genetic differentiation between the species were much lower than those of the Z chromosome (Ф_ST = 0.037–0.058 vs. Ф_ST = 0.166–0.208; Table S4a,b). Pairwise Ф_ST values based on mtDNA, Z-chromosome, and autosomal loci were strongly and significantly correlated (simple Mantel test) with each other (Fig. S4). AMOVA showed that 95.8% and 79.7% of the autosomal and Z chromosome genetic variability, respectively, were due to variability within the populations, while 3.7% and 20.2% of variability were due to interspecies differences (Table S6). Nucleotide diversity values were similar for all populations of Eurasian and American wigeon (0.0037–0.0039 for Z chromosome loci and 0.0077–0.0079 for autosomal loci; Table S3).

Comparing Eurasian and American wigeons and analyzing Z-linked and autosomal markers together as well as separately, PCadapt detected 14 Z-linked loci (7.4%) and 60 autosomal loci (1.9%) as outliers using a false discovery rate (FDR) of 0.01, and 17 Z-linked loci (9.0%) and 88 autosomal loci (2.7%) as outliers using FDR of 0.05. BAYESCAN identified six Z-linked loci (3.2%) and 58 autosomal loci (1.9%) as outliers using an FDR of 0.01 and nine Z-linked loci (4.8%) and 83 autosomal loci (2.7%) as outliers using FDR of 0.05 (Figs 5, S5). Notably, all Z-linked loci and the vast majority of autosomal loci identified by BAYESCAN were a subset of those detected by PCadapt, indicating a strong consensus. Consequently, for subsequent analyses, we prioritized the outlier loci identified by BAYESCAN due to its more conservative and model-based Bayesian framework, which provides a higher level of confidence that the detected signals are true signatures of selection rather than artifacts of population structure. All autosomal and Z-chromosome outliers identified by BAYESCAN were estimated to be under diversifying selection. Z-linked outliers had significantly higher estimates of d_XY than non-outliers (0.0097 vs. 0.0044, two-tailed t-test P = 0.004), whereas d_XY at autosomal outliers was almost the same as non-outliers (0.0078 vs. 0.0079, two-tailed t-test P = 0.93). Among nine Z-outliers, three demonstrated fixed differences and diagnostic SNPs, although one putative hybrid had both SNP variants in its genome (Fig. 6). The other six Z outliers showed significant allele frequency differences with Ф_ST values varying from 0.86 to 0.94; only one to three individuals, mainly American wigeons, had the SNP variant characteristic of the sister species, and the putative hybrid was heterozygous. Autosomal outliers were found to be located on chromosomes 1–12, 14, 15, 19, 22, 24, and 27 and had Ф_ST from 0.11 to 0.68 (Fig. S5). There were no fixed allele differences and species diagnostic SNPs at autosomal outlier loci. Alignment of Z chromosome outliers and autosomal outliers with Ф_ST > 0.5 to the mallard genome assembly ZJU1.0 (GCA_015476345.1) revealed that outliers mostly resided in introns of different protein coding genes (Tables S7–S8).

Diversity estimates across mitochondrial DNA haplotypes and nuclear ddRADseq-loci for American and Eurasian wigeons (Tables 1, S3) were similar to those of other duck species. Thus, autosomal and Z chromosome nucleotide diversity varied from 0.00570 to 0.00678 and from 0.00238 to 0.0038, respectively; for mitochondrial DNA nucleotide diversity was in the range of 0.0020 to 0.0120 in other species of dabbling ducks (McCracken et al. 2001; Kulikova et al. 2005; Peters et al. 2007, 2016; Lavretsky et al. 2015, 2019). Population genetic structure within species was not prominent, and there was very low genetic differentiation between populations based on autosomal and Z-chromosome markers (Table S4). The same results were obtained with other methods: PCoA, ADMIXTURE and STRUCTURE did not resolve population structure within these species (Figs 3, 4, S2) American wigeon populations were also undifferentiated based on mtDNA analysis. The overall lack of genetic structure is therefore widespread, which is a well-known phenomenon in migratory dabbling ducks (Kulikova et al. 2005; Flint et al. 2009; Kraus et al. 2011; Peters et al. 2014; Kulikova et al. 2019). Wigeons, as with mallards for example, exhibit considerable population connectivity and relatively high gene flow across Eurasia and North America. Natal dispersal, common wintering and breeding grounds, and distant annual migrations covering thousands of kilometers including transhemispheric movements contribute to the redistribution of ducks among different geographic regions (Ostapenko et al. 1997), exchange of migrants, and thus high gene flow and weak population structure.

However, mtDNA variability revealed subtle population structure in Eurasian wigeon. Samples from Siberia and the NA Atlantic Flyway were differentiated from the samples collected in Western Beringia, Russian Far East, and NA Pacific Flyway (Ф_ST = 0.14–0.39; Table S4c), which altogether with AMOVA results supported the presence of some barrier to mtDNA gene flow between Siberia and North Asia. We believe it could be explained by common wintering grounds for these groups of populations. It is well known that formation of pairs in waterfowl usually takes place during wintering (Rohwer and Anderson 1988). Wintering areas of European and Siberian wigeon populations overlap from the north of Spain through France, Germany and eastward to the Ural Mountains, while North Asian populations, i.e., Russian Far East and Western Beringia, usually winter in the south of Korean Peninsula and Japan (Bhushan et al. 1993; Ostapenko et al. 1997). Obviously, Eurasian wigeons from the NA Atlantic and Pacific Flyways are migrants from Europe and North Asia, respectively.

Both mtDNA and nuclear DNA confirm that American and Eurasian wigeons are genetically differentiated from each other. Each of two deeply divergent mtDNA haplogroups (Fig. 2) was predominantly found in one of the species as described earlier (Peters et al. 2014). Limited intermixing of mtDNA haplotypes (three of 207, 1.45%, wigeons carried a haplotype of the other species) and deep divergence of mtDNA haplogroups support reciprocal monophyly of mtDNA with occasional intermixing of haplotypes between species due to long-distance dispersal and hybridization. Multiple ddRAD-seq genotypes were different between species. Wigeons were assigned unambiguously to their species-specific clusters on the basis of autosomal and Z-chromosome ddRAD-seq genotypes equally well. We found limited evidence of admixture, suggesting hybridization between these species: two Eurasian wigeons from Alaska and California showed evidence of admixture from American wigeon at autosomal loci, and a wigeon from Alaska was also admixed at Z-linked loci (Figs 3, 4, S2). This Alaska wigeon was probably a F1 hybrid with assignment probabilities to the parental species of around 50% at both autosomal and Z-chromosome markers. A wigeon from California was presumably a F2 hybrid with 74–81% and 26–19% assignment to parental species at autosomal loci (Figs 4, S2). Both hybrids shared mtDNA haplotypes with Eurasian wigeon. One Eurasian wigeon from Western Beringia and two American wigeons from Alaska, USA, and Saskatchewan, Canada, carrying mtDNA haplotypes of the other species (Fig. 2) grouped with conspecifics on the basis of ddRAD-seq genotypes. Hence, four out of five hybrids of different generational crosses had mtDNA haplotypes of Eurasian wigeon. This is in concordance with the shortage-of-mates hypothesis (Rohwer et al. 2022) which assumes: [1] a deficit of females or strongly male-biased sex ratio in wigeon flocks (Bellrose et al. 1961); [2] infrequent wintering of Eurasian wigeon in flocks of American wigeon, on the west coast of North America in particular (Edgell 1984); [3] pair formation in wigeons on wintering grounds; and [4] pairing of males American wigeon with females Eurasian wigeon being much more likely than vice versa because this should result in higher number of hybrids with Eurasian wigeon mtDNA haplotypes. However, Rohwer et al. (2022) obtained the opposite results with 10 of 11 of F1 generation hybrids were sired by Eurasian wigeon, which they interpreted as a possible preference by female American wigeons for male Eurasian wigeon due to the more aggressive behavior of the latter and therefore more effective protection of females from other males.

The previous study of divergence between American and Eurasian wigeon based on sequencing mtDNA control region and 20 nuclear introns (Peters et al. 2014) revealed prominent mito-nuclear discordance. Divergence of nuclear DNA (Ф_ST = 0.046) was much lower than that of mtDNA (Ф_ST = 0.812), and the one Z-sex chromosome locus, chromo-helicase-DNA binding protein 1, CHD1Z, had an intermediate value of (Ф_ST = 0.415). Our estimates reconfirmed values obtained previously by Peters at al. (2014), and thus confirm mito-nuclear discordance within these species: nuclear autosomal DNA Ф_ST = 0.039 vs mtDNA Ф_ST = 0.879 and Z-loci Ф_ST = 0.192. Among the multiple causes behind the discordance are sex-biased dispersal, historical and modern gene flow, mitochondrial capture, natural selection and incomplete lineage sorting (Toews and Brelsford 2012). Two or more processes can act concurrently so that each case of mito-nuclear discordance is most probably the result of unique evolutionary history. Mito-nuclear discordance between American and Eurasian wigeon is of the most widespread type for speciating taxa when there is strong differentiation in mtDNA, moderate differentiation in sex chromosomes, and quite low differentiation in autosomes (reviewed in Irwin 2018). Similar patterns of genomic variability were observed in several pairs of recently diverged dabbling duck species (Peters et al. 2014; Lavretsky et al. 2015; Kulikova et al. 2022). High mutation rate, rapid sorting of haploid mtDNA, which has ¼ the effective population size of autosomes, and lack of recombination contribute to the rapid fixation of mitochondrial mutations via genetic drift. These characteristics also result in positive selection being much faster and more effective in producing divergent species-specific mtDNA haplogroups (Hill 2020). Another potential cause of higher genetic differentiation at mtDNA is that introgression might be hindered as a result of lower viability of hybrid females in comparison to hybrid males (i.e., Haldane’s rule; Haldane 1922). Although it is well established that hybrid fertility tends to decrease with increasing genetic distance—implying that Haldane’s Rule typically does not apply to closely related species (Ottenburghs et al. 2016)—there is evidence that Haldane’s Rule has played a role in another closely related pair of dabbling ducks: the Mallard and the American Black Duck. A postmating isolating mechanism between the two species resulted in a lower proportion of F1 females in hybrid clutches compared to non-hybrid crosses (Kirby et al. 2004). Lower viability of the heterogametic sex can restrict or even prevent gene flow in mtDNA and could also contribute to reduced gene flow in Z-chromosome linked loci (Carling and Brumfield 2008). Therefore, the pattern of mito-nuclear discordance we observe in “northern” wigeons is in concordance with Haldane’s rule in addition to male-mediated gene flow, with stronger drift and faster selection acting on the Z chromosome.

Eurasian and American wigeon had five times higher differentiated Z chromosomes than autosomes (Ф_STZ/Ф_STA = 4.97), whereas genetic diversity of Z-linked loci was almost one half of genetic diversity of autosomal loci (Table 1, S3). Thus, the ratios of Z-linked relative diversity (π_Z/π_A) were 0.442 and 0.443 for Eurasian and American wigeons, respectively, which are far below the expected ratio of 0.75 under equilibrium conditions, and below the mean π_Z/π_A of 0.61 observed for 33 bird species reviewed in Irwin (2018). This ratio is also below the theoretical minimum of 0.56 expected in cases of extreme variance in male reproductive success with a reduction in the effective number of mating males compared to mating females (Charlesworth 2001). Although we do not have much empirical data on male reproductive success in ducks, the evidence we do have suggests that the variance is probably modest (Peters et al. 2003; Kreisinger et al. 2010). Ducks are seasonally monogamous and despite forced extrapair copulations (FEPC) by males being common in waterfowl they barely result in skewed reproduction: only 4.2% of gadwall offspring had genotypes consistent with extrapair fertilizations (Peters et al. 2003), and 9.3% of all offspring were due to extra-pair paternity (EPP) in mallard (Kreisinger et al. 2010). Such frequency of successful FEPCs is not enough to shift the π_Z/π_A ratio even to its theoretical minimum under EPP. Thus, selection might be a better explanation of the observed difference between nucleotide diversities observed on the Z chromosome and autosomes. Specifically, both background selection and selective sweeps can reduce the effective population size by eliminating neutral genetic variation that is physically linked to regions under selection. Given that recombination is generally lower for the Z-chromosome, the effects of linked selection on its genetic diversity are expected to be more pronounced. Low d_XY and high Ф_ST for the Z chromosome compared to autosomes (Table 1) support this inference. Selection in a common ancestor could also have reduced the variation in Z-chromosome with subsequent selection in daughter taxa reducing the nucleotide diversity even more (Cruickshank and Hahn 2014). Thus, the Z chromosome has potentially undergone recurrent selection, where a selective sweep reduced variation within populations both before and after a population split. Such a pattern has been observed in many recently diverged species of birds, and it is more pronounced in species with larger population sizes in which selection would be more effective (Irwin 2018; McCallum et al. 2024; Schield et al. 2025). Indeed, population sizes of American and Eurasian wigeon number in the millions, and therefore, even weak selection can overpower the effects of genetic drift.

The percent of outlier loci under selection was slightly higher among Z-linked loci (4.8%) than among autosomal loci (2.7%). The difference between Z and autosomal loci was much less pronounced in wigeons than in species of the mallard group (Lavretsky et al. 2015; Kulikova et al. 2022). A moderately higher number of outliers on the Z-chromosome could be due to Z-chromosome demography and inheritance that can produce higher neutral divergence even in the complete absence of selection. Additional data are needed to better test the roles of neutral demography and selection on these species. Whole-genome resequencing (WGS) could provide greater precision and clarity for making these inferences.

Nevertheless, Z-linked outliers had elevated absolute divergence (d_XY), depressed within-species nucleotide diversity (π), and high Ф_ST values (mean Ф_ST = 0.93) in comparison to Z non-outlier loci, whereas the absolute divergence of autosomal outlier and non-outlier loci was nearly identical. Furthermore, within-species nucleotide diversity of outlier loci was slightly depressed only for the Eurasian wigeon, and Ф_ST values (mean Ф_ST for autosomal outliers = 0.37) were much lower than for Z-outliers. The signal of strong divergent selection across Z chromosome is also supported by the findings of three outlier loci with species diagnostic SNPs and six outliers with SNPs close to be species-specific with just one to three individuals having the SNP variant characteristic of the sister species (Fig. 6). At the same time there were no fixed SNP variants in autosomal outlier loci. These differences between Z-linked and autosomal outliers further support strong directional selection influencing Z-linked outliers. However, ddRAD loci constitute just 0.0003–0.0004% of the whole genome and, therefore, there could be autosomal loci under strong positive selection that were undetected by ddRAD-seq.

It is noteworthy to mention that alignment of nine Z-chromosome outliers to the mallard genome reveled four of them residing within introns of protein-coding genes and two Z outliers located in introns of long non-coding RNA genes (Table S7). In particular, two of the species diagnostic SNP variants are located in the intron of Janus Kinase 2 (JAK2), a non-receptor tyrosine kinase gene, which is an essential component of signal transduction of the class II cytokine receptors associated with growth and reproduction traits in chickens (Liu et al. 2010). As for long non-coding RNA (lncRNA), they are known to regulate various aspects of cell differentiation and development and affect many other physiological processes in plants and animals (reviewed in Mattick et al. 2023). Among 14 of the most significant autosomal outliers, three are located in lncRNA, and nine outlier loci reside in introns of various signaling proteins, like proteins of ion channels, receptors etc. (Table S8). Despite differences in plumage color between American and Eurasian wigeon, especially nuptial males’ plumage, and association of such traits with sex chromosomes in a number of studies (Toews et al. 2016; Campagna et al. 2017; Wang et al. 2020; Nikelski et al. 2024), no Z-chromosomal or autosomal outlier loci were aligned within proximity of any pigmentation genes (Tables S7, S8). However, the small part of the genome we could analyze with ddRAD-seq could easily have contributed to us missing those associations. Furthermore, plumage traits are typically under control of many gene networks including various regulatory ones (Davoodi et al. 2022; reviewed in Price-Waldman and Stoddard 2021; Kulikova 2021), which makes associations between genes and plumage more challenging to detect.