Genomic signatures of environmental selection despite near‐panmixia in summer flounder

Abstract Rapid environmental change is altering the selective pressures experienced by marine species. While adaptation to local environmental conditions depends on a balance between dispersal and natural selection across the seascape, the spatial scale of adaptation and the relative importance of mechanisms maintaining adaptation in the ocean are not well understood. Here, using population assignment tests, Approximate Bayesian Computation (ABC), and genome scans with double‐digest restriction‐site associated DNA sequencing data, we evaluated population structure and locus–environment associations in a commercially important species, summer flounder (Paralichthys dentatus), along the U.S. east coast. Based on 1,137 single nucleotide polymorphisms across 232 individuals spanning nearly 1,900 km, we found no indication of population structure across Cape Hatteras, North Carolina (FST = 0.0014) or of isolation by distance along the coast using individual relatedness. ABC estimated the probability of dispersal across the biogeographic break at Cape Hatteras to be high (95% credible interval: 7%–50% migration). However, we found 15 loci whose allele frequencies were associated with at least one of four environmental variables. Of those, 11 were correlated with bottom temperature. For summer flounder, our results suggest continued fisheries management as a single population and identify likely response mechanisms to climate change. Broadly speaking, our findings suggest that spatial balancing selection can manifest in adaptive divergence on regional scales in marine fish despite high dispersal, and that these conditions likely result in the widespread distribution of adaptive alleles and a high potential for future genetic adaptation in response to changing environmental conditions. In the context of a rapidly changing world, a landscape genomics perspective offers a useful approach for understanding the causes and consequences of genetic differentiation.

The means by which species respond to changing environments are determined by the potential for dispersal, the strength of selection, and the scale of environmental heterogeneity (Levins, 1968). Thus, understanding marine species' responses to environmental change requires knowledge of how phenotypic and genotypic spatial variations are shaped by the dispersal-selection balance across the seascape.
Our understanding of the spatial scale of genetic mixing in marine species has been undergoing a paradigm shift. The majority of marine species have historically been viewed as highly connected, well-mixed populations by virtue of their high dispersal rates, lengthy pelagic larval durations, and the assumption that larvae act as relatively passive particles (Siegel, Kinlan, Gaylord, & Gaines, 2003). Despite the positive correlation between pelagic duration and dispersal distance (Shanks, 2009), evidence has emerged that many wide-ranging marine species exhibit fine-scale population substructure at neutral loci (Benestan et al., 2015;Therkildsen et al., 2013) or are locally adapted to their surrounding conditions at loci experiencing natural selection (Conover, Clarke, Munch, & Wagner, 2006;Sanford & Kelly, 2011;Therkildsen et al., 2013). These lines of evidence suggest that marine populations are not as homogenous as previously thought (Hedgecock, 1986). For instance, species occupying a wide geographic area characterized by an environmental cline can become genetically matched to their local environment if dispersal is low enough or, alternatively, if selection is strong enough. Thus, spatial local adaptation (Kawecki & Ebert, 2004) or spatial balancing selection (Hedrick, Ginevan, & Ewing, 1976;Sotka, 2012;Whitlock, 2015), which is sometimes called spatially varying selection (Bernatchez, 2016), can result in genetic divergence at particular loci between populations occupying different habitats. This is important because, populations harboring adaptive polymorphisms across a heterogeneous landscape can serve as repositories of genetic variation, contributing potentially beneficial alleles to other populations along the cline, and enabling species persistence as environmental conditions shift.
Summer flounder (Paralichthys dentatus) support an economically important commercial and recreational fishery and inhabit estuarine and continental shelf waters characterized by environmental differences from Nova Scotia, Canada to Florida, USA. In particular, the warm Gulf Stream flows northward from Florida, hugging the coast until it spirals offshore at Cape Hatteras, North Carolina, creating a steep thermal gradient. Summer flounder exhibit a seasonal inshore-offshore migration pattern. Adults and juveniles reside in shallow coastal and estuarine waters during the summer months before moving offshore in the fall and early winter, when spawning occurs. Both adult and juvenile fish remain on the continental shelf throughout the winter and then return to coastal and estuarine waters in the late spring to early summer (O'Brien, Burnett, & Mayo, 1993;Packer et al., 1999). Overfishing of summer flounder resulted in a sharp population decline in the late 1980s and early 1990s, prompting severe fishing restrictions and conservation measures to improve stock abundance (Terceiro, 2011). Today, the summer flounder stock has been successfully rebuilt and is managed as a single population with state allocations of fisheries quota based on the 1980s population distribution. The center of summer flounder abundance lies offshore of the Mid-Atlantic states, between Cape Cod, Massachusetts and Cape Hatteras, North Carolina (Packer et al., 1999). Recent studies have suggested that summer flounder biomass has shifted within the Mid-Atlantic, due either to changing climate (Nye, Link, Hare, & Overholtz, 2009;Pinsky & Fogarty, 2012) or as a response to fishing pressure (Bell, Richardson, Hare, Lynch, & Fratantoni, 2014). Regardless of mechanism, shifts in summer flounder biomass have management consequences because of a growing mismatch between state allocations in quota and the geographic distribution of summer flounder biomass. Even though some fisheries can track shifting species, the disparity between the locations of fish and fishermen is likely to have socioeconomic consequences (Pinsky & Fogarty, 2012). The situation may be further complicated if population genetic structure exists, as state allocations would be affecting different spawning stocks.
Past studies on population structure in summer flounder have found differences in allozymes, morphology, meristic traits, and timing of shallow water ingress between individuals caught north and south of Cape Hatteras, suggesting the existence of multiple spawning stocks, or subpopulations (Able, Matheson, Morse, Fahay, & Shepherd, 1990;Able et al., 2011;Burke, Monaghan, & Yokoyama, 2000;Kraus & Musick, 2000;Van Housen, 1984). However, another study found no genetic differentiation between supposed subpopulations after examining mitochondrial DNA (Jones & Quattro, 1999), so summer flounder population substructure remains equivocal.
With the increasing ease of sampling hundreds to thousands of loci across the genome, reduced representation sequencing in summer flounder can provide a more robust estimate of population structure, as well as assess genetic variation across spatially divergent environmental conditions.
In this study, we use double-digest restriction-site-associated DNA (ddRAD) sequencing to generate 1,137 single nucleotide polymorphisms (SNPs) genotyped in 232 adult summer flounder individuals captured along the U.S. east coast. Using these data, we address two questions: (i) Do summer flounder exhibit population structure and (ii) are there particular loci that suggest locally divergent environmental selection? Our results contribute to a growing number of studies suggesting that a balance between multiple evolutionary forces shapes the genetic makeup of marine populations and that understanding these forces is important for the persistence and management of species under future climate change.

| Sample collection
Regional bottom trawl surveys in 2013-2014 collected adult summer flounder from Massachusetts to Florida (Table 1 and  were classified as mature, but information on maturity was not available for the SEAMAP nor NCDNR individuals. Previous published estimates of median length of maturity for summer flounder range from 28 to 32 cm for females and 24.9 to 28.9 cm for males (Morse, 1981;O'Brien et al., 1993;Wenner et al., 1990). Based on these estimates, some SEAMAP and NCDNR individuals may have been immature, but mature and immature summer flounder move offshore and overwinter together (Packer et al., 1999). For all samples (n = 306), a small muscle tissue plug and/or fin clip was removed and preserved in 95% ethanol. Tissue and fin clips were taken upon capture for NEFSC FBTS and NCDNR specimens.
SEAMAP specimens were shipped frozen to Rutgers University, where tissue and fin clip samples were taken after thawing. All tissue and fin clip samples were held at -20°C until DNA extraction.

| DNA extraction
To extract whole genomic DNA from adult summer flounder, ≤ 20 mg of muscle tissue or fin was lysed, and DNA was washed and eluted using DNeasy 96 Blood & Tissue Kits (QIAGEN; Hilden, Germany) and the manufacturer's recommended protocols. DNA extracts were visualized on 2% agarose gels to assess quality and were subsequently quantified using PicoGreen (Thermo Fisher Scientific, Waltham, MA, USA) and a SpectraMax M3 Microplate Reader (Molecular Devices; Sunnyvale, CA, USA).
TA B L E 1 Sampling year, season, source, latitudinal range, sample size, and range of total lengths (TL) for each collection of adult summer flounder. The range of total lengths reflects available data for each collection

| Bioinformatics and genotyping
To distinguish between pooled libraries, sequenced reads were demultiplexed by Illumina index using a Python script adapted from FASTX Barcode Splitter in the FASTX toolkit programs (Gordon, 2011). To determine reads from particular individuals, sequenced reads were further demultiplexed by barcode and cleaned using process_radtags in STACKS v.1.29 (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013). Next, sequences were run through the dDocent v. 2.14 pipeline (Puritz, Hollenbeck, & Gold, 2014), an analysis software for ddRAD sequencing data that trims, assembles, maps, and calls SNPs. In brief, using default dDocent settings unless otherwise noted, reads were trimmed for quality using TrimGalore! De novo single-end assembly of reference sequences was performed with Rainbow using alleles with a minimum within-individual coverage level of 2.5 and a minimum occurrence in 15 individuals.
Reference sequences with ≥90% similarity were clustered together with CD-HIT. Individual quality-trimmed reads were then mapped to the reference sequences using BWA, and SNPs were identified with FreeBayes.
Variant SNPs that were successfully genotyped in at least 50% of individuals with a minor allele count of at least three were retained in the analysis. To remove individuals that did not sequence well, any genotypes with fewer than three reads were recoded as "missing." About 15% of individuals with the most missing data, corresponding to individuals with more than 77% missing data, were discarded. This filter left 241 of 280 individuals for further analysis. Data were then restricted to variants occurring in 95% of individuals with a minor allele frequency (MAF) of 0.05 and a minimum mean depth of 20. Further filtering was conducted using a filtering script distributed with dDocent that filters reads based on allele balance, locus quality, mapping quality, depth, and strand representation. First, variable sites were removed if the average allele balance at heterozygous genotypes was less than 25%, the ratio of quality score to depth was less than 0.2, and the ratio between the mean mapping quality of the alternate and reference alleles was less than 0.25 or more than 1.75. The average depth and standard deviation were then calculated across all individuals for all remaining sites. Sites with depths greater than the average depth plus one standard deviation were removed if the quality score was less than double the depth. At last, to filter out sites with unreasonably high mean depth, we calculated mean depth and kept the lower 95% of the distribution, corresponding to a mean depth of 88 or less. Indels were then removed, and only the first SNP at each contig was retained to help ensure an unlinked dataset. Hence, we will hereafter refer to these SNPs as loci.
From initial investigations of the dataset, we discovered 10 very closely related summer flounder based on Rousset's distance that were caught in the same tow. No other tow exhibited a similar degree of genetic similarity. These fish were not siblings based on analysis in COLONY v.2.0.6.3 (Jones & Wang, 2010), but could be the result of cross-contamination in the trawl net.
We excluded 9 of 10 of these individuals for the analyses presented in the paper, but also tested whether our findings were sensitive to this choice by performing the analyses with all individuals included. Even with all individuals included, we found no genomewide population structure and the majority of candidate loci were associated with bottom temperature ( Figure S1 and Table S2). Overall, our filtering resulted in 1,137 loci across 232 individuals for analysis.

| Identifying population structure
We used three approaches to investigate spatial genetic structure in our summer flounder dataset: principal component analysis (PCA) followed by two analyses using Bayesian clustering algorithms implemented in the STRUCTURE and Geneland programs. We performed PCA using the adegenet v. 2.0.0 (Jombart, 2008) and ade4 v.1.7-2 packages (Chessel, Dufour, & Thioulouse, 2004) in R version 3.2.2 (R Core Team 2015), and any missing data (0.57%) were interpolated using mean allele frequency.
Next, we used STRUCTURE v.2.3.4, a Bayesian clustering algorithm that is not dependent on location information (Pritchard, Stephens, & Donnelly, 2000), to identify the number of putative populations. In STRUCTURE, we used a burn-in of 10,000 iterations followed by an additional 200,000 Markov chain Monte Carlo (MCMC) steps assuming admixture and correlated allele frequency models with prior sampling location information (Hubisz, Falush, Stephens, & Pritchard, 2009). We selected the admixture and correlated allele frequency models because they are considered to be the best when subtle population structure is expected (Falush, Stephens, & Pritchard, 2003). We initially analyzed all available 1,137 loci, and then subsequently analyzed the 1,005 that remained after removing loci not in Hardy-Weinberg proportions (HWP; p < 0.01, exact test). Loci not in HWP were identified with the pegas v. 0.10 package (Paradis, 2010) in R. We ran 10 replicates of K from 1 to 5, where K is the number of population clusters, for each set of loci and checked for parameter stabilization using the burn-in and MCMC lengths above.
To determine the optimal K, the 10 STRUCTURE replicates for each K were packaged and input into STRUCTURE HARVESTER, a web-based program that facilitates visualization of STRUCTURE output and selection of the optimal K (Earl & VonHoldt, 2012). The mean likelihood values (L(K)), as well as ΔK, a second order rate of change of L(K) with respect to K, were used to determine the optimal number of populations. The ΔK method of determining the optimal K corresponds to the true value of K more often than the L(K) method, but it cannot evaluate scenarios when K = 1, limiting its usefulness when a single population exists (Evanno, Regnaut, & Goudet, 2005).
To visualize the percent assignment of each individual into the putative populations, output from STRUCTURE HARVESTER was fed into CLUMPP, a program that averages the mean of the posteriors from multiple runs of the optimal K (Jakobsson & Rosenberg, 2007).
As different Bayesian clustering methods can produce different results (Frantz, Cellina, Krier, Schley, & Burke, 2009), we also used the Geneland v.4.0.7 package in R (Guillot, Mortier, & Estoup, 2005) to avoid drawing conclusions with potentially important conservation implications on a single analysis. Geneland is a Bayesian clustering program that utilizes geo-referenced multi-locus genetic data to determine population structure. Latitude and longitude were converted to UTM and standardized to zone 18. With the full dataset of 1,137 loci in 232 individuals, we used 100,000 MCMC iterations thinned to every 100th iteration, allowed the number of K clusters to vary between 1 and 10 and applied the spatial method with the uncorrelated frequency model. We then performed an additional analysis using the correlated frequency model and examined how the posterior distribution changed (as suggested by the Geneland manual).

| Isolation by distance
Pairwise genetic differences among 232 individuals using 1,137 loci were calculated using Rousset's distance (Rousset, 2000) in SPAGEDI (Hardy & Vekemans, 2002): where Rousset's distance (a r ) is the difference between the probability of identity of two genes within an individual (Q w ) and the probability of identity of two genes at some geographic distance r (Q r ),  Hatteras (POPONE and POPTWO, respectively) were drawn from log-uniform distributions from 100 to 100,000. Dispersal between these two populations (DISP) was measured as the probability of an individual dispersing to the other population, and this probability was drawn from a uniform distribution from 0 to 0.5. As low frequency SNPs were discarded during the bioinformatics process, the observed data had a MAF cutoff of 0.043. To match that of the observed data, simulated SNPs with a MAF of < 0.043 were also removed. We then downsampled the simulated loci to 663 and computed the joint SFS for each simulation.

| Estimating dispersal
There is still debate about how many summary statistics are appropriate for ABC analysis. Using too few summary statistics is likely to bias ABC estimates (Beaumont et al., 2002;Marjoram, Molitor, Plagnol, & Tavare, 2003), but highly dimensional datasets can also reduce the accuracy of estimation (Blum, 2010). We chose to reduce the dimensionality of our data with PCA. To avoid computational limits, the initial 10,000 simulations were used to define PCA axes in R using the prcomp function and the remaining 90,000 simulations and the observed data were projected onto these PCA axes. We examined the robustness of ABC parameter estimates based on the number of summary statistics used and found that the number of summary statistics employed had little effect on dispersal parameter estimation (Table S3). As we wished to maximize the proportion of variance explained in the data, we present results that retained all 10,000 principal components as summary statistics in the ABC analysis.
Approximate Bayesian Computation model selection was performed using the abc package (Csilléry, François, & Blum, 2012) in R.
Using simple rejection sampling, we accepted the 500 (0.5%) simulations with the shortest Euclidean distance between the simulated and observed summary statistics (the PCA-projected SFS). Using the density function in R, we plotted the posterior probability distributions comprised of these retained simulations for each parameter.
The POPONE and POPTWO parameters were log 10 -transformed for plotting. This allowed us to estimate our parameters of interest (POPONE, POPTWO and DISP) and to determine the degree to which summer flounder were panmictic.

| Environmental associations using BayEnv 2.0 and redundancy analysis
To examine the possibility of locally divergent selection in summer flounder across the species range, we used two methods to identify correlations between allele frequencies and each of four We used the same 232 summer flounder as for the ABC analysis, All 1,137 loci were examined separately for associations with standardized environmental variables using 500,000 iterations for each locus. BayEnv 2.0 was run 10 times to check for consistency between independent runs, as MCMC methods can be sensitive to initial conditions (Blair, Granka, & Feldman, 2014 (Oksanen et al., 2016) in R on standardized environmental data and a centered allele frequency dataset containing only a single allele at each of 1,137 loci across our 232 summer flounder adults. Similar to PCA, any missing data were interpolated using mean allele frequency. Potential outlier loci were identified as loci with scores ± 3 standard deviations from the mean axis score for each of the first three constrained axes. Potential outlier loci were then tested for association with environmental variables by regressing allele frequencies against each environmental variable and calculating the correlation coefficient. We set the cutoff for a significant relationship at p-value < 0.001. To test the significance of locus-environment associations identified by RDA, we randomized the environmental data 10,000 times and identified outliers and associations with the environment as described above. In addition, we relaxed the assumption of linearity by fitting generalized additive models (GAMs) between the RDA outlier loci and each environmental variable. We looked for significant locus-environment associations (p-value < 0.001) in the same way as for the linear models.
We used BLAST annotation to explore the gene function of contigs containing environmentally associated loci from both the BayEnv 2.0 and RDA analyses using the megablast and blastx algorithms with the nucleotide collection (nr/nt) and nonredundant protein sequences (nr) databases, respectively. We also used the megablast algorithm against the Paralichthys olivaceus genome.

| Identifying population structure
The PCA of 1,137 loci sampled from 232 summer flounder adult individuals revealed no distinctive population structure along the U.S. east coast between Massachusetts and Florida ( Figure 2). The first and second principal components explained 1.42% and 1.18% of the variance, respectively.
Analysis in STRUCTURE revealed minimal population structure, with no clear geographic boundaries. After testing K = 1 to K = 5, both the mean likelihood values (L(K)) and the Evanno method (ΔK) suggested two clusters (K = 2) regardless of the set of loci (all 1,137 or 1,005 passing HWP filters) used in the analysis (Figures 3 and S2).
We initially used all available 1,137 loci to assign 232 individuals to populations. We did this to avoid prematurely removing loci potentially contributing to population structure. We then removed 132 loci not in HWP and performed an additional STRUCTURE analysis using the 1,005 loci that passed HWP filters, as is more typical of population assignment analyses. Even though K = 2 using the full dataset of 1,137 loci and the subset of 1,005 loci that passed HWP filters, we interpreted these STRUCTURE analyses as a lack of population structure in summer flounder because individuals did not group into separate populations, but were instead admixed at approximately the same proportions.
Geneland analysis using the uncorrelated frequency model indicated greatest posterior density at one cluster (K = 1). The correlated frequency model in Geneland is better at detecting subtle population structure, but is more prone to instability and sensitive to departure from model assumptions. The correlated frequency model had greatest posterior density at K = 10, but clustering was not geographically related. In comparison with the uncorrelated frequency model, we interpreted the correlated frequency model results as model overfitting due to the limitations of the correlated frequency model. Thus, we conclude that Geneland analyses also suggest a single population in summer flounder ( Figure S3).

| Isolation by distance
We found no detectable relationship between geographic distance and genetic distance, as measured by Rousset's distance (a r ), across 232 individuals genotyped at 1,137 loci (r = 0.026; Mantel test p > 0.86; Figure 4). Distant summer flounder were as genetically similar as those nearby to each other.

| Estimating dispersal
Approximate Bayesian Computation analysis revealed that effective population size both north and south of Cape Hatteras (POPONE & POPTWO) and the dispersal rate between these two groups (DISP) are unlikely to be small ( Figure 5 and

| Environmental associations
Our summer flounder samples were distributed over a wide geographic area, characterized by gradients of environmental variables ( Figures S4 and S5). BayEnv 2.0 analyses revealed that 14 of 1,137 loci (Table 3; Figure S6) had strong positive associations with at least one environmental variable (median BF > 3). In total, bottom temperature was correlated with 10 loci, the greatest number for any environmental variable tested. Distance along the coast was correlated with nine loci, and both depth and bottom salinity with six loci.
Analysis of randomized environmental data in BayEnv 2.0 identified that our null expectation should be five loci associated with depth, four with bottom salinity, two with bottom temperature, and one with distance. For each environmental variable, we found the empirical cumulative observed distribution of BFs to be significantly different than the permuted distribution (two-sample K-S test; p < 0.00001; Figure S7). We also binned BFs to determine if the observed number of BFs > 3 was greater than would be predicted by chance alone for each environmental variable ( Figure S8). There were significant differences between the number of BFs > 3 for distance and for bot-  Figure   S9).
RDA identified 23 outlier loci, of which four were significantly (p < 0.001) associated with one or more environmental variables ( Figure 6 and Table 4). Five locus-environment associations were detected: bottom salinity was correlated with two loci, and bottom temperature, distance and depth were each associated with one.
The GAMs identified the same five locus-environment associations ( Figure S10). Three of the four spatially divergent RDA candidate loci also emerged as candidates in the BayEnv 2.0 analysis (contigs 8420, 35399 and 54409). Although the total number of RDA outliers we identified was to be expected under a null model (p = 0.362; Figure S11A), the number significantly associated with an environmental variable was not. We identified five significant locus-environment associations using RDA and regression, and this was unlikely to occur under a null model (p = 0.044; Figure   S11B).
Using BayEnv 2.0 and RDA, we identified a total of 15 loci whose allele frequencies were regionally differentiated and strongly associated with environmental variables (Table S4).
BLAST annotation of these 15 candidate loci using the megablast  and blastx algorithms revealed several strong matches (e-value ≤ 1E-4) with fish species (Table S5). Six of the candidate loci with strong matches aligned closely with the olive flounder (Paralichthys olivaceus) genome and have putative functions associated with calcium binding, metal binding, demethylation, transcription, and transmembrane proteins. BLAST annotation against the Paralichthys olivaceus genome also returned many annotated genes that were either embedded within or adjacent to candidate loci (Table S5).

| D ISCUSS I ON
Marine populations were often assumed to be homogenous because there are few apparent barriers to dispersal, but recent evidence instead suggests that dispersal is often highly constrained.
Using 1,137 loci genotyped in 232 adult summer flounder along nearly the entire U.S. east coast, we found that these fish comprise a single ecological population with no indication of isolation by distance. However, summer flounder reside across wide environmental gradients, and we also found genomic signatures suggestive of divergent environmental selection on regional scales, despite sequencing only a small part of the genome. In particular, we found 15 loci whose allele frequencies were associated with at least one environmental variable, and 11 (73%) of them were correlated with bottom temperature.

| Population differentiation
Renewed interest in understanding summer flounder population biology in light of shifts in biomass (Bell et al., 2014;Nye et al., 2009;Pinsky & Fogarty, 2012) provides strong motivation for the application of molecular techniques to fisheries-relevant research. From PCA, STRUCTURE, and Geneland analyses, we conclude that summer flounder are a genetically homogenous population across much of the genome. Even if some of our sampled individuals were immature, we expect population structure would be detected because juveniles and adults move offshore together during the spawning season (Packer et al., 1999). A naïve interpretation of our STRUCTURE analyses would suggest multiple populations, but inconsistencies in the number of clusters inferred by different programs can occur (Frantz et al., 2009) and it is important to interpret results with biological relevance in mind (Meirmans, 2015). In addition, while the correlated allele frequency models implemented in STRUCTURE and Geneland can help distinguish between closely related populations, this also makes them more prone to overestimating K (Falush et al., 2003;Pritchard et al., 2000). Using ABC analysis to examine population differentiation from another perspective, we estimated the dispersal rate across Cape Hatteras, North Carolina. Cape Hatteras has been suggested as an important dispersal barrier in this species (Able et al., 2011;Burke et al., 2000;Kraus & Musick, 2000;Van Housen, 1984). However, the ABC analysis suggested that the probability of an individual dispersing across this putative barrier was not low, and in fact, was likely quite high (7% to 50% probability per generation). Dispersal rates of ~10% or greater are often suggested as sufficient to synchronize population dynamics (Hastings, 1993), and a population can be defined under an ecological paradigm by focusing on demographic connectivity (Waples & Gaggiotti, 2006).
Thus, the dispersal rates that we estimated would imply that summer flounder north and south of Cape Hatteras comprise one population from an ecological and population dynamics perspective.
Previous studies investigating population structure in summer flounder using genetic techniques (Jones & Quattro, 1999;Van Housen, 1984) or phenotypic traits (Able et al., 1990(Able et al., , 2011Burke et al., 2000;Kraus & Musick, 2000) have disagreed on the existence of subpopulations. Our genetic analyses with 1,137 SNPs had higher statistical power than earlier mtDNA analyses (Jones & Quattro), but reached the same general conclusion that summer flounder exhibit extensive gene flow across the entire U.S. east coast. Earlier observations of differences in the timing of larval ingress (Able et al., 2011), morphology, and meristics (Burke et al., 2000) suggest that summer flounder may respond plastically to environmental differences across the species range, resulting in differences in phenotype despite extensive gene flow. An alternative hypothesis is that differences in phenotype are controlled or linked to loci under selection, potentially some of the loci identified by this study that were associated with an environmental variable. A productive avenue for future research would be to employ common garden studies to determine the underlying basis of phenotypic variation in this species.

| Evidence for spatially divergent selection
Despite technological advances in an era of genomics, identifying adaptive divergence in natural populations is still challenging (Bernatchez, 2016). Yet, reduced genome representation techniques like ddRAD sequencing can facilitate insight into the demographic and evolutionary processes that shape the genomes of nonmodel organisms (McKinney, Larson, Seeb, & Seeb, 2017). Given that this study surveyed ~ 0.04% of the summer flounder genome (based on the 546 Mb olive flounder genome; Shao et al., 2017), it may have missed many adaptive polymorphisms (Lowry et al., 2017). In addition, recent research suggests that environmental adaptation affects a small fraction of the genome (Bay et al., 2017), and that adaptive traits are often controlled by many loci of small effect (Le Corre & Kremer, 2012;Rockman, 2012;Yeaman, 2015). These two factors can further complicate the identification of loci involved in adaptation. However, despite limited ability to detect polygenic adaptation and loci with small effect sizes, locus-environment associations are still useful for uncovering candidate loci of relatively large effect. Our study likely missed adaptive polymorphisms due to limited genome sampling and limited power to detect small effect loci, but we still detected a number of important candidate loci that appear to be affected by spatially varying natural selection. Furthermore, the lack of overall genetic population structure suggests that the species is locally divergent only at particular loci and that our environmentally associated loci are not the result of underlying neutral population structure. Our findings, which are suggestive of adaptive divergence in summer flounder, are consistent with a growing number of other studies in marine species that have documented adaptive polymorphisms across the species range despite high gene flow (De Wit & Palumbi, 2013;Gagnaire, Normandeau, Côté, Hansen, & Bernatchez, 2012;Gleason & Burton, 2016;Pespeni, Chan, Menge, & Palumbi, 2013;Sandoval-Castillo, Robinson, Hart, Strain, & Beheregaray, 2018;Schmidt & Rand, 1999;Therkildsen et al., 2013).
Even though high gene flow can introduce maladaptive alleles to local populations and reduce the degree to which they can locally adapt, genetic differentiation at particular loci can still arise if selection is strong enough (Hedgecock, 1986;Marshall, Monro, Bode, Keough, & Swearer, 2010;Slatkin, 1987). Selection can more effectively filter out maladaptive alleles or increase the frequency of beneficial alleles when the effective population size is large and thus genetic drift is weaker (Lanfear, Kokko, & Eyre-Walker, 2014;Ohta, 1992), as was suggested for summer flounder by our ABC analysis.
The process of strong selection on a local scale following widespread dispersal can result in the maintenance of adaptive polymorphisms.
These polymorphisms are maintained through spatial balancing (also called spatially varying) selection and have been termed spatially balanced polymorphisms (Sotka, 2005(Sotka, , 2012. In an alternative way, low dispersal or active habitat selection can also drive genetic differentiation across space (Kawecki & Ebert, 2004), which is often known as local adaptation. Although the scale of gene flow and selection differ between local adaptation and spatial balancing selection, they lie on a continuum and both result in genetic differences between local sites (Sanford & Kelly, 2011). In the case of summer flounder, selection following dispersal in each generation appears to be the most likely process maintaining genetic differentiation at certain loci. Evidence for locally adapted polymorphisms maintained by spatial balancing selection in summer flounder can be strengthened by comparing allele frequencies in larvae to those in adults, which provides a promising avenue for future research. If allele frequencies at a set of loci are homogenous in larvae, but divergent in adults, this would suggest that genetic differentiation in adults is maintained by spatial balancing selection.
Processes other than spatial balancing selection could result in polymorphisms throughout the genome, but we believe these are less likely in the summer flounder system. One factor that can drive differentiation at particular loci is sex-specific genetic markers. Phenotypic sex in summer flounder is genetically determined through male heterogamety, or a XX/XY sex-determining system (Colburn, 2008). We did not find any SNPs in our dataset that were present in all males and completely absent in females, as would be expected of male-specific markers in a male heterogeneous species. However, our ability to identify such markers in our dataset is somewhat limited. Studies explicitly seeking to identify sex-specific genetic markers typically aim to understand genetic sex determination and the evolution of sex-determining systems (Carmichael et al., 2013;Fowler & Buonaccorsi, 2016;Gamble & Zarkower, 2014). As a result, these studies undertake bioinformatics steps to limit RAD tags to those found in only one phenotypic sex and not the other (which we did not do). In addition, ~20% of individuals were not confidently sexed in our dataset, and we did not aim to explicitly identify sex-linked markers. Despite these caveats, it seems unlikely that any sex-linked markers were retained in our final dataset. Demographic processes like range expansion and reproductive sweepstakes can also result in spatial differences at a few loci.
Recent range expansion can bring alleles into newly colonized areas, resulting in allele frequency gradients due to genetic drift as the expansion front moves across the landscape. This process is termed allelic surfing and can result in patterns at neutral loci that resemble signatures of selection (Excoffier & Ray, 2008;Hallatschek, Hersen, Ramanathan, & Nelson, 2007). However, allelic surfing is more likely to happen with low gene flow (Excoffier & Ray, 2008;Klopfstein, Currat, & Excoffier, 2006), which makes it unlikely in summer flounder because of high dispersal estimates and no recent history of range expansion. High reproductive variability in broadcast spawners, such as summer flounder, can also create temporal patterns of genetic differentiation (Hedgecock & Pudovkin, 2011). This process is mostly applicable to benthic adult populations and can result in transient "chaotic patchiness," or chaotic patterns of genetic divergence across space in each time step. While adult summer flounder are capable of movement, our candidate loci were associated with environmental gradients rather than exhibiting chaotic patterns typical of sweepstakes reproduction. Instead, we suggest that spatially varying selection due to environmental heterogeneity along the coastline is maintaining adaptive divergence at particular loci in summer flounder, but further investigation is necessary to strengthen the evidence for this.
Environmental conditions can be important selective agents that shape the genotypic and phenotypic composition of local populations (Sanford & Kelly, 2011). Water temperature is perhaps the most important abiotic factor defining a marine organism's habitat (Angilletta, 2009). For many ectotherms, the temperature of their surrounding environment controls a variety of physiological processes, including aerobic scope (Pörtner & Knust, 2007) and developmental (O'Connor et al., 2007) and metabolic rates (Johnston & Dunn, 1987), which in turn can affect growth and survival. In a classic example, Powers and Schulte (1998)

| Potential for adaptation in a changing world
Within the context of global environmental change, spatial balancing selection can allow for all phenotypes and their corresponding genotypes to be present in populations across the species range. Selection acts on the local scale to filter for phenotypes and genotypes that are best suited to local conditions. As local conditions shift over time, different phenotypes and their underlying genotypes may become more or less favorable depending on current conditions. The ability of selection to match phenotype with the current environmental conditions results in a high potential for adaptation as the environment changes (Sanford & Kelly, 2011). Of course, adaptive potential also depends on the underlying genetic variation in a species. Historical studies of exploited fish stocks have shown that fishery-induced population declines can lower effective population size (N e ) and erode standing genetic variation (Hauser, Adcock, Smith, Ramiréz, & Carvalho, 2002;Hutchinson, van Oosterhout, Rogers, & Carvalho, 2003;Pinsky & Palumbi, 2014). If existing genetic variation and effective population size are low in summer flounder due to historical fishing pressure, adaptive evolution in response to changing environments may be limited. Future studies investigating the relationship between standing genetic diversity and the effectiveness of selection at matching phenotypes and genotypes to local environmental conditions would contribute toward a better understanding of evolutionary potential in response to environmental change.
In conclusion, summer flounder from Massachusetts to Florida comprise an effectively panmictic population across much of their genome. We detected genomic signatures suggestive of spatially divergent environmental selection at a few select loci, with many of these divergent loci associated with bottom temperature. These loci are good candidates for further investigation into their functional role in adaptation to changing environmental conditions, and they appear to be maintained by selection rather than limited gene flow. Understanding the relevance and balance between dispersal potential, the strength of selection, and the scale of environmental heterogeneity can help us to understand the adaptive potential and persistence of summer flounder and other important exploited species under future climate conditions.
Stuart assisted in learning and troubleshooting the ddRADseq protocol and bioinformatics.

CO N FLI C T O F I NTE R E S T
None declared.