Genotyping-by-sequencing illuminates high levels of divergence among sympatric forms of coregonines in the Laurentian Great Lakes

Effective resource management depends on our ability to partition diversity into biologically meaningful units. Recent evolutionary divergence, however, can often lead to ambiguity in morphological and genetic differentiation, complicating the delineation of valid conservation units. Such is the case with the “coregonine problem,” where recent post-glacial radiations of coregonines into lacustrine habitats resulted in the evolution of numerous species flocks, often with ambiguous taxonomy. The application of genomics methods is beginning to shed light on this problem and the evolutionary mechanisms underlying divergence in these ecologically and economically important fishes. Here, we used restriction site-associated DNA (RAD) sequencing to examine genetic diversity and differentiation among sympatric species in the Coregonus artedi complex in the Apostle Islands of Lake Superior, the largest lake in the Laurentian Great Lakes. Using 29,068 SNPs, we were not only able to clearly distinguish the three most common forms for the first time, but putative hybrids and potentially mis-identified specimens as well. Assignment rates to form with our RAD data were 93-100% with the only mis-assignments arising from putative F1 hybrids, an improvement from 62-77% using microsatellites. Estimates of pairwise differentiation (FST: 0.045-0.056) were large given the detection of hybrids, suggesting that hybridization among forms may not be successful beyond the F1 state. We also used a newly built C. artedi linkage map to look for islands of adaptive genetic divergence among forms and found widespread differentiation across the genome, a pattern indicative of long-term drift, suggesting that these forms have been reproductively isolated for a substantial amount of time. The results of this study provide valuable information that can be applied to develop well-informed management strategies and stress the importance of re-evaluating conservation units with genomic tools to ensure they accurately reflect species diversity.


Introduction
50% of loci, and (3) removing loci with a minor allele count less than 3. In addition, since all salmonids 213 including members of the C. artedi species complex have experienced a recent genome duplication, 214 putatively paralogous loci were identified with the program HDPlot (McKinney, Waples, Seeb, & Seeb, 215 2017), and any loci with heterozygosity greater than 0.55 or a read ratio deviation greater than 5 and less 216 than -5 were removed. Finally, loci on the same RAD tag may be linked so only the SNP with the highest 217 minor allele frequency on each tag was included in the final dataset. All file format conversions were 218 performed using PGDSpider v2.1.1.5 (Lischer & Excoffier, 2011). 219

Microsatellite amplification and genotyping 220
We used the methods described in Stott et al. (in press) to genotype 12 microsatellites developed for 221 coregonines (Bernatchez, 1996 193bp) size ranges, therefore Sfo8 is considered two loci. Fragment analysis was performed using a 227 Genetic Analyzer 3.0 (Life Technologies), and genotypes were assigned at each locus using GeneMapper 228 3.7 (Life Technologies). We used Genepop v4 (Rousset, 2008) to conduct exact tests for deviations from 229 Hardy-Weinberg and linkage equilibrium (α = 0.01). Three loci were removed for being out of  Weinberg Equilibrium in all three forms (Bwf1, Cocl23, CoclLav6), and no loci showed significant 231 linkage disequilibrium. 232

Genetic differentiation and diversity 233
To ensure that our putative form designations were appropriate, we used two approaches to assess genetic 234 similarity among individuals. First, we conducted a principal component analysis (PCA) in the R package 235 'adegenet' (Jombart, 2008) for both the RAD and microsatellite datasets. Next, we estimated the number 236 of ancestral populations, K, contributing to contemporary genetic clustering for the RAD dataset using the 237 program ADMIXTURE v1.3 (Alexander, Novembre, & Lange, 2009). We tested K from 1-5 with 238 ADMIXTURE's cross-validation procedure and a k-fold of 10 (parameter flag: --cv=10) to examine 239 support for each K. A Q-score of less than 70% was used to identify putative hybrids following similar 240 thresholds applied in the literature ( PCA and ADMIXTURE analysis with RAD genotypes revealed that three individuals originally 245 identified as hoyi fell within the cluster of kiyi samples and, given the discrete clustering of kiyi and the 246 strong possibility that these three hoyi samples were misidentified in the field, we removed these 247 individuals from both RAD and microsatellite datasets for all further analyses to prevent bias in estimates 248 of diversity, differentiation, and Ne. 249 We calculated a variety of summary statistics for the groups comprised of the five forms (artedi -ART, 250 hoyi -HOY, kiyi -KIY, nigripinnis -NIG, and zenithicus -ZEN) including percentage of total observed 251 alleles, observed and expected heterozygosity, and an inbreeding coefficient (GIS). Summary statistics for 252 each locus and populations were calculated using both microsatellite and RAD datasets in GenoDive 253 v2.0b23 (Meirmans & Van Tienderen, 2004). Genetic differentiation among all forms was estimated 254 across all loci in the RAD dataset with pairwise FST (Weir & Cockerham, 1984) in Genepop and tested 255 using exact tests (Goudet et al., 1996; Raymond and Rousset, 1995; alpha = 0.01) in Arlequin (Excoffier 256 & Lischer, 2010). We also calculated locus-specific overall and pairwise-FST values (Weir & Cockerham, 257 1984) in Genepop using a dataset that included the three forms with n>10 (ART, HOY, and KIY). To 258 compare genetic differentiation between RAD and microsatellite datasets, we used GenoDive to estimate 259 standardized pairwise genetic differentiation, G'ST (Hedrick, 2005), which employs an additional 260 correction for bias from sampling a limited number of populations (Meirmans & Hedrick, 2011). 261 The rate at which individuals were able to be assigned back to their form with both RAD and 262 microsatellite datasets was tested using population assignment in GenoDive for all forms with n>10. 263 Assignment was performed by calculating the home likelihood (Lh) that an individual genotype is from a 264 specific group given the allele frequencies (Paetkau, Calvert, Stirling, & Strobeck, 1995)  hybridization among forms will result in introgression and decreased differentiation, we modeled the 286 impacts of hybridization on genetic differentiation with three scenarios: two populations with a starting 287 level of differentiation equal to 1) approximately the level we observed amongst cisco forms with our 288 RAD dataset (FST ≈0.05), 2) approximately twice the amount observed (FST ≈0.10), and 3) approximately 289 four times the amount observed (FST ≈0.20). Preliminary levels of differentiation were achieved with low 290 migration rates (m) set over 1,000 generations (m=0.001, m=0.0005, and m=0.0001, respectively) and 291 followed by 5, 10, or 15 generations of one of three higher migration (i.e. hybridization) rates (m=0.01, 292 0.05, and 0.10) for every initial level of differentiation. Each unique parameter combination was run using 293 a dataset of 1,000 loci and replicated 50 times. Estimates of FST for each replicate were generated in the R 294 was set at FST ≈0.05 using the same method described above, and we allowed populations to hybridize for 300 2, 5, or 10 generations with m=0.05, which was equal to roughly the proportion of hybrids observed in 301 our wild populations between ART-HOY and HOY-KIY. Since no putative ART-KIY hybrids were 302 observed in our RAD dataset, we chose to implement a one-dimensional stepping stone model of 303 migration (Kimura & Weiss, 1964) among simulated populations. Each unique parameter set for these 304 simulations was run using 5,000 loci, a genepop file was output containing a subset of 100 randomly 305 selected individuals per population (50 females/50 males), and ADMIXTURE was used to generate Q-306 scores for each individual. Hybrids were identified using the same Q < 0.70 threshold applied above. 307

Differentiation across the genome 308
We examined genetic differentiation across the genome by pairing our data with the artedi linkage map to 10,000. Contiguous windows that contained at least two loci and exceed the 99 th percentile of the 321 distribution after 10,000 bootstrap replicates were classified as putative islands of divergence. We 322 investigated genomic differentiation using four locus-specific metrics: overall FST and pairwise FST 323 between each of the three putative forms with n>10 (ART-HOY, ART-KIY, and HOY-KIY). We then 324 plotted overall FST across the genome and constructed a bubble plot to visualize the number of significant 325 windows on each LG for each comparison. We found over 50 significant windows for each comparison 326 (see results). Conducting an in-depth investigation of all significant windows was not feasible, therefore 327 we isolated in-depth analysis to one LG (Cart21) that contained the most significant windows in the 328

dataset. 329
To investigate this highly differentiated LG, we aligned consensus sequences for all loci from Cart21 to 330 chromosome Ssa05 in Atlantic salmon (Salmo salar), which is syntenic to Cart21 in artedi (Blumstein,  conducted with BLASTN. The best alignment for each locus was retained, and all alignments had e-333 values < 1e-58. We then visualized the relationship between recombination and physical distance using 334 alignments to Ssa05 and information from the artedi linkage map to determine whether this highly 335 differentiated region is characterized by lower recombination. We also obtained annotation information 336 from the Atlantic salmon genome to determine whether genes of interest were co-located with areas of 337 high divergence. Finally, we plotted the allele frequencies of the 10 SNPs on Cart21 with the highest 338 overall FST to investigate whether these SNPs show consistent patterns of population structure. 339

Results 340
Sequencing and genotyping 341 A total of 137 individuals were RAD sequenced producing more than 455 million reads and an average of 342 3,346,457 reads per sample. After filtering, 119 individuals with representatives from all five putative 343 cisco forms in Lake Superior were genotyped at 29,068 loci (Table 1). More than half of these loci 344 (n=15,348 loci) were also placed on the linkage map. Since both RAD sequencing and microsatellite 345 amplification were performed on the same hoyi and kiyi samples, microsatellite genotypes used in our 346 analyses were restricted to the same individuals that successfully genotyped with our RAD loci. An 347 additional 30 artedi from Stockton Island were genotyped at the 9 microsatellite loci. Paired-end 348 assemblies for each locus are available from Blumstein (2019). 349

Genetic differentiation and diversity 350
PCA showed a sharp contrast in resolution between marker sets with microsatellites producing one large 351 cluster of overlapped forms across the first two principal components and the RAD dataset producing KIY clusters (Fig. 2). The sixth sample fell within the ART cluster, and like the three hoyi in the KIY 358 cluster, possibly represents a misidentified specimen. Unlike the NIG samples, which suggest the 359 possibility for a distinct cluster with the addition of more specimens, the four ZEN samples closely 360 associated with either the HOY cluster (n=1) or the KIY cluster (n=3). Low representation of both NIG 361 and ZEN in our RAD dataset reduces our ability to draw strong conclusions based on these PCA results, 362 so both groups were unaltered for estimates of diversity, inbreeding and differentiation and dropped for 363 the remaining analyses. 364 The most supported number of ancestral populations (K) estimated using the cross-validation procedure in 365 ADMIXTURE was two (Fig. 3). Examining additional Ks for significant sub-structuring among forms 366 generated results that corroborated those from the PCA. When K=2, the ART cluster splits from HOY, 367 KIY, and ZEN. Individuals in the NIG cluster exhibited mixed ancestry between the two major groups as 368 seen on PC1 of the PCA. When K=3, the major genetic ancestries differentiate the ART, HOY, and KIY 369 clusters as seen on PC2. The three putatively misidentified hoyi first noted in the PCA all had Q estimates 370 of 100% for the KIY cluster and were removed from further analyses. Additional Ks did not differentiate 371 either the NIG or ZEN cluster but begin to differentiate small subsets of individuals within groups, which 372 was likely statistical noise.  (Table 1). 375 Inbreeding coefficients were not substantially different from zero in both datasets. The largest GIS was 376 measured in ZEN from only four samples (0.279) and the rest were between -0.028-0.196. All estimates 377 of pairwise genetic differentiation among forms with n>10 (ART, HOY, KIY) with the RAD dataset were 378 significant ( Table 2). The magnitude of genetic differentiation followed similar trends observed in the 379 PCA, with ART being slightly more differentiated from deepwater forms (FST=0.049-0.056). This pattern 380 remained the same with a standardized measure of genetic differentiation in the RAD and microsatellite 381 datasets (Table 3). All pairwise comparisons in both datasets were significant with higher overall values 382 of G'ST generated using microsatellites (0.110-0.122) than SNPs (0.060-0.75). See Tables S1 and S2 for 383 locus-specific summary statistics. 384 Population self-assignment rate using the microsatellite dataset ranged from 61.9-76.7% with artedi 385 exhibiting the highest likelihood of being assigned back to the ART group (Table 1). Assignment rate 386 with the RAD dataset was 100% for the KIY group and 98.3% and 92.5% for the ART and HOY groups 387 (respectively), a result of one putative artedi assigning to HOY and two putative hoyi individuals 388 assigning to KIY. In the ADMIXTURE analysis, these three individuals exhibited relatively high Q 389 estimates for ancestry to the populations to which they were assigned (artedi, QHOY = 67.4%; hoyi, QKIY= 390 64.2% and 56.2%). In the PCA generated from the same data, these individuals were oriented between the 391 main clusters. In the microsatellite dataset with the same hoyi and kiyi samples, neither of the two 392 potentially misclassified hoyi assigned to the HOY group, with one being assigned to ART and one to 393 KIY. Assignment scores are reported in supplementary tables S3-S5. 394 Estimates of Ne with the microsatellite dataset ranged from 73 in HOY to 659 in KIY (Table 1) Waples & Do 2010). For the RAD dataset, estimates of Ne were generated with loci that were placed on 400 the linkage map and ranged from 1,701-2,126 with confidence intervals within 10% of these values. 401 We documented a relatively high level of hybridization in our empirical dataset, with ART-HOY hybrids 402 comprising 8.6% of the total number of sampled artedi and hoyi, HOY-KIY hybrids comprising 4.3% of 403 sampled hoyi and kiyi , and ART-KIY hybrids comprising 0% (Fig. 3, Table S6). Additionally, most of 404 these appear to be F1 hybrids with similar contributions from two genetic groups (Fig. 3). The fact that 405 we observed frequent hybridization coupled with relatively high genetic differentiation was puzzling and 406 prompted us to conduct two types of simulations to investigate how hybridization (i.e. migration) can 407 influence genetic differentiation. Simulated hybridization over a 15-generation period resulted in declines 408 in genetic differentiation among populations for all tested levels of migration (Fig. 4). When the 409 migration rate was set to 5 or 10%, levels of differentiation more than halved after only five generations 410 for all initial levels of FST. A migration rate of 1% resulted in steadily declining genetic differentiation but 411 only began to approach a halved FST after 15 generations. 412 Simulations of stepping stone migration between three populations and subsequent ADMIXTURE 413 analysis resulted in an overall pattern very similar to that observed in our empirical data (Figs 3, S1; 414 Table S6). The major goal of these simulations was to determine whether individuals that we observed in LG=10, SD=10, range=0-36).
LG Cart21 displayed the highest number of differentiated windows, 431 prompting us to conduct an in-depth investigation of this LG to investigate patterns and potential drivers 432 of divergence (Fig. 6). 433 Cart21 contained 36 significantly differentiated windows; with the largest number of windows found for 434 the ART-KIY comparison (18), followed by the overall FST comparison (13), ART-HOY comparison (4), 435 and HOY-KIY comparison (1). We were able to place 351 loci on Cart21, and 12 of these loci displayed 436 overall FST values > 0.3 (Fig. 6a). The largest cluster of high-FST loci was found between 0-10 cM on the 437 linkage map. This region appears to be characterized by relatively low recombination, as loci found in the 438 first 10 cM of Cart21 span about 25 megabases of the Atlantic salmon genome (Fig. 6b). Alignments to 439 the Atlantic salmon genome were possible for 151 loci on Cart21, and these alignments revealed that the 440 highest FST loci were found between positions 15 million and 25 million on Ssa05 (Fig. 6c). Some of 441 these loci were found in genes with functions that include cell signaling and membrane transport. 442 However, there are over 2,000 genes on Ssa05, making it difficult to reach any robust conclusions about 443 the functional significance of our loci. Allele frequencies at the high-FST were generally the most 444 diverged between ART and the other two forms, KIY in particular (Fig. 6d). 445

Discussion 446
The resolution of genetic structure in recently diverged species complexes has proved challenging with 447 traditional genetics methods, prompting the reevaluation with genomics of a taxonomically uncertain 448 species complex in the Laurentian Great Lakes. Using 29,068 SNPs to examine differentiation and 449 diversity of sympatric coregonines in the Apostle Islands of Lake Superior, we were able to 450 unambiguously assign individuals to the three major forms as well as identify putative F1 hybrids and that hybridization beyond the F1 state is not successful. This is further supported by the discovery of 458 widespread differentiation between forms across the genome, indicating that much of the divergence 459 observed has been driven by long term reproductive isolation and drift. 460

Hypotheses for high genetic differentiation among forms 461
Genetic differentiation of the three primary cisco forms in our study (artedi, hoyi, kiyi) was relatively high 462 compared to previous research in cisco using allozymes, mtDNA, microsatellites, AFLPs, and RAD data, 463 which has largely suggested that forms are not frequently diverged within lakes (but see Stott  Provincial Park (Ontario, Canada) using RAD data, which were much lower than the FST values of ~0.05 467 observed in our study. Unfortunately, genetic data for cisco in the Great Lakes is relatively sparse, 468 however, a previous study using mtDNA and microsatellites did not find evidence of differentiation 469 among forms (Turgeon & Bernatchez, 2003). Results from the microsatellites genotyped in the current 470 study are similar and indicate that these markers are unable to differentiate species in Lake Superior, even 471 though genetic structure was relatively high according to the RAD data. Interestingly, our estimates of 472 genetic divergence among forms more closely mirror two studies in lake whitefish and European However, it is also possible that reproductive isolation among cisco forms is more complete, reducing the 485 potential for introgression to erode divergence among forms. 486 Very little is known about the spawning biology of forms outside of artedi, although our observation of 487 relatively frequent F1 hybrids between ART-HOY and HOY-KIY suggests that there is a least some 488 overlap in reproductive timing among the three forms. It also notable that we did not observe any putative 489 hybrids between ART-KIY. These three cisco forms are encountered in different depths in Lake Superior, sperm performance in backcrosses. Our results combined with those from lake whitefish provide evidence 506 that hybrid backcrosses may experience dramatically reduced fitness in cisco. However, future research is 507 necessary to empirically test this hypothesis and investigate potential mechanisms for genomic 508 incompatibilities in hybrid backcrosses. 509

Genetic differentiation across the genome 510
Comparison of patterns of genomic divergence found in our study with previous research suggests that 511 diversification of cisco forms in the Great Lakes is likely polygenic and that these forms have been 512 isolated without gene flow for a relatively long period of time. We identified over 100 significantly 513 differentiated genomic windows in our study, and genetic differentiation among forms was consistently 514 However, as populations drift apart, high levels of divergence accumulate across the genome, making it 546 difficult to differentiate genomic islands involved in the original divergence from accumulated genetic 547 drift in populations that have been isolated for a long time period (Via, 2012). This pattern was observed 548 in sockeye salmon (Oncorhynchus nerka), where an island of divergence that is highly visible across 549 multiple drainages in Alaska was partially obscured in the Copper River, which has experienced lower 550 gene flow and therefore much more genetic drift than the other systems (Larson et al., 2019). The high 551 and relatively homogenous patterns of genetic differentiation observed in our study suggest that genomic 552 islands involved in diversification may be obscured in a similar fashion. 553

Conservation implications 554
We documented high genetic differentiation among the three major cisco forms in Lake Superior and, 555 based on this information, we suggest that separate conservation units could be constructed for each form. 556 This strategy differs from the current conservation paradigm for cisco, which recommends that the entire 557 C. artedi species complex be considered C. artedi sensu lato (translation, in the broad sense) and that 558 units of conservation should be designed to preserve environments that have facilitated the evolution of 559 different forms (i.e. lakes) rather than on forms at a larger scale (Turgeon & Bernatchez, 2003 see Stott et al., in press for a single lake example). Our results suggest that "last generation" markers may 565 be insufficient for capturing differentiation in evolutionarily young species, such as cisco, and highlight 566 the utility of genomic data for designating conservation units in these species. However, it is important to 567 note that we only surveyed animals from one lake, and a larger genotyping effort across the Great Lakes 568 is necessary to accurately inform conservation units for Great Lakes ciscoes. Additionally, our ability to 569 draw conclusions regarding the relationships of the rare forms nigripinnis and zenithicus to the three 570 major forms was limited by low sample size, therefore, genotyping additional contemporary and/or 571 historic specimens may help resolve the placement these forms. 572 The potential of genomic data to revolutionize construction of conservation units has been frequently 573 discussed (reviewed in Allendorf et al., 2010;Funk et al., 2012), and many studies have found that 574 genomic data provides increased resolution for delineating population structure compared to last 575 generation markers, such as microsatellites (Hodel et  suggest that conservation units constructed based on data from last generation markers that found either 582 low or no significant structure among groups of animals with different phenotypes could potentially be re-583 evaluated with genomic tools to ensure within-species diversity is being adequately conserved. 584

Conclusions and future directions 585
Our study provides the first evidence that cisco forms within the Great Lakes are genetically 586 differentiated. We documented high genetic differentiation among the three major forms in Lake 587 Superior, and highly differentiated markers were distributed across the genome, with islands of 588 divergence found on nearly every linkage group. Additionally, we identified putative F1 hybrids but no 589 hybrid backcrosses, suggesting that fitness breakdown of backcrosses may aid in maintaining 590 differentiation of these forms. The results of this study provide the foundation for a new understanding of 591 the ecology and evolution of the C. artedi species complex within the Great Lakes. The ability to 592 differentiate forms with genomics provides researchers with a powerful tool for ground truthing 593 morphological phenotypes and identifying cisco species at any life stage. In particular, the ability to 594 identify larval ciscoes will allow researchers to estimate recruitment, which is vital for management and 595 conservation, and will also significantly improve our understanding of early life history characteristics 596 and reproductive dynamics in this species. Finally, our results suggest that some management units 597 created using last generation markers may not adequately reflect species diversity and could be re-598 evaluated with genomic data. 599 Tables 600 Table 1. Sample statistics, diversity, and effective population size estimates. N is the number of individuals 601 successfully genotyped, A is the percentage of total sampled alleles found in each group, Ho/He are 602 observed/expected heterozygosity, GIS is the inbreeding coefficient, Assignment is the percentage of individuals 603 that were correctly assigned to their population of origin in a leave-one-out test, and Ne is effective population size 604 calculated using the LDNE method and reported with 95% confidence intervals. 605 606 *Three genotyped hoyi samples are suspected to be misidentified kiyi from the RAD-based PCA and ADMIXTURE 607 analysis and were removed to prevent bias in estimates of diversity and Ne. 608 **The one artedi that assigned to HOY and the two hoyi that assigned to KIY appear to be hybrids. Outside of these 609 putative hybrids, assignment to the ART and HOY groups was 100%.    Table 1. 674 Supplementary data: 675 Table S1. Summary statistics for 26,789 loci genotyped with RAD sequencing. The "LG" and "cM" 676 columns denote map location, columns ending in "AF" denote population allele frequencies, "HO" is 677 observed heterozygosity, HE is expected heterozygosity, the "Sequence P1 column" is the sequence from 678 the P1 read for each RAD tag, and the "Sequence PE" is the sequence obtained from paired-end 679 assemblies. Form abbreviations are identical to those in Table 1. 680 Table S2. Summary statistics for 9 microsatellite loci included in this study. Abbreviations are total 681 number of alleles per locus (A), effective number of alleles per locus (Neff), observed heterozygosity 682 (HO), and expected heterozygosity (HE). FST (Weir and Cockerham, 1984) was estimated in Genepop and 683 all other statistics were estimated in Genodive. 684 Table S3. Initial assignment scores from Genodive for all sampled individuals using RAD data. 685 Table S4. Assignment scores from Genodive when tests were limited to the three major forms using RAD 686 data. 687 Table S5. Assignment scores from Genodive for the three major forms using microsatellite data. 688 Table S6. Type, number, and proportion of hybrids observed in empirical data and simulated populations. 689 Data simulated in EASYPOP were comprised of 5,000 biallelic loci from three populations (P1-P3) that 690 experienced stepping-stone migration (m=0.05) for 2, 5, or 10 generations (G2, G5, G10) after 1,000 691 generations of low migration (m=0.001) to approximate an FST similar to that observed between our three stepping-stone migration. Simulations were run with random mating of 1,000 females and 1,000 males in 695 each population using 5,000 biallelic loci, and preliminary conditions that produced a similar level of 696 differentiation observed in our RAD dataset among forms (FST≈0.05; 1,000 generations with an island 697 migration rate of 0.001). Each ADMIXTURE plot represents a random subset of 50 males and 50 females 698 from each population. 699 Fig S2. Overall FST for all loci that could be placed on the linkage map. Each dot is a marker and red lines 700 indicate genomic windows that were significantly differentiated from the rest of the genome according to 701 kernel smoothing analysis. 702 Each dot is a marker and red lines indicate genomic windows that were significantly differentiated from 704 the rest of the genome according to kernel smoothing analysis. Form abbreviations are in Table 1. dot is a marker and red lines indicate genomic windows that were significantly differentiated from the rest 707 of the genome according to kernel smoothing analysis. Form abbreviations are in Table 1. dot is a marker and red lines indicate genomic windows that were significantly differentiated from the rest 710 of the genome according to kernel smoothing analysis. Form abbreviations are in Table 1.