Neutral and adaptive genomic signatures of rapid poleward range expansion
J. Swaegers1, J. Mergeay1,2,A. Van Geystelen3,4, L. Therry1,5, M.H.D. Larmuseau3,4,6, R. Stoks1
1Laboratory of Aquatic Ecology, Evolution and Conservation, University of Leuven, Deberiotstraat 32, 3000 Leuven, Belgium
2Research Institute for Nature and Forest, Gaverstraat 4, B-9500 Geraardsbergen, Belgium
3Laboratory of Socioecology and Social Evolution, Department of Biology,University of Leuven,Naamsestraat 59, 3000 Leuven, Belgium
4Laboratory of Forensic Genetics and Molecular Archaeology, University of Leuven, Kapucijnenvoer 33, 3000 Leuven, Belgium
5Station d’Ecologie Expérimentale du CNRS à Moulis, USR 2936, 09200 Moulis, France.
6Department of Genetics, University of Leicester, University Road, LE1 7RH, Leicester, UK
Keywords: Collateral evolution, countergradient variation, genome scan, Odonata, polygenic multilocus outlier detection method, range expansion
Corresponding author: JanneSwaegers, University of Leuven,Deberiotstraat 32, 3000 Leuven, Belgium,+3216324575,
Running title:Signature of evolution during range expansion
Abstract
Many species are expanding their range polewards and this has been associated with rapid phenotypic change.Yet, it is unclear to what extent this reflects rapid genetic adaptation or neutral processes associated with range expansion, or selection linked to the new thermal conditions encountered.To disentangle these alternatives,we studied the genomic signature of range expansionin the damselflyCoenagrionscitulumusing 4950 newly developed genomic SNPs and linked this to the rapidly evolved phenotypic differences between core and (newly established) edge populations. Most edge populations were genetically clearly differentiated from the core populations and all were differentiated from each other indicating independent range expansion events.In addition,evidence forgenetic drift in the edge populations, and strong evidence for adaptive genetic variation in association with the range expansion was detected. Weidentifiedone SNP under consistent selection in four of the five edge populations and showed that the allele increasing in frequency is associated with increased flight performance. This indicates collateral, non-neutral evolutionary changes in independent edge populations driven by the range expansion process. Wealso detected a genomic signature of adaptation to the newly encountered thermal regimes, reflecting a pattern of countergradient variation. The latter signature wasidentified ata single SNP as well asin a set of covarying SNPs using a polygenic multilocus approach to detect selection.Overall, this studyhighlights how a strategic geographic sampling design and the integration of genomic, phenotypic and environmental data can identify anddisentangle the neutral and adaptive processes that are simultaneously operating during range expansions.
Introduction
Global warming is imposing strong selection, leading to rapid changes in traits such as phenology and thermal resistance (Waltheret al.2002; Hoffman & Sgro 2011; Stokset al.2014). Yet, it is largelyunknownwhether these changes are genetically based and therefore reflect evolution, or are driven by phenotypic plasticity.Plasticity is increasingly reported as a cause of phenotypic responses to global warming, while solid evidence forunderlying genetic change is still scarce (Merilä 2012; Merilä & Hendry 2014).While most studies have focused mainly onin situresponses to global warming, many species are showing rapid poleward range expansions (Hicklinget al.2006; Chenet al.2011).These poleward range expansions are also associated with rapid trait changes, including increased dispersal ability (Hill et al. 2011) and tolerance for newly encountered thermal regimes (Stoks et al. 2014), which are potentially plastic(Hassallet al.2014). Genomic approaches provide a powerful, yet largely unexplored way to identify genetic climate-driven trait changes (Reusch & Wood 2008; Franks & Hoffmann 2012). The few studies on genomic signatures of climate-driven range expansion did not link genetic changes to phenotypes and did not unravel the evolutionary processes involved (Merilä & Hendry 2014).
A key challenge when studyingevolutionary change during range expansionsis to disentangle neutral from non-neutral processes(Klopfsteinet al.2006), which both may be strong and underlie trait evolution (Leinonenet al.2013). Neutral founder effects can randomly change allele frequencies (Excoffier & Ray 2008; Hallatscheket al.2007) and low frequency alleles can thereby ‘surf’ on the wave of expansion leading to high frequencies at the range front (Klopfsteinet al.2006; Hallatschek & Nelson 2008). Empirical studies have indeed documented a reduction in genetic diversity driven by founder effects at expanding range fronts(Wattset al.2010; Garrowayet al.2011; Swaegerset al.2013). Allelic surfing can thereby generate false positive signatures of selection (Excoffier & Ray 2008; Excoffieret al.2009).However, the direction of allele frequency change driven by allelic surfing is expectedto be random across loci and independent among different axes of expansion (Excoffier & Ray 2008; Françoiset al.2010). Hence, when allele frequency changes occur at a locus in the same direction across independently established edge populations, this supports a role of collateral evolution (convergent evolution of a locus that was polymorphic in the shared ancestral population; Stern 2013), rather than genetic drift (Hohenloheet al.2011; Whiteet al.2013).
It is expectedthat the non-neutral range expansion processgenerates strong selective pressures through (i) spatial sorting (Shineet al.2011)where individuals that establish new edge populations aredominated by those with high dispersal ability and/or (ii) natural selection whenevolution occurs due to fitness benefits of being amongst the earliest colonists(Travis & Dytham 2002). Therefore, the range expansion process is predicted to result in rapid evolution of a broad set of traits, including aspects of life history, morphology, physiology and behaviour,which are selected toward values that increase the rate ofspread (Phillipset al.2010; Shineet al.2011). Accordingly, empirical studies have documented thatpopulations at the edge of a moving frontrapidly evolved a faster life history (Phillips 2009;Therryet al.2014b), morphological changes related to ahigher dispersal capacity (Hillet al.2011; Therryet al.2014a,c),and adjustments in their immune function(Therryet al.2014a; Brownet al.2015).It remains, however, unknown howrapid trait changes during range expansions are mediated byallele frequency changes.To our knowledge, only Buckleyet al.(2012) have investigated the signatures of such evolutionary changes at the genomic level in the context of a natural range expansion. They reported in thebutterflyAriciaagestisthata portion (5%) of the genotyped markers showed signatures of selection associated with range expansion,although the identity of these genes remained unknown and could not be matched to phenotypic evolution.
Coupling genomic data with phenotypic data could partly circumvent this problem (Stinchcombe & Hoekstra 2008; LepaisBacles 2014), and at the same time provide evidence for the genetic basis of rapid phenotypic change. However, this strategy has, to our knowledge, not yet been undertaken in the context of natural range expansions(see also Merilä & Hendry 2014). This strategy may be challenging as the above-mentionedphenotypic traits commonly associated with range expansion are typically polygenic (Roff & Fairbairn 2001). In these cases phenotypic change is expected to occur through subtle shifts in allele frequencies of many loci (Hancocket al.2010; Pritchard & Di Rienzo 2010) making evolutionary changes more difficult to detect, especially when genetic drift canmask these changes. Moreover,convergent phenotypic changes (such as those expected at the expansion front along independent axes) can be mediated by different molecular-genetic pathways (Kauttetal. 2012; Rodaetal. 2013; Soria-Carrascoetal. 2014). Alsowithin the more general context of invasions, few studieshave provided molecular evidence of adaptation(Chownet al.2014). A notable exception is the study by Whiteet al.(2013) reporting in the invasive bank vole (Myodesglareolus) recurrent directional changes in allele frequencies of loci related to immunity and behaviour along replicate expansion axes.
Besides therange expansionprocess, poleward moving organisms also experience strong selection imposed by thenovel thermal regimes they encounter (BlanckenhornDemont, 2004; Shamaet al.2011; Therryet al.2014a-d)and disentangling both non-neutral processes may be challenging (Frichotet al.2015).Countergradient variation, whereby genetic influences buffer environmental influences (Conoveret al.2009), may thereby generate similar trait changes as the expansion process. For example, flight ability may both increase in response to selection driven by the process of range expansion (Travis & Dytham 2002; Shineet al.2011) as well as by antagonistic selection to counter potential reduction in flight ability at the suboptimal thermal regimes at higher latitudes (Therryet al.2014c).
Within 20 years,the poleward moving damselflyCoenagrionscitulumhas established populations more than 200 km north-eastward from its historical range in France along a broad expansion front in Belgium, Germany, The Netherlands and Great Britain (Swaegerset al.2013).Due to the overlap in thermal regimes between core and edge populations (Table 1) the study system allows to disentangle the effects of thermal regimes and the range expansion process.The range expansion has been associated with the expected increase in several flight-related traits (Therryet al.2014a-d). Specifically, increased relative investment in flight muscle mass has been detected in field-collected edge animals (Therryet al.2014c), as well as in edge animals reared in a common garden experiment (Therryet al.2014a), suggesting a genetic basis. The finding that flight muscle mass positively covaried with flight endurance in this species (Therryet al.2014d) indicates evolution of higher flight ability during range expansion. Also the newly encountered thermal regimes have been associated with phenotypic changes inC. scitulum.Body mass and flight muscle ratio increased with the number of degree days available for larval development while, consistent with a scenario of countergradient variation, body mass and fat content decreased with the mean air temperature during the adult flight period (Therryet al.2014a,d). During the range expansion, Swaegerset al.(2013) discovered mild founder effects on the basis of genetic structure at 12 microsatellite markers and found that newly established populations along the broad expansion front are genetically more differentiated from each other than populations from the core of the species’ distribution, indicative of independent and replicated colonization and expansion events. This allows distinguishing non-neutral signatures of convergent genetic shifts from genetic drift along different axesof expansion.
Here, we aim to (i) identify the genomic signature of the range expansion processand (ii) differentiate it from signatures generated by neutral processes and non-neutral thermally-induced selection and (iii) link it to the well-documented rapid phenotypic changes inCoenagrionscitulum. To this end we applied both a single-locus and multilocus polygenic outlier approach (Bourretat al.2014) based on a genome scan of ca. 5000 SNPs on individualssampled across independent range expansions(Swaegerset al.2013) where thermal regimes were partly decoupled from the spatial edge-core axis (Therryet al.214a; Frichotet al.2015).
Methods
Study system
Males ofC. scitulumwere collected from five core populations in France and five edge populations in Germany and the Netherlands during the single flight period of the species in June-July 2010-2011 (Fig. 1, Table 1, see also Therryet al.2014d). We only studied males because no differences in dispersal rates are apparent between sexes inCoenagrionspecies (Conradet al.2002),hence no different evolutionary responses in flight ability during range expansion are to be expected between sexes. We defined core populations as those populations situated within the historical continuous distribution of the species (based on Dommanget (1994), Dijkstra & Lewington (2006) and Boudot (2013, unpublished data). Edge populations were defined as populations at the expanding range front and sampled less than five years after founding.Note that as we were interested in the signal of the range expansion process, stable edge populations are not considered as edge populations. The locations of the study populations as well as the entire species distribution are shown in Fig. 1. A phylogeographic study ofC. scitulumrevealed that the origin of the studied edge populations lies in Western Europe (Swaegerset al.2014).To assess the influence of latitude-associated thermal regimes in shaping genomic signals we will explicitly consider the two thermal regimes that influence larval and adult phenotypes in the study species. Damselflies have a complex life cycle with an aquatic larval stage where growth occurs and a terrestrial adult stage where dispersal occurs (Stoks & Cordoba-Aguilar 2012). To represent the larval thermal regimes we use the number of larval degree days available for larval growth(the modelled cumulative number of degrees with water temperatures above 8°C) and to represent the adult thermal regimes we use the mean air temperature during theadult flight period (for details of the calculation, see appendix C in Therryet al.2014d).The variance in thermal regime variables explained by population status (edge/core) was low (larval degree days:R2=19%; adult flight period temperature:R2=0.9%) allowing sufficient separation of the explanatory variables.Each male was phenotyped for ecologically relevant traits related to range expansion (flight endurance, flight muscle ratio, body mass and fat content) and the resulting patterns were described in Therryet al.(2014d).
DNA extractionandgenotyping
DNA from 190 individuals (15-20 males per population, Table 1) was extracted using DNeasy Extraction kit (Qiagen Inc., Valencia, CA, USA) with a final concentration of on average 102 ± 54 (standard deviation) ng/µl measured using a Nanodrop 2000 spectrometer (Isogen Life Science, Belgium). The integrity of the yielded DNA was assessed on 1.5% agarose gels. Extracted DNA was sent to the Cornell Institute for Genomic Diversity to conduct Genotyping-by-sequencing of reduced representation libraries cut with the enzymePstI (CTGCAG) (Elshireet al.2011; Narumet al.2013). This technique constructs reduced representation libraries for the Illumina sequencing platform. Similar to RAD sequencing (Hohenloheet al.2010), DNA is digested using the restriction enzymes. We chose the enzymePstI as it does not cut frequently so that higher coverage is gained and hence greater reliability when genotypes are called. Subsequently, the fragmented DNA was ligated to a barcoded adaptor and a second adaptor with appropriate sticky ends. Apart from DNA samples, also two controls which did not contain DNA were randomly added to the plates. Genomic fragments were amplified using long primers that matched the adaptors. Next, the constructed library was diluted and sequenced on an Illumina HiSeq 2000. For more details on the procedure we refer to Elshireet al.(2011).
The obtained raw sequences were analysed with the Universal Enabled Analysis Kit (UNEAK) pipeline, implemented inTasselv3.0 (Bradburyet al.2007).Firstly, the reads were all trimmed to 64bp to remove the barcodes. Identical reads were clustered into tags and pairwise alignment was subsequently performed on the tag. Tag-pairs with a match of one base pair were considered as candidate SNPs. With an error tolerance rate of 0.03, only reciprocal pairs of tags were retained as SNPs (Luet al.2013).Genotypes were filtered using default parameters inTasselv3.0. Only SNPs with less than 20% missing data and individuals with more than 90% genotyped SNPs were retained.
Geographic patterns in genetic diversity and structure
Using the large set of identified genomic SNPs we assessed the effect of range expansion on genetic diversity and structure. We therefore calculated the mean expected heterozygosity (He) and mean number of polymorphic sites per population and pairwise genetic differentiation (usingFST) inArlequin 3.5(ExcoffierLischer 2010).Patterns of genetic structure were further explored using Principal Components Analysis (PCA) with the adegenet package (Jombart 2008, function dudi.pca) in R (R Development Core Team, 2014).In order to compare our results with those obtained by Swaegerset al.(2013) on microsatellites andto account for population structure in subsequent analyses,we also analysed the genetic structure usingStructure (Pritchardet al.2000).Correlated allele frequencies were assumed and the admixture model simulations were run with prior information on the site. Runs were performed using a burn-in period of 20,000 replicates and a sampling period of 100,000 replicates. Simulations were run for a number of clusters (K) ranging from one to ten, and 20 iterations were performed for each K. The lnlikelihoods from these iterations were averaged for each K and the most likely K was selected with the delta K method(Evannoet al.2005)usingStructureHarvester(Earl & von Holdt 2011).
Single-locus outlier detection
We screened for putative adaptive loci associated with range expansion that are shared between the different edge populations, thereby disentangling genetic drift from collateral evolutionary change.We searched for signatures of selection associated with range expansion using a parallel outlier detection approach (Hohenloheet al.2010). For every edge population we assessed the number of significant outlier loci using latent fixed mixed modelling (LFMM, Frichotet al.2013) where we analysed the five core populations together with every edge population separately with the SNP matrix as dependent variable.FST-outliers repeatedly detected across separate pairwise comparisons can provide stronger evidence for the selective significance of that locus opposed to a signature of genetic drift (Hohenloheet al.2011). The status of the population (edge or core) was included in the model as a dummy variable and we assumed one or two latent factors matching the groups identified by theStructureanalysis (for the analysis involving population DSS we set one latent factor and for the remaining analyses two latent factors). Compared to other methods, LFMM reduces the type I error in situations with complex underlying genetic structure (Frichotet al.2013; de Villemereuilet al.2014; Firchotet al.2015). We used the R package qvalue v1.43 (Dabneyet al.2011) to correct for multiple testing with a false-discovery rate of 0.05. Loci shared among the analyses per edge population that had frequency shifts in the same direction were assumed to be potential candidate SNPs of selection during range expansion.
To assess genomic changeslinked with countergradient variation, we specifically evaluate whether latitude-associated larval and thermal regimes also shaped allele frequencies, and whether this contributed to the outlier loci associated with the range expansion that were detected by the parallel outlier detection approach. Here, we used LFMM using the global dataset with either larval degree days or adult flight period temperature (Therryet al.2014d, Table 1) as environmental independent variables.
Polygenic multilocus outlier detection
In addition to the outlier detection method, which focuses on single locus detection, we also applied a polygenic multilocus approach (Bourretet al., 2014).Within this multivariate genetic framework markers are considered as multilocus entities that influence phenotypic traits in a polygenic manner.This method reduces the number of individual markers to a limited number of PC axes representing sets of independently, covarying groups of markers on the basis of their genotypic distribution. In a first step we extracted PC axes and performed a varimax rotation, the number of extracted PC axes thereby equals the number of markers. In a second step, we only retained the significant PC axes.PC axes were considered significant when their associated eigenvalue was larger than the 99th value of the distribution of eigenvalues derived from 1000 random data sets (this isthe parallel analysis criterion of Horn 1965).For this reduced set of significant PC axes; we determined for each PC axis which SNP markers were significantly correlated with it using a significance level of 0.001. Given the number of loci and individuals, this corresponds to an absolute loading weight (correlation between a given SNP marker and PC axis) higher than 0.24 on a given PC axis.