A chromosomal inversion may facilitate adaptation despite periodic gene flow in a freshwater fish

Abstract Differences in genomic architecture between populations, such as chromosomal inversions, may play an important role in facilitating adaptation despite opportunities for gene flow. One system where chromosomal inversions may be important for eco‐evolutionary dynamics is in freshwater fishes, which often live in heterogenous environments characterized by varying levels of connectivity and varying opportunities for gene flow. In the present study, reduced representation sequencing was used to study possible adaptation in n = 345 walleye (Sander vitreus) from three North American waterbodies: Cedar Bluff Reservoir (Kansas, USA), Lake Manitoba (Manitoba, Canada), and Lake Winnipeg (Manitoba, Canada). Haplotype and outlier‐based tests revealed a putative chromosomal inversion that contained three expressed genes and was nearly fixed in walleye assigned to Lake Winnipeg. These patterns exist despite the potential for high gene flow between these proximate Canadian lakes, suggesting that the inversion may be important for facilitating adaptive divergence between the two lakes despite gene flow. However, a specific adaptive role for the putative inversion could not be tested with the present data. Our study illuminates the importance of genomic architecture consistent with local adaptation in freshwater fishes. Furthermore, our results provide additional evidence that inversions may facilitate local adaptation in many organisms that inhabit connected but heterogenous environments.


| INTRODUC TI ON
Differences in genomic architecture such as structural variants may help to facilitate local adaptation by suppressing recombination and protecting co-adapted alleles, and are therefore key mechanisms of eco-evolutionary processes (Dorant et al., 2020;Shi et al., 2021;Tigano & Friesen, 2016;Wellenreuther & Bernatchez, 2018). One type of structural variant, chromosomal inversions, are taxonomically widespread reversed regions of DNA that have been linked to both adaptation and speciation (reviewed in Wellenreuther & Bernatchez, 2018). For instance, chromosomal inversions have had important eco-evolutionary consequences in sunflowers (Helianthus spp.) (Barb et al., 2014), fruit flies (Drosophila melanogaster) (Corbett-Detig & Hartl, 2012), swallowtail butterflies (Papilio polytes) (Nishikawa et al., 2015), and birds (reviewed in Hooper & Price, 2017). Inversions can effectively protect large chromosomal regions containing hundreds of genes from recombination, facilitating local adaptation even in the face of high gene flow (Wellenreuther & Bernatchez, 2018). Most inversions found in previous studies have been relatively large because sequencing methods can more easily detect large inversions (Wellenreuther & Bernatchez, 2018).
Freshwater fishes often live in heterogenous environments characterized by varying levels of connectivity and therefore, varying opportunities for gene flow (Griffiths, 2015;Mushet et al., 2019).
Gene flow tends to act in opposition to local adaptation (Lenormand, 2002) but nevertheless, evidence for local adaptation despite gene flow has been observed among different fishes that use freshwater habitats (see Fraser et al., 2011 for a review of salmonids; Shi et al., 2021). While genomic architecture has been proposed to facilitate local adaptation (Tigano & Friesen, 2016), few studies test for this possibility in freshwater fishes (Shi et al., 2021). Freshwater habitats contribute a disproportionate richness in species diversity to global species diversity, given that only 0.01% of the planet's total surface is fresh water (Balian et al., 2008). Fishes are a significant proportion of overall freshwater species diversity (approx. 10% of overall freshwater species and 70% of vertebrate freshwater species are fish; Balian et al., 2008). Therefore, an understanding of the mechanisms of local adaptation in freshwater fishes can contribute substantially to our overall understanding of local adaptation in heterogenous environments.
Walleye (Sander vitreus, but see Bruner, 2021) is a freshwater fish with a native range from the Northwest Territories (Canada) to Alabama (USA), approximately 4,600 kilometers apart (Hartman, 2009), and support important commercial and recreational fisheries (Fisheries & Oceans Canada, 2019. The economic importance of walleye underscores the need for research supporting the sustainable management of walleye fisheries, while its wide native distribution provides an opportunity to study the genomic basis of adaptation to a variety of freshwater environments. In addition, northern populations of walleye (such as those in Manitoba, Canada) were suspected to have been recently isolated into smaller lakes following the drainage of glacial Lake Agassiz approximately 7,000 years ago (Rempel & Smith, 1998;Stepien et al., 2009). This event may have resulted in walleye populations that remain closely related at neutral genetic variants but show evidence of divergence at adaptive genetic variants, such as in chromosomal inversions.
We sought to study how heterogenous environments with varying connectivity may have shaped contemporary walleye populations, and the potential for genomic architecture to facilitate adaptation despite opportunities for gene flow. Walleye from Lake Manitoba and Lake Winnipeg in Manitoba (Canada) were analyzed and compared with an entirely stocked outgroup of unknown origin from Cedar Bluff Reservoir in Kansas (USA) (Figure 1). Walleye in Lake Manitoba and Lake Winnipeg live in largely separated waterbodies with a shared river-Lake Manitoba drains eastward into Lake Winnipeg via the Fairford River, Lake St. Martin, and the Dauphin River ( Figure 1). As a mobile species that undergoes broad seasonal migrations, historic gene flow may have been extensive between these Manitoba waterbodies (Backhouse-James & Docker, 2012;Munaweera Arachchilage et al., 2021;Thorstensen et al., 2020;Turner et al., 2021). In addition, historical flooding in 1882, 1902, 1904, 2011, and 2014 brought significantly increased waterflow eastward from Lake Manitoba toward Lake Winnipeg (Ahmari et al., 2016), and such flooding could have facilitated increased fish movement between the waterbodies. Last, walleye K E Y W O R D S adaptive variation, genomic architecture, haplotype, Sander vitreus, selection, walleye

T A X O N O M Y C L A S S I F I C A T I O N
Conservation genetics; Ecological economics; Ecological genetics; Genetics; Genomics; Population genetics fry were opportunistically stocked in Lake Winnipeg from Lake Manitoba (Manitoba, Canada) between 1917(Manitoba Government, 2020. While intensive stocking complicated analyses of walleye ancestry in walleye of Minnesota and Wisconsin (USA) (Bootsma et al., 2021), stocking of walleye fry, which likely experience significant mortality, had little to no effect in South Dakota and Missouri (USA) (Fielder, 1992;Koppelman et al., 1992). Therefore, the genetic impacts of stocking in this system are likely minimal as stocking was generally opportunistic and used only fry.
Here, we employed Rapture (Ali et al., 2016), a reduced representation approach, to genotype n = 345 individuals from the three waterbodies described above. Population structure and demographic F I G U R E 1 Map of waterbodies and ancestry found for walleye (Sander vitreus) included in the present study. Cedar Bluff Reservoir, Kansas, USA, represents walleye in an entirely stocked population of unknown origin. Lake Manitoba and Lake Winnipeg (Manitoba, Canada) walleye represent native populations, although Lake Winnipeg has been stocked by Lake Manitoba walleye throughout the 20th century. In addition, floods in 1882, 1902, 1904, 2011, and 2014 have carried water from Lake Manitoba to Lake Winnipeg, and may have contributed to walleye gene flow. Pie charts represent estimated ancestry proportions with Admixture on n = 345 walleye total from K = 3 groups. Each pie chart corresponds to a sampling site, with the Red River, Matheson Island, and the Dauphin River in Lake Winnipeg, Swan Creek Hatchery in Lake Manitoba, and Cedar Bluff Reservoir in Kansas history were explored to provide context for signatures of selection.
Simulations with EASYPOP were used to estimate historical migration rates between Lake Manitoba and Lake Winnipeg (Balloux, 2001). Candidate regions of adaptive variation were identified with cross-population extended haplotype homozygosity (XP-EHH) (Gautier et al., 2017;Sabeti et al., 2007), and outlier tests were used to confirm haplotype-based results. In addition, mRNA transcript expression in one candidate chromosomal inversion was assessed for Lake Winnipeg walleye to show the potential for the functional significance of genomic architecture in this system. The present study illuminates the importance of genomic architecture for facilitating local adaptation in freshwater fishes, and provides additional evidence that inversions may facilitate local adaptation in many organisms that inhabit connected but heterogenous environments.

| Sample collection
Walleye anal fin clips from Lake Winnipeg individuals were collected Raw reads were processed with STACKS v2.3 (Catchen et al., 2011(Catchen et al., , 2013, and mapped to the yellow perch (Perca flavescens) genome with bwa v0.7.17 using bwa mem (Feron et al., 2020;). Mapping percentages of walleye data to the yellow perch genome were checked with samtools flagstat .
Single nucleotide polymorphisms (SNPs) were called with one random SNP per locus and following best practices with RAD data (Linck & Battey, 2019), then subsequently filtered for paralogy with HDPlot (McKinney et al., 2017). 72,149 SNPs were initially called from mapped reads, and 63,882 used with pcadapt (Luu et al., 2017;Privé et al., 2020) were retained after filtering for paralogy with HDPlot. Data for population structure and demographic history were pruned for linkage with PLINK v1.9 (Purcell et al., 2007) because we reasoned that some linkage blocks may extend beyond 5 kb, the approximate distance between target loci for the RAPTURE panel used here. 46,342 SNPs were retained after pruning for linkage. Because relatedness can introduce bias in genetic studies, pairwise relatedness was checked using the method of moments as implemented by PLINK in SNPRelate v1.22 (Zheng et al., 2012).

| Population genetics
Admixture v1.3.0 was run over K values of one through six, with 10-fold cross validation and 1,000 bootstrap replicates for parameter standard errors (Alexander et al., 2009). Pophelper v2.3.0 in R was used to visualize and organize Admixture results (Francis, 2017).
For population assignments, individuals were considered assigned to a cluster when their Q-values were >0.85 for that cluster. Only assigned individuals were used for population differentiation and demographic reconstruction.
Hierfstat v0.5-7 was used to find β WT as a measure of population-specific differentiation from the entire pool and Weir & Cockerham's pairwise F ST between populations (Goudet, 2005;Weir & Cockerham, 1984;Weir & Goudet, 2017). Ninety-five percent confidence intervals were generated for F ST and β WT over 1,000 bootstrapped iterations each. In addition, observed heterozygosity (H O ), gene diversity (H S ), and inbreeding coefficients (F IS , 95% confidence intervals generated over 1,000 bootstrapped iterations) for each population were calculated.
NeEstimator v2.1 was used to estimate effective population size in 95% confidence intervals for each population using the linkage disequilibrium method and only comparing SNPs on different chromosomes (Do et al., 2014). Here, datasets were filtered for 10% missing data per population using population assignments with vcftools, pulled from the pruned SNP dataset. EASYPOP v2.0.1 (Balloux, 2001) was used to test different scenarios of historic gene flow that may have led to contemporary F ST between Lake Winnipeg and Lake Manitoba. Migration rates varied between 0, 0.0001, 0.001, 0.002, 0.003, 0.004, 0.005, 0.01, 0.02, and 0.03 individuals per generation, with migration rates held constant within each simulation run. Population sizes of 4,500 individuals in Lake Winnipeg and 350 individuals in Lake Manitoba were assumed constant with equal males and females, based on estimates of N e (see results). In EASYPOP, 1,000 simulated loci with two allelic states were used with free recombination allowed, using a KAM mutation model with mutation rate µ of 3.28x10 −9 based on divergence between the blackfin icefish (Chaenocephalus aceratus) and dragonfish (Parachaenichthys charcoti) (Kim et al., 2019). Maximal variability (i.e., randomly assigned alleles) was chosen for the initial population.
Each simulation was run over 930 generations, based on the 4.3- year generation time estimated in Franckowiak et al. (2009) and an approximate lower bound of 4,000 radiocarbon years since Lake Winnipeg had a reduced northern lakebed (Lewis et al., 2002). This reduced northern lakebed would have prevented walleye migration between populations via the Dauphin River and Lake St. Martin in the north basin ( Figure 1). 100 replicate runs of EASYPOP were used for each migration rate tested, with 1,000 simulation runs total.

For each run, Weir & Cockerham's pairwise F ST was measured with
Hierfsat. These F ST values were compared to the F ST observed between Lake Winnipeg and Lake Manitoba walleye populations using linkage disequilibrium-pruned SNPs.

| Linkage disequilibrium decay
The extent to which SNPs in the RAPTURE panel represented linkage disequilibrium (LD) among walleye populations was analyzed via LD decay with PopLDdecay v3.40 (Zhang et al., 2019). Within each population run, SNPs were only included for LD decay analysis with a maximum of 10% missing data and minor allele frequency ≥0.05 because rare alleles lead to increased variance in r 2 (Remington et al., 2001). LD was modeled with a nonlinear least-squares approach to estimate expected r 2 between SNPs (Hill & Weir, 1988;Remington et al., 2001). From the nonlinear regression, the distances' half-decay where estimated LD fell to half of its maximum estimated value, and of moderate LD where r 2 ≤ 0.20 were estimated for each population.

| Signatures of selection
Signatures of selection, where one genomic region has approached or reached fixation in one population compared to another, were explored with XP-EHH using the R package rehh v3.2.1, with the markers previously described as phased and imputed after filtering for 90% present data (Gautier et al., 2017;Sabeti et al., 2007). Data were unpolarized for analyzing XP-EHH because ancestral and derived alleles were unknown. Individuals with phased and imputed SNPs were split into the three assigned populations using Admixture results (Q>0.85 per individual), and extended haplotype homozygosity (EHH) was analyzed for each population separately (scan_hh).
Three pairwise comparisons between each population were used to identify XP-EHH using false discovery rate corrected p-values (q-values). Significant candidate regions of selection were identified by analyzing 100 kilobase (kb) windows overlapping by 10 kb, in which at least three significant SNPs (q < 0.05) showing XP-EHH were found. For visualization, only significant SNPs within candidate regions were highlighted (thus a significant SNP outside a candidate region would be de-emphasized), using the R package ggman v0.99.0 with relative SNP positions to reflect distances between points in the reduced representation data (https://rdrr.io/githu b/veera -dr/ ggman/). Results for XP-EHH were divided into three pairwise tests between each assigned population (Cedar Bluff Reservoir, Lake Manitoba, and Lake Winnipeg), and significant results for each pairwise test were further categorized into which population showed elevated XP-EHH scores (e.g., in a pairwise test between Lake Manitoba and Lake Winnipeg, certain SNPs were significant among haplotypes in Lake Manitoba versus others significant among haplotypes in Lake Winnipeg).
The program pcadapt v4.3.3 was used to explore potential signatures of local adaptation with outlier tests using non-paralogous SNPs unpruned for linkage disequilibrium (Luu et al., 2017;Privé et al., 2020). Here, a scree plot showed that the optimal choice for principal components (PC) was K = 2 based on Cattell's rule, an observation confirmed by a lack of discernible population structure in K > 2 principal components (Cattell, 1966) (Figure S1). Significance for individual SNPs was determined using q values implemented in the R package qvalue v2.20.0, where a SNP was accepted as significant at q < 0.05 (Storey et al., 2020). The PC associated with a significant SNP was retrieved with get.pc, which enabled us to link significant SNPs with different axes of population structure. pcadapt results were separated into datasets of those significant along PC1 or PC2, corresponding to latitude and longitude, respectively.
Hierfstat v0.5-7 was used to provide supporting evidence for selection at individual SNPs using F' ST with basic.stats on the same set of SNPs used with pcadapt. Here, F' ST was used as a sample size corrected F ST (via a sample size correction to gene diversity among samples) and measured in pairwise comparisons between each of three populations defined by Admixture. In addition, absolute allele frequency differences were analyzed between populations. VCFtools was used to calculate allele frequencies for each population. Minor alleles in Lake Winnipeg-assigned individuals were used for comparisons across populations.

| Putative inversion analysis
To assess evidence that a region of interest that we identified (chromosome 8, position 15,260,000-15,900,000; see results for more detail) was a putative inversion, we conducted the following analyses as suggested by Huang et al. (2020) and Shi et al. (2021). First, to test whether the region showed elevated linkage disequilibrium or LD (r 2 ), we calculated r 2 using PLINK v1.9 for SNPs on chromosome 8 (--r2 --ld-window-r2 0.05 --ld-window 999999999 --ld-window-kb 30000) and generated a LD heatmap (Chang et al., 2015;Purcell et al., 2007). Second, we conducted regional PCA using SNPs within the region to identify whether individuals were grouped into three clusters along PC1, which is characteristic of chromosomal inversions. The discreteness of the clustering was calculated as the proportion of the between-cluster sum of squares over the total using the R function kmeans in adegenet (Jombart, 2008). Third, we compared heterozygosity (the proportion of heterozygotes) among the three identified PCA clusters using Wilcoxon tests (α = 0.05) to further confirm that individuals in the middle PCA cluster was heterozygous for the putative inversion. Lastly, we calculated genotype frequencies of the putative inversion in all populations and assumed that the more derived inversion arrangement would likely have had lower heterozygosity given its relatively recent origin compared to the ancestral state (Knief et al., 2016;Laayouni et al., 2003;Twyford & Friedman, 2015) but see Matschiner et al. (2022). Individual assignments for the three identified PCA clusters for the putative inversion were then visualized with respect to population assignments and site collected, to assess the frequency of the inversion across the different waterbodies studied.
To provide evidence of potential functional significance for the identified putative inversion on chromosome 8, gene expression data was collected for genes within the inversion. RNA-seq data was used from n = 48 walleye in three sites in Lake Winnipeg, using previously published reads (NCBI SRA database accession #PRJNA596986) (Thorstensen et al., 2020). Three genes in the yellow perch genome were within the putative inversion on chromosome 8: PDHX, EHF, and LRRC4C. Therefore, expression of these genes was assessed with counts per million (CPM) values from the walleye transcriptome-aligned data Patro et al., 2017;Soneson et al., 2015). Differential gene expression was not tested for because the putative inversion was nearly fixed in the Lake Winnipeg-assigned walleye; differential expression would thus be expected between populations and not within populations.

| Population genetics
Admixture identified K = 2 populations based on the lowest cross validation error between K = 1 through 6, separating Canadian (i.e., Lake Winnipeg and Lake Manitoba) and Kansas (i.e., Cedar Bluff Reservoir) populations of walleye (Figure 1 and S3). However, cross validation error was similar at K = 3 and population structure delineated the three different water bodies used in the present study ( Figure 1); individuals were thus assigned to populations based on K = 3. Using ancestry coefficients Q > 0.85 for population assignments, 60 individuals were assigned to the Kansas population, 67 individuals to the Lake Manitoba population, 189 individuals to a Lake Winnipeg population, and 29 to no population. Notably, 27 out of the 29 unassigned individuals were found in Lake Winnipeg (the 2 remaining in Lake Manitoba), and of those, 18 were caught in the northern Dauphin River, 6 in the central Matheson Channel, and 2 in the southern Red River. In addition, out of the 67 individuals assigned to the Lake Manitoba population, 20 were sampled from the Dauphin River in Lake Winnipeg.

Pairwise F ST showed the greatest differentiation between each
Canadian lake and Cedar Bluff Reservoir, and moderate differentiation between Lake Manitoba and Lake Winnipeg (Figure 2).
Population-specific differentiation from the overall pool, β WT , was lowest for Cedar Bluff Reservoir ( Figure S4).

| Linkage disequilibrium decay
For Lake Winnipeg, half-decay was at approximately 2.4 kb, while moderate LD (r 2 ≤ 0.20) was at approximately 2.8 kb. For Lake Manitoba, half-decay was at approximately 4.2 kb, while moderate LD was at approximately 5.7 kb. For Cedar Bluff Reservoir, halfdecay was at about 1.9 kb, while moderate LD was at approximately 2.6 kb. And for unassigned individuals, half-decay was at approximately 5.6 kb, while moderate LD was at about 8.0 kb. Individual pairwise r 2 , distance, and nonlinear regressions were visualized for each population in Figures S6 and S7.

| Signatures of selection
Among XP-EHH scores between waterbodies, the greatest number of SNPs (45) were significant between Lake Winnipeg and Lake Manitoba (Table 1, Figure 3). There were 10 candidate regions under selection in XP-EHH between Lake Winnipeg and Lake Manitoba, 5 candidate regions between Lake Winnipeg and Cedar Bluff Reservoir, and 0 candidate regions between Lake Manitoba and Cedar Bluff Reservoir (Table   S1). Prominent among candidate regions was a section of chromosome 8 between approximately 15.26 and 15.90 Mb in Lake Winnipeg with high differentiation compared to both Lake Manitoba and Cedar Bluff Reservoir. In this chromosomal region specifically, 7 SNPs showed elevated XP-EHH in Lake Winnipeg relative to Lake Manitoba, out of 11 SNPs within candidate regions total in that comparison (Table 1).  The region between 15.26 and 15.90 Mb along chromosome 8 was notable for its candidate regions identified in XP-EHH that were unique to Lake Winnipeg (Table S1)

| Putative inversion analysis
We found strong evidence of the genomic region 15.26-15.90 Mb on chromosome 8 as a putative chromosomal inversion ( Figure 5). This region demonstrated elevated LD with a mean r 2 of 0.30 and values up to 0.96 (Figure 5b), which appeared as an outlier overall on chromosome 8 (Figure S10). Regional PCA showed that individuals were clustered into three distinct groups along PC1 with a discreteness of 0.9626 (Figure 5c) and the middle PCA cluster displaying significantly higher heterozygosity than the other two clusters (Figure 5d and S11). There were also significant differences in heterozygosity between the two homokaryotypes ( Figure 5d) and the arrangement with lower heterozygosity (cluster 2) was assumed to be the derived inverted type (but see Matschiner et al., 2022). This putative inversion (cluster 2) was nearly fixed in the Red River and Matheson Island sites of Lake Winnipeg (toward the channel and south basin;

| DISCUSS ION
We used a reduced representation approach to genotype n = 345 walleye from three different waterbodies in North America at 46,342 genetic markers. Walleye of the Canadian lakes diverged from each other despite three scenarios for gene flow: a connecting river, stocking of fry throughout the 20th century, and flooding in 1882, 1902, 1904, 2011, and 2014 (and  Mb. An inversion may have played a role in adaptive divergence between Lake Winnipeg and Lake Manitoba walleye (Shi et al., 2021;Wellenreuther & Bernatchez, 2018), and three expressed genes (i.e., mRNA transcripts) within it indicate at least the potential for functional significance of the putative inversion.

| Moderate differentiation between two connected lakes
Despite moderate signals of differentiation between the Canadian lakes, admixture analysis and population assignment identified 11% of individuals sampled from primarily the north basin of Lake TA B L E 1 Number of single nucleotide polymorphisms (SNPs) significant in pairwise comparisons with cross-population extended haplotype homozygosity (XP-EHH)

Number of SNPs (number in candidate regions) Total SNPs (number in candidate regions)
Lake Winnipeg vs. Lake Manitoba 15 (11) (0) 17 (5) Note: Three waterbodies were compared: Lake Winnipeg and Lake Manitoba in Manitoba, Canada, and Cedar Bluff Reservoir in Kansas, USA. In each pairwise comparison, significant SNPs (q < 0.05) are presented in the same order as waterbodies. In parentheticals are SNPs within candidate regions of selection defined by 100 kilobase regions of the genome overlapping by 10 kilobases, in which ≥3 SNPs total were significant.

F I G U R E 3
Signatures of selection inferred from cross-population extended haplotype homozygosity (XP-EHH) between each population of walleye (Sander vitreus) analyzed in the present study, along with absolute allele frequency differences between two waterbodies. In (a), Cedar Bluff Reservoir (Kansas, USA) represents an entirely stocked population, while Lake Manitoba and Lake Winnipeg (Manitoba, Canada) each represent native populations with possible gene flow. Walleye DNA was aligned to the yellow perch (Perca flavescens) reference genome, and chromosome numbers represent yellow perch chromosomes. Unplaced scaffolds are not visualized. Single nucleotide polymorphisms (SNPs) were considered significant with q < 0.05 and only if the SNPs were in a candidate region defined by 100 kilobase regions of the genome overlapping by 10 kilobases, in which ≥3 SNPs total were significant. Because XP-EHH focuses on haplotypes of unusual length and high frequency in one population compared to another, the population in which haplotypes are long and are at or approaching fixation can be identified. Here, the population in which XP-EHH was significant is identified with colors, where blue represents Lake Winnipeg, red represents Lake Manitoba, and yellow represents Cedar Bluff Reservoir. In (b), absolute allele frequency differences were calculated between Lake Manitoba and Lake Winnipeg-assigned walleye. PC1 and PC2 refer to significant SNPs from principal component 1 or 2 from pcadapt (see Figure 4). The genes pyruvate dehydrogenase protein X component (PDHX), ETS homologous factor (EHF), and leucine-rich repeat-containing protein 4C (LRRC4C) were within the candidate region for selection in Lake Winnipeg on chromosome 8. The candidate regions for selection from XP-EHH are emphasized with gray background, which together were identified as the putative chromosomal inversion. Arrows below the yellow regions highlighting each gene indicate the direction of the gene in the genome Winnipeg as originating from Lake Manitoba. This gene flow into Lake Winnipeg may have been associated with stocking programs from Lake Manitoba to Lake Winnipeg in the 1900s and early 2000s, while floods may have facilitated movement in 2011 and 2014 (Ahmari et al., 2016). However, individuals possibly admixed between Lake Manitoba and Lake Winnipeg were likely the direct descendants (i.e., which were attributed to Pleistocene lineages (Bootsma et al., 2021). Walleye population differentiation between Lake Winnipeg and Lake Manitoba was consistent with certain other walleye populations within watersheds (e.g., Grand River, Ontario, Canada vs other Lake Erie stocks in Euclide et al., 2021; Sarah Lake vs. Lake Koronis, Minnesota, USA in Bootsma et al., 2021). Retreating glacial refugia dispersed walleye into different watersheds approximately 7,000-10,000 years ago (Bootsma et al., 2021;Stepien et al., 2009).
Therefore, the population structure identified in the present study represents differentiation following geologically recent changes to the overall landscape. This differentiation is consistent with high spawning site fidelity previously observed in walleye (Horrall, 1981;

F I G U R E 4
Principal components analysis of walleye (Sander vitreus) from three waterbodies across North America. Cedar Bluff Reservoir (Kansas, USA) represent an entirely stocked population of walleye of unknown origin. Lake Winnipeg and Lake Manitoba (Manitoba, Canada) represent native populations of walleye with possible gene flow. Principal components analysis was done using pcadapt, while population assignments were done using Admixture with K = 3 groups. Individuals were assigned to a population with Q > 0.85 for a population, and individuals with maximum Q ≤ 0.85 were considered unassigned. Significant SNPs were distinguished according to the first two principal components in which they most strongly contributed variation, where principal component 1 (PC1) corresponded to latitudinal differences among populations, while principal component 2 (PC2) was most strongly characterized by variation between Lake Winnipeg and Cedar Bluff Reservoir relative to Lake Manitoba. Only SNPs with q-values <0.05 were highlighted for visualization. Walleye DNA was aligned to the yellow perch (Perca flavescens) reference genome. Unplaced scaffolds were not visualized, and chromosome numbers refer to chromosomes in the yellow perch genome Jennings et al., 1996;Stepien & Faber, 1998;Stepien et al., 2009

Given widespread stocking and potamodromous movements in
walleye, it may be surprising that the species showed as much population differentiation as has been observed between waterbodies in multiple watersheds (Bootsma et al., 2021;Euclide et al., 2021;Munaweera Arachchilage et al., 2021;Raby et al., 2018;Turner et al., 2021). While stocking has led to extensive homogenization in some walleye populations, the lack of gene flow from stocking in the present data and among other waterbodies implicates potential mechanisms by which gene flow is prevented (Bootsma et al., 2021). Possible mechanism include environmental interactions with fitness. Here, environmental factors impede gene flow because local adaptation may select against migrants (Hendry, 2004;Sexton et al., 2014). In a comparison between environmental mechanisms of isolation and isolation by distance, isolation by environment was more common in all taxa, including vertebrates (Sexton et al., 2014). As such, environmental differences between source and stocked waterbodies may have precluded gene flow in different walleye systems. Genomic architecture is another possible mechanism that can maintain differentiation despite stocking. Here, physical linkage or chromosomal rearrangements have been hypothesized to play an important role in local adaptation in the face of gene flow (Tigano & Friesen, 2016). Walleye in waterbodies with historical pulses of gene flow or recent stocking may thus maintain differentiation because of genomic architecture facilitating local adaptation (Tigano & Friesen, 2016). Moreover, the hypotheses that environment or genomic architecture may have maintained differentiation despite opportunities for gene flow are not mutually exclusive. Because genomic architecture has played a role in local adaptation in several freshwater fishes (Shi et al., 2021;Wellenreuther & Bernatchez, 2018), and environment is a key variable in local adaptation (Coop et al., 2010;Forester et al., 2018), these two mechanisms may often work in tandem.

| A chromosomal inversion may facilitate adaptive divergence
A putative chromosomal inversion was identified on chromosome 8. Based on PCA, absolute allele frequency differences, and F ' ST , this 640 Kb region on chromosome 8 showed evidence as an outlier region between Lake Manitoba and Lake Winnipeg. For example, the top 10 outlier SNPs (by -log 10 q-value) out of 119 significant outliers (q < 0.05) along PC 2 in a pcadapt analysis are all from this small outlier region (Luu et al., 2017;Privé et al., 2020). Similarly, the top 10 SNPs by absolute allele frequency difference between Lake Manitoba and Lake Winnipeg were also from this small outlier region. As a haplotype-based test, XP-EHH provided evidence that this region was selected for not in Lake Manitoba, but in Lake Winnipeg instead (Sabeti et al., 2007). The extent of linkage in this region was far beyond overall distances of half and moderate LD decay in Lake Winnipeg or Lake Manitoba walleye, as well. This chromosome 8 region was also a consistent outlier in scans for adaptive loci in a different study on walleye in Wisconsin and Minnesota (USA) (Bootsma et al., 2021). However, SNP data were too sparse to interrogate the genomic architecture of this outlier region. While a selective sweep from recent positive selection may underlie such a small chromosomal region of high differentiation (Sabeti et al., 2007), these results were also consistent with some difference in genomic architecture in the observed populations. In the present study, a dense panel of SNPs provided sufficient resolution to investigate the possibility that this outlier region was a chromosomal inversion. However, specific inversion breakpoints are unknown, and may be within the approximate region described here because recombination suppression extends beyond inversion breakpoints (Stevison et al., 2011).
Linkage disequilibrium, regional PCA, discreteness of clustering, and heterozygosity provided evidence for a putative chromosomal inversion at higher frequency in Lake Winnipeg relative to both other waterbodies. This putative inversion was nearly fixed in Lake Winnipeg-assigned walleye, and at intermediate frequen- in Lake Winnipeg (the majority of which were found in the northern Dauphin River), the putative inversion was nearly fixed for the Lake Winnipeg genotype. Therefore, this putative inversion may play a role in walleye adaptation to Lake Winnipeg, but a specific adaptive role for the inversion could not be tested with the present data.
Three genes were within the putative inversion: PDHX, EHF, and LRRC4C. All three were expressed in Lake Winnipeg walleye, indicating at least some gene transcription activity within the putative inversion. While functional information from model species exists for each of these genes, their specific functional significance in Lake Winnipeg walleye is unknown. The specific life stage or environmental context in which the three genes within the putative inversion may have an effect on walleye biology is similarly unknown. Moreover, genes adjacent to the inversion may have functional consequences, and possibly in another group of walleye instead (Matschiner et al., 2022). However, several other characteristics of Lake Winnipeg walleye are consistent with signatures of selection unique to the waterbody. Observations of recentlyevolved dwarf walleye morphotypes possibly due to fishing pressure suggest ongoing selection within Lake Winnipeg (Moles et al., 2010;Sheppard et al., 2018). Additionally, there are patterns of low gonadal investment despite high lipid concentrations in Lake Winnipeg walleye relative to walleye in other lakes, indicating a possibly unusual phenotype in walleye of Lake Winnipeg relative to others (Moles et al., 2008). Walleye habitats are different between Lake Manitoba and Lake Winnipeg, as well. Lake Winnipeg is one of the largest lakes in the world based on surface area, but is relatively shallow compared with other large freshwater lakes (Brunskill et al., 1980;ECCC & MARD, 2020). However, Lake Manitoba is shallower still, with a 4.9 m mean depth (compared with a mean depth for Lake Winnipeg of 13.3 m and 9 m in the north and south basins, respectively) (ECCC & MARD, 2020; Rawson, 1952). Drainage from the more saline Lake Manitoba increases salinity in Lake Winnipeg's north basin, but not in Lake Winnipeg's south basin (Brunskill et al., 1980;ECCC & MARD, 2020;ECCC, 2011). Consequently, possible local adaptation in Lake Winnipeg may be an evolutionary response to environmental differences in the waterbody.

| CON CLUS ION
The identification of a putative inversion contributes to a broad literature on chromosomal inversions in diverse taxa (see Wellenreuther & Bernatchez, 2018), although information on inversions is relatively scarce for obligate freshwater fishes (Arostegui et al., 2019;Penso-Dolfin et al., 2020;Roesti et al., 2015;Shi et al., 2021). In addition, many observed inversions have been >1 Mb in length, possibly because sequencing methods may be biased for detecting large inversions (Wellenreuther & Bernatchez, 2018). The present data show a small putative inversion in an obligate freshwater fish discovered via reduced representation sequencing (i.e., Rapture). That this putative inversion may have played a role in recent divergence between fish of two connected lakes, and that its frequency among likely admixed fish is consistent with a possible adaptive role in Lake Winnipeg walleye, indicates a potential importance for chromosomal inversions in freshwater fishes more generally. Inversions have been previously related to environmental adaptation, such as in freshwater emerald shiner (Notropis atherinoides) with high gene flow (Shi et al., 2021;Wellenreuther & Bernatchez, 2018). Heterogeneous habitats with opportunities for gene flow are common for freshwater fishes (Griffiths, 2015;Mushet et al., 2019), and chromosomal inversions, along with other genomic architectures, may facilitate local adaptation in many organisms that live in such connected environments.

ACK N OWLED G EM ENTS
We thank Dr. Jeff Long for contributing walleye fry from Swan Creek Hatchery and for valuable discussions of Manitoba walleye ecology, diversity, and stocking. We thank Colin Charles, Colin Kovachik, Doug Leroux, Nicole Turner, Mike Gaudry, Sarah Glowa, and Emily Barker, who provided fin clips used for walleye DNA from Lake Winnipeg. We thank Evelien de Greef for guidance on synteny

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

O PE N R E S E A RCH BA D G E S
This article has been awarded Open Data, Open Materials Badges.
All materials and data are publicly accessible via the Open Science