Selection on outlier loci and their association with adaptive phenotypes in Littorina saxatilis contact zones.

doi: 10.1111/jeb.12564

Selection on outlier loci and their association with adaptive phenotypes in Littorina saxatilis contact zones J. HOLLANDER*†, J. GALINDO*‡ & R. K. BUTLIN*§ *Department of Animal and Plant Sciences, University of Sheffield, Sheffield, UK †Department of Biology, Aquatic Ecology, Lund University, Lund, Sweden ‡Departamento de Bioquımica, Xenetica e Inmunoloxıa, Facultade de Bioloxıa, Universidade de Vigo, Vigo, Spain §Sven Loven Centre for Marine Sciences – Tj€arn€o, University of Gothenburg, Str€omstad, Sweden

Keywords:

Abstract

amplified fragment length polymorphism; association analysis; gene flow; hybrid zone; Littorina saxatilis; speciation.

A fundamental issue in speciation research is to evaluate phenotypic variation and the genomics driving the evolution of reproductive isolation between sister taxa. Above all, hybrid zones are excellent study systems for researchers to examine the association of genetic differentiation, phenotypic variation and the strength of selection. We investigated two contact zones in the marine gastropod Littorina saxatilis and utilized landmark-based geometric morphometric analysis together with amplified fragment length polymorphism (AFLP) markers to assess phenotypic and genomic divergence between ecotypes under divergent selection. From genetic markers, we calculated the cline width, linkage disequilibrium and the average effective selection on a locus. Additionally, we conducted an association analysis linking the outlier loci and phenotypic variation between ecotypes and show that a proportion of outlier loci are associated with key adaptive phenotypic traits.

Introduction The population genomic approach (Luikart et al., 2001; Butlin, 2010) has been widely used in recent years to seek genomic regions implicated in local adaptation or reproductive isolation. Typically, these studies identify a few percentage of loci that are ‘outliers’, that is loci with higher genetic differentiation than expected based on the mean divergence and some assumptions about population history and gene flow (Nosil et al., 2009). The latest studies have many markers and so high genomic resolution (e.g. Hohenlohe et al., 2010; Martin et al., 2013), but nevertheless, most analyses lack estimates of the strength of selection on outlier loci and do not associate outlier loci with phenotypes under selection (Seehausen et al., 2014). Potentially, both of these limitations can be overcome by combining outlier analysis with the analysis of contact zones. Hybrid zones provide excellent models for the study of speciation and its genetic architecture (Hewitt, Correspondence: Johan Hollander, Department of Biology, Aquatic Ecology, Lund University, Lund, Sweden. Tel.: +46 46 222 34 73; fax: +46 46 222 45 36; e-mail: [email protected]

1988; Barton & Hewitt, 1989; Bridle & Ritchie, 2001; Buggs, 2007). They are characterized by clinal changes in phenotypic traits and in allele frequencies that are maintained by a balance between gene flow and selection. Intrinsic sources of selection, such as genetic incompatibilities accumulated during a period of allopatry, have often been emphasized, but very similar clinal patterns are expected across abrupt environmental transitions (Kruuk et al., 1999). Only when the environmental transition occurs over a distance much greater than the characteristic scale (r/√s, where r is the standard deviation of parent–offspring distances and s is the selection per locus) are clines expected to be independent of dispersal (Barton & Hewitt, 1985). In practice, intrinsic and extrinsic sources of selection are difficult to distinguish and both are expected to contribute to most hybrid zones (Bierne et al., 2011). Within a hybrid zone, dispersal can generate linkage disequilibria, even between loci that are not physically linked. This linkage disequilibrium can be used to estimate gene flow which, with the width of a cline, can be used to estimate the effective selection experienced by a locus (Barton & Hewitt, 1985; Barton & Gale, 1993).

© 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. JOURNAL OF EVOLUTIONARY BIOLOGY © 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

1

2

J. HOLLANDER ET AL.

In several cases, population genomic analyses have been combined with genetic mapping of adaptive traits to associate outlier loci with phenotypes and so with sources of selection (e.g. Rogers & Bernatchez, 2005; Via & West, 2008; Hohenlohe et al., 2010). However, this approach is limited to systems where the analysis of controlled crosses is possible. An alternative is to exploit the association between alleles that influence phenotypic traits and linked marker alleles that persists following population admixture or is maintained by gene flow into hybrid zones. Although this ‘admixture analysis’ has been used frequently in human genetics (Hindorff et al., 2014), it has only recently been applied in the context of local adaptation and speciation. Malek et al. (2012) analysed admixture following recent hybridization between two divergent stickleback ecotypes within a lake, whereas Lindtke et al. (2013) used a hybrid zone between Populus species to associate loci under selection with adaptive phenotypes. Littorina saxatilis is a marine gastropod that has been studied mainly in Sweden, Britain and Spain with the focus on local adaptation and the evolution of reproductive barriers (reviewed in Johannesson et al., 2010). Within each region, specialized ecotypes occur in contrasting habitats dominated either by crab predation or by hydrodynamic wave force. In Sweden, ecotypes living on boulder shores with abundant crabs (crab or S-ecotype [sheltered]) have a large size, thick shell and a small aperture to withstand crab attacks. The other ecotype in rocky headland environments (wave or E-ecotype [exposed]), which lack crabs, is two to four times smaller than the crab ecotype, thin shelled and has a large aperture to allow an increase in foot size for improved adhesion on the rock. Overall divergent selection between the environments is strong. Estimates of survival by Janson (1983) suggest relative fitness of the crab ecotype in the wave environment of 0.28 and of the wave ecotype in the crab environment of 0.42. As these estimates include only a part of the life cycle and only the survival component of fitness, they must underestimate the overall strength of selection. Between the boulder shore and the smooth rock, the two ecotypes meet at abrupt habitat boundaries, with 5–20 m of intermediate habitat. Despite a tendency to mate assortatively (Hollander et al., 2005), the ecotypes interbreed, generating narrow hybrid zones containing individuals of intermediate phenotype and genotype (Panova et al., 2006) which might be at an advantage in intermediate habitats (Janson, 1983; Johannesson, 2003). The hybrid zones studied here are unlikely to be the result of secondary contact after a long period of spatial isolation. A model of in situ divergence of ecotypes was favoured over a model of separate colonization (Butlin et al., 2014), and background levels of genetic differentiation between ecotypes are low compared with spatial genetic structure (Panova et al., 2006). It is not possible to exclude periods of geographic separation in the past or

contributions of intrinsic incompatibility to selection against hybrids (Bierne et al., 2011, 2013), but we expect the greatest contribution to selection to come from differential adaptation. Although total selection against each ecotype when transplanted to the alternative habitat is strong (Janson, 1983), we expect this to be spread across multiple traits each determined polygenically such that direct selection on each locus is weak. We assume a model with divergent selection that changes abruptly at or near the transition from rocky shore to boulder field resulting in clines determined by selection–dispersal balance. We further assume that selection is not strong enough relative to recombination to generate strong linkage disequilibrium and a stepped cline (Szymura & Barton, 1991; Kruuk et al., 1999). Under these assumptions, some linkage disequilibrium is generated by dispersal, but it is still possible for cline positions and widths to vary among traits (as in the grasshopper Chorthippus parallelus, for example; Butlin et al., 1991). Previous studies of British and Spanish populations of L. saxatilis with amplified fragment length polymorphism (AFLP) markers found a small group (3–5%) of strongly differentiated loci between the two ecotypes (crab vs. wave) (Wilding et al., 2001; Galindo et al., 2009). They considered the differentiated loci to mark areas of the genome under selection for local adaptation or genomic incompatibilities between the ecotypes. These loci showed narrow clines and a small increase in linkage disequilibrium at the cline centre in Britain (Grahame et al., 2006). Also, the genetic clines corresponded closely in position with the morphological transition in both hybrid zones (Grahame et al., 2006; Galindo et al., 2013), but no formal analysis of selection strength and association was conducted in a single study. Given an estimate of gene flow between ecotypes, selection on one of the most differentiated loci has been estimated to be in the region of s = 103 (Wood et al., 2008; Butlin et al., 2014). Here, we examine two transects across environmental transitions in Swedish L. saxatilis that include parental ecotypes and hybrids. We performed a landmark-based geometric morphometric analysis to capture shell shape and size variation, and we used AFLP markers to measure genetic variation. We use these data to quantify (i) the level of heterogeneous genomic divergence between ecotypes, (ii) cline widths and linkage disequilibrium, and so the strength of selection per locus and (iii) the association between differentiated loci and phenotypes.

Materials and methods Sampling We used the same individuals as in Panova et al. (2006), sampled during spring 2003 from two islands in the Swedish Koster archipelago (Rams€ o 58°500 N, 11°040 E and Ramsholmen 58°510 N, 11°030 E). On

© 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12564 JOURNAL OF EVOLUTIONARY BIOLOGY © 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

Selection in Littorina saxatilis contact zones

each island, samples were taken from five points on a transect running from a boulder shore with the crab [S] ecotype to a rocky shore with the wave [E] ecotype. Each transect was centred in the intermediate zone (‘hybrid’ [H]), where the two habitats meet. Then, two samples of each ecotype were collected (2E, 1E, 1S, 2S) at 9 and 27 m from the centre and towards each habitat (see Panova et al., 2006 for details). Sample size was 49 or 50 for each point, on each island. Quantifying morphological shell shape variation The landmark-based geometric morphometric technique was used to assess shell shape and centroid size differences between ecotypes (Bookstein, 1991; Rohlf & Marcus, 1993; Adams et al., 2004). We captured eight anatomical landmarks (tpsDig Slice, 1994) positioned over the shell to describe morphological traits associated with selection (Hollander & Butlin, 2010) and used the software TPSRELW (Rohlf, 2003) to assess shape variation among groups. For a complete description of the procedure, see Hollander et al. (2005). The landmark data set has been submitted to DRYAD: doi: 10.5061/dryad.2jh32. DNA extraction and AFLP analysis Genomic DNA was extracted from the foot tissue with the same protocol as in Wilding et al. (2001). AFLP analysis was performed following the methodology described in Butlin et al. (2014). Briefly, DNA was digested using EcoRI and PstI for 3 h at 37 °C and then ligated to Eco and Pst adaptors for 16 h at 16 °C. Preselective and selective amplifications used the same primers as Butlin et al. (2014). The 12 selective primer combinations were analysed on an ABI 3730 sequencer (Applied Biosystems) with a GENESCAN 500 ROX size standard (Applied Biosystems, Foster City, CA, USA), supplemented with 585- and 685bp ROX-labelled marker bands generated in house. Positive and negative controls were run along with the samples. GENEMAPPER v.3.7 (Applied Biosystems) was used to analyse the AFLP profiles, including fragments (> 50 RFU, relative fluorescent units) between 50 and 685 bp. In addition, between 12% and 14% of the samples were replicated for each primer combination and mean genotyping error rate per locus (see Pompanon et al. 2005) was 8.2 0.4%. The AFLP data set has been submitted to DRYAD: doi: 10.5061/dryad.2jh32. Data treatment We created a relatedness matrix from the AFLP data using AFLP-SURV v1.0 (Vekemans et al., 2002) and used a multidimensional scaling (MDS) plot to detect outlier individuals, using the statistical package R (R Core Team, 2013). Typically, these were individuals with a low proportion of ‘band present’ phenotypes due to weak amplification. After discarding these outliers, we repeated the

3

AFLP-SURV

analysis with the new trimmed data set (464 individuals) and performed a new MDS analysis to check that no further outlier individuals were revealed. Exploring loci under selection

To detect outlier loci, we used the methodology developed by Beaumont & Nichols (1996) with the software package DFDIST (http://www.maths.bris.ac.uk/~mamab/ stuff/). We chose this simple method for comparability with previous analyses of AFLPs in L. saxatilis (Wilding et al., 2001; Galindo et al., 2009) and because our objective was simply to provide an objective cut-off for differentiated loci rather than a strong test for departure from neutrality. The procedure uses an empirical estimate of average FST, calculated from our AFLP data set, and performs coalescent simulations of a classical island model (two islands in our case, as we were primarily interested in reciprocal gene exchange between adjacent populations of the two ecotypes). The method produces a simulated neutral distribution of FST, contingent on heterozygosity, estimates the quantiles of the distribution and identifies loci departing significantly from the neutral distribution as putative outliers influenced by divergent selection. We analysed Rams€ o and Ramsholmen samples independently, and we used only the samples most distant from the contact zone (samples 2E and 2S). For the DFDIST simulations, Ne was set to 1000 (following Wilding et al., 2001); we used the 30% trimmed mean FST (see Caballero et al., 2008) as the expected neutral FST and simulated 100 000 loci. We considered outliers with P-values greater than 0.975. Linkage disequilibrium We measured pairwise linkage disequilibria (D) between outlier loci within each sample using the method developed by Li et al. (2007) for dominant loci, implemented in a custom R script (R Core Team, 2013) (as in Jackson et al., 2012). This method assumes Hardy–Weinberg equilibrium (HWE), an assumption that may be violated as a result of gene flow into a hybrid zone. However, allozyme data have not revealed any consistent pattern of departure from HWE in Swedish L. saxatilis populations (Johannesson & Tatarenkov, 1997). We altered the sign of D so that positive values indicated association between alleles at high frequency within the same ecotype, rather than between presence (or absence) alleles (Grahame et al., 2006; Jackson et al., 2012). Cline analysis and strength of selection in the hybrid zone The five samples from each contact zone did not provide sufficient information to fit clines to each outlier locus independently. Therefore, we were unable to test


4

J. HOLLANDER ET AL.

whether clines at different loci are coincident and/or concordant. In line with the assumptions outlined in the Introduction, our analyses treat clines as independent even though, of necessity, we estimated the average cline width and centre across outlier loci rather than cline parameters for each locus separately. Using the RPARALLEL function in GENSTAT v.14.1 (VSN International Ltd.), we fitted a tanh function with common centre and width but distinct end frequencies for each locus: pi;x ¼ pi;WAVE þ dpi f0:5ð1 þ tan hð2ðx cÞ=wÞÞg; where c is the common centre, w is the common width, pi,WAVE is the allele frequency for locus i in the wave ecotype, dpi is the allele frequency difference between ecotypes, x is distance along the transect and pi,x is the allele frequency for locus i at distance x. Using this width to estimate selection requires an estimate of dispersal, and we obtained this from the linkage disequilibrium estimates. D for a pair of loci is expected to depend on dispersal, the rate of change in allele frequency at each locus and the recombination rate between loci (Barton & Gale, 1993; eqn 1): D¼

r2 dp1 dp2 ; r dx dx

where r is the recombination rate, p1 and p2 are allele frequencies at two loci and x is the distance along the transect. This relationship does not require fixed differences between parental populations. Therefore, we regressed D for each locus pair and sample on the prod 1 dp2 uct of allele frequency slopes dp ; at that cline posidx dx tion based on the cline fits: the slope for locus i being the slope of the fitted tan h curve at position x times the difference in allele frequency between the crab and wave morphs. The regression was constrained to pass through the origin as we expect D = 0 when either slope is zero (as in Jackson et al., 2012). Note that the common cline width still resulted in different rates of allele frequency change for each locus because the frequencies at either end of the cline were fitted separately for each locus and that the cline widths relative to sample spacing meant that some linkage disequilibria may still be generated by dispersal in the extreme samples. We estimated the harmonic mean recombination rate (r) following the procedure described in Porter et al. (1997) and Macholan et al. (2007). We used r0 = 0.001, the minimum distance between loci, chromosome number of 2n = 34 (Rol anAlvarez et al., 1996) and the average chromosome length R/C = 0.5 (i.e. one chiasma per chromosome), which gave r = 0.396. Given the slope of the regression for D (b), the rate of dispersal ( can be estimated from Barton & Gale’s (1993) expectation, assuming that the samples were taken after dispersal: r¼

rffiffiffiffiffiffiffiffiffiffiffi br 1þr

The average effective selection s* (Barton & Gale, 1993; Kruuk et al., 1999) on differentiated loci in the centre of the zone was calculated from the estimated cline width (w) and dispersal rate (r), assuming a step change in the environment: r 2 s ¼ 3 : w Association analysis for shell shape variation The linkage disequilibrium generated by dispersal in the contact zone is also expected to result in associations between marker loci and QTL for phenotypic traits that differ between ecotypes (Lindtke et al., 2013). This will be detectable as correlations between marker genotypes and phenotypes within samples, over and above the gross association due to coincident clines. These associations will be strongest for marker loci close to QTL because the strongest linkage disequilibria are generated where linkage is tight (assuming that cline widths and centres are similar). Therefore, we tested whether any of the outlier loci was associated with snail size (measured as centroid size) or either of the first two relative warps describing shell shape (RW1, RW2). First, we fitted a model of variation in the phenotypic trait across the five samples for each island (GLM with term ‘sample’), and we then tested whether adding the effect of AFLP phenotype (two levels: band present or absent) at a specific outlier locus improved the model fit. If an outlier locus is closely linked to a QTL for a trait such that higher than average linkage disequilibrium is maintained between these loci, the individuals with the band-present phenotype (++ or + genotype) are expected to differ in trait value, within samples, from those with the band-absent phenotype ( genotype). We used a Benjamini & Hochberg (1995) FDR correction, across loci, to adjust significance tests.

Results The transects at Rams€ o and Ramsholmen islands, including the 10 samples of L. saxatilis, provided 1167 putative amplified fragment length polymorphism loci for 464 snails. The multidimensional scaling plot (Fig. 1) illustrates two well-separated groups, mainly on the second dimension, representing the islands of Rams€ o and Ramsholmen, as well as parallel differentiation between ecotypes, mainly on the first dimension (although ecotypes on Ramsholmen were less strongly differentiated). AFLP-SURV gave a mean FST between the extreme samples on Rams€ o of 0.0468 and 0.0351 on Ramsholmen, confirming the lower overall differentiation on the latter island. The DFDIST analysis between ecotypes detected 57 outlier loci in Rams€ o and 69 in Ramsholmen, with seven loci shared between islands (Fig. S1).



Ramsö

●

0.4

●

●

● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ●● ● ●●● ● ●● ● ●●●●● ● ● ●● ● ● ● ● ● ● ● ●● ●● ● ●● ●●● ●●● ● ● ● ● ● ● ● ● ●● ●●● ● ●●●●● ●● ●●● ● ● ●●● ●● ● ● ● ● ● ●●●●●●● ●●● ●● ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●● ●●●●● ● ● ●●●●●● ● ● ●●● ●● ● ● ●● ●● ● ● ● ●● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ●

–0.2

0.0

0.2

●

Dimension 2

●

2S 1S H 1E 2E

● ●

−0.4

Ramsholmen 2S 1S H 1E 2E −0.4

−0.2

0.0 Dimension 1

0.2

0.4

Fig. 1 Multidimensional scaling of the relatedness matrix based on amplified fragment length polymorphism genotypes for the ten samples from the two islands of Rams€ o and Ramsholmen. Samples 2S and 2E collected at 27 m from the hybrid samples (H) and 1S and 1E at 9 m.

The estimated common width of the cline for outlier loci at Rams€ o was 10.4 1.25 m shifted 3.22 0.50 m from the habitat boundary towards the S habitat. The equivalent estimates for Ramsholmen were a width of 23.8 2.84 m and a centre shifted 2.94 0.78 m from the habitat boundary, towards the E habitat (see Supporting Information for an alternative analysis using a hybrid index based on the AFLP data). The generally low differentiation between ecotypes, even at outlier loci (Fig. S1), limited our ability to make accurate estimates of linkage disequilibrium. The unweighted regression of pairwise linkage disequilibria, D, on the product of slopes was not significant in either transect. However, confidence in the estimated D values was variable, and therefore, we used weighted regression, with the likelihood ratio from the test for D 6¼ 0 as the weight. This regression was found to be significant for Rams€ o samples (b = 8.07 2.52, F1,3896 = 10.24, P < 0.001) but not for Ramsholmen, probably because the wider cline generated shallower allele frequency gradients at the centre (Fig. S2). Using the regression slope from Rams€ o provided an estimate of dispersal rate, r = 1.51 0.85 m gen0.5. Combined with the cline width estimate, this leads to an estimated mean effective selection on outlier loci of s* = 0.06 (95%CI: 0.005–0.32). On the assumption that the same dispersal rate is appropriate for the Ramsholmen population, the estimated effective selection per locus on that island was s* = 0.008. To provide a visual comparison between clines at outlier loci and phenotypic clines, we rescaled allele frequencies, RW1, RW2 and centroid size from zero to

5

one and then plotted them with the tanh curve representing the common cline centre and width for AFLP loci for each island (Fig. 2; note that rescaling was used only for the figure, not in any analyses [see Supporting Information for an alternative analysis using the STRUCTURE software]). Clines in centroid size followed allele frequency clines closely for both islands, but for the shape variable RW1, the cline appears to be shifted towards the wave/E habitat (opposite to the displacement for the allele frequency cline) on Rams€ o, whereas on Ramsholmen, the cline centres correspond, but samples taken close to the habitat transition had more extreme shapes than those taken further away. The second shape axis, RW2, showed patterns similar to RW1 on both Rams€ o and Ramsholmen, although overall differentiation was lower relative to within-sample variation, particularly on Ramsholmen. We tested for association between outlier genotypes and phenotypes (centroid size, RW1, RW2). This would be expected if divergence in outlier loci is due to their proximity to QTL for these adaptive traits. No significant associations were found for any of the phenotypic traits in Rams€ o samples, possibly because the narrower clines reduce the power of this analysis. In Ramsholmen samples, three loci were significantly associated with centroid size and one with RW2 (Table 1).

Discussion A major aim in evolutionary biology is to elucidate the mechanisms underlying the origin of reproductive barriers between closely related species. Hybrid zones provide excellent ‘natural laboratories’ in which to tackle this problem (Hewitt, 1988). In the current genomic era, they provide a way to link genetic differentiation with estimates of selection and the relevant adaptive phenotypes, critical objectives in determining the genetic basis of speciation (Seehausen et al., 2014). Here, we have illustrated this potential with an analysis of differentiation across two contact zones in the marine gastropod, L. saxatilis, providing estimates of lifetime dispersal and of the mean effective selection per locus. Significant associations, within samples, between some outlier loci and adaptive traits that differ between ecotypes suggest that the selection operating on these loci arises, at least in part, from close linkage to QTL underlying phenotypic differentiation. In recent years, L. saxatilis has become a valuable model in the marine environment to study the evolution of reproductive barriers, local adaptation and ecological speciation (Rolan-Alvarez et al., 2004; Johannesson et al., 2010; Faria et al., 2014). The species has been studied mainly in three regions of Europe: Sweden, Britain and Spain, where the ecotypes appear to have arisen separately (Butlin et al., 2014). Reciprocal transplant experiments reveal severe selection against misplaced individuals in the pure parental


6

J. HOLLANDER ET AL.

(a)

(b)

Fig. 2 Clinal patterns for amplified fragment length polymorphism outlier loci and phenotypes at (a) Rams€ o and (b) Ramsholmen. The vertical axis is scaled so that the mean for the extreme wave ecotype sample (at 27 m) is zero and for the extreme crab ecotype sample (at +27 m) is 1.0 for each trait. Error bars are 95% confidence intervals (symbols have been jittered to prevent overlap).

Table 1 The five strongest associations between amplified fragment length polymorphism loci and phenotypic traits, for each island. Centroid size

F

P

€ (residual d.f. = 186, 57 outlier loci tested) Ramso Sample (d.f. = 4) 77.465 < 0.001 Locus (d.f. = 1) 216 AGA_AAC 5.491 0.0202 466 AGA_AAG 4.703 0.0314 337 AGA_AAG 4.313 0.0392 117 AAC_AGC 3.696 0.0561 439 AGG_AAC 3.057 0.082 Ramsholmen (residual d.f. = 164, 69 outlier loci tested) Sample (d.f. = 4) 76.196 < 0.001 Locus (d.f. = 1) 80 AGA_ACC 14.159 0.0002* 97 AAC_ACC 11.937 0.0007* 315 AGA_ACT 10.774 0.0013* 134 AGA_ATG 9.769 0.0021 337 AGA_AAG 6.914 0.0094

RW1

F

P

RW2

F

P

141 114 439 340 444

AGA_ATG AGA_AGC AGG_AAC AGA_AGC AGG_AAC

75.572 7.801 5.998 5.042 4.924 4.807

< 0.001 0.0058 0.0153 0.0259 0.0277 0.0296

133 AAC_AGC 672 AGG_AAC 115 AGG_AAC 274 AAC_AGC 94 AGA_AAC

53.07 5.391 4.733 3.999 3.044 2.49

< 0.001 0.0213 0.0308 0.047 0.0827 0.1163

653 672 199 672 337

AGA_ACC AGG_AAC AGA_AAC AGA_AAC AGA_AAG

40.873 9.516 9.102 7.833 7.816 4.224

< 0.001 0.0024 0.003 0.0057 0.0058 0.0414

199 149 191 658 102

23.044 12.625 6.497 4.649 3.72 3.114

< 0.001 0.0005* 0.0117 0.0325 0.0555 0.0795

AGA_AAC AGA_ACC AAC_AGC AGA_AAC AGA_ACT

*Significant at FDR = 0.05. © 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12564 JOURNAL OF EVOLUTIONARY BIOLOGY © 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY


habitats (Janson, 1983; Rolan-Alvarez, 2007). Previous genome scans have detected a few percentage of markers to be influenced by selection (Wilding et al., 2001; Galindo et al., 2009, 2013), and a rough estimate of the strength of selection on one of these markers has been derived from the observed divergence in allele frequency and an estimate of gene flow (Wood et al., 2008; Butlin et al., 2014). Galindo et al. (2013) reported covariance between outlier loci and adaptive phenotypic traits in L. saxatilis from shores of NW Spain. However, they did not separate the association due to clinal variation from the association due to linkage, and so this study is not comparable to our association analysis and does not address the source of selection on outlier loci. No estimate of the effective selection per locus, across a sample of loci, or of the linkage between outlier loci and adaptive traits has been made to date. In this study, we discovered 5–6% of loci to be more differentiated than expected from mutation-drift-gene flow balance between the wave and the crab ecotypes of L. saxatilis in Sweden. This is within the range of previous reports for L. saxatilis from England and Spain (Wilding et al., 2001; Galindo et al., 2009, 2013) and for other organismal groups (Nosil et al., 2009). We found only seven of these outlier loci to be shared between islands (12%; Fig. S1). Galindo et al. (2013) found about 11% shared outliers between pairs of localities in NW Spain using DFDIST and a 95% cut-off, whereas Wilding et al. (2001), using a more stringent cut-off, found the majority of outliers to be shared between sites that were further apart than the islands sampled here (although connected by more continuous habitat). Our current result is in line with the majority of studies reviewed by Nosil et al. (2009), where only a minority of outliers is repeated across multiple comparisons of the same environmental contrast. This may be because adaptation has a different genetic basis in each locality or is due to independent parallel mutations in the same genes (Faria et al., 2014), but it could also reflect the high false-positive and false-negative rates of genome scans (Butlin, 2010). Evidence for the repeated use of the same genes during parallel and convergent evolution (Conte et al., 2012) suggests that the latter explanation is also likely in this hybrid zone, but further work on outlier loci is necessary before this can be verified. Panova et al. (2006) studied the same individuals for neutral microsatellite makers, and they found a similar pattern compared to our multidimensional scaling plot (Fig. 1), islands separated in one axis and ecotypes on the other. Although the latter was nonsignificant for microsatellites, reflecting the lack of outlier loci in their analyses, they showed that different alleles were responsible for the neutral genetic differentiation between ecotypes on different islands. Both results point towards a primary origin of this L. saxatilis hybrid

7

zone (see below, Butlin et al., 2014), independently of the genetic basis of the adaptation, which is unknown. In Sweden, the parental habitat of the wave ecotype is smooth rocks exposed to strong hydrodynamic forces, whereas the crab ecotype habitat consists of boulder shores occupied by the green crab (Carcinus maenas). The transition zone from one habitat to the other is typically narrow. Our hybrid zone analysis assumes a step change in selection at or near this point, corresponding approximately to the central sample on each island, and over a short distance, relative to the characteristic scale for this species. Clines are therefore expected to be ‘dispersal dependent’, with width determined by the balance between gene flow and selection, rather than ‘dispersal independent’, where allele frequencies or phenotypic means track the change in environment (Barton & Hewitt, 1985; Kruuk et al., 1999). The phenotypic transition approximately corresponds to the habitat transition on each island (Fig. 2). However, the hybrid zone analysis suggests that the centres of the genetic clines are, on average, displaced towards the S habitat (crab) at Rams€ o and in the opposite direction at Ramsholmen, nearer the E (wave) environment. The average width of the clines also differs between the islands (Fig. 2), being steeper for Rams€ o than for Ramsholmen. This may indicate different movement patterns, perhaps due to populationdensity asymmetries, or a difference in the relative strengths of selection on the two sides of the transition. Such a difference in selection may also explain why, at Ramsholmen, the most distant S sample was not strongly distinct from the hybrid zone sample. Weaker selection is also consistent with the lower overall FST between ecotypes on Ramsholmen (AFLPs, this study; microsatellites, Panova et al., 2006). However, the differences could be due to stochastic changes in cline position and width, a possibility that can only be excluded by resampling of the same contact zones. Mark–recapture studies in Sweden have previously shown that snails can move 3–4 m in a few months (Janson, 1983) (a typical adult lifetime). Our genetic estimate for the standard deviation of parent–offspring distance, a measure of lifetime dispersal, is lower at 1.5 m gen1/2. Genetic estimates of dispersal have advantages over direct estimates because they integrate over the whole lifespan and over multiple generations and do not exclude the tail of the dispersal distribution (Barton & Gale, 1993). Often genetic estimates exceed direct estimates for these reasons (but see Kawakami et al., 2009; for a counter-example), and so our low estimate is surprising (although note that the standard error is large, mainly due to uncertainty in the estimation of D). It may be biased downwards by nonrandom dispersal at the habitat boundary, for which there is some evidence in the UK (Webster, 2014) and Spain (Erlandsson et al., 1998; Cruz et al., 2004), or by


8

J. HOLLANDER ET AL.

size-assortative mating (Hollander et al., 2005; Johannesson et al., 2008). Hybrid zones are maintained in a balance between dispersal and selection, given that the selection regime changes over a small spatial scale relative to dispersal (Haldane, 1948; Fisher, 1950; Endler, 1977; Barton & Hewitt, 1985). This makes it possible to estimate the effective selection per locus from the cline width, if dispersal is known. Here, we obtain estimates of average effective selection per locus of s* = 0.06 or s* = 0.008 for the two zones, with wide confidence intervals. Combined with our dispersal estimate, this suggests a characteristic scale of 6–16 m, in the same range as the estimated width of transitional habitat (Janson, 1983; Johannesson, 2003), which justifies our use of expectations for dispersal-dependent clines (see Barton & Hewitt, 1985, p. 116). However, over this transition, it is possible that selection on different traits or loci changes at different positions resulting in variation in cline centres and widths. Unfortunately, we were not able to test for coincidence of cline centres across outlier loci. Our selection estimates are somewhat lower than estimates obtained by similar approaches in Bombina toads (s* = 0.2; (Szymura & Barton, 1991) or Vandiemenella grasshoppers (s* = 0.2; Kawakami et al., 2009) but comparable to those obtained for Mus (s* = 0.02; Rafauste et al., 2005), for example. This may reflect the likely recent primary origin of the L. saxatilis contacts where selection may be concentrated on a relatively small number of loci under divergent selection in contrast to the widespread intrinsic and extrinsic barriers accumulated in allopatry and expressed in secondary contact zones that form after long periods of divergence. When selection acts directly on few loci, the effective selection on marker loci tends to be low because they are only weakly associated with directly selected loci. Even though our estimate of selection applies to outlier loci only, rather than randomly chosen loci, it is still likely that they are only loosely associated with the direct targets of selection. If our outlier loci are influenced by selection due to their proximity to QTLs for phenotypic traits involved in differential adaptation, then trait–genotype associations should be detectable within samples, maintained by gene flow across the hybrid zones. Despite our limited sample sizes and the low power of dominant markers, such as AFLP, for such an analysis, we did find significant associations between outliers and both centroid size and the RW2 shape axis. In a hybrid zone where selection is strong relative to recombination, linkage disequilibrium may occur even between unlinked loci, particularly those loci that contribute to the barrier to gene flow and have steep, coincident clines. In such a case, the type of association we observe between phenotypes and outliers may not indicate proximity of marker loci to QTL. However, the

evidence here and in previous studies suggests that effective selection per locus is weak in Littorina. Environmental transitions may not be coincident for all traits under selection. Observed disequilibria are also weak. For these reasons, we expect selection to act independently on different loci, although data for more loci and with finer spatial sampling will be needed to test this assumption. Few studies have taken this approach to date (Malek et al., 2012; Lindtke et al., 2013), but it clearly has considerable promise, particularly with new marker types (e.g. Hohenlohe et al., 2010). If our assumptions are correct, the approach here provides further evidence that these outlier loci are influenced by selection and suggests that the origin of this selection for at least some of the loci is divergent adaptation in size and shape.

Acknowledgments We thank R. Whitlock for writing R-scripts to carry out the LD analyses based on original code provided by Li et al. We are also grateful to Ben Jackson for helping with the LD and MDS analyses. JH was funded by a Marie Curie Intra-European Fellowship (R/120197-111), a Marie Curie European Reintegration Grant (PERG08-GA-2010-276915) and by the Royal Physiographic Society in Lund. RKB is grateful for funding from NERC and the Swedish Research Council VR. JG is currently supported by an ‘Isidro Parga Pondal’ fellowship (Xunta de Galicia). Emilio Rol an-Alvarez and Nicolas Bierne provided valuable comments on an earlier version.

References Adams, D.C., Rohlf, F.J. & Slice, D.E. 2004. Geometric morphometrics: ten years of progress following the ‘revolution’. Ital. J. Zool. 71: 5–16. Barton, N.H. & Gale, K.S. 1993. Genetic analysis of hybrid zones. In: Hybrid Zones and the Evolutionary Process (R.G. Harrison, ed.), pp. 13–45. Oxford University Press, New York. Barton, N.H. & Hewitt, G.M. 1985. Analysis of hybrid zones. Annu. Rev. Ecol. Syst. 16: 113–148. Barton, N.H. & Hewitt, G.M. 1989. Adaptation, speciation and hybrid zones. Nature 341: 497–503. Beaumont, M.A. & Nichols, R.A. 1996. Evaluating loci for use in the genetic analysis of population structure. Proc. R. Soc. Lond. Ser. B 263: 1619–1626. Benjamini, Y. & Hochberg, Y. 1995. Controlling the false discovery rate – a practical and powerful approach to multiple testing. J. R. Stat. Soc. B Methodol. 57: 289–300. Bierne, N., Welch, J., Loire, E., Bonhomme, F. & David, P. 2011. The coupling hypothesis: why genome scans may fail to map local adaptation genes. Mol. Ecol. 20: 2044–2072. Bierne, N., Gagnaire, P.A. & David, P. 2013. The geography of introgression in a patchy environment and the thorn in the side of ecological speciation. Curr. Zool. 59: 72–86. Bookstein, F.L. 1991. Morphometric Tools for Landmark Data: Geometry and Biology. Cambridge University Press, Cambridge.



Bridle, J.R. & Ritchie, M.G. 2001. Assortative mating and the genic view of speciation - Commentary. J. Evol. Biol. 14: 878–879. Buggs, R.J.A. 2007. Empirical study of hybrid zone movement. Heredity 99: 301–312. Butlin, R.K. 2010. Population genomics and speciation. Genetica 138: 409–418. Butlin, R.K., Ritchie, M.G. & Hewitt, G.M. 1991. Comparisons among morphological characters and between localities in the Chorthippus parallelus hybrid zone (Orthoptera, Acrididae). Phil. Trans. R. Soc. Lond. B 334: 297–308. Butlin, R.K., Saura, M., Charrier, G., Jackson, B., Andre, C., Caballero, A. et al. 2014. Parallel evolution of local adaptation and reproductive isolation in the face of gene flow. Evolution 68: 935–949. Caballero, A., Quesada, H. & Rolan-Alvarez, E. 2008. Impact of amplified fragment length polymorphism size homoplasy on the estimation of population genetic diversity and the detection of selective loci. Genetics 179: 539–554. Conte, G.L., Arnegard, M.E., Peichel, C.L. & Schluter, D. 2012. The probability of genetic parallelism and convergence in natural populations. Proc. R. Soc. Lond. Ser. B 279: 5039–5047. Cruz, R., Vilas, C., Mosquera, J. & Garcia, C. 2004. Relative contribution of dispersal and natural selection to the maintenance of a hybrid zone in Littorina. Evolution 58: 2734– 2746. Endler, J.A. 1977. Geographic variation, speciation, and clines. Monogr. Popul. Biol. 10: 1–246. Erlandsson, J., Rolan-Alvarez, E. & Johannesson, K. 1998. Migratory differences between ecotypes of the snail Littorina saxatilis on Galician rocky shores. Evol. Ecol. 12: 913–924. Faria, R., Renaut, S., Galindo, J., Pinho, C., Melo-Ferreira, J., Melo, M. et al. 2014. Advances in Ecological Speciation: an integrative approach. Mol. Ecol. 23: 513–521. Fisher, R.A. 1950. Gene frequencies in a cline determined by selection and diffusion. Biometrics 6: 353–361. Galindo, J., Moran, P. & Rolan-Alvarez, E. 2009. Comparing geographical genetic differentiation between candidate and noncandidate loci for adaptation strengthens support for parallel ecological divergence in the marine snail Littorina saxatilis. Mol. Ecol. 18: 919–930. Galindo, J., Martinez-Fernandez, M., Rodriguez-Ramilo, S.T. & Rolan-Alvarez, E. 2013. The role of local ecology during hybridization at the initial stages of ecological speciation in a marine snail. J. Evol. Biol. 26: 1472–1487. Grahame, J.W., Wilding, C.S. & Butlin, R.K. 2006. Adaptation to a steep environmental gradient and an associated barrier to gene exchange in Littorina saxatilis. Evolution 60: 268–278. Haldane, J.B.S. 1948. The theory of a cline. J. Genet. 48: 277– 284. Hewitt, G.M. 1988. Hybrid zones – natural laboratories for evolutionary studies. Trends Ecol. Evol. 3: 158–167. Hindorff, L.A., MacArthur, J., Morales, J., Junkins, H.A., Hall, P.N., Klemm, A.K. et al. 2014. A Catalog of Published Genome-Wide Association Studies. www.genome.gov/gwastudies. Hohenlohe, P.A., Phillips, P.C. & Cresko, W.A. 2010. Using Population Genomics to Detect Selection in Natural Populations: Key Concepts and Methodological Considerations. Int. J. Plant Sci. 171: 1059–1071. Hollander, J. & Butlin, R.K. 2010. The adaptive value of phenotypic plasticity in two ecotypes of a marine gastropod. BMC Evol. Biol. 10: 333.

9

Hollander, J., Lindegarth, M. & Johannesson, K. 2005. Local adaptation but not geographic separation promotes assortative mating in a snail – support for ecological speciation. Anim. Behav. 70: 1209–1219. Jackson, B., Kawakami, T., Cooper, S., Galindo, J. & Butlin, R. 2012. A genome scan and linkage disequilibrium analysis among chromosomal races of the australian grasshopper Vandiemenella viatica. PLoS ONE 7: e47549. Janson, K. 1983. Selection and migration in two distinct phenotypes of Littorina saxatilis in Sweden. Oecologia 59: 58–61. Johannesson, K. 2003. Evolution in Littorina: ecology matters. J. Sea Res. 49: 107–117. Johannesson, K. & Tatarenkov, A. 1997. Allozyme variation in a snail (Littorina saxatilis) – deconfounding the effects of microhabitat and gene flow. Evolution 51: 402–409. Johannesson, K., Havenhand, J.N., Jonsson, P.R., Lindegarth, M., Sundin, A. & Hollander, J. 2008. Male discrimination of female mucous trails permits assortative mating in a marine snail species. Evolution 62: 3178–3184. Johannesson, K., Panova, M., Kemppainen, P., Andre, C., Rolan-Alvarez, E. & Butlin, R.K. 2010. Repeated evolution of reproductive isolation in a marine snail: unveiling mechanisms of speciation. Philos. Trans. R. Soc. Lond., Ser. B 365: 1735–1747. Kawakami, T., Butlin, R.K., Adams, M., Paull, D.J. & Cooper, S.J.B. 2009. Genetic analysis of a chromosomal hybrid zone in the Australian morabine grasshoppers (Vandiemenella, viatica species group). Evolution 63: 139–152. Kruuk, L.E.B., Baird, S.J.E., Gale, K.S. & Barton, N.H. 1999. A comparison of multilocus clines maintained by environmental adaptation or by selection against hybrids. Genetics 153: 1959–1971. Li, Y.C., Li, Y., Wu, S., Han, K., Wang, Z.J., Hou, W. et al. 2007. Estimation of multilocus linkage disequilibria in diploid populations with dominant markers. Genetics 176: 1811– 1821. Lindtke, D., Gonzalez-Martinez, S.C., Macaya-Sanz, D. & Lexer, C. 2013. Admixture mapping of quantitative traits in Populus hybrid zones: power and limitations. Heredity 111: 474–485. Luikart, G., Gielly, L., Excoffier, L., Vigne, J.D., Bouvet, J. & Taberlet, P. 2001. Multiple maternal origins and weak phylogeographic structure in domestic goats. Proc. Natl. Acad. Sci. USA 98: 5927–5932. Macholan, M., Munclinger, P., Sugerkova, M., Dufkova, P., Bimova, B., Bozikova, E. et al. 2007. Genetic analysis of autosomal and X-linked markers across a mouse hybrid zone. Evolution 61: 746–771. Malek, T.B., Boughman, J.W., Dworkin, I. & Peichel, C.L. 2012. Admixture mapping of male nuptial colour and body shape in a recently formed hybrid population of threespine stickleback. Mol. Ecol. 21: 5265–5279. Martin, L.B.B., Fei, Z.J., Giovannoni, J.J. & Rose, J.K.C. 2013. Catalyzing plant science research with RNA-seq. Front. Plant Sci. 4: 66. Nosil, P., Funk, D.J. & Ortiz-Barrientos, D. 2009. Divergent selection and heterogeneous genomic divergence. Mol. Ecol. 18: 375–402. Panova, M., Hollander, J. & Johannesson, K. 2006. Site-specific genetic divergence in parallel hybrid zones suggests nonallopatric evolution of reproductive barriers. Mol. Ecol. 15: 4021–4031.


10

J. HOLLANDER ET AL.

Pompanon, F., Bonin, A., Bellemain, E. & Taberlet, P. 2005. Genotyping errors: causes, consequences and solutions. Nat. Rev. Genet 6: 847–859. Porter, A.H., Wenger, R., Geiger, H., Scholl, A. & Shapiro, A.M. 1997. The Pontia daplidice-edusa hybrid zone in northwestern Italy. Evolution 51: 1561–1573. R Core Team 2013. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. Rafauste, N., Orth, A., Belkhir, K., Senet, D., Smadja, C., Baird, S.J.E. et al. 2005. Inferences of selection and migration in the Danish house mouse hybrid zone. Biol. J. Linn. Soc. B 84: 593–616. Rogers, S.M. & Bernatchez, L. 2005. Integrating QTL mapping and genome scans towards the characterization of candidate loci under parallel selection in the lake whitefish (Coregonus clupeaformis). Mol. Ecol. 14: 351–361. Rohlf, F.J. 2003. tpsRelw, Relative Warp Analysis, 1.36. Department of Ecology and Evolution, State University of New York, Stony Brook, NY. Rohlf, F.J. & Marcus, L.F. 1993. A revolution in morphometrics. Trends Ecol. Evol. 8: 129–132. Rolan-Alvarez, E. 2007. Sympatric speciation as a by-product of ecological adaptation in the galician Littorina saxatilis hybrid zone. J. Molluscan Stud. 73: 1–10. Rolan-Alvarez, E., Carballo, M., Galindo, J., Moran, P., Fernandez, B., Caballero, A. et al. 2004. Nonallopatric and parallel origin of local reproductive barriers between two snail ecotypes. Mol. Ecol. 13: 3415–3424. Rol an-Alvarez, E., Bu~ no, I. & Gosalvez, J. 1996. Sex is determined by sex chromosomes in Littorina saxatilis (Olivi) (Gastropoda: Prosobranchia). Hereditas 124: 261–267. Seehausen, O., Butlin, R.K., Keller, I., Wagner, C.E., Boughman, J.W., Hohenlohe, P.A. et al. 2014. Genomics and the origin of species. Nat. Rev. Genet 15: 176–192. Slice, D.E. 1994. DS-DIGIT: Basic Digitizing Software. Department of Ecology and Evolution, State University of New York, Stony Brook, NY. Szymura, J.M. & Barton, N.H. 1991. The genetic-structure of the hybrid zone between the fire-bellied toads Bombina bom-

bina and B. variegata – comparisons between transects and between loci. Evolution 45: 237–261. Vekemans, X., Beauwens, T., Lemaire, M. & Roldan-Ruiz, I. 2002. Data from amplified fragment length polymorphism (AFLP) markers show indication of size homoplasy and of a relationship between degree of homoplasy and fragment size. Mol. Ecol. 11: 139–151. Via, S. & West, J. 2008. The genetic mosaic suggests a new role for hitchhiking in ecological speciation. Mol. Ecol. 17: 4334–4345. Webster, S. 2014. Selection for local adaptation in Littorina saxatilis. PhD Thesis, University of Sheffield, UK. Wilding, C.S., Butlin, R.K. & Grahame, J. 2001. Differential gene exchange between parapatric morphs of Littorina saxatilis detected using AFLP markers. J. Evol. Biol. 14: 611–619. Wood, H.M., Grahame, J.W., Humphray, S., Rogers, J. & Butlin, R.K. 2008. Sequence differentiation in regions identified by a genome scan for local adaptation. Mol. Ecol. 17: 3123–3135.

Supporting information Additional Supporting Information may be found in the online version of this article: Figure S1 Correlation of FST values between Rams€ o and Ramsholmen calculated with DFDIST for each AFLP locus (1167 loci). Figure S2 The relationship between linkage disequilibrium and the product of allele frequency slopes for the two transects: (A) Rams€ o, (B) Ramsholmen. Table S1 Parameters of best-fitting unimodal cline models for the hybrid index. Data deposited at Dryad: doi:10.5061/dryad.2jh32 Received 2 May 2014; revised 24 November 2014; accepted 25 November 2014


Detection of outlier loci and their utility for fisheries management.

Outlier Loci and Selection Signatures of Simple Sequence Repeats (SSRs) in Flax (Linum usitatissimum L.).

Lipids in digestive gland of Littorina saxatilis rudis (Maton) and in daughter sporocysts of Microphallus similis (Jäg. 1900).

Positive Selection on Loci Associated with Drug and Alcohol Dependence.

Historical Biogeography of the Marine Snail Littorina saxatilis Inferred from Haplotype and Shell Morphology Evolution in NW Spain.

Genome Wide Association Study Identifies New Loci Associated with Undesired Coat Color Phenotypes in Saanen Goats.

Natural and sexual selection on many loci.

Classification of High Intensity Zones of the Lumbar Spine and Their Association with Other Spinal MRI Phenotypes: The Wakayama Spine Study.

Validation of outlier loci through replication in independent data sets: a test on Arabis alpina.

Numerical studies on two-loci selection models with general viabilities.

Adaptive Laboratory Evolution of Antibiotic Resistance Using Different Selection Regimes Lead to Similar Phenotypes and Genotypes.

An Adaptive Association Test for Multiple Phenotypes with GWAS Summary Statistics.

Evidence of selection on splicing-associated loci in human populations and relevance to disease loci mapping.

Ontogenetic selection on hatchery salmon in the wild: natural selection on artificial phenotypes.

The Spatial Mixing of Genomes in Secondary Contact Zones.

Migration-selection balance at multiple Loci and selection on dominance and recombination.

Intense natural selection preceded the invasion of new adaptive zones during the radiation of New World leaf-nosed bats.

Tricarboxylic acid cycle enzymes in the digestive gland of Littorina saxatilis rudis (Maton) and in the daughter sporocysts of Microphallus similis (Jäg.) (Digenea: Microphallidae).

Aerobic glucose metabolism in the digestive gland of Littorina saxatilis rudis (Maton) and in the daughter sporocysts of Microphallus similis (Jäg.).

EigenGWAS: finding loci under selection through genome-wide association studies of eigenvectors in structured populations.

Outlier Loci Detect Intraspecific Biodiversity amongst Spring and Autumn Spawning Herring across Local Scales.

Association of PRPS1 Mutations with Disease Phenotypes.

Variable patterns of obesity and cardiometabolic phenotypes and their association with lifestyle factors in the [email protected] study.

Selection towards different adaptive optima drove the early diversification of locomotor phenotypes in the radiation of Neotropical geophagine cichlids.