Contemporary ancestor? Adaptive divergence from standing genetic variation in Pacific marine threespine stickleback

Populations that have repeatedly colonized novel environments are useful for studying the role of ecology in adaptive divergence – particularly if some individuals persist in the ancestral habitat. Such “contemporary ancestors” can be used to demonstrate the effects of selection by comparing phenotypic and genetic divergence between the derived population and their extant ancestors. However, evolution and demography in these “contemporary ancestors” can complicate inferences about the source (standing genetic variation, de novo mutation) and pace of adaptive divergence. Marine threespine stickleback (Gasterosteus aculeatus) have colonized freshwater environments along the Pacific coast of North America, but have also persisted in the marine environment. To what extent are marine stickleback good proxies of the ancestral condition? We sequenced > 5800 variant loci in over 250 marine stickleback from eight locations extending from Alaska to California, and phenotyped them for platedness and body shape. Pairwise FST varied from 0.02 to 0.18. Stickleback were divided into five genetic clusters, with a single cluster comprising stickleback from Washington to Alaska. Plate number, Eda, body shape, and candidate loci showed evidence of being under selection in the marine environment. Comparisons to a freshwater population demonstrated that candidate loci for freshwater adaptation varied depending on the choice of marine populations. Marine stickleback are structured into phenotypically and genetically distinct populations that have been evolving as freshwater stickleback evolved. This variation complicates their usefulness as proxies of the ancestors of freshwater populations. Lessons from stickleback may be applied to other “contemporary ancestor”-derived population studies.


Background
The ecological theory of adaptive divergence predicts that populations diverge phenotypically and genetically if they reside in distinct environments [1,2], potentially resulting in speciation. Some of the most striking examples of adaptive divergence come from species in which contemporary populations persist in environments putatively occupied by ancestral populations (e.g. [3][4][5][6][7][8][9][10][11]). To the extent that such populations have not undergone evolution, contemporary descendants of populations in ancestral environments can be used as a proxy for the ancestor, with phenotypic and genetic differences between these "contemporary ancestors" and "derived" populations being used to infer the direction, source, and pace of adaptation to derived environments. However, recent changes in the ancestral environment may stimulate evolutionary responses in the contemporary populations that inhabit it, complicating their utility as a proxy.
Standing genetic variation (SGV), defined as the variety of alleles segregating in a population [12,13], is expected to play an important role in parallel evolution. In particular, SGV permits rapid adaptation compared to de novo mutation, and increases the likelihood that the same beneficial allele will be present in different derived populations [14][15][16]. The role of SGV in adaptive divergence is readily measurable: if an allele fixed in the derived population is present in the contemporary ancestor at low frequencies, it likely contributed to adaptation [13]. However, the inference that an allele present in the contemporary ancestral population resulted in adaptation via SGV requires three assumptions. (1) The subset of individuals that originally colonized the derived environment must have contained the rare adaptive allele at some frequency; otherwise it arose from de novo mutation or subsequent gene flow. (2) The contemporary ancestor has undergone little evolution, including gene flow from the derived population (e.g., [17]). (3) The ancestral population has to have been properly characterized (i.e., population structure and allele frequencies associated with SGV). If there are multiple potential ancestral populations, each with a different pool of SGV, inference about the source and pace of evolution in the derived population will vary depending on which putative ancestral population is investigated (e.g. [12]). These assumptions must be verified to characterize accurately the role of SGV during population divergence.
Threespine stickleback (Gasterosteus aculeatus) provide perhaps the best documented examples of adaptation from SGV. Marine threespine stickleback occur widely in the northern hemisphere, including along the Pacific coast of North America from Alaska south to southcentral California. Across the north Pacific coast, much freshwater habitat formed recently (~10,000-20,000 years ago) in association with isostatic rebound following glacial retreat. The subsequent colonization of this habitat by stickleback allows tests of the significance of de novo mutation and SGV for adaptation (e.g. [17][18][19]). For instance, marine stickleback bodies are often covered by > 29 bony lateral plates, but fewer plates (0-10) have evolved in parallel in freshwater populations through selection on a rare marine allele [18]. Despite numerous studies indicating the role of SGV at either a single locus for platedness (Ectodysplasinhereafter Eda) or for multiple loci with unknown phenotypic effects [20][21][22], assumptions about the appropriateness of considering extant marine sticklebacks as representative of the ancestors of freshwater populations remains untested. Despite evident genetic variation in threespine stickleback among geographic clades [23][24][25][26], marine stickleback on the eastern Pacific are largely assumed to constitute a single population (e.g. [27][28][29][30][31]). This assumption is justified by the absence of barriers to gene flow in the marine environment [27], the migratory capacity of marine stickleback [32], the relative "evolutionary stasis" of marine stickleback inferred from the fossil record [28,29], and the low marine population structure reported from several local studies ( [20,33], but see [34]). Nevertheless, substantial evidence indicates local adaptation even in highly migratory marine fishes [35][36][37], and indeed in Baltic Sea threespine stickleback [38][39][40]. If Pacific threespine stickleback constitute a single population, their large population size should limit the effects of genetic drift and high gene flow should offset local adaptation [41]. Given these conditions, marine stickleback should not exhibit local differentiation that would generate regional differences in the initial colonists of freshwater lakes and streams. SGV should be the same along the Pacific coastand all freshwater populations could have evolved from the same initial pool of marine SGV facilitating parallel adaptation. These assumptions require formal testing to elucidate the role of SGV during adaptive divergence.
In this study, we consider phenotypic and genotypic variation of > 200 marine threespine stickleback from eight locations from California to Alaska to test hypotheses about the genetic structure of marine stickleback and its evolutionary consequences. Based on variation in plate phenotypes and genotypes associated with SGV at Eda, three-dimensional body morphology from micro-computed tomography (μCT) scans quantified using geometric morphometrics, and Genotype-by-Sequencing [42], we assess whether marine stickleback constitute a single population. By doing so, we test assumptions about the distribution of SGV in "contemporary ancestors" of freshwater stickleback. Additionally, we test predictions regarding the influence of SGV on the source of adaptation based on genomic sequences for a freshwater population from British Columbia. We specifically test the following null predictions: (1) Marine stickleback populations will not vary in the frequency and content (e.g. private alleles) of SGV. (2) Similarly, marine stickleback will not exhibit genetic population structure, which would otherwise influence the SGV regionally available for selection. (3) Marine populations will exhibit no phenotypic divergence in body shape or platedness [e.g. 18]. (4) If differences in SGV among marine populations occurs, genetic variation will show no evidence of having been shaped by natural selection. (5) If population structure in marine stickleback occurs, geographic proximity to a freshwater population will determine the extent of genetic divergence between marine and freshwater stickleback. (6) Differences in SGV among marine populations, if they occur, will not affect the candidate loci identified as contributing to adaptation in freshwater stickleback.

Methods
Threespine stickleback (n = 383, Table 1) were collected with minnow traps or seines during the summers of 2010 (Brannen Lake, British Columbia, hereafter BCFW), 2012 (Alaska) and 2013 (all other localities). Sampling locations extended along a 21.8 degree latitudinal spread (Table 1, Fig. 1), from California (south to north, CA01, CA02, CA03), through Oregon (OR01, OR02), the Puget Sound area of Washington (WA01), Vancouver Island (BC01) and Alaska (AK01). Locations varied in terms of benthos, freshwater input, and protectionfor instance, CA01 fish were sampled in a slough with freshwater input determined by precipitation, while OR02 were collected near a tidal gate close to the mouth of a river. Other marine species were collected alongside stickleback, such as bay pipefish (Syngnathus leptorhyncus) or smelt (Atherinops/ Atherinopsis sp.). Adults were captured in all localities with the exception of WA01, while OR01 contained a range of age classes. Stickleback were euthanized using buffered tricaine methanesulfonate (MS-222) or Eugenol (clove oil) and preserved in 70% ethanol. Fin clips were preserved in 95% ethanol for later sequencing. All collections were conducted in accordance with CCAC guidelines (AUP AC13-0040) and state/provincial/ national collection and import permits.
Sex was identified using primers developed by [43] that amplify sex-specific alleles at the idh locus. Alleles were visualized in a 2% agarose gel for 367 individuals.

Library preparation and analysis
Reduced representation DNA sequencing was used to generate Single Nucleotide Polymorphisms (SNPs) in order to assess population structure and adaptive divergence. Two hundred nanograms total genomic DNA was extracted per fish in January 2016 using Qiagen  [44]). Thirty to thirty-five individuals were included from each marine location, and 15 from BCFW. After digestion-ligation, fish were pooled into groups of nine. Cleanup and size selection were performed simultaneously using SPRI beads (Beckman Coulter), at a bead ratio of 0.8× and 0.61× for left and right-side cleanup, respectively. This left a fragment range between 250 and 600 bp. Pooled samples were divided into three technical replicates to ameliorate stochastic differences during PCR, and were amplified.
Replicates were pooled and left-side cleaned using SPRI beads. Pooled samples were quantified using a 2200 TapeStation (Agilent) and Qubit (Thermofisher) dsDNA high sensitivity assay. Equal volumes of each 2 nM pooled sample were pooled to make the final library. Library preparation followed the Illumina protocols for the Illumina NextSeq 500 Mid-Output kits with version 2 chemistry. Two sequencing runs were completed on the Illumina Next Seq 500 using 150 cycles and different final library concentrations -the first at 1.8 pM final concentration, the second at 1.1 pM final concentration. A 20% PhiX spike-in was used for both to compensate for the low diversity nature of the library. Results from the two sequencing runs were merged. Sequenced reads were cleaned and processed using Stacks v.1.35 [45,46]. Reads were de-multiplexed and cleaned using process_radtags, rescuing barcodes if the correction of a single sequencing error made them identifiable. GSnap [47] was used to align reads to the stickleback reference genome (Ensembl release 72 [48]), allowing for five mismatches with soft masking disabled. SNP calls were corrected in rxstacks using a bounded SNP model with an upper error rate of 0.1. Stringent filtering criteria were applied to the data, with slightly different filtering criteria used to address different questions. In order to determine population structure, each sampling site was assumed to constitute a distinct population. The filtering criteria for this "marine site" data set included: i) log likelihood threshold > − 60; ii) sequenced in more than 75% of individuals; iii) in 6 of 8 populations; iv) with minimum 4× coverage; v) minimum minor allele frequency of 2%, and vi) F IS > − 0.3. Individuals were included if they retained > 10,000 RADloci after cleaning (n = 239, Table 1). After population structure was determined, Adegent-recognized clusters rather than sampling location were used in the second Stacks run. This "marine cluster" data set used the same filtering criteria as above, with the exception that the variant needed to be present in all clusters. Finally, a "marine-freshwater" data set was run, which treated each sampling site as a distinct population and additionally included a freshwater population (see below). In this case the variant needed to be present in all eight marine populations and in the freshwater population. For each data set, all population statistics except for F ST were calculated using the populations module of Stacks.

Population genetic structure
Pairwise global and per-locus F ST were calculated using the Weir and Cockerham [49] adaptation implemented in hierfstat v0.04-22 [50] in R [51]; pairwise global F ST were tested for significance (> 0) using 999 permutations in GenoDive [52]. Discriminant Analysis of Principal Components (DAPC) [53] was used on the "marine site" data to assess population structure in the marine environment using Adegenet v2.0.1 [54], as it has shown to perform better than Structure under a stepping-stone model of dispersal [53]. The optimal number of Principal Components (PCs) to retain was calculated using both xvalDAPC and a-scores, which gave similar answers. The optimal number of clusters was assigned based on the lowest Bayesian Information Criteria (BIC) score using k-means clustering. As several possible k clusters had similarly low BIC scores, analyses were run and compared using 3 to 8 clusters.
An analysis of molecular variance (AMOVA), implemented in poppr v2.3.0 [55,56] using the Ade4 package [57], was used to determine the proportions of genetic variance among versus within sampling sites or Adegenetrecognized clusters. Missing values were replaced with the average frequency for a locus; ignoring missing values did not alter overall patterns. To explore the possibility of cryptic population structure, each sampling site was further analysed individually using Stacks and Adegenet.
The distance between each sampling site was measured as distance along the coast (km) using Google Maps. Distances were measured to or from the mouth of each bay. Neighbouring localities were separated by 242-479 km, except for BC01-AK01, which were separated by approximately 2500 km of coastline. The location of WA01 in Puget Sound resulted in all locations south of Washington being closer to BC01 than they were to WA01. Genetic distance was calculated using the pairwise global Weir and Cockerham F ST measures from hierfstat. Geographic and genetic distance matrices were compared using a Mantel test from the Adegenet package with 999 replications to determine Isolation-by-Distance (IBD).
Population statistics were also estimated for the "marine cluster" data set in Stacks using the optimal Adegenet-identified clusters.
A phylogenetic network was calculated using SNPhylo [58] and visualized using FigTree v.1.4.3 [59]. The "marine sites" data set was used, but SNPhylo additionally filtered loci based on linkage disequilibrium. Individuals were colour-coded according to their recognized genetic cluster. A hierarchical clustering tree was additionally constructed using BayPass v.2.1 [60].

Platedness
Plate variation among populations was assessed using plate number and Eda genotype. Adult stickleback (i.e. fish > 30 mm standard length (SL)) (n = 281, Table 1, Additional file 1: Table S2) were stained in Alizarin red. Plate number, including keel, was counted on both sides of the body and summed. Low-plated individuals without keel (LPNK) were defined as individuals with < 20 anterior plates. Partially-plated keeled (PPK) fish had 21-59 plates, including at least one plate at the caudal keel. Fully-plated keeled (FPK) fish had ≥ 60 plates. Additionally, some low-plated fish had a keel (LPK) and were defined as having < 20 anterior plates plus additional plates at the caudal keel. Partially plated stickleback that lacked a keel (PPNK) had > 20 anterior plates but had no plates at the caudal keel. Individuals were also genotyped at the Stn382 locus [18] as this microsatellite is linked to an indel in intron 1 of the Eda gene, yielding a 218 bp "fully-plated" allele (C) or a 158 bp "low-plated" allele (L) [61]. Genotyping followed the protocol of [43]. Individuals were genotyped as LL (homozygous for the low-plated allele), CL (heterozygous), or CC (homozygous for the fully-plated allele). This approach allowed juveniles (< 30 mm SL) to be included in the analysis (total n = 361), and provided genetic information at a locus with known adaptive significance that was not recovered from sequencing. Hardy-Weinberg equilibrium (HWE) was assessed for Stn382 for each marine site using a goodness-of-fit Chi-squared test.

Morphometrics
Phenotypic variation among populations was further assessed using morphometric analysis. Stickleback > 30 mm SL that had been preserved with relatively little bending (n = 272) were straightened, and spines and fins held flat against the body using plastic wrap. μCT scanning at a resolution of 20 μm was conducted in a standardized fashion for all individuals using a Scanco μCT35 instrument (Scanco AG). Three-dimensional images were generated from the anterior point of the premaxilla to the posterior tip of the pelvic spine, using standardized isosurface thresholds in Amira 5.4 (FEI Visualization Sciences Group). Fifty-five landmarks were plotted on the left side of each fish (Additional file 1: Table S1, Fig. 2) and raw landmark scores were exported to MorphoJ v1.06a [62] for further analyses. A prior study had removed the operculum on the left side of all AK01 stickleback, so landmarks were plotted on their right sides. Data were first transformed to remove differences associated with isometric scaling, rotation and translation using Procrustes superimposition. Residuals from a within-marine site multivariate regression on centroid size were estimated and used in all subsequent analyses. Principal Components Analysis (PCA) determined the major axes of phenotypic variation. Canonical Variate Analysis (CVA) was used to determine Procrustes distances using sex and marine site as categorical variables, although BC01 had only one female and OR01 included one individual of unknown sex. The significance of Procrustes distances among pairwise comparisons of marine site-sex combinations was determined based on 10,000 permutations and a corrected α of 0.0005. Discriminant Function Analysis (DFA) was used to determine the likelihood that individuals could be reassigned to their site of origin, given their phenotypes. For this analysis the effect of sex on reassignment success was not assessed.

Selection on phenotypic variation
Patterns of phenotypic variation among populations may be the result of genetic drift, natural selection, or phenotypic plasticity. One way to rule out neutral evolutionary processes is to compare estimates of phenotypic divergence (P ST ) with neutral expectations based in part on observed neutral genetic differentiation. Observed phenotypic divergence was estimated as: where σ 2 B and σ 2 W were the between-and withinpopulation components of variance, respectively, for plate count and the first four PCs from the morphometric analysis (as per [63,64]). Variance components were estimated for all marine sampling sites together (global P ST ) and pairwise using lme4 [65], with population as a random effect. Genetic divergence at the Stn382 locus for Eda (F STQ ) was estimated using the Weir and Cockerham method in Genepop V4 [66]. Neutral genetic divergence (F ST ) was estimated in hierfstat using non-genic SNPs identified from our data set using Biomart [67]. Non-genic SNPs may still be linked to loci under selection, so this approach provides a conservative estimate of neutrality.
Selection was inferred based on two methods. The first assessed the association between P ST -F ST and F STQ -F ST using Mantel tests. This measure is based on the expectation that phenotypic or QTL divergence will be uncorrelated with neutral genetic divergenceby extension implicating selection to explain such patterns. The second test involved Whitlock and Guillaume's [68] method using the R-code from [69]. In brief, the expected between-population variance component for a neutral phenotype was estimated by using observed non-genic F ST and observed within-population variance component for the phenotype: As per [69], the distribution of neutral σ 2 B was estimated by generating a χ 2 distribution with six degrees of freedom (one less the number of sampling sites excluding Washington), and multiplying a randomly drawn value from this distribution by σ 2 B . From this new distribution expected neutral σ 2 B were drawn 10,000 times and used to create a distribution of neutral P ST -F ST . The observed P ST -F ST was then compared to this distribution and the quantile of the neutral distribution that lay beyond the observed value was used as the probability of the observed outcome in the absence of selection, p. Under the expectation of no selection, p is > 0; selection is evidenced if p = 0. F STQ -F ST values was also compared to the neutral P ST -F ST , as per [40].
To quantify the relation between phenotype and Stn382 genotype, a generalized linear model (GLM) was fit to the data in R, using the glm routine, as per [40], with plate number as the dependent variable, genotype as a fixed effect, and using a log-link function with a quasi-Poisson error distribution. Furthermore, a Mantel test was used to estimate the correlation between pairwise F STQ and P ST measures.

Selection on genetic variation in the ocean
Under the assumption that marine stickleback populations have a shared history, the covariance matrix of population allele frequencies (Ω) was estimated in BayPass [60] using the "marine site" data. From this a hierarchical clustering tree [60] was generated, assuming no gene flow. A covariate-free genome scan was then performed to identify outlier loci putatively under selection, using per-locus measures of differentiation (XtX). The simulate.baypass function was used to estimate the posterior predictive distribution of XtX using a pseudo-observed data set (POD) [60]. Any loci in the "marine site" data set with XtX values above the POD-estimated threshold were scored as outlier loci potentially under selection. Genic outliers were identified using BioMart [67].

Marine-freshwater genetic divergence
To assess the extent to which the choice of putative "contemporary ancestor" affected inference about adaptation to fresh water, pairwise F ST was estimated for each marine-freshwater pair using hierfstat. In BayPass, covariance matrices, PODs and XtX thresholds were estimated for each marine-freshwater comparison. Outlier loci were examined to determine if the same outliers were being consistently recovered irrespective of the origin of the marine fish.

Sequencing results
Over 192 million reads passed initial filters (Additional file 1: Table S3) in the "marine sites" data set. Two hundred eighty-two RAD-loci were excluded due to excess heterozygosity. Filtering minor alleles at a threshold of 2% reduced the number of retained loci by 32%. After filtering, between 230,010 (AK01) and 426,018 (OR01)  Table S1 for identity of landmarks loci were retained for each site sampled, generating between 1877 (AK01) and 5204 (OR01) SNPs (Table 2), for a total of 6655 variant loci.

Standing genetic variation
All marine samples exhibited SGV, ranging from an average of 0.82% (AK01) to 1.22% (OR01) of the total SNPs genotyped in a given population; however, the pool of SGV varied from California to Alaska (Table 2). Stickleback from each marine location contained multiple private alleles (alleles found only at that location) (Fig. 3) and were polymorphic for a portion of the variant loci (loci that were polymorphic in at least one marine site). Polymorphism among variant loci varied from 57% (AK01) to 80% (OR01). For variant loci, the average frequency of the major allele (present in > 50% of all sequenced stickleback) ranged from 88% (OR02) to 93% (AK01) (Additional file 1: Figure S1), suggesting that the frequencies of SGV also differed among locations. Heterozygosity ranged from 0.10 (AK01) to 0.17 (OR02) for variant loci. Population-level average F IS over all variant loci varied from 0.032 (CA01) to 0.064 (OR01) ( Table 2). Although F IS was close to 0 for most loci, it approached 1 for a few loci (Additional file 1: Figure S2).
Significant population genetic structure was detected. The best supported number of clusters from the eight marine locations sampled was five (BIC = 1379, Fig. 4, Additional file 1: Table S4). The five clusters were, from south to north, CA01, CA02, CA03-OR01, OR02, and WA01-BC01-AK01. The CA03-OR01 cluster also contained seven individuals from OR02 and a single individual from AK01; otherwise individuals clustered with others from their sampling locality. The possibility of a single genetic cluster was as well-supported as ten genetic clusters (BIC = 1392). Three to six clusters had BIC values that differed little from the best-supported model. Altering the number of putative clusters revealed different population structures (Additional file 1: Table S4), with WA01 and AK01 continuing to cluster until k = 13. Cryptic population structure was not evidenced for any sampling locality (Additional file 1:   Table S5). Even when WA01, BC01, and AK01 were included in a single analysis, k = 1 was the best supported cluster (BIC = 436.9). However, k = 2 (BIC = 438.2) and k = 3 (BIC = 440.2) still separated individuals by locality. The genetic variance partitioned between clusters by AMOVA was low (19%) compared to within clusters (81%). When clusters were not considered, 33% of variation occurred between marine sites. A Mantel test of pairwise geographic distances and pairwise Weir and Cockerham F ST was non-significant when all sites were included (r = − 0.2, p = 0.8). This was largely driven by the extreme distance between the WA01, BC01, and AK01 sites. If AK01 was excluded from the analysis the association between geographic and genetic distance was weakly significant (r = 0.5, p = 0.02) (Additional file 1: Figure S3).
The phylogenetic network largely agreed with Adegenet assignment (Fig. 5). The network revealed greater intermixing of groups than did Adegenet, but WA01, BC01, and AK01 still largely clustered together and comprised a separate lineage from most southern  Table 1 for label meanings  Table S2).

Morphometrics
Phenotypes varied extensively among sites. The first eight Principal Components explained 71% of all phenotypic variance (Additional file 1: Figure S5). The first two Canonical Variates (CVs) explained 57% of the variation among combinations of site and sex, and the first four CVs explained 77% (Additional file 1: Table S8). CV1, after accounting for differences in centroid size, revealed that BC01 fish had narrow, streamlined bodies with dorsal and ventral landmarks both shifting inward relative to the consensus fish (Figs. 7 and 8). Californian fish were grouped close together on CV1 and had squatter, less streamlined bodies with dorsal and ventral landmarks shifted away from one another relative to the consensus. AK01, OR01 and OR02 had intermediate phenotypes between BC01 and California. CV2 showed a gradual transition from CA01 to BC01, but AK01 was clearly distinct from all other sites along this axis. AK01 showed substantial dorsolateral and anterior-posterior constriction of the body relative to all other sites (Figs. 7 and 8).
The sexes differed morphologically at all sites except AK01 (Additional file 1: Table S9), but the sexes still largely grouped together according to sampling location. Fig. 4 Adegenet-identified clusters for k = 5. Inset shows hypothetical range of each cluster. Note that the cluster identified as CA03, OR01 contains one AK01 and seven OR02 individuals. See Table 1 for label meanings CA02 and OR02 were exceptions, with males from both sites clustering with OR01 males. Similarly, CA01 and CA02 females had morphologies that were not significantly distinct.

P ST -F ST and F STQ -F ST comparisons
P ST was estimated as 0.46 for platedness, 0.60 for PC1, 0.29 for PC2, 0.23 for PC3, and 0.09 for PC4. Fig. 5 a Result of the phylogenetic analysis using SNPhylo. Individuals are colour-coded according to their five Adegenet-recognized clusters (Fig. 4). b Hierarchical clustering tree using BayPass. See Table 1 for label meanings Fig. 6 The frequency of different Eda genotypes using the Stn382 marker, for each sampling site. "The North" refers to samples from Washington, British Columbia, and Alaska. CC = homozygous for the fully-plated allele. LL = homozygous for the low-plated allele. CL = heterozygous F STQ was 0.60. Plate P ST and F STQ greatly exceeded the range of the neutral P ST -F ST distribution (p = 0 for both), as did P ST for PC1 (p = 0). P ST for PC2 was marginally significant but within the tail of the neutral distribution (p = 0.002), while P ST for PC3 (p = 0.02) and PC4 (p = 0.7) were well within the neutral distributions (Fig. 9).
Of all of the Mantel tests between P ST , F ST , F STQ , and distance, only one was significant and two were marginally significant with corrected α = 0.05/18 = 0.0027: PC1 P ST correlated positively with geographic distance, plate count P ST correlated positively with F STQ , and, surprisingly, PC1 P ST correlated positively with F STQ (Table 4, Additional file 1: Figure S6 and S7). The relation between plate P ST and F STQ was substantiated with a generalized linear model which showed a decrease in plate number with number of L alleles (null deviance = 4813 with 277 d.f., residual deviance = 349 with 275 d.f., p < 0.0001, Hosmer and Lemeshow goodness of fit text:

Selection on genetic variation in the ocean
102 of 6655 loci from the "marine sites" data set were flagged as outliers, using an XtX threshold of 15.52 (Additional file 1: Figure S8). Although variant loci were sequenced across all 21 chromosomes and an additional 65 scaffolds, outliers were only detected on 16 chromosomes and two scaffolds. Of these, 8 of 36 (22%) variant loci on scaffold 37 were outliers, followed by 13 of 178 (7%) on linkage group (LG) XXI, 11 of 184 (6%) on LGXIX, and 18 of 498 (4%) on LGIV. Of these 102 outliers, 20 were located within 16 genes (Table 5), although none of these genes have been previously studied in stickleback.

Genetic differentiation and outlier analysis for marinefreshwater comparisons
The "marine-freshwater" data set identified 132,415 loci, of which 1912 were variant. As expected, marine-freshwater divergence was high and in all but one instance was > 0.25 ("very great" differentiation).  Table S11). The maximum difference in pairwise per-locus F ST estimates for a single locus was 0.991, and the average difference was 0.19. 314 (17%) SNPs had a minimum F ST estimate of little genetic differentiation (< 0.05) in at least one marine-freshwater contrast, but great genetic differentiation (> 0.25) in another.
97 of 1912 loci were identified as outliers in at least one marine-freshwater comparison, using XtX thresholds of 4.9 (CA03 -FW) to 6.0 (CA01 -FW). Seventeen outlier loci resided in 13 genes (Table 6). Only 3 of 97 loci were flagged as outliers in all eight marine-freshwater comparisons, while nearly a third (31%) were flagged in only a single comparison. Five of the genic loci were outliers in seven comparisons (tub, S100P, and three novel genes), and five were outliers in only one (including a different S100P locus).
Surprisingly, 12 of the loci flagged as outliers in the "marine sites" data set were also flagged as outliers in at least one marine-freshwater comparison (Table 6), including a novel gene on LGII (n = 7 comparisons), and csnk1g2a (n = 3 comparisons). The gene OVGP1 was also flagged as containing outlier loci in both analyses, but different loci were flagged in each case.  Table 1 for label meanings

Marine stickleback exhibit between-population genetic variation
The significance of SGV for parallel evolution depends on its occurrence in the ancestral populations. Marine stickleback harbour SGV, and it differs between populations. Differences in gene expression between two marine BC populations suggested this possibility [10], but here it has been quantified across an extensive latitudinal range. The extent of SGV, 0.8-1.2% of all sequenced loci, is intermediate to that reported from other studies [34,70]. Nucleotide diversity varied from 0.0016 to 0.0027, consistent with results from Alaska (0.0022 and 0.0025 [20]) and slightly lower than that reported from Oregon (0.003-0.0036 [34]). All marine locations harboured some degree of private alleles, even after ignoring minor alleles at < 2% frequency, suggesting that not just frequencies of SGV but also content of SGV can differ from site to site. Furthermore, the superiorly. a CV1 for a BC01-type body shape; b CV1 for a CA01-type body shape; c CV2 for an AK01-type body shape; d CV2 for a BC01-type body shape. Light blue wireframe shows the consensus morphology, while dark blue shows the conformational change best-known example of SGV, Eda, was present at varying frequencies between populations and is likely under selection in the marine environment. Such variation in the content and frequency of SGV, in turn, led to compelling evidence for population genetic structure. Marine threespine stickleback showed substantial population genetic structure along the Pacific coast of North America. Although F ST values (average F ST = 0.088) were generally lower than those reported for marine-freshwater divergence (e.g. [20,34,71]), they were higher than those reported for other marine stickleback populations along the North American Pacific coast (two Alaskan populations: F ST = 0.0076 [20]; three Oregonian populations: F ST = 0.007 [34]). However, they align with studies from Europe [39,40,72].
Five genetic clusters were identified for the eight sampled localities, although structuring was hierarchical, with some clusters more genetically diverged than others. The most widespread cluster occupied > 2700 km of coastline from Washington to Alaska. Marine stickleback from this genetic cluster have been well characterized, with genetic divergence reported to be low between populations separated by up to 1000 km [20,33,73]. Such low structuring between proximate marine populations has led to a basic assumption in stickleback literature that marine stickleback exhibit little genetic or phenotypic diversity globally (e.g. [27-29, 64, 74])an assumption supported by the iconic image of distinct freshwater stickleback forms radiating from a single fully-plated marine stickleback type (e.g. [75]). In contrast, results from a broad range suggest that such generalizations should be restricted to the northern genetic clusterand even it contains morphological and genetic differentiation that could be adaptively significant.
The southern genetic clusters were sequentially separated by a few hundred kilometres, well within the migratory ability of marine stickleback [76][77][78]. IBD was evident only after removing AK01 from the dataset, suggesting that limited migration could explain patterns of divergence between the southern genetic clusters. However, IBD needs to be interpreted with caution, as geographic distance was correlated with latitude, and latitudinal variation can be associated with environmental clines [79]. Whatever the causes that shape genetic variation between stickleback populations, the distribution of SGV among different marine populations affects inference about the source and pace of selection in the freshwater environment, and complicates attempts to uncover loci that are under selection in derived populations.

Marine stickleback exhibit between-population phenotypic variation
Marine stickleback are generally considered to be fullyplated (e.g. [27]), yet the Eda genotype for low-platedness has an ancient marine origin [18]. The low-plated allele has been hypothesized to exist in the marine environment as SGV only when transported from the freshwater environment [17] or when masked by marine modifying alleles [18]. If the low-plated allele exists at low frequencies in the ocean, behaviours that facilitate the movement of low-plated marine stickleback into fresh water could also account for the consistent colonization of rare low-plated stickleback in lakes and streams [80]. We found substantial variation in the frequency of the low-plated allele, to the point that it was the major allele in some Californian and Oregonian populations, but was absent from BC01. This finding is consistent with other records of low-plated marine stickleback in California ( [81], but see [18,82]) and high frequencies of low-platedness in European marine stickleback [61,74,[83][84][85]. However, it contrasts with other Pacific and Atlantic North American studies that focussed on northern sites [86][87][88][89][90][91].
The observed SGV at Eda could affect the rate at which adaptation to lakes occurred in the past. Indeed, it may explain why reduced plate size has evolved in some freshwater populations, rather than reduced plate number, as an alternative strategy that may have been required in the absence of SGV at Eda [92,93]. Such variation at Eda is a particularly striking reminder that the function of full-platedness in marine stickleback remains unknown (for an example of its possible use, demonstrated in freshwater populations, see [94]).
Marine stickleback body shape also varied extensively along the coast. Californian populations tended to have squatter body shapes that appeared to be less streamlined than their northern counterparts. The functional Table 4 The observed correlation and p-values (p) for Mantel tests between geographic distance, neutral genetic distance (F ST ), phenotypic distance (P STfor plates or the first four Principal Components (PCs) of body shape), or genetic distance at Eda (F STQ )

Comparison
Observed  LG Linkage Group. XtX provides an average if the number of SNPs is > 1. The XtX threshold, above which a SNP was considered an outlier, was 15.52  Table 6 Pairwise XtX values and outlier locus information for each of the eight marine-freshwater comparisons.
LG is linkage group, although scaffolds are also included. Position refers to the nucleotide position along the linkage group (Continued) LG significance of these differences requires testingbut it is interesting that the streamlined fish were from a single genetic cluster, while the squat Californian fish exhibited significant population structuring between neighbouring localities. This potentially indicates extensive migration along the northern coast that is not mirrored in the south. Jamniczky et al. [95] reported considerable morphological divergence between neighbouring sampling sites in British Columbiagroups presumably with little to no genetic divergence, implicating plasticity as a driver of morphological variation. Morris et al. [79] similarly reported variation in vertebral number and standard length with latitude. However, two related analyses suggest that selection may also play a role in shaping phenotypic diversity. Based on DeFaveri and Merilä's [40] method, Pacific coast stickleback exhibited selection for platedness similar to that observed for Baltic Sea stickleback. There was also suggestive evidence for selection on PC1 of body shape, which largely corresponded to CV1more streamlined northern fish, more squat southern fish. The role of plasticity in affecting these results remains to be determined. PC1 was marginally associated with Eda genotype. To our knowledge, this is the first study to demonstrate pleiotropic or linked effects of Eda on body shape in the marine environment (for marine-freshwater or freshwater-only evidence for pleiotropy, see [96]; for other forms of phenotypes associated with Eda see [80,[97][98][99][100][101]), which could result in low-platedness being an indirect target of selection in some marine habitats. It is possible, for instance, that a pleiotropic relationship exists between body shape, Eda, and thermal tolerancea possibility suggested by the relationship between low plate frequency and latitude in anadromous populations of Europe [83]. This is one of several possible explanations that requires formal testing.
Although Eda is the best-characterized example of SGV, the outlier analysis revealed other potential candidate genes for selection in the marine environment. Extensive differentiation was evidenced at some loci. Local adaptation despite gene flow has been found in other marine fish populations [35,[102][103][104], including European marine stickleback [72]; but the extent to which stickleback south of Washington exhibited population structure was unanticipated.
Smaller bodies and distinct body shapes tend to evolve in freshwater stickleback populations [105,106], often with significant correlations between morphology and freshwater biotic and abiotic factors [107]. Given the morphological variation among marine stickleback, elucidating whether freshwater morphology is the result of plasticity, SGV, or de novo mutation requires informed decisions about what constitutes the marine ancestor.

Inferring the source and pace of adaptation in freshwater stickleback
The choice of marine stickleback affected inference concerning the source and pace of adaptation in one freshwater population from Vancouver Island. Few outlier loci were consistently recovered in all marine-freshwater CA01 CA02 etc. identifies the marine group that was compared to the Brannen Lake British Columbia freshwater population. The number in brackets () after each population identifier is the XtX threshold value used for that marine-freshwater comparison above which a locus was flagged as an outlier. If the locus was an outlier in that comparison its XtX value is given; otherwise the cell is blank. Total refers to the number of comparisons for which that locus was flagged as an outlier. If Marine = YES that same locus was also flagged as an outlier in the marine environment comparisons; had only a single marine population been used in this study, the outliers reported would differ depending on which marine population had been chosen. Many studies involve comparisons between geographically proximate marine and freshwater stickleback pairs (e.g. [21,22]), presumably to account for the possibility of population structure in the marine environment. Yet F ST was the lowest when the freshwater population was paired with a geographically distant population from northern Oregon, a finding that is difficult to reconcile with the assumption that the nearest marine population is the most suitable ancestral type. The occurrence of Japanese mtDNA haplotypes in Haida Gwaii lake populations that are not present in Haida Gwaii marine populations [108,109] suggests that this may not be an isolated incident. Similarly, several studies have noted, but not explained, the fact that northern marine stickleback are genetically more similar to southern than northern freshwater stickleback [17,21,22]. Nearly one third of outlier loci were only outliers for a single marine-freshwater comparison. Thus, one's choice of marine population could produce spurious inferences about the loci under selection, or miss true candidate genes. Clearly more information is needed going forward about the relationship between marine and freshwater stickleback.
Inferring the role of SGV during adaptive divergence in ancestral-derived comparisons requires three conditions to be met, for which we have provided varying evidence. The first, that colonists had to contain the same variants as present in the ancestral population, could not be tested here. Second, the contemporary ancestor has to have been relatively evolutionarily static. Both phenotypic and genetic evidence demonstrates that this condition does not hold for marine stickleback. Contemporary marine populations have diverged phenotypically in a manner beyond neutral evolutionary expectations, and genetically in a way partially explicable by selection. This means that contemporary marine threespine stickleback populations are genetically and phenotypically distinct from their own ancestorsand it was these ancestors that also originally colonized lakes and streams along the coast. Thus the term "contemporary ancestor" is a misnomer, as contemporary marine threespine stickleback populations do not reflect the ancestral condition. Interpretations of freshwater stickleback evolution need to be tempered by marine stickleback evolutionary history.
Third, the ancestral population has to have been properly characterized. Eastern Pacific stickleback exhibit some population genetic structure, although consistent with other reports [20,33] stickleback north of Oregon constitute a single genetic population. This means that along much of the coast freshwater environments were likely colonized by distinct marine stickleback populations, which differed in SGV frequency and content. Furthermore, it is likely that marine stickleback have exhibited range contractions and expansions along the southern and northern coasts throughout their evolutionary history, most recently in the north after the last glacial retreat [110]. This means that there is no a priori reason to expect that a marine population currently proximate to freshwater populations are descendants of the ancestors of those freshwater populations. Population structuring and evolutionary history thus changes our understanding of the ancestral condition of marine stickleback, and requires that we carefully consider the use of contemporary marine populations when addressing evolutionary questions.

Conclusions
Studies that compare marine and freshwater stickleback may need to adjust their methodologies in light of marine stickleback variation. The typical image in textbooks is of a single marine stickleback form from which numerous freshwater forms radiate [75]. Our data suggests that there is phenotypic and genetic variation in marine threespine stickleback which likely impacts freshwater stickleback diversificationto assume a single marine form is no longer tenable. Furthermore, assuming ancestral status for the marine population most geographically proximate to the freshwater population of interest is problematic, unless it can be directly demonstrated (e.g. [111]). So where does this leave the comparative method? One possibility would be to reconstruct the genotype of the ancestor to all eastern Pacific marine sticklebackbut this would ignore the important role that local variation has played in freshwater stickleback evolution. Another possibility could be to conduct larger-scale geographic sampling than has heretofore been done, of both marine and freshwater forms, in order to determine a more thorough evolutionary history of this species. Then, having taken evolutionary relationships into account, comparisons can be made using better-justified "contemporary ancestors".
"Contemporary ancestors" are used in a number of systems [3][4][5][6][7][8][9][10][11] for addressing evolutionary questions. They are particularly useful for determining the role of SGV during evolution, and for identifying the alleles involved in adaptation to new environments. Clearly, care must be exercised in characterizing these proxies of the ancestral form, as unaccounted population structure and current evolution can lead to spurious interpretations of adaptation. Whether the lessons from stickleback apply to other species with smaller geographic distributions or more limited opportunities for gene flow waits to be seen.

Additional file
Additional file 1: Additional tables and figures. This document contains information supporting the main text, including: the list of landmarks used for 3D morphometrics, population-specific details on sex and platedness, additional information pertaining to methodology, the distributions of major allele frequencies and FIS per population, isolation-by-distance analyses, details regarding Adegenet-recognized clusters, CVA and DFA morphometrics results, additional Mantel tests, global pairwise FST values for each marinefreshwater comparison, and results pertaining to the outlier analysis.