DOI: 10.5061/dryad.4xgxd25d4
title: Low levels of hybridization between domestic and wild Mallards wintering in the Lower Mississippi Flyway

description: 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.
description: <em>Sample collection and DNA
 extraction.</em> 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. <em>Mitochondrial DNA Sequencing and
 Analyses</em> 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 <em>TBD</em>). 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). <em>Nuclear DNA ddRADseq library
 preparation, Sequencing, and Bioinformatics</em>
 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
 <em>TBD</em>). We used the
 <em>ddRADparser.py</em> 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 &gt;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.
description: 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 &gt;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) &gt;
 0.5. We conducted all analyses without <em>a priori</em>
 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 &lt;0.0001. We
 ran separate ADMIXTURE analyses that included all possible samples and
 another excluding Khaki Campbell mallards. Each analysis was run for
 <em>K</em> populations of 1 through 5, and with 100 iterations
 per each value of <em>K</em>. The optimum
 <em>K</em> in each analysis was based on the average of
 cross-validation errors across the iterations per <em>K</em>
 value; however, we examined additional values of <em>K</em> 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 <em>K</em> value, and to determine the
 robustness of assignments of individuals to populations at each
 <em>K</em> 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 <em>K</em> value and per sample assignment probabilities
 (Q estimates; the log-likelihood of group assignment) were based on CLUMPP
 analyses of all 100 replicates per <em>K</em> value. In
 addition to the above replicates, standard deviations were calculated
 under the optimum <em>K</em> 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.