Invasion genetics of the silver carp Hypophthalmichthys molitrix across North America: Differentiation of fronts, introgression, and eDNA metabarcode detection

In the 1970s, the introduced silver carp Hypophthalmichthys molitrix (which is indigenous to eastern Asia) escaped from southern U.S. aquaculture to spread throughout the Mississippi River basin, and since has steadily moved northward. This large, prolific filter-feeder reduces food availability for other fishes. It now has reached the threshold of the Laurentian Great Lakes, where it likely will significantly impact food chains and fisheries. Our study evaluates population genetic variability and differentiation of the silver carp using 10 nuclear DNA microsatellite loci, and sequences of two mitochondrial genes–cytochrome b and cytochrome c oxidase subunit 1, along with the nuclear ribosomal protein S7 gene intron 1. We analyze population samples from: two primary Great Lakes’ invasion fronts (at the Illinois River outside of Chicago, IL in Lake Michigan and in the Wabash River, which leads into the Maumee River and western Lake Erie), the original establishment “core” in the Lower Mississippi River, and expansion areas in the Upper Mississippi and Missouri rivers. We analyze and compare our results with bighead and other invasive carps, and cyprinid relatives. Results reveal that the silver carp invasion possesses moderate levels of genetic diversity, with more mtDNA haplotypes and unique microsatellite alleles in the “core” Lower Mississippi River population, which also diverges the most. The two invasion fronts also significantly genetically differ. About 3% of individuals (including all populations except the Illinois River) contain a unique and very divergent mtDNA haplotype, which likely stems from historic introgression in Asia with female largescale silver carp H. harmandi. The nuclear microsatellites and S7 sequences of the introgressed individuals do not differ from silver carp and are very distant from bighead carp. These sequence variation data are employed to design and evaluate a targeted high-throughput metabarcoding sequence assay that identifies and distinguishes among species of invasive carps (i.e., silver, bighead, grass, black, and common carps, along with goldfish), as well as native cyprinids, using cytochrome b. Our assay further differentiates among selected silver carp haplotypes (including between H. molitrix and H. harmandi), for use in population genetics and future analyses of spread pathways. We test and evaluate this assay on environmental (e)DNA water samples from 48 bait shops in the Great Lakes’ region (along the Lake Erie, Lake St. Clair, and Wabash River watersheds), using positive and negative controls and custom bioinformatic processing. Test results discern silver carp eDNA in four of the shops–three in Lake Erie and one in the Wabash River watershed–and bighead carp from one of the same Lake Erie venues, suggesting that retailers (who often source from established southerly populations) comprise another introduction vector. Our overall findings thus provide key population genetic and phylogenetic data for understanding and tracing introductions, vectors, and spread pathways for silver carp, their variants, and their relatives.


Introduction
Discerning the population genetic trajectories of nonindigenous species invasions can enhance our overall understanding of the evolutionary adaptations and changing ecological community dynamics governing today's ecosystems [1][2][3]. Establishments by invasive species comprise accidental experiments to ground-truth evolutionary and ecological theory with reality [4][5][6]. The genetic variation of invasive populations often influences their relative success and persistence in new habitats, including colonizing, reproducing, spreading, and overcoming biotic resistance [2,7,8].
Traditional invasion theory predicted that most introduced populations would be characterized by low genetic diversity due to founder effect, which would limit adaptive potential [9][10][11]. More recently, some invasions founded by large numbers of introduced propagules, multiple events, and multiple sources have been found to possess as high or even higher levels of population genetic variability than those in native regions, due to admixture [12][13][14]. For example, ballast water introductions from Eurasian sources into the North American Laurentian Great Lakes of the zebra mussel Dreissena polymorpha (Pallas, 1771) [15,16], the quagga mussel D. rostriformis (Deshayes, 1838) [15,16], and the round goby Neogobius melanostomus (Pallas, 1814) [5,[17][18][19] all possessed relatively high population genetic diversities and significant population divergences across their introduced ranges, on par with those of native populations. These invasions stemmed from large numbers of introduced propagules and multiple introduction events, which have been traced to several founding sources [5,[15][16][17][18][19]. Such genetic variability may enhance adaptive potential of invasions (see [11]).
The "leading edge" hypothesis postulates that expansion populations at an invasion's front (s) should possess less genetic variability than those at the invasion's core (i.e., the original successful founding area) [20,21]. According to that hypothesis, individuals at the front(s) likely would be adapted for dispersal and high reproductive output [22], as they would experience low population density, greater resource availability, and less competition, thereby enhancing their relative reproductive success [23,24]. This might lead to population genetic differences between the expansion fronts versus longer-established core areas, due to drift and/or selection. However, studies of these long-term trends across the spatial and temporal courses of invasions are relatively rare in the literature [5,6,22].
Closely related species, which may interact competitively, sometimes are introduced together or in close succession, as in the mid-1980s cases of zebra and quagga mussels in the Great Lakes [25][26][27]. The quagga mussel typically overtakes the zebra mussel over time, with populations becoming dominated by the former except in shallow areas [28,29]. In another Great Lakes' example, the Eurasian round goby and the tubenose Proterorhinus semilunaris

Invasion ecology of the silver carp and its relatives
Silver and bighead carps are prolific filter-feeders that significantly alter both the numbers of plankters and their community compositions, reducing food for sport and commercial fishes [ [62][63][64][65]. Silver and bighead carps often swim just below the surface, and can travel in large schools (either as single species or together) [66]. In their native ranges, they mature at age 4-8 years, but as early as 2 years old in North America, with each female laying up to 5 million eggs annually [67]. Their reproductive season is longer in the Mississippi River system than in their native habitats [68]. The eggs are buoyant, hatch in one day, and the larvae readily disperse with currents [66]. Adults live to 20 years, reaching >1 m and 27.3 kg [66]. Modeling studies suggest that even a few reproductive individuals could establish successful populations in the Great Lakes [69].
The silver carp's current North American geographic range is depicted in Fig 1 [70]. Throughout the Mississippi River drainage, its populations extend upstream and downstream of 23 locks and dams (three in the Arkansas River, seven in the Illinois River, eight in the Mississippi River, and five in the Ohio River). At present, there are two potential man-made impediments to silver and bighead carps reaching the Great Lakes basin, of which the first is an electric barrier in the Chicago Area Waterway System that separates the Illinois River from Lake Michigan [71]. This "barrier" frequently is breached by small fishes and larger ones traveling in the wakes of large boats [71]. In 2016, a 2.7 km long, 2.3 m-tall earthen berm was completed in Eagle Marsh at Fort Wayne, IN, between the Wabash and Maumee Rivers (the latter leads into Lake Erie) [72]. This wetland area frequently experienced flooding and connection between the two watersheds, and previously was separated by just a chain link fence, which smaller fishes (and young carps) could readily swim through [72]. The present investigation analyzes populations at these two invasion front areas (i.e., the Illinois and Wabash Rivers), in comparison to genetic variation throughout much of the range. The issue of silver and bighead carps potentially entering and establishing in the Great Lakes is of great concern to managers, the commercial and sport fishing communities, ecologists, and many of the general public [73,74].

Research objectives
Little is known about the population genetic and genomic patterns of the silver carp in North America (e.g., [75,76]). Thus, our objectives are to evaluate genetic diversity and its distribution, and to test for possible population structure. Population samples include the seat of its original establishment in the Mississippi River basin, expansion areas in the Upper Mississippi and the Missouri rivers, and the two most likely invasion fronts leading to the Great Lakes-the Illinois and Wabash rivers. We test invasion genetics theory, including possible founder effect and correspondence to the leading edge hypothesis at the front areas. We further analyze whether population genetic divergence and differentiation occurs across the range, including between the fronts. We evaluate genetic variation using 10 nuclear (n)DNA microsatellite (μsat) loci and sequences from two mitochondrial (mt) genes, which include 992 nucleotides (nts) of cytochrome b (cytb) and 549 nts of the cytochrome oxidase subunit I (COI). We compare representative subsets encompassing all mt and μsat variants with nDNA sequence variation from the single copy ribosomal protein S7 gene intron 1. We relate our results to other genetic studies of silver carp, its relatives, and other exotic species.
Based on those baseline genetic variability data, we design a targeted cytb metabarcode assay to discern and identify silver carp, bighead carp, and their cyprinid relatives, in order to facilitate early detection and tracking of species and possible hybrids at all life stages, using Illumina MiSeq high-throughput sequencing (HTS). The assay is designed to also simultaneously distinguish among silver carp haplotype variants, for use in population comparisons. We test our assay's application for discerning and identifying taxa and possible introduction vectors using eDNA-containing water sampled from 48 bait shops in the Lake Erie, Lake St. Clair, and Wabash River watersheds. Overall, our investigation aims to provide and apply population genetic, evolutionary phylogenetic analysis, and systematic biological knowledge of genetic variation in silver carp and its cyprinid relatives to interpret and understand invasion dynamics, pathways, and progression.

Sample collection
Tissue samples from silver carp individuals were collected by us, our laboratory members, federal or state agencies, and university collaborators (see Acknowledgments). Collections were made under their state or federal collection permits, using their protocols, and under the University of Toledo's IACUC protocol #205400, "Genetic studies for fishery management" to CAS. Samples included 11 collection areas (Table 1), representing much of the silver carp's North American range (Fig 1). Bighead carp also were sampled, along with black, common, and grass carps to provide comparative data from related species. Morphological characters, including gill-raker structure, were used to distinguish silver from bighead carps and discern possible hybrids (see [33]). Samples were labeled, stored in 95% EtOH, and archived at CAS' G3 lab in the NOAA Pacific Marine Environmental Laboratory.
Genetic diversity measures, including the number of alleles per locus (N A ), observed heterozygosity (H O ), and allelic richness (A R ; adjusted for sample size using rarefaction) were calculated in FSTAT v2.9.3.2 [87], their standard errors (±S.E) with EXCEL (Microsoft, Redmond, WA), and the significance values (p<a) adjusted with standard Bonferroni correction. [85]. Significant differences in H O and A R were evaluated using paired two-tailed t-tests in R v3.2.1 [88]. Number of private alleles (N PA ) per locus, i.e., those appearing unique in a population sample, were identified with CONVERT v1.31 [89]. The percentage of private alleles (P PA ) was the number of private alleles in a given sample divided by its total number of alleles, using rarefaction representation in ADZEv1.0 to adjust for sample size disparity [90]. COLONY v2.0.6.1 was used to test for the respective presences of full and half siblings in the samples [91]. The possible influence of selection was evaluated using the outlier method of Beaumont and Nichols 1996 [92] in LOSITAN [93].
Pairwise genetic divergences were calculated between all population samples, and separately between the two invasion fronts at the Illinois and Wabash Rivers with the F ST analog θ ST [94] in FSTAT, which is regarded as appropriate for analyses of high gene flow species, small sample sizes, and unknown number of subpopulations [95][96][97], and to facilitate comparisons with other studies. Since F-statistic estimates assume a normally distributed data set and may be influenced by sample size [98], we additionally conducted pairwise exact tests of differentiation (χ 2 ) in GENEPOP, using MCMC chains of 10,000, 1000 batches, and 10,000 iterations. Probabilities for both types of pairwise genetic divergences were assessed with sequential Bonferroni correction [99].
Microsatellite DNA variation for all individuals that possessed a highly divergent mtDNA haplotype ("H") was statistically compared against a group that included all of the other haplotypes (see Results and Table 1), to better match the size of the comparative datasets. We also ran the μsat variation of "H" group against the data set containing all other individuals (i.e., the entire data set).

DNA sequence analyses
We evaluated mtDNA sequence variation in the cytb and COI regions, both separately and together as concatenated sequences (Table 1), adding available complete gene sequences from NIH GenBank (https://www.ncbi.nlm.nih.gov/genbank/; detailed in S1 Table). Haplotype relationships are depicted with TCS haploptype networks [107] in POPART (http://popart.otago. ac.nz), in reference to four bighead carp haplotypes. ARLEQUIN calculated numbers of haplotypes (N H ), haplotypic diversity (h), numbers and proportions of private haplotypes (N PH and P PH ) per sample, and pairwise divergences using θ ST and exact tests [108]. All of the highly divergent ("H") haplotype silver carp individuals and representatives of all other haplotypes and populations also were sequenced for the nDNA S7 intron 1 (see Results and Table 1). The purpose was to determine whether their nDNA genome also was genetically divergent.
Relative percentage of mtDNA concantanted mtDNA haplotypes per population was illustrated using a stacked bar graph, drawn with R. Evolutionary relationships among silver carp sequence variants were evaluated with Bayesian phylogenetic trees in MRBAYES v3.2.6 [109], in comparison with bighead carp, grass carp, black carp, and common carp, with the latter as the selected outgroup. These used the GTR + I + Γ model of substitution, a relaxed molecular clock, the rate at which the variance of the effective branch length increased over time, and a default prior exponential distribution rate = 10.0. The MCMC was run for 100,000 generations to calculate support values for branch nodes. Consensus trees were visualized with FIGTREE v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

Design and testing a diagnostic assay for silver, bighead, and other invasive carps
Using the mtDNA sequence variation discerned for the cytb gene, we developed primers from conserved aligned sequence regions (see Results) that bracketed a diagnostic region of sequence variation, which differentiated among silver carp (and many of its haplotypes), bighead carp, and other invasive carp species for use in a targeted metabarcode HTS assay. The purpose was to facilitate detection, species identification, population variation, and quick, accurate differentiation in field samples.
Our assay was used to test for and evaluate the possible cryptic presence of silver carp (including haplotype variants), bighead carp, and other invasive cyprinid species in retail bait shops, as a pilot study for a larger ongoing project. We evaluated water samples and live bait from 48 bait shops located in the Lake Erie, Lake St. Clair, and Wabash River watersheds in 2016 and 2017. Approximately two dozen bait fishes were purchased from each shop, which were immediately sacrificed in the parking lots using our IACUC protocol. The water containing the fishes was drained into a sterile container, placed on ice, and then frozen at -80˚C in the laboratory.

Conducting the targeted HTS assay
We briefly thawed the water samples on ice, and centrifuged 250 ml from each at 4500 rpm and 4˚C, for 45 min. The supernatant was discarded and the pellet, containing extra and intracellular DNA and debris, was resuspended in 95% EtOH and stored at -20˚C. DNA was extracted within 24 hr with Qiagen DNeasy kits following the manufacturer's protocol, except that two AW 1 and 2 buffer washes were performed. The Zymo OneStep PCR Inhibitor Removal Kit (Zymo Genetics, Seattle, WA) was used to alleviate PCR inhibition.
Primers contained spacer inserts to increase the diversity of the HTS libraries [110]. 25μl PCRs comprised 1X Radiant TAQ Reaction Buffer (Alkali Scientific Inc., Ft. Lauderdale, FL), 3mM MgCl2, 1mM total dNTPs, 0.6mM of each primer (with spacer inserts and an Illumina, MiSeq sequencing primer tail), and 1.25 units of Radiant TAQ polymerase. Conditions were 2 min at 95˚C, followed by 40 cycles of 95˚C for 45 sec, 56˚C for 30 sec, and 72˚C for 45 sec, capped by 3 min at 72˚C. The final PCR used the prior step's column-cleaned product as its template (2μl) and incorporated Nextera paired end indices (Illumina, kit FC-121-1011), which included P5 and P7 adaptor sequences for binding the prepared library onto the flowcell's surface. Adding unique combinations of forward and reverse indices permitted multiple samples to be pooled together in each lane. Each reaction contained no template controls (NTCs), of which only those free of NTC amplification were used for libraries (see below). After column clean up, products were sized and quantified by us on a 2100 Bioanalyzer (Agilent Technologies), and concentrations of pooled products were measured with a Qubit fluorometer (Invitrogen). Pooled samples were run with 2X 300 base pair V3 chemistry. An additional 40-50% PhiX DNA spike-in was added to improve data quality of low nt diversity samples [110]. Each HTS run contained positive control mixtures of equal volumes of DNA extractions from 10 marine fish species (not present in the Great Lakes), which had been previously Sanger sequenced, to estimate sequencing error and help validate results [111]. Illumina MiSeq was performed by Ohio State University's Molecular and Cellular Imaging Center in Wooster, OH, USA (http://mcic.osu.edu/).

High-throughput sequencing data analyses
Custom bioinformatic scripts in PERL v5.26.1 trimmed primers and removed any sequences having a different spacer insert (which might result from index hopping if a few sequences with a different index became incorporated into an HTS library [112]), or were of the wrong length (i.e., short primer dimer sequences in the FASTQ files of the HTS runs; [113]). Trimmed sequences then were merged, dereplicated (grouped by 100% sequence similarity), and any chimeras eliminated with DADA2 [114], which employed a denoising algorithm with a Poisson model to calculate error probabilities at each nt position. This removed variant nts having a lower percentage representation than the predicted threshold. The resultant screened sequences are termed amplicon sequence variants (ASVs). Unexpected ASVs remaining in the positive controls after merging in DADA2 would be attributable to PCR and/or sequencing error, or sequence-to-sample mis-assignment (via index hopping). The greatest percent representation of any single erroneous ASV was set as the cutoff. If an ASV in a sample occurred at a lower frequency than this cutoff, then it was removed. Remaining ASVs were queried using the Basic Local Alignment Search Tool (BLAST; https://www.ncbi.nlm.nih.gov/) from the command line against a custom database containing cytb sequences of all fishes native to the Great Lakes basin [115] and non-native species that have been documented or have been predicted to possibly be introduced in the future [116]. All BLAST results with the lowest e value (best match) per species were summarized with a custom PERL script that also grouped together multiple detections in a sample. Samples with positive silver carp detections (see Results), were re-sequenced on a separate run using spacer inserts that did not share a forward or reverse index, markedly reducing possibility of undetectable index hops.

Microsatellite diversity, divergence, and population structure
Results indicated no linkage disequilibrium, null alleles, or selection. A single locus (Hmo-B4) that did not conform to Hardy-Weinberg equilibrium expectations was removed from further analyses. There were no differences in genetic diversity and θ ST values when population samples were pooled regionally (e.g., all samples from the Lower Mississippi together, N = 60) versus separately for individual sampling locations ( Table 2; Fig 1). Thus, our results are presented with regional pooling to best depict overall geographic patterns across the range, encompassing larger sample sizes. Four pairs of full silver carp siblings were identified by COLONY, whose inclusion or exclusion did not significantly alter the θ ST values; we presented the complete data here. One individual from three of the sibling pairs was collected in the Illinois River, and both individuals of the fourth pair also occurred there. Of the three located in separate populations, the other siblings respectively were found in the Lower Mississippi River, Upper Mississippi River, and the Missouri River. A total of 509 pairs of half siblings were identified, constituting 0.5% of the possible pairs and representing all possible combinations of population areas. Statistically fewer of the half sibling pairs (t-test p = 0.008 � ) co-occurred within the same population (19.8±4.85) versus in different populations (41.0±4.07), indicating high dispersal.
A total of 70 μsat alleles were discerned for silver carp (Table 2), with the most occurring at the invasion fronts, including 57 in the Illinois River and 58 in the Wabash River, versus 54 for the Lower Mississippi River population. Observed heterozygosity was 0.60 across all populations, ranging from 0.57 in the Illinois River, to 0.60 in the Wabash River and the Upper and Lower Mississippi rivers; these values did not significantly differ. Allelic richness appeared highest in the Wabash River front population at 6.3 and next highest in the Lower Mississippi River at 5.9, ranging down to 5.7 in the Illinois and Upper Mississippi rivers. No significant differences in allelic richness occurred among the populations. A total of 15 private alleles were identified, constituting 21% of the overall alleles (Table 2). Of these, the Illinois River invasion front possessed the most (seven, totaling 12% of its alleles), followed by the Lower Mississippi River (four at 7%), the Wabash River (three at 5%), and the Missouri River (one at 2%).
Significant genetic divergences distinguished all but just two pairs of populations, assessed by exact tests, indicating overall population structure ( Table 3). The two exceptions were similarities between the Upper Mississippi River versus populations in the Missouri and the Illinois rivers. Analyses using θ ST divergences revealed four significant differences, between the Lower Mississippi River population versus the Missouri, Illinois, and Wabash rivers, as well as between the Upper Mississippi and Wabash rivers (Table 3A). The two invasion front populations, the Illinois and Wabash rivers, significantly diverged in all analyses (Table 3B).
Silver carp that had a highly divergent mtDNA haplotype (designated "H"; see below section, "Mitochondrial and nuclear DNA sequence haplotypes") did not significantly differ in μsat composition from other silver carp (including analyses against all other individuals and for a subgroup that included representatives of all other haplotypes (see Table 1), with θ ST and exact tests.
The 3-d FCA (Fig 2) explained 88.54% of the data, distinguishing among all populations ( Table 3). All appeared widely separated, with the most divergence among the Upper Mississippi, Illinois, and Wabash river samples. STRUCTURE and STRUCTURE HAR-VESTER best-supported K = 2 population groups (not shown), for which distribution plots did not reveal appreciable structure (i.e., all individuals and samples showed approximately equal assignment to both). GENECLASS2 results for the entire data set also discerned overall mixed assignments. However, the GENECLASS2 assignment test between the two invasion front populations, the Illinois versus the Wabash rivers, indicated that 62% and 98% of their individuals respectively self-assigned. Individuals from the Wabash River population thus showed very high self-assignment.

Mitochondrial and nuclear DNA sequence haplotypes
In total, we analyzed 233 mtDNA cytb sequences, along with 248 COI, and 226 concatenated (COI and cytb) sequences for silver carp (Tables 1 and S1). Eight silver carp haplotypes occurred in North America (lettered "A-H" on Fig 3). We analyzed two other cytb haplotypes from GenBank, one from the Yangtze River in China (Accession AF051866, lettered "Q") and another from the Black River in Russia (AB198974, lettered "R"), which were not available for our sequencing. Genetic relationships among the haplotypes are depicted in Fig 4. The most haplotypes (seven) occurred in the Lower Mississippi River population (which lacked "F" alone) (Fig 3, Tables 2, 4 and 5). Population samples overwhelmingly were dominated by the common haplotypes "A" and "B" (Fig 3), which together comprised >90%, occurred in all populations, and differed by just a single nt (Fig 4). Haplotype "B" was less abundant in the Missouri River than in the other populations. All populations but the Illinois River additionally contained haplotype "H". Haplotypes "C", "D", and "E" were found in the Lower Mississippi River population alone; these were singletons that evolutionary analyses depicted as most closely related to common haplotype "A" (Figs 4 and 5A). Haplotypes "D" and "E" shared a closer relationship with each other. Haplotype "F" (a single individual) was discerned in the Upper Mississippi River population alone (Figs 3 and 4).
The invasion front populations of the Illinois and Wabash Rivers contained just three haplotypes each, with the Illinois River containing "G" (as well as a greater abundance of "B") and the Wabash River with "H" (Figs 3 and 4). The sole differentiations at the population level in mtDNA pairwise divergence tests were between the Illinois River versus the Missouri and the Lower Mississippi rivers, respectively. (Table 4). Haplotypic diversity across all of our samples was 0.54 and did not significantly differ among populations (Table 2), being highest in the Lower Mississippi River (0.61±0.04), followed by the Upper Mississippi River, and lowest in the Missouri River (0.48±0.07). Most silver carp haplotypes differed by 1-9 nt substitutions from haplotype "A", with the exception of "H" (whose divergence was much greater; Fig 4). Haplotypes exclusively found outside of North America, i.e., "R" and "Q", were closest in sequence relationship to our rare "F" and "G" haplotypes (Fig 4).
The haplotype networks revealed a distinctly divergent haplotype "H", which occurred in all populations except for the Illinois River, and was 47 nt substitutions from the closest silver carp haplotype and 76 from the bighead carp (Fig 4B). The "H" haplotype constituted a basal branch outside of the silver carp clade on the Bayesian phylogenetic tree using concatenated mtDNA sequences, with 1.00 posterior probability support (Fig 5A). Individuals with the "H" mtDNA haplotype did not vary from other silver carp in their nDNA S7 intron sequences, and lacked a common nDNA genotype or pattern, as shown on the Bayesian phylogenetic tree (Fig 5B). There was no differentiation of those with "H" from the 14 other silver carp individuals examined, which included all other mtDNA haplotypes ("A-G"), i.e., four randomly selected ones with "A", four with "B", the singletons "C", "D", "E", and "F", and both "G" individuals ( Table 1). Many of them shared identical sequences for the S7 intron with other "H" individuals and other silver carp haplotypes, while others differed in their S7 sequences alone (S2 Table). For example, 12 individuals possessed the same S7-"1" sequence, including four that had "H", two with haplotype "A", and one each with haplotypes "B-G". Individuals possessing sequences S7-"2" and -"3" each contained one individual of "A", whereas three individuals with S7-"4" had "B". A silver carp having the S7-"5" sequence Silver carp invasion genetics contained "G", and one with the S7-"6" sequence contained "H" (GenBank Accessions MH938833-43; Fig 5B and S2 Table). Note that one of the "H" individuals did not amplify for S7, and five of the six were analyzed. Thus, the mtDNA "H" variant showed no similar nDNA divergence, in either S7 sequences or μsat allelic frequencies (detailed in above section, "Microsatellite diversity, divergence, and population structure").
The Bayesian evolutionary trees using mtDNA ( Fig 5A) and nDNA sequences (Fig 5B) both showed strong support for the respective silver carp and bighead carp species clades, which each had posterior probabilities of 1.00. There were 21 nt nDNA S7 substitutions separating the silver and bighead carp species (Fig 5B).

Targeted HTS assay for silver and other carps
The invasive carp assay that we designed, developed, and tested here used the forward primer 5'-TGATGAAAYTTYGGMTCYCTHCTAGG-3' and reverse primer 5'-AARAAGAATGATGC YCCRTTRGC-3' to amplify a 135 nt region of the cytb gene that begins at nt114 (Table 5). This region can discern all invasive carp species and many native cyprinids, as well as selected  Table). Intra-specifically, it differentiates among silver carp haplotypes "B", "D", "E", and the highly divergent "H", along with "A/C" and "F/G" (with the slash indicating that "A" and "C" group together, as well as "F" with "G").
PCRs detected eDNA in water from 48 bait shops sampled in the Lake Erie, Lake St. Clair, and Wabash River watersheds in summers 2016 and 2017. We obtained 11,538,299 reads for these libraries (mean per eDNA water sample = 122,109±5,982 reads), of which 7,998,256 (mean = 81,615±3,021) were verifiable; i.e., they contained both primers, the correct spacer (indicating it was not an index hop), and were the correct length (were not attributable to primer dimer; S4A Table). Of those, 6,279,133 (mean = 66,302±2,876) were successfully merged using DADA2 (S4A Table). Following bioinformatic filtering (including removal of ASVs below the cutoff for the positive controls), our assay identified eDNA of several native cyprinids, including emerald shiner Notropis atherinoides Rafinesque, 1818, rosyface shiner N. rubellus (Agassiz, 1850), spottail shiner N. hudsonius (Clinton, 1824), fathead minnow Pimephales promelas Rafinesque, 1820, and golden shiner Notemigonus crysoleucas Rafinesque, 1820. Among invasive cyprinids, goldfish eDNA was detected in seven samples, grass carp in four, and common carp in five. Silver carp eDNA occurred in four shops after re-sequencing the libraries that had positive detections. Bighead carp eDNA was found in a single shop (which also contained silver carp eDNA), but whose ASV frequency (0.07% of reads) fell below the cutoff (0.19%). However, its haplotype matched our bighead carp haplotype "L", which differed from the silver carp "A" and "B" haplotypes, by 14 and 15 nts respectively (S3 Table). This large number of nt differences (>10% of the sequence length) indicated that the eDNA presence of bighead carp was not likely attributable to error.
Of the four bait shop samples that were positive for silver carp, two contained haplotype "B" alone (0.19-0.78% of reads; S4B Table). One store near central Lake Erie possessed both haplotypes "A/C" (2.72% of reads) and "B" (2.02%). A single shop located near the Wabash River, where silver carp are well-established in the wild [72], contained two previously undiscerned "novel" silver carp haplotypes (labeled "N1", GenBank Accession MK205185, 0.39% of reads and "N2", MK205186, 0.57% of reads), which each differed from haplotype "B" by single nts (S3 Table). No other silver carp sequences were recovered from that shop. Neither silver nor bighead carps were identified with morphological examination of fish samples purchased from the 48 bait stores. NTCs did not amplify and no detections occurred in the positive controls. Thus, it is highly unlikely that eDNA detection was due to contamination or error.

Genetic diversity patterns, founder effect, and the "leading edge" hypothesis
Moderate and consistent levels of genetic diversity (measured by heterozygosity and allelic richness) characterize silver carp populations tested here across much of their invasive North American range, based on both μsats and DNA sequences. Genetic diversity at the two invasion front populations did not significantly differ from longer-established core populations, refuting the leading edge hypothesis. Moreover, the two invasive front populations significantly genetically diverged, according to μsat results.
Our findings indicate that silver carp populations possess genetic diversity levels (A R = 6.1 ±0.07, H O = 0.60±0.14) that appear comparable to those of other North American invasive fish populations (see [6]). For example, the introduced Eurasian ruffe Gymnocephalus cernua (Linnaeus, 1758) in the Great Lakes had much lower allelic richness (A R = 3.09±0.91) and similar heterozygosity (H O = 0.59±0.03), based on comparable numbers of μsat loci; its lower allelic richness is attributable to a likely single source introduction and founder effect [6]. For neutral     Silver carp invasion genetics loci, such as μsats, allelic richness generally decreases more than does heterozygosity with founder effect and bottlenecks, since rare alleles are more readily affected by drift than are common ones [117,118]. In comparison, allelic richness for the round goby's invasion of the Incongruence with mtDNA tree implies a historical mtDNA introgression, possibly with the closely related largescale silver carp. SVC S7-"1" contained 12 individuals with haplotypes "A-H" (two of haplotype "A", one individual each with the "B", "C", "D", "E", "F", and "G" haplotypes, and four having haplotype "H"), S7-"2" (one individual with haplotype "A"), S7-"3" (one individual with haplotype "A"), S7-"4" (three individuals with haplotype "B"), S7-"5" (one individual with haplotype "G"), and S7-"6" (one individual with haplotype "H"). Note that one of the six individuals having the mt "H" haplotype did not amplify for S7, and thus was not included. Posterior probability values presented in decimals at nodes. https://doi.org/10.1371/journal.pone.0203012.g005

Silver carp invasion genetics
Great Lakes (A R = 9.01±0.92) was higher and its heterozygosity levels were similar (H O = 0.60 ±0.05) [5] to those determined here for silver carp, again from a similar number of μsat loci. Population genetic data for silver carp from its native range are scarce, rendering it difficult to directly compare levels of genetic diversity and structure. Farrington et al. 2017 [75] described μsat allelic richness of silver carp in North America at 4.7±0.1 versus 6.1±0.2 in its native range, whose values for the former were lower than those discerned here in our results (mean 6.1±0.7). Their mean heterozygosity values in North America appeared slightly greater (0.65) than ours (0.60). These variations likely reflect the respective loci used and the population samples analyzed by them versus by us. Results of both studies indicate that the silver carp showed slight to no evidence of founder effect or bottlenecks in its North American invasion and expansion [75].
Li et al. 2011 [76] reported 18 unique mtDNA haplotypes (based on concatenated COI and control region (D-loop) sequences) for 94 silver carp individuals from the Illinois and Mississippi rivers, which were undetected in their samples from the native range (N = 121; Yangtze, Pearl, and Amur rivers). Their haplotypic diversity was reported to be lower in the Mississippi River basin (0.74) than in the native range (0.95). Since their sequences are not in GenBank or other publicly accessible data bases, we cannot directly compare our data with theirs. Our haplotypic diversity was 0.61 in the Lower Mississippi River, ranging down to 0.48 in the Missouri River, and averaged 0.54 across all of our samples. Our results and those from other studies suggest little to no founder effect for the introduced silver carp in North America.
Some invasive fish populations have been successful in new areas despite reduced genetic diversity levels, in comparison to native ones. For example, the Indo-Pacific lionfishes Pterois volitans (Linnaeus, 1758) and P. miles (Bennett, 1828), which are widely invasive and prolific in the western Atlantic Ocean and the Caribbean Sea, respectively had nine mtDNA control region haplotypes and just a single haplotype in their introduced ranges, versus 36 and 38 in their native ranges [119]. Sacramento pikeminnow Ptychocheilus grandis (Ayres, 1854) populations have been highly successful invaders of coastal Californian rivers despite 49.6% reduction in allelic richness due to founding by just four individuals [120]. Independent, allopatric introductions of ruffe into the upper Great Lakes and Bassenthwaite Lake, England showed significant founder effects in μsat allelic richness and heterozygosity levels, with low effective population sizes, and declines in genetic diversity over their 30-year time courses [6]. Chen et al. 2012 [121] discerned reduced μsat variation in grass carp introductions on three separate continents in comparison to the native range, with founder effect and bottlenecks were indicated for two of the introductions. The exception was North America, where rapid expansion into the large Mississippi River basin and likely multiple introduction events was hypothesized to have counteracted bottlenecks [121]. Our results suggest that similar rapid expansion may have influenced the genetic diversity of silver carp populations.
Roman and Darling 2007 [12] described that just 37% of introduced populations experienced a significant loss of genetic diversity. For example, most introductions of wakame algae Undaria pinnatifida (Harvey, 1873) contained high genetic diversity worldwide, likely stemming from admixture of multiple native strains [122]. Likewise, the ballast water introductions of dreissenid mussels and round goby in the Great Lakes were large and genetically diverse, involving several likely source populations and introduction events, and showed no founder effects [5,7,[15][16][18][19]. Such ballast water introductions likely involved up to hundreds or thousands of individual propagules being introduced at a time [12,123]. These contrast with the more moderate diversity of the silver carp invasion.
The leading edge hypothesis predicts that species with rapidly expanding ranges would possess lower genetic diversity in peripheral populations [20,124]. This hypothesis has been supported for expansions of several native species from glacial refugia after the Ice Ages [20], including μsat and mtDNA sequence analyses of yellow perch Perca flavescens (Mitchill, 1814) populations [125]. The leading edge hypothesis also was supported at three invasive fronts in Ireland of the bank vole Myodes glareolus (Schreber, 1780), which showed reduced single nucleotide polymorphisms (SNPs), no accumulation of deleterious alleles, and possible selection for traits involved in immunity and behavior [126]. Lower diversity at the leading edge was implicated in North American expansion areas of the vampire bat Desmodus rotundus (É. Geoffroy, 1810) [21].
The leading edge hypothesis has not been supported for some invasions, such as the round goby, whose genetic diversity levels did not differ across the Great Lakes, yet whose populations have remained divergent and respectively stable over their 25 year history [5]. The present study similarly discerns that silver carp does not significantly vary in genetic diversity levels across its invasive range-including at the two fronts tested.
When invasive species spread via jump dispersal, new population areas frequently possess reduced genetic diversity (i.e., founder effect), as was seen with transport of invasive dreissenid mussels from the Great Lakes to smaller lakes to the west (via overland trailered boats) [15,16] and with the European green crab's Carcinus maenas (Kinnaeus, 1758) spread from the Atlantic to the Pacific coasts of North America [127]. This has not been the case for silver carp here, whose populations appear to have slowly and steadily progressed northward, without "leading" edge drop in genetic diversity. Our results indicate continued maintenance of large population sizes and regular gene flow.
The considerable geographic distances separating silver carp full and half siblings in our study support high dispersal, including within a single generation. Tagging studies of silver carp established in the Wabash River showed that some individuals traveled widely (up to 409 km), while others remained close to the tagging location [128]. This might be an adaptive "strategy" for an interplay between dispersal and site fidelity among various offspring, whose genetic and gene expression components [36] merit future investigation.
Our targeted metabarcoding HTS assay could be applied en masse to analyzing eDNA from water or entire plankton samples from net tows (i.e., discerning early life stages of eggs and larvae) to detect new occurrences and dispersal patterns by managers. A targeted assay designed by our laboratory to be specific for the two species of invasive dreissenid mussels, distinguishing between them and many of their haplotypes, has been successfully deployed to simultaneously analyze tens of thousands of veliger larvae to species and population levels from net tows, as well as eDNA water samples [129].

Genetic divergence patterns and population relationships
Significant population genetic differentiation sometimes occurs across the range of invasive species, which often is attributed to high diversity and multiple founding sources, including the round goby [19], whose population patterns in the Great Lakes have remained stable over 25 years [5]. The ruffe exhibited limited genetic diversity and some population divergence across its 30 year history in the upper Great Lakes, revealing founder effect attributable to a single European source area; it has declined in effective population size and numbers in recent years [6]. In comparison, the zebra mussel's genetic composition has changed over its three decades in the Hudson River, experiencing several population turnovers and recolonizations [15,16], yet has remained consistent in Lake Erie throughout [15,129]. Some invasions possessed little population genetic structure across their ranges; for example, colonization by Virginia's warbler Oreothlypis virginiae (Baird, 1860) of the Black Hills in South Dakota showed no differentiation, implicating high gene flow and dispersal [130].
In the present investigation, the silver carp displays modest levels of inter-population divergences across its invasive North American range. Population differences are most pronounced between the "core" area of original introduction versus expansion areas. Overall, the silver carp invasion appears to have progressed steadily, without bottlenecks or marked changes in genetic composition. Most alleles are distributed throughout the populations, with some significant differences between the invasion fronts of the Illinois and Wabash Rivers. The latter shows especially high self-assignment, which merits evaluation of genomic expression.

Genetic history of the invasion, revealed by mtDNA
Interpreting the genetic compositions, diversity levels, and phylogeographic patterns of introduced populations and their sources facilitates understanding how they colonize, expand their ranges (see [118-122, 125-127, 129-131]), and adapt to new environments [2,36]. The silver carp's origins are believed to trace to China, based on aquaculture records summarized by Kolar et al. 2007 [132] and genetic work by Li et al. 2011 [76]. Relation to the latter's results suggest that origins of the two highest frequency haplotypes found in our study ("A" and "B"), which dominated all of our North American populations, likely trace to the Amur and/or Yangtze Rivers [76]. In our results, "A" was more common in the Wabash River invasion front and "B" at the Illinois River invasion front, accounting for the significant divergence between them.
Sequences of cytb haplotype "Q" from the native Yangtze River in China and "R" from an introduction to the Black River in Russia (both from GenBank) are closest to our rare "F" (Upper Mississippi River) and "G" (Lower Misssissippi River) haplotypes. The singleton haplotypes "C", "D", and "E" all occurred in our Lower Mississippi River "core" population. These results suggest that the original population escape of silver carp into the Lower Mississippi River region migrated to found the expansion areas, and that there likely have not been appreciable additional introductions from other overseas sources.

Introgression with largescale silver carp and bighead carp
In comparison to Farrington et al. 2017 [75], who sequenced the entire mitochondrion of 30 silver carp individuals from North America and two from their native range, we sequenced 225 individuals for~10% of some of the most variable portiion of the mt genome, finding eight haplotypes compared to their six, including a much higher frequency of the highly divergent "H" haplotype (for which they uncovered just a single individual) ( Table 5). We included Farringon et al.'s 2017 [75] sequences in our analyes, except for four of their haplotypes, which either were outside of the cytb and COI gene regions that we sequenced, or had incomplete sequences. With our greater sample size, we discerned that the "H" haplotype is found throughout most of the North American range.
"H" occurred in~3% of the individuals, which may represent a historic introgression of the silver carp H. molitrix with the closely related largescale silver carp species, H. harmandi Sauvage, 1884. Our mtDNA sequence for "H" matches a GenBank haplotype (EU315941) from the Yangtze River. The largescale silver carp H. harmandi is morphologically distinguishable in possessing larger scales and a deeper body than does H. molitrix [133]. The largescale silver carp is native to Hainan Island and Vietnam, where silver carp also has been extensively introduced [133]. The two species have widely hybridized in China and Vietnam [110]. Our findings discern that the H. harmandi haplotype "H" is widespread, occurring in the Lower Mississippi, Upper Mississippi, Missouri, and Wabash rivers, as well as in China.
The S7 sequences and μsat allelic frequencies possessed by the "H" haplotype individuals did not show nDNA differentiation. The "H" haplotype thus appears to have resulted from mtDNA introgression, likely of female H. harmandi interbreeding with male H. molitrix in Asia. Introgressed native silver carp likely were present as a small percentage of founders in the North American introduction.
There are 73 nt substitutions separating haplotype "H" from the most similar bighead carp cytb sequence. In comparison, silver carp haplotype "Q", which is closest to "H", diverges by 37 nts. Moreover, bighead carp also is very divergent in its S7 sequence.
Studies have reported that wild silver and bighead carp do not widely interbreed in their native range, having developed pre-zygotic reproductive barriers in sympatry [132,133]. However, the two species often hybridize in aquaculture [77] and hybrids occur across their North American ranges [33,34]. Lamer et al. 2015 [34] employed 57 SNPs and one mtDNA polymorphism to determine that 44.7% of Hypophthalmichthys spp. individuals in the Mississippi River Basin have interbred or backcrossed, suggestive of a hybrid swarm. The flow of maternal DNA was found to be biased from silver carp to bighead carp in both the F1's and bighead carp backcrosses [34]. Greater interspecific recombination in females has been shown by microsatellite genetic linkage mapping of two silver and bighead carp crosses by Guo et al. 2013 [134]. Across all of our samples (which we selected to avoid morphological hybrids; see Methods), we identified just a single individual with bighead carp mtDNA, either due to misidentification or hybridization. Thus, it is not likely that our results were influenced by hybridization between bighead and silver carps.
Given our objective to analyze the population genetic relationships of the silver carp in its invasive range, it is relevant to consider whether and how genetic components from bighead and possibly largescale silver carp may influence them in the future. It is possible that hybridization has acted to enhance invasive success. Our targeted HTS assay can be applied to track these mtDNA genotypes en masse, while genomic techniques can be used to evaluate the relative fitness of variants.

Application of our targeted silver carp and cyprinid HTS assay
Our assay to detect and identify silver carp and its variants, as well as other invasive carps, and native cyprinid species, has potential wide application. Similar metabarcoding eDNA assays that utilized degenerate primers have been shown to be highly sensitive, detecting down to tens of copies of a target marker [110,135]. In a similar but much more limited study of six retailers in the Great Lakes, Mahon et al. 2014 [136] detected white perch Morone americana with metabarcoding in one bait shop in Michigan, in which only a single individual was observed in the tank. Previously published methods to detect invasive carps have relied on quantitative (q)PCR (e.g., [137]), cannot differentiate haplotypes, and therefore do not provide any information on population genetic diversity.
We did not physically find silver carp in any of the bait stores surveyed, based on morphological examination of purchased samples. However, we discovered silver carp eDNA in the tanks of four shops, indicating likely presence at low densities. In addition, eDNA of bighead carp was identified from one of those shops (albeit below the cutoff frequency based on positive controls). We also found eDNA of other invasive cyprinids. Our eDNA findings highlight the probability that silver and bighead carps are prevalent in the retail trade in the Great Lakes watershed, where they likely blend in with other "minnows" for sale. Releases of live bait by anglers may comprise an important vector for their introduction, meriting further investigation. Silver carp in the Mississippi River did not display patterns associated with jump dispersal [15,16]. It is therefore not likely that bait release has played a significant role in its spread in the region to date, but might in its potential spread to the Great Lakes.
The eDNA from the store near the Wabash River contained two silver carp haplotypes that were not recovered in our traditional Sanger sequencing and were not referenced in GenBank (search on November 30, 2018). Although those haplotypes were at low frequency, their possibility of being erroneous would stem from the highly unlikely situation of an undetected index hop co-occurring with two separate sequencing errors that were greater than the probabilities calculated from positive controls. This is unlikely due to use of our custom pipeline that removed index hops and a de-noising algorithm that corrected for sequencing error (see Methods). Both haplotypes have single transitional nt mutations, which are common mtDNA cytb mutations [138]. These haplotypes may have been undetected by traditional sampling and sequencing due to sample size limitations. This shop is located within the invasive range of silver carp. Thus, it is possible that these new eDNA haplotypes are present in the silver carp Wabash River population.
Further refinement of bioinformatic pipelines and use of PCR replicates would enhance the veracity of rare haplotypes discerned with targeted eDNA HTS assays [139,140]. In fact, another targeted HTS assay discerned several novel zebra mussel haplotypes via mass sequencing tens of thousands of veliger larvae, which is beyond the scope of traditional individual population genetic studies [129]. Additional eDNA and plankton sampling with our assay should be undertaken in order to verify these silver carp haplotypes (and others), and their respective distributions.
Our assay has potential widespread use for screening bait fishes, as well as with water and plankton samples from the field. The latter will allow inexpensive and accurate identification of eggs and larvae to species and population [129], which often lack diagnostic morphological characters and cannot be visually discerned. Moreover, surveying water and plankton samples using this approach will allow their relative species proportions to be evaluated, as well as yield population genetic variability statistics for ecological and biogeographic comparisons (see [129]).

Conclusions
The present investigation resolves phylogeographic and genetic diversity patterns across the invasive North American range of silver carp, significantly advancing from prior knowledge using larger sample sizes and a combined population genetic approach. We discern consistent genetic diversity levels across the invasion's expansion, including at the fronts leading to the Great Lakes, which significantly differ in genetic composition. We discover and describe a widespread and highly differentiated ancient mtDNA haplotype, which likely originated from historic introgression with the largescale silver carp in Asia, predating the North American introduction. The introgressed individuals do not differ from silver carp in nDNA μsats or gene sequences, suggesting that female largescale silver carp reproduced with male silver carp.
There is great interest in preventing entry of silver carp into the Great Lakes system using man-made barriers and the use of genetic tools to detect them at likely introduction points, expansions, or via other possible introduction vectors. Our finding of eDNA from silver carp in bait stores in the Lake Erie and Wabash River watersheds points to another likely introduction vector, as these retailers often source from southerly areas where silver carp are wellestablished.
Our targeted HTS assay is designed to discern and detect multiple divergent haplotypes, including the introgressed haplotype "H", as well as distinguish and identify related cyprinid species. It can be used on water samples or with plankton tows, which will yield abundant population genetic information, including at early life stages. At egg and larval dispersal stages, new invasions and expansions may be detected early, when it is more feasible to eradicate or control them. The ability to track changes in species composition and population genetic variation, rather than just single species presence or absence, should significantly enhance our understanding of the successes of invasive carps in colonizing new areas. If silver carp spread into the Great Lakes via natural dispersal, populations likely will possess similar genetic diversity, yet may subtly vary in composition, allowing tracking. The fact that the two invasive fronts analyzed here differ in genetic composition, as well as from the longer-established core area of the lower Mississippi River, shows the influence of genetic drift and possible adaptation. Elucidating these genetic patterns will significantly increase ecological understanding of the relative successes and adaptations of invasions, and aid the rapid identification of new populations and range areas.
Supporting information S1 Table. Additional samples, GenBank Accession numbers, species, population locations, gene regions, and haplotypes (Hap) included in our mtDNA sequence analyses. Table. S7 nuclear ribosomal protein gene, intron 1, partial sequence haplotypes of silver carp. Nucleotide (nt) positions and substitutions, based on alignment to GenBank Accession AY325777 S7-"1" reference sequence, with dots indicating congruence with S7-"1". Dashes denote sequence deletion (indels). Variants are GenBank Accessions MH938813-43. (DOCX) S3 Table. MtDNA sequence alignment of the invasive carp HTS assay region. Alignment shows silver carp (SVC) haplotypes "A-H", and two "novel" haplotypes "N1" and "N2" recovered with the targeted HTS assay (GenBank Accessions: MK205185-6), bighead carp (BHC), and other invasive cyprinid sequences. Nucleotide positions (above sequence) are based on the complete cytochrome b gene, with those differing from silver carp haplotype "A" shown and dots denoting homology. (DOCX) S4 Table. Summary of targeted assay high throughput sequencing run output for targeted invasive carp HTS assay from 48 bait shops in the Lake Erie, Lake St. Clair, and Wabash River watersheds. Raw sequence reads, trimmed reads (had both primers and were the correct length), the number and percent that DADA2 successfully merged for all samples and those having sequences that matched silver carp (% per haplotype, "A/C", "B", "N1", "N2"