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 dXY, 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.

"freshwater" alleles arose millions of years ago, even though the freshwater populations studied arose much more recently. While it has been hypothesized that other genes in the stickleback genome may share these patterns, large-scale surveys of genomic variation have been unable to test this prediction directly.
Here, we use new sequencing technologies to survey DNA sequence variation across the stickleback genome for patterns like those at the eda gene. We find that nearly every region of the genome associated with marine-freshwater genetic differences shares this pattern to some degree. Moreover, many of these regions are as old or older than eda, stretching back over 10 million years in the past and perhaps even predating the species we now call the threespine stickleback. We conclude that natural selection has maintained this variation over geological timescales and that the same alleles we observe in freshwater stickleback today are the descendants of those under selection in ancient, now-extinct freshwater habitats. Our findings highlight the need to understand evolution on macroevolutionary timescales to understand and predict adaptation happening in the present day.
The mode and tempo of adaptive evolution depend on the sources of genetic variation affecting fitness (Wright 1932;Orr 2005). While new mutation is the ultimate origin of all genetic variation, recent studies of adaptation in the wild have documented adaptive genetic variation that was either segregating in the ancestral population as standing genetic variation (SGV) Domingues et al. 2012;Schrider and Kern 2017), or introgressed from a separate population or species (Huerta-Sánchez et al. 2014;Fontaine et al. 2015). The use of SGV during evolution appears particularly important when dramatic responses to selection occur on ecological timescales, in dozens of generations or fewer . When environments change rapidly, SGV can propel rapid evolution in ecologically relevant traits even in populations of long-lived organisms like Darwin's finches (Grant and Grant 2002), monkeyflowers (Wright et al. 2013), and threespine stickleback fish (Colosimo et al. 2005).
The contribution of SGV to rapid divergence has important consequences for our understanding of evolutionary genetics. Existing genetic variants have evolutionary histories that are often unknown, but which may none-the-less have significant impacts on subsequent adaptation (Kirkpatrick and Barton 2006;Wright et al. 2013). The abundance, genomic distribution, and fitness effects of SGV are themselves the products of evolution (Charlesworth et al. 1993;Colosimo et al. 2005;Kirkpatrick and Barton 2006;Linnen et al. 2009;Stankowski and Streisfeld 2015), and their unknown history raises fascinating questions for the genetics of adaptation in the wild. When did adaptive variants orig-inally arise? How are they structured, across both geography and the genome? Which evolutionary forces shaped their current distribution and how might this history channel future evolutionary change?
Answers to these questions are critical for our understanding of the importance of SGV in nature, as well as our ability to predict the paths available to adaptation on ecological timescales (Wright et al. 2013). Biologists are beginning to probe evolutionary histories of SGV using genome-wide sequence variation across multiple individuals in numerous populations (Pease et al. 2016), but this level of inference has been unavailable for most natural systems because of methodological limitations that remove phase information (e.g., pool-seq: Schlotterer et al. 2014) or produce very short reads (e.g., RAD-seq: Davey et al. 2011).
Here, we investigate the structure and evolutionary history of divergent SGV by modifying the original sheared RAD-seq method to generate ß700 bp haplotypes at tens of thousands of loci sampled across the stickleback genome. This approach allows us to accurately measure sequence variation and estimate divergence times across the genome. By collecting more detailed sequence information at each RAD locus, this approach also provides more accurate estimates of polymorphism and divergence at each locus, and with far smaller sample sizes, compared to traditional short-read methods (Nei 1987 chapters 10 and 13;Wakeley 2009;Cruickshank and Hahn 2014).
SGV has long been postulated to be critical to adaptation in stickleback, and several recent population genomic studies have supported this hypothesis (Hohenlohe et al. 2010;Jones et al. 2012;Feulner et al. 2013;Roesti et al. 2015;Samuk et al. 2017). Marine stickleback have repeatedly colonized freshwater lakes and streams (Bell and Foster 1994b;Jones et al. 2012;Wund et al. 2016), and adaptive divergence in isolated freshwater habitats is highly parallel at the phenotypic (Colosimo et al. 2004;Cresko et al. 2004) and genomic levels (Hohenlohe et al. 2010;Jones et al. 2012; but see Stuart et al. 2017). In addition, analyses of haplotype variation at the genes eda (Colosimo et al. 2005;Roesti et al. 2014) and atp1a1 (Roesti et al. 2014) present two clear results: separate freshwater populations share common "freshwater" haplotypes that are identical-by-descent (IBD), and sequence divergence between the major marine and freshwater haplogroups suggests their ancient origins--perhaps over two million years ago in the case of eda (Colosimo et al. 2005). While intriguing, it is not clear whether the deep evolutionary histories of these loci are rarities or representative of more widespread ancient history across the genome.
To address fundamental questions of genealogical relationships and molecular evolution in stickleback, we utilize the new RAD-seq haplotyping approach to assay genome-wide variation associated with adaptive divergence in two young freshwater ponds, which formed during the end-Pleistocene glacial retreat (c. 12,000 years ago: Francis et al. 1986;Cresko et al. 2004, Fig. 1). In addition, we generated a de novo genome assembly of the sister taxon ninespine stickleback (Pungitius pungitius), allowing us to estimate divergence times for genealogies across the genome. Our results clearly demonstrate that the previous findings of deep evolutionary history based upon candidate loci are not unique but in fact the rule. A suite of adaptive variation structured into distinct marine and freshwater haplotypes that evolved over millions of years structures a deep pool of SGV fueling repeated and rapid evolution in stickleback.

SAMPLE COLLECTION
Wild threespine stickleback were collected from Rabbit Slough (N 61.5595, W 149.2583, n = 5 fish), Boot Lake (N 61.7167, W 149.1167, n = 5 fish), and Bear Paw Lake (N 61.6139, W 149.7539, n = 4 fish) (Fig. 1A). Rabbit Slough is an offshoot of the Knik Arm of Cook Inlet and is known to be populated by anadromous populations of stickleback that are stereotypically oceanic in phenotype and genotype (Cresko et al. 2004;Hohenlohe et al. 2010). Boot Lake and Bear Paw Lake are both shallow lakes formed during the end-Pleistocene glacial retreat. Fish were collected in the summers of 2009 (Rabbit Slough), 2010 (Bear Paw Lake), and 2014 (Boot Lake) using wire minnow traps and euthanized in situ with Tricaine solution. Euthanized fish were immediately fixed in 95% ethanol and shipped to the Cresko Laboratory at the University of Oregon (Eugene, OR, USA). DNA was extracted from fin clips preserved in 95% ethanol using either Qiagen DNeasy spin column extraction kits or Ampure magnetic beads (Beckman Coulter, Inc.) following manufacturer's instructions. Yields averaged 1-2 μg DNA per extraction (ß30 mg tissue). Treatment of animals followed protocols approved by the University of Oregon Institutional Animal Care and Use Committee (IACUC).

SEQUENCING STRATEGY AND RATIONALE
We designed our sequencing to maximize detection of DNA sequence variation and divergence, with the ultimate goal being the estimation of absolute divergence times of marine and freshwater haplogroups. Previous work by us and others using short sequence reads provided clear evidence of changes in relative frequencies of alleles across stickleback populations (Hohenlohe et al. 2010;Roesti et al. 2014;Lescak et al. 2015;Roesti et al. 2015), but could not sufficiently address questions of haplotype ages. We therefore designed a RAD sequencing approach to (1) accurately estimate sequence diversity within and divergence between threespine stickleback ecotypes and (2) recover sufficient RAD loci that map unambiguously to an outgroup genome sequence from the ninespine stickleback that we could confidently compare diversity within threespine stickleback to divergence from the ninespine stickleback. We chose the ninespine stickleback as an outgroup because the threespine-ninespine split is sufficiently old that lineage sorting should be nearly complete (>>4N e , assuming N e <10 6 , Aldenhoven et al. 2010) yet recent enough to facilitate sequence mapping between species (Rastas et al. 2015).
To achieve our aims, we designed a sequencing method to produce phased haplotypes of ß700 bp at each RAD locus ( Fig. 1B-D) and to sample the genome densely enough to identify signatures of selection after the likely dropout of RAD loci without clear mapping in the ninespine stickleback genome. We used the single-digest, sheared RAD approach to limit biases in our estimates of sequence diversity. RAD-seq has known biases due to mutations in restriction sites causing allele dropout (Arnold et al. 2013;Gautier et al. 2013), the potential for which increases with increasing sequence divergence and leads to underestimates of genetic diversity. Diversity estimates are, however, substantially more accurate with sheared RAD-seq compared to other RAD-seq approaches (e.g., double-digest RAD-seq: Peterson et al. 2012). Importantly for the coalescent analyses we present here, such allele dropout is unlikely to affect estimates of overall divergence across the clade of alleles. When in the rare cases it does, the bias is toward underestimation of the divergence age (Arnold et al. 2013), which would make our findings of deep divergence even more striking.
Our sequencing design facilitated accurate inference of sequence variation even with smaller population samples than are typical among population genomic studies. While allele frequency-based statistics like F ST have particularly high variance with small sample sizes (Willing et al. 2012), our study is fortunate to be built upon numerous properly powered, previous population genomic studies in stickleback, including extensive previous work in the three populations we analyze here. The genome-wide patterns of F ST we observed using our new approach closely matched multiple previous studies (Hohenlohe et al. 2010;Jones et al. 2012, Fig. 2A). Because of this extensive body of previous work, we relied on F ST only to draw inference of larger genomic regions containing tens or hundreds of RAD loci. Instead, as stated above, the focus of this work is to extend these previous findings by addressing the ages of allelic divergence. We therefore do not expect the higher variance associated with smaller sample size to qualitatively influence our results. Importantly, estimation of sequence diversity (π) (Nei 1987) and divergence (d XY ) (Nei 1987;Cruickshank and Hahn 2014) at a given locus improves greatly with increases in sequence length. Using equations (10.9) and (13.83) from Nei (1987; Box 1 in Cruickshank and Hahn 2014), the predicted sampling variances in both π and d XY using 700 bp sequences in five individuals are lower than those obtained using standard 100 bp sequences at any sample size (Fig. S1). Therefore, not only is this novel application of RAD-seq ideally suited for our questions, our findings show that this approach may significantly decrease the necessary sample size, and thus resource expenditure, for many population genomic studies.

LIBRARY PREPARATION
To identify sufficient sequence variation at a RAD locus, and to simplify downstream sequence processing and analysis, we took advantage of longer sequencing reads available on newer Illumina platforms and the phase information captured by paired-end sequencing. We generated RAD libraries from these samples using the single-digest sheared RAD protocol from Baird et al. (Baird et al. 2008) with the following specifications and adjustments: 1 μg of genomic DNA per fish was digested with the restriction enzyme PstI-HF (New England Biolabs), followed by ligation to P1 Illumina adaptors with 6 bp inline barcodes. Ligated samples were multiplexed and sheared by sonication in a Bioruptor (Diagenode). To ensure that most of our paired-end reads would overlap unambiguously and produce longer contiguous sequences, we selected a narrow fragment size range of 425-475 bp. The remainder of the protocol was per Baird et al. (Baird et al. 2008). All fish were sequenced on an Illumina HiSeq 2500 using paired-end 250 bp sequencing reads at the University of Oregon's Genomics and Cell Characterization Core Facility (GC3F).

SEQUENCE PROCESSING
Raw Illumina sequence reads were demultiplexed, cleaned, and processed primarily using the Stacks v1.46 pipeline (Catchen et al. , 2013. Paired-end reads were demultiplexed with process shortreads and cleaned using process radtags using default criteria (throughout this article, names of scripts, programs, functions, and command-line arguments will appear in italics). Overlapping read pairs were then merged with fastq-join (Aronesty 2011). Pairs that failed to merge were removed from further analysis. To retain the majority of the sequence data for analysis in Stacks and still maintain adequate contig lengths, merged contigs were trimmed to 350 bp and all contigs shorter than 350 bp were discarded. We aligned these contigs to the stickleback reference genome Glazer et al. 2015) using bbmap v35.69 with the most sensitive alignment settings ("vslow = t"; http://jgi.doe.gov/data-and-tools/bbtools/) and required that contigs mapped uniquely to the reference. We then used the pstacks, cstacks, and sstacks components of the Stacks pipeline to identify RAD-tags and call SNPs using the maximum likelihood algorithm implemented in pstacks, create a catalog of RAD tags across individuals, and match tags across individuals. All data were then passed through the Stacks error correction module rxstacks to prune unlikely haplotypes. We ran the Stacks component program populations on the final dataset to filter loci genotyped in fewer than four individuals in each population and to create output files for sequence analysis. We use the naming conventions of Baird et al. (Baird et al. 2008): A "RAD tag" refers to sequence generated from a single end of a restriction site and the pair of RAD tags sequenced at a restriction site comprises a "RAD locus" (Fig. 1D).
We used the program phase v2.1 (Stephens et al. 2001;Stephens and Scheet 2005) to phase pairs of RAD tags originating from the same restriction site. We coded haplotypes present at each RAD tag, which often contain multiple SNPs, into multiallelic genotypes. This both simplified and reduced computing time for the phasing process. We also performed coalescent simulations to generate, "cut," and rephase haplotypes to demonstrate the high accuracy of this method using sequences and sample sizes similar to those in this study (Fig. S2). Custom Python scripts automated this process and are included as supplementary files. We required that each individual had at least one sequenced haplotype at each tag for phasing to be attempted. If a sample had called genotypes at only one tag in the pair, the sample was removed from further analysis of that locus. The resultant phased haplotypes were uniformly 696 bp in length [(350 bp x 2)--4 bp PstI overlap] with 690 potential variable sites (6 bp PstI motif was invariant) and were used to generate sequence alignments for import into BEAST.
We recovered a total of 236,787 RAD tags after filtering, mapping to 151,813 PstI restriction sites. At 84,974 restriction sites, we recovered and successfully phased adjacent RAD tags (169,948 RAD tags) into single RAD loci. RAD tags with no variable sites were simply concatenated to the adjacent tag to form a single locus. We retained these 84,974 RAD loci for our analysis. For population genetic analyses, inclusion of singleton (i.e. unpaired) RAD tags did not qualitatively change our results. We chose to restrict genealogical analyses to loci of uniform length and to use the same set of loci in analyses of polymorphism and gene tree topologies.

NINESPINE STICKLEBACK GENOME ASSEMBLY
To estimate T MRCA of threespine stickleback RAD alleles, we used the ninespine stickleback (Pungitius pungitius) as an outgroup. RAD sequence analysis, however, relies on the presence of homologous restriction sites among sampled individuals and results in null alleles when mutations occur within a restriction site (Arnold et al. 2013). Because this probability increases with greater evolutionary distance among sampled sequences, we elected to use RAD-seq to only estimate sequence variation within the threespine stickleback. We then generated a contig-level de novo ninespine stickleback genome assembly from a single ninespine stickleback individual from St. Lawrence Island, Alaska (collected by J. Postlethwait) using DISCOVAR de novo revision 52488 (https://software.broadinstitute.org/software/discovar/blog/). We used this single ninespine stickleback haplotype to estimate threespine-ninespine sequence divergence and time calibrate coalescence times within the threespine stickleback. DISCOVAR de novo requires a single shotgun library of paired-end 250-bp sequence reads from short-insert-length DNA fragments. High molecular weight genomic DNA was extracted from an ethanolpreserved fin clip by proteinase K digestion followed by DNA extraction with Ampure magnetic beads. Purified genomic DNA was mechanically sheared by sonication and size selected to a range of 200-800 bp by gel electrophoresis and extraction. We selected this fragment range to agree with the recommendations for de novo assembly using DISCOVAR de novo. This library was sequenced on a single lane of an Illumina HiSeq2500 at the University of Oregon's Genomics and Cell Characterization Core Facility (GC3F: https://gc3f.uoregon.edu/). We assembled the draft ninespine stickleback genome using DISCOVAR de novo. Raw sequence read pairs were first quality filtered and adaptor sequence contamination removed using the program process shortreads, which is included in the Stacks analysis pipeline. We ran the genome assembly on the University of Oregon's Applied Computational Instrument for Scientific Synthesis (ACISS: http://aciss-computing.uoregon.edu).

ALIGNMENT OF RAD TAGS TO THE NINESPINE ASSEMBLY
We included the single ninespine stickleback haplotype into our sequence analyses by aligning a single-phased threespine stickleback RAD haplotype from each locus to the ninespine genome assembly. For those that aligned uniquely (59,254 RAD loci), we used a custom Python script to parse the alignment fields of the output BAM file (Li et al. 2009) and reconstruct the ninespine haplotype by introducing threespine-ninespine substitutions into the threespine RAD locus sequence. The final dataset consists of 57,992 RAD loci that mapped to the 21 threespine stickleback chromosomes and aligned uniquely to the ninespine assembly.

COMMON ANCESTOR
Allelic divergence can occur by multiple modes of lineage sorting during adaptation. To identify patterns of lineage sorting associated with freshwater colonization, we analyzed gene tree topologies at all RAD loci using BEAST v. 1.7 (Drummond and Rambaut 2007;Drummond et al. 2012). We chose BEAST because it coestimates tree topologies and node ages for sequenced genomic loci. BEAST does not explicitly model natural selection, and this may affect divergence time estimates in genomic regions influenced by selection. However, other methods developed to estimate the age of adaptive alleles make assumptions that are likely not relevant to the evolutionary histories we infer here. First, some models assume a recent origin of an adaptive allele compared to adjacent genomic variation (Peter et al. 2012;Ormond et al. 2016)--the opposite of what we describe here--so that measures of variation at linked sites and the decay of linkage disequilibrium can be used to estimate when a sweep began. Selection in the stickleback populations we study likely acted on SGV, as has been supported by previous studies, and we hypothesize that this SGV may be quite old. Therefore, adaptive alleles already existed on distinct haplotype backgrounds, which masks the differences between selected and linked neutral sites.
Second, a recent model developed to infer ages of standing genetic variants assumes that the variant was evolving neutrally at some point during its trajectory through a population (Peter et al. 2012). This assumption is unlikely for many of the variants we detect here, except in the very distant past and for those variants that have evolved recently arose in genomic regions already heavily influenced by selection. Rather, the patterns of haplotype variation we observed in the genomic regions that differentiate marine and freshwater populations reflect long-term maintenance and isolation of separate haplogroups that mimics population structure and even speciation, with selective sweeps being important but constituting a small minority of the time these haplotypes have segregating in the stickleback metapopulation. For all of these reasons, we therefore chose to estimate tree topologies and divergence times with BEAST, which makes minimal assumptions regarding specific evolutionary processes.
We used blanket parameters and priors for BEAST analyses across all RAD loci. Markov chain Monte Carlo (MCMC) runs of 1,000,000 states were specified, and trees logged every 100 states. We used a coalescent tree prior and the GTR+ substitution model with four rate categories and uniform priors for all substitution rates. We identified evidence of lineage sorting by using the program treeannotator v1.7.5 to select the maximum clade credibility (MCC) tree for each RAD locus and the is.monophyletic() function included in the R package "ape" v3.0 (Paradis et al. 2004;Popescu et al. 2012). We determined for each MCC tree whether tips originating from marine (RS) or freshwater (BL + BP) formed monophyletic clades.
To convert node ages estimated in BEAST into divergence times, in years, we assumed a 15 million-year divergence time between threespine and ninespine stickleback at each RAD locus (Aldenhoven et al. 2010). The T MRCA of all alleles in each gene tree was set at 15 Mya and each node age of interest was converted into years relative to the total height of the tree. Additionally, to use the ninespine stickleback as an outgroup, we required that threespine stickleback haplotypes at a RAD locus were monophyletic to the exclusion of the ninespine haplotype. Doing so reduced our analysis to 49,672 RAD loci for analyses included in Fig. 4. We used medians of the posterior distributions as point estimates of T MRCA for each RAD locus. Because of the somewhat limited information from any single RAD locus, and because the facts of the genealogical process mean that the true T MRCA at any locus likely differs from the 15 My estimate (Kingman 1982a, b;Tajima 1983), we do not rely heavily on T MRCA estimates at individual RAD loci. Rather, we use these estimates to understand broad patterns of ancestry throughout the threespine stickleback genome.
We determined T MRCA outlier genomic regions by permuting and kernel smoothing the genomic distribution of T MRCA estimates using the same window sizes as we present in the main text. Windows where the smoothed T MRCA exceeded 99.9% of permuted windows were considered outliers. This method controls for the local density of RAD loci (poorly sampled regions will have larger confidence bands) and the size of the windows used.

SEQUENCE DIVERSITY AND HAPLOTYPE NETWORKS
We quantified sequence diversity within and among populations and sequence divergence between populations using R v3 (R Core Team 2016). We used the R package "ape" (Paradis et al. 2004) to compute pairwise distance matrices for all alleles at each RAD locus and used these matrices to calculate the average pairwise nucleotide distances, π, within and among populations along with d XY , the average pairwise distance between two sequences using only across-population comparisons (Nei 1987). We also calculated the haplotype-based F ST from Hudson et al. (1992) implemented in the R package "PopGenome" v2.2.4 (Pfeifer et al. 2014). We used permutation tests written in R to identify differences in variation within-and between-habitat type at divergent RAD loci versus the genome-wide distributions. Mann-Whitney-Wilcoxon tests implemented in R were used to identify variation in genome-wide diversity among populations and habitat types.
We constructed haplotype networks of the RAD loci at eda and atp1a1 using the infinite sites model with the function haploNet() in the R package "pegas" (Paradis 2010 Glazer et al. (2015)]). The eda network spans exon 2 and portions of introns 1 and 3 of eda (chr4: 12,808,396-12,809,030).

CODE AVAILABILITY
Scripts used to phase RAD-tags, summarize gene trees, calculate population genetic statistics, and produce figures and statistics presented in article are available at https://github.com/ thomnelson/ancient-divergence. Scripts for processing raw sequence data are available from the authors upon request.

DATA AVAILABILITY
Raw sequence data supporting these findings are available on NCBI at PRJNA429207. The Pungitius pungitius draft assembly is available at PRJNA429208. The final datasets needed to reproduce the figures and statistics presented in the article are available at https://github.com/thomnelson/ancientdivergence.

Results and Discussion
Parallel adaptation to freshwater environments has been a major theme of stickleback evolutionary history (Bell and Foster 1994a). Stereotypical morphological changes to, for example, bony armor (Colosimo et al. 2004) and craniofacial structures (Kimmel et al. 2005) presumably reflect adaptation to similar selective regimes (Reimchen 1994;Arnegard et al. 2014). These phenotypic changes are accompanied by parallel genomic divergence (Hohenlohe et al. 2010;Jones et al. 2012), which involves large regions spanning many megabases (Schluter and Conte 2009;Roesti et al. 2014), including multiple chromosomal inversions . The leading hypothesis for the genetics of parallel divergence in stickleback posits that distinct freshwater-adaptive haplotypes that are identical-by-descent (IBD) are shared among freshwater populations due to historical gene flow between marine and freshwater populations (Schluter and Conte 2009). We tested for the presence of these haplotypes directly and at a genomic scale.

GENOME-WIDE SUITE OF HAPLOTYPES
Our sequencing strategy produced 57,992 RAD loci, each with 690 potential variable sites, present across the three threespine stickleback populations and aligned to the ninespine stickleback genome assembly. These data comprise over 40 Mb of sequence, or nearly 10% of the threespine stickleback genome (9.5% of 419 Mb assigned to chromosomes) Glazer et al. 2015). All loci we recovered were polymorphic and we observed a median of seven segregating sites per locus (range: 2-155, Fig. S3, Table S1). By including haplotypes from all three populations in these genealogical analyses, we were able to jointly calculate population genetic statistics (F ST , π, d XY ) and identify patterns of identity-by-descent (IBD) among populations, which we defined as haplotypes from two populations forming a monophyletic group to the exclusion of the third population.
We find that parallel population genomic divergence in the two freshwater pond populations consistently involved haplotypes that were identical-by-descent (IBD) among both freshwater populations (Fig. 2). Background F ST between populations was moderate, ranging from 0.139-0.226 (F ST(RS-BL) = 0.139, F ST(RS-BP) = 0.194, F ST(BL-BP) = 0.226; two-sided Mann-Whitney test for all pairwise comparisons: P ࣘ 1 × 10 −10 ). Genome-wide differentiation was highest between the freshwater populations, supporting the hypothesis of independent colonization and in agreement with previous morphological (Francis et al. 1986) and genetic data (Hohenlohe et al. 2010). The genomic distributions of marine-freshwater F ST were similar to those previously reported (Hohenlohe et al. 2010), including outlier regions over a broad span of chromosome 4 in which the eda gene is embedded (orange triangle in Fig. 2A), and three regions now known to be associated with chromosomal inversions on chromosomes 1, 11, and 21 . The gene atp1a1 (green triangle in Fig. 2A) is contained within inv1. As expected, we found distinct haplogroups associated with marine and freshwater habitats at both eda and atp1a1 (Fig. 3, insets). Strikingly, this finding of habitat-specific haplogroups was not at all unique to these well-studied genes or chromosomal inversions. The two isolated freshwater populations shared IBD haplotypes within all common marine-freshwater F ST peaks even though IBD was rare elsewhere (Fig. 2B). Furthermore, we observed a separate clade of haplotypes representing the marine RS population at the majority (1129 of 2172, 52%) of RAD loci showing freshwater IBD. The result was a genome-wide pattern of reciprocal monophyly between marine and freshwater haplotypes. Notably, this is the same genealogical structure previously reported at eda (Colosimo et al. 2005;Roesti et al. 2014) and atp1a1 (Roesti et al. 2014), demonstrating that these loci are but a small part of a genome-wide suite of genetic variation sharing similar habitat-specific evolutionary histories, and the previous documentation of their genealogies was a harbinger of a much more extensive pattern across the genome revealed here. Here-after, we refer collectively to this class of RAD loci as "divergent loci."

DIVERGENCE INVOLVES ANCIENT ALLELIC ORIGINS
Because the genealogical structure of divergence across the genome mirrors that at eda and atp1a1, we asked whether levels of sequence variation and divergence also showed consistent genomic patterns. At all RAD loci we therefore calculated π within each population, as well as in the combined freshwater populations, and d XY between marine and freshwater habitat types. Genome-wide diversity was similar across populations and habitat types (mean π RS = 0.0032, π BL = 0.0034, π BP = 0.0026, π FW = 0.0038) and comparable to previous estimates (Hohenlohe et al. 2010). Likewise, genome-wide d XY among habitat types was modest (0.0049) when compared to π across all populations (π = 0.0042, two-sided Mann-Whitney test: P ࣘ 1 × 10 −10 ; Fig. S4). Among divergent loci, however, we observed reductions in diversity in both habitats (mean π RS-divergent = 0.0012, π FW-divergent = 0.0016, two-sided permutation test: P ࣘ 1 × 10 −4 , Fig. 3), indicating natural selection in both habitats. Sequence divergence associated with reciprocal monophyly was striking, averaging nearly three times the genome-wide mean (mean d XY-divergent = 0.0124). This divergence ranged more than an order of magnitude (0.0013-0.0442), from substantially lower than the genome-wide average to ten times greater than the average. These findings indicate that much of the genetic variation underlying adaptive divergence is vastly older than the diverging freshwater populations themselves. Not only was adaptive variation standing and structured by habitat, but it has been segregating and accumulating for millennia.
These data clearly support the hypothesis of Schluter and Conte (2009) of ancient haplotypes "transported" among freshwater populations. Much of the divergence we observed was ancient in origin, with levels of sequence divergence at some RAD loci exceeding that observed at eda (Fig. 3, gold line) and suggestive of divergence times of at least two million years ago (Colosimo et al. 2005). Our observation that sequence variation was consistently reduced in both habitat types emphasizes that alternative haplotypes at these loci are likely selected for in the marine population as well as the freshwater. These alternative fitness optima-driven by divergent ecologies-provide a favorable landscape for the maintenance of variation (Charlesworth et al. 1997;Lenormand 2002), but also lead to a more potent barrier to gene flow among freshwater populations if there are fitness consequences in the marine habitat for stickleback carrying freshwater-adaptive variation. Conditional fitness effects through genetic interactions (i.e., dominance or epistasis: Otto and Bourguet 1999;Phillips 2008)  Adaptive divergence between marine and freshwater stickleback genomes is likely ongoing, with recently derived alleles arising on already highly divergent genomic backgrounds. We found reciprocal monophyly associated with a spectrum of sequence divergence, including a substantial fraction of divergent loci (11.0%, 124/1129) with d XY below the genome-wide average. Thus, ongoing marine-freshwater ecological divergence may continue to yield additional marine-freshwater genomic divergence. Moreover, while this younger variation is shared between the freshwater populations in this study, and localizes to genomic regions of divergence shared globally , some adaptive variants may be distributed only locally (e.g., limited to southern Alaska or the eastern Pacific basin). Global surveys of shared variation have been performed ), but future work in this system should quantify the distributions of locally or regionally limited genomic variation involved in ecological divergence, because regional pools of variation may contribute substantially to stickleback genomic and phenotypic diversity (Stuart et al. 2017).

OLD AS THE THREESPINE STICKLEBACK SPECIES
Sequence divergence provides an important relative, but ultimately incomplete, evolutionary timescale. To more directly compare the timescales of ecological adaptation and genomic evolution, we translated patterns of sequence variation into the time to the most recent common ancestor (T MRCA ) of allelic variation, in years. We find that the divergence of marine and freshwater haplotypes has been ongoing for millions of years and extends back to the split with the ninespine stickleback lineage (Fig. 4B). Genome-wide variation averaged 4.1 MY old, and T MRCA for the vast majority of RAD loci was under 5 MY old. In contrast, divergence times at habitat-associated loci averaged 6.4 MYA and, amazingly, the most ancient 10% (118 of 1129) are each estimated at over 10 MY old.
This deep genomic divergence not only underscores that local adaptation to marine and freshwater habitats has been occurring throughout the history of the threespine stickleback lineage--for which there is evidence in the fossil record going back 10 million years (Bell et al. 1985)--but it also demonstrates that at least some of the variation fueling those ancient events has persisted until the present day. In fact, the most highly divergent regions of the threespine stickleback genome also had the lowest rates of monophyly of threespine stickleback haplotypes (Fig. S6). This is consistent with marine-freshwater allelic divergence occurring near to or before the divergence of the threespine and ninespine lineages. It is tempting to think that marine-freshwater ecological divergence, which has evolved convergently in both species and involves similar phenotypes, may involve alleles that are IBD and arose before the split of these species. However, QTL analyses of ecologically relevant armor traits point either to independent genetic architectures in the two species (Shapiro et al. 2009) or to a locus associated with repeated mutational events among threespine stickleback populations (Chan et al. 2010;Shikano et al. 2013). Nevertheless, in some genomic regions marine and freshwater threespine stickleback are as divergent as threespine and ninespine stickleback, which are classified into separate genera.

LONG-TERM DIVERGENCE MAINTAINS LINKED VARIATION AND PROMOTES GENOMIC STRUCTURAL EVOLUTION
Adaptive divergence has impacted the history of the stickleback genome as a whole (Fig. 4C). We identified 32.6 Mb, or 7.5% of the genome, as having elevated T MRCA (gray boxes in Fig. 4C; two-sided permutation test of smoothed genomic intervals, P ࣘ 0.001). Outside of the nonrecombining portion of the sex chromosome (chr. 19), the oldest regions of the stickleback genome were those enriched for divergent loci. Patterns of ancient ancestry closely mirrored recent divergence in allele frequencies ( Fig. 2A) and it appears that historical and contemporary marine-freshwater divergence has impacted ancestry across much of the length of some chromosomes. Chromosome 4, for example, contains at least three broad peaks in T MRCA and a total of 5.9 Mb identified as genome-wide outliers (two-sided permutation test, P ࣘ 0.001). This chromosome has been of particular interest because of its association with a number of phenotypes (Colosimo et al. 2004;Miller et al. 2014), including fitness ). We found the major-effect armor plate locus eda comprised a local peak (mean T MRCA = 6.4 MYA) nested within a large region of deep ancestry spanning 8.1 Mb. Moreover, at least two other peaks distal to eda, centered at 21.4 Mb and 26.6 Mb, were also several million years older than the genomic average at 6.8 MYA and 7.0 MYA, respectively.
Intriguingly, genomic regions of elevated T MRCA remained outliers even after removing marine-freshwater relative divergence outlier loci (as measured by F ST : Fig. S5). We estimated that 7.5% of the genome had increased T MRCA even though only 1.9% of RAD loci (1129 of 57,992) were classified as divergent based on marine-freshwater reciprocal monophyly. When we removed these loci, along with loci with elevated marine-freshwater differentiation (F ST > 0.5), many of the regions in which they resided were still T MRCA outliers. It is possible that the remainder of this old variation is neutral with respect to fitness. However, we identified divergence outliers based on only a single axis of divergence: the marine-freshwater axis. Throughout the entire species range, populations are locally experiencing multiple axes of divergence, including lake-stream and benthic-limnetic axes (McKinnon and Rundle 2002), that often shares a common genomic architecture (Deagle et al. 2012;Roesti et al. 2015). Our data may indicate underlying similarities in selection regimes. Alternatively, this colocalized ancient variation may represent the accumulation of adaptive divergence along multiple axes in the same genomic regions, whether or not the underlying adaptive variants are the same. Aspects of the genomic architecture--such as gene density, local recombination rates, and segregating structural variants-may in part govern where in the genome adaptive divergence can occur (Feulner et al. 2013;Roesti et al. 2013;Aeschbacher et al. 2017;Samuk et al. 2017). Multiple axes of divergence may therefore act synergistically to maintain genomic variation across the stickleback metapopulation. Nevertheless, much of the ancient variation we observe may in fact itself be neutral, having been maintained by close linkage to loci under divergent selection between the marine and freshwater habitats (Barton and Bengtsson 1986;Charlesworth et al. 1997). Indeed, the broadest peaks of T MRCA we observe occur in genomic regions with low rates of recombination in other stickleback populations (Roesti et al. 2013;Glazer et al. 2015), which would extend the size of the linked region affected by divergent selection. On ecological timescales, low recombination rates in stickleback are thought to promote divergence by making locally adapted genomic regions resistant to gene flow (Roesti et al. 2013). Our results potentially extend the inferred impact of recombination rate variation on genomic variation to timescales that are 1000-fold longer, maintaining both multimillion-year-old adaptive variation and large stores of linked genetic variation. Future modeling efforts will be needed to explore the range of population genetic parameter values (e.g., selection coefficients, migration rates, and recombination rates) required to produce the extent of divergence we see here.
Lastly, our findings further support the hypothesis that known chromosomal inversions maintain globally distributed, multilocus haplotypes. The three chromosomal inversions known to be associated with marine-freshwater divergence Roesti et al. 2015) (inv1, inv11, and inv21; yellow bars in Fig. 4C) all showed sharp spikes in T MRCA . Genomic signatures of these inversions are distributed throughout the species range, including coastal marine-freshwater population pairs in the Pacific and Atlantic basins ) and inland lake-stream pairs in Switzerland (Roesti et al. 2015). Despite our limited geographic sampling, our finding that all three of these inversions are over six million years old is further evidence of single, ancient origins of each, followed by their spread across the species range. Each inversion contained a high density of divergent RAD loci (inv1: 64% of loci divergent; inv11: 60%; inv21: 71%) but we also identified regions within these inversions in which haplotypes from marine or freshwater habitats, or both, were not monophyletic. inv1 and inv11 both contained two regions separated by loci in which neither habitat type was monophyletic; inv21, the largest of the three, contained ten such regions. Additionally, T MRCA and F ST decreased sharply to background levels outside of the inversions, demonstrating the potential for gene flow and recombination to homogenize variation in these regions. We interpret this as evidence that these inversions help maintain linkage disequilibrium among multiple divergently adaptive variants in regions susceptible to homogenization (Kirkpatrick and Barton 2006;Guerrero et al. 2012). The presence of these inversions in addition to divergence in regions of generally low recombination (Glazer et al. 2015), therefore, further supports the hypothesis that the recombinational landscape can influence where in the genome adaptive divergence can occur (Roesti et al. 2013;Samuk et al. 2017) and emphasizes the degree to which gene flow among divergently adapted stickleback populations has impacted global genomic diversity.

Conclusions
Selection operating on two very different timescales-the ecological and the geological-has shaped genomic patterns of SGV in the threespine stickleback. On ecological timescales, selection drives phenotypic divergence in decades or millennia by sorting SGV across geography and throughout the genome (Hendry et al. 2002;Hohenlohe et al. 2010;Lescak et al. 2015;Roesti et al. 2015). Our findings show that persistent ecological diversity and continual local adaptation of stickleback has set the stage for long-term divergent selection and for the accumulation and maintenance of adaptive variation over geological timescales. Some of the genetic variants fueling contemporary, rapid adaptation may even have been present--and under selection--since before the threespine-ninespine stickleback lineages split. The genomic architecture of ecological adaptation in one population is therefore the product of millions of years of evolution taking place across multiple populations, many of which are now extinct. These findings underscore the need to understand macroevolutionary his-tories of genetic variation when studying microevolutionary processes, and vice versa. AUTHOR CONTRIBUTIONS T.C.N. and W.A.C. conceived of the project and designed sampling, sequencing, and analysis. T.C.N. prepared sequencing libraries, wrote software, and performed data analysis. T.C.N. and W.A.C. wrote the article.
Associate Editor: J. Slate

Supporting Information
Additional Supporting Information may be found in the online version of this article at the publisher's website: Figure S1. Longer sequences reduce variance in estimates of sequence diversity and divergence. Figure S2. Accurate phasing of RAD loci even at low population-level sampling. Figure S3. RAD-seq effectively samples genome-wide sequence diversity. Figure S4. Relative (FST) and absolute (dXY) sequence divergence are positively correlated genome-wide in two instances of marine-freshwater divergence. Figure S5. TMRCA outlier regions remain outliers after removing highly differentiated RAD loci. Figure S6. Marine-freshwater divergence in threespine sticklebacks is associated with reduced signal of monophyly of threespine haplotypes. Table S1. Sequencing summary for threespine stickleback samples. Table S2. Genome assembly statistics for Pungitius pungitius.