Abstract

Numerous studies addressing the impact of habitat fragmentation on genetic diversity have been performed. In this study, we analyze the effects of a seemingly nonpermeable matrix on the population structure of the forest-dwelling butterfly Pararge aegeria in geographically isolated oases at the northern margin of the Sahara desert using microsatellites, morphological characters, and species distribution modeling. Results from all analyses are mostly congruent and reveal 1) a split between European and North African populations, 2) rather low divergence between populations from the eastern and western part of North Africa (Morocco vs. Tunisia), 3) a lack of differentiation between the oasis and Atlas Mountain populations, 4) as well as among the oasis populations, and 5) no reduction of genetic variability in oasis populations. However, one exception to this general trend resulted from the analyses of wing shape; wings of butterflies from oases are more elongated compared with those from the other habitats. This pattern of phenotypic divergence may suggest a recent colonization of the oasis habitats by individuals, which might be accompanied by a rather dispersive behavior. Species distribution modeling suggests a fairly recent reexpansion of the species’ climatic niche starting in the Holocene at about 6000 before present. The combined results indicate a rather recent colonization of the oases by highly mobile individuals from genetically diverse founder populations. The colonization was likely followed by the expansion and persistence of these founder populations under relatively stable environmental conditions. This, together with low rates of gene flow, likely prevented differentiation of populations via drift and led to the maintenance of high genetic diversity.

Anthropogenic habitat fragmentation is one of the major drivers of biodiversity loss. Fragmentation causes the decline of habitat size and limits gene flow between populations in isolated habitat patches. The effects of fragmentation differ depending on the location within the species’ distribution, patch size, and distances among adjacent patches. Small populations, with large geographic distances among them, as well as populations at the distribution edge suffer most severely and often display strong population divergence and reduced genetic diversity (Lesica and Allendorf 1995; Allendorf and Luikart 2007). Reduced genetic diversity in isolated populations has been shown to negatively affect individual fitness and thus population viability (Saccheri et al. 1998; Reed and Frankham 2003; Melbourne and Hastings 2008; Vandewoestijne et al. 2008).

Most research on habitat fragmentation in terrestrial systems has been performed in regions with strong human impact (Keller et al. 2004, 2005). These systems are usually characterized by semipermeable and highly heterogeneous landscapes, which often include stepping stones and corridors that aid in maintaining at least low levels of gene flow (Cook et al. 2002; Hanski and Gaggiotti 2004). In contrast, many naturally fragmented systems consist of habitat patches embedded in a matrix of highly unsuitable environment. Habitat patches in these systems are often separated by much larger geographic distances and have been isolated over long time scales. Although being of particular interest for the study of island biogeography, very few studies on the effects of fragmentation have been performed in naturally fragmented landscapes.

A unique and well-suited model system to analyze naturally fragmented populations is represented by forest species living in oases within deserts. These rare conditions can be found in the oases of the northwestern Sahara desert. These oases are strictly isolated from each other by large areas of stone desert since the final desiccation of the Sahara approximately 6000 years ago (de Noblet-Ducoudré and Prentice 2000). Forest-like habitats within these oases are lined up along several vadis (rivers only temporarily carrying water) extending from the southern High Atlas Mountains of Morocco into the Sahara desert.

Butterfly species of the subfamily Satyrinae are good models to study evolutionary processes, including effects of population isolation. Many insular and geographically restricted endemics are known for this subfamily (Cesaroni et al. 1994; Dapporto 2010; Schmitt and Besold 2010,; Thomson 2011). The Speckled Wood butterfly Pararge aegeria (Linnaeus 1758) is a particularly suitable model species to study the impact of population isolation because it is restricted to shady and moist forest habitats, whereas arid conditions are completely unfavorable (Wiklund and Persson 1983; Shreeve 1984; Merckx et al. 2003; Stevens 2004). Pararge aegeria is widespread over the Western Palaearctic region, ranging from North Africa to Scandinavia and from Iberia to the Urals. The species is common in the Central and Southern European broad-leaved forests, but also occurs in fragmented agricultural landscapes along hedgerows (Dover and Sparks 2000; Merckx et al. 2003; Vandewoestijne and Van Dyck 2010). Populations of the Western Mediterranean area show significant genotypic and phenotypic variation on large as well as regional scales (Weingartner et al. 2006; Dapporto 2010). South of the Atlas Mountains, the butterfly species is strictly confined to oases where it occurs in shady areas with a dense cover of olive trees and date palms (personal observations).

We sampled 14 oasis populations of P. aegeria stretching 250 km along an east–west transect from the Todra Valley to the Guir River vadi in the transitional zone between the High Atlas and the northern Sahara. Additionally, 5 populations from the Atlas Mountains, 2 from Tunisia, and 5 from Europe (3 from France and 2 from Belgium) were analyzed (Table 1). We used 6 microsatellite loci as our genetic marker system to calculate parameters of genetic differentiation, genetic variability and to generate estimates of effective population sizes (Ne). This noncoding marker system is highly polymorphic and fast evolving and hence is well suited for intraspecific analyses of genetic divergence (Selkoe and Toonen 2006). In contrast, the morphological characters (wing shape and genitalia phenotypes) analyzed here are known to be frequently affected by natural (and sexual) selection and therefore might provide information on local selective pressures. Genitalia can evolve relatively independently from environmental factors and thus are often highly correlated with neutral genetic divergence (Cesaroni et al. 1994; Dapporto 2010; Dapporto et al. 2011b; Dincǎ et al. 2011, but see Dapporto et al. 2011b). Genitalia shape is involved in reproductive isolation and therefore can be under strong sexual selection (discussed in Shapiro and Porter 1989). In contrast, wing shape is more likely influenced by natural selection (Berwaerts et al. 2002). Variations in these traits may therefore mirror contrasting behavior and life-histories (Van Dyck et al. 1997; Berwaerts et al. 2002) and are used in this study to detect different past and present selection regimes (Dennis and Shreeve 1989; Joyce et al. 2009). In addition to genetic and morphometric analyses, we employ a modeling approach to understand past and present species distribution (and possible range contraction and expansion events). Specifically, we aim to answer the following questions using a combination of the above-mentioned complementary analyses:

Table 1

Sampling localities of Pararge aegeria. Given are the name of the sampling location, with gps coordinates, and numbers of collected individuals analyzed for microsatellites (µsats), wing shape (W) and genitalia morphology (G), as well as sampling dates

 Name of locality Site Gps coordinates N (µsats) N (W) N (G) Sampling dates 
Oases (Morocco) Gorges du Todra 31°32ʹN; 5°34ʹW 21 15 13 7-3-2009; 14-3-2010 
South Goulmima 31°36ʹN; 4°57ʹW 30 29 24 13-3-2010 
Gaouz 31°41ʹN; 4°56ʹW 25 25 19 13-3-2010 
Taltfraout 31°48ʹN; 4°57ʹW 24 23 18 21-5-2008; 6-3-2010 
Tardirhoust 31°53ʹN; 4°55ʹW 22 20 13 6-3-2010 
Tarda 31°49ʹN; 4°36ʹW 28 25 21 20-5-2008; 7-3-2010 
Errachidia 31°58ʹN; 4°27ʹW 29 30 24 7-3-2010 
Ifri 32°04ʹN; 4°23ʹW 30 30 22 10-3-2010 
Source Bleu de Meski 31°51ʹN; 4°16ʹW 23 17 15 22-5-2008; 7-3-2010 
Boudnib 10 31°56ʹN; 3°36ʹW — 25 15 21-5-2008; 8-3-2010 
Sahli 11 31°57ʹN; 3°32ʹW 21 21 18 21-5-2008; 8-3-2010 
Zrigat 12 31°38ʹN; 4°12ʹW 21 — — 23-5-2008; 10-3-2010 
Ba Touroug 13 31°33ʹN; 4°42ʹW 30 30 28 12-3-2010 
Mellaab 14 31°33ʹN; 4°50ʹW 24 20 34 12-3-2010 
Atlas Mountains (Morocco) Imlil 15 31°12ʹN; 7°58ʹW 26 — — 13-5-2008 
Vallée d’Ourika 16 31°14ʹN; 7°37ʹW 17 — — 14-5-2008 
Cascade d’Ouzout 17 32°04ʹN; 6°40ʹW 11 — — 16-5-2008 
Kasba Tadla 18 32°51ʹN; 5°56ʹW 10 — — 9-5-2008 
Azrou 19 33°27ʹN; 5°13ʹW — 12 18-5-2008 
Tunisia Ain Draham 20 36°43ʹN; 9°11ʹE 20 31 28 5-2010 
Béjà 21 36°46ʹN; 8°41ʹE 20 25 21 5-2010 
Europe Gerompont (Belgium) 22 50°39ʹN; 4°53ʹE 20 20 20 8-2007 
 Meerdaal (Belgium) 23 50°48ʹN; 4°43ʹE 20 — — 8-2007 
 Corgens (France) 24 46°22ʹN; 5°06ʹE 20 — — 8-2007 
 Bois Fevrier (France) 25 46°11ʹN; 4°56ʹE 20 — — 8-2007 
 Le Muy, Var (France) 26 43°28ʹN; 6°34ʹE — — 20 4-2006 
 Name of locality Site Gps coordinates N (µsats) N (W) N (G) Sampling dates 
Oases (Morocco) Gorges du Todra 31°32ʹN; 5°34ʹW 21 15 13 7-3-2009; 14-3-2010 
South Goulmima 31°36ʹN; 4°57ʹW 30 29 24 13-3-2010 
Gaouz 31°41ʹN; 4°56ʹW 25 25 19 13-3-2010 
Taltfraout 31°48ʹN; 4°57ʹW 24 23 18 21-5-2008; 6-3-2010 
Tardirhoust 31°53ʹN; 4°55ʹW 22 20 13 6-3-2010 
Tarda 31°49ʹN; 4°36ʹW 28 25 21 20-5-2008; 7-3-2010 
Errachidia 31°58ʹN; 4°27ʹW 29 30 24 7-3-2010 
Ifri 32°04ʹN; 4°23ʹW 30 30 22 10-3-2010 
Source Bleu de Meski 31°51ʹN; 4°16ʹW 23 17 15 22-5-2008; 7-3-2010 
Boudnib 10 31°56ʹN; 3°36ʹW — 25 15 21-5-2008; 8-3-2010 
Sahli 11 31°57ʹN; 3°32ʹW 21 21 18 21-5-2008; 8-3-2010 
Zrigat 12 31°38ʹN; 4°12ʹW 21 — — 23-5-2008; 10-3-2010 
Ba Touroug 13 31°33ʹN; 4°42ʹW 30 30 28 12-3-2010 
Mellaab 14 31°33ʹN; 4°50ʹW 24 20 34 12-3-2010 
Atlas Mountains (Morocco) Imlil 15 31°12ʹN; 7°58ʹW 26 — — 13-5-2008 
Vallée d’Ourika 16 31°14ʹN; 7°37ʹW 17 — — 14-5-2008 
Cascade d’Ouzout 17 32°04ʹN; 6°40ʹW 11 — — 16-5-2008 
Kasba Tadla 18 32°51ʹN; 5°56ʹW 10 — — 9-5-2008 
Azrou 19 33°27ʹN; 5°13ʹW — 12 18-5-2008 
Tunisia Ain Draham 20 36°43ʹN; 9°11ʹE 20 31 28 5-2010 
Béjà 21 36°46ʹN; 8°41ʹE 20 25 21 5-2010 
Europe Gerompont (Belgium) 22 50°39ʹN; 4°53ʹE 20 20 20 8-2007 
 Meerdaal (Belgium) 23 50°48ʹN; 4°43ʹE 20 — — 8-2007 
 Corgens (France) 24 46°22ʹN; 5°06ʹE 20 — — 8-2007 
 Bois Fevrier (France) 25 46°11ʹN; 4°56ʹE 20 — — 8-2007 
 Le Muy, Var (France) 26 43°28ʹN; 6°34ʹE — — 20 4-2006 
  • Did geographical isolation and variation in ecological conditions lead to the evolution of distinct genetic and morphologic lineages?

  • Are oasis populations genetically and morphologically differentiated and impoverished due to their inherently limited size and isolation?

  • When were the Saharan oases colonized by P. aegeria? Did climatic conditions favor this invasion process?

Materials and Methods

Sampling

A total of 441 P. aegeria individuals were collected in North Africa (Morocco, Tunisia) and Europe between 2006 and 2010. The sampled individuals represent 14 oasis populations, 5 populations from the Atlas Mountains, 2 populations from Tunisia, 2 populations from Belgium, and 3 populations from France (Figure 1). Genetic data were obtained for 24 of these populations. Details about sampling locations are given in Table 1. Immediately after netting, individuals were killed and stored in absolute ethanol or dried in the sun and stored in envelopes until further processing. Details on samples from Europe (excluding Le Muy) are given in Vandewoestijne and Van Dyck (2011).

Figure 1.

Geographic locations of the 26 sampling sites of Pararge aegeria over North Africa and Europe (Morocco, Tunisia, Belgium, France). The detailed inlay shows 14 sampled populations of the oases at the southern edge of the Atlas Mountains. Numbers of localities coincide with Table 1. Elevations more than 500 m asl are shaded in light gray.

Figure 1.

Geographic locations of the 26 sampling sites of Pararge aegeria over North Africa and Europe (Morocco, Tunisia, Belgium, France). The detailed inlay shows 14 sampled populations of the oases at the southern edge of the Atlas Mountains. Numbers of localities coincide with Table 1. Elevations more than 500 m asl are shaded in light gray.

Molecular Analyses

DNA was isolated from the thorax using the DNeasy Blood and Tissue kit (Qiagen, Hilden, Germany) following the manufacturer’s protocol for tissue samples. Six polymorphic microsatellite loci (Helsen et al. 2010) were amplified. These microsatellites were already successfully used previously for P. aegeria (Vandewoestijne and Van Dyck 2010) and show high levels of polymorphism. Most of these loci are characterized by insertions and deletions of single base pairs in the flanking regions; these substitutions of single base pairs lead to deviations from the repeat motif (occurrence of 1-nucleotide steps, supplementary appendix 1 online) as described in Helsen et al. (2010). Sequencing studies showed that such deviations from the microsatellite motif occur at nonnegligible rates in the flanking regions of microsatellites (Angers and Bernatchez 1997; Grimaldi and Crouau-Roy 1997) and were frequently observed for microsatellites in Lepidoptera. This was explained by high similarities in flanking regions between different microsatellites within the same species and/or the lack of conserved flanking regions (Meglécz et al. 2007).

The forward primers were labelled with the fluorescent dyes 6FAM, 6TAMN, HEX (IDTdna) using multiplex PCR following a modified protocol provided by Helsen et al. (2010). Although 8 primer pairs were originally included in multiplex PCR, only 6 amplified successfully for North African samples. Multiplex amplifications were performed with QIAGEN Multiplex PCR kit (Qiagen) using a Veriti Thermal Cycler (Applied Biosystems). PCR conditions were as follows: 95 °C for 15min, followed by 32 cycles (94 °C for 30 s, 60 °C for 90 s, 72 °C for 90 s) and 72 °C for 30min, terminating at 25 °C. Amplified fragments were visualized with an ABI PRISM 3100 Genetic Analyser (Applied Biosystems). Allele sizes were scored against the internal standard ROX-400SD (Applied Biosystems) using GeneMapper 3.5 (Applied Biosystems).

The program Micro-Checker (Van Oosterhout et al. 2004) was used to test for distortion of the genotypic data as a result of stutter bands, large allele dropouts, or null alleles (Selkoe and Toonen 2006). Allelic richness (AR) was calculated using FSTAT (Goudet 1995). Nonhierarchical and hierarchical analyses of genetic variance (analyses of molecular variance [AMOVAs]), calculations of observed (Ho) and expected (He) heterozygosities, tests for Hardy–Weinberg equilibrium, linkage disequilibrium, and isolation-by-distance analyses were performed with Arlequin 3.1 (Excoffier et al. 2005).

AMOVAs were carried out using 2 different approaches: the conventional F-statistics based on allele frequencies and the infinite alleles model and the microsatellite-specific R-statistics based on allele lengths and the stepwise mutation model (Slatkin 1995). The latter assumes that alleles similar in length are more likely to share a more recent common ancestor. Both approaches comprised 3 hierarchical levels: within populations, among populations within regions, and among regions. Differentiation between populations was estimated as pairwise RST values.

Isolation-by-distance was tested using Mantel tests performed on matrices of pairwise geographic distances (km) and linearized pairwise RST values, that is, RST/(1 − RST). Pairwise geographic distances were calculated with GeoDist (Heidenreich A, unpublished data). Mantel tests were performed using the isolation-by-distance web service version 3.23 (Jensen et al. 2005); hereby, isolation-by-distance tests were performed on different geographical units to adjust for swamping effects large geographic distances can have on the method. Mantel tests were performed using 30 000 randomizations. Isolation by distance is one of the most common patterns observed in natural populations not disrupted by natural or anthropogenic barriers and hence can be seen as a null model in studies of population divergence. If genetic differences among specific populations are higher than the isolation-by-distance model predicts barriers to gene flow often can be invoked. If genetic differentiation is less pronounced than predicted, the effects of genetic drift are mediated by large population sizes or gene flow.

Analyses of population structure were performed using a Bayesian approach as implemented in STRUCTURE 3.1 (Pritchard et al. 2000). For each run, burn-in and simulation lengths were 100 000 and 300 000, respectively. The calculations were run under the admixture model with correlated gene frequencies. To obtain information about recent bottlenecks, each sample was tested using the Wilcoxon test as implemented in BOTTLENECK (Cornuet and Luikart 1997). This test is based on the expectation that a recent reduction in effective population size leads to a reduction in both allele numbers and gene diversity. As allele number decreases faster than gene diversity in a recently bottlenecked population, the observed gene diversity is higher than the expected equilibrium gene diversity (Luikart and Cornuet 1998).

Additionally, effective population sizes Ne were estimated with ONeSAMP (Tallmon et al. 2008). Ne is an important population genetic parameter because it indicates a population’s potential to evolve via random processes. Generally, populations with small Ne are more susceptible to genetic drift; consequently geographically separated populations with small Ne will differentiate faster and loose rare alleles quicker than larger populations.

Wing Shape Morphometrics

Digital images of the ventral side of the forewings of 386 P. aegeria individuals from 16 populations (13 oases, 2 from Tunisia, and 1 from Belgium; Table 1) were taken with a digital camera (Canon M1, 50-mm lens) under standardized light conditions (Häuser et al. 2005). Seven landmarks on wing veins (crossing points or at the wing margin) of the forewing were used to characterize the wing shape and the relation between wing width (distance between landmarks 5 and 7) and wing length (distance between landmarks 1 and 5) (Figure 2). Landmarks were digitized with the program TpsDig 2.16 (Rohlf 2010a). The order of specimens was randomized using TpsUtil 1.46 (Rohlf 2010b). We applied a generalized least square Procrustes fit analysis using TpsRelw 1.49 (Rohlf 2010c) to obtain relative warps (similar to principle components [PCs]), which can be used as shape variables. The Levene’s test of homogeneity of variances was used to verify if PC scores showed different variances among populations. A full cross-validation test was conducted to identify variances among oasis and non-oasis populations using U tests. To evaluate whether different populations can be distinguished on the basis of different wing shape, we performed a discriminant function analysis of PC scores using population membership as a priori grouping variable. Using PC scores in discriminant analyses allows for the detection of shape variation among groups. A full cross-validation test was applied to evaluate the power of the discriminant functions and the possibility to blindly attribute each specimen to the original group. Statistic analyses were conducted with the SPSS 17 software package (IBM, Somers, NY).

Figure 2.

(a) Forewing of Pararge aegeria with selected landmarks used to calculate total wing shape, wing length (between 5 and 7), and wing width (between 1 and 5). (b) Valva and (c) aedeagus with selected landmarks (black dots) and sliding semilandmarks (white dots).

Figure 2.

(a) Forewing of Pararge aegeria with selected landmarks used to calculate total wing shape, wing length (between 5 and 7), and wing width (between 1 and 5). (b) Valva and (c) aedeagus with selected landmarks (black dots) and sliding semilandmarks (white dots).

Wing Coloration Analyses

The same digital pictures as used for wing shape analyses were transformed, splitting the photographs into 3 spectral colors (red, green, and blue) using the program Matlab R2007a 7.4 (The Mathworks Inc.). These 3 colors were merged into a mean gray value. A wing polygon bordered by the landmarks 1, 5, 6, and 7 was used to curtail a defined and identical wing area for each individual to calculate the degree of melanization. The total sum of pixels within this polygon was counted for each picture. A threshold value was selected being sensitive enough to classify certain proportions of the pixels being white and black. All pixels falling below this threshold were recorded as black, all others as white. The ratio divided by the total number of pixels (in %) was analyzed for the wing of each individual. Furthermore, we calculated the same ratio (in %) using each spectral range (red, green, and blue) separately. We performed discriminant function analysis by using the percentages of each threshold value for each color as predictors and population membership as a priori grouping variable.

Genitalia Morphology

The 9th abdominal segment of male butterflies is divided into a dorsal tegumen and a ventral vinculum forming a structure for the attachment of genital parts and a pair of lateral clasping organs (valvae). Males also have a median tubular organ (called the aedeagus), which is directly involved in the insemination of females. Lateral sections of the valvae and the aedeagus showed significant differences among Mediterranean populations of P. aegeria in previous analyses and hence are suitable and sensitive characters for detecting intraspecific differentiation (Dapporto 2010). Male genitalia from 359 individuals belonging to 18 populations (13 oasis populations, 2 from Tunisia, each one from the Atlas Mountains, Belgium, and France; Table 1) were dissected using standard procedures; abdomina were boiled in 10% caustic potash. Genitalia were cleaned and the left valva and the aedeagus removed. Tegumen and valva were mounted on Euparal between microscope slides and cover slips. Aedeagi were included in dimethyl hydantoin formaldehyde resin without using cover slips in order to avoid possible deformations. Lateral views of aegeadi and valvae were photographed with a digital camera mounted on a dissecting microscope. We applied a combination of landmarks and sliding semilandmarks (Bookstein 1991). Three landmarks on the outlines of valvae and 5 landmarks defining the shape of aedeagi were considered (Figure 2b). Landmark choice was based on easy and secure identification of homologous structures across individuals (type II and type III landmarks). An additional set of 38 points for valvae and 23 points for aedeagi were defined as sliding semilandmarks using TpsUtil. Although landmarks represent fixed homologue points, sliding semilandmarks provide an opportunity to digitize outlines when few landmarks are available. Each point is allowed to slide along a tangential direction so that the contour is homologue, whereas each individual point does not need to be (Perez et al. 2006). Landmarks and sliding semilandmarks were digitized with TpsDig and the order of specimens randomized using TpsUtil. We applied a generalized least square Procrustes fit analysis using TpsRelw to obtain relative warps (PCs). Only PCs explaining more than 1% of the total variance were considered in successive analyses (Dapporto 2010). We performed discriminant function analysis by using the aedeagus and valva PCs as predictors and population membership as a priori grouping variable.

Species Distribution Modeling

In order to assess the current and past potential distribution of P. aegeria, we obtained a comprehensive set of species records from our own field surveys and from scientific collections linked to the Global Biodiversity Information Facility (www.gbif.org). Including only specimen-based information with a spatial precision of coordinates meeting the resolution of the climate data (see below), we excluded all raster data. As the remaining records showed a high degree of spatial clustering caused by locally varying sampling efforts rather than by actual preference of the species to specific habitats, we applied a spatial filter leaving 2753 randomly selected records with a minimum distance to the nearest neighbor of 0.092° (median = 0.425°). The minimum distance corresponds to a separation of approximately 2 grid cells of the climate data.

A set of 19 bioclimatic variables (Busby 1991; Beaumont et al. 2005), representing average conditions between 1950 and 2000 with a spatial resolution of 0.041°, was obtained from the worldclim database (www.worldclim.org; Hijmans et al. 2005a). In order to assess the species’ past potential distribution, we downloaded 2 different sets of bioclimatic variables for the last glacial maximum (LGM) 21 000 before present (BP) derived from the Community Climate System Model version 3 (CCSM; Otto-Bliesner et al. 2006) and the Model for Interdisciplinary Research on Climate version 3.2 (MIROC; Hasumi and Emori 2004) and spatially downscaled by R.J. Hijmans for worldclim. Using the same downscaling procedure as for the LGM dataset described by Peterson and Nyári (2008), we downscaled further CCSM and MIROC simulations for 6000 BP using the respective original files downloaded from the Paleoclimatic Modelling Intercomparison Project 2 (http://pimp2.lsce.ipsl.fr/).

A variety of different modeling algorithms for Species Distribution Model (SDM) applications is currently available. Herein, we selected a simple profiling technique termed BIOCLIM (Busby 1991) assessing a multivariate rectangular environmental space occupied by the target species. BIOCLIM does not allow the incorporation of complex interactions between variables, but is at the same time not sensitive to high intercorrelations among predictors. Furthermore, BIOCLIM can be used to visualize the spatial distribution and position of areas within the bioclimatic envelope spanned by the species records in an easily interpretable form.

We believe that BIOCLIM may outperform other more complex algorithms in this specific case because the study species virtually occurs all over the study area. Presence/absence or presence/pseudo-absence methods are likely to fail in effectively discriminating between environmental conditions at the species records and random background conditions. BIOCLIM models were computed in Cran R version 2.14 using the dismo package (Hijmans et al. 2011). Past bioclimatic conditions unavailable within the study area under current conditions were quantified using Multivariate Environmental Similarity Surfaces (MESS; Elith et al. 2010).

Results

Molecular Data

The 6 microsatellite loci showed only marginal evidence of null alleles or large allele dropouts for the populations analyzed (Locus 11: 3 of 21; Locus 16: 2 of 21; Locus 17: 3 of 21; Locus 19: 2 of 21; Locus 2: 4 of 21; and Locus 7: 3 of 21). All microsatellite loci were polymorphic throughout all 21 populations from North Africa. Locus-specific allele frequencies for each population are given in the Supplementary Appendix 1 online. The mean value of AR was 6.75 (± 0.47 standard deviation [SD]) and 7.72 (± 0.31 SD) calculated for a reduced set of 4 microsatellites (data for 2 of the 6 loci were unavailable for the European samples). Measures of heterozygosity were also high: He = 79.7% (± 6.3 SD) (4 loci: 84.3% ± 1.6 SD) and Ho = 61.6% (± 0.9 SD) (4 loci: 60.1% ± 0.8 SD) (Table 2). Population-based diversity estimates are given in the Supplementary Appendix 3 online. Heterozygosity deficiency was detected for the following populations: Gaouz, Taltfraout, Ifri, Sahli, Zrigat, Corgens, and Bois Fevrier (P < 0.05). No heterozygostiy excess was detected. Calculation of effective population sizes yielded fairly similar estimates for all populations ranging from 97 to 156 individuals (values of each population are given in Supplementary Appendix 3 online). Neither the calculations of potential population bottlenecks nor the estimates of effective population sizes revealed significant differences between oases and other populations (Mann–Whitney U tests, all P > 0.05). Genetic diversity was not significantly different among oases, the Atlas Mountains and Tunisia (Kruskal–Wallis test, P > 0.05). No difference in genetic diversity was found between North Africa and Europe either (Kruskal–Wallis test, P > 0.05).

Table 2

Means and SDs of parameters of genetic diversity analyzed for geographic groups of Pararge aegeria: allelic richness (AR), percentage of expected heterozygosity (He), and percentage of observed heterozygostiy (Ho). Furthermore, the ratio of wing length to wing width is given as parameter “wing design.” Calculations were performed using 4 and 6 microsatellite loci. Differences between North Africa and Europe were tested by using Kruskal–Wallis ANOVA

 Oases Atlas Mountains Tunisia North Africa Belgium France P value 
AR_4  7.78 (±0.42) 7.38 (±0.83)  7.99 (±0.38)  7.72 (±0.31)  6.88 (±0.14)  7.73 (±1.08) ns 
AR_6  6.85 (±0.24) 6.24 (±0.38)  7.17 (±0.11)  6.75 (±0.47) — — ns 
He_4 84.6 (±1.6) 82.6 (±3.9) 85.8 (±2.5) 84.3 (±1.6) 76.4 (±0.3) 84.1 (±6.7) ns 
He_6 73.7 (±14.3) 79.3 (±11.5) 86.3 (±1.2) 79.7 (±6.3) — — ns 
Ho_4 59.2 (±18.1) 60.7 (±6.0) 60.2 (±3.1) 60.1 (±0.8) 72.7 (±1.5) 78.0 (±0.8) ns 
Ho_6 60.7 (±18.4) 61.9 (±8.8) 62.36 (±4.86) 61.6 (±0.9) — — ns 
Wing design  1.23 (±0.03) —  1.18 (±1.19)  1.22 (±0.03)  1.17 (±0.04) — <0.05 
 Oases Atlas Mountains Tunisia North Africa Belgium France P value 
AR_4  7.78 (±0.42) 7.38 (±0.83)  7.99 (±0.38)  7.72 (±0.31)  6.88 (±0.14)  7.73 (±1.08) ns 
AR_6  6.85 (±0.24) 6.24 (±0.38)  7.17 (±0.11)  6.75 (±0.47) — — ns 
He_4 84.6 (±1.6) 82.6 (±3.9) 85.8 (±2.5) 84.3 (±1.6) 76.4 (±0.3) 84.1 (±6.7) ns 
He_6 73.7 (±14.3) 79.3 (±11.5) 86.3 (±1.2) 79.7 (±6.3) — — ns 
Ho_4 59.2 (±18.1) 60.7 (±6.0) 60.2 (±3.1) 60.1 (±0.8) 72.7 (±1.5) 78.0 (±0.8) ns 
Ho_6 60.7 (±18.4) 61.9 (±8.8) 62.36 (±4.86) 61.6 (±0.9) — — ns 
Wing design  1.23 (±0.03) —  1.18 (±1.19)  1.22 (±0.03)  1.17 (±0.04) — <0.05 

ns, not significant.

Nonhierarchical AMOVAs detected a moderate genetic variance within individuals throughout all populations analyzed (12.59, RIS: 0.171, P < 0.0001). The major proportion of molecular variance was detected among populations (17.33, RST: 0.191, P < 0.0001). Analyses for each geographic group separately (Oases, Atlas Mountains, Tunisia, Morocco and Tunisia, Belgium, France, and Central Europe) yielded no significant genetic differentiation among the North African populations, among any of the oasis populations, or among the populations from the Atlas Mountains and between Morocco and Tunisia (Table 3). Results obtained from F-statistics are given in the Supplementary Appendix 4 online.

Table 3

Nonhierarchical analyses of molecular variance on Pararge aegeria populations. Hereby, 3 levels are taken into account: genetic variance among populations, among individuals within populations, and within individuals. Given are values of genetic variance with respective R-statistics metric (in parenthesis). For the AMOVA analyses, we predefined the following regions: Oases, Atlas Mountains, Tunisia, Morocco + Tunisia, Belgium, France, and Central Europe

 Among populations (RSTAmong individuals within populations (RISWithin individuals 
Oases 0.30590 (0.00565) 4.04940 (0.07520) 49.79808 
Atlas Mountains −0.48498 (−0.01129) 6.95335 (0.16010) 36.47692 
Tunisia −0.74915 (−0.01229) 6.66813 (0.10804) 55.05263 
Morocco + Tunisia 0.24942 (0.00469) 4.76885 (0.09004) 48.19277 
Belgium −0.82953 (−0.00872) 14.88065 (0.15516) 81.02381 
France 7.96330 (0.07297**) 3.10781 (0.03072) 98.05319 
Central Europe 2.71967 (0.02682) 8.65922 (0.08775) 90.01685 
All 17.32697 (0.19075***) 12.58912 (0.17126***) 60.92063 
 Among populations (RSTAmong individuals within populations (RISWithin individuals 
Oases 0.30590 (0.00565) 4.04940 (0.07520) 49.79808 
Atlas Mountains −0.48498 (−0.01129) 6.95335 (0.16010) 36.47692 
Tunisia −0.74915 (−0.01229) 6.66813 (0.10804) 55.05263 
Morocco + Tunisia 0.24942 (0.00469) 4.76885 (0.09004) 48.19277 
Belgium −0.82953 (−0.00872) 14.88065 (0.15516) 81.02381 
France 7.96330 (0.07297**) 3.10781 (0.03072) 98.05319 
Central Europe 2.71967 (0.02682) 8.65922 (0.08775) 90.01685 
All 17.32697 (0.19075***) 12.58912 (0.17126***) 60.92063 

*P < 0.05; **P < 0.01; ***P < 0.001.

The hierarchical analysis of molecular variance indicated a shallow split between Morocco and Tunisia (RCT: 0.025, P < 0.05). However, considerable differentiation was detected between European and African populations (54.95, RCT: 0.426, P < 0.0001) (Table 4). Respective values obtained for F-statistics are fairly similar and are provided in the Supplementary Appendix 5 online.

Table 4

Hierarchical variance analyses of Pararge aegeria populations. Hereby, 4 levels of genetic variance were taken into account: genetic variance found among predefined groups, among populations within groups, among individuals within populations, and genetic variance found within individuals. Predefined geographic groups are in accordance with the groups given in Table 3. Given are variance values with respective R-statistics (in parenthesis)

Groups Among groups (RCTAmong populations within groups (RSCAmong individuals within populations (RISWithin individuals 
Oases vs. Atlas Mountains −0.73165 (0.00395) 0.26710 (−0.01093) 12.89644 (0.19132***) 54.51061 
Morocco vs. Tunisia 1.76532 (0.02528*) −0.02416 (−0.00035) 13.41835 (0.19704***) 54.68072 
North Africa vs. Central Europe 54.95278 (0.42576***) 0.60815 (0.00821***) 12.58912 (0.17126***) 60.92063 
Belgium vs. France −1.61731 (−0.01604) 3.79824 (0.03707*) 8.65922 (0.08775*) 90.01685 
Groups Among groups (RCTAmong populations within groups (RSCAmong individuals within populations (RISWithin individuals 
Oases vs. Atlas Mountains −0.73165 (0.00395) 0.26710 (−0.01093) 12.89644 (0.19132***) 54.51061 
Morocco vs. Tunisia 1.76532 (0.02528*) −0.02416 (−0.00035) 13.41835 (0.19704***) 54.68072 
North Africa vs. Central Europe 54.95278 (0.42576***) 0.60815 (0.00821***) 12.58912 (0.17126***) 60.92063 
Belgium vs. France −1.61731 (−0.01604) 3.79824 (0.03707*) 8.65922 (0.08775*) 90.01685 

*P < 0.05; **P < 0.01; ***P < 0.001.

Population- (neighbor joining tree, not shown) and individual-based (Figure 3) genetic analyses supported a significant split between Europe and North Africa (Figure 3a) and absence of differentiation within North Africa (Figure 3b, 3c). No significant isolation by distance was detected for North Africa excluding the oases (N = 7 populations), z = 199.0357, r = 0.3400, P = 0.1455, including the oases (N = 20), z = 1240.9531, r = 0.221, P = 0.082, Morocco including the oases (N = 18), z = −89.7969, r = −0.2281, P = 0.9146, the oases populations themselves (N = 13), z = 88.9228, r = 0.3088, P = 0.0700, or Europe (N = 4), z = 44.1515, r = −0.0723, P = 0.7104. However, significant results were obtained for Morocco excluding the oases populations (N = 5), z = −1.4181, r = 0.5748, P = 0.0084 and for the complete set including all samples (N = 24), z = 77696.6009, r = 0.8098, P < 0.0000 (Table 5). Respective graphs are given in Supplementary Appendix 2 online.

Table 5

Results obtained from isolation-by-distance analyses performed separately for different geographic groups

Group Number of populations Z r P 
All populations 24 77696.6009 0.8098 <0.0001 
North Africa (with oases) 20 1240.9531 0.2201 0.0820 
North Africa (without oases) 199.0357 0.3400 0.1455 
Morocco (with oases) 18 −89.7969 −0.2281 0.9146 
Morocco (without oases) −1.4182 0.5748 0.0084 
Oases 13 88.9228 0.3088 0.0700 
Europe 44.1515 −0.0723 0.7104 
Group Number of populations Z r P 
All populations 24 77696.6009 0.8098 <0.0001 
North Africa (with oases) 20 1240.9531 0.2201 0.0820 
North Africa (without oases) 199.0357 0.3400 0.1455 
Morocco (with oases) 18 −89.7969 −0.2281 0.9146 
Morocco (without oases) −1.4182 0.5748 0.0084 
Oases 13 88.9228 0.3088 0.0700 
Europe 44.1515 −0.0723 0.7104 
Figure 3.

Bayesian analysis performed on populations of Pararge aegeria using the STRUCTURE software (Pritchard et al. 2000). (a) Performing the analysis for K = 2 including all samples available (dividing Europe from North Africa), (b) performing the analyses for K = 2 including all individuals collected in North Africa (western vs. eastern), and (c) performing the analysis for K = 2 with all samples from the western part of North Africa (Oasis populations vs. populations from the Atlas Mountains).

Figure 3.

Bayesian analysis performed on populations of Pararge aegeria using the STRUCTURE software (Pritchard et al. 2000). (a) Performing the analysis for K = 2 including all samples available (dividing Europe from North Africa), (b) performing the analyses for K = 2 including all individuals collected in North Africa (western vs. eastern), and (c) performing the analysis for K = 2 with all samples from the western part of North Africa (Oasis populations vs. populations from the Atlas Mountains).

Wing Morphometrics

The morphometric analyses of the wings resulted in 12 PCs explaining a cumulative variance of 100%. Levene’s test revealed that only PC1 contained slightly significant different variances separating populations (Levene statistic = 1.949, P = 0.014). However, oasis populations did not show lower variances compared with the remaining Moroccan and Tunisian populations (Mann–Whitney U = 23.0, P = 0.578). Discriminant analysis extracted 6 functions showing a significant Wilks’s lambda (P < 0.05). These functions cumulatively explained 86.9% of the variance. The cross-validation test only assigned a low percentage of individuals (29.5%) to their original groups. Oasis populations showed particularly low percentages of correctly assigned individuals (17.4%) compared with the Atlas and Tunisian populations (35.7%) and the Belgian population (74.6%). Discriminant functions 1 and 2 explained the largest fraction of variance (52.3%, Wilks’s lambda = 0.36, P < 0.001; and 9.5%, Wilks’s lambda = 0.127, P < 0.001) and were mainly correlated with PC1 and PC5. Visual inspection of the plotted scores for these variables and thin-plate splines revealed that European and the Atlas and Tunisian specimens have lower PC1 values with respect to most oasis population. Thin-plate spline grids reveal that PC1 is linked to wing elongation, with oasis populations showing more elongated wings (Figure 4a). The ratio of wing width to wing length shows slight evidence that individuals living in oases have more elongated wings (oases, ratio wing length/wing width, 1.23 [± 0.03]) if compared with individuals from Europe (1.17 [± 0.04]) or Tunisia (1.18 [± 1.19]) (P < 0.05) (Table 2). Population-specific values are given in the Supplementary Appendix 3 online.

Figure 4.

Graphic representation of geometric morphometric variables mostly correlated with discriminant functions among areas. Variations in shape along both axes are shown in thin-plate spline deformation grids. (a) Forewing: the dashed lines only represent an indication for wing shape, but has no morphometric values; (b) male genitalia: numbered squares represent means ± SD values, populations are numbered as in Table 1, black squares, European populations; gray squares, non-oasis population from North Africa; white squares, oasis populations.

Figure 4.

Graphic representation of geometric morphometric variables mostly correlated with discriminant functions among areas. Variations in shape along both axes are shown in thin-plate spline deformation grids. (a) Forewing: the dashed lines only represent an indication for wing shape, but has no morphometric values; (b) male genitalia: numbered squares represent means ± SD values, populations are numbered as in Table 1, black squares, European populations; gray squares, non-oasis population from North Africa; white squares, oasis populations.

Wing Coloration

By using 20 different threshold values for each color, we obtained 60 color variables; 51 of them showed variances different from zero. These variables were included in a discriminant function analysis, which extracted 3 functions showing a significant Wilks’s lambda (P < 0.05). These functions cumulatively explained 79.0% of the variance. The cross-validation test assigned only a low percentage of individuals (20.1%) to their original group. The oasis and the other North African populations showed low percentages of correctly assigned individuals (15.8% and 21.0%); however, the oasis population 6 showed a rather high correct assignment (52.2%). As for wing shape, the Belgium population showed a high correct assignment value (88.2%). Functions 1 and 2 explained the largest fraction of variance (48.1%, Wilks’s lambda = 0.23, P < 0.001; and 19.6%, Wilks’s lambda = 0.426, P < 0.001). The inspection of values from variables entered in the discriminant analysis revealed that individuals from Europe tend to be darker compared with individuals from North Africa.

Genitalia Morphometrics

The morphometric analyses of aedeagi and valvae resulted in 14 and 11 PCs, respectively, explaining more than 1% of the total variance, resulting in a cumulative variance of 92.5% and 93.8%. Levene’s test revealed that only PC1 and PC11 showed slightly significant different variances between populations (PC1 Levene statistic = 1.859, P = 0.024; PC11 Levene statistic = 1.743, P = 0.042). However, oasis populations did not show lower variances compared with the Atlas and Tunisian populations (Mann–Whitney U = 23.0, P = 0.578). Discriminant analysis extracted 3 functions showing a significant Wilks’s lambda (P < 0.05). These functions cumulatively explained 80.1% of the variance. The cross-validation test assigned only a low percentage of individuals (18.3%) to their original group. Oasis populations showed particularly low percentages of correctly assigned individuals (10.5±6.7% SD) compared with non-oasis populations of North Africa (35.8%) and the European populations (94.1%). Functions 1 and 2, represented in Figure 4b, explained the largest fraction of variance (48.1%, Wilks’s lambda = 0.23, P < 0.001; and 19.6%, Wilks’s lambda = 0.426, P < 0.001). Aedeagus PC1 and PC2 were mostly correlated with the first 2 discriminant functions. Visual inspection of the plotted PC scores and thin-plate splines distinguished European from North African populations, mainly by the presence of a larger cornutus in the latter ones. Within the North African populations, no distinction was found between populations from oases, Morocco, and Tunisia (Figure 4b).

Species Distribution Modeling

The BIOCLIM SDM developed under current conditions correctly identified major parts of Europe and North Africa as adequately situated within the bioclimatic envelope defined by the species records (Figure 5). Current distribution positions the oases populations at the southernmost margin of the species’ distribution range. The northern part of North Africa yielded permanently suitable climatic conditions over the past 21000 years for P. aegeria. The species was modeled to be widely distributed over major parts of the North African region with a southern distribution limit at the northern edge of the Sahara desert during the Pleistocene. The model suggested that during the Holocene (6000 BP) the species retracted northwards. Under this scenario P. aegeria disappeared from the Saharan oases, but likely persisted in the northern part of North Africa, where climatic conditions remained stable. Starting at about 6000 years BP, rising humidity at the northern edge of the Sahara has resulted in the southwards expansion of P. aegeria and the recolonization of the oases studied here.

Figure 5.

Spatial distribution of the bioclimatic envelope of Pararge aegeria derived from a BIOCLIM SDM under current and past conditions (6000 and 21 000 years BP) assuming 2 different climate scenarios (CCSM and MIROC). Areas with bioclimatic conditions during the past, which are currently not realized within the shown geographic extent, are indicated as extrapolation areas (MESS).

Figure 5.

Spatial distribution of the bioclimatic envelope of Pararge aegeria derived from a BIOCLIM SDM under current and past conditions (6000 and 21 000 years BP) assuming 2 different climate scenarios (CCSM and MIROC). Areas with bioclimatic conditions during the past, which are currently not realized within the shown geographic extent, are indicated as extrapolation areas (MESS).

Discussion

Given the ecologic characteristics of P. aegeria, we expected to detect significant genetic and morphologic differentiation among oasis populations as a result of strong isolation of habitat patches in the unfavorable matrix of the Saharan desert. However, our findings strongly contrasted with these initial predictions: The only major split found in our data was detected between North African and European populations; isolated oasis populations maintain similar intraspecific genetic variability compared with the interconnected populations from the Atlas Mountains, Tunisia, and Europe and no significant genetic differentiation was found among the oasis populations. In addition no differentiation was detected over the entire North African region. In contrast to the genetic homogeneity and the lack of differentiation in wing color and genitalia morphology among the oasis populations, wing design differed from all other geographic groups. Indeed, oases butterflies are characterized by a more elongated wing design. The 2 unexpected findings: 1) the lack of differentiation among isolated oases populations and 2) the high level of genetic diversity within oasis populations are discussed in further detail below.

Lack of Intraspecific Differentiation

Our data show strong genetic and morphologic cohesiveness of populations from both oases and even across all of Morocco and Tunisia. Two nonexclusive scenarios can explain the lack of differentiation for the isolated oasis populations: 1) a recent colonization from interconnected source populations from the North and hence a limited divergence time and 2) ongoing gene flow among local oasis populations. In addition, the oases provide relatively stable habitat conditions limiting population size fluctuations and therefore the influence of bottlenecks and similar random processes.

The SDMs suggest that P. aegeria was widely distributed over the North African region during the LGM (21 000 BP). The species presumably disappeared from its southernmost distribution range during the Holocene (6000 BP), but subsequently reexpanded southwards leading to the recolonization of the Saharan oases (Figure 5). Hence, the models indicate that the oases were recently colonized from more northern areas (i.e., Atlas Mountains).

Interestingly, oases populations, although being similar in their genetic makeup and in the studied coloration and genitalia phenotypes, show a slightly elongated wing design compared with individuals from Europe and Tunisia, which distinguishes them from the remaining studied populations. Previous studies on the wing design in P. aegeria revealed that elongated wings can be associated with increased dispersal ability (   Van Dyck et al. 1997; Berwaerts et al. 2002  ). Therefore, the detected phenotypic uniqueness in wing design of oases populations in this butterfly species might be caused by selection towards increased dispersal behavior most probably during the colonization process of the oases habitats, when only the most mobile individuals were probably capable of reaching these habitat islands.

Population genetic and phylogeographic studies on P. aegeria support the capability of long-distance dispersal in the wake of habitat expansion, as for example, shown for the colonization of Europe from North Africa (Weingartner et al. 2006) and the rapid colonization of the agricultural landscape in some parts of Europe (Vandewoestijne and Van Dyck 2010). The high dispersal ability and the results of species distribution modeling support the hypothesis of a fairly recently colonization of the oases.

Similar to the genetic data, wing coloration and genitalia morphology showed no significant differentiation, but relatively high within-population variation across North Africa. This was unexpected given the diverse ecological environments the populations were found in (e.g., desert oases vs. mountain populations). However, these findings are congruent with previous studies on P. aegeria in Europe revealing high variability of this species for both morphological and life-history traits (Nylin et al. 1989; Van Dyck and Wiklund 2002; Kemp et al. 2006).

The lack of differentiation in phenotypic characters might be explained by similar selective pressures due to high similarity in climatic and ecological conditions across oasis habitats. Yet, one would expect that neutral characters would diverge in isolated populations as a result of drift, especially when considering the relatively small effective population sizes in oases (97–156 individuals).

It appears likely that the lack of genetic differentiation is the result of a fairly recent invasion of the oases by a large number of founder individuals. The lack of isolation- by-distance for the oases populations further adds evidence for this scenario. However, frequent migration of P. aegeria individuals among oases through passive migration events (e.g., caused by sand storms) is rather unlikely, given the extremely unfavorable climatic conditions outside the oases and the large distances among oasis habitats. Yet, only few migrants are needed to prevent differentiation when strongly divergent selective pressures are missing. Additionally, similar environmental conditions across the oases since the recolonization might have prevented any morphological differentiation of these isolated populations.

The lack of genetic differentiation between these populations is unlikely to be the result of an insufficient resolution of the selected microsatellite markers; the same loci showed low but significant genetic differentiation along a latitudinal transect in Europe and even at the landscape level (Vandewoestijne and Van Dyck 2010). Furthermore, data obtained with the same set of microsatellites also indicate the presence of 4 genetic clusters within the P. aegeria species complex: the first cluster is formed by the North European subspecies P. aegeria tircis, whereas the more southern P. a. aegeria is divided into 3 lineages: North Africa, Iberia, and other southern European regions (unpublished data). The southern European and North African lineages have been supported in this study.

High Intraspecific Variability

Small areas located at the margin of a species distribution often do not host permanent populations, but represent sinks within the metapopulation network (Dennis et al. 2010). In such cases, however, ephemeral populations are usually reduced in size and are characterized by low genetic variability as a result of repeated founder events. Conversely, genetic diversity within oasis populations was high and similar compared with populations from the Atlas Mountains and Tunisia. Furthermore, the genetic diversity found in the populations of North Africa is similar to that detected in European populations (Vandewoestijne and Van Dyck 2010), and relatively high compared with other butterfly species (Meglécz and Solignac 1998; Nève and Meglécz 2000; Meglécz et al. 2007; Habel et al. 2008; Sinama et al. 2011). Estimates of effective population sizes showed similar population sizes across all oases. Furthermore, no significant population size differences were detected between oases and mountain populations; also, no recent bottlenecks were detected for any of the oasis populations. Similarly, a study by Shaibi and Moritz (2010) revealed relatively high genetic diversity (and lack of genetic differentiation) in isolated oasis populations of the honey bee Apis mellifera. Although additional factors such as introgression from domestic stocks altered the genetic patterns in some oases, the authors concluded that oases populations generally were large enough to ensure local survival. A similar genetic pattern was found for desert populations of the green toad, Bufo boulengeri (Nevo et al. 1975). Here, the authors suggested that the high genetic diversity is an adaptive strategy in heterogeneous environments. Overall, persisting stable environmental conditions in oases appear to provide good conditions for the maintenance of viable populations and high genetic diversity despite their isolation from other populations.

Conclusion

During the late Pleistocene, P. aegeria occurred over major parts of North Africa including the northern edge of the Sahara desert. The species retracted northwards during the early Holocene but reexpanded southwards in later stages of the Holocene. In the northern part of North Africa, P. aegeria persisted during phases of northward retractions and occurred in interconnected populations. This resulted in a lack of genetic differentiation between the western and eastern margins of the North African distribution range, which still persists today. The oases on the southern margin of the distribution were presumably colonized fairly recently by immigrants showing a rather high dispersal behavior, being mirrored in the slightly elongated wing design found for individuals of oases populations, but not in populations of Europe or Tunisia. Habitat stability over the last millennia probably prevented strong population fluctuations and concomitant loss of genetic diversity. This, together with low rates of ongoing gene flow likely prevented the divergence of local populations.

Supplementary Material

Supplementary material can be found at http://www.jhered.oxfordjournals.org/.

Funding

Fonds National de la Recherche Luxembourg (FNR); German Academic Exchange Service (DAAD); Natural History Museum Luxembourg.

Acknowledgments

J.C.H. was a DAAD postdoctoral fellow and S.V. was a FRS-FNRS postdoctoral fellow. We thank Joachim Besold for statistical help.

References

Allendorf
FW
Luikart
G.
2007
Conservation and the genetics of populations
 .
Oxford
:
Blackwell Publishing
..
Angers
B
Bernatchez
L.
1997
.
Complex evolution of a salmonid micro satellite locus and its consequences in inferring allelic divergence from size information.
Mol Biol Evol
 .
14
:
230
238
.
Beaumont
LJ
Hughes
L
Poulsen
M
.
2005
.
Predicting species distributions: use of climatic parameters in BIOCLIM and its impact on predictions of species’ current and future distributions
.
Ecol Model.
 
186
:
250
269
.
Berwaerts
K
Van Dyck
H
Aerts
P.
2002
.
Does flight morphology relate to flight performance? An experimental test with the butterfly Pararge aegeria
.
Func Ecol.
 
16
:
484
491
.
Bookstein
FL.
1991
Morphometric tools for landmark data: geometry and biology
 .
Cambridge (UK), New York
:
Cambridge University Press
..
Busby
JR.
1991
.
BIOCLIM - a bioclimatic analysis and prediction system
. In:
Margules
CR
Austin
MP
, editors.
Nature conservation: cost effective biological surveys and data analysis
 .
Melbourne
:
CSIRO
.. p.
64
68
.
Cesaroni
D
Lucarelli
M
Allori
P
Russo
F
Sbordoni
V.
1994
.
Patterns of evolution and multidimensional systematics in graylings (Lepidoptera, Hipparchia)
.
Biol J Linn Soc.
 
52
:
101
−–
111
.
Cook
WM
Lane
KT
Foster
BL
Holt
D.
2002
.
Island theory, matrix effects, and species richness patterns in habitat fragments
.
Ecol Lett.
 
5
:
619
623
.
Cornuet
JM
Luikart
G.
1997
.
Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data.
Genetics
 .
144
:
2001
2014
.
Dapporto
L.
2010
.
Satyrinae butterflies from Sardinia and Corsica show a kaleidoscopic intraspecific biogeography (Lepidoptera, Nymphlidae)
.
Biol J Linn Soc.
 
100
:
195
212
.
Dapporto
L
Habel
JC
Dennis
RLH
Schmitt
T.
2011
.
The biogeography of the western Mediterranean: elucidating contradictory distribution patterns of differentiation in Maniola jurtina (Lepidoptera, Nymphalidae)
.
Biol J Linn Soc.
 
103
:
571
577
.
Dapporto
L
Schmitt
T
Vila
V
Scalercio
S
Biermann
H
Dinca
V
Gayubo
SF
González
JA
Lo Cascio
P
Dennis
RLH
.
2011
.
The phylogeographic island disequilibrium: evidence for macroecology dynamics in Mediterranean butterflies
.
J Biogeogr.
 
38
:
922
932
.
De Noblet-Ducoudré
N
Prentice
CM.
2000
.
Mid-Holocene greening of the Sahara: first results of the GAIM 6000 year BP Experiment with two asynchronously coupled atmosphere/biome models
.
Clim Dynam.
 
16
:
643
659
.
Dennis
RLH
Shreeve
TG.
1989
.
Butterfly wing morphology variation in the British Isles. The influence of climate, behavioural posture and the hostplant-habitat
.
Biol J Linn Soc.
 
38
:
323
348
.
Dennis
RLH
Dapporto
L
Sparks
TH
Williams
SR
Greatorex-Davies
JN
Asher
J
Roy
DB.
2010
.
Turnover and trends in butterfly communities on two British tidal islands: stochastic influences and deterministic factors
.
J Biogeogr.
 
37
:
2291
2304
.
Dincǎ
V
Dapporto
L
Vila
R.
2011
.
A combined genetic-morphometric analysis unravels the complex biogeographical history of Polyommatus icarus and Polyommatus celina common blue butterflies.
Mol Ecol
 .
20
:
3921
3935
.
Dover
J
Sparks
T.
2000
.
A review of the ecology of butterflies in British hedgerows
.
J Env Manag.
 
60
:
51
63
.
Elith
J
Kearney
M
Phillips
S.
2010
.
The art of modelling range-shifting species
.
Meth Ecol Evol.
 
1
:
330
342
.
Excoffier
L
Laval
G
Schneider
S.
2005
.
Arlequin (version 3.0): an integrated software package for population genetics data analysis.
Evol Bioinform Online
 .
1
:
47
50
.
Goudet
J.
1995
.
FSTAT (version 1.2): a computer program to calculate F-statistics
.
Heredity.
 
86
:
485
486
.
Grimaldi
MC
Crouau-Roy
B.
1997
.
Microsatellite allelic homoplasy due to variable flanking sequences.
J Mol Evol
 .
44
:
336
340
.
Habel
JC
Finger
A
Meyer
M
Schmitt
T
Assmann
T.
2008
.
Polymorphic microsatellite loci in the endangered butterfly Lycaena helle (Lepidoptera: Lycaenidae)
.
Eur J Entomol.
 
105
:
361
362
.
Hanski
IA
Gaggiotti
OE
, editors.
2004
Ecology, genetics, and evolution of metapopulations
 .
San Diego (CA)
:
Academic Press
.,
696p
.
Hasumi
H
Emori
S
, editors.
2004
.
K-1 coupled GCM (MIROC) description. K-1 Technical Report No. 1.
Japan:
Center for Climate System Research, University of Tokyo
..
Häuser
CL
Holstein
J
Steiner
A.
2005
.
Digital imaging of butterflies and other lepidoptera – more or less “flat” objects?
In:
Häuser
CL
Steiner
A
Holstein
J
Scoble
MJ
, editors. Digital imaging of biological type specimens – a manual of best practice. Stuttgart (Germany): ENBI. p. 209.
Helsen
P
Vandewoestijne
S
Van Dongen
S
Matthysen
E.
2010
.
Isolation and characterization of ten polymorphic microsatellite markers from the Speckled Wood butterfly (Pararge aegeria, Nymphalidae)
.
Mol Ecol Res.
 
10
:
232
236
.
Hijmans
RJ
Cameron
SE
Parra
JL
Jones
PG
Jarvis
A.
2005
.
Very high resolution interpolated climate surfaces for global land areas
.
Intern J Climatol.
 
25
:
1965
1978
.
Hijmans
RJ
Guarino
L
Jarvis
A
O’Brien
R
Mathur
P
Bussink
C
Cruz
M
Barrantes
I
Rojas
E.
2005
DIVA-GIS version 5.2 manual. Available from: http://www.diva-gis.org/docs/DIVA-GIS5_manual.pdf
Hijmans
RJ
Phillips
S
Leathwick
J
Elith
J.
2011
dismo: species distribution modeling. R package version 0.7-13
. Available from: http://CRAN.R-project.org/package=dismo
Jensen
JL
Bohonak
AJ
Kelley
ST.
2005
.
Isolation by distance, web service.
BMC Genet
 .
6
:
13
.
Joyce
DA
Dennis
RLH
Bryant
SR
Shreeve
TG
Ready
JS
Pullin
AS.
2009
.
Do taxonomic divisions reflect genetic differentiation? A comparison of morphological and genetic data in Coenonympha tullia (Muller), Satyrinae
.
Biol J Linn Soc.
 
97
:
314
327
.
Keller
I
Excoffier
L
Largiadèr
CR.
2005
.
Estimation of effective population size and detection of a recent population decline coinciding with habitat fragmentation in a ground beetle.
J Evol Biol
 .
18
:
90
100
.
Keller
I
Nentwig
W
Largiader
CR.
2004
.
Recent habitat fragmentation due to roads can lead to significant genetic differentiation in an abundant flightless ground beetle.
Mol Ecol
 .
13
:
2983
2994
.
Kemp
DJ
Wiklund
C
Van Dyck
H.
2006
.
Contest behaviour in the speckled wood butterfly (Pararge aegeria): seasonal phenotypic plasticity and the functional significance of flight performance
.
Behav Ecol Sociobiol.
 
59
:
403
411
.
Ladle
RJ
Whittaker
RJ.
2010
Conservation biogeography
 .
Sussex (UK): Wiley-Blackwell
..
p. 301
.
Lesica
P
Allendorf
FW.
1995
.
When peripheral populations are valuable for conservation
.
Conserv Biol.
 
9
:
753
760
.
Luikart
G
Cornuet
JM.
1998
.
Empirical evaluation of a test for identifying recently bottlenecked populations from allele frequency data
.
Conserv Biol.
 
12
:
228
237
.
Meglécz
E
Anderson
SJ
Bourguet
D
Butcher
R
Caldas
A
Cassel-Lundhagen
A
d’Acier
AC
Dawson
DA
Faure
N
Fauvelot
C
et al
2007
.
Microsatellite flanking region similarities among different loci within insect species.
Insect Mol Biol
 .
16
:
175
185
.
Meglécz
E
Solignac
M.
1998
.
Microsatellite Loci for Parnassius mnemosyne (Lepidoptera)
.
Hereditas.
 
128
:
179
180
.
Melbourne
BA
Hastings
A.
2008
.
Extinction risk depends strongly on factors contributing to stochasticity.
Nature
 .
454
:
100
103
.
Merckx
T
Van Dyck
H
Karlsson
B
Leimar
O.
2003
.
The evolution of movements and behaviour at boundaries in different landscapes: a common arena experiment with butterflies.
Proc R Soc Ser B Biol Sci
 .
270
:
1815
1821
.
Nève
G
Meglécz
E.
2000
.
Microsatellite frequencies in different taxa.
Trends Ecol Evol
 .
15
:
376
377
.
Nevo
E
Dessauer
HC
Chuang
KC.
1975
.
Genetic variation as a test of natural selection.
Proc Natl Acad Sci USA
 .
72
:
2145
2149
.
Nylin
S
Wickman
PO
Wiklund
C.
1989
.
Seasonal plasticity in growth and development of the Speckled Wood Butterfly, Pararge aegeria (Satyrinae)
.
Biol J Linn Soc.
 
38
:
155
171
.
Otto-Bliesner
BL
Brady
EC
Clauzet
G
Tomas
R
Levis
S
Kothavala
Z.
2006
.
Last Glacial Maximum and Holocene climate in CCSM3
.
J Climate.
 
19
:
2526
2544
.
Perez
SI
Bernal
V
Gonzalez
PN.
2006
.
Differences between sliding semi-landmark methods in geometric morphometrics, with an application to human craniofacial and dental variation.
J Anat
 .
208
:
769
784
.
Peterson
AT
Nyári
AS.
2008
.
Ecological niche conservatism and Pleistocene refugia in the Thrush-like Mourner, Schiffornis sp., in the neotropics.
Evolution
 .
62
:
173
183
.
Pritchard
JK
Stephens
M
Donnelly
P.
2000
.
Inference of population structure using multilocus genotype data.
Genetics
 .
155
:
945
959
.
Reed
DH
Frankham
R.
2003
.
Correlation between fitness and genetic diversity
.
Conserv Biol
 
17
:
230
237
.
Rohlf
FJ.
2010
TpsDig, digitize landmarks and outlines, version 2.16.
 
Department of Ecology and Evolution, State University of New York at Stony Brook
..
Rohlf
FJ.
2010
TpsUtil, version 1.46.
 
Department of Ecology and Evolution, State University of New York at Stony Brook
..
Rohlf
FJ.
2010
TpsRelw, version 1.49.
 
Department of Ecology and Evolution, State University of New York at Stony Brook
..
Saccheri
I
Kuusaari
M
Kankare
M
Vikman
V
Fortelius
W
Hanski
I.
1998
.
Inbreeding and extinction in a butterfly metapopulation
.
Nature.
 
392
:
491
494
.
Schmitt
T
Besold
J.
2010
.
Up-slope movements and large scale expansions: the taxonomy and biogeography of the Coenonympha arcania-darwiniana-gardetta butterfly species complex
.
Zool J Linn Soc.
 
159
:
890
904
.
Selkoe
KA
Toonen
RJ
.
2006
Microsatellites for ecologists: a practical guide to using and evaluating microsatellite markers.
Ecol Lett
 .
9
:
615
629
.
Shaibi
T
Moritz
RFA.
2010
.
10,000 years in isolation? Honeybees (Apis mellifera) in Saharan oases
.
Conserv Genet.
 
11
:
2085
2089
.
Shapiro
AM
Porter
AH.
1989
.
The lock-and-key hypothesis: evolutionary and biosystematic interpretation of insect genitalia
.
Ann Rev Entomol.
 
34
:
231
245
.
Shreeve
TG.
1984
.
Habitat selection, mate location, and microclimatic constraints on the activity of the speckled wood butterfly Pararge aegeria
.
Oikos.
 
42
:
371
377
.
Sinama
M
Dubut
V
Costedoat
C
Gilles
A
Junker
M
Malausa
T
Martin
J-F
Nève
G
Pech
N
Schmitt
T.
et al
2011
.
. Challenge for microsatellite development in Lepidoptera: Euphydryas aurinia Rottemburg, 1775 (Nymphalidae) as a case study.
Eur J Entomol.
 
108
:
261
266
.
Slatkin
M.
1995
.
A measure of population subdivision based on microsatellite allele frequencies.
Genetics
 .
139
:
457
462
.
Stevens
DJ
.
2004
Pupal development temperature alters adult phenotype in the speckled wood butterfly
,
Pararge aegeria.
  J Thermal Biol.
29
:
205
210
.
Tallmon
DA
Koyuk
A
Luikart
GH
Beaumont
MA.
2008
.
ONeSAMP: a program to estimate effective population size using approximate Bayesian computation
.
Mol Ecol Res.
 
8
:
299
301
.
Thomson
G.
2011
The Meadow Brown butterflies. A study in genetics, morphology and evolution
 .
Waterbeck (Scotland)
:
Privately published by George Thomson
..
Vandewoestijne
S
Schtickzelle
N
Baguette
M.
2008
.
Positive correlation between genetic diversity and fitness in a large, well-connected metapopulation.
BMC Biol
 .
6
:
46
.
Vandewoestijne
S
Van Dyck
H
.
2010
.
Population genetic differences along a latitudinal cline between original and recently colonized habitat in a butterfly.
PLoS ONE
 .
5
:
e13810
.
Vandewoestijne
S
Van Dyck
H.
2011
.
Flight morphology along a latitudinal gradient in a butterfly: do geographic clines differ between agricultural and woodland landscapes?
Ecography.
 
34
:
876
886
.
Van Dyck
H
Matthysen
E
Dhondt
AA
.
1997
The effect of wing colour on male behavioural strategies in the speckled wood butterfly
.
Animal Behav
 .
53
:
39
51
.
Van Dyck
H
Wiklund
C.
2002
.
Seasonal butterfly design: morphological plasticity among three developmental pathways relative to sex, flight and thermoregulation
.
J Evol Biol
 
15
:
216
225
.
Van Oosterhout
C
Hutchinson
WF
Wills
DPM
Shipley
P.
2004
.
MICRO-CHECKER (version 2.2.3): software for identifying and correcting genotyping errors in microsatellite data
.
Mol Ecol Notes.
 
4
:
535
538
.
Weingartner
E
Wahlberg
N
Nylin
S.
2006
.
Speciation in Pararge (Satyrinae: Nymphalidae) butterflies – North Africa is the source of ancestral populations of all Pararge species
.
Sys Entomol.
 
31
:
621
632
.
Wiklund
C.
2003
.
Sexual selection and the evolution of butterfly mating systems
. In:
Boggs
CL
Watt
BW
Ehrlich
PR
, editors.
Butterflies - ecology and evolution taking flight
 . Chicago (IL): University of Chicago Press. p.
67
78
.
Wiklund
C
Persson
A.
1983
.
Fecundity, and the relation of egg weight variation of offspring fitness in the speckled wood butterfly Pararge aegeria, or why don′t butterfly females lay more eggs?
Oikos.
 
40
:
53
63
.

Author notes

Corresponding Editor: Adalgisa Caccone