Genetic differentiation and population structure of "northern" wigeons (Anseriformes: Anatidae: Mareca americana, M. penelope)

757 Genetic differentiation and population structure of “northern” wigeons (Anseriformes: Anatidae: Mareca americana, M. penelope) Irina V. Kulikova1, Philip Lavretsky2, Kevin G. McCracken3,4,5, Yury N. Zhuravlev1, Irina L. Miroshnichenko1, Andrew B. Correll6, Jeffrey L. Peters6 1 Federal Scientific Center of the East Asia Terrestrial Biodiversity, Far Eastern Branch of the Russian Academy of Sciences, 159 100-let Vladivostoku avenue, Vladivostok 690022, Russia 2 Department of Biological Sciences, University of Texas at El Paso, 500 West University Avenue, El Paso, Texas 79968, USA 3 Department of Biology, University of Miami, 1301 Memorial Drive, Coral Gables, Florida 33124, USA 4 Department of Marine Biology and Ecology at the Rosenstiel School of Marine, Atmospheric, and Earth Science, University of Miami, 4600 Rickenbacker Causeway, Miami, Florida 33149, USA 5 Human Genetics and Genomics at the Miller School of Medicine, 1501 NW 10th Avenue, Biomedical Research Building, Miami, Florida 33136, USA 6 Department of Biological Sciences, Wright State University, 3640 Colonel Glenn Hwy, Dayton, Ohio 45435, USA Corresponding author: Irina V. Kulikova ([email protected]) Academic editor Martin Päckert | Received 6 August 2025 | Accepted 21 November 2025 | Published 3 December 2025 Citation: Kulikova IV, Lavretsky P, McCracken KG, Zhuravlev YuN, Miroshnichenko IL, Correll AB, Peters JL (2025) Genetic differentiation and population structure of “northern” wigeons (Anseriformes: Anatidae: Mareca americana, M. penelope). Vertebrate Zoology 75: 757–772. https://doi.org/ 10.3897/vz.75.e167908 Abstract Eurasian wigeon (Mareca penelope) and American wigeon (Mareca americana) are sister species with diagnosable differences mostly in male plumage. They breed in the Palearctic and Nearctic, respectively, but due to transoceanic migrations come in contact in North America, Western Europe, and north-eastern Asia, where they occasionally hybridize. To estimate genomic divergence and study their population structure we sequenced mitochondrial (mt) DNA control region and obtained 3092 autosomal and 189 Z chromosome loci from double-digest restriction associated DNA sequencing (ddRAD-seq). Consistent with previous work with few nuclear loci, we observed discordant patterns between mtDNA and nuclear DNA divergence. Deeply divergent species-specific mtDNA haplogroups contrasted with low autosomal differentiation and moderate Z-sex chromosome divergence. Meanwhile, Z-linked differentiation (ФST = 0.192) between taxa was five times higher than differentiation of autosomal loci (ФST = 0.0386), with four fixed and eight nearly fixed differences in SNPs discovered in three and six Z-linked outlier loci, respectively. No species-specific SNP variants were found among 83 autosomal outlier loci. This elevated Z-chromosome differentiation is most likely the result of selection that has been important in speciation. The lack of population genetic structure within Eurasian wigeon and American wigeon supports the common notion that migratory waterfowl have high dispersal ability that contributes to strong genetic connectivity between geographic populations. Keywords Diagnostic single nucleotide polymorphisms, elevated Z-sex chromosome divergence, mito-nuclear discordance, population genomics, population structure, wigeon Vertebrate Zoology 75, 2025, 757–772 | DOI 10.3897/vz.75.e167908 Copyright Irina V. Kulikova et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Kulikova IV et al.: Population genomics of “northern” wigeons 758 Introduction Eurasian wigeon (Mareca penelope) and American wigeon (Mareca americana) compose the “northern” wigeons (Peters et al. 2014), and are closely related waterfowl species that breed widely across the Palearctic and Nearctic, respectively. Their South American counterpart, the Chiloe wigeon (Mareca sibilatrix) is the sister species to the American wigeon, and collectively, these three species comprise the “wigeon clade” (Johnson and Sorenson 1999; Gonzalez et al. 2009). Eurasian and American wigeons are highly migratory species with vagrants making transhemispheric exchange and often occurring with flocks of other species when wintering, migrating, and even breeding (Krechmar and Кondratiev 2006; Withrow 2023). Eurasian wigeons are observed regularly in Canada and US, mainly within the Pacific Flyway and fewer within the Atlantic Flyway (Edgell 1984; Newton and Dale 1996; Campbell and Ryder 2013), and some have even been documented in Mexico and South America (Williams and Beadle 2003; Ramírez-Albores et al. 2021). American wigeons are occasionally observed in Western Europe and Northeast Asia (Makatsch 1980; Madge and Burn 1988; Mackinnon and Phillipps 2000; Votier et al. 2003; Lee et al. 2005; Krechmar and Кondratiev 2006; Matyushkov and Zdorikov 2018). Phenotypic hybrids of different generational crosses are not uncommon (Randler 2001) and have been reported in Europe, North America, and Eastern Asia (Merrifield 1993; Nechaev and Gorchakov 1995; Randler 2001; McCarthy 2006). Female hybrids are recorded much less frequently than males, and most probably overlooked (Gillham and Gillham 1996) due to the phenotypic similarity of females of these species. Both species are highly migratory with well-defined migratory pathways. American wigeon migrates along all four North American flyways: Pacific, Central, Mississippi, and Atlantic (Bellrose 1980), and it is one of the most numerous dabbling ducks in North America, with an estimated breeding population of 2.7 million birds (BirdLife International 2021a). Banding records of Eurasian wigeon (Ostapenko et al. 1997) confirm five biogeographic populations in Eurasia: Icelandic, European, Western Siberian, Eastern Siberian, and the Russian Far East with high migrant exchange between populations (5 to 30 %). The global population of Eurasian wigeon is estimated to number around 2.65–3.59 million individuals (BirdLife International 2021b). High levels of migration between continents likely leads to gene flow between populations within continents and, to some extent, between continents and across species boundaries. Indeed, studies focused on mitochondrial (mt) DNA control region sequences of Eurasian wigeon recovered weak phylogeographic structure, implying strong maternal gene flow and connectivity between distant populations in the Palearctic (Kulikova et al. 2019). American and Eurasian wigeon have diagnosable differences in male plumage and deeply divergent mtDNA haplogroups (Peters et al. 2005; Peters et al. 2014). However, some Eurasian wigeons share mtDNA haplotypes with American wigeons, likely resulting from introgressive hybridization (Peters et al. 2005; Peters et al. 2014; Kulikova et al. 2019). Moreover, comparing mtDNA sequences with those from 20 nuclear introns revealed prominent mito-nuclear discordance in which nuDNA differentiation (ФST = 0.046) was much lower than mtDNA divergence (ФST = 0.812; Peters et al. 2014). Such mito-nuclear discordance has been widely observed in the literature and across all tree of life taxa from protozoans to birds and mammals (DeRaad et al. 2023). In the previous study, the only nuclear locus with high ФST of 0.415 was the Z-sex chromosome linked intron 19 of chromo-helicase-DNA binding protein gene 1, CHD1Z (Peters et al. 2014). In birds, which have a ZW sex chromosome system (females are ZW), the pattern of higher divergence of Z chromosome in comparison to autosomes in sister species pairs is well known and widely discussed (reviewed in Irwin 2018). Sex chromosomes likely play a very important role in speciation because they house genes related to sexual selection and reproductive isolation (Charlesworth et al. 1987; Presgraves 2008). For instance, some Z-chromosome linked loci are linked to sexually selected plumage traits involved in mate selection (Toews et al. 2016; Campagna et al. 2017). A higher rate of differentiation of the Z chromosome relative to autosomes is also commonly explained by faster Z evolution due to genetic drift acting faster because of reduced effective population size of sex chromosomes (¾ that of autosomes). Lower recombination rates and the large Z effect disproportionately affect the fitness of hybrids (Irwin 2018; Payseur et al. 2018). In dabbling ducks, higher divergence of Z-linked loci in contrast to autosomal loci is reported for sexually dichromatic mallard (Anas platyrhynchos) and its monochromatic close relatives: Mexican duck (A. diazi), mottled duck (A. fulvigula), American black duck (A. rubripes), Chinese spot-billed duck (A. zonorhyncha; Lavretsky et al. 2015, 2019; Kulikova et al. 2022). Similarly, high Z-chromosome divergence was found between the sexually monochromatic grey teal Anas gracilis and the dichromatic chestnut teal A. castanea (Dhami et al. 2016), as well as between gadwall Anas strepera and falcated duck A. falcata (Dhami et al. 2018). The objective of this study was to test for genome-wide genetic differentiation between the American and Eurasian wigeon and among populations of each species. We sampled both species widely across the Palearctic and Nearctic and applied ddRAD-seq (double digest restriction associated DNA sequencing). We combined ddRADseq data with mtDNA control region variability to answer the following questions: (a) are American and Eurasian wigeon populations genetically structured across the Holarctic, (b) is there congruence between different types of markers: mtDNA, sex chromosomes and autosomes, (c) is the Z chromosome more divergent than the autosomes, and (d) what is the number and distribution of divergent loci between these species, especially on the Z chromosome. Vertebrate Zoology 75, 2025, 757–772 759 Methods Mitochondrial DNA: sample collection, sequencing and analyses 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 PV787184PV787239; 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 (dXY), 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. Nuclear DNA ddRAD-seq: sample collection, library preparation, sequencing and bioinformatics 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 Figure 1. Approximate breeding distributions of the “northern” wigeons (modified from Peters et al. 2005) with inset photographs of males of both species. Closed circles indicate sampling locations for Mareca penelope, open circles indicate sampling locations for M. americana. Wintering distributions are not shown, although both species are seasonal migrants, and several individuals were sampled from migrating and over-wintering populations (Tables S1–S2). Photos: Eurasian wigeon (M. penelope), Graeme Travers / Pexels; American wigeon (M. americana), Bryan Hanson / Unsplash. Kulikova IV et al.: Population genomics of “northern” wigeons 760 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-RADseq/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. Nuclear population structure 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 (dXY), 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 dXY, 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 Vertebrate Zoology 75, 2025, 757–772 761 diversity divided by an estimate of 1.1 for substitution rate ratio of Z vs. autosomes) to autosome diversity (πZ/πA). Outlier analysis and tests of selection We calculated the pairwise per locus values of ФST, dXY, 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. Results Genetic diversity and differentiation – mtDNA 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. Nuclear species diversity and differentiation 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 Table 1. Nucleotide diversity (π), absolute divergence (dXY), and relative divergence (ФST) in Eurasian wigeon, Mareca penelope (M.p.) and American wigeon, M. americana (M.a.) calculated with introgressed haplotypes included. Loci M.p. π M.a. π dXY ФST mtDNA 0.00278 0.00509 0.0328 0.880 A-loci 0.00792 0.00791 0.0082 0.039 Z-loci 0.00385 0.00386 0.0048 0.192 Kulikova IV et al.: Population genomics of “northern” wigeons 762 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 (R2 = 0.974). Eurasian and American wigeon had similar autosomal and Z-chromosome nucleotide diversities (Table 1). For autosomes absolute divergence (dXY) closely approximated the nucleotide diversities while for Z chromosome dXY 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 (dXY = 0.0082) than for the Z chromosome (dXY = 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 Figure 2. Haplotype network of mitochondrial DNA based on mtDNA control region sequences (659–661 bp) obtained from 205 Eurasian and American wigeons, Mareca penelope (M.p.) and M. americana (M.a.). Vertebrate Zoology 75, 2025, 757–772 763 Figure 3. Scatter plots of the first two principal coordinates for (a) 3092 autosomal and (b) 189 Z-linked ddRAD-seq loci for Eurasian wigeon, Mareca penelope (M.p.) and American wigeon, M. americana (M.a.). Figure 4. ADMIXTURE assignment probabilities for Eurasian wigeon, Mareca penelope (M.p.) and American wigeon, M. americana (M.a.) for (a) 3092 autosomal and (b) 189 Z-linked ddRAD-seq loci for K population values of 2 and 3 and one randomly chosen SNP from each locus. I – M.p. Siberia, II – M.p. Russian Far East, III – M.p. Western Beringia, IV – M.p. NA Pacific Flyway, V – M.p. NA Atlantic Flyway, VI – M.a. NA Pacific Flyway, VII – M.a. NA Central Flyway, VIII – M.a. NA Atlantic Flyway. Kulikova IV et al.: Population genomics of “northern” wigeons 764 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. Nuclear population differentiation 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). Outlier loci 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 BAYESFigure 5. Z chromosome Manhattan plot with significant outliers identified by the BAYESCAN analysis between Eurasian and American wigeons shown as diamonds (FDR = 0.01) and circles (FDR = 0.05). Vertebrate Zoology 75, 2025, 757–772 765 CAN were estimated to be under diversifying selection. Z-linked outliers had significantly higher estimates of dXY than non-outliers (0.0097 vs. 0.0044, two-tailed t-test P = 0.004), whereas dXY 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). Discussion Population genetic structure of “northern” wigeons 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 Figure 6. Haplotype networks for outlier SNPs on Z-chromosome (SNP positions are shown above each network, Z-loci names and SNP variants below each network) of Eurasian wigeon (dark fill), American wigeon (light fill), and American x Eurasian wigeon hybrid (white fill). Each haplotype network includes one (females) or two (males) alleles per individual. Kulikova IV et al.: Population genomics of “northern” wigeons 772 Supplementary Material 1 Figures S1–S5 Authors: Kulikova IV, Lavretsky P, McCracken KG, Zhuravlev YuN, Miroshnichenko IL, Correll AB, Peters JL (2025) Data type: .docx Explanation notes: Figure S1. Cross-validation errors for autosomal and Z-chromosome ADMIXTURE results for K 1–10. — Figure S2. STRUCTURE assignment probabilities for Eurasian wigeon (Mareca penelope) and American wigeon (M. americana) for 3092 autosomal and 189 Z-linked ddRAD-seq loci with number of populations K of 2–4. — Figure S3. Delta K values for autosomal and Z-chromosome STRUCTURE results for K = 1–8. — Figure S4. Comparisons of pairwise ΦST values based on autosomal ddRAD loci, Z-linked ddRAD loci and mtDNA. — Figure S5. Manhattan plots of autosomal chromosomes with significant outliers as red circles (FDR = 0.01) and yellow circles (FDR = 0.05). Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/ odbl/1.0). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/vz.74.e167908.suppl1 Supplementary Material 2 Tables S1–S8 Authors: Kulikova IV, Lavretsky P, McCracken KG, Zhuravlev YuN, Miroshnichenko IL, Correll AB, Peters JL (2025) Data type: .xlsx Explanation notes: Table S1. Sample information on “population”, location, NCBI accession number for mtDNA dataset. — Table S2. Sample information on “population” and location for ddRAD-seq dataset. — Table S3. Nucleotide diversity (dX), gene diversity (H) and Tajima's D values in populations of Eurasian Wigeon (M.p.) and American Wigeon (M.a.). Far East – Russian Far East, West Ber – Western Beringia, Pac NA – Pacific Flyway (North America), Atl NA – Atlantic Flyway (North America), Cent NA – Central Flyway (North America). — Table S4. A-loci pairwise values of nucleotide divergence (dxy, above diagobal) and gene flow (Fst, below diagonal) between populations of Eurasian wigeon (M.p.) and Americam wigeon (M.a) based on (a) 3092 ddRAD autosomal loci; (b) 189 ddRAD Z-loci; (c) 661 bp mtDNA. — Table S5. Numbers of ddRAD loci aligned to Mallard chromosomes (reference genome GCA_015476345.1). — Table S6. AMOVA analysis of Eurasian and American wigeons. — Table S7. BLASTn searches for the Z chromosome loci detected at the outlier SNP analysis. — Table S8. BLASTn searches for the autosomal loci detected at the outlier SNP analysis. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/ odbl/1.0). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/vz.74.e167908.suppl2