Complications of estimating hatchery introgression in the face of rapid divergence: A case study in brook trout (Salvelinus fontinalis)

Fish stocking has been utilized for over a century to offset extirpations or declines in abundance of many native species. These historical declines and hatchery contributions have led to uncertainty surrounding whether many contemporary populations are native, introgressed with hatchery sources, or entirely of hatchery origin. Such uncertainty is problematic for the conservation of native biodiversity as it hampers management agencies’ ability to prioritize the conservation of indigenous locally adapted populations. Fortunately, genetic and genomic tools have allowed researchers to investigate these questions, often through the use of clustering or assignment approaches that are predicated on identifiable and consistent divergence between native populations and hatchery sources. Here, we apply these methods to genomic data from 643 brook trout (Salvelinus fontinalis) originating from 13 wild populations and an exogenous hatchery strain to investigate the extent of historical extirpations, hatchery contributions, and processes affecting population structure in a small area of the previously unglaciated Driftless Area of Wisconsin, USA. Results from these analyses suggest that wild populations in this region are genetically distinct even at small spatial scales, lack hydrologically associated population structure, rarely exchange gene flow, and have small effective population sizes. Further, wild populations are substantially diverged from known hatchery strains and show minimal evidence of introgression in clustering analyses. However, we demonstrate through empirically informed simulations that distinct wild populations may potentially be hatchery-founded and have since diverged through rapid genetic drift. Collectively, the lack of hydrological population structure and potential for rapid drift in the Driftless Area suggest that many native populations may have been historically extirpated and refounded by stocking events, and that commonly used genomic clustering methods and their associated model selection criteria may result in underestimation of hatchery introgression in the face of rapid drift.

Methods
Sampling

Fin clips were collected from 13 wild brook trout populations in Wisconsin and the St. Croix Falls hatchery strain with a target sample size of 50 fish per population.  All wild samples were collected via backpack electrofishing from 2015-2019, St. Croix Falls hatchery strain samples were collected from broodstock maintained at the St. Croix Falls State Fish Hatchery in 2019, and all fin clips were preserved in 95% ethanol prior to DNA extraction.  Our primary objective for these collections was to examine putative fine scale population structure of seven populations originating from the Pine River watershed in the Driftless Area of Wisconsin.  However, given the diverse stocking sources used in the Wisconsin, we also included four additional populations (West Branch of Mill Creek, Lowery Creek, South Fork of the Hay River, and Lawrence Creek) that have historically or currently serve as wild brood or translocation sources.  In addition, we also sampled two additional populations near the Pine River watershed in the Driftless Area.

Library preparation

Total genomic DNA was extracted using the Promega Wizard Genomic DNA purification kit (Promega Corporation, Madison, Wisconsin).  DNA concentrations were estimated using Quant-iT dsDNA PicoGreen assays (ThermoFisher Scientific, Waltham, Massachusetts) and a BioTek Synergy HTX Plate Reader (Agilent Technologies, Santa Clara, California).  Approximately 100ng of each sample was digested in 2µl reactions with restriction enzyme SbfI (New England Biolabs, Ipswich, Massachusetts).  Barcodes were then ligated to digested samples and sheared to fragment lengths of 300-500bp using 12-14 30-second cycles of a Q500 sonicator (Qsonica, Newtown, Connecticut).  These fragments were subsequently bound to Dynabeads M-280 Streptavidin magnetic beads (ThermoFisher Scientific, Waltham, Massachusetts) and washed to remove non-target sequences.  We then size-selected these fragments using AMPure XP beads (Beckman Coulter, Brea, California) targeting an insert size of 250bp and enriched our libraries following the NEBNext Ultra DNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, Massachusetts) using plate-specific i7 barcodes and 12 cycles of PCR.  Enriched and size-selected libraries were visually confirmed on a 2% agarose E-Gel (ThermoFisher Scientific, Waltham, Massachusetts) before a final AMPure XP bead purification.  Final libraries were quantified using Quant-iT dsDNA PicoGreen assays, pooled to equimolar concentrations, and sequenced on an Illumina Novaseq 6000 using 1.4 lanes of 2x150bp S4 chemistry (Novogene, Sacramento, California).

Sequencing, genotyping, and filtering

Raw reads were processed using stacks v2.3 (Catchen et al. 2011, 2013).  Briefly, sequences were demultiplexed by i5 and i7 barcodes, filtered for SbfI restriction sites, and trimmed to lengths of 140bp (i.e., -e SbfI -c -q -r -t 140 --filter_illumina --bestrad).  Individual-specific reads passing these filters were aligned with ustacks and parameters previously shown to avoid under- and over-merging loci from salmonid RAD datasets (i.e., --disable-gapped -m 3 -M 5 -H --max_locus_stacks 4 --model_type bounded --bound_high 0.05).  A catalog of loci was then built using all resultant stacks and the cstacks module (-n 3 -p 6 --disable_gapped).  Individual-specific stacks were then matched against the locus catalog using the sstacks module (--disable_gapped).  Paired end reads were then imported and all data were transposed by locus using tsv2bam.  Paired end contigs were assembled and SNPs were called with the gstacks module.  Resultant individual genotypes were exported in variant call format for loci genotyped in at least 30% of individuals using the populations module (-r 0.3).  We further used vcftools v0.1.16 to filter this dataset by sequentially removing loci with minor allele counts less than three, loci that were not successfully genotyped in ≥70% of individuals in each population, and individuals that were not genotyped at ≥50% of loci (--mac 3, --exclude-positions, and --remove, respectively; Danecek et al. 2011).  Paralogous loci, resultant from a recent whole genome duplication event in salmonids, were identified and removed using HDPlot with absolute read deviations greater than 8 and observed heterozygosity less than 0.5.  Finally, as some of our downstream analyses assume independence of loci, we further filtered this dataset by retaining only one SNP per contig and chose to retain the SNP with the highest minor allele frequency in these instances.