Evolutionary history and patterns of differentiation among European Maniola butterflies ( Lepidoptera : Satyrinae )

Phylogenetic relationships of European Maniola butterflies are reconstructed using molecular sequences from two regions of the mitochondrial DNA, cytochrome oxidase subunit I (COI) and cytochrome b (Cytb). A total of 988 base pairs (486 for Cytb, and 502 for COI) were aligned for 15 individuals of Maniola and an outgroup species. The phylogenetic tree obtained through Bayesian inference analysis of the combined data sets shows evidence that the island endemic M. chia is indistinguishable from M. jurtina on the basis of the mtDNA genes studied. Net nucleotide divergence between M. jurtina and M. chia is 0.4%, but 2% between the M. jurtina and the M. nurag clade. A phenetically distinct entity of individuals from Sardinia appears to be a hybrid between M. nurag and M. jurtina. The southern and northern European ecotypes of M. jurtina, which differ in the summer aestivation period of the southern type, are not structured genetically at the level of coding mtDNA genes. Divergence time between M. nurag and M. jurtina was estimated to be 1.1 to 1.2 million years. Speciation most likely took place in the early Pleistocene as a consequence of the isolation of Sardinia, when the sea reflooded the Mediterranean basin after the Messinian crisis (about 5 million years ago).


INTRODUCTION
The Palaearctic genus Maniola (Satyrinae: Maniolini) includes one of the most widespread and common European butterflies, the Meadow Brown, Maniola jurtina (L.), and six less widespread and largely vicariant species among which three of Europe's rarest butterflies, the island endemics M. chia (Thomson, 1987) on the Greek island of Chios, M. cypricola (Graves, 1928) on Cyprus, and M. nurag (Ghiliani, 1852) on Sardinia.Endemic species often replace more widespread congeners.On Sardinia, however, M. jurtina is sympatric with the endemic M. nurag.M. jurtina is replaced by M. telmessia (Zeller, 1847) in southern and western Turkey, and eastwards to the Bosporus, by M. halicarnassus (Thomson, 1990) in the Bodrum peninsula and the Aegean island of Nissiros, and by M. megala (Oberthür, 1909) in southern Turkey and eastwards to Iran.
Maniola jurtina is morphologically and ecologically the most variable species in the genus (Ford, 1945;Thomson, 1971Thomson, , 1973;;Brakefield, 1982).Intraspecific morphological variation includes pronounced wing-pattern polymorphism, as well as ecological variability related to latitude and altitude (Shreeve, 1989;Goulson, 1993a, b).In the Southern areas of the Palaearctic, adult females perform a diapause during the hottest part of the summer, with a concomitantly delayed ovarian maturation (Masetti & Scali, 1972;García-Barros, 1987), while in the Northern Palaearctic eggs are deposited immediately after fertilization.
A priori, such a high degree of ecological variability within a single species could suggest the existence of cryptic species (Thomson, 1971;Goulson, 1993a), as has been shown for a number of butterflies (e.g.Beneš et al., 2003;Hebert et al., 2004).In fact, a distinct morphotype within the genus Maniola has recently been found on Sardinia (Grill et al., 2004).These individuals occur in areas, where M. nurag and M. jurtina are sympatric, and have a wing pattern that appears intermediate between the other two species.In addition, they are also distinct in genitalia morphology (Grill et al., 2004).The Greek island endemic, M. chia, on the other hand, is identical to M. jurtina in wing patterns and can only be distinguished by the male genitalia (Grill et al., 2004), making it apparently more similar to M. jurtina than the Sardinian morphotype.
Despite their ubiquitousness, Maniola butterflies have as yet remained little studied in phylogenetic contexts.The most recent phylogeny of European satyrids includes only two individuals of M. jurtina from France, the other species in the genus are still missing (Martin et al., 2000).
Here, we determined sequences of the mitochondrial genes cytochrome oxidase subunit I (COI) and cytochrome b (Cytb) to infer patterns of genetic differentiation among European Maniola butterflies.The three main aims of this study are (1) to clarify the systematic position of M. chia within the genus Maniola, (2) to investigate if the Sardinian morphotypes are genetically differentiated from the Maniola with which they occur in sympatry, (3) to test our hypothesis that the ecological variability of M. jurtina might result from the existence of various cryptic Maniola species.Further, we estimate the time scale of the speciation events between M. jurtina and the island endemics M. nurag and M. chia.

Sample collection
We sampled a total of 16 individuals, comprising specimens of the following morphotypes: the distinct Sardinian individuals (further termed "intermediate"), each of the two different ecological races of M. jurtina and its morphologically most similar sister species, M. chia and M. nurag (Table 1).These three species are morphologically very similar and can often only be distinguished on the basis of the shape of the male genitalia (Grill et al., 2004), while the remainder of the species in this genus can readily be classified on the basis of wing patterns.Individuals were collected from Sardinian and continental populations of M. jurtina, representing the two different ecological races, a M. chia population from Chios and a M. nurag population, as well as four individuals of the morphologically distinct group from Sardinia (Table 1).Pyronia cecilia was collected to serve as an outgroup in the phylogenetic reconstruction.

Molecular analysis
DNA was extracted from the legs of freshly frozen butterflies, using the following protocol, modified from Aljanabi & Martinez (1997).Homogenization was carried out in the presence of 200 µl of buffer: 50 mM Tris-HCl pH 8.0, 20 mM EDTA, 2% SDS.The material was grounded and 2.5 µl (20 mg/ml) proteinase-K was added.The mix was incubated at 37°C for 1-2 h and 50 µl of NaCl (360 g NaCl / 1 l H2O) was added to the samples.This mixture was then microcentrifuged for 15 min at 14,000 rpm.The supernatant was transferred to new eppendorf tubes and the DNA precipitated by using 0.5 ml of ice-cold 99% ethanol.This mixture was again microcentrifuged for 30 min at 14,000 rpm.This was followed by the first 70% ethanol wash.After another microcentrifugion at maximum speed at 4°C for 15 min the fluid was removed and the DNA was dried in a vacuum desiccator for 15-30 min at room temperature.The pellet was resuspended in 20 µl of HPLC-purified dH2O.The samples were incubated at 37°C for 30 min and then stored at -20°C.
All PCR reactions were performed in a 50 µl reaction volume.COI PCR reaction was prepared according to the following protocol: 33.9 µl H2O, 5 µl PCR buffer, 5 µl dNTP's, 1 µl of MgCl (50 mg/ml), 2 µl BSA (10 mg/ml), 0.5 µl of each primer, 0.4 µl of Taq polymerase, and 2 µl of DNA.The cycling profile for COI was 94°C for 3 min, 35 cycles of 94°C for 45 s, 46°C for 30 s and 72°C for 1 min, and a final extension period of 72°C for 5 min.For Cytb, 25.6 µl H2O, 5 µl PCR buffer, 10 µl dNTP's, 3 µl of MgCl (50 mg/ml), 3 µl BSA (10 mg/ml), 0.5 µl of each primer, 0.4 µl of Taq polymerase, and 2 µl of DNA were used.The cycling profile was 94°C for 3 min, 38 cycles of 94°C for 1 min, 51°C for 1 min and 72°C for 1 min, and a final extension period of 72°C for 5 min (Palumbi, 1996).The primers used to amplify the COI fragment were C1-J-1718: GGA GGA TTT GGA AAT TGA TTA GTT CC and C1-N-2191: CCC GGT AAA ATT AAA ATA TAA ACT TC, and for the Cytb CB-J-10933: TAT GTA CTA CCA TGA GGA CAA ATA TC and CB-N-11367: ATT ACA CCT CCT AAT TTA TTA GGA AT (Simon et al., 1994).Purified products were sequenced using an Applied Biosystems ABI 3700 automated sequencer in Laboratories of the Academische Medisch Centrum, Universiteit Van Amsterdam.Since all sequenced fragments were of the same size, alignment was trivial.BLAST searches were conducted on all sequences to check for possible contamination.

Phylogenetic analysis
Homogeneity of partitions between the two data sets (COI and Cytb) was tested by analyzing the incongruence length difference (ILD; Farris et al., 1994) between trees obtained separately for each data set, using PAUP* 4.0.b10(Swofford, 2003).Maximum Parsimony analyses were performed on the whole data set using PAUP, with heuristic search option (TBR, random sequence addition, MaxTrees set to 500) with 1,000 randomaddition replicates.Partitioned Bayesian analyses were performed on the combined data set using MrBayes version 3.0b4 (Huelsenbeck & Ronquist, 2001).For each partition (COI and Cytb), we used Modeltest 3.06 (Posada & Crandall, 1998) to assess the best-fit substitution model, through hierarchical likelihood ratio tests.The asymptote of the fluctuating likelihood values of the Bayesian trees (or burnin period) was determined in preliminary runs.We ran four Metropolis-coupled chains in one run of 500,000 generations, and sampled one tree every 1000 cycles after the burnin period had passed.The sampled trees were used to generate a majority rule tree and the support for the nodes of this tree was given by posterior probability estimates for each clade.The hypothesis of equal mutation rates and substitution rates among lineages was tested by performing a Likelihood Ratio test comparing the likelihood scores of Maximum Likelihood trees either unconstrained or with a molecular clock enforced in the two Phylip 3.6b package programs (Felsenstein, 2004), DNAML and DNAMLK (less timeconsuming than under PAUP*).Net distances between groups and average distances within groups were calculated using Mega 3 (Kumar et al., 2004), with the Tamura-Nei model of substitution.

RESULTS
A total of 988 base pairs (486 for Cytb, and 502 for COI) were aligned for 15 individuals of Maniola and one outgroup species, Pyronia cecilia (Table 1).Among the 988 nucleotides of our dataset, 32 sites were variable and 29 were parsimony-informative (without considering the outgroup species).The total number of distinct haplotypes was 12.We computed Bayesian inferences using the following prior probability parameters determined by Modeltest: for Cytb, best-fit model of sequence evolution was the Hasegawa-Kishino-Yano model of substitution, rates "gamma"; for COI, the best-fit model of sequence evolution was the Tamura-Nei model of substitution, rates "equal".After 50,000 cycles the chain converged to stable likelihood values.Four hundred and fifty trees were sampled to reconstruct the majority rule tree with posterior probability estimates.The Maximum Parsimony analysis resulted in three equally parsimonious trees (140 steps; CI = 0.8929; RI = 0.9045), which only slightly differed in the level of resolution of the M. nurag clade (trees not shown).The 50% majority-rule consensus tree (see bootstrap values in Fig. 1) resulted in a tree similar to the Bayesian inference reconstruction, with all nodes in The two studied genes (COI and Cytb) showed differences in the best-fit model of evolution, as attested by Modeltest.Homogeneity between the two partitions was rejected when tested for incongruence length differences (P = 0.0222).However, since the two genes are physically linked in the mitochondrion, we considered them a part of one single and consistent evolutionary unit and tested the molecular clock hypothesis on the whole data set.Trees obtained without any constraint were congruent with those obtained under the molecular clock model of sequence evolution.Values of -Ln (likelihood) for "molecular clock-unconstrained" vs. "molecular clockconstrained trees" were 1987.211and 1973.499,respectively (dl = 14, P = 0.471), supporting the use of the molecular clock model.Net nucleotide divergence for COI and Cytb between M. jurtina (including intermediates) and M. chia was 0.4%.Net nucleotide divergence for COI and Cytb between the M. nurag clade (including intermediates) and the M. jurtina clade (including intermediates and M. chia individuals) was 2.0%.The withingroup divergence of the two clades was 0.9% for the jurtina-clade and 0.1% for the nurag-clade.Average nucleotide divergence between the northern and southern ecotype of M. jurtina is 0.8%, whereas the within-group divergence is respectively 1.4% (northern individuals) and 0.4% (southern individuals).Based on the assumption of a mutation rate of 1.1% to 1.2% per million years in arthropod mitochondrial DNA (Brower, 1994;Gaunt & Miles, 2002), which translates to about 2.3% sequence divergence per million years, divergence time between M. nurag and M. jurtina was estimated to be between 1.1 and 1.2 million years.

Genetic identity of the Sardinian morpho-species
The individuals that are morphologically intermediate between M. nurag and M. jurtina branch either in the jurtina-clade or in the nurag-clade, suggesting that this third group of Sardinian Maniola is the result of hybridisation between M. nurag and M. jurtina.This result confirms the pattern obtained from allozyme data, in which intermediate individuals clustered together on a separate branch, equally distant from the jurtina-clade and from the nurag-clade (Grill et al., in prep.).

Status of the different ecotypes of M. jurtina
The two ecological races of M. jurtina cluster in the same clade and are not differentiated from each other [i.e., the divergence between the two individuals from northern Europe is larger than the mean divergence between individuals from southern and northern Europe (Fig. 1)].Therefore, and despite the important ecological and morphological differences between northern and southern European populations, M. jurtina is probably not a complex of cryptic species, or if it were, such species would have diverged in the late Pleistocene.In this case, differentiation can only be detected by analysing more variable genes.The fact that genetic distances within the jurtina-clade are larger than within the nurag-clade may reflect the larger geographic area across which the individuals in the jurtina-clade were sampled (while all individuals in the nurag-clade came from Sardinia).

Genetic identity of M. chia
Our results show that the Maniola individuals from the island of Chios are not more differentiated from the continental M. jurtina than one M. jurtina individual from another.Despite the morphological differences in the genitalia of M. chia, on the basis of which the species was originally described (Thomson, 1987), it is genetically very close to M. jurtina.Moreover, the two species are indistinguishable by their wing patterns.Also ecologically, M. chia and M. jurtina are very similar.They occur in the same kind of habitats, use the same nectar plants, and both species are able to undergo a summer-diapause (A.Grill, pers. observ.).
As the island of Chios is only a few hundred meters from the Turkish mainland, a distance that can easily be overcome by M. jurtina (e.g.Brakefield, 1982;Conradt et al., 2000), a certain degree of gene-flow between mainland and island populations can not be excluded.

Evolutionary history of M. nurag
The estimated divergence time for the split between M. nurag and M. jurtina is in congruence with divergence time estimates based on allozyme data (Grill et al., submitted) as well as with data available on other endemic Sardinian taxa (Fig. 2).
The time scale of these speciation events can be explained with regard to the geological history of the island of Sardinia.About 5 million years ago the Mediterranean basin was almost entirely desiccated, a period referred to as the Messinian salinity crisis (Hsü et al., 1977;Por, 1989).This happened after tectonic action lifted the Strait of Gibraltar above sea level, so that the Mediterranean Sea dried out and became a desert basin with salt lakes at the bottom (Grove & Rackham, 2003).When the connection to the Atlantic was re-established, the massive influx of water brought a dramatic change in marine flora and fauna, and at the same time induced the differentiation of many island species, which then became isolated from the continent.A number of Sardinian lineages in various organismal groups, for example butterflies (Marchi et al., 1996), cave-beetles (Caccone & Sbordoni, 2001) and lizards (Oliviero et al., 1998) differentiated after the Messinian crisis.Until well into the Pleistocene, severe sea level oscillations (up to 1000 m) continuously reconnected Sardinia with northern Italy and southern France, creating land bridges between the Tyrrhenian islands and the continent (Arias et al., 1980).
During the last glacial maximum (20 000 years ago), Sardinia was connected with neighbouring Corsica.Although the two islands share a number of endemic butterflies (Biermann, 1998), M. nurag is absent from Corsica, while M. jurtina flies there widespread and common.Our divergence time estimate suggests that M. nurag is the result of vicariance and differentiated when Sardinia was disconnected from the continental mainland after a rise of the sea level.During the last glacial maximum, Sardinia was reconnected with Corsica.In that period, M. jurtina could have colonized Sardinia.During the postglacial warming, M. nurag retreated to higher altitudes, whereas the Sardinian M. jurtina remained mostly at sea level, resulting in the distribution pattern we observe today, where M. nurag is restricted to areas above 500 m while M. jurtina has its largest populations at sea level, and is only occasionally found at altitudes above 1000 m.
In the zone where the distributions of both species overlap, occasional hybridisation between the two species occurs, resulting in morphologically intermediate individuals, like those investigated here.Similar hybridisation events have been detected between Papilio machaon and P. hospiton, a papilionid butterfly endemic to Sardo-Corsica (Cianchi et al., 2003).

CONCLUSIONS
(I) The distinct Sardinian morphotype is a hybrid between M. nurag and M. jurtina.As these individuals are very low frequent in the populations (<5%), we suppose that they are in the majority F1 hybrids demonstrating (i) low levels of fertility, or (ii) backcross reproduction with one of the parent species.
(II) In the light of these results, it seems appropriate to include the Maniola individuals coming from Chios within the taxonomic entity of M. jurtina.Considering the probably limited gene flow between the Maniola from Chios and their continental conspecifics, however, it is justifiable to treat the former as an "evolutionary significant unit" in nature conservation contexts.
(III) Both Maniola nurag and M. jurtina are "true" species, which seem to have diverged in the early Pleistocene about 1 million years ago.This divergence is likely to be the consequence of the isolation of Sardinia after to the Messinian crisis.However, the two species are not reproductively isolated, as attested by the existence of hybrids.Ecotypes within M. jurtina are not structured genetically at the level of coding mtDNA genes.
Needless to say, further genetic work will be necessary to fully understand the evolutionary history of the species in the genus Maniola, as well as their position in the phylogeny of the Maniolini.

Fig. 1 .
Fig. 1.Phylogenetic relationships of the meadow brown butterfly, Maniola jurtina, and its endemic sister taxa, M. nurag and M. chia.The "intermediate" individuals are the result of hybridisation between M. jurtina and M. nurag.The phylogenetic tree was obtained through Bayesian inference analysis of the combined data set, using Bayesian estimates of branch lengths.Numbers adjacent to nodes give Bayesian posterior probabilities of nodes greater than 50% (in bold) and Maximum Parsimony bootstrap values in the 50% majority-rule consensus tree (in italics).

Fig. 2 .
Fig. 2. Schematic presentation of time scale of the geophysical events in relation to speciation time estimates of Sardinian taxa inferred from molecular data (see main text for references).

TABLE 1 .
List of specimens, with indications of species distribution, samling locality, sex of individuals, voucher number and accession numbers for Cytb and COI.Abbreviations for sampling localities are as follows: NL = Netherlands, SA = Sardinia, Italy, CH = Chios, Greece.