Genomic differentiation in Pacific cod using Pool-Seq

Abstract

Patterns of genetic differentiation across the genome can provide insight into selective forces driving adaptation. We used pooled whole genome sequencing, gene annotation, and environmental covariates to evaluate patterns of genomic differentiation and to investigate mechanisms responsible for divergence among proximate Pacific cod (Gadus macrocephalus) populations from the Bering Sea and Aleutian Islands and more distant Washington Coast cod. Samples were taken from eight spawning locations, three of which were replicated to estimate consistency in allele frequency estimation. A kernel smoothing moving weighted average of relative divergence (F _ST) identified 11 genomic islands of differentiation between the Aleutian Islands and Bering Sea samples. In some islands of differentiation, there was also elevated absolute divergence (d _XY) and evidence for selection, despite proximity and potential for gene flow. Similar levels of absolute divergence (d _XY) but roughly double the relative divergence (F _ST) were observed between the distant Bering Sea and Washington Coast samples. Islands of differentiation were much smaller than the four large inversions among Atlantic cod ecotypes. Islands of differentiation between the Bering Sea and Aleutian Island were associated with SNPs from five vision system genes, which can be associated with feeding, predator avoidance, orientation, and socialization. We hypothesize that islands of differentiation between Pacific cod from the Bering Sea and Aleutian Islands provide evidence for adaptive differentiation despite gene flow in this commercially important marine species.

Keywords: Pacific cod; Pool‐Seq; genomic differentiation; genomics; local adaptation; selection; selection–migration balance.

FIGURE 1

Spawning locations included in this study were named after nearby landmarks. There were three collections from the Bering Sea (Pervenets, Zhemchug, and Pribilof), three from the Aleutian Islands (Near Island, Kiska Island, and Adak Island), as well as a collection from Kodiak Island and the Washington Coast.

FIGURE 2

Principal components analysis for all pools and all (1,944,870) SNP loci. One principal component was the optimal choice to represent the data, which explained 13.3% of the variance; plots are shown along a single axis only.

FIGURE 3

Weighted average F _ST for all SNPs aligned to GadMor2 in 180 kbp windows that overlap with a step size of 20,000 bp between the Bering Sea and Aleutian Islands superpools (upper panel), and the Bering Sea and Washington Coast superpools (lower panel). F _ST outlier windows are colored in red shades. Islands of differentiation identified in Table 2 are highlighted with gray bars. The island on linkage group 12 is not filled because it was not present in the EBS‐WA comparison.

FIGURE 4

Density plot of normalized F _ST within and outside F _ST outlier windows for the Bering Sea vs. Aleutian Islands comparison (EBS‐AI) and the Bering Sea vs. Washington comparison (EBS‐WA).

FIGURE 5

Weighted average F _ST and d _XY (the latter multiplied by 10 for visualization) in the six linkage groups containing F _ST outlier regions: linkage groups 2, 8, 12, 16, 18, and 22, labeled corresponding to F _ST outlier regions in Table 2. Shaded vertical lines represent outlier regions within which linkage groups with F _ST outliers are greater than 0.030 and pairwise d _XY is greater than the genome average (0.000539 × 10, horizontal black line). Blue points represent the upper 5% quantiles for d _XY, and pink dots indicate F _ST outlier windows. Orange triangles represent locations of annotated genes within F _ST outlier windows, which are labeled corresponding to the gene number in Table A5, and vision gene names are listed.

FIGURE 6

Generalized additive model (GAM) fitted between depth (meters) and near‐bottom optical depth data in the Aleutian Islands (AI) and eastern Bering Sea from summer (June–July) 2006–2018. Panels show: (a) GAM fitted mean (line) ± 1 standard error (shading), (b) AI GAM residuals by longitude, mean residual based on LOESS regression (line) ± 1 standard error (shading), and position of Samalga Pass (~169°29′W). Point and line colors denote region (AI and EBS).