Genetic and phenotypic divergence in an island bird: isolation by distance, by colonization or by adaptation?

Discerning the relative roles of adaptive and nonadaptive processes in generating differences among populations and species, as well as how these processes interact, is a fundamental aim in biology. Both genetic and phenotypic divergence across populations can be the product of limited dispersal and gradual genetic drift across populations (isolation by distance), of colonization history and founder effects (isolation by colonization) or of adaptation to different environments preventing migration between populations (isolation by adaptation). Here, we attempt to differentiate between these processes using island populations of Berthelot's pipit (Anthus berthelotii), a passerine bird endemic to three Atlantic archipelagos. Using microsatellite markers and approximate Bayesian computation, we reveal that the northward colonization of this species ca. 8500 years ago resulted in genetic bottlenecks in the colonized archipelagos. We then show that high levels of genetic structure exist across archipelagos and that these are consistent with a pattern of isolation by colonization, but not with isolation by distance or adaptation. Finally, we show that substantial morphological divergence also exists and that this is strongly concordant with patterns of genetic structure and bottleneck history, but not with environmental differences or geographic distance. Overall, our data suggest that founder effects are responsible for both genetic and phenotypic changes across archipelagos. Our findings provide a rare example of how founder effects can persist over evolutionary timescales and suggest that they may play an important role in the early stages of speciation.


Introduction
Understanding the processes that drive differentiation among populations is a major aim in evolutionary biology. Most population genetic studies have been carried out using neutral genetic markers, at which variation across populations has been viewed as the product of mutation, migration and genetic drift (Wright 1943;Hartl & Clark 2007). As migration tends to occur more commonly between neighbouring populations, a relationship between the geographic and genetic distance of pairs of populationsisolation by distanceis expected in many wild systems (Wright 1943(Wright , 1946. In contrast, most phenotypic differences across populations are assumed to have arisen in response to differential selection regimes, and drift is considered a null hypothesis against which selection can be tested (Orr 1998).
Recently, several authors have questioned whether patterns of neutral genetic differentiation can be understood in terms of migration and drift alone and have suggested that selection can play a role (e.g. Nosil et al. 2009;Orsini et al. 2013a). In cases where populations exist along environmental gradients, individuals are likely to be adapted to local conditions. This local adaptation may act as a barrier to migration, and subsequent drift will result in neutral genetic structure across populations. This phenomenon, known as 'isolation by adaptation', has received both theoretical and empirical support, and a recent review suggests that it may be substantially more common than previously realized (Orsini et al. 2013b). Importantly, just as selection may structure neutral DNA, genetic drift may shape populations at the phenotypic level (e.g. Carson 1968;Kolbe et al. 2012). However, because drift is usually considered a null hypothesis, its role in shaping differences among populations and species remains poorly understood (Marie Curie SPECIATION Network 2012).
A further related factor that must be considered in studies of population differentiation is colonization history. The founding of a population by a small number of individuals can result in founder effectsgenetic and/or phenotypic divergence from the parent population. This can lead to a pattern of 'isolation by colonization', whereby differences between populations primarily reflect colonization history, rather than contemporary patterns of migration, gradual genetic drift or selection. Laboratory and field experiments designed to simulate founder events have shown that foundermediated genotypic and phenotypic divergence can (but do not necessarily) occur and persist in the presence of selection (Moya et al. 1995;Travisano et al. 1995;Rundle 2003;Kolbe et al. 2012). Cases where recent natural or human-mediated introductions have been observed have produced mixed results, with the extent of founder divergence depending on the timing, severity and extent of bottlenecks (Carson 1968;Brinkmann et al. 1998;Coyne & Orr 2004;Frantz et al. 2009;Kolbe et al. 2012). Several recent studies have used modelling approaches to study isolation by colonization, again with mixed results. In humans, patterns of neutral genetic differentiation are largely the product of serial founder effects that arose as small populations colonized out of Africa (Prugnolle et al. 2005a). These founder effects have also shaped patterns of morphological variation and resulted in structure at disease resistance genes, both of which are traits under selection (Prugnolle et al. 2005b;Manica et al. 2007). In silvereyes (Zosterops lateralis), a recent island colonizer, sequential founder events produced genetic divergence among populations (Clegg et al. 2002a;Estoup & Clegg 2003). At the morphological level, however, selection was found to be the primary force driving differences among populations (Clegg et al. 2002b(Clegg et al. , 2008. The presence of isolation by distance can be detected by correlating pairwise levels of genetic structure (e.g. F ST ) with pairwise geographic distances between populations, usually using Mantel and partial Mantel tests (Rousset 1997). Similarly, isolation by adaptation can be tested by correlating pairwise F ST values with a matrix of ecological or environmental distance (Nosil et al. 2009). However, it will only be possible to differentiate between isolation by distance and isolation by adaptation using this approach when ecological and geographic distances between populations are not strongly correlated. Under isolation by colonization, no relationship is expected between genetic distance and either geographic or ecological distance (Orsini et al. 2013b); instead, genetic structure is expected to reflect the colonization and bottleneck history of the populations being studied. Patterns of phenotypic divergence (e.g. Q ST ) can be studied in the same way as neutral genetic data. However, there is a lack of studies simultaneously examining isolation by distance, by colonization and by adaptation at the genetic and morphological levels (Orsini et al. 2013b). Therefore, while it is clear that these may all occur, their relative importance and how these processes interact are not well understood.
In this study, we investigated patterns of genetic and phenotypic divergence in Berthelot's pipit (Anthus berthelotii), a small passerine endemic to the Canary Islands, Madeira and Selvagens archipelagos in the North Fig. 1 The location of Berthelot's pipit populations used in the present study, and (inset) the location of the three archipelagos in relation to the mainland (CI = Canary Islands, SG = Selvagem Grande; reproduced from Illera et al. 2007). Atlantic (Fig. 1). Previous analyses of neutral and functional genetic diversity across Berthelot's pipit populations have revealed very low levels of variation and a recent northward colonization across its range, from the Canary Islands to Madeira and the Selvagens Spurgin et al. 2011). We use a panel of 21 autosomal microsatellite markers and population genetic simulations to test whether the patterns of neutral genetic divergence across the pipit populations are better explained by isolation by distance, by adaptation or by colonization. We then compare patterns of morphological and neutral genetic variation to determine the evolutionary forces driving divergence at the phenotypic level.

Sampling, morphological and molecular methods
Representative samples (ca. 30 individuals per population) were obtained from all 13 Berthelot's pipit populations in April 2005 (Selvagens), January to March 2006 (Canary Islands) and September 2006 (Madeiran Islands). Birds were caught from multiple localities encompassing the entire range of each population, ensuring that population-level comparisons could be made (see below). Sampled birds from all populations were from low-lying coastal regions, with the exception of birds caught on the mountain of El Teide on Tenerife. These birds exist on a plateau ca. 2000 m above sea level and are separated from the rest of the island by a band of dense forest, which the pipit does not inhabit. For this reason, we considered El Teide a separate thirteenth population. While Berthelot's pipits do exist on alpine plateaus in other islands, such as La Palma and Madeira, these populations are much smaller and were not targeted for trapping. We caught all birds using spring traps baited with Tenebrio molitor larvae, fitted each with a uniquely numbered aluminium ring from the relevant national authority, or with a plastic colour ring when these were not available, and took seven morphological measurements (mass, tarsus length, wing length, head length and bill length, width and height; see Illera et al. 2007 for details). Blood samples (ca. 40 ll) were taken by brachial venipuncture and stored at room temperature in screw-topped, rubber-sealed, microfuge tubes containing 800 ll of absolute ethanol.
Genomic DNA was extracted using a modified salt extraction method (Bruford et al. 1992;Richardson et al. 2001) and diluted to a concentration of 10-50 ng/ll. We genotyped all individuals at a panel of 21 microsatellite loci (Primmer et al. 1995;Martinez et al. 1999;Dawson et al. 2010Dawson et al. , 2012Dawson et al. , 2013 (Table S2, Supporting information). Loci were confirmed unique by comparing their sequences using BLASTN software (Altschul et al. 1997).
PCRs were carried out in 2-ll reactions using a method based on Kenta et al. (2008). Briefly, 1 ll (10-50 ng) of genomic DNA was added to each PCR well and the liquid evaporated. To this we added 1 ll primer mix (containing all forward and reverse primers in the multiplex reaction at 0.2 lM) and 1 ll 2x QIAGEN Multiplex PCR Master Mix. Samples were overlaid with mineral oil before PCR. An initial hot-start denaturing phase of 95°C for 15 min was followed by 40 cycles of 94°C for 30 s, 56°C for 90 s and 72°C for 60 s. A final hold of 60°C for 30 min completed the reaction. PCR products were diluted 1 in 400 and separated on an ABI 3730 DNA analyser. We determined allele sizes using GENEMAPPER version 3.7 (Applied Biosystems).

Genetic analyses
Unless stated otherwise, all statistics and plots were generated in R version 2.15.2 (R Development Core Team 2012). For all variables used in parametric tests, we tested for normality using Shapiro-Wilks tests and used transformations where appropriate.
For each microsatellite locus and population, Hardy-Weinberg equilibrium and linkage disequilibrium were tested using GENEPOP version 4.1 (Raymond & Rousset 1995) and null allele frequencies were estimated using CERVUS version 3.0.3 (Marshall et al. 1998). Allele frequencies and expected heterozygosity were calculated using FSTAT version 2.9 (Goudet 1995). Allelic richness was calculated after controlling for differences in sample size, using a rarefaction approach implemented in HP rare (Kalinowski 2004(Kalinowski , 2005. We tested for differences in genetic diversity across archipelagos using linear mixed models, implemented in the lme4 and nlme packages in R (Bates et al. 2011;Pinheiro et al. 2013). Allelic richness and heterozygosity were included as response variables in separate models, with archipelago as a fixed factor and locus identity as a random factor. To assess whether population identity had a significant effect on genetic diversity above and beyond the effect of archipelago identity, we then ran a second mixed model including the same variables, plus population identity as a random factor, and compared the two models using likelihood ratio tests.
To visualize the level of genetic structure across islands and archipelagos, we performed a principal component analysis (PCA) of the 21 microsatellite loci, using the adegenet package in R (Jombart 2008). We also used STRUCTURE (Pritchard et al. 2000) to identify the number of distinct genetic clusters (K) in our data set. For this, we used a model allowing admixture and correlated gene frequencies and carried out four independent runs for each value of K = 1-13. We accepted the value of K with the highest average 'log probability of data' as the most likely number of clusters. For each run, we ran 500 000 steps, with a burn-in of 10 000 steps.
We used two approaches to detect whether populations had undergone genetic bottlenecks. First, we calculated the degree of heterozygosity excess, which occurs due to the loss of rare alleles shortly after bottlenecks (Cornuet & Luikart 1996), in the program BOTTLE-NECK (Piry et al. 1999). We used a two-phase mutation model, with 95% stepwise and 5% nonstepwise mutations. The probability of heterozygosity excess was then calculated using Wilcoxon tests. Second, we calculated Garza and Williamson's index (M) by dividing the number of alleles in a population (k) by the range in allele size (r) (Garza & Williamson 2001). This was modified to M = k/(r + 1) to avoid dividing by zero in monomorphic populations (Excoffier et al. 2005).
We used the program DIY-ABC (Cornuet et al. 2008) to reconstruct the most likely demographic history of Berthelot's pipit. We constructed the four most plausible model colonization scenarios based on a northward expansion across the three archipelagos ). These involved two colonization routes from the Canary Islands, one with a sequential colonization to Selvagens then to Madeira and the other with a simultaneous colonization of Selvagens and Madeiraeach with and without bottlenecks associated with the colonization events. The models described an ancestral N e , as well as the contemporary N e and the colonization time of each archipelago. For the two scenarios involving bottlenecks, we also included the number of founders for each archipelago and the duration of the bottleneck. Priors were given uniform distributions, informed where possible by knowledge of the past and present distribution of Berthelot's pipit (Table S2, Supporting information). We used the default settings provided by the program as prior information on microsatellite mutation rates. We simulated one million data sets for each scenario and calculated the following summary statistics for observed and simulated data sets: mean number of alleles per locus, mean gene diversity, mean size variance and mean M ratio within each population, as well as pairwise F ST across each pair of populations. The posterior probability of each scenario was estimated using a logistic regression approach implemented in DIY-ABC. We evaluated confidence in our choice of scenario by comparing the match between simulated and observed summary statistics and by calculating type I and II error rates using the standard procedures available in DIY-ABC (Cornuet et al. 2008).

Isolation by distance, colonization and adaptation
We tested whether patterns of neutral genetic structure were the product of isolation by distance, by adaptation or by colonization, using a rationale largely based on that of Nosil et al. (2009) and Orsini et al. (2013b). Under a scenario of limited dispersal, a linear relationship between pairwise genetic and geographic distance matrices is expected. Under isolation by adaptation, we expect patterns of genetic structure to reflect ecological differences across islands, irrespective of geographic distance. Finally, under isolation by colonization, we expect patterns of genetic structure to mirror the colonization history of the pipit, with highest levels of genetic structure between populations that have been separately colonized by a small number of founders (i.e. the most bottlenecked populations).
We calculated population-level pairwise genetic differentiation as F ST /(1ÀF ST ) (Slatkin 1995) in ARLEQUIN version 3.5 (Excoffier & Lischer 2010). Geographic distance was calculated as the closest linear distance between pairs of islands, estimated using Google Earth (http://earth.google.com). As a proxy for ecological distance, we used a selection of climate variables. For each island, we extracted data on average monthly precipitation, average temperature in August (the hottest month) and average temperature in January (the coldest month; see Appendix S1 for sources, Supporting information). Our rationale for choosing these variables was based on the extent to which they differed within and among archipelagos, and on their expected relationship with divergence at the morphological traits we had measured. Temperature is the most important predictor of body size in birds globally (Olson et al. 2009), and under isolation by adaptation, we therefore predict that birds on islands with similar temperature profiles will be of similar size. Annual rainfall has been shown, via food availability to result in natural selection on beak morphology in birds (Grant & Grant 1996), and has been linked to arthropod abundance in the Canary Islands (Illera & Diaz 2006). We standardized climate variables (v) using the following equation: Using the three standardized climate variables for each island, we then calculated Euclidean distances between all pairs of islands, thus obtaining a single distance matrix (hereafter referred to as 'environmental distance').
To obtain a distance matrix that reflects the relative severity of the combined bottleneck effects between two populations, we modified Garza & Williamson's (2001) M ratio as follows: M ðpopulation pairÞ¼À logðM ðpop1Þ ÂM ðpop2Þ Þ Thus, M (population pair) is highest when both populations are bottlenecked (and M is low) and is therefore hereafter termed 'bottleneck severity'.
We tested separately whether genetic distance was related to geographic distance, bottleneck severity and environmental distance using Mantel tests, implemented in the Ecodist package in R (Goslee & Urban 2007). Importantly, environmental and geographic distance in this system are not correlated (Mantel test, R = 0.05, P = 0.3), making it possible to assess the relative importance of these variables. To control for the potentially confounding effect of hierarchical structure (i.e. of archipelago-level effects) on relationships between distance matrices, we ran partial Mantel tests, including a coding variable that identified which archipelagos were included in each pairwise comparison e.g. Canary Islands-Canary Islands, Canary Islands-Madeira, etc).

Morphometric analyses
To visualize the level of morphological structure across archipelagos, we conducted a PCA of six of the seven morphological measurements (excluding mass, which fluctuates too much over short timescales in small passerines) and plotted the first two principal components (PC1 and PC2). We tested how morphological variation was partitioned within and across archipelagos; for each PC as the dependent variable, we first ran a linear model, including age, sex and archipelago identity as fixed factors, and we then ran a mixed model, including the same variables, plus population identity as a random factor. To assess whether island identity had a significant effect on morphology above and beyond the effect of archipelago identity, we compared the linear mixed models to the linear models using likelihood ratio tests.
We tested whether morphological differences between populations were the product of isolation by distance, by colonization or by adaptation using the same rationale as for the genetic data. We calculated Euclidean distances in morphology between all pairs of populations, using the mean values per population of the six morphometric traits (i.e. we obtained a single matrix of Euclidean distances in six-dimensional space). Morphometric traits were standardized in the same manner as environmental variables (equation 2). We tested whether morphometric distance was related to geographic distance, bottleneck severity and environmental distance using Mantel tests and partial Mantel tests (accounting for archipelago-level differences). We also directly tested whether morphometric and genetic distance matrices were related, again using Mantel and partial Mantel tests.

Results
In total, 371 individuals were genotyped at 21 unique microsatellite loci (sample sizes for each population provided in Table 1). No loci were in Hardy-Weinberg disequilibrium in more than two of the 13 populations typed, no loci had estimated null allele frequencies greater than 0.2, and no loci were in linkage disequilibrium after Bonferroni correction for multiple tests. Allele frequencies are given in Table S3 (Supporting information). Allelic richness and heterozygosity both differed significantly across archipelagos (Fig. 2), with highest levels of variation found in the Canary Islands and lowest levels in Selvagens (F = 27.9; P < 0.001). Including population identity as a random factor did not significantly improve the model fit for either allelic richness or heterozygosity models (likelihood ratio tests, both P-values >0.99), suggesting that differences in genetic diversity were predominantly across, compared to within, archipelagos (Fig. 2). The two tests employed to detect genetic bottlenecks yielded concordant results: there was no evidence for a bottleneck in any of the Canary Island populations, whereas there was a bottleneck signal in the Selvagens, and in all three Madeiran Islands (Table 1).
A PCA of the 21 microsatellite loci revealed high levels of structure across, but not within, archipelagos (Fig. 3A). Likewise, Bayesian clustering revealed K = 3 as the most likely number of genetic clusters, and these clusters corresponded to the three archipelagos ( Figure  S1, Supporting information). There was an overall relationship between pairwise genetic and geographic distance matrices (Mantel test, R = 0.39, P = 0.007), suggesting a pattern of isolation by distance; however,  Fig. 4A). Genetic distance was strongly related to bottleneck severity (R = 0.86, P < 0.001), and this relationship remained highly significant after controlling for archipelago effects (R = 0.89, P < 0.001; Fig. 4B). Genetic distance was not related to environmental distance, regardless of whether archipelago effects were accounted for (R = À0.14, P = 0.89 and R = À0.14, P = 0.90 respectively; Fig. 4C). The most highly supported demographic scenario in DIY-ABC analyses involved simultaneous colonization of Madeira and Selvagens both directly from the Canary Islands, with bottlenecks associated with the founder events (scenario 1, posterior probability = 0.96; Table 2). Posterior estimates from scenario 1 suggested moderate contemporary N e in all three archipelagos and that Madeira and Selvagens were colonized from the Canary Islands between ca. 1000 and 26 000 years ago (Table 3). While estimates of the bottleneck N e were very broad (4-49) given the prior distribution (1-50), the simulations did suggest that Madeira and Selvagens remained at a low N e for a long time following colonization, with the length of the bottleneck estimated at between 200 and 12 000 years (Table 3). DIY-ABC was unable to estimate the colonization time from the mainland to the Canary Islands with any degree of accuracy (25 000-2 000 000 years ago), most likely because of a lack of samples from the ancestor of Berthelot's pipit. Using scenario 1, all but one of our simulated summary statistics were within the range or our observed summary statistics at the P > 0.05 level (the M ratio in the Canary Islands was the only exception to this). Rates of type I and II error were low, 0.044 and 0.040, respectively.
Variation across archipelagos for each of the six morphological traits is shown in Figure S2 (Supporting information). As with the genetic data, morphometric variation was greater across than within archipelagos. The first two components from a PCA of the morphological traits explained 83% and 12% of the total variation, respectively. PC1 was moderately negatively correlated with all six original variables and thus represents body size. PC2 was highly negatively correlated with bill width and height and thus represents bill shape. We detected highly significant differences across archipelagos in both PC1 (F = 47.6, P < 0.001) and PC2 (F = 46.8, P < 0.001). Madeiran birds were associated with increased body size (decreased PC1) with no change in bill shape, whereas on Selvagens birds exhibited decreased overall size and increased bill height and width (decreased PC2) (Fig. 3B). Archipelago-level differences in both PC1 and PC2 exceeded population-level differences (Fig. 3C); however, including population identity as a random factor did result in a significantly improved model for PC1 (likelihood ratio test, P < 0.001), but not for PC2 (P = 0.07). This effect on PC1 was caused by an increased body size in the mountainous population on El Teide, in Tenerife (Fig. 3C).

Discussion
Here, we show that these latter two archipelagos were colonized independently from the Canary Islands between 1000 and 26 000 years ago. Both basic population genetic analyses and Bayesian simulations indicate that the patterns of structure observed across the island populations are the result of these dispersal events, with genetic bottlenecks having occurred in the two most recently colonized archipelagos. Similarly, explicit tests for isolation by distance, colonization and adaptation revealed that only colonization history was a significant predictor of genetic structure across populations. We also found that patterns of morphological variation closely mirror patterns of structure at the neutral microsatellite loci and are consistent with the bottleneck history of the pipit, but not with a scenario of selection or gradual genetic drift. Together, our data suggest that founder effects have been the predominant force shaping genetic and morphological variation in this species. Revisiting our previous study ) with additional loci and new analytical methods has changed our interpretation of the Berthelot's pipit population history. The most striking result was the detection of genetic bottlenecks in the recently colonized populations (Table 1), which were not observed when five microsatellite loci were used . Simulations have highlighted that the two most commonly used methods for detecting genetic bottlenecks (Cornuet & Luikart 1996;Garza & Williamson 2001) are sensitive to type II error, depending on the number of loci used and the variability of these loci (Williamson-Natesan 2005) our findings here provide a striking example of how this can occur, and how it can influence our interpretations. In the light of this, that we did not detect a bottleneck in the Canary Islands does not mean that one did not occur there either, but it does suggest that if a bottleneck did occur, it was older and/ or less severe than in the Selvagens or Madeiran populations, which is in accordance with our colonization scenario.
The other result that contrasts with our previous study is that the apparent pattern of isolation by distance was no longer observed when archipelago-level effects were controlled for (Fig. 4A). This suggests that patterns of neutral genetic structure across Berthelot's pipit populations are in fact not the simple product of limited migration and gradual genetic drift. We also suggest that our patterns of neutral genetic structure are not the product of environmental differences limiting migration across archipelagos. Our evidence for this is the following: (i) we found no relationship between divergence at our selected environmental variables and genetic distance (Fig. 4C); (ii) our own observations indicate that ecological differences (e.g. food availability) are greater within compared to across archipelagos; and (iii) we do not expect a relationship between  genetic structure and bottleneck history (Fig. 4B) if adaptation is the primary driver of differences among populations. It can be very difficult to distinguish between true isolation by distance, isolation by colonization and isolation by adaptation, and we echo recent advice that more care needs to be taken in differentiating between these phenomena (Meirmans 2012;Orsini et al. 2013b).
In contrast to what one would expect from selection (Clegg et al. 2002b(Clegg et al. , 2008, the divergence observed at different morphological traits strongly reflected patterns of neutral genetic divergence (Fig. 3). As with our genetic data, most of the differences in morphology were across archipelagos, despite there being greater variation in ecological parameters within populations, or among populations within archipelagos, than there is among archipelagos. Morphological divergence was not related to geographic or environmental differences between populations, but was strongly related to the bottleneck history of the pipit (Fig. 4). While it will never be possible to rule out that some environmental or ecological variables that we have not considered, such as diet, are driving morphological (and genetic) differences across archipelagos, the most parsimonious explanation is that isolation by colonization has been the predominant force driving both genetic and phenotypic divergence in Berthelot's pipits. Importantly, the relationship between morphometric divergence and bottleneck history was strongly significant even after controlling for archipelago-level effects, suggesting that differences in morphology within, as well as among, archipelagos were largely driven by neutral processes.
An exception to this is the increased body size in the mountain population on El Teide, on the island of Tenerife (Fig. 3C). This population is situated >2000 m above sea level, and temperatures frequently drop below freezing in winter, contrasting starkly with the other populations, which are predominantly low-lying, subtropical and coastal. In passerine birds, body size and altitude are positively related, as suggested by Bergmann's rule (Blackburn & Ruggiero 2001). Thus, it seems probably that the colder weather on El Teide has selected for larger body size in Berthelot's pipits. We note, however, that overall differences in body size due to founder effects observed across the archipelagos far outweigh the differences between El Teide and the other Canarian populations (Fig. 3C).
Morphological changes in island birds constitute some of the classic examples of natural selection in the wild and can occur over remarkably short timescales (Grant 1986). It is somewhat surprising, therefore, that the apparent founder effects observed across Berthelot's pipit populations have not yet been obscured by natural selection. This may be because i) with the exception of El Teide (and the other, smaller isolated mountain populations on Madeira and La Palma, which we did not sample), this species occupies uniform habitats with benign and stable climate. Additionally, very few other insectivorous bird species can be found in the same habitat as Berthelot's pipits (Illera 2007), suggesting that interspecific competition is low. However, although natural selection does not appear to have been the primary driver of morphological change in Berthelot's pipits, it is almost certain that natural selection and sexual selection have acted on other genotypic and phenotypic traits than the ones we measured here. Indeed, we have previously shown that these pipits face different levels of pathogen-mediated selection across populations (Spurgin et al. 2012) and that selection is involved in differentiation at major histocompatibility complex (MHC) genes (Spurgin et al. 2011). This species therefore provides an excellent opportunity to study how selection and drift interact at the genotypic and phenotypic levels.
In addition to a lack of (or weak) selection on morphological change, it is thanks to a fortuitous evolutionary time frame and a lack of gene flow across archipelagos ) that we have been able to observe founder effects so clearly in this study. In very recently separated populations and experimental studies, in which it is possible to directly observe founder events, it is unclear how long any observed founder effects will persist (Kolbe et al. 2012). On the other hand, in populations or species that have been separated for very long periods of time, a combination of selection, mutation and drift will have produced Table 3 Mean and quantiles of posterior distribution samples for parameters estimated in DIY-ABC analyses of Berthelot's pipit colonization history. Parameters were estimated using the best-supported scenario (scenario 1, see Table 2) of an initial colonization of the Canary Islands, followed by simultaneous colonization from the Canary Islands directly to both Madeira and Selvagens, with bottlenecks associated with each colonization event differences so great that it becomes difficult to disentangle the signature of founder effects from other evolutionary forces (e.g. Hoeck et al. 2010;Aleixandre et al. 2013). The fleeting nature of founder effects may be one reason that their prevalence and importance in evolution are poorly understood (Templeton 2008). By studying moderately diverged populations, we have been able to show that founder effects can indeed persist over evolutionary timescales, and our results suggest that under certain circumstances, founder effects may play an important role in the initial stages of differentiation and speciation. Further targeted research using systems where founder speciation is predicted to play a role will help us to better understand its evolutionary importance.

Data accessibility
Microsatellite, morphological and climate data sets have been posted on Dryad (doi:10.5061/dryad.n1506).

Supporting information
Additional supporting information may be found in the online version of this article.
Appendix 1 Sources for climate variables used to calculate 'environmental distance' across Berthelot's pipit populations.
Table S1 Demographic scenarios used for approximate Bayesian computation (ABC) analyses (see Fig. 1 for visual representation).

Table S2
Twenty-one microsatellite loci characterised in the Berthelot's pipit.