The genetic basis for ecological adaptation of the Atlantic herring revealed by genome sequencing

Abstract

Ecological adaptation is of major relevance to speciation and sustainable population management, but the underlying genetic factors are typically hard to study in natural populations due to genetic differentiation caused by natural selection being confounded with genetic drift in subdivided populations. Here, we use whole genome population sequencing of Atlantic and Baltic herring to reveal the underlying genetic architecture at an unprecedented detailed resolution for both adaptation to a new niche environment and timing of reproduction. We identify almost 500 independent loci associated with a recent niche expansion from marine (Atlantic Ocean) to brackish waters (Baltic Sea), and more than 100 independent loci showing genetic differentiation between spring- and autumn-spawning populations irrespective of geographic origin. Our results show that both coding and non-coding changes contribute to adaptation. Haplotype blocks, often spanning multiple genes and maintained by selection, are associated with genetic differentiation.

Keywords: Atlantic herring; evolutionary biology; genetic adaptation; genomics; natural selection.

Figure 1.. Demographic history and phylogeny.

(A) Geographic location of samples. The salinity of the surface water in different areas is indicated schematically. Autumn spawners are marked with an asterisk. (B) Demographic history. Black circles indicate effective population size over time estimated by diCal (Sheehan et al., 2013); estimates are averages from four arbitrarily chosen genomic regions. The grey field is confidence interval ( ± 2 sd), while light grey lines show the underlying estimates from each genomic region. (C) Neighbor-joining phylogenetic tree. The evolutionary distance between Atlantic and Pacific herring was calculated using mtDNA cytochrome B sequences; right panel, zoom-in on the cluster of Atlantic and Baltic herring populations. Colour codes for sampling locations are the same as in Figure 1A. (D) Global distribution of F_ST –values based on 19 populations of Atlantic and Baltic herring. The inset illustrates the tail of the distribution. The mean and median of this distribution are indicated. To reduce the F_ST sampling variance, we only used SNPs with ≥30x coverage in each population. DOI: http://dx.doi.org/10.7554/eLife.12081.003

Figure 1—figure supplement 1.. Population genetics and Q-Q plot.

(A) Average decay of linkage disequilibrium estimated from 30 genomic regions of 50 kb each, covering a total of 40,000 SNPs. (B) Comparison between observed and expected distribution of minor allele frequencies in Atlantic herring. (C) Genome-wide distribution of Tajima’s D values, calculated in 10 kb windows (mean =-0.57 ± 0.01). (D) Q-Q plots for the contrasts Atlantic vs. Baltic and spring- vs. autumn-spawners. The inflation factor is slightly higher than 1 for both these distributions which is expected considering that the average genetic differentiation among populations is not zero (mean F_ST = 0.038, Figure 1D). (E) Tajima’s D values, calculated in 10 kb windows for scaffold 1420 for 16 individually sequenced Atlantic and Baltic herring. The shaded region highlights the putative selective sweep region present in spring-spawning herring (see Figure 4E). DOI: http://dx.doi.org/10.7554/eLife.12081.004

Figure 2.. Genome assembly and annotation.

(A) Phylogeny of ray-finned fishes (Actinopterygii) from the Devonian to the present, time-calibrated to the geological time scale based on Near et al. (2012). Geological abbreviations: C (Carboniferous), CZ (Cenozoic), D (Devonian), J (Jurassic), K (Cretaceous), Ng (Neogene), P (Permian), Pg (paleogene) and Tr (Triassic). Dating of the specific rounds of whole genome duplication is based on Glasauer and Neuhauss (2014). Abbreviations: Ts3R (teleost-specific third round) and Ss4R (salmonid-specific fourth round) of duplication. The number of species with a genome assembly available is marked within parentheses after their group’s name. Atlantic herring belongs to Clupeiformes, the order indicated in red letters. (B) Orthologous gene families across four fish genomes (C. harengus, D. rerio, L. chalumnae and G. morhua). DOI: http://dx.doi.org/10.7554/eLife.12081.005

Figure 2—figure supplement 1.. Schematic overview of the annotation pipeline.

Program versions used are detailed in Materials and methods. DOI: http://dx.doi.org/10.7554/eLife.12081.006

Figure 2—figure supplement 2.. Density plot of the Annotation Edit Distance (AED) score distribution for gene builds rc4 and rc5.

The refined rc5 (red stroke) shows a distribution closer to 0 compared to the evidence build (black stroke), thereby supporting a second pass of ab initio gene modeling that includes in rc5 those genes lacking support from evidence models alone. DOI: http://dx.doi.org/10.7554/eLife.12081.007

Figure 2—figure supplement 3.. Overall read length histogram for the five synthetic long reads (SLR) libraries.

The shape of the histogram looks similar to other published length histograms (Voskoboynik et al., 2013). The sequence provider has filtered out SLRs shorter than 1,500 bp from the dataset. The red bar indicates the frequency of reads longer than the 99.9% percentile of the data (i.e. 10,887 bp). DOI: http://dx.doi.org/10.7554/eLife.12081.008

Figure 2—figure supplement 4.. Read coverage of the assembly with synthetic long reads (SLRs) is uneven and not Poisson-shaped.

The average coverage for SLRs is ~5-6X. In shotgun sequencing, a Poisson distribution of coverage would be expected (Lander and Waterman, 1988). A Poisson distribution (shown as red dots) fitting our data (λ=5) indicates that this would result in most bases having a coverage of 5 and less than 1% of bases not to be covered at all. The observed distribution is quite different (blue bars), with most bases (36.7% of non-N-bases) not being covered at all and a large variance of coverage across the genome. DOI: http://dx.doi.org/10.7554/eLife.12081.009

Figure 2—figure supplement 5.. Phylogeny of endogenous retroviruses (ERVs).

Phylogeny of 62 Clupea harengus ERVs representing 150 detected ERVs, in relation to 92 reference sequences (Genbank and Repbase) rooted on the Cer1 retrotransposon. The phylogeny was constructed using FastTree2 (Price et al., 2010). Sequences and annotations for herring ERVs are presented in Supplementary file 2. DOI: http://dx.doi.org/10.7554/eLife.12081.010

Figure 3.. Genetic differentiation between Atlantic and Baltic herring.

(A) Manhattan plot of significance values testing for allele frequency differences between pools of herring from marine waters (Kattegat, Skagerrak, Atlantic Ocean) versus the brackish Baltic Sea. Lower panel, corresponding plot for scaffold 218 only; both P- and F_ST-values are shown. (B) Neighbor-joining phylogenetic tree based on all SNPs showing genetic differentiation in this comparison (p<10^-10). (C) Comparison of allele frequencies in five strongly differentiated regions. The major allele in the AB1 sample (Atlantic Ocean) was used as reference at each SNP. Lower panel, neighbor-joining tree based on haplotypes formed by 128 differentiated SNPs from scaffold 218. (D) Heat map showing copy number variation partially overlapping the HCE gene. Orientation of transcription is marked with an arrow; the position of SNPs significant in the χ² test is indicated by stars. Population samples and salinity at sampling locations are indicated to the right; abbreviations are explained in Table 2. (E) Strong genetic differentiation between Atlantic and Baltic herring in a region downstream of SLC12A3; statistical significance based on the χ² test is indicated. DOI: http://dx.doi.org/10.7554/eLife.12081.013

Figure 3—figure supplement 1.. Comparison of allele frequencies estimated using pooled whole genome sequencing or by individual genotyping using a SNP chip.

Abbreviations for localities are given in Table 2. DOI: http://dx.doi.org/10.7554/eLife.12081.014

Figure 3—figure supplement 2.. Additional neighbor-joining trees for the contrast Atlantic versus Baltic.

Trees based on individual SNP genotype data from four highly differentiated regions. DOI: http://dx.doi.org/10.7554/eLife.12081.015

Figure 4.. Genetic differentiation between spring- and autumn-spawning herring.

(A) Manhattan plot of significance values testing for allele frequency differences. (B) Neighbor-joining phylogenetic tree based on all SNPs showing genetic differentiation in this comparison (p<10^-10). (C) Comparisons of allele frequencies in four strongly differentiated regions. The major allele in the AB1 sample (Atlantic Ocean) was set as reference at each SNP. Scaffolds 190 and 1420 have been merged in this plot since it was obvious that they were overlapping. *The signal in scaffold s1440 is present ~27 kb upstream of SOX11 and ~46 kb downstream of DCDC2/ALLC. (D) Neighbor-joining tree based on haplotypes formed by 70 differentiated SNPs from scaffold 190/1420; same populations as in Figure 4C. (E) Plot of average heterozygosity, per SNP in 5 kb windows, across scaffold 1420 indicating a selective sweep among spring-spawners in the region marked with vertical hatched lines. Autumn-spawning populations are marked by an asterisk. DOI: http://dx.doi.org/10.7554/eLife.12081.016

Figure 4—figure supplement 1.. Analysis of deviations from Hardy-Weinberg equilibrium using the F_IT statistic in spring- (BÄV), summer- (BÄS) and autumn- (BÄH) spawners from the same locality (Gävle).

The result is based on individual genotyping, using a SNP-chip, of 1500 SNPs that show strong genetic differentiation between spring- and autumn-spawners. The data are presented as a violin plot that shows a box plot in the middle and a rotated kernel density plot on each side. The central rectangle of the box plot spans the first to third quartiles of the distribution, and the ‘whiskers’ above and below the box show the maximum and minimum estimates. The line inside the rectangle shows the median. Mean F_IT and its standard error is shown at the bottom; n represents the number of loci that was polymorphic in that particular population. DOI: http://dx.doi.org/10.7554/eLife.12081.017

Figure 4—figure supplement 2.. Additional neighbor-joining trees for the contrast between spring- and autumn spawning herring.

Trees based on individual SNP genotype data from three highly differentiated regions. DOI: http://dx.doi.org/10.7554/eLife.12081.018

Figure 5.. Nucleotide diversity within and between samples with different ecological adaptations as regards

(A) salinity and (B) spawning time. For each contrast 30 strongly differentiated regions of the genome and 30 control regions showing no significant differentiation were used. The nucleotide diversity within and between populations for the control regions was estimated around 0.3% consistent with the genome average whereas diversity in differentiated regions was significantly higher. BK=Baltic herring, Kalix; AB=Atlantic herring, Bergen; BÄH=autumn-spawning Baltic herring from Gävle; BÄV, spring-spawning Baltic herring from Gävle; see Table 1. The data are presented as box plots; the central rectangle spans the first to third quartiles of the distribution, and the ‘whiskers’ above and below the box show the maximum and minimum estimates. The line inside the rectangle shows the median. DOI: http://dx.doi.org/10.7554/eLife.12081.019

Figure 6.. Analysis of delta allele frequency (dAF) for different categories of SNPs.

(A) dAF calculated for the contrast marine vs. brackish water. (B) dAF calculated for the contrast spring- vs. autumn-spawning. The black line represents the total number of SNPs in each dAF bin and coloured lines represent M values of different SNP types. M values were calculated by comparing the frequency of SNPs in a given annotation category in a specific bin with the corresponding frequency across all bins. DOI: http://dx.doi.org/10.7554/eLife.12081.020