Transcriptomic analyses of juvenile Striped Bass (Morone saxatilis) exposed to chronic and acute temperature change.

Striped Bass are economically important, migratory fishes, which occur across a wide range of latitudes. Given their wide-ranging nature, Striped Bass can cope with a broad range of environmental temperatures, yet the mechanisms underlying this ability have not been thoroughly described. Heat shock proteins (HSPs) are inducible molecular chaperones, which help mitigate protein damage resulting from increased temperatures. The importance of HSPs has been demonstrated in a number of fish species, but their role in Striped Bass is poorly understood. This study characterizes changes in gene expression in juvenile Striped Bass, following acute and chronic temperature change. Fish were acclimated to one of three temperatures (15, 25 or 30°C) and sampled at one of two treatments (control or after CTmax), following which we assessed differential gene expression and gene ontology in muscle. It is clear from our differential expression analyses that acclimation to warm temperatures elicits more robust changes to gene expression, compared to acute temperature increases. Our differential expression analyses also revealed induction of many different heat shock proteins, including hsp70, hsp90, hsp40 and other small HSPs, after both acute and chronic temperature increase in white muscle. Furthermore, the most consistent gene ontology pattern that emerged following both acclimation and CTmax was upregulation of transcripts involved in "protein folding", which also include heat shock proteins. Gene ontology analyses also suggest changes to other processes after acclimation, including decreased growth pathways and changes to DNA methylation. Overall, these data suggest that HSPs likely play a major role in the Striped Bass's ability to tolerate warm waters.

Reproductive ecology of muskellunge (Esox masquinongy), an introduced predator, in the lower Wolastoq/Saint John River, New Brunswick, Canada.

Introduced predators can have harmful top-down effects on their newly colonized system through competition with and direct predation on native species. Following an initial introduction of muskellunge in Lac Frontière, Québec in the 1970s at the headwaters of the Wolastoq/Saint John River, the species rapidly migrated downstream, expanding its range by ~500 km over ~20 years. Despite this expansive colonization and concern over possible threats to native species, little is known about the basic ecology of muskellunge in this system. The last downstream barrier is the hydroelectric facility, Mactaquac Generating Station (MGS), 150 km upstream of the sea. While there are no downstream fish passage facilities at MGS, adult muskellunge have been recorded downstream. In this study, muskellunge (n = 23) were surgically tagged with very-high-frequency (VHF) radio or combined acoustic radio telemetry (CART) tags and tracked over two spawning seasons. We sought to determine if there was a reproducing population downstream of MGS and tracked Tagged muskellunge over two spawning seasons. We tracked fish to locate and confirm spawning sites, and followed up with egg and/or juvenile sampling surveys. Tagged muskellunge (90%) moved upstream towards the MGS during the spawning period in each year (2016 and 2017), where they remained throughout the entire spawning period. No spawning or nursery sites were confirmed near MGS, but in 2016 three distinct spawning locations and six distinct nursery sites were confirmed 10-12 km downstream amongst a chain of flooded islands. In 2016, eggs, sac-fry and juveniles were collected and confirmed as muskellunge by genetic sequencing, providing the first empirical observation of successful spawning downstream of MGS.

Comparing mixed models and Random Forest association tests using naturalgwas and a Striped Bass SNP data set.

In this study, we used the phenotype simulation package naturalgwas to test the performance of Zhao's Random Forest method in comparison to an uncorrected Random Forest test, latent factor mixed models (LFMM), genome-wide efficient mixed models (GEMMA), and confounder adjusted linear regression (CATE). We created 400 sets of phenotypes, corresponding to five effect sizes and two, five, 15, or 30 causal loci, simulated from two empirical data sets containing SNPs from Striped Bass representing three and 13 populations. All association methods were evaluated for their ability to detect genotype-phenotype associations based on power, false discovery rates, and number of false positives. Genomic inflation was highest for uncorrected Random Forest and LFMM tests and lowest for Gemma and Zhao's Random Forest. All association tests had similar power to detect causal loci, and Zhao's Random Forest had the lowest false discovery rate in all scenarios. To measure the performance of association tests in small data sets with few loci surrounding a causal gene we also ran analyses again after removing causal loci from each data set. All association tests were only able to find true positives, defined as loci located within 30 kbp of a causal locus, in 3%-18% of simulations. In contrast, at least one false positive was found in 17%-44% of simulations. Zhao's Random Forest again identified the fewest false positives of all association tests studied. The ability to test the power of association tests for individual empirical data sets can be an extremely useful first step when designing a GWAS study.

Considering Fish as Recipients of Ecosystem Services Provides a Framework to Formally Link Baseline, Development, and Post-operational Monitoring Programs and Improve Aquatic Impact Assessments for Large Scale Developments.

In most countries, major development projects must satisfy an Environmental Impact Assessment (EIA) process that considers positive and negative aspects to determine if it meets environmental standards and appropriately mitigates or offsets negative impacts on the values being considered. The benefits of before-after-control-impact monitoring designs have been widely known for more than 30 years, but most development assessments fail to effectively link pre- and post-development monitoring in a meaningful way. Fish are a common component of EIA evaluation for both socioeconomic and scientific reasons. The Ecosystem Services (ES) concept was developed to describe the ecosystem attributes that benefit humans, and it offers the opportunity to develop a framework for EIA that is centred around the needs of and benefits from fish. Focusing an environmental monitoring framework on the critical needs of fish could serve to better align risk, development, and monitoring assessment processes. We define the ES that fish provide in the context of two common ES frameworks. To allow for linkages between environmental assessment and the ES concept, we describe critical ecosystem functions from a fish perspective to highlight potential monitoring targets that relate to fish abundance, diversity, health, and habitat. Finally, we suggest how this framing of a monitoring process can be used to better align aquatic monitoring programs across pre-development, development, and post-operational monitoring programs.

Correction: The sockeye salmon genome, transcriptome, and analyses identifying population defining regions of the genome.

[This corrects the article DOI: 10.1371/journal.pone.0240935.].

Ancient DNA SNP-panel data suggests stability in bluefin tuna genetic diversity despite centuries of fluctuating catches in the eastern Atlantic and Mediterranean.

Atlantic bluefin tuna (Thunnus thynnus; BFT) abundance was depleted in the late 20th and early 21st century due to overfishing. Historical catch records further indicate that the abundance of BFT in the Mediterranean has been fluctuating since at least the 16th century. Here we build upon previous work on ancient DNA of BFT in the Mediterranean by comparing contemporary (2009-2012) specimens with archival (1911-1926) and archaeological (2nd century BCE-15th century CE) specimens that represent population states prior to these two major periods of exploitation, respectively. We successfully genotyped and analysed 259 contemporary and 123 historical (91 archival and 32 archaeological) specimens at 92 SNP loci that were selected for their ability to differentiate contemporary populations or their association with core biological functions. We found no evidence of genetic bottlenecks, inbreeding or population restructuring between temporal sample groups that might explain what has driven catch fluctuations since the 16th century. We also detected a putative adaptive response, involving the cytoskeletal protein synemin which may be related to muscle stress. However, these results require further investigation with more extensive genome-wide data to rule out demographic changes due to overfishing, and other natural and anthropogenic factors, in addition to elucidating the adaptive drivers related to these.

Life-stage-dependent supergene haplotype frequencies and metapopulation neutral genetic patterns of Atlantic cod, Gadus morhua, from Canada's Northern cod stock region and adjacent areas.

Among highly migratory fish species, nursery areas occupied by juveniles often differ from adult habitats. To better understand the spatial dynamics of Canada's Northern cod stock, juveniles caught off the east coast of Newfoundland and Labrador were compared to adults from the same region as well as individuals from other areas in Atlantic Canada using double-digest restriction site-associated DNA sequencing-derived single nucleotide polymorphisms. A reduced proportion of homozygotes with a chromosomal inversion located in linkage group 1 (LG1) was detected between juvenile and adult samples in the Northern cod stock region, potentially indicating age-dependent habitat use or ontogenetic selection for attributes associated with the many genes located in LG1. No selectively neutral genetic differences were found between samples from the Northern cod stock; nevertheless, significant differences were found between some of these samples and cod collected from St. Pierre Bank, Bay of Fundy, Browns Bank and the southern Scotian Shelf. Clustering analysis of variants at neutral loci provided evidence for three major genetic units: (a) the Newfoundland Atlantic Coast, (b) eastern and southern Gulf of St. Lawrence and Burgeo Bank and (c) the Bay of Fundy, Browns Bank and southern Scotian Shelf. Both adaptive and neutral population structure within the Northern cod stock should be considered by managers to promote demographic rebuilding of the stock.

The sockeye salmon genome, transcriptome, and analyses identifying population defining regions of the genome.

Sockeye salmon (Oncorhynchus nerka) is a commercially and culturally important species to the people that live along the northern Pacific Ocean coast. There are two main sockeye salmon ecotypes-the ocean-going (anadromous) ecotype and the fresh-water ecotype known as kokanee. The goal of this study was to better understand the population structure of sockeye salmon and identify possible genomic differences among populations and between the two ecotypes. In pursuit of this goal, we generated the first reference sockeye salmon genome assembly and an RNA-seq transcriptome data set to better annotate features of the assembly. Resequenced whole-genomes of 140 sockeye salmon and kokanee were analyzed to understand population structure and identify genomic differences between ecotypes. Three distinct geographic and genetic groups were identified from analyses of the resequencing data. Nucleotide variants in an immunoglobulin heavy chain variable gene cluster on chromosome 26 were found to differentiate the northwestern group from the southern and upper Columbia River groups. Several candidate genes were found to be associated with the kokanee ecotype. Many of these genes were related to ammonia tolerance or vision. Finally, the sex chromosomes of this species were better characterized, and an alternative sex-determination mechanism was identified in a subset of upper Columbia River kokanee.

Genomic population structure of Striped Bass (Morone saxatilis) from the Gulf of St. Lawrence to Cape Fear River.

Striped Bass, Morone saxatilis (Walbaum, 1792), is an anadromous fish species that supports fisheries throughout North America and is native to the North American Atlantic Coast. Due to long coastal migrations that span multiple jurisdictions, a detailed understanding of population genomics is required to untangle demographic patterns, understand local adaptation, and characterize population movements. This study used 1,256 single nucleotide polymorphism (SNP) loci to investigate genetic structure of 477 Striped Bass sampled from 15 locations spanning the North American Atlantic coast from the Gulf of St. Lawrence, Canada, to the Cape Fear River, United States. We found striking differences in neutral divergence among Canadian sites, which were isolated from each other and US populations, compared with US populations that were much less isolated. Our SNP dataset was able to assign 99% of Striped Bass back to six reporting groups, a 39% improvement over previous genetic markers. Using this method, we found (a) evidence of admixture within Saint John River, indicating that migrants from the United States and from Shubenacadie River occasionally spawn in the Saint John River; (b) Striped Bass collected in the Mira River, Cape Breton, Canada, were found to be of both Miramichi River and US origin; (c) juveniles in the newly restored Kennebec River population had small and nonsignificant differences from the Hudson River; and (d) tributaries within the Chesapeake Bay showed a mixture of homogeny and small differences among each other. This study introduces new hypotheses about the dynamic zoogeography of Striped Bass at its northern range and has important implications for the local and international management of this species.

The effects of taxonomy, diet, and ecology on the microbiota of riverine macroinvertebrates.

Freshwater macroinvertebrates play key ecological roles in riverine food webs, such as the transfer of nutrients to consumers and decomposition of organic matter. Although local habitat quality drives macroinvertebrate diversity and abundance, little is known about their microbiota. In most animals, the microbiota provides benefits, such as increasing the rate at which nutrients are metabolized, facilitating immune system development, and defending against pathogenic attack. Our objectives were to identify the bacteria within aquatic invertebrates and determine whether their composition varied with taxonomy, habitat, diet, and time of sample collection. In 2016 and 2017, we collected 264 aquatic invertebrates from the mainstem Saint John (Wolastoq) River in New Brunswick, Canada, representing 15 orders. We then amplified the V3-V4 hypervariable region of the 16S rRNA gene within each individual, which revealed nearly 20,000 bacterial operational taxonomic units (OTUs). The microbiota across all aquatic invertebrates were dominated by Proteobacteria (69.25% of the total sequence reads), but they differed significantly in beta diversity, both among host invertebrate taxa (genus-, family-, and order-levels) and temporally. In contrast to previous work, we observed no microbiota differences among functional feeding groups or traditional feeding habits, and neither water velocity nor microhabitat type structured microbiota variability. Our findings suggest that host invertebrate taxonomy was the most important factor in modulating the composition of the microbiota, likely through a combination of vertical and horizontal bacterial transmission, and evolutionary processes. This is one of the most comprehensive studies of freshwater invertebrate microbiota to date, and it underscores the need for future studies of invertebrate microbiota evolution and linkages to environmental bacteria and physico-chemical conditions.

Mercury bioaccumulation in aquatic biota along a salinity gradient in the Saint John River estuary.

Although estuaries are critical habitats for many aquatic species, the spatial trends of toxic methylmercury (MeHg) in biota from fresh to marine waters are poorly understood. Our objective was to determine if MeHg concentrations in biota changed along a salinity gradient in an estuary. Fourspine Stickleback (Apeltes quadracus), invertebrates (snails, amphipods, and chironomids), sediments, and water were collected from ten sites along the Saint John River estuary, New Brunswick, Canada in 2015 and 2016, with salinities ranging from 0.06 to 6.96. Total mercury (proxy for MeHg) was measured in whole fish and MeHg was measured in a subset of fish, pooled invertebrates, sediments, and water. Stable sulfur (δ34S), carbon (δ13C), and nitrogen (δ15N) isotope values were measured to assess energy sources (S, C) and relative trophic level (N). There were increases in biotic δ13C and δ34S from fresh to more saline sites and these measures were correlated with salinity. Though aqueous MeHg was higher at the freshwater than more saline sites, only chironomid MeHg increased significantly with salinity. In the Saint John River estuary, there was little evidence that MeHg and its associated risks increased along a salinity gradient.

RAD-Seq Reveals Patterns of Additive Polygenic Variation Caused by Spatially-Varying Selection in the American Eel (Anguilla rostrata).

The American Eel (Anguilla rostrata) has an exceptional life cycle characterized by panmictic reproduction at the species scale, random dispersal, and selection in a highly heterogeneous habitat extending from subtropical to subarctic latitudes. The genetic consequences of spatially-varying selection in this species have been investigated for decades, revealing subtle clines in allele frequency at a few loci that contrast with complete panmixia on the vast majority of the genome. Because reproduction homogenizes allele frequencies every generation, sampling size, and genomic coverage are critical to reach sufficient power to detect selected loci in this context. Here, we used a total of 710 individuals from 12 sites and 12,098 high-quality single nucleotide polymorphisms to re-evaluate the extent to which local selection affects the spatial distribution of genetic diversity in this species. We used environmental association methods to identify markers under spatially-varying selection, which indicated that selection affects ∼1.5% of the genome. We then evaluated the extent to which candidate markers collectively vary with environmental factors using additive polygenic scores. We found significant correlations between polygenic scores and latitude, longitude and temperature which are consistent with polygenic selection acting against maladapted genotypes in different habitats occupied by eels throughout their range of distribution. Gene functions associated with outlier markers were significantly enriched for the insulin signaling pathway, indicating that the trade-offs inherent to occupying such a large distribution range involve the regulation of metabolism. Overall, this study highlights the potential of the additive polygenic scores approach in detecting selective effects in a complex environment.

The evolution of the major histocompatibility complex in upstream versus downstream river populations of the longnose dace.

Populations in upstream versus downstream river locations can be exposed to vastly different environmental and ecological conditions and can thus harbor different genetic resources due to selection and neutral processes. An interesting question is how upstream-downstream directionality in rivers affects the evolution of immune response genes. We used next-generation amplicon sequencing to identify eight alleles of the major histocompatibility complex (MHC) class II ß exon 2 in the cyprinid longnose dace (Rhinichthys cataractae) from three rivers in Alberta, upstream and downstream of municipal and agricultural areas along contaminant gradients. We used these data to test for directional and balancing selection on the MHC. We also genotyped microsatellite loci to examine neutral population processes in this system. We found evidence for balancing selection on the MHC in the form of increased nonsynonymous variation relative to neutral expectations, and selection occurred at more amino acid residues upstream than downstream in two rivers. We found this pattern despite no population structure or isolation by distance, based on microsatellite data, at these sites. Overall, our results suggest that MHC evolution is driven by upstream-downstream directionality in fish inhabiting this system.

Draft genome of the American Eel (Anguilla rostrata).

Freshwater eels (Anguilla sp.) have large economic, cultural, ecological and aesthetic importance worldwide, but they suffered more than 90% decline in global stocks over the past few decades. Proper genetic resources, such as sequenced, assembled and annotated genomes, are essential to help plan sustainable recoveries by identifying physiological, biochemical and genetic mechanisms that caused the declines or that may lead to recoveries. Here, we present the first sequenced genome of the American eel. This genome contained 305 043 contigs (N50 = 7397) and 79 209 scaffolds (N50 = 86 641) for a total size of 1.41 Gb, which is in the middle of the range of previous estimations for this species. In addition, protein-coding regions, including introns and flanking regions, are very well represented in the genome, as 95.2% of the 458 core eukaryotic genes and 98.8% of the 248 ultra-conserved subset were represented in the assembly and a total of 26 564 genes were annotated for future functional genomics studies. We performed a candidate gene analysis to compare three genes among all three freshwater eel species and, congruent with the phylogenetic relationships, Japanese eel (A. japanica) exhibited the most divergence. Overall, the sequenced genome presented in this study is a crucial addition to the presently available genetic tools to help guide future conservation efforts of freshwater eels.

Regional variation of gene regulation associated with storage lipid metabolism in American glass eels (Anguilla rostrata).

Variation in gene regulation may be involved in the differences observed for life history traits within species. American eel (Anguilla rostrata) is well known to harbor distinct ecotypes within a single panmictic population. We examined the expression of genes involved in the regulation of appetite as well as lipid and glycogen among glass eels migrating to different locations on the Canadian east coast and captured at two different periods of upstream migration. Gene expression levels of three reference and five candidate genes were analyzed by real-time PCR with Taqman probes in recently captured wild glass eels. All gene transcripts were detected in glass eels. Of the five candidate genes, bile salt activated and triacylglycerol lipases were respectively 7.65 and 3.25 times more expressed in glass eels from the St. Lawrence estuary than in those from Nova Scotia, and there was no effect related to the two-week difference in capture date. These two genes explained 82.41% of the dissimilarity between the two rivers. In contrast, glycogen phosphorylase, ghrelin, and leptin receptor genes showed no significant differences in gene transcription. These results confirmed at the molecular level an observation that was recently made at the phenotypic level that glass eels from the St. Lawrence estuary have a greater capacity to use lipid reserves to sustain their metabolic needs. These observations add to the body of evidence supporting the hypothesis that regional phenotypic variation observed in American eel is determined early in life and that part of this variation is likely under genetic control.

RAD sequencing highlights polygenic discrimination of habitat ecotypes in the panmictic American eel.

The two primary ways that species respond to heterogeneous environments is through local adaptation and phenotypic plasticity. The American eel (Anguilla rostrata) presents a paradox; despite inhabiting drastically different environments [1], the species is panmictic [2, 3]. Spawning takes place only in the southern Sargasso Sea in the Atlantic Ocean [1]. Then, the planktonic larvae (leptocephali) disperse to rearing locations from Cuba to Greenland, and juveniles colonize either freshwater or brackish/saltwater habitats, where they spend 3-25 years before returning to the Sargasso Sea to spawn as a panmictic species. Depending on rearing habitat, individuals exhibit drastically different ecotypes [4-6]. In particular, individuals rearing in freshwater tend to grow slowly and mature older and are more likely to be female in comparison to individuals that rear in brackish/saltwater [4, 6]. The hypothesis that phenotypic plasticity alone can account for all of the differences was not supported by three independent controlled experiments [7-10]. Here, we present a genome-wide association study that demonstrates a polygenic basis that discriminates these habitat-specific ecotypes belonging to the same panmictic population. We found that 331 co-varying loci out of 42,424 initially considered were associated with the divergent ecotypes, allowing a reclassification of 89.6%. These 331 SNPs are associated with 101 genes that represent vascular and morphological development, calcium ion regulation, growth and transcription factors, and olfactory receptors. Our results are consistent with divergent natural selection of phenotypes and/or genotype-dependent habitat choice by individuals that results in these genetic differences between habitats, occurring every generation anew in this panmictic species.

Characterization of MHC class IIB for four endangered Australian freshwater fishes obtained from ecologically divergent populations.

Genetic diversity is an essential aspect of species viability, and assessments of neutral genetic diversity are regularly implemented in captive breeding and conservation programs. Despite their importance, information from adaptive markers is rarely included in such programs. A promising marker of significance in fitness and adaptive potential is the major histocompatibility complex (MHC), a key component of the adaptive immune system. Populations of Australian freshwater fishes are generally declining in numbers due to human impacts and the introduction of exotic species, a scenario of particular concern for members of the family Percichthyidae, several of which are listed as nationally vulnerable or endangered, and hence subject to management plans, captive breeding, and restoration plans. We used a next-generation sequencing approach to characterize the MHC IIB locus and provide a conservative description of its levels of diversity in four endangered percichthyids: Gadopsis marmoratus, Macquaria australasica, Nannoperca australis, and Nannoperca obscura. Evidence is presented for a duplicated MHC IIB locus, positively selected sites and recombination of MHC alleles. Relatively moderate levels of diversity were detected in the four species, as well as in different ecotypes within each species. Phylogenetic analyses revealed genus specific clustering of alleles and no allele sharing among species. There were also no shared alleles observed between two ecotypes within G. marmoratus and within M. australasica, which might be indicative of ecologically-driven divergence and/or long divergence times. This represents the first characterization and assessment of MHC diversity for Percichthyidae, and also for Australian freshwater fishes in general, providing key genetic resources for a vertebrate group of increasing conservation concern.

High-throughput SNPs for all: genotyping-in-thousands.

Understanding the genetic structure of species is essential for conservation. It is only with this information that managers, academics, user groups and land-use planners can understand the spatial scale of migration and local adaptation, source-sink dynamics and effective population size. Such information is essential for a multitude of applications including delineating management units, balancing management priorities, discovering cryptic species and implementing captive breeding programmes. Species can range from locally adapted by hundreds of metres (Pavey et al. ) to complete species panmixia (Côté et al. ). Even more remarkable is that this essential information can be obtained without fully sequenced or annotated genomes, but from mere (putatively) nonfunctional variants. First with allozymes, then microsatellites and now SNPs, this neutral genetic variation carries a wealth of information about migration and drift. For many of us, it may be somewhat difficult to remember our understanding of species conservation before the widespread usage of these useful tools. However most species on earth have yet to give us that 'peek under the curtain'. With the current diversity on earth estimated to be nearly 9 million species (Mora et al. ), we have a long way to go for a comprehensive meta-phylogeographic understanding. A method presented in this issue by Campbell and colleagues (Campbell et al. ) is a tool that will accelerate the pace in this area. Genotyping-in-thousands (GT-seq) leverages recent advancements in sequencing technology to save many hours and dollars over previous methods to generate this important neutral genetic information.

Neutral and selective processes shape MHC gene diversity and expression in stocked brook charr populations (Salvelinus fontinalis).

The capacity of an individual to battle infection is an important fitness determinant in wild vertebrate populations. The major histocompatibility complex (MHC) genes are crucial for a host's adaptive immune system to detect pathogens. However, anthropogenic activities may disrupt natural cycles of co-evolution between hosts and pathogens. In this study, we investigated the dynamic sequence and expression variation of host parasite interactions in brook charr (Salvelinus fontinalis) in a context of past human disturbance via population supplementation from domestic individuals. To do so, we developed a new method to examine selection shaping MHC diversity within and between populations and found a complex interplay between neutral and selective processes that varied between lakes that were investigated. We provided evidence for a lower introgression rate of domestic alleles and found that parasite infection increased with domestic genomic background of individuals. We also documented an association between individual MHC alleles and parasite taxa. Finally, longer cis-regulatory minisatellites were positively correlated with MHC II down-regulation and domestic admixture, suggesting that inadvertent selection during domestication resulted in a lower immune response capacity, through a trade-off between growth and immunity, which explained the negative selection of domestic alleles at least under certain circumstances.

The genetic architecture of reproductive isolation during speciation-with-gene-flow in lake whitefish species pairs assessed by RAD sequencing.

During speciation-with-gene-flow, effective migration varies across the genome as a function of several factors, including proximity of selected loci, recombination rate, strength of selection, and number of selected loci. Genome scans may provide better empirical understanding of the genome-wide patterns of genetic differentiation, especially if the variance due to the previously mentioned factors is partitioned. In North American lake whitefish (Coregonus clupeaformis), glacial lineages that diverged in allopatry about 60,000 years ago and came into contact 12,000 years ago have independently evolved in several lakes into two sympatric species pairs (a normal benthic and a dwarf limnetic). Variable degrees of reproductive isolation between species pairs across lakes offer a continuum of genetic and phenotypic divergence associated with adaptation to distinct ecological niches. To disentangle the complex array of genetically based barriers that locally reduce the effective migration rate between whitefish species pairs, we compared genome-wide patterns of divergence across five lakes distributed along this divergence continuum. Using restriction site associated DNA (RAD) sequencing, we combined genetic mapping and population genetics approaches to identify genomic regions resistant to introgression and derive empirical measures of the barrier strength as a function of recombination distance. We found that the size of the genomic islands of differentiation was influenced by the joint effects of linkage disequilibrium maintained by selection on many loci, the strength of ecological niche divergence, as well as demographic characteristics unique to each lake. Partial parallelism in divergent genomic regions likely reflected the combined effects of polygenic adaptation from standing variation and independent changes in the genetic architecture of postzygotic isolation. This study illustrates how integrating genetic mapping and population genomics of multiple sympatric species pairs provide a window on the speciation-with-gene-flow mechanism.