A replicated climate change field experiment reveals rapid evolutionary response in an ecologically important soil invertebrate.

Whether species can respond evolutionarily to current climate change is crucial for the persistence of many species. Yet, very few studies have examined genetic responses to climate change in manipulated experiments carried out in natural field conditions. We examined the evolutionary response to climate change in a common annelid worm using a controlled replicated experiment where climatic conditions were manipulated in a natural setting. Analyzing the transcribed genome of 15 local populations, we found that about 12% of the genetic polymorphisms exhibit differences in allele frequencies associated to changes in soil temperature and soil moisture. This shows an evolutionary response to realistic climate change happening over short-time scale, and calls for incorporating evolution into models predicting future response of species to climate change. It also shows that designed climate change experiments coupled with genome sequencing offer great potential to test for the occurrence (or lack) of an evolutionary response.

Introduction

It is well‐documented that species have adapted evolutionarily to climatic changes over geological time scales (Hewitt, 2000; Davis & Shaw, 2001). However, current climate change occurs at a much faster rate than in the past as inferred from Quaternary climate records (IPCC, 2014) and this is likely to impact the future of numerous species (Walther et al., 2002; Parmesan & Yohe, 2003; Deutsch et al., 2008). Possible biotic responses to altered environmental conditions include migration, phenotypic plasticity, and genetic adaptations in response to the natural selection pressure exerted by the altering conditions (Schilthuizen & Kellermann, 2014). Phenotypic responses to climate change are known from species in both terrestrial and aquatic systems (Portner & Knust, 2007; Deutsch et al., 2008), but there are obvious limits to how much a genotype can respond plastically (Overgaard et al., 2014). Genetic adaptations are thus crucial for the long‐term persistence especially for species of low vagility (Franks & Hoffmann, 2012).

Most of the studies published so far use natural climatic (latitudinal or altitudinal) gradients (Bradshaw & Holzapfel, 2001; Balanyá et al., 2006), or the evolution of invasive species within a novel environment (Moran & Alexander, 2014) as proxy for predicted climate change rather than controlled replicated experiments (Merilä, 2012). However, very few studies have documented genetic responses to climate change in manipulated experiments in the natural environment of an organism.

Field‐scale climate experiments, in which the abiotic environment is manipulated to mimic predicted future climate, are useful tools to study the effects of changes in abiotic factors, such as temperature, precipitation, and CO2 levels (Mikkelsen et al., 2008). However, to our knowledge, these powerful randomized and replicated designs have rarely been used for studying either evolutionary potential or actual evolutionary genetic responses to changes in climatic conditions (Lau et al., 2014). Here, we show that these designs coupled with genome surveys at the population level offer a great potential to test for the occurrence (or lack) of an evolutionary response. We illustrate this by applying this approach to a non‐model soil invertebrate that is ecologically important for terrestrial ecosystems.

The terrestrial annelid genus, Chamaedrilus [previously referred to as a single species named Cognettia sphagnetorum Vedjovsky, 1878 (Clitellata; Enchytraeidae) (Martinsson et al., 2015)] is widespread in northern and western parts of Europe (Briones et al., 2007). Chamaedrilus chlorophilus Friend 1913 (Martinsson et al., 2015) is abundant in dry heathland of Denmark and other parts of continental Europe (Schmelz & Collado, 2010; Holmstrup et al., 2012). This species reproduces primarily asexually by fragmentation probably with 2–4 generations per year, depending on the environmental conditions (Augustsson & Rundgren, 1998; Schmelz & Collado, 2010). It reaches a size of 10–20 mm and is confined to the upper soil layers that are rich in dead organic matter. Chamaedrilus spp. have very limited dispersal abilities (Sjögren et al., 1995) and therefore is ideal for studies of genetic changes in open field experiments. Chamaedrilus spp. forms a substantial fraction of the soil animal biomass in temperate organic rich soils and plays an important role in the decomposition process and nutrient mineralization of these soils by consuming microorganisms and dead organic matter (Williams & Griffiths, 1989; Maraldo & Holmstrup, 2010; Larsen et al., 2011). Enchytraeids, and in particular C. chlorophilus, are sensitive to drought (Maraldo et al., 2009, 2010) and increased soil temperatures are also expected to influence abundance of enchytraeids (Briones et al., 1997, 2007). Previous field experiments have shown that increased intensity of summer drought causes transient bottlenecks in populations of C. chlorophilus followed by recovery to normal population densities within months (Maraldo & Holmstrup, 2010). Such recurring drought periods could foster genetic adaptation to occur due to natural selection processes. However, no studies have been conducted to test if observable ecological responses – in population dynamics or change of species abundance – are also accompanied by evolutionary changes in these populations.

Next‐generation sequencing technologies (NGS) offer the opportunity to investigate within‐species genetic variation in absence of prior genomic resources. Here, we survey genetic variation in the expressed fraction of the genome of C. chlorophilus using transcriptome sequencing of samples from an experimental design where realistic future climate conditions were enforced and replicated in a natural setting since 2005 (Mikkelsen et al., 2008). We investigate whether the manipulated environmental changes left a detectable evolutionary signature in these replicate populations. We document rapid adaptation to the climatic environment manipulations. Our findings pave the way for using molecular tools to detect adaptation in a variety of non‐model organisms and will impact on our understanding of the adaptation of organisms to climate change.

Materials and methods

Experimental design and sampling of individuals

The experimental site is situated at Brandbjerg, Denmark (55°53′N 11°58′E) on a nutrient‐poor sandy deposit with a heath/grassland ecosystem dominated by a grass (Deschampsia flexuosa (L.) Trin) and an evergreen dwarf shrub (Calluna vulgaris (L.) Hull). The mean annual air temperature is 8.0 °C and the mean annual precipitation is 613 mm. The experiment was initiated in October 2005 and included single and all possible combinations of elevated CO2 concentration, increased temperature (T), drought in late spring/early summer (D), and untreated controls for reference (A, for ambient).

The drought treatment was applied with waterproof curtains automatically pulled over the vegetation during rain events to prevent precipitation to reach the ground and retracted when the rain stopped. Since drought was applied by reducing precipitation, the resulting soil water content in drought plots was variable from year to year. The CO2 was distributed by a Free Air Carbon Enrichment (FACE) system with a target concentration of 510 ppm. CO2 fumigation started 30 min after sunrise and ended 30 min before sunset all year round, except during periods with full snow cover of the vegetation. The temperature enhancement was achieved by ‘passive night‐time warming’, where a light scaffolding (0.5 m height) carrying a curtain reflects the outgoing infrared radiation. The curtains were automatically pulled over the vegetation at sunset and retracted at sunrise. In case of rain or heavy winds (>7 m s−1) during the night the curtains were automatically retreated to avoid hydrological disturbance or damage to the curtains (for further details, see Mikkelsen et al., 2008).

The eight treatments were placed in pairwise octagons of 6.8 m across receiving ambient and elevated CO2 respectively. Each octagon was divided into four ‘slices’ (9.1 m2 per plot) to provide all eight treatment combinations (Fig. 1a). Each combination was replicated six times (total of 48 ‘slices’, hereafter termed as ‘plots’) with six octagons at ambient CO2 and six receiving elevated CO2 (Fig. 1b). The distances between the octagons were at least 17 m, to prevent contamination of ambient CO2 plots by the adjacent FACE system. Soil temperature (5 cm depth) and soil moisture (TDR measurements; 0–20 cm depth) were logged on an hourly basis throughout the experimental period (Mikkelsen et al., 2008). Based on these data two covariates were calculated for each plot: the average temperature at 5 cm depth (Ts5) and the average soil water content (0–20 cm depth) during the last week of drought periods when soil moisture was at its lowest (SWCdrought). The soil temperature covariate was calculated as the average of hourly measurements of Ts5 over the 7‐year experiment. The soil moisture covariate was calculated as the average SWCdrought across the 7 years.

Figure 1

(a) Overview of the experimental manipulation within one octagon (6.8 m width) comprising four slices (hereafter plots). Note that half of the octagons where also submitted to an elevated CO 2 treatment (not shown on the figure). A denotes the ‘Ambient’ (control) condition, while T and D denote, respectively, plots where temperature and drought were manipulated separately (T, D) or jointly (TD). (b) Overview of the design and spatial location of samples listed in Table 1. Octagons 2, 4, 5, 8, 10, and 12 received elevated (CO 2) and are circled in bold. Solid lines show the boardwalks and the locations of the two meteorological stations are marked as M1 and M2. The two rectangles represent buildings that house computers, control systems and field laboratories.

The increased temperature treatment resulted in 25–75 additional degree‐days annually in the soil (5 cm depth) and somewhat more in air and on the soil surface. Warmed plots had the same temperature as ambient plots until the treatments were started in October 2005 (Fig. S1).

During the 7‐year field experiment, abundance of C. chlorophilus was repeatedly reduced in drought treated plots whereas experimental warming did not have detectable effect (Maraldo et al., 2008; Andresen et al., 2011; Larsen et al., 2011; Holmstrup et al., 2015). Thus, low soil water content reduced the local population density, but populations recovered before the drought manipulation of the following year.

Plots for our study were sampled in November 2012, i.e. after 7 years of treatments. Among the set of plots available we sampled 15 plots each representing five independent replicates of three environmental conditions: ambient (no manipulation), drought, or a combination of increased temperature, drought, and increased CO2. Hereafter, we consistently refer to these plots as ‘populations’ (Table 1). Within each of these 15 populations, two soil samples were collected. Samples were collected using a soil corer (inner diameter 5.5 cm, depth 9 cm), divided into three layers of 3 cm each and kept in a cold room at 5 degrees C for 5 days. Enchytraeids were subsequently extracted using Baerman technique (for further details, see Maraldo et al., 2008). The enchytraeids of these two samples were pooled and 30 C. chlorophilus individuals were morphologically identified using a dissecting microscope (Olympus) at 60× magnification and pooled for RNA extraction and sequencing. Chamaedrilus chlorophilus was the only representative of this genus (Holmstrup et al., 2015) and was therefore easily distinguishable from other species present at the study site. Thus, it was not necessary to mount specimens on slides before sampling (species identity and homogeneity was confirmed a posteriori by barcoding, see below).

Table 1

List of populations, spatial location on the design and experimental treatment associated with each sampled population

Population id	Octogon	Treatment
F1‐P1	1	Ambient
F1‐P2	1	Drought
F3‐P1	3	Drought
F3‐P4	3	Ambient
F4‐P2	4	TDCO2
F5‐P3	5	TDCO2
F6‐P2	6	Ambient
F7‐P2	7	Drought
F8‐P3	8	TDCO2
F9‐P1	9	Ambient
F9‐P4	9	Drought
F10‐P2	10	TDCO2
F11‐P1	11	Drought
F11‐P4	11	Ambient
F12‐P3	12	TDCO2

Octogon number refers to spatial position depicted on Fig. 1.

RNA extraction, cDNA library construction, and sequencing

For each sample, RNA was extracted from the pooled 30 individuals following (Gayral et al., 2011). This was followed by quantification by the Invitrogen (Carlsbad, CA, USA) Quant‐iT™ RiboGreen® RNA Reagent based on the manufacturer's protocol and verification of RNA quality by Agilent (Santa Clara, CA, USA) 2100 Bioanalyzer profile. Synthesis of cDNA and construction of libraries were done with the Illumina, Inc. (San Diego, CA, USA) TruSeq RNA Sample Preparation v2 Kit. Fragment size of selected cDNA was between 300 and 600 pb. Libraries were indexed and sequenced using two flow cells of a Illumina, Inc. HiSeq 2000 sequencer with the 2 × 100 cycles, paired‐end, indexed protocol provided by the Montpellier Genomix Platform (http://www.mgx.cnrs.fr).

Bioinformatic pipeline for contig assembly, read mapping, and polymorphism detection

After filtering for read quality, we used a bioinformatics pipeline developed specifically for de novo assembly of transcriptome data from non‐model organisms, read mapping and Single Nucleotide Polymorphism (SNP) detection. Reads were assembled to contigs using a combination of the programs Abyss v1.1. (Simpson et al., 2009) and Cap3 (Huang & Madan, 1999), following the protocol described in (Cahais et al., 2012). Then reads were mapped to contigs using the BWA program (Li & Durbin, 2009) with default options. A position was called as polymorphic (SNP) when the frequency of the second most frequent base was above 0.05 across all samples. This threshold implies that each detected variant was carried by at least ~15 distinct reads, thus conservatively avoiding any impact of sequencing errors.

Control for contamination of the samples

We first checked that the samples in our transcriptome data comprised a single species. This is imperative to ensure that we are tracking genetic changes within species and not merely changes in abundance of morphologically identical but genetically distinct entities.

To do so we used highly expressed contigs of our assembly that are homologs of well‐established barcoding/phylogenetic markers for animals. We used published sequences of Cox‐1 and 18S ribosomal genes from various annelids, including sequences deposited as C. sphagnetorum (and recently re assigned to Chamaedrilus chlorophilus), as well as unpublished Cox‐1 sequences kindly made available by Rüdiger Schmelz (ECT Oekotoxikologie GmbH, Flörsheim am Main, Germany). We aligned these to all contigs displaying strong homology to Cox‐1 as a check for potential contamination by other Chamaedrilus species that had been wrongly identified at the pooling stage.

Accuracy of the pool based estimates of SNP frequency: simulated pooling from datasets with individual SNP calling

We validated the SNP calling approach on pooled samples by testing the approach on independent datasets in which individuals were separately sequenced and genotyped. The three datasets examined here (marine tunicate Ciona intestinalis A, butterfly Melitaea cinxia, and freshwater snail Physa acuta) were obtained using the same RNA‐seq protocols and the same bioinformatics pipeline for transcriptome assembly and read mapping as in the current study. Each dataset consisted of 10 individuals.

In each of the three species, we compared the ‘pool’ (individual labels ignored) and ‘separate’ (individual label aware) methods of SNP calling and allele frequency estimation. The ‘pool’ approach was identical to that followed in this study. The ‘separate’ approach was that introduced by (Gayral et al., 2011) and applied by (Romiguier et al., 2014). The two approaches were convergent, with the same two alleles being called in >99% of cases, and allele frequencies being well‐correlated (r 2>0.6). This was so in spite of a rather imbalanced contribution of the distinct individuals. For instance, in C. intestinalis, in 36% of the SNPs, the among‐individuals minimal number of reads was below one‐tenth of the median, and in 16% of the SNPs the maximal number of reads was above 10 times the median. Then we defined two arbitrary groups of five individuals each and calculated the difference in allele frequency between the two groups, Δall_freq again comparing the two approaches (Table S1). A high correlation was obtained between Δall_freq[separate] and Δall_freq[pool] (r 2 ~ 0.8, Fig. S2). The variance across SNPs of Δall_freq was actually found to be higher with separate than pooled individuals, suggesting that the pooled approach is conservative and unlikely to yield an excess of spurious SNPs. This control demonstrates a good agreement between the pooling and labeling strategies, and suggests that pooling is conservative regarding the detection of allele frequency variation between samples even at relatively low per individual coverage, consistent with the conclusions of Gautier et al. (2013) and Schlötterer et al. (2014).

Statistical analysis of SNP frequency data

Counts for each allele at each SNP position were modeled using a logistic regression framework assuming that counts of reads per alleles are binomially distributed with a mean that can be a function of a continuous covariate (either the temperature or the soil moisture covariate). This accommodates naturally the unequal coverage per SNP position in each pool. Likelihood ratio tests comparing models with and without the continuous covariate provides formal tests for the association of each SNP with the underlying environmental gradient. We inspected visually the empirical distribution of P‐values. In particular, we checked visually that the empirical distribution of P‐values in the range (0.5–1) was uniform. This attests that the test we used is well‐behaved (not too conservative or liberal) and we subsequently used the false discovery rate estimator of (Storey & Tibshirani, 2003) to estimate the proportion of false positives (q‐value) given an individual P‐value threshold.

We assessed the robustness of the logistic regression framework used to analyze differences in SNPs counts in each sample. Our framework assumes that counts of reads per allele at each SNP are binomially distributed. This is only strictly true if each individual in the pool contributes exactly the same number of reads obtained through sequencing. Even in the absence of any allele‐specific bias in gene expression, we expect the number of reads contributed to a pool by each individual to be at best Poisson distributed with identical means. This will yield estimates of SNP frequency within pools with inflated variance relative to the nominal variance expected if counts were strictly binomial.

In order to further check that our likelihood ratio test for the effect of an environmental covariate is statistically sound, we used a permutation‐based approach to obtain a P‐value for each test. Such permutation‐based test also guards against spurious association merely created by experimental plot closer to each other having also more similar covariates and more similar SNPs. As this procedure is computer intensive, we re‐analyzed a random subset of 500 SNPs in our dataset. For each SNP, the environmental covariate values were permuted 10 000 times and, for each permuted dataset, the logistic regression was applied and a value for the likelihood ratio test statistic was obtained. The empirical distribution of the likelihood ratio test statistic obtained by permutation was then used to obtain a permutation‐based P‐value for the actual observed likelihood ratio test statistic on each SNP. The P‐value was calculated as the fraction of permutations that yielded a likelihood ratio test statistic as high or higher than the observed likelihood ratio test. We then calculated the correlation between the asymptotic P‐value and the permutation‐based P‐value. We obtained a strong correlation (r 2 = 0.86) between both types of P‐values.

Annotation and ranking of contigs

We used homology searches to annotate contigs containing an open reading frame (ORF). ORFs were submitted to similarity search by BLASTP against the NR protein database and GO functional annotations were obtained using the blast2GO program. GO terms were grouped in broader GO categories using the generic GOSlim vocabulary (v.1.2) available from www.geneontology.org

We ranked SNP containing contigs for their overall significance of association of SNP(s) frequency and either the drought or the temperature covariate. A score for overall significance was computed for each covariate on every contig containing SNPs as Z = − (log(P 1) + … + log(P ))/m; where P's are P‐value obtained on each of m SNPs using the logistic regression detailed above.

Results

Data generation and species barcoding

For each of the 15 samples (1 sample = 30 individuals pooled), 24–68 millions of 100 bp reads per sample were generated, of which 93.9–96.1% passed quality control filtering. Reads were assembled into 593 143 contigs, of which 19 122 included a predicted ORF of length above 200 bp and sequenced at average depth (coverage) above 20× per pool (Fig. 2a). We use that stringent cut‐off for coverage of contigs in order to have both the ability to detect SNP and to obtain good estimates of allele frequency at the SNP positions within each pool.

Figure 2

(a) Summary statistics of the transcriptome assembly: genome‐wide distribution of variation in genetic diversity from contig to contig within each pool. For each pool, genetic diversity is measured, within each contig, as the number of variable (SNP) position detected per 1000 nucleotides surveyed. Distributions are reported as smoothed histograms (green: pools sampled in ambient plots, orange: drought; red: drought + temperature + CO 2 plots). Inset graph: (smoothed) empirical distribution of contig length for all contigs (gray), and contigs containing open reading frames (blue). (b) Effect of experimental treatments on two covariates measuring the abiotic environment experienced by the 15 pools used for this study: mean temperature at 5 cm depth and mean soil water content during drought period at 0–20 cm of depth. Each plot is represented by a dot with a color reflecting the type of plot sampled (green: ambient, orange: drought treatment, red: drought + temperature + CO 2).

The taxon Cognettia sphagnetorum has recently been subdivided into four nominal species and its generic name changed into Chamaedrilus (Martinsson et al., 2015). We checked that all individuals in our study were coming from a single species.

To do so, we used the fraction of our sequencing (RNA‐seq) data that corresponds to mitochondrial gene Cox‐1 used to do barcoding; i.e. operationally define species. We mapped all reads obtained from RNA‐seq against a series of reference Cox‐1 sequences of annelids (published sequences available in July 2013 and a series of sequences provided by Rüdiger Schmelz, personal communication, see methods). We found that all reads (across all samples typically 99.9% and always more than 98%) that had any similarity to Cox‐1 are virtually identical at the nucleotide level to a single reference sequence. Cox‐1 is a highly expressed gene (>10 000× coverage in all samples analyzed). This gives us enormous power to detect contamination even by a single individual from a divergent species. Assuming conservatively a minimum coverage of 10 000 in any of our RNA seq samples (1 sample = 30 individuals pooled), even if only 1 individual of 30 in the sample originates from a different species, then we expect conservatively 300 reads should contain obvious sequences differences. The probability of missing just by chance such a divergent sequence is astronomically low (

Furthermore, we realigned the Cox‐1 reference sequence we used to barcode our samples – labeled JN259763 – with a set of Cox‐1 sequences published more recently: http://datadryad.org/resource/doi:10.5061/dryad.6mh29

Redoing a phylogenetic tree based on this sequence alignment shows unambiguously that our reference sequence is placed within the so‐called ‘sphagnetorum clade C’ in Martinsson & Erséus (2013, Figure 2) and renamed Chamaedrilus chlorophilus (Friend 1913) more recently (Fig. S3). Note that all other clades exhibit at least 7–10% divergence for that gene.

SNP diversity detected among populations

Reads were mapped to contigs and base counts were performed to characterize allelic variation across pools at each position of each contig. A total of 21 781 coding single‐nucleotide polymorphisms (SNPs) were called within 2629 distinct ORF‐containing contigs.

The detected amount of genetic variation, as measured by the density of SNPs in the contigs, was similar among the 15 samples irrespective of their environment of origin (Fig. 2a). Our pool‐based sequencing approach precludes estimating the numbers of distinct genotypes present in each population and no systematic difference in amounts of genetic diversity among populations from different treatment was found; but we reveal ample nucleotide variation in expressed regions of the genome upon which selection can act.

Association between SNPs and environmental variation

To characterize the abiotic environment of each plot studied, two covariates were calculated based on hourly soil temperature and moisture data. These covariates revealed substantial variation in mean temperature and drought level within the soil corresponding to the same climate treatments (Fig. 2b). The existence of heterogeneity between replicates of the same treatment motivates the use of these two underlying covariates instead of mere treatments effect in our SNP analysis model (see methods). So instead of trying to use treatment as a factor to explain variation in SNPs among populations, we chose to use directly these two covariates to characterize the abiotic environment of each population sampled as they describe the effect of the experimental manipulation carried out for both temperature and humidity in the soil. Our main goal was to test whether the variation in SNPs frequencies revealed between populations does associate statistically with these covariates. Here, temperature and soil moisture are taken as proxies for abiotic environmental variation that potentially creates spatially varying natural selection – through differential viability or fecundity of genotypes – and in turn triggers a detectable evolutionary response. If a covariate measures an environmental variation that does not trigger any evolutionary response, we expect SNPs frequency to vary randomly from population to population (through mere genetic drift) and no statistical effect of that covariate. If, however, C. chlorophilus populations responded evolutionarily to an environmental gradient quantified by a covariate, then a statistical association across samples between local SNPs frequency and the environmental covariate should be detected (Fig. 3).

Figure 3

Testing for micro‐evolutionary response to climatic gradients. Two examples of variation in allele frequency at individual SNPs along an environmental gradient measured through (a) the mean temperature at 5 cm depth or (b) mean soil water content during drought period (0–20 cm depth). Dots indicate actual observed SNP frequencies in ambient (green), drought (orange), and drought + CO 2 + temperature plots (red) and the predicted frequencies (black) using a logistic regression model including the covariate.

Every SNP detected in coding sequence (CDS), 5′UTR and 3′UTR regions of contigs containing an ORF was analyzed using a logistic regression in which the response variable (number of reads of a given allele) in each sample is modeled as a function of an environmental covariate (see methods). A P‐value for the test of environmental association is obtained for each SNP and a false discovery rate (FDR) analysis was conducted using the empirical distribution of P‐values (Fig. 4) to account for multiple testing (n = 21 781 SNPs tested). From that analysis, we estimate that 2700 SNPs from 1192 contigs are significantly affected by soil temperature, and 1208 SNPs from 664 contigs significantly affected by soil moisture at an expected FDR of 1%, thus demonstrating a pervasive effect of the environment on variation in allele frequency in this 7‐year long experiment.

Figure 4

Empirical distribution of P‐values for the association of SNPs with (a) the temperature (b) the drought covariate (right). Red line depicts the expected distribution of P‐values under the null hypothesis of no effect of the underlying environmental covariate on SNP frequency. Inset: Venn diagrams describing the amount of overlap between contigs carrying SNPs with significant association in either the 5′UTR, the CDS or the 3′UTR region (counts report the number of contigs in each category).

Although we could recover homologs from two other annelid draft genomes (Capitella teleta and Hellobdella robusta) for ca. 8000 contigs, only 802 contigs were homolog to functionally annotated sequences yielding 3266 gene ontology (GO) terms (Fig. S4). A list of contigs that could be annotated by homology and contain SNPs exhibiting the highest association with environmental covariates are given Table S2. Unfortunately, the majority of contigs containing the strongest signal could not be annotated. For instance the first two contigs that could be annotated and harbor SNPs with the strongest signal of covariation with the temperature covariate are only ranked 7 and 31 (Table S2).

We also computed the amplitude of allele frequency, defined as the difference between maximum and minimum values among the 15 analyze populations, as a proxy for of the overall amount of evolution across the experiment. Although we cannot assume that populations were initially harboring the same allele frequency, we take this statistic as our best (albeit indirect) measure of amount of evolutionary change during the course of this evolutionary experiment. If the detected SNPs were associated with environmental covariates ‘just by chance’, we would expect the distribution of these amplitudes to be the same between environment‐responding and environment‐non‐responding SNPs. Here we report a strikingly different picture (Fig. 5): SNPs whose frequency variation is associated with specific environmental covariates varied dramatically in frequency from sample to sample relative to SNPs for which no association was detected.

Figure 5

Distribution of amplitude in the minor allele frequency at SNPs along a gradient defined through either (a) temperature or (b) drought. SNPs with a low associated FDR (q < 0.001) and thus responding to selection are categorized as rejecting the null (HA, in red) while SNPs displaying no association with the environmental gradient when tested (P > 0.5) are categorized as stemming from the null distribution (H0, in gray).

Discussion

We found that ~12% of SNPs with differences in allele frequencies among populations can be associated to measurable changes in the abiotic environment. We document clear genetic differences among populations (Figs 3 and 5) subject to mean temperatures increased by only 0.5 °C over the 7‐year period in which soil temperature was manipulated. Importantly, such differences arose in a randomized and replicated experiment. We emphasize that the genetic differences we detect cannot not merely be driven by more intense bottlenecks in some environmental conditions as all pools exhibited extremely similar patterns of genome‐wide genetic diversity (Fig. 2a). Furthermore, although environmental effects can affect gene expression levels, pools of individuals were kept in the same lab environment prior to RNA extraction and RNA seq. This explains probably why we did not detect any significant effect of individual treatment on total gene expression level among populations (not shown). Note also that the statistical model we used – logistic regression with binomial link for SNP counts – naturally takes into account differences among population in gene expression levels (whether these are significantly different between treatment or not) and thus does not bias our statistical test for evolutionary response for a given SNP.

We argue in more detail below why the genetic differences we observe represent a genuine evolutionary change induced by selection and cannot be merely attributed to confounding of spatial and environmental proximity. Firstly, our approach is expected to be more robust than FST‐based approaches that focus on SNPs that are outliers displaying extreme genetic differentiation relative to the rest of the genome as candidate for local adaptation. Such methods make rather strong assumptions regarding the (unknown) underlying population structure to model allelic frequencies and are very prone to yield false positives (Bierne et al., 2013; Savolainen et al., 2013). In contrast, within our randomized experimental design, testing for association between variation in SNPs frequency and environment is a lot more straightforward because we are not dealing with a naturally occurring environmental gradient where typically populations that are geographically close are more likely to share identical genes (through isolation by distance) and closer environments (Coop et al., 2010).

Secondly, in our design, migration may occur between populations. However, migration events, if any have occurred, will merely dilute the signal of covariation between SNP frequency and the environmental covariate. Any findings of association between SNPs and environment variation cannot be spuriously generated by uncontrolled migrations between plots (these are unlikely given the dispersal abilities of our study organism, but cannot be ruled out).

Lastly, we do not have temporal series in our experiment. One concern might therefore be that the association we detect between SNPs and environment manipulations are merely caused by past population differentiations predating our climatic manipulation and not evolutionary changes. However, if this were the case, we would expect zero correlation between our likelihood ratio test P‐values and those obtained via permutations of environmental covariate. We observe a strong correlation (r 2 = 0.86) between both types of P‐values. This attests that the associations between SNPs and environment are robust findings and not artifacts of previous population differentiations. Even if the absence of a temporal sampling in our study precludes the direct estimation of the strength of selection underlying the adaptive changes of SNPs frequencies (Olson‐Manning et al., 2012), the evolutionary response we detect in our design is genuine.

The two covariates used for soil temperature and moisture exhibit some level of correlation (Fig. 2b). This, and the clonal nature of the populations we studied preclude us from inferring precisely what abiotic factor(s) and which SNPs are driving the micro‐evolutionary changes detected between populations. Similarly, the fact that annelids are understudied in genomics precludes parsing precisely which biological functions are associated with the genes exhibiting the observed evolutionary response. These are inherent limitations of our study system.

Can species adapt to climatic changes at a rate that matches the rate of climate change? This pressing question is pivotal for understanding the response of species to current climatic changes and future distribution of biodiversity. Studies documenting convincingly genetic change mediated by climate change selection remain scarce (Merilä, 2012). A handful of studies have used an experimental evolution approach to examine evolutionary response in model organisms propagated under lab conditions by typically altering either CO2 (Collins & Bell, 2004) or temperature within a realistic range given current climate change models. To our knowledge, the present study provides the first analysis of the evolutionary response within populations in a replicated experiment where climate change conditions are enforced under natural conditions.

Our approach demonstrates that the sequencing of samples from a replicated manipulated experiment provides sufficient power to measure the extent to which populations respond genetically to environmental perturbations under field conditions. This was achieved in a species lacking any prior genomic resources, which indicates that our approach can be carried out in a wide range of nonmodel organisms with limited capacity for dispersal – provided that a sufficient number of generations can be reached.

Incorporating evolutionary response into climate change modeling is pivotal for correctly predicting the response and future distribution of plant and animal species. Climate change models incorporating evolutionary responses often yield very different predictions for future species distributions compared to models that do not incorporate evolution (Kearney et al., 2009). Two scenarios may be foreseen: either no adaptation occurs and genetic diversity decreases through repeated bottlenecks ultimately increasing the risk of extinction in a changing environment, or, adaptive genetic changes do occur and may pave the way to evolutionary rescue. We did not observe any environment specific loss of genetic diversity in our study species, and the genetic differences we observe between environments therefore most likely represent adaptive changes to a future climate. If this is the case, and if our findings apply to many other species, the good news might be that models neglecting evolution when predicting species distributions according to climate change projections are perhaps too conservative and overestimating the risk of extinction. As more studies like this accumulate, comparing their outcomes will help assessing how pervasive evolutionary response are to climate change and ultimately reveal which biological factors limit the potential for evolutionary response and rescue.

Conflict of interest

Authors declare that they have no conflict of interest.

Author contributions

TB, BE, JGS, and MH conceived & designed the experiment. NF, DS, and SS performed RNA extraction, cDNA library construction & sequencing. TNM, KSL, and CB designed the field experiment and provided physical data. TB, AB, NG, and NC analyzed the data. TB, BE, and MH wrote the article with contributions from all other co‐authors.

Data availability

Contig sequences and allele counts data for the analysis of the SNP allele frequencies are deposited in the Dryad archive. http://dx.doi.org/10.5061/dryad.00v07

Table S1. Overview of number of SNPs called by two different methods: pools consisting of 10 individuals versus separate individuals.

Table S2. ranking and annotation of the 10 highest ranking annotated contigs.

Figure S1. Comparison of degree‐days accumulated in experimental plots in warming (T) versus ambient (A) plots.

Figure S2. Individual versus pool based SNP calling.

Figure S3. Maximum likelihood phylogeny of the Cox‐1 gene in the focal study species and known related genera, including the reference sequence (JN259763) used to barcode our samples.

Figure S4. Proportion of GO SLIM categories for all annotated contigs, the subset of contigs containing SNPs and tested for potential association, and the subset with at least SNP significantly associated with the meanTS5 covariate (FDR <1%).

Click here for additional data file.