Size reductions and genomic changes associated with harvesting within two generations in wild walleye populations

The extent and rate of harvest-induced genetic changes in natural populations may impact population productivity, recovery and persistence. While there is substantial evidence for phenotypic changes in harvested fishes, knowledge of genetic change in the wild remains limited, as phenotypic and genetic data are seldom considered in tandem, and the number of generations needed for genetic changes to occur is not well understood. We quantified changes in size-at-age, sex-specific changes in body size, and genomic metrics in three harvested walleye (Sander vitreus) populations and a fourth reference population with low harvest levels over a 15-year period in Mistassini Lake, Quebec. We also collected Traditional Ecological Knowledge (TEK) surrounding concerns about these populations over time. Using ∼9000 SNPs, genomic metrics included changes in population structure, neutral genomic diversity, effective population size and signatures of selection. TEK revealed concerns about overall reductions in body size and number of fish caught. Smaller body size, smaller size-at-age, changing population structure (population differentiation within one river and homogenization between two others), and signatures of selection between historical and contemporary samples reflected coupled phenotypic and genomic change in the three harvested populations in both sexes, while no change occurred in the reference population. Sex-specific analyses revealed differences in both body size and genomic metrics but were inconclusive about whether one sex was disproportionately affected. Our results support that harvest-induced genetic changes can arise within 1-2.5 generations in long-lived wild fishes, demonstrating the need to investigate concerns about harvest-induced evolution quickly once they have been raised.


Introduction
fishing may drive genetic changes in the life history (e.g. body size, size-at-age, by sex), and genetic 78 characteristics (e.g. population structure, genetic diversity and composition) of wild populations. 79 While Western Scientific Methods (WSM) are most often used to inform fisheries management, 80 inclusion of Traditional Ecological Knowledge (TEK) has become an integral complement to scientific 81 knowledge for wildlife management and community-based conservation (Berkes et al. 2000;Fraser et al. 82 2006;Polfus et al. 2016;Polfus et al. 2014). TEK is defined as the "cumulative body of knowledge, 83 practice and belief, evolving by adaptive processes and handed down through generations by cultural 84 transmission, about the relationship of living beings (including humans) with one another and with their 85 environment" (Berkes et al. 2000). Importantly, TEK provides extensive location-specific knowledge, can 86 detect changes in wildlife more quickly than WSM (Huntington 2011), and often provides increased 87 knowledge of environmental linkages (Chapman 2007;Drew 2005). 88 Walleye (Sander vitreus) are important for commercial, sport and Indigenous subsistence 89 fisheries across North America (Bozek et al. 2011;Hansen et al. 2015;Scott & Crossman 1979). 90 Mistassini Lake in northern Quebec, Canada, is the province's largest natural lake (161 km long, 2 335 91 km 2 , 183 m maximum depth), is in Grand Council of the Crees land, Eeyou Istchee, and is considered to 92 be largely pristine (minimal mining, forestry, development; no known invasive species) (Fraser et al. 93 2006;Marin et al. 2017). The motivation behind this study was observations by Cree elders and fishers of 94 reduced body size and catch rates in walleye populations in three of Mistassini Lake's southern tributaries 95 that are close to the community, and a desire by the community to determine if management actions were 96 needed. We also studied a fourth river at the northeastern tip of the lake, where the population was 97 perceived to be largely unaffected by fishing until very recently (~2015, TEK, see methods). Subsistence 98 harvest takes place on the rivers during spawning in the spring, and walleye from different rivers 99 comprise a mixed-population fishery in the lake during the summer, both recreationally and for 100 subsistence (Table S1). However, recreational non-Cree fishers are only permitted to fish below the 51 st 101 parallel when they are without a Cree guide (Figure 1), fishing by Cree appears to be mostly in the south 102 (Table 2), and the genetically-distinct populations that contribute most to the mixed summer fishery are 103 those from the rivers of concern (Dupont et al. 2007). Documented catch by non-Cree fishers without a 104 guide has not increased between 1997 and 2015 (Table S1), but we do not have data on direct or latent 105 mortalities due to local fishing derbies. In addition, the human population and the number of households 106 in Mistissini almost doubled between 1997 and 2016 (Table S1). Cumulatively, this information indicates 107 an increase in fishing pressure in the southern rivers. 108 Using tissue samples and body size measurements collected in 2002/03 (Dupont et al. 2007) 109 ("historical") and between 2015-2017 ("contemporary"), we tested the general hypothesis that harvesting 110 over a period of 1-2.5 generations (based on ages of spawners in the southern rivers of Mistassini Lake 111 (supplementary data, Dupont et al. 2007)) was sufficient to cause coupled phenotypic and genetic changes 112 in wild walleye populations. Specifically, we predicted that, in association with recent, increased fishing 113 effort in Mistassini Lake, the following should be evident within the southern, harvested rivers but not in 114 the northern river with limited harvesting, when comparing contemporary versus historic samples. 1) 115 Reduced body size (total length and mass). 2) Reduced size-at-age. 3) Changes to population structure 116 such as collapsing/homogenization of between-river population structure. 4) Reductions in genetic 117 diversity and effective population size. 5) Signatures of selection, with putatively selected loci related to 118 growth, body size and/or maturation. 6) Greater reductions in body size, size-at-age and stronger 119 signatures of selection in females than in males, as a sexually dimorphic species with larger females than 120 males. As one of the relatively few studies incorporating genomic and phenotypic data in wild 121 populations to date, and the first to show rapid genetic change in a long-lived species, this study could be 122 used to inform population genomics parameters and monitoring practices for the sustainable harvest and 123 management of other similar long-lived species. 124 125

Materials and Methods 126
Fishing pressure, Traditional Ecological Knowledge 127 Currently, there is no mechanism in place for indigenous fisheries to report the number of fish 128 caught in Mistassini Lake. Thus, to establish trends in fishing pressure, fish abundance and body size, we 129 conducted semi-directed interviews as in Fraser et al. (2006) during February and July of 2018 with 17 130 elders and fisherman (30 -79 years of age, with 13 respondents > 40 years) ( Table 2). Answers were not 131 used for questions where respondents explicitly stated a lack of knowledge as per Gagnon and Berteaux 132 (2009), and the frequency of respondents for a given answer has been provided, using the total number of 133 respondents for that question as the denominator. In addition, we obtained census numbers for all people in 134 the community close to the lake and the number of fish caught by non-Cree fishers for a subset of years 135 (Table S1). Ethics certificate 30008247. 136 Rivers included in the study were Chalifour, Icon and Perch in the south and Takwa in the north. 137 (contemporary) by us (see Table 1 for sample sizes). Sampling was collaborative with subsistence fishers 145 for 2015 -2017. Walleye were captured via angling using the same lures and a combination of boats and 146 shore fishing, from the same locations within rivers, for both historical and contemporary sampling (Table  147 1). Catch-per-unit-effort was not available for historic samples or collaborative sampling and is therefore 148 not included here for contemporary sampling. After capture, fish were immediately placed in freshwater 149 baths with aerators. From each walleye, we collected total and fork length (TL ± 1 mm), wet mass (± 50 g), 150 sex (M, F, U (unknown, either spawned out or premature)) and a tissue sample for genetics; otoliths were 151 collected from a random subsample. Live walleye were returned to the water near the location of capture. 152 Opercular bones were collected for aging for historic samples (Table 1) To determine the best model, we used backward step-wise model selection and AIC (Akaike 1974). 173 Significance was detected at an alpha of 0.05 and all multi-comparison P-values were adjusted using the 174 false discovery rate (FDR) method for 64 planned contrasts (Tables S3 and S5 2003 within rivers (Dupont et al. 2007), they were combined for all analyses, and denoted as 2002/03. In 179 addition, in 2017 we were unable to collect any female walleye from Perch River, nor a sufficient number 180 of female walleye in Icon River to be able to use them in length/mass models (see Table 1  To test our prediction of reduction in size-at-age in the southern populations relative to the 185 northern one through time, we used a Bayesian hierarchical regression model. We used Bayesian as 186 opposed to frequentist modelling to account for possible bias due to sampling gear, as well as small and 187 variable sample sizes for aging structures across rivers. In addition, only total length was modeled 188 because mass data were not included in the historical dataset containing age information. We assumed 189 walleye total length (TL) was normally distributed such that: respectively. Mean total length for the ith walleye was then modelled using linear regression: 195 196 = 0 + 1 + 2 + 3 + 4 + 5 × 197

198
We used vague normal priors for all coefficients and modelled hyperpriors for age and sex by river. 199 Location (southern rivers vs northern river), and history (contemporary vs historic samples) were coded 200 as categorical variables. 201 The Bayesian model was run using JAGS version 4.3.0 (Plummer 2003) in R, using rjags and 202 run.jags (Denwood 2016). We described the posterior distribution for the model using four MCMC 203 chains. Starting parameter values for each chain were jittered. Each chain took 20,000 samples of the 204 posterior, thinned at a rate of 50. The adaption period was 1000 iterations, and a burn-in rate of 50% was 205 used, for a total chain length of 2,050,00. We evaluated MCMC chain convergence by visual inspection 206 of trace-plots to assess mixing. Additionally, we ensured that each parameter had effective samples sizes 207 >1,000 and that they passed the Gelman-Rubin diagnostic test with potential scale reduction factors 208 (PSRF) <1.1 suggesting convergence on a common posterior mode (Gelman et al. 2013). 209 210 Sequencing 211 DNA was extracted using a modified Qiagen blood and tissue kit protocol (Qiagen Inc., Valencia, 212 CA) (see Table 1 for sample sizes) and was sequenced using individual-based genotyping-by-sequencing 213 (GBS). Libraries for Ion Proton GBS were prepared using the procedure described by Masher et al. 214 (2013) at IBIS, Université Laval, Québec, Canada, with modifications described in Abed et al. (2018). 215 Libraries were prepared for sequencing using an Ion CHEF, Hi-Q reagents and P1 V3 chips 216 (ThermoFisher), and the sequencing was performed for 300 flows. Enzymes used to cleave the DNA were 217 rare cutter pst1 and frequent cutter msp1. 218 Single-nucleotide polymorphisms (SNPs) were determined from raw sequence reads using the 219 stacks pipeline v1.45 (Catchen et al. 2013), and de novo sequence alignment, on the supercomputer 220 Guillimin from McGill University, managed by Calcul Québec and Compute Canada. Pre-processing of 221 fastq files was completed using fastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to 222 assess read-quality before and after using cutadapt (Martin 2011) to trim any remaining adapters and 223 remove sequences <50bp. Our stacks parameter optimization method was similar to Mastretta-Yanes et 224 al. (2015), but we did not estimate error rate because we did not have enough positive controls to do so. 225 Final stacks parameters included default settings with the following custom options: within 226 process_radtags, 80bp trim-length; ustacks, SNP model, alpha = 0.1 for SNP calls, -m = 7; cstacks, -n = 3; 227 rxstacks, log-likelihood cut-off = -30 for SNP calls; populations, log-likelihood cut-off of -30 for SNP 228 calls, choose single SNP, maf = 0.01, -r = 0.8. We ran populations twice. First, we used the parameters 229 listed here, and generated a blacklist of loci consisting of loci with FIS < -0.3. We then re-ran populations 230 with the same parameters listed here using this blacklist, with each -p=6/8, 5/7 and 4/6. No negative 231 controls produced stacks, and all positive controls assigned to the correct populations using Discriminant 232 Analysis of Principal Components (DAPC) analysis in Adegenet (Jombart 2008;Jombart et al. 2010). 233 After quality trimming and filtering, an average of 8457 high-quality SNPs were used to estimate 234 population structure, genetic diversity and effective population size (Ne). 235 236

Population structure 237
To test our prediction that harvesting in the southern populations would change genetic 238 population structure, potentially homogenizing structure between southern rivers over time, we assessed 239 population structure using DAPC, ADMIXTURE (Alexander et al. 2009), and genetic distance (FST) 240 (using GenoDive and 999 permutations, Meirmans & Van Tienderen 2004;Weir & Cockerham 1984). 241 The optimal number of principal components (PCs) to retain for DAPC was determined using the xval 242 procedure, using n/3 (recommended by the manual) as the maximum number of PCs allowable, and 500 243 replicates. The population grouping that best fit the data was assessed using Bayesian Information 244 Criterion (BIC) for DAPC, while for ADMXTURE analysis, we used cross validation (CV) and 500 245 bootstrap replications. Both analyses were completed at least four times using different bootstrap values 246 and numbers of replicates to ensure results were stable. Based on DAPC, it was clear that 28 individuals, 247 primarily sampled in Icon and Chalifour Rivers, were from different unsampled genetic source 248 populations. We removed these individuals, re-ran the populations module of stacks, and conducted all 249 subsequent analyses using this reduced dataset. 250 251

Removing loci potentially under selection 252
Global outlier loci (loci putatively under selection) were detected using PCAdapt (Luu et al. 253 2017), using the scree plot method to determine the best number of PCs (K) to retain (Jackson 1993), and 254 Mahalanobis distance with alpha = 0.1 to determine outliers. PCAdapt does not use pre-defined 255 population structure, but instead ascertains structure based on PCA. The program detects outliers based on 256 how they relate to the structure of populations on the PCA (i.e., the distance between a point and a 257 distribution). After removing the outlier loci from the dataset, the effect of linkage disequilibrium (LD) on 258 population structure was assessed by finding markers that were in LD (r 2 = 0.7) using plink v1.9 (Chang 259 et al. 2015). We removed these markers (n = 507) and re-analyzed population structure with DAPC. Since 260 linked loci had no effect on structure, they were retained for all subsequent analyses. Genetic diversity, 261 FST and Ne analyses were completed both including and excluding global outlier lociwhile there was 262 little difference between results including or excluding outlier loci for genetic diversity and Ne, the 263 magnitude of FST was greater with outlier loci included, and the conclusion changed in one case for FST. 264 Thus, we have included only results for neutral loci here for these two metrics. 265 266

Genetic diversity and Ne 267
To test the prediction that genetic diversity and Ne were reduced over time in southern 268 populations, genetic diversity estimates were obtained using the populations model of stacks, and per-269 generation Ne (5-7 year generation time) was estimated using the linkage disequilibrium method in 270 NeEstimator v2.01 (Do et al. 2014). Ne estimates and confidence intervals were corrected for linkage by 271 correcting for chromosome number according to Waples et al. (2016). 272

Signatures of selection 274
To test the prediction that signatures of selection would be most evident between timepoints 275 within southern and not northern river(s), and that putatively selected loci would be associated with 276 relevant biological processes, analyses to determine outlier loci were conducted with PCAdapt using the 277 method described above. Analyses were conducted both for sexes combined and separately: (i) for all 278 populations combined, (ii) for southern rivers only, and (iii) within each population. For all analyses, 279 except when all populations and years were included, the populations module of stacks was re-run 280 including only the populations and/or sexes that were being contrasted, specifying that loci had to be 281 present in both populations (see Table S7 for sample sizes and numbers of loci in each analysis). For sex-282 based analysis of the full dataset including all rivers and years, loci were required to be present in 6 of 8 283

populations. 284
To determine possible functions for outlier loci, for all within-river historic-contemporary 285 contrasts for which there were outlier loci, FASTA files were blasted, mapped and annotated using 286 blast2go (Götz et al. 2008). Default parameters were used with the following custom choices: proprietary 287 cloudblast, fast-blast, UniProtKB/Swiss-Prot (Swissprot_v5) database, blast e-value 1.0E-5, 10 blast hits, 288 filtered GO by taxonomy taking only matches to animals (Metazoa).  (Table S1). In addition, while there is currently no data on the number of fishers in the community nor the 300 proportion of the population that fishes, the majority of fishers (16/17) fished in the southern area of 301 concern in the lake, more than double than in all other areas of the lake except Takwa River (where 9/17 302 fishers fished) ( Table 2). 303 304 Body size at spawning and size-at-age 305

Body size 306
Our prediction that body size would be reduced in southern populations within a 1-2.5-generation 307 period was supported. The main effects, year, river and sex all had a significant effect on total length and 308 mass (Table 3), and the best-fit models for each TL and mass each contained a three-way interaction 309 between year, sex and river. AIC was >10 better for the full model for TL and >6 better for mass. Both  Figure 2); lack of significance may relate to low female sample size in this river ( Table  316 1). Indeed, the trend for Chalifour females was particularly evident when contrasted to Takwa (the 317 reference northern river), where mean sizes of fish were consistent across all sampling years. Finally, a 318 sex-bias for more males than females being captured at spawning sites was consistent for all sampling 319 years and for all rivers, including Takwa ( Figure S1). 320 321

Size-at-age 322
Our prediction that fish in the southern rivers would be smaller for their age through time was 323 supported, although the reduction in size was small (See Table 4 for posterior means and 95% credible 324 intervals). Overall, fish were larger in the southern than the northern river(s), over all rivers fish were 29.4 325 mm larger contemporaneously than they were historically, and males were 48.2 mm smaller than females 326 on average. Lastly, and the term that tested our hypothesis, fish in the southern rivers were 13.7 mm 327 smaller relative to fish in the north in contemporary relative to historical samples. 328 All parameter estimates passed convergence checks. Each parameter had effective samples sizes 329 >1,000 and passed the Gelman-Rubin diagnostic test with potential scale reduction factors (PSRF) <1.1, 330 suggesting convergence on a common posterior mode (Gelman et al. 2013).  (Table 5), and then merged as a single metapopulation in 2015. At a within-population 344 level, Chalifour River was differentiated between timepoints, Icon did not diverge between timepoints, 345 Perch diverged marginally between timepoints, and Takwa did not diverge between timepoints. Given 346 that k = 3 was identified as the best structure overall, subsequent genetic diversity and Ne analyses were 347 conducted using a meta-population structure for Icon-Perch, with Chalifour and Takwa identified 348 separately, but years were defined as separate populations for temporal analysis. show that all populations remained large (i.e., Ne in the thousands). 362

Signatures of selection 364
Our prediction that signatures of selection would be present between historic and contemporary 365 timepoints within southern rivers but not the northern river, with putatively selected loci related to 366 growth, body size and/or maturation, was supported. Eleven to 263 loci were outliers (0.17-2.83%) 367 (Figures 4a -b). Outliers were found in the global PCA that included all rivers and both timepoints, for Lastly, parallel outliers existed between the southern rivers ( Figure 4d); of note, more outliers were in 379 common between Icon-Perch F and Chalifour M (18 outliers) than between Icon-Perch F and Chalifour F 380

(0) or Chalifour M and F (3). 381
Sex-specific analyses provided greater resolution than with sexes combined (Figures 4c and d). 382 For example, outlier SNPs maintained population structure between timepoints in more of the southern 383 rivers in sex-specific rather than combined analyses (Table S7, Figures S2-7). In addition, there was no 384 population structure in Icon-Perch males, but there was in females (Figures 4c and S12-13). Lastly, the 385 number of outliers in common between rivers differed between the sexes (Figure 4d), though this may be 386 due in part to the number of individuals sequenced. 387 Blasts were conducted for southern (F/M combined and separately), Chalifour (F/M combined 388 and separately), and Icon-Perch F (no outliers were found for Icon-Perch M or Takwa) Rivers. Southern 389 F/M combined and M had no blast hits. Otherwise, between 2 and 6 alleles were annotated for each blast. 390 Because annotations were completed against all mapped Metazoa, at level 2 go-annotation, functional 391 annotations included many different biological processes, molecular functions and cellular components. 392 Three relevant processes indicated were growth, metabolism and developmental process (Table S8).  (Table 4), genomic change evidenced by changing genomic population structure (Figure 1, 399 Table 5) and putative signatures of selection within rivers (Figure 4), both with sexes combined and 400 separately. These changes were present in the southern rivers most-impacted by increased fishing pressure 401 by Cree and non-Cree fishers alike (Tables 2 and S1), and not in the northern river where there were 402 fewer boats and fishers. Importantly, not only is fishing pressure greatest in the south, but southern fish 403 from the affected spawning runs remain close to those spawning runs in the summer mixed population 404 fishery (Dupont et al. 2007). A genetic drift hypothesis could posit that the observed genomic changes are 405 stochastic. However, all of the observed populations had large Ne, making it unlikely that drift caused the 406 phenotypic or genetic changes observed. The difference between neutral and adaptive genomic results, 407 however, illustrates the capacity for genetically-large populations in nature to rapidly respond to harvest-408 induced selective pressures. moreover, size-at-age was reduced in the south over time (Table 4). In fact, size reductions in all southern 412 rivers were likely underestimates of the true change. Namely, 2016 monitoring was largely collaborative; 413 approximately 48% of all sampled walleye (216 of 446) were harvested and donated by fishers. Donated 414 2016 walleye were 639 g (stderr ± 21 mm) and 424 mm (stderr ± 4 mm) on average compared to 603 g 415 (stderr ± 20 g) and 410 mm (stderr ± 4 mm) from our caught and released 2016 walleye (note that 416 sampling was not collaborative in this way in 2002/03). Lastly, given that Perch River males were 417 smaller than females, it is less likely that they would be subject to size-selective harvesting. In sum, these 418 results support the idea that fishers often target larger fish, and this type of size-selective harvesting has 419 been documented to lead to the evolution of smaller body size (Heino et al. 2015;Hutchings 2005;Swain 420 et al. 2007). 421 Genomic change occurred between timepoints in the southern Rivers but not in the northern river. 422

Population structure was homogenized over time between Icon and Perch Rivers. Signatures of selection 423
were evident within southern Rivers: rivers were genetically differentiated between years (Table 5), 424 outlier loci maintained that structure (Figure 4, Table S7), and exploratory analysis revealed relevant 425 biological functions associated with a small number of those outlier loci (Table S8). In addition, although 426 the extent of parallelism in events of natural selection is variable (Oke et al. 2017), parallel outlier loci 427 were detected between timepoints in the southern rivers ( Figure 4c). However, genomic change was 428 clearly nascent. Neutral genetic diversity did not change between timepoints. Differences between the 429 preferred population structures in ADMIXTURE were small (Figure 1). FST within rivers between years 430 and between Icon and Perch 2015 were small (Table 5), and scree plots for outlier locus detection showed 431 weak structure in two cases ( Figures S9 and S12). Nonetheless, results were generally consistent when 432 sexes were analyzed together and separately; although congruent with reductions in body size between 433 2002/03 and 2015 in Perch River, genetic structure was present between timepoints in Icon-Perch females 434 but not males (Figure 4). were different historically could have moved to a different spawning location in later years (Bigrigg 440 2008), though Mistassini Lake populations generally have strong spawning site fidelity (Dupont et al. 441 2007). Another possibility is that individuals from Icon River could be using Perch River to spawn much 442 more now than historically, either replacing genotypes that have been fished out, or increasing gene flow 443 substantially (Allendorf et al. 2008). Although the observed neutral and putatively selective genomic 444 change in Icon-Perch is rapid, it is not without precedent (3 generations or less, Chebib et al. 2016;van 445 Wijk et al. 2013), and even though these are genetically large populations (Figure 3a), rapid adaptation is 446 possible via soft sweeps (Hermisson 2005;Messer & Petrov 2013). In sum, nascent genomic change 447 occurred within a 12-year period (genomic samples were 2003 and 2015) within the southern most-448 harvested rivers, which represents 1-2.5 generations maximum. 449 Our data are consistent with rapid genetic changes to population structure and genetically-based 450 phenotypes due to harvest, but alternative explanations must be explored. The observed reduction in size-451 at-age in the south was small relative to the overall change in body size between historical and 452 contemporary samples, and may be due to a difference in aging structures used (historical using opercula, 453 contemporary using otoliths) (Faust & Scholten 2017). However, ages calculated using opercula and 454 otoliths have been highly correlated in walleye (Geisler 2012), and opercula have been validated to the 455 age of 16 in walleye (94% of aged fish in Mistassini were < 16) (Faust & Scholten 2017). 456 Another alternative could be that the large body size changes are entirely plastic due to changes in 457 the environment or a habitat shift unrelated to fishing. However, climate change is expected to warm the 458 Mistassini region; as a cold oligotrophic lake, Mistassini is not ideal habitat for walleye, which prefer 459 mesotrophic lakes (Kitchell et al. 1977;Niemuth et al. 1972). Climate warming is expected to increase 460 the growing season length for walleye, and thus in the absence of fishing an increase in body size is a 461 more likely response with climate warming than a decrease. Regarding plasticity, growing degree day 462 (GDD) was shown to account for 96% of the variation in length of immature walleye over 416 463 populations in Ontario and Quebec, though variation in growth associated with food availability was also 464 evident (Neuheimer & Taggart 2007;Venturelli et al. 2010). Thus, although it is unlikely that all 465 observed changes are due to selection, smaller body size at spawning and smaller size-at-age could 466 indicate that fish are selectively growing slower (Enberg et al. 2012). 467 Under variable recruitment (Bozek et al. 2011;Hansen et al. 1998), fish captured in 2002 and468 2003 could represent distinct, strong year classes, biasing estimates of mean size and contributing to 469 temporal genetic differences. However, there was no genetic differentiation within rivers between 2002 470 and 2003 (Dupont et al. 2007), Takwa River walleye had consistent body size in all years sampled, and 471 we found a significant reduction in size-at-age in the southern rivers. 472 We estimated the selective pressure required to generate the observed changes in body size within 473 1-2.5 generations to assess whether they were biologically plausible using the breeders equation (R = h 2 S, 474 where R = response to selection, h 2 = heritability, S = selection differential). Using averages for 11-year 475 old walleye in the south for each 2002 (515 mm) and 2015 (424 mm) (a 13-year interval), R = -7 mm per 476 year. Given realistic h 2 estimates (0.3) (Law 2000;Nussle et al. 2009), if observed changes were entirely 477 due to selection, S would need to be -23 mm per year. Our Bayesian model held all variables constant 478 when assessing the size change associated with harvest, and in so doing found that the size change 479 attributed to selection was smaller than what is shown here. Thus, it is very likely that some of the 480 observed body size change was plastic and/or due to stochasticity in addition to selection pressure. 481

Conclusions and management implications 482
We have presented coupled phenotypic and genetic evidence consistent with fisheries-induced 483 genetic changes within 1-2.5 generations in wild walleye populations; with rare comprehensive evidence, 484 this study sets a precedent for the timeframe needed for investigating concerns regarding harvest-induced 485 evolution in fisheries. Furthermore, sex-specific dynamics for both body size and genomics herein 486 highlight the importance of collecting sex-specific data. 487 Our study illustrates the power of integrating life history and genomics methods for conservation 488 in order to understand the factors affecting population change (Bernatchez et al. 2017), of integrating 489 these with TEK, and of iterative population monitoring practices (Flanagan et al. 2018); i.e., this study 490 would not have been possible without the historic data. Considerations for Cree management could 491 include that observed phenotypic and genetic changes may cause reduced productivity (Allendorf et al. 492 2008;Hutchings 2005), and that genomic change is clearly nascent here. Depending on the severity of 493 harvest (which is not precisely known in this case) and life-history (Audzijonyte & Kuparinen 2016), 494 fisheries-induced changes may be reversed in 9 generations (Conover et al. 2009) or less (Feiner et al. 495 2015) if fishing is halted.  Table 1: Details of sample sizes for walleye caught in each tributary of Mistassini Lake for each sex and 513 year for each analysis, as well as whether samples were caught from a boat or on shore. boat 51(M), 22(F), 76(U) ǂTakwa and Chalifour spawning grounds are only accessible by boat; thus, while we do not have record of exactly how fishing was conducted, boat can be inferred ^Individuals with unknown sex were not used in body size models, and nor were categories with any fewer than 8 observations *We extracted DNA from 371 walleye (historical n = 173 from 2003, contemporary n = 198 from 2015), and that Icon and Perch were analyzed as a single population 515 Damage to fish after being handled 1 Fishing derby's causing many dead fish 2 Pike are after walleye eggs 2 Pollution in lake, dirty water 2 No other concerns, or they will be fine 8 520 521      contrast (i.e., the number of SNPs used was unique for each sex and population). B) Percentage of SNPs that were outliers in each pairwise contrast. C) Number of outlier SNPs associated with each PCaxis. The star ( ) denotes which PC axis separated years. Where there is no star on a bar population structure between years was not maintained by outlier loci. See Table S7 and Figures S2 -S16 for further explanations and detailed descriptions of how outliers on each PC axis maintained the observed population structure. Where no bars are shown (i.e., for Takwa and Icon-Perch m in % outlier SNPs and SNPs associated with each PCaxis) there was no population structure and thus no outlier loci between years. D) Outlier loci overlap between the southern populations.