Adaptation to Low Salinity Promotes Genomic Divergence in Atlantic Cod (Gadus morhua L.).

How genomic selection enables species to adapt to divergent environments is a fundamental question in ecology and evolution. We investigated the genomic signatures of local adaptation in Atlantic cod (Gadus morhua L.) along a natural salinity gradient, ranging from 35‰ in the North Sea to 7‰ within the Baltic Sea. By utilizing a 12 K SNPchip, we simultaneously assessed neutral and adaptive genetic divergence across the Atlantic cod genome. Combining outlier analyses with a landscape genomic approach, we identified a set of directionally selected loci that are strongly correlated with habitat differences in salinity, oxygen, and temperature. Our results show that discrete regions within the Atlantic cod genome are subject to directional selection and associated with adaptation to the local environmental conditions in the Baltic- and the North Sea, indicating divergence hitchhiking and the presence of genomic islands of divergence. We report a suite of outlier single nucleotide polymorphisms within or closely located to genes associated with osmoregulation, as well as genes known to play important roles in the hydration and development of oocytes. These genes are likely to have key functions within a general osmoregulatory framework and are important for the survival of eggs and larvae, contributing to the buildup of reproductive isolation between the low-salinity adapted Baltic cod and the adjacent cod populations. Hence, our data suggest that adaptive responses to the environmental conditions in the Baltic Sea may contribute to a strong and effective reproductive barrier, and that Baltic cod can be viewed as an example of ongoing speciation.

Introduction

Teleost fishes occupy a range of aquatic habitats from freshwater environments to extreme marine environments with salinities well above 40‰. In order to spawn in marine habitats, teleosts have eggs adapted to a hyper-osmotic environment, resulting in less permeable and highly hydrated eggs that compensate for water outflow until osmoregulatory organs develop (Fyhn et al. 1999; Finn and Kristoffersen 2007). After the egg stage, osmoregulation primarily takes place in the intestines (Whittamore 2011; Laverty and Skadhauge 2012) and through salt excretion from specialized ionocyte cells in the gill epithelium (Evans 2008; Hiroi and McCormick 2012). Adaptation to specific osmotic conditions thus requires a wide range of molecular and physiological modifications in order to maintain water homeostasis both at the egg stage and in adult fish. Several studies have demonstrated population-specific acclimatization to different salinity environments through differential gene expression (see, e.g., Boutet et al. 2006; Larsen et al. 2011; Papakostas et al. 2014). Although many genes involved in the physiological response to variation in salinity have been identified, the genetic basis of broad salinity tolerance remains unclear (Tine et al. 2014). In addition to the characteristic physical environment (salinity, temperature, and oxygen), Baltic cod also experiences a different parasite load compared with Kattegat and North Sea cod populations (Mellergaard and Lang 1999; Perdiguero-Alonso et al. 2008). However, the role of parasites in promoting genetic divergence has mixed support (Karvonen and Seehausen 2012).

The Baltic Sea is one of the world’s largest semienclosed brackish seas, originating approximately 8,000 years ago (Zillén et al. 2008). Since then, the Baltic Sea has been colonized by both freshwater and marine teleosts (Ojaveer and Kalejs 2005). Limited water exchange and almost no tidal flow allow for stable salinity conditions, ranging from nearly fresh water in the northern Baltic Sea to around 30‰ at the border to the North Sea (fig. 1). Fossil records suggest that all contemporary marine teleost fish species have descended from a freshwater ancestry (Romer and Grove 1935; Griffith 1994; Long 1995; Fyhn et al. 1999), reflected by their hypo-osmotic state at the adult stage and their eggs and sperm in seawater (Griffith 1987; Fyhn et al. 1999; Finn and Kristoffersen 2007). The adaptation of neutral egg buoyancy toward those salinity levels found in the marine environment could arguably impede the successful colonization of less saline environments, as hyper-osmotic eggs would sink to the bottom, requiring an even higher degree of oocyte hydration, which is the case for Atlantic cod adapted to the Baltic environment (Nissling et al. 1994). Despite such considerations, some marine fishes with pelagic eggs, such as sprat, plaice and Atlantic cod, have colonized and successfully reproduce in low saline waters such as the Baltic Sea. Indeed, Atlantic cod colonized the central Baltic Sea some 8–6,000 YBP (Schmölcke 2006) when the surface salinity was 12–14‰. Since then, surface salinity has gradually decreased to around 7‰ today and around 14‰ at the spawning depth for Atlantic cod (50–100 m; Ignatius et al. 1981). Hence, the selection pressure to adapt to the low saline waters is a major force influencing the spawning success of cod in the Baltic Sea (Westin and Nissling 1991; Nissling et al. 1994), limited by factors such as egg buoyancy, sperm motility, and general osmoregulation. Importantly, Baltic cod spawns in a different season, compared with the North Sea and Western Baltic cod stocks, likely promoting reproductive isolation (Brander 1993; Wieland 2000; Brander 2005). The adaptation of Atlantic cod to low saline conditions is thus likely of relatively recent evolutionary origin, even though we cannot exclude the possibility that the founders of the Baltic population were already adapted to low salinity. This provides an excellent opportunity to study the genomic architecture behind salinity adaptation in a natural environment. As the ecological adaptation of Atlantic cod to a low-saline environment may contribute to reduced gene flow and thereby promote population divergence (cf. Nosil 2012), investigation of the genomic architecture of Baltic cod may give insights into ecological speciation in nature, and especially the genetic link between adaptation and reproductive isolation (cf. Seehausen et al. 2014).

F

Sampling locations of Atlantic cod specimens. Samples were obtained between 2002 and 2008 using trawl. Only mature specimens were selected for genetic analysis. See table 1 for a detailed description of the samples. Average surface salinity is denoted in ‰.

Table 1

Location and Details of the Atlantic Cod Samples Included in This Study, Combined with Basic Population Genetic Parameters

Population ID	Location	Sampling Time	Lat.	Long.	Spawning	Sample Size	Avg. Call Rate	# Polymorphic Loci	Ho (SD)	He (SD)
North Sea	North Sea	Mar 2002	N55.60	E05.85	Yes	42	0.983	8,659	0.362 (0.149)	0.361 (0.138)
Kattegat	Kattegat	Feb 2004	N56.90	E12.15	Yes	48	0.995	8,690	0.360 (0.145)	0.360 (0.135)
Öresund	Öresund	Feb 2003	N55.95	E12.70	Yes	48	0.995	8,715	0.363 (0.146)	0.363 (0.134)
Baltic	Baltic Sea					56	0.993	8,517	0.344 (0.155)	0.345 (0.147)
Baltic	Bornholm 05	June 2005	N55.50	E16.00	Yes	8	0.989	7,605	0.381 (0.188)	0.382 (0.135)
Baltic	Bornholm 04	June 2004	N55.59	E16.30	Yes	8	0.993	7,389	0.379 (0.192)	0.378 (0.137)
Baltic	Gotland	Nov 2008	N57.15	E18.78	No	20	0.994	8,267	0.359 (0.168)	0.356 (0.144)
Baltic	Öland	Nov 2008	N57.40	E17.00	No	20	0.994	8,301	0.353 (0.163)	0.357 (0.143)

Note.—All sampled individuals were adults. Estimates of observed (Ho) and expected heterozygosity (He) were calculated using ARLEQUIN 3.5.1.3. SD, standard deviation. Latitude and longitude values are given in degrees and minutes. The two samples from Bornholm are denoted as Bornholm 04 and 05 based on the sampling year.

Atlantic cod has high fecundity, pelagic eggs and larvae and large effective population sizes (Ward et al. 1994; Poulsen et al. 2006; Therkildsen et al. 2010), suggesting that selection, rather than genetic drift, is the main driver for genetic differentiation (Allendorf et al. 2010). Dispersal and gene flow could however prevent genetic structuring caused by local adaptation. Nevertheless, low levels of genetic structuring may not necessarily extend to genomic regions of functional importance (Weir et al. 2005; Nosil et al. 2009; Orsini et al. 2013) and highly differentiated genomic regions have been found in marine fish populations (Cano et al. 2008; Nielsen, Hemmer-Hansen, Larsen, et al. 2009; Nielsen, Hemmer-Hansen, Poulsen, et al. 2009; Hohenlohe et al. 2010; André et al. 2011; Shimada et al. 2011; Jones et al. 2012; Hemmer-Hansen et al. 2013). It has been shown in several cases that adaptive genetic variation may be selectively maintained, despite homogenizing effects of gene flow (see, e.g., Cano et al. 2008; Hauser and Carvalho 2008; Nielsen, Hemmer-Hansen, Larsen, et al. 2009; Bradbury et al. 2013; Defaveri, Jonsson, et al. 2013; Teacher et al. 2013; reviewed in Orsini et al. 2013). To identify genomic regions under selection and to get a better understanding of the evolutionary processes that generate such regions are important as the sheer number of genes and the number of potential traits evolving together are related to the size of the “islands of divergence” (see, e.g., Via 2009; Feder and Nosil 2010; Feder, Gejji, et al. 2012; Nosil and Feder 2012; Via et al. 2012; Flaxman et al. 2013).

So far, the genomic studies targeting the Baltic cod have used a relatively restricted number of genetic markers (e.g., Nielsen et al. 2003, 2012; Poulsen et al. 2006; Nielsen, Hemmer-Hansen, Poulsen, et al. 2009; Bradbury et al. 2010, 2013; Larsen et al. 2011; Hemmer-Hansen et al. 2013), limiting inference of the genomic architecture underlying salinity adaptation. Here, we investigate divergence patterns using 8,809 polymorphic single nucleotide polymorphisms (SNPs) distributed throughout the entire genome of Atlantic cod (based on an Illumina 12 K SNPchip described in The Cod SNP Consortium, in preparation), providing the first individually genotyped whole genome-wide approach to date on this species. We used different statistical outlier approaches on four different populations across the steep salinity gradient from the fully marine North Sea to the low saline Baltic Sea, identifying SNPs and genomic regions under selection. These genomic regions were also identified by a landscape genomic analysis, corroborating that certain genomic regions with underlying SNPs are influenced by defined environmental variables.

Materials and Methods

Sample Collection

Fin clips or spleen from adult Atlantic cod specimens (n = 194) were collected from seven localities and stored on ethanol or RNAlater (Qiagen). Sample sizes per location ranged from 8 to 48 individuals (fig. 1 and table 1). All individuals were collected during spawning, except for the Öland and the Gotland samples that presumably belong to the Bornholm basin spawning population. The samples collected within the Baltic Sea were treated as a one single population (see below).

DNA Extraction and Genotyping

DNA was extracted using the E.Z.N.A Tissue DNA kit (Omega Bio-Tek) and normalized to 100 ng/μl using a NanoDrop DN1000 instrument (Thermo Scientific Inc.) prior to genotyping, using a 12 K Illumina SNP-chip (The Cod SNP Consortium, in preparation). Out of the 10,913 SNPs on the SNPchip, 8,809 SNPs were polymorphic, showed good clustering/separation, had a call rate of greater than 95% and also showed Mendelian segregation among family individuals (data not shown). Each individual SNP locus was manually inspected and clusters were adjusted if necessary, using the Genome Studio 2011.1 software (Illumina). Of these, 262 SNPs are previously published (Moen et al. 2008; Hubert et al. 2010), 648 SNPs target 313 selected candidate genes (initially selected for a family-based study on growth and maturity, hence these markers could be treated as random genes in this study), 1,554 SNPs are nonsynonymous whereas the remaining 6,345 SNPs are randomly distributed throughout the 23 different LGs in the Atlantic cod genome (the source of each SNP is listed in supplementary table S4, Supplementary Material online). All SNPs are referred to by their rs# or ss#, available in dbSNP (www.ncbi.nlm.nih.gov/SNP, last accessed May 1, 2015). The LGs are numbered as described by Hubert et al. (2010) and the order of SNPs is based on preliminary linkage data (Lien S, unpublished data).

Population Genetics

Within each population, estimates of observed (Ho) and expected heterozygosity (He) were calculated in ARLEQUIN 3.5.1.3 (Excoffier and Lischer 2010), whereas allele frequencies were calculated for all SNPs in all populations using Convert 1.31 (Glaubitz 2004). Departure from HWE was tested locus by locus in each population using the exact test by Guo and Thompson (1992) implemented in ARLEQUIN with 100,000 iterations and a Markov Chain of 1 million. We corrected for multiple testing by computing the q value for each locus, using the QVALUE package (Storey 2002) in R (R Core Team 2012), using a q value of 0.05 as a threshold for significance.

Based on results from the outlier analyses (see below), each SNP was categorized either as an outlier or as a neutral marker. Consequently, three data sets were constructed and used in the FST and STRUCTURE analyses. The first data set contains all 8,809 SNP loci and is referred to as the “full data set.” The second data set contains 6,913 loci that were not under selection in any of the outlier tests (prior to corrections for multiple testing), and is referred to as the “neutral data set.” The third data set, the “outlier data set,” contains 233 physically unlinked outlier SNPs from all 23 LGs. Hence, outlier SNPs in high LD (r > 0.8) with each other were represented by a single tag SNP in this data set, selected using PLINK v1.07 (Purcell et al. 2007). To minimize the use of false positives in this data set, only loci detected to be under directional selection by two or more outlier tests, after corrections for multiple testing, were included.

Locus specific FST values for all pairwise population comparisons and the weighted average FST values (Weir and Cockerham 1984) between all population pairs for all data sets were calculated in ARLEQUIN, using 10,000 permutations for significant testing. We used the Bayesian clustering model implemented in the program STRUCTURE v2.3.4 (Pritchard et al. 2000) to identify major genetic clusters (see supplementary text S3, Supplementary Material online, for details).

Linkage Disequilibrium

We evaluated the presence of LD among all 8,809 SNPs, using a data set consisting of all populations combined, calculating both inter- and intrachromosomal LD between all SNPs. LD was quantified in PLINK with the r 2 estimate, reporting all r values for the intrachromosomal calculations and r > 0.2 for the interchromosomal calculations.

Outlier Detection

Roesti suggested that uninformative markers at low frequency be excluded from outlier analyses to increase the power of genome scans. The initial selection of SNPs on the Illumina SNP chip, based on a sequencing of eight individuals, resulted in a relatively high minor allele frequency of most of the SNPs, hence no further SNPs were excluded prior to the outlier analyses. There are several challenges in outlier studies, including detection of false positives, false negatives, and complications due to underlying population structuring. For this reason, three independent methods were used to identify candidate loci under selection.

First, we used a Bayesian regression approach implemented in BAYESCAN v2.1 (Foll and Gaggiotti 2008) which, based on FST coefficients, measures the discordance between global and population-specific allele frequencies. The degree of differentiation, based on FST, is decomposed into a locus-specific component (α), shared by all populations, and a population-specific component (β), shared by all loci. Selection is assumed when alpha is necessary to explain the observed pattern of diversity. To control for variation in the BF distribution caused by randomness and imperfection in each independent run of BAYESCAN, the median value of ten independent runs was calculated for each SNP. We carried out the simulations on the global data set, including all four populations, using stringent criteria, assuming selection to be 10%. The FDR was set to 0.01. We also calculated empirical P values, as described in Lotterhos and Whitlock (2014), based on 3,137 neutral SNPs that are more than 5 kb away from any annotated gene. We report both the median log10 values of the posterior odds (PO) as well as the q values and the empirical P values, only accepting log10(PO) above 2.0 and 1.0 as “decisive” and “strong” evidence for selection (according to Jeffreys 1961) and empirical P values < 0.01 as significant.

Second, we used the FDIST2 approach by Beaumont and Nichols (1996) implemented in the software LOSITAN (Antao et al. 2008). When using this approach, comparisons are made of FST values in relation to heterozygosity of individual loci, based on a neutral distribution, generated by means of coalescence simulations in a symmetric island migration model at mutation–drift equilibrium. We carried out the simulations on all four populations together under the Infinite Allele Method with 1 million simulations, a confidence interval of 0.99 with an FDR of 0.01, using the neutral mean FST option and forcing mean FST option. We report that the median value of ten independent runs of LOSITAN. FDIST2 is known to be robust to a wide range of nonequilibrium conditions, but it can be sensitive to demographic variations among populations as well as hierarchical genetic structure that may result in the detection of false outliers (Storz 2005). We corrected for multiple testing by computing the q value for each locus, using the QVALUE package (Storey 2002) in R (R Core Team 2012) and calculated empirical P values, as described above. We report both the q values and the empirical P values, using 0.01 as a threshold for significance.

Third, outlier tests were performed pairwise between all four population pairs, as overall divergence based on global FST values may not detect candidates that are under selection in only some of the populations (Vitalis et al. 2001) and also to identify which populations that are under directional selection for each locus. By performing pairwise comparisons of populations, we omit some of the methodological weakness associated with population structure/demographic processes in the data sets. BAYESCAN was used as this is considered to be the most conservative test (Narum and Hess 2011). We used the same stringent criteria as in the global analyses, reporting both the median log10(PO) of ten independent BAYESCAN runs as well as the q values and the empirical P values, only accepting log10(PO) above 2.0 and 1.0 as decisive and strong evidence for selection (according to Jeffreys 1961) and empirical P values < 0.01 as significant.

Environmental Association

Correlations between allele frequencies and the environmental variables, salinity (‰), temperature (°C) and oxygen concentration, all at surface and at spawning depth, were tested using the Bayesian method by Coop et al. (2010) implemented in the software BAYENV 2.0 (Günther and Coop 2013). The environmental variables were retrieved from ICES and Helcom, as detailed in table 3. For the combined Baltic population, data from the proposed spawning area of the Bornholm Basin were used in the analyses. The rationale behind the method by Coop is to estimate a neutral covariance matrix based on a large set of control loci. Second, a test for covariance between the environmental variables and the population specific allele frequencies at each SNP is performed, using the neutral covariance matrix as a reference (null model) to control for shared population history and gene flow. For each test, a BF is calculated based on the ratio of the posterior probabilities between the two models. A high BF indicates a support for the alternative model where the environmental variable has a linear effect on the locus of interest. As a basis for the null model, we used 3,137 neutral SNP loci that are more than 5 kb away from any annotated gene as control loci. BAYENV 2.0 was run for 500,000 Markov chain Monte Carlo (MCMC) iterations and the covariance matrix was created using an average of all matrices estimated by the program (output every 500 iteration). Next, all SNPs were tested for correlation with the environmental variables using 500,000 MCMC iterations each. Following Blair , to control for variation in the BF distribution caused by randomness in the MCMC algorithm and imperfection in the null model, 32 independent runs were carried out. The starting point for the algorithm (random seed) was drawn randomly for each run. Based on the BAYENV 2.0 results, q values were calculated according to de Villemereuil et al. (2014). Empirical p values were calculated as described above. SNPs showing median log10(BF) > 2 indicating “decisive evidence” and median log10(BF) > 1 indicating “strong evidence” according to Jeffreys (1961) were regarded as significant results when also significant at q < 0.01.

Table 3

Environmental Conditions at the Atlantic Cod Sampling Locations

Population ID	Location	Spawning Month	Spawning Depth (m)	Salinity Surface	Salinity Spawning Depth	Temp. Surface	Temp. Spawning Depth	Oxygen Surface	Oxygen Spawning Depth	Source
North Sea	North Sea	Feb.–April	40–60	34.7	34.8	6.0	6.1	6.8	6.7	ICES
Kattegat	Kattegat	Feb.–April	60–100	24.2	34.1	3.5	6.3	8.0	6.5	Helcom
Öresund	Öresund	Feb.–April	40–53	16.4	30.9	3.4	6.4	8.5	6.2	Helcom
Baltic	Bornholm 05	June–Aug.	50–82	7.4	13.6	13.6	6.2	7.2	4.2	Helcom
Baltic	Bornholm 04	June–Aug.	50–90	7.4	13.6	13.6	5.7	7.4	4.5	Helcom
Baltic	Gotland	June–Aug.	50–90	6.9	8.6	9.3	4.4	7.3	4.1	Helcom
Baltic	Öland	June–Aug.	50–82	6.7	8.4	8.4	4.8	7.6	3.4	Helcom

Note.—Salinity, temperature, and oxygen values were obtained online from ICES or Helcom. Salinity is denoted in ‰, temperature in °C and oxygen in concentration. The two samples from Bornholm are denoted as Bornholm 04 and 05 based on the sampling year.

SNP Annotation

The original SNP chip design was based on the ATLCOD1A genome assembly (www.codgenome.no, last accessed May 1, 2015) that has a different coordinate system than the genome assembly for which Ensembl annotation is available (ATLCOD1C). Therefore, the location of each SNP was determined by aligning 200 bp of flanking sequence to the ATLCOD1C assembly (Star et al. 2011) using BWA aln -o 1 -n 0.04, seeding disabled (Li and Durbin 2009). Only unique hits were selected. The locations of annotated genes were downloaded from the Ensembl database Gadus_morhua.gadMor1.73 using Biomart (Kasprzyk 2011). The distance of the flanking sequence of the SNPs from genes (exons) was determined using BEDclosest including the options -t “first” and -d (BEDtools v2.16.2; Quinlan and Hall 2010). The protein transcripts of Ensembl genes that were associated with the SNPs through this approach were annotated with BLAST2GO (Conesa et al. 2005) using public database b2g_sep13. Protein transcripts were aligned to the refseq_protein data using the BLASTP algorithm allowing a maximum of 20 hits with a minimum e value of 1 E-3, as implemented in BLAST2GO. Apart from setting the evidence code weight of IEA (electronic annotation evidence) to 1, default weights were used. Annotation was augmented using the Annex function in BLAST2GO.

Results

SNP-Set and Samples

Based on a 12 K SNPchip, we analyzed a set of 8,809 polymorphic SNPs (see Materials and Methods for details), distributed over all 23 linkage groups (LGs) with an average distance of 94,000 bp between SNPs, based on a genome size of 830 Mb (Star et al. 2011). A total of 5,565 SNPs were located within 5,000 bp of 4,522 Ensembl annotated genes (supplementary table S4, Supplementary Material online). A final number of 194 individuals of Atlantic cod (individual call rate > 95%) from seven different localities were investigated (fig. 1, table 1). As the four Baltic localities likely comprised individuals from the same breeding population, these data were pooled into a single sample, which was justified by assignment testing (see Materials and Methods and supplementary text S1, Supplementary Material online) and further supported by tests for deviation from Hardy–Weinberg equilibrium (HWE) (supplementary text S2, Supplementary Material online). Only one SNP locus was significantly out of HWE after false discovery rate (FDR) correction (q < 0.05) in the pooled Baltic sample, indicating no sign of a Wahlund effect. The number of polymorphic loci, observed- and expected heterozygosity (Ho and He) were similar in all populations (table 1).

Population Genetic Structuring

In addition to the full data set, SNPs were categorized either as outliers or as neutral, based on outlier analyses (see Outlier Detection and Environmental Association section). In the full- and neutral data set, all average pairwise FST comparisons were statistically significant except between Kattegat and Öresund, indicating little genetic differentiation between these two populations (table 2). In comparison, based on 233 unlinked SNPs from the outlier data set, all pairwise FST comparisons were significant, indicating some genetic differentiation also between the Kattegat and Öresund populations.

Table 2

Pairwise FST Values among Atlantic Cod Populations, Using a Full-, Neutral-, and Outlier Data Set (See Text for Details)

	North Sea	Kattegat	Öresund	Baltic
Full data set (8,809 SNPs)
North Sea	—	0.01595	0.01467	0.07095
Kattegat		—	0.00026	0.05140
Öresund			—	0.04824
Baltic				—
Neutral data set (6,913 SNPs)
North Sea	—	0.00962	0.00996	0.04840
Kattegat		—	0.00007	0.03138
Öresund			—	0.03061
Baltic				—
Outlier data set (233 SNPs)
North Sea	—	0.12065	0.10278	0.40130
Kattegat		—	0.00825	0.32685
Öresund			—	0.29778
Baltic				—

Note.—Values in bold are significant values (P values < 0.05), calculated using ARLEQUIN 3.5.1.3.

Locus-specific FST values for all pairwise population comparisons indicate a genome-wide pattern of high FST when comparisons are made to the Baltic Sea, whereas few discrete regions of the genome show elevated FST values in the remaining comparisons (fig. 2). A total of 40 private alleles (supplementary table S2, Supplementary Material online) were detected among the populations, 38 of which were rare in frequency (<0.04). The highest number of private alleles (23) was observed in the North Sea population and the lowest number (3) was observed in the Baltic population. The two private alleles of higher frequency (0.179 and 0.232) were found in the North Sea population, in complete linkage disequilibrium (LD) with each other (ss1712302407 and ss1712304757, located 241 bp apart, in LG3). In total, 334 fixed allele frequencies were detected in any population; however, only 84 of these showed frequency differences of more than 0.1 in any pairwise comparison. Out of these, 59 were unique to the Baltic population (19 potentially under selection) whereas 11, 2 and 1 were unique to the North Sea, Kattegat and Öresund populations, respectively. In addition, 59 outlier SNPs were close to fixation (frequency > 0.95; supplementary table S2, Supplementary Material online).

F

Pairwise Fst values between Atlantic cod populations using 8,809 SNPs in 23 LGs. SNPs are ordered according to LG and position within LGs based on the results of a preliminary SNP linkage map (Lien S, unpublished data). LG nomenclature follows that of Hubert et al. (2010). Median Fst estimates between the population pairs are denoted in parentheses.

Bayesian cluster analyses (STRUCTURE analysis; supplementary text S3, Supplementary Material online) strongly support a separation between the Baltic population and the other populations for all three data sets (supplementary fig. S1, Supplementary Material online). At higher K values, the Kattegat and the Öresund populations, representing the transition zone between the North Sea and the Baltic Sea, form an intermediate group with more similarity to the North Sea than to the Baltic population (supplementary fig. S2, Supplementary Material online). The highest delta K value was detected for K = 2 for all three data sets, indicating a deep-rooted structuring in the data and that the Baltic cod is separated by a larger genomic distance to all the other populations. The identified clusters correspond well with the overall population structure estimated by pairwise FST comparisons (table 2).

LD was evaluated among all 8,809 SNPs, independent of LG localization to detect both inter- and intrachromosomal LD. As expected, we find SNPs with high r values (>0.75) within all LGs (supplementary table S3, Supplementary Material online). The distributions of the intrachromosomal r values among the different LGs show that SNPs with high r values often reside within the same scaffolds and hence are physically close. However, for LG2, 7 and 12, an extensive number of SNPs with high LD (fig. 3) covers multiple scaffolds including 15 and 5 scaffolds in two close but separated blocks in LG2, 14 scaffolds in LG7, and 31 scaffolds in LG12 (supplementary table S3, Supplementary Material online). The LD pattern (fig. 3) also shows less distinct blocks of elevated r values in LG4, 10, and 17. Interestingly, relatively high levels of r values (>0.3) also occur between SNPs on different LGs, particularly between LG1 and 2 and between LG2 and 4 (supplementary table S3, Supplementary Material online).

F

Global outlier pattern and LD in Atlantic cod. The extent of pairwise LD among loci within 23 LGs, measured as r2, is estimated using all individuals and populations. A color bar on the right hand side denotes the strength of the r2 values. A global outlier pattern (underneath each LD plot) based on the same SNP data is visualized using a moving average of the median log10(BF), calculated using the global BAYESCAN outlier analyses (see text for details). The SNPs are plotted according to LG and their respective position within the LGs, based on the results of a preliminary SNP linkage map (Lien S, unpublished data). LG nomenclature follows Hubert et al. (2010).

Outlier Detection and Environmental Association

We used two outlier detection approaches, BAYESCAN and LOSITAN, combined with a method that explores the correlation between environmental variables and allele frequencies (BAYENV 2.0) to uncover various aspects of the genomic regions under selection (see Materials and Methods for details).

In the global outlier analyses comprising all four populations, we identified 123 (1.4%) and 114 (1.3%) SNPs as candidates for divergent selection (log10(Bayes factor [BF]) > 1, empirical P < 0.01), using BAYESCAN and LOSITAN, respectively. Pairwise outlier analyses using BAYESCAN, revealed 228 (2.6%) SNPs as candidates for divergent selection (supplementary fig. S3 and table S4, Supplementary Material online). However, the pattern and the number of outliers in each pairwise comparison varied (fig. 4b–g). In total, 266 SNPs (3.0%) were detected as candidates for divergent selection using global or pairwise outlier test. LG2 and 12 had the highest proportion of SNPs potentially under selection (36 and 98, respectively; supplementary table S4, Supplementary Material online). These two LGs contain large LD blocks, where most, but not all, of these outliers reside (fig. 3).

F

Global and pairwise outlier patterns in Atlantic cod. Manhattan plots of global and pairwise outlier analyses based on median log10(BF) from ten replicate runs of BAYESCAN. The SNPs are plotted according to LG and their respective position within the LGs along the X axis based on the results of a preliminary SNP linkage map (Lien S, unpublished data). LG nomenclature follows Hubert et al. (2010). The dotted line at log10(BF) = 1 indicates “strong association” and the solid line at 2 indicates “decisive association” according to Jeffreys (1961).

Correlations between allele frequencies and salinity (‰), temperature (°C) and oxygen concentration at surface and at spawning depth (table 3) were investigated using a landscape genomic approach (BAYENV 2.0). Significant association (median log10(BF) > 1 and q < 0.01) with at least one environmental variable was detected in 234 (2.6%) of the 8,809 examined SNP loci, distributed across all LGs except LG6 (supplementary fig. S3 and table S4, Supplementary Material online). The pattern and the number of SNPs with correlation to each environmental variable varied (fig. 5a–f). An association with salinity and oxygen level at spawning depth is clearly observed at a genome-wide scale, and in particular toward the end of LG2, whereas temperature shows a weak correlation in LG12 (fig. 5). There is a strong correlation between SNPs associated with salinity and oxygen, whereas SNPs associated with temperature are not strongly correlated with salinity or oxygen SNPs (fig. 5).

F

SNPs associated with environmental variables in Atlantic cod. Manhattan plots of SNP association with salinity, temperature, and oxygen level (all at surface and at spawning depth, see table 1) based on median log10(BF) from 32 independent runs of BAYENV. The SNPs are plotted according to LG and their respective position within the LGs along the X axis, based on the results of a preliminary SNP linkage map (Lien S, unpublished data). LG nomenclature follows Hubert et al. (2010). The dotted line at log10(BF) = 1 indicates strong association and the solid line at 2 indicates decisive association according to Jeffreys (1961).

In total, 326 SNPs (3.7%) were candidates for divergent selection or were associated with an environmental variable (supplementary fig. S3, Supplementary Material online). The overlap between the outlier and the landscape genomic data sets reveals 174 SNPs with significant evidence of selection in both the outlier approach (global or pairwise comparisons) and the landscape genomic approach, distributed on 102 different scaffolds on most LGs, except for LG6, and 23 (supplementary table S4, Supplementary Material online). Of these 174 SNPs, 129 loci were located in or within 5 kb of an annotated gene, of which 61 were located in exons and 48 were nonsynonymous substitutions causing amino acid changes (supplementary table S4, Supplementary Material online). Further, 123 of the 174 outlier SNPs reside in five LGs (LG1, 2, 4, 8, and 12; supplementary fig. S3, Supplementary Material online). In LGs where outliers are detected, they do not seem to be clustered within the LGs, except for the outliers in LG2, 12 and to some extent in LG4, 8, 10, 17, and 18 (fig. 3).

Identifying Genes and Genomic Regions under Selection

Three LGs (LG2, 7, and 12) have substantial long-distance LD (fig. 3) and the LD pattern is similar within the different populations (data not shown). In LG2, 51 SNPs were candidates for selection, out of which 32 were located in two tightly linked regions (consisting of 22 and 10 SNPs, respectively) toward the end of the LG (fig. 3). All SNPs that were candidates for selection on this LG were associated with oxygen and salinity (primarily at spawning depth). Further, all of these outliers were detected as pairwise outliers relative to the Baltic Sea, but not in any of the other pairwise comparisons, except for an SNP (ss1712297946) close to a C-type lectin gene (CLEC16A) toward the beginning of the LG. Combined, these results indicate that the outliers on LG2 are of particular relevance for the adaptation to the low saline conditions in the Baltic Sea. In LG7, one single outlier SNP (ss1712305126) with decisive evidence of selection was located within an exon in an Na+/K+ ATPase gene (ATP1A1) and was not associated with the large LD region present on LG7 (fig. 3). In LG12 the long-distance LD area contains 37 decisive outliers, but notably only detected in the pairwise comparison between North Sea and Kattegat/North Sea and Öresund. No decisive pairwise outliers were detected in any of the comparisons with the Baltic population (supplementary table S4, Supplementary Material online). However, the majority of outliers detected on LG12 appears to be under a less stringent selection regime as the number of outliers rise to 103, using log10(BF) > 1 and empirical P < 0.01 as a threshold for significance. Similarly, there is a clear selection pattern within the linked region in LG12, associated with temperature—but not with salinity or oxygen (fig. 5 and supplementary table S4, Supplementary Material online).

In addition, we find three less distinct LD regions (LG4, 10, and 17) that also harbor several outlier candidates. In LG4, the LD region contains 11 of the 15 decisive outlier SNPs. All of these SNPs showed signs of selection in the pairwise comparisons with the Baltic population, but not with any of the other populations and were all decisively associated with salinity and oxygen, primarily at spawning depth (fig. 5 and supplementary table S4, Supplementary Material online). Hence, there are strong indications that also the genomic region on LG4 has relevance for adaptation to the low saline conditions in the Baltic Sea. Notably, SNPs within the LD block on LG4 are linked to SNPs on the LD blocks on LG2 (0.3 < r2 < 0.4; supplementary table S3, Supplementary Material online), possibly indicating a common selection pressure on the two genomic regions. In LG10, the main outlier peak lies just outside the LD area, with some less significant outlier SNPs residing within the LD area (fig. 3). All SNPs under selection on LG10 were associated with salinity and oxygen at spawning depth (supplementary table S4, Supplementary Material online). In LG17, seven SNPs displayed decisive evidence of selection and four of these SNPs, associated with oxygen at surface, were within the LD region (fig. 3).

A total of 51 SNPs, within 5 kb of 43 annotated genes (see Discussion) were decisively associated with salinity and oxygen and were specific to the Baltic Sea, relative to the other three sampling locations (table 4). These pairwise outliers are likely to identify the strongest genomic differentiation, shaped by the environment in the Baltic Sea. Surrounding these strong outliers, other less pronounced outlier SNPs were also detected (figs. 4 and 5 and supplementary table S4, Supplementary Material online).

Table 4

SNPs Closer than 5 kb to an Annotated Gene with Decisive Association (i.e., median log10(BF) > 2 and Empirical P < 0.01, see text for details) to Salinity and Oxygen at Spawning Depth that Are Also Strong Pairwise Outliers across the Salinity Barrier

					BAYENV 2.0 Log10(BF)		BAYESCAN Log10(PO)
Locus	LG	Distance to Exon	HGNC	Gene ID	Salinity Depth	Oxygen Depth	Baltic Öresund	Baltic Kattegat	Baltic North Sea	North Sea Kattegat	North Sea Öresund	Öresund Kattegat
ss1712302658	1	0	TSHB	ENSGMOG00000011010	2.78	2.62	2.75	2.38	0.37	−1.02	−1.03	−1.02
ss1712302659	1	244	TSHB	ENSGMOG00000011010	2.72	2.57	2.64	2.22	0.38	−1.02	−1.02	−1.00
ss1712299135	1	21 (I)	GRIP2A	ENSGMOG00000017983	2.58	2.71	−0.37	1.32	2.18	−1.02	−0.73	−0.99
ss1712303317	1	0 (NS)	FBLN2	ENSGMOG00000017143	2.77	2.68	2.92	2.27	1.30	−1.01	−1.03	−1.00
ss1712303323	1	0 (NS)	COL2A1	ENSGMOG00000009692	4.36	4.34	>4	>4	>4	−1.02	−1.02	−1.00
ss1712300221	1	0 (NS)	—	ENSGMOG00000013560	4.98	4.94	>4	>4	>4	−1.04	−1.03	−1.00
ss1712304198	1	0 (NS)	—	ENSGMOG00000013560	2.80	2.88	2.92	2.55	3.22	−1.01	−1.01	−1.00
ss1712304199	1	0 (NS)	—	ENSGMOG00000013560	2.88	2.94	2.96	2.82	3.40	−1.02	−1.02	−1.03
ss1712304200	1	0 (NS)	—	ENSGMOG00000013560	4.64	4.59	>4	>4	>4	−1.03	−1.03	−1.02
ss1712304720	2	0 (NS)	ZFAND2A	ENSGMOG00000004883	4.11	4.14	>4	3.70	>4	−1.02	−1.02	−1.02
ss1712296776	2	832 (I)	TMEM235	ENSGMOG00000013183	2.71	2.62	−0.88	1.82	−0.31	−1.02	−1.00	−0.84
ss1712298606	2	120 (I)	TMEM104	ENSGMOG00000018126	3.40	3.36	1.86	1.78	1.51	−1.02	−1.02	−1.00
ss1712302705	2	621	SSTR3	ENSGMOG00000013330	2.68	2.63	−0.86	1.48	−0.18	−1.01	−1.00	−0.91
ss1712302707	2	1,058	SSTR3	ENSGMOG00000013330	3.38	3.38	−0.91	1.31	0.27	−1.02	−0.86	−0.88
rs119055764	2	0	SLC25A39	ENSGMOG00000009149	4.31	4.31	0.05	3.00	2.27	−1.02	−0.67	−0.89
ss1712302819	2	337	PRL	ENSGMOG00000017998	4.55	4.39	3.70	>4	2.82	−1.03	−1.01	−1.02
ss1712299603	2	202	PDGFAB	ENSGMOG00000013650	3.89	3.80	0.21	2.92	1.30	−1.01	−0.98	−0.92
ss1712298457	2	4,924	PACSIN1A	ENSGMOG00000000364	4.39	4.39	0.20	3.11	2.42	−1.02	−0.67	−0.88
ss1712304723	2	0 (NS)	MPG	ENSGMOG00000004225	2.83	2.97	3.40	3.55	>4	−0.98	−0.97	−1.00
ss1712300659	2	0 (NS)	KCNH6	ENSGMOG00000004695	3.73	3.90	−0.89	0.63	1.29	−1.02	0.11	−0.96
ss1712296497	2	1,603 (I)	IQCK	ENSGMOG00000018312	4.62	4.69	0.36	3.22	3.00	−1.02	−0.12	−0.93
ss1712296037	2	82	HEATR2	ENSGMOG00000004347	3.73	3.62	0.06	3.00	1.02	−1.02	−1.01	−0.90
ss1712304682	2	0 (NS)	GRN	ENSGMOG00000000761	3.16	3.18	−0.36	1.84	1.17	−1.02	−0.95	−0.96
ss1712298365	2	109 (I)	COL1A1A	ENSGMOG00000008472	5.75	5.61	−0.05	1.80	1.11	−0.94	−0.94	−0.96
ss1712296495	2	7	CDK6	ENSGMOG00000017917	4.30	4.28	0.20	3.16	2.15	−1.02	−0.81	−0.89
ss1712298100	2	0	CDC27	ENSGMOG00000013717	4.60	4.32	0.62	2.30	0.57	−1.01	−1.03	−0.97
ss1712301880	2	1,421	APOL4	ENSGMOG00000009453	2.51	2.55	0.92	2.17	1.80	−1.02	−1.02	−0.99
ss1712298561	2	43 (I)	ADCY9	ENSGMOG00000014968	4.47	4.53	0.20	3.31	2.60	−1.02	−0.40	−0.88
ss1712295765	3	11 (I)	SERPING1	ENSGMOG00000013214	3.28	3.21	1.41	2.66	1.55	−1.02	−1.02	−1.00
ss1712299811	3	0	LITAF	ENSGMOG00000014597	2.52	2.79	3.22	1.78	>4	0.70	−0.65	−1.01
ss1712304088	3	0 (NS)	LITAF	ENSGMOG00000014597	2.48	2.71	2.96	1.77	>4	−0.11	−0.81	−1.01
ss1712304766	3	0 (NS)	DCHS2	ENSGMOG00000000621	3.08	3.17	0.91	1.81	2.49	−1.02	−0.93	−1.00
ss1712304767	3	0 (NS)	DCHS2	ENSGMOG00000000621	2.79	2.93	0.11	1.43	2.37	−1.00	−0.81	−1.00
ss1712300044	4	3,371	XKR6	ENSGMOG00000011568	3.30	3.22	2.03	1.96	1.23	−1.01	−1.01	−1.00
ss1712297689	4	387 (I)	PSD4	ENSGMOG00000014165	3.83	3.84	2.18	2.16	2.36	−1.03	−1.02	−1.01
ss1712301665	4	4,026 (I)	MGAT5	ENSGMOG00000018588	2.62	2.50	1.90	1.11	−0.16	−1.00	−1.02	−1.01
ss1712297584	4	56	CCNT2B	ENSGMOG00000006401	3.31	3.23	2.20	1.97	1.22	−1.02	−1.01	−1.02
ss1712305126	7	0	ATP1A1	ENSGMOG00000005261	4.12	4.00	1.34	1.71	1.18	−1.01	−1.02	−0.99
ss1712303106	8	4,358	ZP2L1	ENSGMOG00000010293	3.66	3.66	>4	>4	3.40	−1.04	−1.02	−1.02
ss1712303105	8	2,840	PER2	ENSGMOG00000010256	3.45	3.35	>4	>4	2.54	−1.01	−1.00	−1.00
ss1712299931	8	219	MCOLN3	ENSGMOG00000011255	2.91	2.80	2.01	>4	1.58	−1.02	−1.01	−0.99
ss1712304112	8	0 (NS)	MCOLN3	ENSGMOG00000011255	2.64	2.44	1.63	3.22	−0.44	−1.02	−1.03	−0.99
ss1712302864	8	1,986	KISS1R	ENSGMOG00000011215	3.84	4.06	2.60	2.35	>4	−0.46	−0.58	−1.01
ss1712299109	9	0	PSMA1	ENSGMOG00000013869	4.02	4.11	>4	>4	>4	−0.98	−0.94	−1.01
ss1712303592	9	0 (NS)	ITPR2	ENSGMOG00000003496	3.58	3.51	2.13	0.58	0.91	−1.02	−1.02	−0.99
ss1712302672	11	0 (NS)	TLR18	ENSGMOG00000003793	2.64	2.56	1.78	2.12	0.69	−1.02	−1.01	−1.01
ss1712303453	14	0 (NS)	RNF128	ENSGMOG00000008481	3.12	3.07	>4	2.47	2.62	−1.03	−1.01	−0.98
ss1712301982	15	0	FBXW4	ENSGMOG00000015043	2.73	2.78	1.68	3.05	2.89	−1.02	−0.99	−1.00
ss1712299551	16	4,571	KCNJ5	ENSGMOG00000013278	3.50	3.63	2.08	0.55	2.25	−0.86	−1.01	−0.99
ss1712304231	21	0 (NS)	SLC22A16	ENSGMOG00000002642	3.57	3.73	0.87	1.32	2.60	−0.93	−0.81	−0.99
ss1712302811	21	188	ESR1	ENSGMOG00000014898	2.45	2.74	0.44	0.47	>4	0.45	0.65	−1.01

Note.—I, intron; NS, nonsynonymous. Distance to gene is in base pairs. Values in the BAYENV2 and BAYESCAN columns have negative values written in red color and log10(BF/PO) > 1 values are in bold text. HGNC, gene names according to the HUGO gene nomenclature committee. Gene ID refers to the annotated Atlantic cod genome, available at www.codgenome.no and Ensembl.

Discussion

We identified a set of directionally selected loci with allele frequencies strongly correlated with habitat differences in salinity, oxygen, and temperature. These loci reside within several discrete regions within the Atlantic cod genome, suggesting that these regions are influenced by divergence hitchhiking and indicate presence of genomic islands of divergence.

Patterns of Genomic Diversification Are Caused by Adaptation

Neutral and selective evolutionary forces shape the genetic makeup among populations and it is important to disentangle these effects. Using the set of presumed neutral SNPs we show that the Baltic population is divergent relative to the other three populations (table 2 and supplementary fig. S1, Supplementary Material online). Although effects of selection on some of these SNPs cannot be totally excluded, our results are in line with microsatellite data, showing low levels of genetic differentiation among North Sea, Kattegat and Öresund populations compared with the Baltic population (Nielsen et al. 2003) which could possibly be a secondary effect of local adaptation to the Baltic conditions. The Baltic cod population is on the margin of its ecological and geographical range (Johannesson and André 2006), suggesting that the observed divergence between the Baltic population and the other populations could be partly caused by decreased population density at its distribution margin (Bridle and Vines 2007; Orsini et al. 2013). This is less likely however, as the effective population size estimates of Baltic cod still remain relatively high (Poulsen et al. 2006). Poulsen et al. (2006) estimated Ne for Baltic and North Sea cod using temporal analyses of allele frequencies in microsatellites, and found Ne in the Baltic population to be only slightly smaller than in the North Sea population (both in the thousands), but still large enough to not loose its evolutionary potential (Franklin 1980). Poulsen et al. (2006) concluded by urging more attention to selection than to bottleneck effects on genetic composition of Baltic cod. Hence, there are strong indications that Ne in the Baltic population is sufficiently large and that bottleneck effects should not be a major issue, supporting the view that demographic processes alone are not the ultimate driver for the outlier patterns that we observe. Rather, a scenario where adaptation drives neutral structure and “isolation by adaptation” or “isolation by environment“ is more likely (Orsini et al. 2013; Wang and Bradburd 2014).

Outlier Detection and Methodological Considerations

The consistency of the different approaches for outlier detection (supplementary fig. S3, Supplementary Material online) and the functional relevance of the underlying genes strongly suggest that the majority of the identified loci and their associated genomic regions is subject to divergent selection (see however Bierne et al. 2013). Nevertheless, it is inherently difficult to define a biologically meaningful cutoff value in the various outlier tests, and identified outliers remain candidates. Here, we attempted to achieve robust conclusions by combining cutoff values based on log10(BF/PO) values (Jeffreys 1961) with q values and empirical P values (for a detailed discussion, see supplementary text S4, Supplementary Material online).

It has been suggested that outlier tests may have high false positive rates due to the effects of population structure and bottlenecks (see, e.g., Narum and Hess 2011; de Villemereuil et al. 2014; Lotterhos and Whitlock 2014). One way to reduce the effect of population structure is to perform outlier analyses between pairs of populations (cf. Vitalis et al. 2001). Another way to compensate for population structuring is by calculating empirical P values based on a set of putatively neutral SNP markers (determined a priori) to create a null distribution to test all SNP markers against (cf. Lotterhos and Whitlock 2014). We employed both of these methods to minimize the effects of population structuring in the data (supplementary text S4, Supplementary Material online).

Recent simulation studies have suggested that correlation-based approaches in many instances outperform more traditional population genomic approaches in accurately identifying loci under divergent selection (de Mita et al. 2013; de Villemereuil et al. 2014; Lotterhos and Whitlock 2014). By also analyzing our data using BAYENV 2.0, which accounts for demographic signals, we were able to disentangle the genetic signals caused by selection rather than demography with more confidence (details in supplementary text S4, Supplementary Material online).

Genomic Divergence Implies Ongoing Speciation

In the process of population divergence and speciation, heterogeneous genomic divergence can be formed (Nosil et al. 2009; Smadja and Butlin 2011), observed as elevated levels of divergence in selected regions, surrounded by physically linked loci through divergence hitchhiking (cf. Charlesworth et al. 1997; Via 2009, 2012). Divergent selection may further promote extrinsic reproductive isolation where migrants between different environments have reduced fitness (Feder, Egan, et al. 2012). For Atlantic cod, ecological adaptation seems to restrict gene flow across the salinity barrier, and in the Baltic Sea population we observe a genome-wide pattern of divergence, with SNPs that are candidates for selection in all LGs except LG23 (fig. 4b–d). This is consistent with the theory of isolation by adaptation and the later stages of ecological speciation where gene flow is small or nonexistent (Feder, Egan, et al. 2012). Further, most of these outlier SNPs cluster in genomic regions in high LD (fig. 3 and supplementary table S3, Supplementary Material online) suggesting that mutations in close proximity (high LD) hitchhike with the selected genomic region, increasing the size of the local genomic islands (Smith and Haigh 1974; Feder, Egan, et al. 2012; Feder, Gejji, et al. 2012). However, recombination will reduce LD across the genome over time, resulting in larger patterns of genome-wide divergence, although heterogeneity among regions may still be present due to varying degree of selection and recombination. Eventually, ecologically favored alleles will predominate in one habitat and neutral- and universally favored alleles will potentially be present in all habitats (see, e.g., Savolainen et al. 2013). At this point, diagnostically fixed differences can be observed between populations (Feder, Egan, et al. 2012; Seehausen et al. 2014). In our data set, fixed allele frequencies were detected in 19 SNPs potentially under selection in the Baltic population and 59 SNPs were nearly fixed (frequency > 0.95), suggesting strong ongoing diversification. These results are comparable to those in other fish species such as three-spine sticklebacks (Gasterosteus aculeatus) (Roesti et al. 2012, 2014) and lake whitefish (Coregonus clupeaformis) (Gagnaire et al. 2013) for which similar patterns of divergence have been described.

Among the North Sea, Kattegat, and Öresund populations, where obvious environmental differences are small, we observe a pattern consistent with the early stages of genetic divergence, where a substantial amount of gene flow is still prevalent, and one would expect the presence of relatively few but potentially large genomic islands of divergence (Feder, Egan, et al. 2012). In these populations, we detected regions that are subject to divergent selection in few LGs (fig. 4e–g) compared with the Baltic population and most of these SNPs reside in LG12, in high LD with each other (fig. 3 and supplementary table S3, Supplementary Material online). The finding of selection in relatively few LGs agrees with previous findings in Atlantic cod (Bradbury et al. 2010, 2013; Hemmer-Hansen et al. 2013; Karlsen et al. 2013). This pattern of early divergence has also been shown in the studies of pea aphids (Acyrthosiphon pisum pisum), African malaria mosquitoes (Anopheles gambiae), and Heliconius butterflies (Turner et al. 2005; Via and West 2008; Via 2009; Nadeau et al. 2012; Via et al. 2012). In some instances, studies of the same species under different environmental conditions have arrived at different results regarding the distribution of the divergence pattern (see, e.g., Hohenlohe et al. 2011; Roesti et al. 2012; Gagnaire et al. 2013). This seems to be the case also in our investigation, where the observed divergence patterns for the Baltic population comparisons (fig. 4b–d) differ markedly to that observed in the North Sea, Kattegat, and Öresund populations (fig. 4e–g). In some species strong signatures of selection have been shown for chromosome inversions, although few studies have identified the actual target genes for selection (Kirkpatrick and Kern 2012). Nevertheless, our current data cannot accurately resolve if the clustering of outliers in LD regions is due to inversions.

Genomic Regions under Divergent Selection in the Baltic Sea

Studying the genetic changes that contribute to reproductive isolation in partly reproductively isolated populations (before they become confounded by additional genetic differences after speciation is complete) may reveal important aspects of the speciation process (Via 2009). Our sampling design enabled us to do so and our results indicate that several discrete regions of the Atlantic cod genome are candidates for directional selection and most are associated with adaptation to the local conditions in the Baltic Sea. It is likely that salinity and oxygen levels are strong evolutionary forces for Baltic cod (figs. 4 and 5) and our results agree with findings that natural selection can shape population structure on short spatial scales, despite the high dispersal capacity of marine organisms (see, e.g., Gaggiotti et al. 2009; Bradbury et al. 2010, 2013; André et al. 2011; Lamichhaney et al. 2012; Nielsen et al. 2012; Corander et al. 2013; Defaveri, Shikano, et al. 2013; Hemmer-Hansen et al. 2013; Vandamme et al. 2014).

During speciation in the presence of gene flow, the establishment and maintenance of genomic regions that sufficiently can resist gene flow is only likely if divergent selection (or sexual selection) is strong and hence the initial barriers to gene flow are likely to evolve quickly (Via 2001; Hendry et al. 2007). As a result of such strong selection, the genomic regions causing the reproductive isolation become particularly distinctive relative to the remaining genome, facilitating its discovery in empirical analyses (Via 2009). We detect several distinct genomic regions in which a number of SNPs in high LD appear as candidates for selection. Several of these regions have been associated with environment variables in other studies. For instance, we identify 32 SNPs in linked regions toward the end of LG2 as candidates for selection. The same regions have previously been associated with variation in temperature (Bradbury ). Nevertheless, in our study these SNPs appear more associated with salinity and oxygen (at spawning depth) than with temperature (fig. 5 and supplementary table S4, Supplementary Material online). Similarly, one SNP located 337 bp from the prolactin (PRL, in LG2) gene, lying outside the region defined by Bradbury , was also associated with temperature but more strongly so with salinity and oxygen (supplementary table S4, Supplementary Material online). The gene for hemoglobin beta-1 (Hb-β1) is also located in LG2 (in contig_03031) and identified to be 17.8 cM away from the outlier region by Bradbury . It has been proposed by Andersen , but see also Sick (1961), that the Baltic cod is adapted to the environmental hypoxia and temperature in the Baltic Sea by possessing the high-affinity Val-Ala form of Hb-β1. In the same scaffold as Hb-β1 an outlier in a nonsynonymous SNP (ss1712304723) was identified, which is associated with salinity (at surface and spawning depth) and oxygen level (at spawning depth). If this defined outlier SNP is linked to Hb-β1, it may reflect the hemoglobin polymorphism in the Baltic Sea.

Moreover, we detected 17 SNPs as candidates for selection in LG1. The same region has been associated with a migratory ecotype by Hemmer-Hansen and by Karlsen between migratory Northeast arctic cod and stationary Norwegian coastal cod. This region contains the well-known pantophysin gene polymorphism (PanI, in scaffold_09065) that has been used to determine individuals as either stationary or migratory (Fevolden and Pogson 1997; Pogson and Fevolden 2003). This marker is known to be nearly monomorphic for the coastal type (PanIA) in North Sea and Baltic populations (Case et al. 2005) which was also the case for the closest genotyped SNP in this study (ss1712302787, 5,148 bp away; supplementary table S2, Supplementary Material online). Notably, within LG1, all pairwise outlier SNPs relative to the Baltic Sea, and loci environmentally associated with salinity and oxygen (table 4), are located within the same scaffolds (supplementary table S4, Supplementary Material online) as Hemmer-Hansen describes as associated with a migratory ecotype. The fact that we detect outliers in the same region that has previously been identified as associated with a migratory ecotype, but also between the Baltic- and the North Sea (Hemmer-Hansen et al. 2013), in the Baltic population, might indicate that this region is not selected by migratory behavior per se.

Furthermore, in LG12, we observe an outlier pattern driven by the North Sea/Kattegat comparison and to some extent also the North Sea/Öresund comparison (fig. 4). The landscape genomic analysis suggests a temperature association within the long-range LD region on LG12 (figs. 3 and 5) that corresponds to the temperature-associated region identified by Bradbury . This correlation is not as robust as the correlations found for salinity and oxygen level at the other LGs, as decisive association was only detected in eight SNPs (using q < 0.01 as a cutoff). These results suggest that LG12 is not central in the adaptation to the Baltic Sea environment, but rather the selection in this genomic region is driven by adaptation to temperature variation, not encountered by the Baltic Sea cod. In the tightly linked outlier block on LG12, the allele frequencies within the Baltic Sea and within the North Sea are similar (supplementary table S2, Supplementary Material online) possibly reflecting that the detected outlier pattern either predates the Baltic split or is subject to similar selection pressure in both areas. In Hemmer-Hansen it is evident that SNPs in this region shows elevated FST values in a Norwegian coastal cod/North Sea comparison, which may indicate that it is the North Sea that is driving this outlier pattern. However, it is also possible that the Öresund and Kattegat populations are experiencing specific selection pressure in this genomic region or share a common ancestry with the Norwegian coastal cod populations.

Apart from outliers in long LD regions, we also identify single SNPs as significant candidates for selection, associated with at least one environmental variable (figs. 4 and 5 and supplementary table S4, Supplementary Material online) indicating that several smaller genomic regions also play a role in the adaptation to the different environments under study. As expected, we did not identify all previously known outliers regions in Atlantic cod; however, we did identify some novel regions under selection. For instance, in LG7 we detected a single decisive outlier SNP, located within an exon in the Na+/K+-ATPase gene (ATP1A1); see discussion on physiological adaptation to low salinity. Nevertheless, no SNPs under selection were detected in the large LD region on LG7 (fig. 3) which was previously described as associated with temperature and behavior (Bradbury et al. 2010; Hemmer-Hansen et al. 2013).

Bradbury identified strong interchromosomal LD (ILD) among the large LD regions in LG2, 7 and 12, supporting a hypothesis of divergence hitchhiking, transitioning to genome hitchhiking with reproductive isolation. Even though we find some evidence of high ILD between SNPs among the same regions, the most notable findings of ILD are between outlier SNPs on LG1 and 2 and between outlier SNPs on LG2 and 4 (supplementary table S3, Supplementary Material online). Altogether, our results correspond well with earlier work in Atlantic cod: Bradbury and Hemmer-Hansen separately identified outlier SNPs in LG1, 2, 7, and 12 whereas a later study by Bradbury identified outliers in four additional LGs (LG6, 8, 13, and 22).

Outlier SNPs Indicate Genes Associated with Adaptation to Low Salinity at the Egg Stage

In the low saline environment in the Baltic Sea, cod eggs need to be highly hydrated to maintain buoyancy and avoid sinking into the deeper anoxic water layers. As a consequence, Baltic cod eggs have a neutral egg buoyancy at 14.5‰ (Nissling et al. 1994; Nissling and Westin 1997) compared with 33.0‰ for fully marine cod (Thorsen et al. 1996) and Belt Sea and Skagerrak cod (19.6‰ and 26.6‰; Nissling and Westin 1997).

The molecular and cellular mechanisms behind the extreme hydration of the oocytes are poorly understood, but recent studies have revealed an important role of aquaporins in oocyte hydration (Fabra et al. 2005; Raldúa et al. 2008; Cerdà 2009; Cerdà et al. 2013). The discovery of a teleost-specific aquaporin, AQP1ab, which is specifically expressed in the oocyte (Fabra et al. 2005, 2006), and the fact that hydration is inhibited when AQP1ab is blocked (Fabra et al. 2005, 2006) suggests a role for this AQP channel in hydration of oocytes. Further details are given in supplementary text S5, Supplementary Material online. Two outlier SNPs close to a tandem duplicate of AQP1a (on LG8) were associated with salinity and oxygen and were also pairwise outliers between the Baltic and the Kattegat populations, indicating that this gene could be a key player in the adaptation to low saline conditions in Atlantic cod.

We also detected nine SNPs as candidates for selection in solute carrier protein (SLC) genes (supplementary table S4, Supplementary Material online), which is significantly more than expected by chance alone (Fisher’s exact test; P < 0.0001), belonging to seven SLC families (SLC 1, 4, 9, 22, 25, 29, and 35) where five families (SLC 4, 9, 25, 29, and 35) have previously been detected in the ovary transcripts of striped bass (Reading et al. 2012). SLCs consist of a wide range of membrane proteins that control cellular influx and efflux of solutes such as inorganic cations and anions, salts, metals, amino acids, fatty acids, and lipids (Hediger et al. 2004; He et al. 2009; Schlessinger et al. 2010). The characterization of SLC gene expression in growing oocytes and during ovarian maturation has direct importance for understanding oogenesis in teleosts (Reading et al. 2012), especially in the context of egg buoyancy in pelagophilic marine species. It is likely that a repertoire of the SLC genes in Atlantic cod acts in concert with aquaporin genes during oogenesis, supported by Bobe showing an upregulation of SLC26 and AQP4 in rainbow trout (Oncorhynchus mykiss) during ovarian maturation. For a detailed discussion of the SLC genes, see supplementary text S5, Supplementary Material online.

During oocyte maturation, the hydrolysis of yolk proteins generates a pool of free amino acids that enables water influx into the oocyte during oocyte hydration (Finn and Fyhn 2010). Central to this process are the vitellogenins, and indeed, we do find a SNP in vitellogenin (VTG6), a lipid transporter, as a candidate for selection (supplementary table S4, Supplementary Material online) associated with both salinity at spawning depth and temperature. One of the outlier SNPs with strongest association to salinity and oxygen level at spawning depth and a decisive pairwise outlier in all Baltic comparisons was localized close to the zona pellucida glycoprotein-2 gene (ZP2L1) on LG8 (table 4). ZP2L1 shows ovary-specific expression in zebrafish (Danio rerio) (Mold et al. 2001) and is liable for species-restricted binding of sperm to unfertilized eggs (Wassarman 2008). Further details are given in supplementary text S5, Supplementary Material online.

Outlier SNPs Indicate Physiological Adaptation to Low Salinity

After the egg stage, osmoregulation primarily takes place in gills, intestines and kidneys, in which the enzyme Na+/K+-ATPase, an enzyme that is fundamental to osmoregulation and ion exchange, is abundantly expressed. The role of this enzyme within osmoregulatory organs in teleosts is to actively pump K+ in to the cell and Na+ out of the cell across a concentration gradient (Geering 1990) and aiding in other transport processes associated with osmoregulatory function (McCormick 1995). We detected two different Na+/K+-ATPase genes (ATP1A1 and ATP1B1) as candidates for selection (see supplementary text S6, Supplementary Material online). The ATP1A1 gene (in LG7) is decisively associated with salinity and oxygen at spawning depth and surface temperature and is also a pairwise outlier only across the salinity barrier (table 4), highlighting its presumptive osmoregulatory role. This gene has been shown to have significant different gene expression in gill tissue between North Sea- and Baltic cod, suggesting adaption to local salinity conditions (Larsen et al. 2011).

A SNP (ss1712302819) located 337 bp away from the PRL gene, localized on LG2, is one of the SNPs with strongest environmental correlation and is a pairwise outlier in all Baltic comparisons (table 4). Prolactin is central in ion uptake and in regulating water and ion permeability in osmoregulatory surfaces in freshwater and euryhaline fish species (Manzon 2002; Sakamoto and McCormick 2006) and is known to contribute to the Na+/K+-ATPase regulation in for example gilthead seabream (Sparus aurata) (Mancera et al. 2002). It is widely accepted that PRL is involved in freshwater adaptation (Hirano 1986; Manzon 2002), antagonized by growth hormone (GH) and insulin-like growth factor-1 (IGF1)—known to promote acclimatization to seawater in many teleosts (see, e.g., Sakamoto and McCormick 2006). Here, we find significant association between salinity and/or oxygen level within or close to the genes GHRHR, GFBP1, and IGFBP3 (see supplementary text S6, Supplementary Material online), indicating central roles in salinity adaptation.

Conclusions

By combining outlier and landscape genomic analyses, we identified a set of directionally selected loci that are strongly correlated with habitat differences in salinity and oxygen in the Baltic Sea. Relative to the Baltic Sea, we detect an outlier pattern consistent with later stages of ecological divergence (Feder, Egan, et al. 2012) where most of the outlier SNPs cluster in extended genomic regions in high LD, suggesting divergence hitchhiking and the presence of genomic islands of divergence. Fixed and nearly fixed allele frequencies were detected in the Baltic population at a low frequency, suggesting that the process of diversification is ongoing. Candidate SNPs for selection was detected in all LGs, but LG2 clearly stands out with the highest number of outliers with respect to the Baltic Sea and salinity and oxygen association. A large portion of the outlier SNPs detected reside within ecologically important genes affecting egg buoyancy and general osmoregulation, and are thus likely to constitute an evolutionary response to the ecological conditions in the Baltic Sea. Such an adaptive response may contribute to a strong and effective reproductive barrier, leading to ecological speciation in the Baltic cod.

Data Access

Previously published SNPs are referred to by their rs-numbers whereas novel SNPs are referred by their ss-numbers, available in dbSNP (www.ncbi.nlm.nih.gov/SNP, last accessed May 1, 2015).

Supplementary Material

Supplementary texts S1–S6, tables S1–S4, and figures S1–S4 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).