Introduction

Despite the fact that clonal reproduction is a priori expected to be reproductively advantageous (Maynard-Smith, 1978; Bell, 1982), obligate asexual organisms are evolutionarily short lived and have been traditionally regarded as evolutionary dead ends (Bell, 1982; Judson and Normark, 1996). It has been suggested that asexual organisms accumulate deleterious mutations more quickly than sexuals (Muller, 1964; Kondrashov, 1988) and that co-evolutionary arm races with parasites, predators and competitors will favor sexual reproduction (Lively, 2010). One way by which asexuals may temporarily escape these forces is through the observed differences in the distribution of closely related sexual and asexual lineages.

Asexuals tend to occupy more isolated habitats, have broader distributions and are often found in marginal habitats, relative to their sexual counterparts; a pattern that may provide insight into potential mechanisms by which asexual lineages may persist, at least in the short term. Thus, a number of hypotheses have been proposed to explain the observed patterns of distribution of sexual and asexual lineages: (1) the superior colonizing ability of parthenogens could explain their larger geographical ranges, given that every single asexual individual is able to found a new population, because finding a mate is not necessary to ensure reproduction (Cuéllar, 1994); (2) the interactions with parasites, predators and competitors will influence the distribution of asexual forms, and thus parthenogens will be more common in sparsely inhabited habitats, where these interactions are rare, and where environmental factors are dominating (Glesener and Tilman, 1978; Hamilton et al., 1990); (3) the lower fitness of asexual/bisexual hybrids might keep populations separated (Lynch, 1984); (4) heterozygosity assurance provided by clonal reproduction will favor asexual lineages in subdivided metapopulations (Vrijenhoek, 1985; Haag and Ebert, 2004). Most of these models are based on the idea that parthenogens may persist in the long term if they escape from negative interactions with either sexual lineages or biological enemies (predators and parasites).

Ischnura hastata (Zygoptera, Coenagrionidae) is the only known example of parthenogenetic reproduction within the insect order Odonata, and represents a typical example of geographic parthenogenesis: sexual populations of this species are widely distributed in the American continent, whereas only female, parthenogenetic populations have been described only in the Azores archipelago (Cordero Rivera et al., 2005). This system allows us the opportunity to study the mechanisms that maintain parthenogenesis and how they have arisen in a previously unidentified group. It also allows us to investigate the influence of geographic isolation on the evolution of parthenogenetic populations, and to contribute further to our understanding of geographic parthenogenesis.

Previous molecular genetics work strongly suggests that parthenogenetic I. hastata females are diploid and that the number of chromosomes of the unfertilized eggs is restored by an apomictic mechanism of parthenogenesis; although these results are not entirely conclusive without further cytological and karyotype analyses (see Lorenzo-Carballa and Cordero-Rivera, 2009). No evidence for bacterial involvement in the parthenogenesis has been found, which suggests that tychoparthenogenesis (that is, a small capacity for parthenogenesis in otherwise sexually reproducing females) could occur in the sexual populations of this species (Lorenzo-Carballa and Cordero-Rivera, 2009). I. hastata has been inferred to be a migratory species (Corbet, 1999) based on its widespread distribution in the American continent. It is found in temporal or permanent water bodies, which suggests that the sexual demes may have an island-like structure subject to frequent extinction and/or recolonization, as described for other odonate species (Gibbons et al., 2002; Purse et al., 2003; De Block et al., 2005; Rouquette and Thompson, 2007). This potential for long distance dispersal should favor population mixing, leading to low levels of genetic subdivision.

In order to gain insight in the origin and patterns of distribution of parthenogenetic I. hastata lineages, we have studied the population structure of this species over a wide geographic area. Here, we use microsatellite loci, mitochondrial and nuclear DNA sequence data to examine the population genetic structure of sexual and parthenogenetic populations of I. hastata. Our results suggest that the North-American demes of this species are strongly interconnected by dispersal so they can be considered a single population, and that the asexual population found in the Azores archipelago is the result of a demographic expansion after a recent single colonization event.

Materials and methods

Sample collection and DNA extraction

Adult individuals of I. hastata were collected from 14 locations: six parthenogenetic populations (Pico, Terceira, Faial, Flores, São Miguel and Corvo), covering most of the species’ range in the Azores islands and eight sexual populations spread across North America (Arkansas; Kentucky; Texas; Florida, Lake County; Florida, Vero Beach; México, Las Fuentes and México, Vicam) and the Caribbean (Cuba) (Figure 1; Table 1). Genomic DNA from all of these samples was extracted using a CTAB protocol (modified from Doyle and Doyle, 1987) or the NucleoSpin Tissue kit (Macherey-Nagel, Düren, Germany), following the manufacturer's instructions.

Figure 1
figure 1

Sampling locations of specimens of Ischnura hastata in North America (a) and the Azores islands (b). Sites are labelled according to abbreviations in Table 1.

Full size image
Table 1 Ischnura hastata: geographic distribution of the analyzed specimens
Full size table

Sequencing and haplotype determination

Nucleotide variation was assessed in one mitochondrial (Cytochrome c Oxidase I, COI) and one nuclear gene (Elongation Factor 1-α, EF1-α). Genomic DNA was PCR amplified using locus-specific primers and conditions (see Supplementary information). Sequencing reactions were carried out bidirectionally either on a MegaBACE 500 sequencer, using the DYEnamic ET Dye Terminator Cycle Sequencing Kit (Amersham Biosciences, Buckinghamshire, UK) or on a 3730xl DNA Analyzer, using the BigDye Terminator Kit (Applied Biosystems, Foster City, CA, USA). Sequencing results were consistent across platforms, and only high quality traces were included in our analyses. For EF1-α data, allelic phases of heterozygous individuals were inferred using PHASE (Stephens et al., 2001). Excluding autapomorphies, all identified haplotypes had posterior probabilities >0.8.

Microsatellite genotyping

Nine polymorphic microsatellite loci (Lorenzo Carballa et al., 2007) were used to screen 144 individuals from sexual and parthenogenetic populations (Table 3). Forward primers were labelled with a fluorescent dye (HEX or 6-FAM, Eurofings MWG Operon, Ebesberg, Germany) for microsatellite typing. PCR conditions were the same as described in Lorenzo Carballa et al. (2007). Automated genotyping was performed on a MegaBACE 500 automated DNA sequencer using the ET-550 Size Standard (Amersham Biosciences). Fragment size was recorded and scored using Genetic Profiler v. 1.2 (Amersham Biosciences).

Genetic diversity analyses

Forward and reverse sequence strands were assembled and edited using SeqManII v. 5.03 (DNAstar, Inc., Madison, WI, USA) and consensus sequences were aligned using Clustal W (Thompson et al., 1994), as implemented in MEGA v. 4.0 (Tamura et al., 2007). The population genetic parameters number of haplotypes (h), number of segregating sites (S), haplotype diversity (H), nucleotide silent diversity (πS) and number of synonymous and non-synonymous substitutions for each locus and population were calculated using DNAsp v. 4.20.2 (Rozas et al., 2003). COI and EF1-α haplotypes are deposited in GenBank with accession numbers JN578143–JN578216.

For the microsatellite data set, allelic richness, observed and expected heterozygosities, and the number of private alleles for each locus and population were calculated using the software GDA (Lewis and Zaykin, 2001). Deviations from Hardy–Weinberg equilibrium for each locus in each population were calculated using Arlequin v. 3.01 (Excoffier et al., 2005). Linkage disequilibrium analyses in the sexual populations were performed using correlation-based tests for linkage disequilibrium (Zaykin et al., 2008).

Phylogenetic and phylogeographic analyses

Genealogical relationships between I. hastata haplotypes were analyzed using statistical parsimony networks (95% probability criterion) as implemented in TCS 1.21 (Clement et al., 2000). Additionally, we used multiple tree-building methods (NJ, Maximum Likelihood and Bayesian Inference) to reconstruct phylogenetic relationships among these haplotypes. The results obtained with these methods were almost identical to the reported networks (data not shown).

For the COI data set, we used the refined nesting algorithm of the nested clade phylogeographic analysis (Templeton et al., 1995) as implemented in Geodis v. 2.5 (Posada et al., 2000) to identify haplotype clades showing significant geographic associations. However, because of the severe limitations of the nested clade phylogeographic analysis to infer the multiple historical processes that may characterize a species’ history (see Knowles and Maddison, 2002; Petit and Grivet, 2002; Knowles, 2004; Panchal and Beaumont, 2007; Petit, 2008a, 2008b; Beaumont et al., 2010) we did not use this analysis to make any inferences about the demographical or historical processes that could have shaped the genetic variation observed in I. hastata.

Population structure

First, we assessed population structure assuming a priori populations (that is, sampling locality=population), using hierarchical analysis of molecular variance (Excoffier et al., 1992) including (i) populations, (ii) populations within reproductive mode and (iii) reproductive mode. Analysis of molecular variance was also implemented to analyze the genetic structure of sexual populations. Analyses were conducted for all loci (COI, EF1-α and microsatellites) in Arlequin v. 3.01 (Excoffier et al., 2005). For the microsatellites we selected a locus-by-locus strategy, which is the most desirable option when missing data exist (Excoffier et al., 2005).

Pairwise Fst values were estimated with Arlequin v. 3.01 (Excoffier et al., 2005), using pairwise distances among EF1-α and COI sequences. The significance of Fst for population comparisons was assessed by 10 000 permutations.

To test for a possible scenario of isolation by distance in the sexual populations, we tested the correlation of genetic differentiation over geographical distances for all pairs of populations. Geographical distances among locations were estimated using the latitude/longitude distance calculator (http://jan.ucc.nau.edu/~cvm/latlongdist.html), and correlations between geographic and genetic distances for both microsatellites and COI were tested using a Mantel permutation procedure (Sokal and Rohlf, 1995) as implemented in IBD v. 3.15 (http://ibdws.sdsu.edu/~ibdws/; Jensen et al., 2005). The significance of the Mantel test was assessed by performing 10 000 random iterations.

To get insight into the substructuring among the sexual populations, the results with predefined populations were then compared with the Bayesian model-based clustering methods implemented in Structure v. 2.1 (Pritchard et al., 2000). We used the admixture model, and the number of clusters, K, was estimated by comparing the log-likelihood ratios in two independent runs for values of K between 1 (panmixia) and 4 (the total number of independent sites sampled for microsatellites). Each run consisted in 106 iterations, with a burn-in period of 105. All summary statistics stabilized before the end of the burn-in period and independent runs always attained the same results.

Patterns of demographic history in sexual populations

We used COI data to infer patterns of demographic history in sexual I. hastata. First, two tests based on the distribution of segregating sites (Tajima's D; Tajima, 1989) and on the haplotype distribution (Fu's FS; Fu, 1997) were applied to the data. Negative values of these tests are expected in populations that have undergone a recent expansion, because rare alleles are more numerous than expected; whereas positive values occur if rare alleles are eliminated from population following genetic bottlenecks (Tajima, 1989). Significance of D and FS statistics was tested by generating 10 000 coalescent simulations in Arlequin v. 3.01 (Excoffier et al., 2005). To further investigate the possibility of demographic expansions, we plotted the distribution of pairwise sequence differences (mismatch distribution) for all COI haplotypes using Arlequin v. 3.01 (Excoffier et al., 2005). Unimodal mismatch distributions are characteristic of haplotypes drawn from a population that has undergone a recent expansion and is the null hypothesis by which the mismatch distribution is assessed. Significant divergences from this model (multimodality) reflect the highly stochastic shape of gene trees under demographic equilibrium (Rogers and Harpending, 1992). All these analyses were carried out in each of the sexual localities sampled, as well as in the entire sexual data set.

As an independent test for demographic events, we used microsatellite data to test for the possibility of recent population size reduction and expansion in the sexual cluster. For bottleneck detection, calculations were performed with the program BOTTLENECK (Cornuet and Luikart, 1997), using the Two-Phase Model of mutation as suggested by Luikart et al. (1998). Wilcoxon's signed rank tests were used to determine whether populations exhibit significant number of loci with gene diversity excess. Signatures of population expansion were examined using the β imbalance index (Kimmel et al., 1998). This statistic is analogous to a mismatch distribution, and relies in two different estimates of the genetic diversity parameter θ (4Neμ). One estimate is based on allele size variance (θV), and the other is estimated from the variance in repeat numbers (θP). The imbalance index should be <1 for expanding populations whose pre-expansion history is stable, but is characteristically >1 in populations with a reduction in size that precedes a detectable growth phase (Kimmel et al., 1998, King et al., 2000). Here, we calculated β1, the difference in the natural logarithm of these estimators, averaged over all loci. This logarithm-based metric was chosen because of its higher sensitivity to detect historic signals of population expansion (King et al., 2000). Confidence intervals of β1 were estimated by bootstrapping over loci using R 2.9.0.

Results

Genetic diversity

Analyses of sequence polymorphisms in sexual and parthenogenetic I. hastata are summarized in Table 2. For EF1-α, 228 sequences were obtained. In the 46 inferred haplotypes 44 out of 459 sites (9.59 %) and two indels were found to be polymorphic. Only one haplotype was found in the parthenogenetic populations, which was also present in the sexual populations. The other 45 haplotypes were found only in the sexual individuals: one haplotype was present in four populations, seven haplotypes were present in two populations and three were found in different individuals from the same population. Each of the remaining (38) was found in single individuals. Nei's nucleotide diversity (π) varied from 0.006 to 0.022 in the sexual populations, with an average of 0.014. The number of haplotypes per location varied from 3 (K) to 10 (VB), and the number of private haplotypes per site fluctuated between 2 (VI) and 7 (VB), which yields an overall haplotype diversity in the sexual populations of 0.9 (ranging from 0.78 to 1).

Table 2 Polymorphism statistics for each gene
Full size table

One-hundred and fifty-eight sexual and parthenogenetic I. hastata individuals were sequenced for COI. In all, 34 out of 550 sites were polymorphic (6.18%) and 28 haplotypes were identified, 2 of which were found only in the parthenogenetic populations. The most common haplotype in the parthenogens was found in 54 individuals, whereas the second, which differs only in 1 nucleotide, was only found in 2 individuals from São Miguel Island. In the sexual populations, 26 haplotypes were found: 1 was shared by all sexual populations and represented 61.85% of the total number of sexual individuals, 1 haplotype was found in 4 populations, 3 haplotypes were found in 2 populations and 2 were found in different individuals of the same population. The rest (19) were each found in single individuals. Nucleotide diversity varied from 0 to 0.01, with an average of 0.003. The number of haplotypes per location fluctuated between 11 (FU) and 1 (CU), and the number of private haplotypes between 7 (FU) and 0 (CU), which yields an average haplotype diversity of 0.51 (ranging from 0 to 0.95).

Polymorphism was analyzed across 9 microsatellite loci in 144 individuals from 10 populations (Table 3). Allelic diversity ranged from 7 to 13 alleles per locus, with 87 alleles identified across all 9 loci. All the alleles present in the parthenogenetic females were also found in the sexual populations, with the exception of one allele at microsatellite locus Ihas16, exclusive of the parthenogenetic populations, that was found in one individual from Pico Island.

Table 3 Microsatellite polymorphisms
Full size table

Given that apomictic parthenogenesis produces offspring that is heterozygous at all loci in which their mothers were heterozygous (Lorenzo-Carballa and Cordero-Rivera, 2009), deviations from Hardy–Weinberg equilibrium might be expected at segregating loci. Accordingly, for variable microsatellite loci (Ihas05, Ihas08, Ihas09, Ihas11, Ihas13 and Ihas16) all parthenogenetic females were heterozygous (Table 3). Variation in parthenogens was only observed at locus Ihas16, with four different alleles that defined the three highly related clones found in the asexual population (Table 4). In contrast, 40 unique multilocus genotypes were found in the sexuals. In all the sexual populations, we observed a deficit of heterozygotes across multiple loci, although this excess of homozygotes was only significant (P<0.05) in the populations from Mexico (FU and VI) and Florida (LA and VB) (Table 3). Cases of linkage disequilibrium were found between loci Ihas01-Ihas05, Ihas05-16 (LA); Ihas01-Ihas09, Ihas05-Ihas13 (VB); Ihas01-Ihas10, Ihas08-Ihas09 (VI); and Ihas01-Ihas09, Ihas01-Ihas11 (FU).

Table 4 Microsatellite multilocus genotypes found in the parthenogenetic Ischnura hastata populations
Full size table

Phylogenetic and phylogeographic analyses

The COI and EF1-α haplotype networks are shown in Figure 2. The main characteristic of these networks is the lack of genetic structure in the sexual populations. For the COI, parthenogenetic I. hastata cluster together, whereas for the EF1-α, there is no separation among sexual and asexual populations, suggesting a relative recent divergence between both lineages. Both loci networks showed a star-like shape, with the most frequent haplotypes in the center, surrounded by several low frequency haplotypes; a pattern typically observed in expanding populations.

Figure 2
figure 2

Mutational haplotype networks based on statistical parsimony, representing genealogical relationships between 28 COI haplotypes (a) and 46 EF1-α haplotypes (b) in sexual and parthenogenetic populations of Ischnura hastata. Black dots represent missing mutational steps or non-sampled haplotypes. Haplotypes connected by a single line differ only in one mutational step. The size of the circles correlates with haplotype frequency. The dashed lines in the COI network indicate the two main clades with significant geographic association (see main text).

Full size image

For the COI data, the nested clade phylogeographic analysis nesting algorithm resolved two two-step clades that showed significant geographic association: clade I, which contained the most common haplotype in sexual populations and its surrounding haplotypes and clade II, which included the haplotypes found in the parthenogenetic populations and their more closely related sexual haplotypes. Parthenogens differ from their closest sexual relative by a single mutational step (Figure 2).

Population structure

I. hastata exhibited significant partitioning of genetic variation. The analysis of molecular variance including sexual and parthenogenetic populations showed that for the nuclear markers almost all genetic variation (∼60–85%) resides within the sampled populations. In contrast, the main source of genetic variation for mtDNA (∼72%) is attributed to haplotype differences between reproductive groups (sexual vs asexual populations; Table 5). Within these groups, both the nuclear and the mitochondrial genome showed very little population structure. A complementary analysis of molecular variance of the sexual populations showed a very similar pattern. At a continental scale, I. hastata can be characterized as a large subdivided population where the vast majority of the observed variation corresponds to variation within the independently sampled populations (Table 5). Fst estimates in North-American (that is, sexual) populations were higher for EF1-α than for any other locus reflecting the existence of many population-exclusive alleles in this marker (Table 6). We found a strong correlation between pairwise Fst values and geographic distances for both the COI (r=0.5508, Mantel test P=0.0062) and the microsatellite data sets (r=0.8416, P=0.1626). However, as we were only able to genotype four of the sexual populations (see Table 3) the latter correlation is not significant.

Table 5 AMOVA and hierarchical analyses for Ischnura hastata populations
Full size table
Table 6 Pairwise Fst values between sexual Ischnura hastata populations based on COI (upper diagonal) and EF1-α (lower diagonal)
Full size table

As for the phylogenetic analyses, the microsatellite-based Bayesian population assignment tests failed to detect any genetic structure among sexual populations: an increase in the number of clusters resulted in nearly the same probability for all the sexual individuals to be assigned to each of the sampled localities (Figure 3).

Figure 3
figure 3

Genetic clustering of sexual Ischnura hastata, based on nine microsatellite loci. Estimated population structure (from K=1 to K=4). Sampled localities (labelled below the figure) are separated by black lines. Each individual is represented by a narrow vertical column with the proportion of the colors indicating the genome proportion derived from each population.

Full size image

Patterns of demographic history in sexual populations

Because demographic events in continental (that is, sexual) populations might help to explain the colonization of the Azores archipelago, we conducted a series of analyses aimed to characterize the recent demographic history of I. hastata in North America. Overall sequencing data (Fu's FS, Tajima's D and mismatch statistics) rejected constant population and support the hypothesis of population expansion. When the data from each locality were analyzed separately, neutrality tests were significantly negative at Vero Beach; and significantly negative FS values were also found in Kentucky and Las Fuentes. With the exception of Vicam, a trend toward negative values of both statistics was found in all populations. According to the mismatch distributions, the hypothesis of a sudden expansion could not be rejected in most populations, thus providing further evidence of population expansion (Table 7). Similarly, the frequency distribution of pairwise differences in the sexual cluster was not significantly different from that expected under a model of demographic expansion, and neutrality tests were also significantly negative (Figure 4).

Table 7 Demographic estimates from COI for each of the sexual localities sampled
Full size table
Figure 4
figure 4

Mismatch distributions for the sexual group of Ischnura hastata. Raggedness index and probability levels are given under the null distribution (model frequency) consistent with a sudden population expansion.

Full size image

The microsatellite data set further confirmed this pattern of population expansion. The observed imbalance index for the sexual cluster (β1 mean: −0.492; range: −0.889 to −0.137) suggests that the sexual population has recently expanded without significant previous bottlenecks. Consistently, our bottleneck analyses did not detect a non-significant excess of heterozygosity across all loci (Wilcoxon's test, P=0.684), suggesting that the sexual population has not suffered significant bottlenecks.

Discussion

Genetic diversity, structure and demography of sexual populations

Our results suggest that in I. hastata, local populations may be so strongly interconnected that they can be better considered as a single ‘patchy population’, where gene flow (that is, dispersal) exerts a strong homogenizing force over wide geographic ranges. Bayesian model-based clustering, phylogenetic and genetic differentiation analyses suggest that at a broad continental scale there is no significant population structure, and that very little genetic differentiation exists between demes. In fact, most of the total genetic variation found across North America (average 97% for microsatellite loci) is found within each of the independently sampled populations. However, some of these local populations are significantly differentiated, suggesting varying amounts of dispersal among localities. Only moderate population structure is detected over distances of 1400 km, but the vast extent of the North-American subcontinent allows significant population differentiation. Thus, demes may currently be panmictic over a regional scale with low or sporadic gene flow over greater distances. This restricted pattern of gene flow over very long distances, combined with high levels of local (or even regional) recruitment, can explain the significant association between genetic differentiation and geographic distance detected in I. hastata populations (see Palumbi, 2003). Yet, other non-equilibrium alternatives cannot be entirely excluded. For example, a very recent disruption of gene flow (for example, Pogson et al., 2001), or an ongoing turnover of local populations where colonizers tend to be from different demes (Wade and McCauley, 1988) could potentially explain the observed geographic patterns of genetic variation. The analyses of the site-frequency distribution in two different markers (that is, mitochondrial COI and nuclear EF1-α) revealed a relatively small, but consistent, excess of rare alleles within local populations. This excess could be the result of sampling admixture (Ptak and Przeworski, 2002; Hammer et al., 2003), as the genetic pool of adults flying in each local population includes both local recruitment and natal dispersal (Corbet, 1999). Accordingly, analyses of genetic variation at nine unlinked microsatellite loci revealed a significant genome-wide reduction of heterozygosity, suggesting that sexual populations are indeed comprising genetically dissimilar individuals from different neighboring locations (Freeland et al., 2000).

This pattern of population structure that we observe in I. hastata is the result not only of the present-day gene flow, but also reflects historical processes. Both slow (mitochondrial COI and nuclear EF1-α) and fast (microsatellites) evolving polymorphic markers have allowed us to infer the demographic history of this species in North America. First, the star-like shape of the nuclear and mitochondrial haplotype networks is indicative of a population expansion occurred during the recent history of the species. Neutrality tests of mitochondrial DNA provided also strong indication for the deviation from mutation-drift equilibrium, and mismatch analysis was unimodal for the sexual cluster, thus providing further indication of population expansion.

The erosion of mitochondrial polymorphism might be also due to natural selection either purifying or adaptive (Bazin et al., 2006), since the lack of recombination in mitochondrial DNA makes it act as a single locus, and thus selection acting anywhere in the mtDNA would exactly mimic a demographic bottleneck/expansion. To distinguish between these two evolutionary processes, we have used the imbalance index and bottleneck tests to infer aspects of population history and demography of sexual I. hastata using the microsatellite data. We found concordance between the results of these analyses and the ones obtained from COI data, suggesting that the sexual populations of this species have suffered a recent expansion, and that the demography of pre-expansion was relatively stable.

Genetic diversity of asexual populations: age and origin of parthenogenetic lineage

Our results confirm that sexual populations of I. hastata are connected by dispersal over a wide geographic area. The detection of a pattern of population expansion in multiple, independent, loci allows us to distinguish it from selection and to attribute it to a past demographic change, which is consistent with the low levels of genetic differentiation currently found across the species’ range in North America. This range expansion is likely to be related to the biology and habitat preferences of I. hastata, a pioneer species usually found in temporary ponds where local extinctions are frequent, and adult individuals are thus forced to emigrate to new suitable habitats. Weather fronts have a major role in the dispersal of odonates, and long distance dispersal being passively transported as part of the aerial plankton is frequent in the case of weak flyers such as I. hastata. In fact, this species has been reported over water far from mainland in the Gulf of Mexico (Corbet, 1999). The possibility then exists that this species has been transported to the Azores archipelago by the surface winds associated with the Gulf Stream. Thus, the origin of the populations of this species in the Azores islands is likely to be the result of a single long distance dispersal event, rather than the result of a human introduction associated with shipping trade between the North-American continent and the Azores archipelago (Cordero Rivera et al., 2005).

As expected if the parthenogenetic populations from the Azores have their origin in the sexual populations that inhabit the American continent, microsatellite alleles and the only EF1-α allele found in the islands are a subset of those found in the sexual populations. In contrast, mtDNA revealed parthenogenetic populations as a differentiated cluster, and none of the two almost identical COI haplotypes that are present in the archipelago are found in any of the sexual populations sampled. This observed discordance between nuclear and mitochondrial markers can be potentially explained by a recent gene exchange between both reproductive modes, although this hypothesis is extremely is unlikely, given the great geographic distance that exists between both lineages. Instead, the observed discordance is better explained as the result of a recent fixation of a single diploid thelytokous lineage in the Azores archipelago (see Lorenzo-Carballa and Cordero-Rivera, 2009).

Because apomixis (that is, thelytokous parthenogenesis in which the eggs do not undergo meiosis) generates instant complete linkage disequilibrium between the mitochondrial and nuclear genome, shortly after the fixation there should be no genetic variation among clones. However, there is at least one very recently derived (that is, S=1) COI haplotype, and one derived microsatellite allele segregating in the population at frequencies of 3.5 and 0.6%, respectively, indicating that the ancestral parthenogenetic lineage has been in the archipelago long enough to start accumulating mutations.

The use of mitochondrial DNA to study sexual–asexual lineage divergence can lead to invalid age estimates, due to the extinction of sexual lineages, an inadequate sampling, or if sexual populations have diverged considerably (see Johnson, 2006; and references therein). In our case, we found a single mutational step between sexual and parthenogenetic lineages suggesting a very recent origin for the parthenogenetic lineages. Population frequencies of the Ihas16 alleles can be used to obtain maximum likelihood estimates of the age of the specific parthenogenetic clones (Slatkin and Rannala, 1997; data not shown) However, these estimates showed wide confidence intervals, and could not reject a very recent origin (<20 generations ago), again indicating the young age of the parthenogenetic lineages.

Despite the great geographic area comprised by the sexual localities sampled in this study, they do not represent the actual distribution range of sexual I. hastata, which has been recorded from Canada to South America. A more extensive sampling effort in the American continent would probably reveal higher genetic diversity in the sexual lineage and could help to better understand metapopulation dynamics in this species. Also, the fact that we have found three highly related microsatellite multilocus genotypes in the asexual populations suggests that genetic variability in the parthenogens could be higher than what has been detected in the current study. Thus, including more samples from the Azores, and increasing the resolution power of the analyses by adding more loci (for example, amplified fragment lenght polymorphisms, AFLPs) will probably reveal more diversity in the asexual populations and could also help to better date the time of divergence among both lineages.

One of the most intriguing questions that arise from this work is whether the I. hastata population found at the Azores arose from the immigration of sexual damselflies, which subsequently evolved into parthenogens, or whether parthenogens originated in America and then dispersed to the archipelago. Thus, asexuality could have arisen once in the Azores, and asexuals subsequently outcompeted and drove to extinction their sexual ancestors, given that the lower intrinsic rate of demographic growth of sexual reproduction would theoretically cause sexual individuals to be replaced by parthenogens, all else being equal (Maynard-Smith, 1978). Alternatively, the capacity for parthenogenetic reproduction could be already present in the ancestral sexual population, and therefore in the individuals that colonized the Azores. In this case, parthenogenetic reproduction could have been advantageous for the colonization of the islands. However, with the data available, it is not possible to answer this question. Future research should focus on the study of the potential for parthenogenetic reproduction in sexual females, to better understand the ecological conditions that could have favored the loss of sex in this species.