Ancient genomic variation underlies repeated ecological adaptation in young stickleback populations

Abstract

Adaptation in the wild often involves standing genetic variation (SGV), which allows rapid responses to selection on ecological timescales. However, we still know little about how the evolutionary histories and genomic distributions of SGV influence local adaptation in natural populations. Here, we address this knowledge gap using the threespine stickleback fish (Gasterosteus aculeatus) as a model. We extend restriction site-associated DNA sequencing (RAD-seq) to produce phased haplotypes approaching 700 base pairs (bp) in length at each of over 50,000 loci across the stickleback genome. Parallel adaptation in two geographically isolated freshwater pond populations consistently involved fixation of haplotypes that are identical-by-descent. In these same genomic regions, sequence divergence between marine and freshwater stickleback, as measured by d_XY , reaches tenfold higher than background levels and genomic variation is structured into distinct marine and freshwater haplogroups. By combining this dataset with a de novo genome assembly of a related species, the ninespine stickleback (Pungitius pungitius), we find that this habitat-associated divergent variation averages six million years old, nearly twice the genome-wide average. The genomic variation that is involved in recent and rapid local adaptation in stickleback has therefore been evolving throughout the 15-million-year history since the two species lineages split. This long history of genomic divergence has maintained large genomic regions of ancient ancestry that include multiple chromosomal inversions and extensive linked variation. These discoveries of ancient genetic variation spread broadly across the genome in stickleback demonstrate how selection on ecological timescales is a result of genome evolution over geological timescales, and vice versa.

Keywords: Adaptation; RAD‐seq; evolutionary genomics; speciation islands.

Figure 1

Stickleback sampling and RAD sequencing to measure haplotype variation. (A) Threespine stickleback sampling locations in this study. Colors represent habitat type: red: marine; blue: freshwater. Number of haploid genome sampled is shown. (B–D) We modified the original RAD‐seq protocol to generate local haplotypes. Colored bars represent polymorphic sites. For a detailed description of haplotype construction, see Methods. (B) Overlapping paired‐end reads are anchored to PstI restriction sites. (C) Paired reads mapping to each half‐site are merged into contigs. Contigs mapping to the same restriction site are identified by alignment to the reference genome. (D) Sequences from each half of a restriction site are phased to generate a single RAD locus. RAD tags in the background represent multiple genotypes used in phasing.

Figure 2

The genealogical structure of parallel genomic divergence. (A) Genome‐wide F_ST for both marine‐freshwater comparisons was kernel‐smoothed using a normally distributed kernel with a window size of 500 kb. Inverted triangles indicate the locations of two genes known to show extensive marine‐freshwater haplotype divergence, eda and atp1a1. Three chromosomal inversions are highlighted in yellow. (B) Lineage sorting patterns were identified from maximum clade credibility trees for each RAD locus. Blue bars: haplotypes from both freshwater populations form a single monophyletic group; red: haplotypes from the marine population form a monophyletic group; black: A RAD locus is structured into reciprocally monophyletic marine and freshwater haplogroups.

Figure 3

Extensive sequence divergence between marine and freshwater haplogroups accompanies reciprocal monophyly. For each reciprocally monophyletic RAD locus, we calculated sequence variation (π) within and sequence divergence between habitat types (d_XY). Each RAD locus is shown as a pair of lines connecting estimates of π and d_XY. Boxplots show distributions across all reciprocally monophyletic RAD loci: Boxes are upper and lower quartiles, including the median; whiskers extend to 1.5 × interquartile range. Dashed lines are the genome‐wide medians. Single RAD loci from within the transcribed regions of eda and atp1a1 are shown as gold and green lines, respectively, and presented as haplotype networks. Dots represent mutational steps. Circle sizes indicate the number of haplotypes and colors indicate population of origin as in Figure 1. Each network = 28 haplotypes.

Figure 4

Marine‐freshwater divergence has evolved over millions of years, affecting large genomic regions. We performed Bayesian estimation of the time to the most recent common ancestor (T_MRCA) of alleles at threespine stickleback RAD loci. We calibrated coalescence times within threespine stickleback by including a de novo genome assembly from the ninespine stickleback (Pungitius pungitius) and setting threespine‐ninespine divergence at 15 million years ago. (A) Maximum clade credibility RAD gene tree representative of the genome‐wide average T_MRCA. Branches within threespine are colored by population of origin. (B) Kernel‐smoothed densities of T_MRCA distributions for all RAD loci containing a monophyletic group of threespine stickleback alleles (light gray) and those structured into reciprocally monophyletic marine and freshwater haplogroups. (C) The genomic distribution of reciprocally monophyletic RAD loci (black, as in Fig. 2) is associated with increased T_MRCA at a genomic scale. T_MRCA outlier windows (those exceeding 99.9% of permuted genomic windows) are shown as gray bars. Genome‐wide T_MRCA was kernel‐smoothed using a normally distributed kernel with a window size of 500 kb. Inverted triangles indicate the locations of eda and atp1a1. Three chromosomal inversions are highlighted in yellow.