Complex introgression among three diverged largemouth bass lineages.

Hybrid zones between diverged lineages offer a unique opportunity to study evolutionary processes related to speciation. Natural and anthropogenic hybridization in the black basses (Micropterus spp.) is well documented, including an extensive intergrade zone between the widespread northern Largemouth Bass (M. salmoides) and the Florida Bass (M. floridanus). Phenotypic surveys have identified an estuarine population of Largemouth Bass (M. salmoides) in the Mobile-Tensaw Delta, with larger relative weight and smaller adult size compared to inland populations, suggesting a potential third lineage of largemouth bass. To determine the evolutionary relationships among these Mobile Delta bass populations, M. salmoides and M. floridanus, putative pure and intergrade populations of all three groups were sampled across the eastern United States. Phylogenetic analyses of 8582 nuclear SNPs derived from genotype-by-sequencing and the ND2 mitochondrial gene determined that Delta bass populations stem from a recently diverged lineage of Largemouth Bass. Using a novel quantitative pipeline, a panel of 73 diagnostic SNPs was developed for the three lineages, evaluated for accuracy, and then used to screen 881 samples from 52 sites for genetic integrity and hybridization on the Agena MassARRAY platform. These results strongly support a redrawing of native ranges for both the intergrade zone and M. floridanus, which has significant implications for current fisheries management. Furthermore, Delta bass ancestry was shown to contribute significantly to the previously described intergrade zone between northern Largemouth Bass and Florida Bass, suggesting a more complex pattern of secondary contact and introgression among these diverged Micropterus lineages.

INTRODUCTION

Hybrid zones are geographical areas of ongoing hybridization between diverged lineages and therefore offer a unique opportunity to study evolutionary processes related to speciation (Barton & Hewitt, 1985; Gompert et al., 2017). The most accepted process by which hybrid zones arise is secondary contact of diverged lineages that have not evolved complete reproductive isolation. This contact may occur due to natural shifts in geological or environmental barriers, or anthropogenic introductions of non‐native species (Largiadèr, 2007; Woodruff, 1973). Nearly all studied hybrid zones are those involving two parental taxa; however, more complex hybrid zones involving multiple taxa are known to occur in nature (Chhatre et al., 2018; Keck & Near, 2009; Nevado et al., 2011; Peñaloza‐Ramírez et al., 2010). When hybridization is thought to be the result of anthropogenic translocations of non‐native species, it can be of particular concern due to loss of biodiversity (McDonald et al., 2008; Rhymer & Simberloff, 1996; Viard et al., 2020).

Interspecific hybridization is relatively common among freshwater fishes. Of over 150 described pairs of hybridizing species in the United States, 20% are found within Centrarchidae, including the black basses (Micropterus spp.; Bolnick, 2009). Natural and anthropogenic hybridization in black basses is well documented, using either morphometrics (Bailey & Hubbs, 1949; Baker et al., 2013) or diagnostic genetic markers (Bagley et al., 2011; Bangs et al., 2018; Barthel et al., 2010; Li et al., 2015; Lutz‐Carrillo et al., 2006; Philipp et al., 1983; Thongda et al., 2020). One of the most notable cases of hybridization in Micropterus is an extensive intergrade zone between the widespread northern Largemouth Bass (M. salmoides) and the Florida Bass (M. floridanus), of which pure populations are thought to be restricted to peninsular Florida. This hybrid zone was initially described using morphology and meristics and was thought to extend through northern Florida, Georgia, and small parts of Alabama and South Carolina (Bailey & Hubbs, 1949). Thirty years later, Philipp et al. (1983) evaluated allele frequencies at two diagnostic allozymes which extended the hybrid zone west to Mississippi and north to Virginia. Their result is unsurprising, for disconnects between morphological and genetic descriptions of hybrid zones are common, as phenotype can often be inaccurate for identifying later‐generation hybrids (Lutz‐Carrillo et al., 2018; Meyer et al., 2017; Taylor et al., 2018). While Philipp et al. (1983) acknowledged the inaccuracy of meristics for detecting hybrids, they also suggested that the apparent hybrid zone expansion was due to stockings of Florida Bass outside their native range for recreational fishing. This latter view has persisted, with a recent review portraying the geographic extent of natural hybridization the same as was first described by Bailey and Hubbs (Taylor et al., 2019). Modern genetic tools are required to accurately characterize the hybrid zone of M. salmoides–M. floridanus and determine whether this zone is primarily shaped by stocking or historical hybridization.

In addition to the high diversity of Micropterus species, the southeastern United States has the greatest aquatic diversity in North America (Warren et al., 1997), which is exemplified by the numerous freshwater and brackish water fish species supported by the Mobile‐Tensaw River Delta in Alabama (Swift et al., 1986; Swingle & Bland, 1974; Swingle et al., 1966). With the confluence of the Tombigbee and Alabama rivers at its northernmost point and Mobile Bay to the south, the Mobile Delta is the second largest delta in the contiguous United States. Due to a rich history of geologic shifts, glaciation, and sea‐level fluctuations, freshwater fishes of the southeastern United States have experienced repeated bouts of habitat isolation and reconnection (Swift et al., 1986). These dramatic changes in recent evolutionary history are likely responsible for the region's high endemicity and hybridization observed in numerous freshwater fish, such as mosquitofishes (Wilk & Horth, 2016), crappies (Travnichek et al., 1996), sunfish (Avise & Saunders, 1984), and bass (Near et al., 2003). Such sea‐level fluctuations in the Late Pliocene are thought to have isolated M. floridanus from M. salmoides on the Florida peninsula (Near & Kim, 2021), resulting in genetically and phenotypically diverged sister taxa.

In the Mobile Delta exists another phenotypically distinct population of largemouth bass. These “Delta” bass have a smaller adult size and higher relative weight compared to other M. salmoides populations, as well as physiological differences that may facilitate their survival in water with elevated salinity (up to 13 ppt; DeVries et al., 2015; Glover et al., 2012, 2013; Tucker, 1985). Preliminary genetic analysis using isozymes and microsatellites showed that Delta bass were genetically similar to M. salmoides, but were unable to conclusively determine whether they were a distinct genetic lineage (DeVries et al., 2015; Hallerman et al., 1986). Given the proximity of the Mobile Delta to the M. salmoides–M. floridanus intergrade zone, a thorough genetic characterization of Delta bass is necessary for determining how this putative lineage contributes to existing hybridization patterns.

The objectives of this study were to employ modern genomic techniques to analyze the M. salmoides–M. floridanus hybrid zone and resolve the phylogenetic relationship of Largemouth Bass from the Mobile‐Tensaw Delta. Specifically, we (1) used thousands of single nucleotide polymorphisms (SNPs) and mitochondrial sequencing to characterize the phylogenetic and population genetic relationships between Delta Largemouth Bass, northern Largemouth Bass, and Florida Bass; (2) developed a SNP assay for accurate and rapid identification of pure and hybrid individuals; and (3) applied this SNP assay to hundreds of samples across the eastern United States in order to characterize the geographic extent of hybridization and assess the role of stocking in driving the extent of Florida Bass introgression.

METHODS

Sample collection and GBS library preparation

Sample collection and genotype‐by‐sequencing (GBS) using the PstI restriction enzyme was previously described in Thongda et al. (2020). For the present study, we used GBS data for 144 largemouth bass samples across 18 sites and captive hatchery populations: Largemouth Bass from the Mobile‐Tensaw Delta (DLB; M. salmoides, N = 29), northern Largemouth Bass (NLB; M. salmoides; N = 42), Florida Bass (FLB; M. floridanus; N = 29), and putative intergrade largemouth bass (ILB; N = 44; Table S1). For phylogenetic analysis, we also included 17 outgroup black bass samples: Spotted Bass (M. punctulatus; n = 2), Coosa Redeye Bass (M. coosae; 5), Shoal Bass (M. cataractae; n = 5), Smallmouth Bass (M. dolomieu; 3), and Guadalupe Bass (M. treculii; n = 2). FLB, NLB, and ILB individuals were identified using previously developed 25–38 SNP panels that are diagnostic for northern Largemouth and Florida Bass (Li et al., 2015; Zhao et al., 2018). We also sampled an additional 786 largemouth bass individuals across 52 populations for extensive population genetic structure analysis and hybrid classification using the diagnostic SNP panel developed in this study (Figure 1; also refer to File S2). For the 881 total samples genotyped with the diagnostic SNP panel (including some of the samples used for GBS), DNA was extracted from fin clips using a simple sodium hydroxide (NaOH) and hydrochloric acid (HCl) tissue digestion (Truett et al., 2000). For mitochondrial gene sequencing, DNA was extracted from fin clips or whole blood samples using E.Z.N.A. DNA Kits (Omega Bio‐Tek).

FIGURE 1

Sampling sites of largemouth bass assayed with 73 SNP panel. Populations are labeled as in Table S2 and represented by pie graphs showing the mean estimated ancestry proportions for three bass lineages based on 73 diagnostic SNPs and the population genetic program STRUCTURE (K = 3). Delta Largemouth Bass (DLB) are pink, northern Largemouth Bass (NLB) are blue, and Florida Bass (FLB) are green. Native ranges of NLB and FLB are shown, as well as the NLB‐FLB hybrid zone as described by Bailey and Hubbs (1949). Map data from ESRI; FLB and NLB native ranges from Taylor et al. (2019)

GBS analysis and SNP calling

Raw GBS reads were demultiplexed, quality filtered, and aligned to a M. floridanus genome reference (N50 of 11,136 bp and 249,768 scaffolds (Zhao et al., 2018)) using STACKS v2.55 (Rochette et al., 2019). GBS reads were trimmed and demultiplexed using the process_radtags program in STACKS, while removing reads with an uncalled base (‐c) and low‐quality scores (‐q). For reference‐based SNP calling, we mapped the demultiplexed reads to the M. floridanus genome using the BWA‐MEM algorithm from BWA (Burrows–Wheeler Aligner) v0.7.17 with default settings (Li & Durbin, 2009). The mapped reads were sorted using the “sort” function in SAMTOOLS v1.6 (Li et al., 2009) and then used to call SNPs with the “gstacks” pipeline and default SNP model (“marukilow”), which uses a Bayesian approach to identify SNPs across all samples for each GBS locus and then genotypes each individual at each SNP. The alpha thresholds for SNP discovery and calling genotypes were left at the defaults of 0.05, and the maximum soft‐clipping level was increased to 0.3 (–max‐clipped 0.3). The “populations” program was used to merge loci that were produced from the same restriction enzyme cut sites (–merge‐sites), perform initial SNP filtering, and export a variant call format (VCF) file. Stringent additional SNP filtering was done with VCFtools v0.1.17 (Danecek et al., 2011), custom Python code, and code adapted from Jon Puritz's laboratory (Puritz et al., 2014) to remove potentially erroneous genotypes due to low coverage, sequencing error, or paralogs.

For population genetic analyses, biallelic SNPs were retained if they had a minor allele frequency greater than 0.05, had a minor allele count greater than 4, were genotyped in at least three population groups (NLB, DLB, FLB, or ILB), and had <15% missing data across all individuals. SNPs that departed Hardy–Weinberg equilibrium with a p‐value cutoff of 0.05 in two or more sampling sites were also removed. SNPs in linkage disequilibrium (r 2 > 0.2) or within the same 1000‐bp window were pruned using BCFtools +prune plugin (Li et al., 2009). For phylogenetic analyses that included outgroup Micropterus species, variable sites were retained if they were genotyped in at least 4 species and 70% of all individuals, were heterozygous in fewer than 70% of samples (–max‐obs‐het), and had a minor allele count greater than 3. The VCF output of stacks was converted to a phylip file using vcf2phylip v2.0 (Ortiz, 2019).

ND2 sequencing

DNA from 58 samples was amplified with previously published primers targeting the mitochondrial NADH subunit 2 (ND2) gene (Kocher et al., 1995). Amplification of the ND2 gene used 25 ml reaction volumes and included 1× Kapa HiFi HotStart ReadyMix (2×; Kapa Biosystems), 0.2 mM forward primer, 0.2 mM reverse primer, and 2 μl of template DNA. All PCRs included negative control reactions (no DNA template). Thermal cycler reactions were conducted with a heated lid and the following conditions: initial denaturation at 95°C for 3 min, followed by 5 cycles at 94°C for 30 s, annealing at 57°C for 30 s, and 72°C for 75 s; an additional 5 cycles of annealing at 56°C for 30 s and 5 cycles with annealing at 55°C for 30 s; followed by 20 cycles at 94°C for 30 s, annealing at 54°C for 30 s and 72°C for 75 s, and a final extension at 72°C for 10 min, followed by a final hold at 4°C (Baker et al., 2013). PCR products were visualized on a 1% agarose gel to check the amplification success and that a single band was observed. Amplicons were then cleaned using ExoSAP‐IT (Thermo Fisher Scientific). Quantity and quality of PCR products were determined with a NanoDrop Spectrophotometer (Thermo Fisher Scientific).

Sanger cycle sequencing of all products was performed by Genewiz, Inc. using the PCR primers and internal ND2 primer MET‐F (Kocher et al., 1995). Consensus sequences of ND2 were obtained by analyzing the forward sequencing trace files using CodonCode Aligner (2019.2.1) and then trimmed to the first 598 bp. The reads were inspected manually for quality and for any discrepancies in base calls. All sequences were aligned in MEGA‐X using MUSCLE to determine the number of unique haplotypes present for subsequent analyses (Edgar, 2004; Kumar et al., 2018). A minimum spanning network for haplotypes was calculated and visualized in the R package pegas v0.14 (Paradis, 2010).

Phylogenetics and population genetics

Maximum likelihood phylogenetic analyses of 100 individuals from pure NLB, FLB, and DLB populations and 17 outgroup Micropterus samples were conducted using 18,492 concatenated variable GBS sites with IQ‐TREE v1.6.12 (Nguyen et al., 2015). ModelFinder as implemented in IQ‐TREE used BIC and determined the best‐fit substitution model that included ascertainment bias correction to be TVM + F + ASC + R3 (Kalyaanamoorthy et al., 2017). Branch supports were assessed with both 1000 ultrafast bootstrap approximation replicates and 1000 bootstrap replicates for the SH‐like approximate likelihood ratio test (Guindon et al., 2010; Hoang et al., 2018). The resulting tree was rooted with the five outgroup bass species. The same analysis was performed on 18,641 SNPs with ILB individuals included. The resulting best trees were plotted with the R package ggtree (Yu et al., 2017).

Population genetic summary statistics were calculated on GBS SNPs for population groups (NLB, DLB, FLB, ILB) and individual sampling sites with at least 5 individuals. Based on equations from Nei and Chesser (1983), observed heterozygosity (Ho), expected heterozygosity (He), overall FST, and FIS were calculated using the basic.stats function in the R package hierfstat v0.5.7 (Goudet & Jombart, 2015). Confidence intervals for population‐specific FIS were determined using the boot.ppfis function in hierfstat with 1000 bootstrap replicates. Pairwise FST following Weir and Cockerham (1984) was calculated using the genet.dist function in hierfstat (Weir & Cockerham, 1984). Variable sites were determined to be “fixed” between groups if one group had an allele frequency ≥0.98 in one group and <0.02 in the other (based on allele frequencies reported by STACKS). This cutoff was chosen to allow for singleton alleles observed in a lineage that may be due to genotyping error.

The model‐based Bayesian clustering method STRUCTURE v2.2.4 (Pickrell & Pritchard, 2012) was used to determine the number of distinct genetic clusters (K), using the admixture model with correlated allele frequencies, no prior location information, and a burn‐in period of 50,000 repetitions followed by 200,000 repetitions. Five replicate analyses were performed on the thinned GBS SNP dataset with values of K = 1–10. Replicates were summarized and visualized using the CLUMPAK server (Kopelman et al., 2015). The ∆K method implemented in STRUCTURE HARVESTER was used to determine an optimal K (Earl & vonHoldt, 2012). Principal component analysis (PCA) was implemented in the R package adegenet v2.1.1 (Jombart & Ahmed, 2011) and visualized with the R package PCAviz v0.3.29 (Novembre et al., 2018). Missing data were filled by randomly drawing an allele based on the overall allele frequency across all individuals in one of the prespecified groups (NLB, DLB, FLB, ILB) using custom R code. As other methods support significant differentiation among these groups, this approach is more appropriate than filling in missing data using the mean or overall allele frequency.

Panel development, testing, and genotyping

Our goal when developing a SNP assay for accurate identification of largemouth bass lineages and hybrids was to find a subset of SNPs that best matched the ancestry results determined by STRUCTURE when using the full GBS SNP dataset. We first used custom Python scripts and the populations.sumstats.tsv output from Stacks to identify 2809 SNPs that were diagnostic between DLB, NLB, and FLB, with the criteria that they had >90% or <10% allele frequency in at least two of the three largemouth bass lineages. PCA was then performed on these diagnostic SNPs for 100 “pure” largemouth bass samples (excluding intergrades) as previously described. To identify SNPs that had the greatest contributions to the PCA and therefore were the most diagnostic, we extracted 549 SNPs in the 90% quantile of loadings for PC1 and PC2. We then randomly sampled 80–100 SNPs from this set for analysis with STRUCTURE and compared the Q‐values for each individual with the Q‐values derived using the full SNP dataset. This analysis was iteratively repeated 300 times. A subset of 127 total SNPs with the fewest samples exhibiting a > 5% difference in Q‐values was retained for marker development.

A MassARRAY System (Agena Bioscience) was used to develop, genotype, and evaluate SNP panels. For each of the 127 diagnostic SNPs, a 201‐bp sequence was extracted from the M. floridanus genome (the SNP and 100 bp flanking on either side) and imputed into the MassARRAY Assay Design Software to design two multiplex assays with a maximum of 60 SNPs per assay. This produced forward, reverse, and extension primer sequences for assays of 52 and 39 SNPs (File S3), which were ordered through IDT (Integrated DNA Technologies, Inc.). We used 91 samples to test the concordance of SNP genotypes generated with the MassARRAY platform and GBS, as well as 70 technical replicates made up of 40 individuals to assess the consistency of genotype calls between runs. Discordant genotypes due to missing data were excluded from this analysis. SNPs that consistently failed in greater than 90% of samples were invariant, or inconsistent between runs were removed, resulting in a final set of 73 SNPs. For all samples run on the MassARRAY platform, amplification and extension reactions were performed using the iPLEX Gold Reagent Kit (Agena Bioscience) according to the manufacturer's protocol (Gabriel et al., 2009). SNP genotypes were called using the MassARRAY Typer 4 analysis software and manually confirmed. This software uses a three‐parameter (mass, peak height, and signal‐to‐noise ratio) model to estimate genotype probabilities. Genotype concordance and SNP panel evaluation was performed with custom Python scripts and Excel.

Some samples used for GBS, as well as 786 additional samples across 52 sites, were genotyped at these 73 SNPs for a total of 881 samples. To determine the FLB, NLB, and DLB ancestry proportions of these samples, STRUCTURE was run using 62 reference individuals and the USEPOPINFO, POPFLAG, and PRCOMP settings. Individuals which had a Q‐value >0.94 from STRUCTURE when using the full GBS SNP dataset were chosen to be references (14–24 per lineage). Other parameters: admixture model, K = 3, correlated allele frequencies, migration prior 0.05, and burn‐in of 20,000 followed by 150,000 MCMC iterations. A subset of 751 samples were also genotyped using a panel of 35 SNPs that had previously been developed for differentiating FLB and NLB, where ancestry proportions were determined by counting the number of FLB or NLB alleles using a custom R script (Li et al., 2015; Zhao et al., 2018).

Genotype simulation and assignment

To evaluate the accuracy of our ancestry assignment using the 73 SNP panel and STRUCTURE, we employed a simulation approach developed by Vähä and Primmer (2006). Individuals used as the pure reference populations in our STRUCTURE runs (14 FLM, 24 DLB, 24 NLB) were used to simulate genotypes from random mating and hybridization with the hybridize R function in adegenet. One hundred genotypes were simulated for each pure lineage, 50 for each type of F1, 150 for each type of F2 through F4 generations assuming neutral admixture, and 300 triple hybrids (considered the product of crossing an F1 from a pair of species with a pure individual of the third species or a backcross between a triple hybrid and a pure species). These simulated genotypes were analyzed using STRUCTURE with the same settings and reference genotypes. The performance of STRUCTURE for hybrid and pure individuals was evaluated based on efficiency (number of individuals in a simulated group that are correctly assigned), accuracy (proportion of an identified group that truly belongs to that group), and performance (efficiency multiplied by accuracy). Finally, the optimal threshold values of Q to assign individuals to the different genotypic categories were determined.

Geographical patterns of Florida Bass introgression

To determine whether patterns of FLB genetic ancestry could be explained by both geographical variables and stocking history, we first calculated the great circle distance in km between our reference FLB population (St. Johns River, FL) and all other sites sampled with our SNP panel using the R package geosphere v1.5‐10 (Hijmans et al., 2016). We then modeled the proportion of FLB ancestry for all 881 individuals as the interaction between distance from St. Johns River, FL, and categorical stocking history, using a fractional logistic regression model implemented in base R stats [glm(FLB ancestry ~ distance * stocking, family = “quasibinomial”)]. Based on published literature and discussions with biologists at state agencies, sites were designated as either “River” or “Reservoir,” where reservoirs were natural or artificial bodies of water with known stocking history (Alford & Jackson, 2009; Bunch et al., 2017; Hargrove et al., 2019). Reservoirs or lakes that had no known history of largemouth bass stocking (e.g., Lake Mattoon, IL) were coded as “River.” Ancestry proportions were those calculated using the 73 SNP diagnostic panel described below. This model was compared with a model based only on distance [glm(FLB ancestry ~ distance, family = “quasibinomial”)], by comparing the reduction in deviance with a chi‐squared test.

RESULTS

GBS assembly

A total of 42,390 biallelic SNPs across 102,685 GBS loci were genotyped in greater than 50% of 144 largemouth bass individuals and at least three of the four groups (NLB, ILB, FLB, DLB). Mean length of GBS loci was 91 bp (SD = 0.02). Further filtering by missing data <15%, HWE, MAF >5%, and thinning by LD reduced the dataset to 8582 SNPs. Average read depth across SNPs per individual ranged from 3.3 to 29.2 (mean = 11.2 ± 4.6).

Population genetics

A PCA of 144 individuals using 8582 SNPs clearly differentiated the three lineages across PCs 1 and 2, with PC1 representing 30% of SNP variation and PC2 representing 11% variation (Figure 2). We calculated population genetic summary statistics on these 8582 SNPs, with samples separated by either lineage group or sampling site. As expected, intergrade samples had the highest observed genome‐wide heterozygosity (0.31), followed by DLB (0.16). NLB had the highest levels of inbreeding (0.105; Table 1). Overall FST when separated by lineage group was 0.426. Pairwise FST between lineages was highest between FLB–NLB (0.749) and lowest between DLB–NLB (0.459). Fixed SNPs were detected between each pair of lineages, with the fewest between DLB–NLB (46 unlinked SNPs) and the most between FLB–NLB (1807; Table 2). There were 532 unlinked SNPs that were fixed between FLB–NLB but not between DLB–FLB, suggesting possible introgression between DLB–FLB. When comparing summary stats between sampling sites, Sugar Lake, MN, had the lowest observed heterozygosity, followed by FLB samples from a captive hatchery population (Table S1). Pairwise FST between sites within a lineage group was generally low (mean pairwise FST 0.051–0.124), with the highest within‐lineage FST observed between Sugar Lake, MN, and Hatchery NLB (0.239; Table 1).

FIGURE 2

Principal component analysis plot of PC1 and PC2 for 144 Micropterus individuals based on 8582 GBS SNPs. Individuals are colored by lineage (DLB = Delta Largemouth Bass, NLB = northern Largemouth Bass, FLB = Florida Bass, ILB = intergrade largemouth bass)

TABLE 1

Population genetics of largemouth bass lineages

Group	Ho	He	Avg. pairwise FST (min–max)	FIS
ILB	0.31	0.32	0.051 (0.0134–0.083)	0.032 (0.027–0.036)
DLB	0.16	0.17	0.031 (0.008–0.05)	0.070 (0.063–0.077)
FLB	0.11	0.13	0.09 (0.065–0.11)	0.079 (0.069–0.089)
NLB	0.12	0.14	0.124 (0.037–0.239)	0.105 (0.097–0.112)

Population genetic summary statistics based on 8582 SNPs for 144 individuals across the three largemouth bass lineages and their intergrades (DBS = Delta Largemouth bass, FLB = Florida Bass, NLB = northern Largemouth bass, ILB = intergrades). Ho: observed heterozygosity averaged across loci; He: expected heterozygosity averaged across loci; FIS: Wright's F‐statistics averaged across loci (Nei & Chesser, 1983).

TABLE 2

Genetic differentiation of largemouth bass lineages

	ILB	DLB	FLB	NLB
DLB	0.203	0	1275 (15%)	46 (0.5%)
FLB	0.329	0.677	0	1807 (21%)
NLB	0.347	0.459	0.749	0

Lower triangle: pairwise F ST between lineages of largemouth bass and intergrades using 8582 GBS SNPs. Upper triangle: number and percentage of unlinked, “fixed” (>0.98 allele frequency) SNPs between lineages.

Phylogeny and inferred ancestry

Maximum likelihood phylogenetic inference using 18,492 GBS SNPs supported DLB as a distinct lineage that is sister to NLB, with 100% support from both the SH‐like approximate likelihood ratio test (SH‐aLRT) and ultrafast bootstrapping (Figure 3). When samples from intergrade populations were included, SH‐aLRT support for DLB as a unique lineage was 99.9% and ultrafast bootstrap support was 96% (Figure S1). For the first 598 bp of the ND2 mitochondrial gene, there were two Delta haplotype, ten FLB haplotypes, and ten NLB haplotypes (Figure 4). There were two fixed ND2 SNPs between DLB‐NLB, 21 between DLB‐FLB, and 19 between FLB‐NLB. Intergrade samples, which had STRUCTURE‐inferred ancestry proportions of at least 6% in two or more lineages, were distributed across both the GBS phylogenetic tree and across mitochondrial haplotypes (Figure S1 and Figure 4). Individuals from some sampling sites, including Sugar Lake, MN, and Demopolis Lake, AL, sequenced well on the MassARRAY platform but did not amplify or sequence well with available ND2 primers (Kocher et al., 1995), suggesting a polymorphism in the ND2 primer recognition sites. Some intergrade individuals appeared to be heterozygous at ND2 sites that were diagnostic between lineages, which may be a sign of heteroplasmy (Piganeau et al., 2004) or sequencing of homologous regions in the nuclear genome (nuMTs; Simone et al., 2011). This included two individuals from West Point Reservoir, GA, one from Lake Eufaula, AL, one from Hatchie River, TN, and one from Sutton Lake, NC. All of these sites have known recent stocking history.

FIGURE 3

Maximum likelihood phylogeny constructed using 18,538 concatenated GBS SNPs and IQ‐TREE (1000 ultrafast bootstrap replicates). Nodes in red indicate 100% ultrafast bootstrap support; nodes in yellow indicate 50%–99% ultrafast bootstrap support. STRUCTURE plots (K = 3) are included and show the ancestry membership proportions (Q‐values) for the three largemouth bass lineages as inferred using 8582 SNPs (DLB = Delta Largemouth Bass, NLB = northern Largemouth Bass, FLB = Florida Bass)

FIGURE 4

(a) Minimum spanning network of mitochondrial haplogroups from 58 individuals, based on the first 529 bp of the ND2 gene. Circle size is proportional to haplotype frequency and colors represent sampling location. (b) STRUCTURE (K = 3) results based on 73 diagnostic SNPs, ordered by mitochondrial haplotype group (DLB = Delta Largemouth Bass, FLB = Florida Bass, NLB = northern Largemouth Bass)

Panel development

From the 127 diagnostic SNPs retained for panel development, the MassARRAY Assay Design software developed primers for 91 SNPs across two multiplexed panels (52 and 39). Detailed information on SNP panels, including SNP ID, genome position, and primer sequences (forward, reverse, and extension), are listed in File S3. After removing SNPs that sequenced poorly or were invariant, 73 SNPs were retained and run on 881 individuals. STRUCTURE results on the 144 GBS individuals using 8582 GBS SNPs, 2809 diagnostic GBS SNPs, and the 73 panel SNPs were highly concordant (Figure 5). We also examined the consistency of MassARRAY genotype calls among 70 technical replicates and found 98.4% of genotypes matched across multiple runs. Analysis of simulated genotypes indicated that a threshold value in the admixture coefficient of Q ≥ 0.94 resulted in adequate separation between pure species and hybrid individuals (Table 3). The genetic assignment of simulated pure individuals from the three lineages showed high levels of performance (98%–100%). For hybrid classes involving two lineages and fewer than four generations of backcrossing, performance values ranged from 99% to 93%. While individuals were classified as NLB–FLB hybrids with 100% accuracy, 6.5% of the simulated NLB–FLB hybrids were misclassified as triple hybrids, showing a small observed DLB Q‐value (<0.1). Triple hybrids had >90% performance. When fourth‐generation backcrossed hybrids (Bx4) were included, efficiency decreased to 82.8%–92.4% for hybrids and accuracy decreased to 58%–78.7% for pure species, as many Bx4 hybrids were assigned as pure (Table S2).

FIGURE 5

STRUCTURE results (K = 3) from three GBS SNP subsets (K = 3). (a) 8582 SNPs, (b) 2809 diagnostic SNPs, and (c) 73 SNPs chosen for SNP panel. NLB = northern Largemouth Bass, DLB = Delta Largemouth Bass, FLB = Florida Bass, ILB = intergrade largemouth bass

TABLE 3

Evaluation of 73 SNP panel on simulated genotypes

Simulated/assigned	DLB	FLB	NLB	DLB‐FLB	DLB‐NLB	NLB‐FLB	Triple hybrid	Total simulated
DLB	100							100
FLB		100						100
NLB			100					100
DLB‐FLB F1				50				50
DLB‐FLB F2:F4				149			1	150
DLB‐FLB BX	2			198				200
DLB‐NLB F1					50			50
DLB‐NLB F2:F4					150			150
DLB‐NLB BX			1		199			200
NLB‐FLB F1						47	3	50
NBB‐FLB F2:F4						131	19	150
NLB‐FLB BX				1	2	193	4	200
Triple Hybrid				4			296	300
Total assigned	102	100	101	402	401	371	323	1800
Efficiency	100.00%	100.00%	100.00%	99.25%	99.75%	92.75%	98.67%
Accuracy	98.04%	100.00%	99.01%	98.76%	99.50%	100.00%	91.64%
Performance	98.04%	100.00%	99.01%	98.02%	99.25%	92.75%	90.42%

Number of simulated individuals (rows), which were assigned to one of three lineages or into a hybrid category (columns). Simulated individuals were assigned based on their STRUCTURE Q‐values, with a threshold of Q ≥ 0.94 for pure individuals. Parameters of efficiency, accuracy, and overall performance of the assignment method are given in percent. Individuals correctly assigned are in bold type. DLB: Delta Largemouth Bass (M. salmoides), NLB: northern Largemouth Bass (M. salmoides), FLB: Florida Bass (M. floridanus); F1–4: first‐fourth generation hybrids, BX: backcrosses.

Geographical patterns of ancestry and hybridization

Population structure analysis of 52 sites using the 73 SNP panel identified populations with a majority of pure DLB individuals extending north of the Mobile Delta, encompassing all tested populations below the Fall Line and outside of peninsular Florida, including those in remote areas with no history of stocking (e.g., Sipsey River; Figure 1, #41). Pure FLB individuals were found far outside of the currently accepted Florida Bass native range, extending into Georgia south of the Fall Line and up the Atlantic coastal plain into North Carolina, while intergrade populations between FLB and NLB were also found outside of the intergrade zone as described by Bailey and Hubbs (1949). Twenty‐six intergrade populations with >6% mean ancestry from all three lineages were identified, primarily in Alabama, Georgia, Tennessee, and Louisiana. Compared to previous estimates of FLB–NLB ancestry using a 35 SNP panel, this 73 SNP panel appears to more accurately represent FLB ancestry in sites with high DLB ancestry. For example, while the 35 SNP panel found that populations from the Mobile Delta were 27%–30% FLB and 70%–73% NLB, the 73 SNP panel identifies them as 95%–98% DLB (Figure S3, File S2).

To evaluate the relative roles of geography and stocking history for determining FLB genetic ancestry, we implemented fractional logistic regression models based on the distance of a site from our reference FLB population and the stocking history of a site (coded as either River or Reservoir). The proportion of FLB ancestry in an individual was best described by a model that included stocking: FLB ancestry ~ distance * stocking (chi‐squared test of reduction in deviance, p‐value = 7.43e‐12), where separate geographic clines were fit for River and Reservoir sites (Figure 6). Based on values predicted by the model for Reservoir sites, stocking has increased the overall geographic extent of FLB intergrades (FLB ancestry >0.06) by 436 km on average. However, the model also indicates that individuals with pure FLB ancestry observed in Georgia below the Fall Line are likely not the result of anthropogenic movement. Both the model and our observations in unstocked River sites support a larger natural distribution of pure FLBs and intergrades than is currently accepted (Taylor et al., 2019).

FIGURE 6

Scatterplot of M. floridanus (FLB) ancestry proportions for 881 individuals, plotted by great circle distance from a reference pure FLB site in St. Johns River, FL. FLB ancestry proportions are based on STRUCTURE results (K = 3) using a panel of 73 diagnostic SNPs. Individuals are colored by stocking history (River = no significant FLB stocking, Reservoir = known stocking history). Lines show the fit of a fractional logistic regression model (FLB ~ distance * stocking)

DISCUSSION

Using a combination of 8582 SNPs and targeted mitochondrial gene sequencing, we identified a complex hybrid zone involving three distinct lineages of largemouth bass (M. floridanus and two from M. salmoides), including the first thorough genetic characterization of a unique lineage of M. salmoides radiating out from the Mobile‐Tensaw Delta. To assess the geographic extent of hybridization, we developed a suite of 73 SNPs that can accurately classify advanced‐generation hybrids between all three lineages. This SNP panel will be a valuable tool for fisheries management and conservation of these economically and ecologically important species, while also facilitating future investigations of introgression and species boundaries. While stocking of Florida Bass in reservoirs for recreational fishing has contributed in part to observed hybridization patterns, the natural geographic extent of both pure M. floridanus and introgressed populations are nonetheless much larger than currently accepted.

Distinct mobile delta lineage of largemouth bass

Support from both nuclear and mitochondrial data indicates that Largemouth Bass from the Mobile‐Tensaw Delta (DLB) are an evolutionarily diverged lineage of M. salmoides. Genetic differentiation between DLB and northern Largemouth Bass (NLB) was high (FST = 0.459), with 46 unlinked nuclear SNPs fixed between the two groups and two fixed SNPs in the first 598 bp of the mitochondrial ND2 gene. While significant, this differentiation is still much lower than that observed between Florida Bass (FLB) and NLB (FST = 0.749, 1807 fixed unlinked nuclear SNPs, 19 fixed ND2 SNPs). These values are higher than those observed using GBS SNPs for different species of Squalius chubs (FST = 0.126–0.414) and comparable to FST between diverged lineages of walleye (0.805; Mendes et al., 2021; Zhao et al., 2020). Our results contribute to a growing scientific consensus around the validity of M. floridanus as a distinct species (Kassler et al., 2002; Near et al., 2003; Near & Kim, 2021). While this study cannot address whether Delta bass should be considered a new species or subspecies of M. salmoides, it does indicate that morphological evaluations aided by genetic assays are warranted.

Maximum likelihood phylogenetic analysis of nuclear loci and a haplotype network analysis of a mitochondrial locus suggest that DLB–NLB diverged more recently than NLB–FLB. Fossil‐calibrated phylogenies of Micropterus using either two mitochondrial genes or 16 nuclear genes dated the divergence of NLB–FLB at the end of the Pliocene (~3.3 Mya; Near et al., 2003; Near & Kim, 2021), suggesting that a NLB population in the Mobile Delta may have diverged later during the Pleistocene, when glaciation and sea‐level changes resulted in fragmented coastal freshwater habitat (Bermingham & Avise, 1986). Once the population was reconnected with the mainland, we hypothesize these Delta Largemouth Bass migrated out of the Delta and met with inland M. salmoides and M. floridanus. The Mobile Delta region is also home to genetically diverged populations of walleye and yellow perch, although the estimated divergence time between these relict populations and their northern counterparts range from 1.7 Mya for walleye to 170,000 ya for perch (Stepien et al., 2015), suggesting that the timing of population divergence in this geographic locale can vary across fish taxa.

Genome‐wide genetic diversity was higher in pure DLB populations compared to pure populations of FLB or NLB (Table 1 and Table S1). As we only had two to four wild populations of each lineage, this may reflect our limited geographical sampling. This result contrasts with walleye and perch diversity, where their Mobile Delta populations had lower genetic diversity than those in the north (Stepien et al., 2015). Elevated genetic diversity is encouraging, as it suggests that the DLB populations have not experienced dramatic population declines and bottlenecks (Nei et al., 1975). Historical or recent introgression with FLB is one potential explanation for higher genetic diversity in DLB populations. Genetic divergence between DLB–FLB is lower than between NLB–FLB (Table 2), supporting this hypothesis. Introgression between the two may have occurred naturally across the Florida panhandle, or due to unsuccessful stocking efforts that occurred in the Mobile Delta sporadically from 1988 to 2000 (Armstrong et al., 2000; Rainer, 2010). Sugar Lake, MN, had the lowest genetic diversity out of all sampled sites, which may be due to a recent genetic bottleneck and range expansion after the last glacial maximum (Dyke & Prest, 1987). Sugar Lake is also the only NLB population with >99% mean inferred NLB ancestry based on STRUCTURE. We hypothesize that the low levels of DLB ancestry (<5%) observed in other “pure” NLB populations, such as those in Illinois, are likely due to ancestral polymorphisms and incomplete lineage sorting. The potential genetic bottleneck experienced by Sugar Lake may have removed some of these shared variants between DLB and NLB from the population.

Multilineage hybrid zone

The hybrid zone between NLB and FLB was first described by Bailey and Hubbs in 1949 using meristics and found to extend only through Georgia and a small part of Alabama (Bailey & Hubbs, 1949). Although subsequent studies using larger sample sizes and biochemical or genetic assays have determined the geographic extent of hybridization to be much larger, these results are often discounted as merely the product of M. floridanus recreational stocking outside of Florida. Using our diagnostic 73 SNP panel, we assayed 52 populations across the east and southeast, representing both water bodies with known stocking history and unaltered sites (Figure 1). We found that stocking did not explain the presence of FLB intergrades observed in North Carolina; however, stocking has likely increased the mean percent of FLB ancestry in those reservoirs. Future studies should incorporate the magnitude, duration, frequency, and timing of stocking in order to assess how these variables may influence the modeling results shown here. Unstocked river populations with many pure FLB fish were observed in Georgia south of the Fall Line, which is also in contrast to the currently accepted native range of M. floridanus. These results strongly support a redrawing of native ranges for both the intergrade zone and M. floridanus. To fully describe the extent of introgression in these lineages, more geographic sampling with the 73 SNP panel is needed, particularly in Mississippi and South Carolina.

One surprising result of our study was the high number of hybrids detected with ancestry contributions from all three lineages (367, 41% of samples analyzed with 73 SNP panel). As our panel was shown to have an 8.4% false discovery rate (Table 3), we can estimate that at least 336 of these triple hybrids are real. Populations with high proportions of triple hybrids (>20%) were primarily in the Atchafalaya Basin in Louisiana, reservoirs of Alabama and Georgia, and two sites in Tennessee. While the low FLB ancestry observed in Louisiana and Tennessee are likely due to stocking (Hargrove et al., 2019), we hypothesize that the triple hybrids in Georgia and Alabama are a result of secondary contact among the three lineages after allopatric divergence. Hybrids between three diverged lineages are more commonly detected in plants, but have been observed in some fishes such as cichlids, suckers, darters, and other micropterids (Keck & Near, 2009; Littrell et al., 2006; McDonald et al., 2008; Nevado et al., 2011).

Whether the result of natural introgression, anthropogenic introductions, or a combination of the two, the large hybrid zone observed between these three diverged lineages provides an excellent opportunity to study the speciation process (Barton & Hewitt, 1985; Gompert et al., 2017). When hybrid zones act as semipermeable barriers to gene flow between parental species, patterns of introgression can vary across the genome (Harrison & Larson, 2014). Differential introgression among loci can be due to stochastic processes (i.e., drift) or selection associated with genotype‐by‐environment interactions or reproductive isolation (Fitzpatrick et al., 2009; Gompert et al., 2012; Taylor et al., 2015). With additional genomic sequencing and a careful geographic sampling design aided by the diagnostic panel produced in this study, largemouth bass can be developed into a new case study for genomic patterns of introgression when multiple lineages are involved.

SNP assay for hybridization and species identification

Characterizing genomic contributions of hybrids is critical for understanding hybrid zone dynamics and informing fisheries management (Fitzpatrick et al., 2015). Empirical and simulation‐based studies have determined that at least 50 ancestry‐informative genetic markers may be required for accurately classifying F2 hybrids and advanced‐generation backcrossed individuals (Fitzpatrick, 2012; Malde et al., 2017). Previous work has resulted in panels of 18 microsatellites or 25–38 SNPs for assessing integrity and hybridization between M. salmoides/M. floridanus; however, these panels did not include samples from the Mobile‐Tensaw Delta during development and therefore may not accurately describe patterns of introgression (Li et al., 2015; Seyoum et al., 2013; Zhao et al., 2018). When comparing results to our 73 SNP panel, we find that M. floridanus ancestry inferred with a previous 35 SNP FLB–NLB assay is often inflated in sites with high DLB ancestry. For example, populations of pure (>94%) DLB ancestry in the Mobile Delta were shown to have >25% FLB ancestry when using the previous assay (Figure S3, File S2).

Designing a diagnostic SNP panel for more than two species is a challenge, as fixed triallelic SNPs are rare and therefore only biallelic SNPs are retained. Unlike a previous NLB–FLB SNP assay, we are not able to simply count the number of FLB and NLB alleles to determine hybrid ancestry proportions. Instead, we developed a rigorous approach for screening diagnostic loci by leveraging the program STRUCTURE and reference genotypes for accurate estimation of ancestry. Our panel of 73 SNPs captures the same ancestry proportions as the full set of 8582 GBS SNPs for all individuals within Q‐values of 0.1. By evaluating our panel on simulated hybrid genotypes, we also determined that a Q‐value cutoff of ≥0.94 was accurate for classifying whether an individual was likely derived from only one lineage. Many other studies that use STRUCTURE to identify pure individuals of a species use an arbitrary Q‐value cutoff, such as ≥0.95 (Lutz‐Carrillo et al., 2006; Thongda et al., 2020) or ≥0.90 (Dakin et al., 2015). When utilizing a reduced number of diagnostic markers in STRUCTURE, simulations should be applied to ascertain performance at various Q‐value thresholds.

Management implications for largemouth bass

As highly popular sportfish species, management of NLB and FLB involves balancing angler‐desirable traits (e.g., size, catchability), stable fisheries, and biodiversity concerns (Young et al., 2006). The phenotypic and demographic impacts of anthropogenic hybridization between the two species have been studied extensively, with some finding negative fitness consequences due to outbreeding depression (Cooke et al., 2001; Cooke & Philipp, 2006; Philipp & Claussen, 1995) and others seeing a potential to enhance managed fisheries (Maceina & Murphy, 1992; Maceina et al., 1988). Many of these early studies failed to take genetic ancestry into account or evaluated hybrids between geographically disparate sites, making it difficult to extrapolate results to naturally occurring intergrades. The discovery of a genetically distinct lineage of Largemouth Bass that may differ in fishery‐relevant traits (e.g., smaller size, slower growth rate, and environmental tolerances) further strengthens the need for accurate genetic tools to determine purity and hybridization. To better inform stocking success, future performance assays should assess genetic ancestry and leverage natural intergrades to provide a complete picture of fitness variation and environmental adaptation across largemouth bass lineages.

Current diagnostic markers in largemouth bass include a set of 35 SNPs designed for the MassARRAY Sequenom (Li et al., 2015; Zhao et al., 2018) and a panel of 18 microsatellites (Barthel et al., 2010; Seyoum et al., 2013). While we have not evaluated our new panel against these microsatellites, they likely face the same issues as the 35 SNP panel for mischaracterizing ancestry. Given the numerous populations with Delta bass ancestry observed in the southeast, genetic markers used by state agencies and evolutionary biologists need to be updated to better characterize the geographic extent of pure and hybrid populations of all three lineages. For example, some sites in Alabama and Georgia that were thought to have roughly the same FLB–NLB ancestry are now shown to differ considerably in their proportions of DLB ancestry. The continued movement and stocking of largemouth bass by state agencies and commercial hatcheries makes this a pressing issue. Genotyping of 5 “pure” northern Largemouth Bass from a popular commercial hatchery found DLB ancestry ranging from 0.02 to 0.17. Without understanding the potential genotype‐by‐environment interactions of these hybrids, DLB alleles may be spread far outside of their native range to detrimental outcomes.

The observation of individuals with DLB mitochondria and ancestry extending outside of the Mobile Delta raises the question of whether phenotypic variation observed in fish from the Delta are due to adaptation or phenotypic plasticity from the environment. While Largemouth Bass populations of the Mobile Delta are currently stable, they may be at particular risk from climate change due to unique environmental challenges in coastal habitats, such as low dissolved oxygen or storm‐mediated shifts in salinity. Microcosm experiments and functional genomics with both inland and coastal Delta bass can help elucidate whether these fish have adapted or acclimated to the estuarine environment.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

Supplementary Material

Click here for additional data file.

File S2

Click here for additional data file.

File S3

Click here for additional data file.