Diversifying selection and color-biased dispersal in the asp viper
Study site and tissue sampling
The study site is located in the French Alps (Mont Blanc massif), between 1’100 and
2’100 m above sea level. We collected a total of 170 samples from blotched (N?=?113)
and non-blotched (N?=?57) snakes between 2006 and 2010. Based on their location, we
grouped the samples into 12 populations (see Figure 4; in red, populations with a proportion of non-blotched individuals higher than 50%).
Figure 4. Distribution and proportion of non-blotched individuals of sampled populations. In
red and green, populations mainly assigned to the first and second genetic cluster
identified by the STRUCTURE analyses (Hubisz et al.57]).
For each captured individual, we collected the coordinates either with a GPS or with
the help of GOOGLE EARTH v5.0 (Google, Mountain View, US), both methods allowing an accuracy of about 10 m.
We collected DNA samples for the genetic analyses from blood, ventral scales, tip
of the tail (stored in 90% ethanol prior to DNA extraction), and/or buccal swab. In
order to avoid duplicate samplings, we took dorsal and head pattern pictures (dorsal
and lateral view) for each specimen. Color morphs of individuals were determined in
the field. Several persons scored the phenotype of snakes simultaneously and none
of the 170 individuals used in the present study had an ambiguous phenotype. Indeed,
blotched (characterized by a zigzag-like pattern) and atypical non-blotched individuals
(light-colored individuals with an absence of zigzag-like pattern) were sufficiently
different to avoid misclassifications (Figure 1).
DNA extraction
We extracted DNA using the QIAGEN DNeasy® kit (QIAGEN, Hombrechtikon, Switzerland)
according to the DNeasy® Blood Tissue Handbook. In order to improve the quality
of extractions, we modified a few steps: (i) overnight lysis was conducted for most
samples, (ii) elution was performed twice using 100 ?l of Buffer AE and finally (iii)
incubation time was extended to 5 min so that a higher concentration of DNA was obtained.
For buccal DNA sampling, the DNA extraction was also following the same protocol,
except that swabs were placed into the DNeasy Mini spin column and centrifuged to
remove the remaining liquid from them between steps 4 and 5 of the protocol provided
by the manufacturer.
Microsatellite analyses
Seven loci from V. aspis (Va-P70, 91, 69, 35, 81, 25, 20; 49]) and two loci from V. berus (Vb-A8 and Vb-D17; 50]) were amplified and scored. We performed PCR in an Eppendorf Mastercycler Gradient
(Vaudaux-Eppendorf AG, Basel, Switzerland), with a final volume of 12 ?l per reaction,
using 2-4 ?l of DNA extraction 1x PCR buffer, 2 mg/ml of Q solution, 0.2 mM dNTPs,
0.5 units of Taq Polymerase, and, depending on the amplified locus, 1.5 – 3.5 mM of MgCl2 (all reagents from QIAGEN) and 0.25 – 0.5 mM of each primer (see Additional file
1 for more details). Cycling conditions included 37 to 40 cycles of 95°C for 30 sec,
50 to 57°C (depending on the microsatellite locus) annealing temperature during 30 sec, 72°C for 45 sec, and a final extension of 72°C for 7 min (see 49] for more details). Then, we genotyped amplified products with an ABI3130xl genetic
analyzer (Applied Biosystems) and visualized with PEAK SCANNERâ„¢ Software v1.0 (Applied
Biosystems).
F-statistics and genetic diversity parameters
First, we tested for the presence of null-alleles, scoring error due to stuttering
and large allele dropout with MICRO-CHECKER v 2.2.3 51]. We calculated the genotypic disequilibrium between loci in each sample based on
10,000 randomizations to check for linked loci. We tested deviations from Hardy–Weinberg
equilibrium (HWE) within samples based on 10,000 randomizations. We estimated the
Wright’s fixation indices for within-population deviation from random mating (FIS), as well as pairwise subpopulation differentiation (FST), following Weir Cockerham 52]. We computed deviations from random mating within populations (FIS) per locus and sample with a bootstrap procedure including 10,000 randomizations.
We estimated expected (HE) and observed (HO) heterozygosities following the methods of Nei Chesser 53] and allelic richness (AR) with FSTAT. In addition, we performed a Mantel test 54] with genetic distance (pairwise FST) as the dependent variable and the distance between sites as explanatory variable.
We carried out permutation tests in order to detect significant differences in allelic
richness, expected (HE) and observed (HO) heterozygosities and FST indices among the populations with a high versus a low amount of non-blotched individuals
(populations 1-4 and 5-12, respectively). We performed all summary statistics and
tests using the software FSTAT Version 2.9.3.2 55]. The critical p-value of 0.05 was adjusted using the Bonferroni correction 56] due to multiple comparisons.
Structure of populations
We used the software STRUCTURE version 2.3.4 57], a Bayesian model-based clustering method 58], to infer population structure and to assign individuals to populations. Based on
allele frequencies, we used this MCMC simulation to assign a membership coefficient
for each individual to each K populations. Ten runs of 600,000 iterations (the first
200,000 considered as burn-in) for K?=?1–12 were performed including all individuals.
Then, we defined the number of clusters that best fits our data set as described in
Evanno et al.35]. This approach compares the rate of change in the log probability of data between
successive K and the corresponding variance of log probabilities.
Unidirectional gene flow among populations
In order to quantify unidirectional migration rates (m) between populations, we used
the software BAYESASS 3.0.3 59]. This Bayesian method relies on the tendency for immigrants to show temporary disequilibrium
in their genotypes relative to the focal population, which allows their identification
as immigrants or offspring of immigrants. After initial runs were conducted with variable
values of deltaA, deltaM and deltaF in order to improve the acceptance levels, the
best values were set to deltaA?=?1, deltaM?=?0.2 and deltaF?=?1. For the final analysis,
we used 3 x 106 MCMC iterations, including a burn-in length of 3 x 105 iterations.
Sex and color-biased dispersal
We tested for sex- and color-biased dispersal in our populations, using GENALEX 6.5 60]. Then, we compared the mean of the corrected assignment index (mAIc; 61]) between sexes (males vs. females), as well as between colors (blotched vs. non-blotched
individuals). With this statistical approach, residents tend to have higher mAIc values
than immigrants.
Comparison between genetic and morphological variation
A commonly used method to estimate population differentiation for a quantitative trait
is a metric called QST, an analog of FST, which calculates the genetic differentiation at neutral genetic markers. Because
QST calculation requires experimental estimates of additive genetic variances and since
we did not estimate additive genetic variance for asp viper coloration, we will refer
in this study to PST (phenotypic or pseudo-QST) rather than QST, as proposed by Saether et al.62]. PST was estimated with two different methods. First, we coded the color as a locus having
only two alleles, blotched individuals being homozygotes coded 11 at this locus and
non-blotched individuals being homozygotes 22. Then, we used FSTAT to calculate the
pairwise PST between populations. Second, we estimated pairwise PST-values as a function of heritability (h2), within- and between-populations phenotypic variances (; 29],63]):
We obtained the within- and between populations variances by extracting the mean squares
(MS) with an ANOVA on color in JMP 10.0 (SAS Institute, Cary, NC, USA). Because the heritability
of the color pattern is not known, we used three different values of heritability
(0.1, 0.5 and 1). We estimated within-population variance without bias by within-population MS, whereas between-population variance is calculated as follows:
where MSb and MSw are the within- and between-population MS and n0 is a weighted average of the sample size for each population comparison estimated
following Sokal Rohlf 64] as:
where ni the number of individuals in the ith population and a is the number of populations to be compared.
To test whether the coloration differentiations between populations evolved by neutral
processes or selection, we compared PST values obtained for individual traits and with the two different methods with the
distribution of pairwise FST values. As in Slavov et al.65], we considered PST values as extreme when exceeding the empirical 95th percentile of the FST distribution (in this study: 0.086). In addition, to check for significant difference
between PST – and FST-values, we performed a Wilcoxon signed rank test 66].
Availability of supporting data
Microsatellite genotypes and phenotypes: Dryad doi: 0.5061/dryad.87478.
Animal ethics
The samples have been taken with the authorizations of the local authorities (Autorisation
préfectorale No. 2009-14; Direction de l’Administration Territoriale et de l’Environnement,
Bureau de l’Environnement et du Développement Durable, Préfecture de la Savoie, 73018
Chambéry, France).