Mallard ducks are a ubiquitous and socio-economically important game bird in North America. Despite their generally abundant midcontinent population, Mallards in eastern North America are declining, which may be partially explained by extensive hybridization with human-released domestically-derived game-farm Mallards. We investigated the genetic composition of Mallards in the middle and lower Mississippi flyway, key wintering regions for the species. We found that nearly 30% of wild Mallards carried mitochondrial haplotypes derived from domestic mallards present in North America, indicating that the individuals had female game-farm Mallard lineage in their past; however nuclear results identified only 4% of the same sample set as putative hybrids. Recovering 30% of samples with OW A mtDNA haplotypes is concordant with general trends across the Mississippi flyway and this percentage was stable across Mallards we sampled a decade apart. The capture and perpetuation of OW A mtDNA haplotypes is likely due to female breeding structure, whereas reversal of the nuclear signal back to wild ancestry is due to sequential backcrossing and lower and/or declining admixture with game-farm Mallards. Future studies of wild ancestry of Mississippi flyway Mallards will benefit from coupling molecular and spatial technology across flyways, seasons, and years to search for potential transitions of Mallard populations with different genetic ancestry, and whether the genetic ancestry is somehow linked to an individual’s natal and subsequent breeding location.

Sample collection and DNA extraction.

A total of 369 and 130 Mallards were sampled from hunter-killed birds during autumn-winters of 2011–12 and 2020–21, respectively. Specifically, samples from the 2011–12 hunting season were collected in west-central Missouri and Mississippi, while we expanded geographical sampling for the 2020–21 hunting season to also include Arkansas and Louisiana (Figure 1A). In the end, we either conducted temporal comparisons between states sampled between decades (i.e., west-central Missouri and Mississippi) or among states within the same time period, and we acknowledge that sample sizes are unequal between the two types of analyses.

DNA was extracted for all 499 samples using a DNeasy Tissue Kit (Qiagen, Valencia California), following manufacturers’ protocols. We ensured DNA quality based on the presence of high molecular weight bands visualized using gel electrophoresis with a 1% agarose gel, and quantified using a Qubit 3 Flourometer (Invitrogen, Carlsbad, CA) having a minimum concentration of 20 ng/µL.

Mitochondrial DNA Sequencing and Analyses

We used primers L78 and H774 to PCR amplify and sequence ~625 bp of the mitochondrial control region (Sorenson & Fleischer 1996; Sorenson et al. 1999), following PCR reaction concentrations and thermocycler conditions described in Lavretsky et al. (2014b). The PCR products were visualized via agarose electrophoresis and then purified using either the QIAquick PCR Purification kit (Qiagen) or ExoSAP-IT (ThermoFisher). Clean PCR products were then sent for Sanger sequencing using the L78 primer on a 3130XL Genetic Analyzer at either the Core Laboratories of Arizona State University School of Life Sciences or the University of Texas at El Paso, Border Biomedical Research Center’s Genomic Analysis Core Facility. Sequences were aligned and edited using Sequencher version 4.8 (Gene Codes) and all sequences were subsequently submitted to GenBank (accession numbers TBD). We note that mtDNA control region sequences for reference wild and domestic mallards from previous studies were included in the analyses (Lavretsky et al. 2014a; Lavretsky et al. 2014b; Lavretsky et al. 2019; Lavretsky et al. 2020).

Relationships among mtDNA haplotypes were assessed through a median-joining network constructed with Network v. 5.0 (Bandelt et al. 1999). Given that mtDNA haplotypes can be permanently captured through maternal inheritance, we were particularly interested in individuals possessing OW A mtDNA haplotypes as this is a proxy to determine individuals in which their lineage at some point included a female captive-bred Mallard (Lavretsky et al. 2020). This was especially informative for early 2011–12 Mallard samples for which we were unable to assess nuclear variation (see below). Consequently, in addition to visualizing haplotype relationships, we examined OW A versus NW B haplotype ratios across space (i.e., 4 states) and time (2011–12 versus 2020–21).

Nuclear DNA ddRADseq library preparation, Sequencing, and Bioinformatics

For the 2020–21 Mallard samples, we followed procedures presented by Lavretsky et al. (2015), but with fragment size selection following Hernandez et al. (2021) to create multiplexed ddRAD-seq fragment libraries. Briefly, we enzymatically fragmented genomic DNA using SbfI and EcoRI restriction enzymes and ligated Illumina TruSeq compatible barcodes that permitted future de-multiplexing. We pooled libraries in equimolar concentrations, and 150 base pair (bp), single-end (SE) sequencing was completed on an Illumina HiSeq X at Novogenetics LTD (Sacramento, CA). Illumina reads were deposited in NCBI’s Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/sra; SRA data TBD).

We used the ddRADparser.py script of the BU ddRAD-seq pipeline (DaCosta and Sorenson 2014) to de-multiplex raw Illumina reads based on perfect barcode/index matches. As with mtDNA, previously published ddRAD raw sequence data generated using the same protocols were included in alignments and subsequent analyses, serving as reference wild Mallard (Lavretsky et al. 2019) and domestic Mallards (Lavretsky et al. 2020). Note that known game-farm Mallards and Khaki Campbells were included as known domestic lineages, with the latter one serving as a genetic signature of a “park” duck, another potential source of hybridization with wild Mallards (also see Lavretsky et al. 2020). All sequences were first trimmed or discarded for poor quality using trimmomatic (Bolger et al. 2014), and then remaining quality reads were aligned to a reference Mallard genome (Huang et al. 2013) using the Burrows Wheeler Aligner v. 07.15 (bwa; Li and Durbin, 2011). Next, samples were sorted and indexed in Samtools v. 1.6 (Li et al., 2009) and combined using the “mpileup” function with the following parameters “-c –A -Q 30 -q 30.” All steps through “mpileup” were automated using a custom in-house Python script (Python scripts available at https://github.com/jonmohl/PopGen). Next, we used VCFtools v. 0.115 (Danecek et al., 2011) to further filter VCF files for any base-pair missing >5% of samples, which also required a minimum base-pair sequencing depth coverage of 5X (i.e., 10X per genotype) and quality per base PHRED scores of ≥30 to be retained.

Relationships among mtDNA haplotypes were assessed through a median-joining network constructed with Network v. 5.0 (Bandelt et al. 1999). Given that mtDNA haplotypes can be permanently captured through maternal inheritance, we were particularly interested in individuals possessing OW A mtDNA haplotypes as this is a proxy to determine individuals in which their lineage at some point included a female captive-bred Mallard (Lavretsky et al. 2020). This was especially informative for early 2011–12 Mallard samples for which we were unable to assess nuclear variation (see below). Consequently, in addition to visualizing haplotype relationships, we examined OW A versus NW B haplotype ratios across space (i.e., 4 states) and time (2011–12 versus 2020–21).

To evaluate nuclear population structure, we used only the autosomal ddRAD-seq bi-allelic single nucleotide polymorphisms (SNPs). Prior to analyses, we used PLINK v. 0.70 (Purcell et al. 2007) to ensure that singletons (i.e., minimum allele frequency [maf] = 0.0038) and any SNP missing >5% of data across samples were excluded in each dataset. Additionally, we identified independent SNPs by conducting pair-wise linkage disequilibrium (LD) tests across ddRAD-seq autosomal SNPs (--indep-pairwise 2 1 0.5) in which 1 of 2 linked SNPs are randomly excluded if we obtained an LD correlation factor (r2) > 0.5. We conducted all analyses without a priori information on population or species identity.

First, we used the PCA function in PLINK to perform a principal component analysis (PCA). Next, we used ADMIXTURE version 1.3 (Alexander et al. 2009, Alexander and Lange 2011) to attain maximum likelihood estimates of population assignments for each individual, with datasets formatted for the ADMIXTURE analyses using PLINK and following steps outlined in Alexander et al. (2012). We ran each ADMIXTURE analysis with a 10-fold cross-validation, incorporating a quasi-Newton algorithm to accelerate convergence (Zhou et al. 2011). Each analysis used a block relaxation algorithm for point estimation and terminated once the change in the log-likelihood of the point estimations increased by <0.0001. We ran separate ADMIXTURE analyses that included all possible samples and another excluding Khaki Campbell mallards. Each analysis was run for K populations of 1 through 5, and with 100 iterations per each value of K. The optimum K in each analysis was based on the average of cross-validation errors across the iterations per K value; however, we examined additional values of K to test for further structural resolution across analyses. We used the R package PopHelper (Francis 2016) to convert ADMIXTURE outputs into CLUMPP input files at each K value, and to determine the robustness of assignments of individuals to populations at each K value with the program CLUMPP version 1.1 (Jakobsson and Rosenberg 2007). In CLUMPP, we employed the Large Greedy algorithm and 1,000 random permutations. Final admixture proportions for each K value and per sample assignment probabilities (Q estimates; the log-likelihood of group assignment) were based on CLUMPP analyses of all 100 replicates per K value. In addition to the above replicates, standard deviations were calculated under the optimum K population value based on 1,000 bootstraps as implemented in the ADMIXTURE program. Doing so permitted us to evaluate how sensitive our assignment probabilities were given our SNP dataset.