Introduction
Cicada barbara Stål is usually found associated with olive, pine and, occasionally, eucalyptus trees and occurs in scattered warm habitats in the Iberian Peninsula, whereas it is widely distributed in northwest Africa (Boulard, Reference Boulard1982; Quartau, Reference Quartau1988). Isozyme studies have revealed that C. barbara is a genetically homogeneous species in Portugal, especially when compared to C. orni, a species with which it occurs in sympatry in some areas of the Iberian Peninsula (Quartau et al., Reference Quartau, Ribeiro, Simões and Crespo2000, Reference Quartau, Ribeiro, Simões and Coelho2001). One potential explanation for this homogeneity is that C. barbara is a relatively recent immigrant to the Iberian Peninsula from North Africa. A North African origin for this species has previously been proposed by Boulard (Reference Boulard1982), who classified specimens from Portugal as a subspecies C. barbara lusitanica, distinct from the North African subspecies C. barbara barbara, based on differences at the male genitalia and female ovipositor.
Traditional recognition of subspecies based on discontinuities in the geographical distribution of phenotypic traits (Mayr & Ashlock (1991) inPhillimore & Owens, Reference Phillimore and Owens2006) has often been challenged by molecular phylogenetic data, especially when the phylogenetic monophyly of a subspecies is not recovered (e.g. Zink, Reference Zink2004; Phillimore & Owens, Reference Phillimore and Owens2006). However, historical introgression or incomplete lineage sorting might reduce the probability of finding reciprocal monophyly in recently divergent populations, inducing an underestimation of the differentiation (Phillimore & Owens, Reference Phillimore and Owens2006). Nevertheless, patterns of genetic variation found within species at present can reflect differentiation due to geographic isolation, as well as the retreats and expansions of populations due to climatic changes (Hewitt, Reference Hewitt1996, Reference Hewitt1999, Reference Hewitt2000, Reference Hewitt2001, Reference Hewitt2004a; Taberlet et al., Reference Taberlet, Fumagalli, Wust-Saucy and Cosson1998; Besnard et al., Reference Besnard, Khadari, Baradat and Bervillé2002a; Gómez et al., Reference Gómez, Vendramin, González-Martínez and Alía2005). During the ice ages, many species, particularly thermophilic species (such as cicadas), could only survive in favourable regions or refugia, so genetic differentiation found in the present day reflects their survival in different refugia combined with genetic drift and founder effects during re-colonization. The Iberian Peninsula and Maghreb have been frequently recognised as a key southern refugium for many European species during the ice ages (the Atlantic-Mediterranean refugial area; sensuSchmitt, Reference Schmitt2007), but the survival of exclusively Mediterranean taxa would also have been constrained by these climatic oscillations. According to Jansson & Dynesius (Reference Jansson and Dynesius2002), in areas with little climatic change over long periods of time (low ORD, Orbitally Forced Range Dynamics, over an entire period of Milankovitch oscillations, about 100 kyr) such as southern Europe and the tropics, cladogenesis and anagenesis may be induced by a wide range of processes. As a result, high intraspecific genetic variation and finely geographically subdivided species are expected to occur, while, in regions with high ORD, climatic shifts will alter the dynamics of genetic diversity directly by subsampling gene pools and indirectly by selecting for vagility and generalism, consequently decreasing genetic variability.
Information on the influence of the ice ages on Mediterranean taxa has become increasingly available (e.g. Hewitt, Reference Hewitt2004b; Schmitt, Reference Schmitt2007), and different geological and geophysical conditions have been identified as possible promoters of the present-day distribution and variability of Mediterranean species. The strait of Gibraltar has separated the Iberian Peninsula from North Africa since the end of the Messinian Salinity Crisis about 5.33 Ma (Hsü et al., Reference Hsü, Montadert, Bernoulli, Bianca, Erickson, Garrison, Kidd, Mèliéres, Müller and Wright1977), dividing the Atlantic-Mediterranean refugial area and representing a barrier for several taxa, resulting in vicariant differentiation (Hewitt, Reference Hewitt2004b; Albert et al., Reference Albert, Zardoya and García-París2007; Pinho et al., Reference Pinho, Harris and Ferrand2007; Schmitt, Reference Schmitt2007). On the other hand, with the increase in the volume of the ice sheets, sea levels dropped and, during the Late Glacial Maximum, were 120–130 m (Yokoyama, et al., Reference Yokoyama, Lambeck, De Deckker, Johnston and Fifield2000) lower than at present. As a consequence, the European and the African continental shelves emerged, the distance between these continents was reduced, and several islands emerged (Collina-Girard, Reference Collina-Girard2001), providing stepping stones for different organisms between the two continents (e.g. snakes (Colubridae) and plants (Asteraceae) (Carranza et al., Reference Carranza, Arnold and Pleguezuelos2006; Ortiz, et al., Reference Ortiz, Tremetsberger, Talavera, Stuessy and García-Castaño2007). In the Iberian Peninsula, the expansion and contraction of ice sheets have been implicated as generators of intraspecific diversity and substructure within Iberian taxa (e.g. Paulo et al., Reference Paulo, Dias, Bruford, Jordan and Nichols2001; Hewitt, Reference Hewitt2004b).
Here, we attempt to improve understanding of the demographic history and population structure of C. barbara in the Iberian Peninsula and North Africa, based on genetic analysis of a section of the cyt b gene. The mtDNA gene cytochrome (cyt) b has not been as widely used in insects as other mitochondrial genes, such as cytochrome oxidase I (Caterino et al., Reference Caterino, Cho and Sperling2000), largely because cyt b studies have reported a lack of phylogenetic signal for comparisons among distant taxa (e.g. Krzywinski et al., Reference Krzywinski, Wilkerson and Besansky2001; Simmons & Weller, Reference Simmons and Weller2001). However, it has been highlighted that cyt b is likely to be more informative within and among very closely related species and populations (Krzywinski et al., Reference Krzywinski, Wilkerson and Besansky2001; Simmons & Weller, Reference Simmons and Weller2001; Drotz, Reference Drotz2003; Monteiro et al., Reference Monteiro, Donnelly, Beard and Costa2004).
We used the cyt b gene to explore the genetic history of this species in order to better understand the current patterns of genetic differentiation and geographical distribution. We also aim to determine if the geographical separation of the Iberian Peninsula and North Africa acts as an effective natural barrier between C. barbara populations and discuss the concept of the differentiation of two subspecies within C. barbara.
Material and methods
Field work was carried out during the summer seasons (July–September) of 1995–2001 in 12 localities from the Iberian Peninsula and three from Morocco (table 1). Male Cicada barbara were identified in the field by their calling songs then collected either by hand or with a sweep net and taken to the lab where they were preserved in liquid nitrogen or in 96% ethanol. Mitochondrial DNA sequences from five specimens per locality (except Toledo: N=2) were analysed (table 1).
Table 1. Sampled populations of Cicada barbara with the distribution of cyt b gene haplotypes (364 bp) and genetic diversity among localities.
N, number of specimens used for mtDNA analysis; h, haploptype diversity; π, nucleotide diversity; d, mean number of pairwise differences.
Molecular analysis
DNA was isolated from a small section of the thoracic muscle (frozen specimens) or from one leg (alcohol preserved) for each cicada, using the methods described in Pinto et al. (Reference Pinto, Quartau, Morgan and Hemingway1998) or using the QIAamp DNA mini kit (QIAgen) following the manufacturer's protocol. Primers to amplify a segment of the mitochondrial cytochrome b gene were designed based on sequences obtained directly from purified mtDNA from 30 C. barbara; these were CB-J-11004: 5' CTTACTTAGGGCTCAAC 3' and CB-N-11351: 5' TAAAAAGTATCATTCTGG 3' with a product size of 365 bp (nomenclature based on Simon et al. (Reference Simon, Frati, Beckenbach, Crespi, Liu and Flook1994), numbers refer to the mitochondrial Drosophila yacuba sequence). PCR reactions comprised 25 μl of 1× Taq buffer, 1.25 mm MgCl2, 0.2 mm of dNTPs, 50 pmol of each primer, 1 unit of Taq polymerase (Invitrogen) and approximately 50 ng of genomic DNA; or, instead, the puReTaq Ready-To-Go PCR beads (Amersham Biosciences, UK) were used, just adding the same primers and template as before. Cycle conditions were as follows: 3 min at 95°C, then 30 cycles of 95°C for 30 s, 50°C for 45 s and 72°C for 1 min, and a final extension step of 72°C for 7 min. Microcon-100 microconcentrators (Amicon) or the Geneclean turbo for PCR kit (QIAgen, BIO101) were used to clean the PCR products prior to direct sequencing.
Sequencing reactions using the primers described above were performed using the ABI Big Dye v. 2 with better buffer (WebScientific Ltd.) and sequences were electrophoresed in a Perkin Elmer 3100 DNA sequencer. Both DNA strands were always sequenced.
DNA sequences were initially analysed with Chromas software version 1.4 (Conor McCarthy, Griffith University, Queensland, Australia), then edited with Editseq and aligned by the CLUSTAL method with MegAlign (DNAstar applications, Lasergene, USA). Transitions, transversions and synonymous versus nonsynonymous substitutions between sequences were analysed in MEGA version 2.1 (Kumar et al., Reference Kumar, Tamura, Jakobsen and Nei2001). Modeltest v. 3.5 (Posada & Crandall, Reference Posada and Crandall1998) was used to search for the most likely model for DNA substitution. The hierarchical likelihood ratio test (hLRT) and the Akaike information criterion (AIC) revealed HKY 85 (Hasegawa et al., Reference Hasegawa, Kishino and Yano1985) as the most likely model of DNA evolution for these data. Therefore, a HKY distance matrix between the sequences was computed using PAUP (4.0 beta version; Swofford, Reference Swofford1999). Genetic diversity, nucleotide variation, minimum spanning network based on the number of pairwise differences between all pairs of haplotypes and AMOVA (analysis of molecular variance) were calculated using the software Arlequin Ver. 2.000 (Schneider et al., Reference Schneider, Roessli and Excoffier2000). A medium joining network based on the maximum likelihood connections between haplotypes was produced by Network 4.1.1.2 (Fluxus Technology Ltd). Pairwise genetic distances based on Slatkin's linearised Fst (Slatkin, Reference Slatkin1995) were also determined using Arlequin. AMOVA was used to test different groupings of populations based on the geographical proximity between sampling sites and the pattern of distribution of the haplotypes. The significance of estimates of the variance components for different hierarchical subdivisions (within populations, within groups of populations and among groups) was tested using 1000 non-parametric permutations (Excoffier et al., Reference Excoffier, Smouse and Quattro1992).
Arlequin was also used to estimate the parameters expected under a sudden demographic expansion and to test the goodness-of-fit of mismatch distributions using 1000 bootstrap replicates. DNAsp4 (Rozas et al., Reference Rozas, Sánchez-DelBarrio, Messeguer and Rozas2003) was used to plot mismatch distributions among haplotypes and to determine the confidence intervals and significance of the raggedness statistics using coalescent simulations based on the neutral infinite-sites model; and, together with Arlequin, to perform neutrality tests in order to corroborate the mismatch distributions statistics (Tajima's (Reference Tajima1989) D-test, Fu & Li's (Reference Fu and Li1993) D* and F*, and Fu's (Reference Fu1997) Fs statistics). Fu's (Reference Fu1997) Fs test statistics have been shown to be more powerful for detecting population growth than the mismatch statistics (Fu, Reference Fu1997; Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002; Mes, Reference Mes2003). Fu's Fs assesses population growth (Fu, Reference Fu1997), recognizing the excess of low frequency alleles in an expanding population compared to the number expected in a stationary one. Tajima's D test assumes that, under neutrality, the number of nucleotide differences between sequences from a random sample should be equal to the number of differences between the polymorphic sites; so, significant negative deviations of Tajima's D from zero might reflect population expansions (Tajima, Reference Tajima1989). Fu & Li's D* and F* statistics were used to compare with Fs tests, as the pattern of significance of these tests may distinguish population growth from the effects of background selection; expansion is expected when Fs is significant and D* and F* are not, while the converse indicates selection (Fu, Reference Fu1997).
The ratio between the number of variable sites (S) and mean pairwise differences among sequences (d) was also calculated (S d−1, expansion coefficient). Expanding populations tend to have sequences with higher frequency of unique sites; thus, large S d−1 values may indicate a recent population expansion, while small values indicate relatively constant long-term population sizes.
Assuming population expansion happened, Arlequin generated the expansion time in mutational units (τ), and the time since expansion was calculated based on the estimators described by Rogers & Harpending (Reference Rogers and Harpending1992), Harpending (Reference Harpending1994) and Rogers (Reference Rogers1995) (τ=2ut, where t is the number of generations elapsed between initial population and current population; and u=2 μk, where μ is the mutation rate per million years and k is the length of the sequence).
Results
Six haplotypes from the 364 bp cyt b gene were detected in the 72 cicadas (fig. 1). All localities analysed from the Iberian Peninsula+Ceuta (North Africa) possessed the same haplotype (HCb1); Seville was one exception, here an additional haplotype was found in one individual (HCb4). Neither haplotype was present in Morocco. Both sites in Morocco shared a dominant haplotype (HCb2); but each locality also revealed other exclusive haplotypes, HCb3 was found in Fès and HCb5 and HCb6 in Meknès.
Fig. 1. Distribution and relative frequency of cyt b haplotypes among populations of Cicada barbara in the Iberia Peninsula and Morocco. The volume of the pie chart representing each haplotype is proportional to the total number of specimens showing that particular haplotype.
Genetic variation
The mean proportion for each nucleotide was 30.3%: A; 39.5%: T; 12.4%: G; and 17.9%: C – demonstrating the expected AT bias for invertebrate mtDNA (Simon et al., Reference Simon, Frati, Beckenbach, Crespi, Liu and Flook1994). Six polymorphic sites were found. Substitutions were all transitions, except one A→T transversion, only present in haplotypes found in Morocco. The mean Ti/Tv ratio between all haplotypes was 1.25. Only one non-synonymous substitution occurred in haplotype CBHCb5, where Ile (Isoleucine) was replaced by Thr (Threonine). The sequences encoded 120 amino acids, with leucine the most common (14.2%), followed by phenylalanine (10.8%) as observed in most insects for this gene (Simmons & Weller, Reference Simmons and Weller2001).
The average number of pairwise differences was 0.576 and the average nucleotide diversity, π, was 0.00159. The mean nucleotide diversity among synonymous sites was 0.00423. Haplotype diversity for all data was 0.2758±0.07 (within Morocco, 0.5333±0.18; Iberian Peninsula+Ceuta, 0.0323±0.03).
Meknès had the highest haplotype (h) and nucleotide (π) diversity (h=0.7 and π=0.0022; table 1). Both Fès and Seville had two haplotypes and the same values for haplotype and nucleotide diversity (h=0.4 and π=0.0011). The mean number of pairwise differences within these localities was 0.4 for Fès and Seville, and 0.8 for Meknès. All remaining localities, including Ceuta, did not show any diversity.
Population structure
A minimum spanning network constructed by hand, based on the number of pairwise nucleotide differences computed using Arlequin and a medium joining network computed by Network based on maximum likelihood, are shown in fig. 2. Even though the number of haplotypes and the connection lengths are low, Moroccan haplotypes group monophyletically.
Fig. 2. Medium joining and minimum spanning networks for Cicada barbara in the Iberia Peninsula and Morocco based on cyt b haplotypes. The volume of the circles representing each haplotype is proportional to the total number of specimens showing that particular haplotype. The colours are representative of the number of specimens per locality. Each connection indicates one mutational step and black dots in the medium joining network represent missing intermediate haplotypes.
HKY distances between all pairwise combinations of haplotypes varied from 0.0028 to 0.0084 (table 2). The average genetic distance among haplotypes within Morocco was 0.42%, within the Iberian Peninsula+Ceuta, 0.28% and between the two groups the average HKY distance was 0.63%.
Table 2. HKY pairwise distances between all pairs of haplotypes based in a 364 bp fragment of the Cyt b gene within Cicada barbara.
Slatkin's linearised Fst revealed significant genetic differentiation (P<0.05) in all pairwise comparisons between Morocco samples (Fès and Mèknes) versus all other populations (table 3), with values ranging from 0.6 to 0.8. All other sample pairwise comparisons revealed zero or non-significant Fst values (samples within the Iberia Peninsula and also between Fès and Meknès).
Table 3. Slatkin's linearised Fst values (Slatkin, Reference Slatkin1995) between all pairs of samples based in a 364 bp fragment of the Cyt b gene within Cicada barbara, below diagonal.
* , significant values (P<0.05).
Several hierarchical analyses of molecular variance were carried out to compare hypotheses of genetic partitioning. The highest ΦCT value, which was statistically significant, was found when samples were separated into two groups (table 4), one including Iberian Peninsula+Ceuta samples and the other including Moroccan samples (ΦCT 0.91). Other partitioning tests with lower ΦCT were also significant. For instance, comparing Iberian Peninsula vs. North Africa (Ceuta and Morocco), a ΦCT value of 0.78 was found. When testing five groups (North Portugal (Foz Côa), South Portugal, South Spain, North Africa (Ceuta) and Morocco), the ΦCT was even lower (0.67).
Table 4. AMOVA results for the two groups of Cicada barbara maximizing the ΦCT values analysed for cyt b gene (Morocco vs Iberian Peninsula and Ceuta).
d.f., degrees of freedom; P, significance of percentage variation and fixation indices estimated from permutation tests (1000 permutations); ΦCT, correlation of random haplotypes within a group of populations relative to that drawn from the entire species;ΦSC, correlation of random haplotypes within populations relative to that within a regional grouping of populations within a region;ΦST, correlation of random haplotypes within a population relative to that from the whole species (Excoffier et al., Reference Excoffier, Smouse and Quattro1992); significant+P<0.05.
Demographic analysis
Since there was a significant genetic differentiation between Morocco and the Iberian Peninsula+Ceuta, the historical demography of these two regions was analysed separately. The distribution of pairwise nucleotide differences for the Iberia Peninsula+Ceuta samples and for the Morocco samples alone were computed and plotted (fig. 3). The Iberian Peninsula+Ceuta group showed a unimodal distribution, while Morocco revealed a slightly bimodal curve, even though both distributions showed the expected pattern for populations that have undergone sudden expansion (table 5). The observed raggedness was high for the Iberia Peninsula+Ceuta and moderate for Morocco, though not significant in both cases. This was most likely due to the small number of haplotypes observed for the first region. The coalescent simulated raggedness, on the other hand, was notably low for the Iberia Peninsula+Ceuta and, again, moderate for Morocco.
Fig. 3. Mismatch distributions within Cicada barbara (—, Exp; –-o–-, Obs).
Table 5. Diversity, neutrality and expansion parameters for Cicada barbara based on a section of cyt b gene.
S, polymorphic sites; d, mean number of pairwise nucleotide differences; CI, confidence interval. Significant values for the neutrality estimations (P<0.05) are accentuated with (+) Estimated mean τ parameter was obtained according to the percentile method from the sudden expansion model (confidence level: α=0.05 based on 1000 replicates).
Both regions showed significantly negative Fu's Fs values (table 5) indicative of expanding populations. In the Morocco samples, Tajima's D was significantly negative, and Fu and Li's D* and F* were non-significant negative, which, again, reflects population growth and no background selection. Iberia Peninsula+Ceuta also had negative values for these statistics, but Tajima's D was not significant contrary to Fu and Li's D* and F*.
The expansion coefficients estimated for each region revealed high values for the two regions analysed; however, Iberia Peninsula+Ceuta had a much higher expansion coefficient than Morocco.
The τ parameter for the Iberian samples was considerably higher than for the Moroccan samples. Cicada nymphs develop underground and the exact time they spend at this stage is still unknown for most species within genus Cicada. In fact, the cicada immature stages may last between two to three years (Quartau, Reference Quartau1995; Boulard & Mondon, Reference Boulard and Mondon1996; Giralda et al., Reference Giralda, Rodriguez and Pintor1998). Considering this, there are two possible scenarios for the expansion time. Considering a mutation rate for the cyt b gene of 2.15% (as for Agabus species, Coleoptera: Drotz, Reference Drotz2003) and one generation every two years, the time of expansion (t) based in the estimators described by Rogers & Harpending (Reference Rogers and Harpending1992) and Harpending (Reference Harpending1994) for Iberian Peninsula cicadas occurred around 191,668 years BP (20,444–191,668 years BP), while Moroccan cicadas are estimated to have expanded about 48,492 years ago (0–113,724 years BP). The estimated time of expansion for the Iberian Peninsula is also at the maximum of its error range, most likely an artefact due to the low number of haplotypes detected. If we consider that these cicadas have a generation time of three years; then, for the Iberian Peninsula, the expansion is estimated to have occurred around 287,502 years BP (30,666–287,502 years BP) and for Morocco about 72,738 years BP (0–170,586 years BP).
Discussion
Population structure and genetic variation
The divergence between Moroccan and the Iberian Peninsula+Ceuta populations is strongly corroborated by the molecular data. The AMOVA results, for the cyt b sequence analysed for C. barbara, revealed a highly significant ΦCT value of 0.91; and the high mean pairwise ΦST value (0.90) implies genetically isolated populations with a low level of gene flow.
The Ceuta population from Spanish North Africa was more similar to the Iberian populations than to the surrounding Moroccan populations. This relationship might suggest that the Strait of Gibraltar does not act as a strict barrier to dispersal. The narrowest distance between the southern-most point of Spain and the northwestern-most point of North Africa is 13 km (the strait ranges from 12.9 km to 43 km), possibly within potential flight distances for cicadas. Moreover, the drop in sea level during the last glaciations, revealing the European and the African continental shelves and several islands in between (Collina-Girard, Reference Collina-Girard2001) would have facilitated the movement of cicadas through the Strait of Gibraltar. Wind dispersal over water has been invoked as one of the explanations for dispersal in New Zealand cicadas (Arensburger et al., Reference Arensburger, Simon and Holsinger2004) and sea-surface dispersal through the Strait of Gibraltar by some coleoptera, including wingless insects, has been postulated (e.g. Palmer & Cambefort, Reference Palmer and Cambefort2000; Horn et al., Reference Horn, Roux-Morabito, Lieutier and Kerdelhue2006). Alternatively, considering that Ceuta is a Spanish territory, a possible human-induced translocation should not be ignored. Trees or shrubs, containing C. barbara specimens (eggs or adults), from the Iberian Peninsula might have been introduced. In contrast, the high altitude of the Rif Mountains (highest point is 2455 m) to the south of Ceuta and its climatic extremes, including snow and frost, might provide a strong obstacle to the natural dispersal of these thermophilic cicadas. There is no record of C. barbara in the Rif Mountains and, considering the distribution of vegetation and climate along the altitude gradient here (e.g. Moore et al., Reference Moore, Fox, Harrouni and El Alami1998), it is highly unlikely that C. barbara occurs above 1000 m. The cultivation of olive trees in these mountains occurs below 500 m (Moore et al., Reference Moore, Fox, Harrouni and El Alami1998). The presence of a rock glacier site in these mountains is also indicative of the severe past climatic conditions. According to Hughes et al. (Reference Hughes, Woodward and Gibbard2006), glacial deposits extend to exceptionally low altitudes for this latitude (ca. 36°N), reaching as low as 750 m on some slopes of the massif. The Rif Mountains are known to hold a high rate of endemism, partly as a result of the long isolation of Holartic taxa in the high elevations of these North African mountain ranges. The HKY genetic distances observed between the haplotypes of C. barbara were relatively low (0.28–0.84%). For the same gene, Drotz (Reference Drotz2003) found genetic distances between 0–3.9% within Agabus species (Coleoptera) and Monteiro et al. (Reference Monteiro, Donnelly, Beard and Costa2004) referred distances of 0.7–1.2% within groups of Triatoma brasiliensis (Hemiptera).
Demographic expansions
Overall, the demographic analysis suggests sudden expansions for both regions. Even though the observed raggedness was high, especially for the Iberian Peninsula, coalescent simulations revealed much smaller values, which is probably due to the low number of haplotypes observed in the Iberian Peninsula. All other estimators analysed, including those usually considered as more robust and reliable, strongly supported sudden expansion models. However, the low number of haplotypes and low nucleotide diversity in the Iberian Peninsula could indicate the occurrence of a strong bottleneck at some time in the past.
The Iberian Peninsula and Ceuta specimens showed a more ancient signature of demographic expansion than in Morocco. These results should, however, be treated with caution because the mutation rate used for cyt b gene might not be appropriate for C. barbara, and the low number of haplotypes detected weakens these tests. Also, the generation time has yet to be accurately determined for this species; and the sample sizes, especially for Morocco, are relatively low.
Cicada barbara can be seen on different kinds of shrubs or trees; but, in the majority of sites, they are found associated with olive and pine. Therefore, it is most likely that the life history of cicadas might be related with the life history of these trees. In fact, Burban & Petit (Reference Burban and Petit2003) found a similar pattern of genetic differentiation within maritime pine (Pinus pinaster). Three different mitotypes were detected in maritime pine using RFLP analysis within the Mediterranean area: one haplotype was exclusively present in the Iberian Peninsula and in Punta Cires (in the north of Morocco very close to Ceuta), one other was present only in Morocco and the third was present in Eastern Europe and Northeast Africa. On the other hand, the cultivation of olive trees is known to have happened since ancient times, and humans have greatly influenced its dispersal (Besnard et al., Reference Besnard, Khadari, Baradat and Bervillé2002a). As such, the genetic structure of olive trees reflects not only the occurrence of refugial zones but also the particular biogeographic conditions of the Mediterranean Basin and the human influence (Besnard et al., Reference Besnard, Khadari, Baradat and Bervillé2002a). Nevertheless, strong genetic differentiation between the eastern and western Mediterranean, as well as for Morocco, was shown for olive cpDNA (Besnard et al., Reference Besnard, Khadari, Baradat and Bervillé2002a); and high levels of genetic differentiation were revealed by both cpDNA and mtDNA, probably due to the contraction of the distribution during the Quaternary glaciations and subsequent re-colonization during Holocene (Besnard et al., Reference Besnard, Khadari, Baradat and Bervillé2002b). Even though no eastern samples were available for this study, these data show that in the Iberia Peninsula and Morocco C. barbara has a similar pattern of genetic differentiation as both pine and olive trees.
Iberia Peninsula expansions
Whether we consider that cicadas have a two or three year generation time, the expansion time for the Iberian Peninsula cicadas coincides with the Mindel-Riss interglacial period (300,000–180,000 years BP). Furthermore, the previous (Mindel) glaciation (410,000–380,000 years BP) correlates with the strong bottleneck we detected for Iberian cicadas. In fact, it is considered that this glaciation marks the absolute maximum extent of continental ice sheets in the Quaternary (Lindner et al., Reference Lindner, Dzierzek, Marciniak and Nitychoruk2003). During the subsequent glaciations (Riss, 180,000–130,000 years BP; and Würm, approximately 70,000–10,000 years BP), C. barbara individuals could have potentially survived in refugia, which persisted in the Iberia Peninsula. Moreover, the last interglacial, the Eemian, which started approximately 126,000 years BP (Sánchez Goñi et al., Reference Sánchez Goñi, Turon, Eynaud, Schackleton and Cayre2000) allowed Mediterranean taxa to expand. For instance, the maritime pine in the Iberian Peninsula seems to have been largely favoured in comparison with other Mediterranean conifers due to its wide ecological range (Gómez et al., Reference Gómez, Vendramin, González-Martínez and Alía2005). During glacial periods, P. pinaster was found near the coast and mountains, and it has been suggested to have a continuous presence at lower altitudes since 30,000 years BP (Carrión & Van Geel (1999) inGómez et al., Reference Gómez, Vendramin, González-Martínez and Alía2005). Nevertheless, C. barbara did not show any geographical structure in Iberia, unlike other Mediterranean species reported (e.g. Paulo et al., Reference Paulo, Dias, Bruford, Jordan and Nichols2001; Jansson & Dynesius, Reference Jansson and Dynesius2002; Hewitt, Reference Hewitt2004b), suggesting high vagility and vulnerability to climate change. The extreme northwest of Africa, such as Ceuta, might, however, have acted as a southern refuge for Iberian cicadas during the most severe climatic conditions. Specimens of C. barbara from the Iberian Peninsula and from Ceuta (North Africa) did not show any genetic differentiation; and, as stated before, access between the two margins of the Strait of Gibraltar would have been easier during the Pleistocene.
Moroccan expansions
The estimated expansion time for Moroccan populations, considering a two year generation time, occurred during the last Pleistocene glaciation, which ended around 14,000 years BP. On the other hand, considering a three year generation time, Moroccan cicadas would have expanded after Riss glaciation. Since the end of the Riss glaciation, during the interglacial, the Mediterranean oak forests expanded greatly in Morocco (Hooghiemstra et al., Reference Hooghiemstra, Stalling, Agwu and Dupont1992), allowing the expansion of other species, e.g. macaques (Modolo et al., Reference Modolo, Salzburger and Martin2005). Moroccan cicadas might have expanded during this phase and retreated when conditions were not so favourable and finally expanded again at the end of the last ice age, when climatic conditions were improving. In fact, the mismatch distribution of Moroccan specimens shows two peaks, and the signal of the expansion is not as clear as for the Iberian Peninsula populations. This might indicate that populations in Morocco were more stable; and, even though they suffered retreats and expansions, there is no evidence of a very strong bottleneck during these ice ages. Also, the high haplotype variability found in these populations suggests a higher stability than in the Iberian Peninsula populations. These results are in agreement with other studies of Mediterranean taxa suggesting that the most southern populations show higher geographical substructure and stability (e.g. Jansson & Dynesius, Reference Jansson and Dynesius2002; Hewitt, Reference Hewitt2004b; Pinho et al., Reference Pinho, Harris and Ferrand2007).
Origin of Iberian cicadas and subspecies classification
The hypothesis of C. barbara being a recent immigrant to the Iberian Peninsula from North Africa (Quartau et al., Reference Quartau, Ribeiro, Simões and Coelho2001) cannot be confirmed or refuted with these data. Ceuta, one of the most northwestern localities in North Africa, did not clarify this issue, as it showed the same genetic pattern as Iberian populations, but remains significant since it potentially shows that the Strait of Gibraltar does not seem to have prevented dispersal of these cicadas in contrast to the Rif Mountains. Although the low genetic variability of the Iberian populations might suggest a very recent colonization, the complete absence of shared haplotypes with Morocco seems to indicate an ancient separation between these regions. Also, the expansion signal in the Iberian Peninsula is much more ancient than in Morocco. It is possible that a sampling artefact might have hidden shared haplotypes; however, these results do not seem to suggest a Moroccan origin of the Iberian cicadas. Nevertheless, we cannot exclude the hypothesis of a North African origin. In fact, Morocco topography, with the Rif Mountains in the north, the Atlas Mountains in the southeast and the Sahara desert further south, favours the isolation of species (Modolo et al., Reference Modolo, Salzburger and Martin2005) and consequently their divergence, as also found with the maritime pine (Burban & Petit, Reference Burban and Petit2003) and the olive, Olea europaea subsp. maroccana (Besnard et al., Reference Besnard, Khadari, Baradat and Bervillé2002a,Reference Besnard, Khadari, Baradat and Bervilléb). Also, during the Quaternary climatic fluctuations, the Sahara Desert expanded and retracted repeatedly (Cox & Moore, Reference Cox and Moore2005), separating a northwestern and a northeastern refuge which favoured the differentiation between populations in northwest and northeast Africa, as observed in some coleoptera, snails, lizards, freshwater turtles, shrews and snakes (e.g. Cosson et al., Reference Cosson, Hutterer, Libois, Sarà, Taberlet and Vogel2005; Fritz et al., Reference Fritz, Barata, Busack, Fritzsch and Castilho2005; Horn et al., Reference Horn, Roux-Morabito, Lieutier and Kerdelhue2006; Carranza et al., Reference Carranza, Arnold and Pleguezuelos2006). It will be important to sample cicadas from other North African regions to clarify this.
The definition of the two subspecies within C. barbara (C. barbara lusitanica and C. barbara barbara) finds some support here, although Moroccan specimens might be also differentiated from the other Northern African individuals as found in other taxa. Species delimitation remains conceptually difficult and for subspecies the subjectivity may be even more pronounced. Davies & Nixon (Reference Davies and Nixon1992) suggested that evolutionary lineages may become different species if at least one character state is fixed in one group and absent in the other. Templeton (Reference Templeton, Schierwater, Streit, Wagner and DeSalle1994) stated that evolutionary lineages must be phylogenetically distinct, show no recent gene flow and reveal ecological or demographic limitations to reproduction. Zink (Reference Zink2004) and Phillimore & Owens (Reference Phillimore and Owens2006) also used monophyly as an indicator of phylogenetic distinctiveness of subspecies. On the other hand, Lee (Reference Lee2004) emphasized the importance of analysing other biological information besides molecular data to define species, and pointed out that it is not possible to associate a specific level of genetic divergence to reproductive isolation. Even so, some of these concepts might suggest that the two groups of cicadas analysed here are differentiating and may be on the way to splitting, as supported by some acoustic differentiation which is involved in mating attraction (see Pinto-Juma et al., this journal).
Acknowledgements
We thank Genage André, Sofia Seabra, Paula Simões, Teresa Rebelo, Mónica Ribeiro and Teresa Fernandes (Faculty of Science, University of Lisbon) for their help with field work. We also thank G. Hewitt and an anonymous reviewer for their comments and suggestions on a previous version of this manuscript. This research was supported by the PhD grant PRAXIS BD/18229/98 from the Fundação para a Ciência e Tecnologia (Lisbon, Portugal).
