Tracing the radiation of Maniola (Nymphalidae) butterflies: new insights from phylogeography hint at one single incompletely differentiated species complex

The use of DNA sequence data often leads to the recognition of cryptic species within putatively well-known taxa. The opposite case, detecting less diversity than originally described, has, however, far more rarely been documented. Maniola jurtina, the Meadow Brown butterfly, occurs all over Europe, whereas all other six species in the genus Maniola are restricted to the Mediterranean area. Among them, three are island endemics on Sardinia, Cyprus, and Chios, respectively. Maniola species are almost indistinguishable morphologically, and hybridization seems to occur occasionally. To clarify species boundaries and diversification history of the genus, we reconstructed the phylogeography and phylogeny of all seven species within Maniola analyzing 138 individuals from across its range using mitochondrial and nuclear genetic markers. Examination of variation in mitochondrial and nuclear DNA surprisingly revealed a case of taxonomic “oversplitting”. The topology of the recovered phylogenetic tree is not consistent with accepted taxonomy, but rather reveals haplotype clades that are incongruent with nominal species boundaries: instead of seven species, we recognized only two major, yet incompletely segregated, lineages. Our results are consistent with the hypothesis that Maniola originated in Africa. We suggest that one lineage dispersed over the Strait of Gibraltar and the Iberian Peninsula to the west of Europe, while the other lineage spreads eastward through Asia Minor and over the Bosporus to Eastern Europe.


Introduction
Taxonomic mishaps such as "oversplitting" (the misinterpretation of individual variants as distinct specific entities) and "lumping" (erroneously grouping several species into a single one) continue to pose problems in modern taxonomy and systematics (Dayrat 2005). In the last decade, DNA barcoding has become a handy tool for resolving species identifications, especially in clades where morphological characterization of putative taxa is weak or inconsistent (Dupuis et al. 2012;Puillandre et al. 2012;Ratnasingham and Hebert 2013).
As a consequence, supposedly well-established taxonomic systems for many groups of organisms have become subject to frequent, and often drastic, change due to ongoing revisions. Newly described cryptic species and novel insights into relationships between taxa frequently overturn traditional systematics, even in comparatively well-known groups such as butterflies (e.g., Dinc a et al. 2011 ;Talavera et al. 2013). In particular, the use of DNA barcoding approaches has led to a very substantial increase of recognized species numbers during the last decade (Hebert et al. 2004;Hausmann et al. 2011).
Meadow Brown butterflies have been in the focus of studies in evolutionary genetics and ecology since decades (Scali 1971;Brakefield 1982a,b;Goulson 1993). The center of these studies has been the widespread M. jurtina. The other species in the genus, especially the island endemics, have been largely neglected since their description, apart from a few anecdotal taxonomic papers (e.g., Jutzeler et al. 1997), some work on ecological aspects of differentiation (Grill et al. 2006a,b) and work on morphometrics in M. jurtina and M. nurag, particularly in genitalia (Dapporto 2010). Maniola jurtina alone has been studied more extensively, but so far, no conclusive picture of phylogeographic patterns has emerged (Schmitt et al. 2005;Dapporto et al. 2011). As taxonomic uncertainties last due to great morphological variation and overlap in distribution areas, we here investigated species boundaries with molecular markers.
Maniola comprises seven described species with a rather peculiar distribution: Maniola jurtina (Linnaeus, 1758) is widely distributed over much of Europe, northwest Africa, the Canary Islands, and eastward to the Ural mountains and NW Iran (Tshikolovets 2011), whereas the other six species in the genus are restricted to the Mediterranean area (Grill et al. 2006a(Grill et al. , 2007. Among them, four species are narrowly endemic to islands and therefore belong to those butterflies of Europe with the smallest ranges: M. cypricola Graves, 1928 on Cyprus, M. chia Thomson, 1987 on the Aegean islands Chios and Inousses, M. halicarnassus Thomson, 1990 on the Aegean island Nisyros and the neighboring Turkish Bodrum peninsula, and M. nurag Ghiliani, 1852 on Sardinia (Grill et al. 2007; Kudrna et al. 2011). The remaining two species inhabit largely overlapping parts of the eastern Mediterranean realm: M. telmessia (Zeller, 1847) is found on many Aegean islands and throughout southern Turkey from the Bosporus eastward to Syria, Lebanon, Israel, Jordan, Iraq, and NW Iran (Hesselbarth et al. 1995;Tshikolovets 2011). The range of M. megala (Oberth€ ur, 1909) is smaller, from the Aegean island Lesbos through southern Turkey as far as the Syrian border (Tshikolovets 2011).
With these unusual distributions, Maniola stands as an example for a genus with one widespread and several regional species, a pattern that is found similarly in a range of Palaearctic animal genera, not only in butterflies but in fact in many groups including vertebrates (e.g., in wall lizards Podarcis: Poulakakis et al. 2005; in birds of the genus Oenanthe: Randler et al. 2012;in butterflies: Dennis et al., 2008 ;Dapporto and Strumia 2008).
Maniola jurtina, M. telmessia, M. halicarnassus, and M. megala are broadly sympatric in Turkey, yet evidence for differential habitat preferences is vague (Hesselbarth et al. 1995). Despite minor variation in emergence times, hybridization has been recorded, for example between M. halicarnassus and M. telmessia (Hesselbarth et al. 1995) as well as M. jurtina and M. nurag on Sardinia (Grill et al. 2006b). All Maniola species strikingly resemble each other in terms of wing patterns and genitalia morphology, making taxonomy a challenging task even for experts (Thomson 1973;Olivier 1993;Grill et al. 2004). In Table 1, a survey of described morphological and ecological differences of the Maniola species is given.
Available studies on phylogeography within Maniola conducted so far never included all species of the genus and never had a geographically large enough taxon sampling (Schmitt et al. 2005;Grill et al. 2006a;Dapporto et al. 2011). So, this study aims at (1) reconstructing the phylogeny of the genus Maniola, (2) comparing it with the currently accepted taxonomy, (3) testing the usability of DNA barcoding for species identifications in Maniola, and (4) investigating whether molecular data reveal information about the existence of refugia in the Mediterranean region and possible expansion routes that have led to the current distribution of the species. What we find for Maniola may indicate that similar outcomes are to be expected for a number of other genera with similarly peculiar distribution patterns.

Material and Methods
Sample collection, DNA extraction, amplification, sequencing, and alignment We used 138 Maniola individuals, from sites across Europe, Anatolia, and northern Africa (Table S2) representing all seven taxonomically described species. Specimens had either been dried or put into 99% ethanol after collection. Samples were then stored at À20°C until DNA extraction. The applied mitochondrial and nuclear genetic markers (cytochrome-c-oxidase [COI], cytochrome-B [CytB], wingless [wgl], and elongation factor 1a [EF1a]) were amplified with varying success (Table 2), as the specimens had been collected between the years 1980 and 2013; older samples performed often worse in PCR than the more recent ones. We additionally included COI sequences from 51 individuals from NCBI GenBank.
As outgroups, we downloaded sequences of closely related satyrine butterfly species from GenBank (see Table  S1), from the subtribes Coenonymphina, Erebiina, Maniolina, Melanargiina, Pronophilina, and Satyrina, to improve support of the phylogeny of the genus. Outgroups were selected from Peña et al. (2006), but as they did not use CytB in their study, we additionally sequenced two specimens of Pyronia cecilia that were at hand in our laboratory to have at least one species available with complete gene sampling.
Total genomic DNA was extracted from two legs, respectively, or, if missing, from thoracic muscles following a standard protocol (DNeasy Blood & Tissue Kit, Qiagen Inc., Valencia, CA). DNA samples were extracted and amplified in a separate room exclusively dedicated to DNA extractions. Primer names, references, primer sequences as well as respective annealing temperature and time are shown in Table 3.
Gene fragments of the mitochondrial (COI, CytB) and nuclear DNA (wgl, EF1a) were amplified using polymerase chain reaction (PCR) in a thermal cycler (Eppendorf Mastercycler pro S vapo.protect). Negative (sterile water)  Hesselbarth et al. (1995), Olivier (1993), Tshikolovets (2011), Tolman and Lewington (2008) and Makris and John (2003 and positive (samples with known genotypes) controls were always used. PCRs were performed in 25 lL volumes containing 1 lL of genomic DNA, 22.5 lL of Red-dyMix â , 0.5 lL of the respective forward and reverse primer, and 1 lL BSA. PCR conditions were optimized for each primer pair. PCR products were visualized on an agarose gel to verify amplification success. Afterward, they were analyzed using an ABI capillary sequencer (3730 DNA analyzer, Applied Biosystems). Sequences contained no gaps or stop codons and were aligned and edited in BioEdit 7.1.3.0 (Hall 1999) for each gene separately. All new sequences have been submitted to GenBank, under accession numbers KP032366 -KP032458 (Cytb), KP032241 -KP032365 (COI), KP032552 -KP032639 (EF1a), and KP032459 -KP032551 (wgl). Accession numbers of outgroup species are provided in Table S1.

Phylogenetic analysis
We used Bayesian inference (BI) to reconstruct phylogenetic trees. The combined genes tree (including 65 individuals, for which all genetic markers worked) was based on an alignment of 2427 base pairs. Trees of COI included 118 sequences with 657 bp length. We used the program JModelTest (v. 2.0) on the Phylemon-Server 2.0 (S anchez et al. 2011) to estimate models of nucleotide substitutions (of three substitution schemes: JC, HKY, and GTR) as judged by the Bayesian information criterion (BIC) for BI. Only the best-fit models were subsequently used for evaluations of tree topology. BI analysis was performed with MrBayes v.3.2, considering the estimated best-fit substitution models (Huelsenbeck and Ronquist 2001). Trees were visualized using Figtree v.1.3.1 (Rambaut 2009).
For BI trees under BIC, the best-fit models were as follows: [GTR + G] for COI, [HKY + I + G] for CytB, [GTR + I + G] for mtDNA and all genes combined, [K80 + I] for Elongation factor and nDNA, and [HKY + I] for wingless.
Genetic variation within described species was estimated as the numbers of variable sites (S), average numbers of nucleotide differences (k), haplotype diversity (h: Nei, M. 1987), and nucleotide diversity (p: Nei and Li 1979) for each gene as well as species, using the software DnaSP v5 (Librado and Rozas 2009). Genetic distances (Kimura-2-distances) were calculated with MEGA 5.1 (Tamura et al. 2011).

Haplotype networks
Additionally, we generated median-joining networks for several datasets using the program Network 4.6.1.1 (http://www.fluxus-engineering.com). The median-joining method is based on a maximum parsimony algorithm that searches for all shortest trees of a particular dataset (Bandelt et al. 1999).  For all species and all genetic markers (Table 4), estimates of haplotype diversity (h) were rather high, ranging from 0.333 to 0.996, whereas estimates of nucleotide diversity (p) were much lower (especially in nuclear genes), ranging from 0.00112 to 0.02804. Interestingly, the endemic species did not show lower haplotype and nucleotide diversity than their mainland congeners, and M. nurag even showed the highest diversity of all as well as very high average numbers of nucleotide differences (k). The available M. megala specimens all had exactly the same haplotype (see Fig. 4), and M. halicarnassus revealed rather low numbers of variable sites. As two lineages could be detected in the phylogenetic trees (see subsection Phylogeny and phylogeography), they were also investigated for their nucleotide diversity: Lineage A (comprising the taxa M. cypricola, M. telmessia, and M. halicarnassus) generally had lower nucleotide diversity than Lineage B (all other nominal species).
To use the barcoding region for species identification, a distinct barcoding gap should existseparating intrafrom interspecific pairwise genetic distances. Among the 118 Maniola COI sequences, no such gap could be found, as intra-and interspecific distances for conventionally delineated species were intermixing. When only the two major genetic lineages (see below) were compared, genetic distances between and within lineages were separated more clearly from each other, although a convincing barcoding gap could still not be found (Fig. 1). Approximately 99.8% of pairwise COI sequence comparisons showed two or more percent genetic distance between these two lineages (interlineage), but only 51% of pairwise comparisons showed negligible (0-1% genetic distance) distances within lineages (intralineage). About 26% of sequence comparisons showed high (≥2%) distances within lineages.

Phylogeny and phylogeography
The topology of the Bayesian Inference tree of all Maniola samples, based on specimens for whom all genetic markers were available plus the selected outgroup species (Fig. 2), is not consistent with current taxonomy. Rather, most nominal species form mixed clades. Nevertheless, recurrent clusters can be recognized: there is one lineage of M. jurtina (from various provenances; mixed with M. chia, M. nurag, and M. megala) (Lineage B in the lower part of the tree diagram); another one containing a few M. jurtina (from Crete), M. nurag, two M. cypricola specimens, and one M. telmessia sequence (upper Lineage B); and one M. telmessia (together with all M. halicarnassus and the majority of M. cypricola sequences) lineage (Lineage A). In the BI tree of COI sequences (Fig. 3), Lineage B forms a single clade, so the two parts of this lineage are both called Lineage B in the BI tree of the combined genes, as Lineage B can be defined as all specimens that do not belong to Lineage A. The "mixed island species clade" within Lineage B (marked in blue in Fig. 2 The tree based exclusively on the COI barcoding sequences shows similar groupings (Fig. 3), but more clearly emphasizes the split of haplotypes into two main lineages: One lineage (Lineage A; light gray) predominantly contains Maniola telmessia, M. halicarnassus, and M. cypricola, whereas the other contains M. jurtina, M. nurag, M. chia, and M. megala (Lineage B; dark gray), although both with some exceptions (black arrows). Lineage A shows a probability value of 76%.
Median-joining networks were calculated for COI, mtDNA, nDNA, and combined genes and resembled the phylogenetic trees largely; but in comparison, they show a clearer resolution of the relationships between haplotypes. Networks of COI, mtDNA, and combined genes were much alike: they showed the same underlying clustering, but the more genes were used the more single haplotypes and mutational steps between clusters were found. Only the COI network is shown (Fig. 4), as it is the most clearly arranged. Groupings are similar as in the phylogenetic trees described above.  Table 4. Summary of molecular diversity indices of COI, CytB, wgl, elongation factor 1a genes, mitochondrial DNA (combined COI and CytB), and nuclear DNA (nDNA); sample size (n), number of haplotypes (no.), number of variable sites (S), average number of nucleotide differences (k), haplotype diversity (h), and nucleotide diversity (p) with standard deviation (SD). Nuclear genes have been doubled for analysis by the program -DNAsp to avoid ambiguous sites, as DNAsp cannot analyze diploide genetic information.

Discussion
Our data indeed show one of theat least up to nowrare cases (cf. Vila et al. 2010), where the differentiation of several nominally described, but rather ambiguously characterized species, did not receive stronger support and higher resolution of "cryptic diversity" by use of genetic methods. Although one might have suspected the existence of even more cryptic species among Maniola because of their geographic distribution around the Mediterranean Sea and its many islands, the opposite was found, viz. according to our data, the whole genus consists of only one, quite variable species.

DNA barcoding, phylogenies, and taxonomic implications
Haplotype as well as nucleotide diversity of the whole Maniola genus is comparable with those of a single widely distributed species (e.g., Jeratthitikul et al. 2013). The high genetic diversity of the island species within Maniola indicates gene flow among the so-called species. Island populations (Frankham 1997), and isolated populations in general (Cassel and Tammaru 2003), typically have less genetic variation than continental populations due to bottlenecks experienced by small founder populations and subsequent inbreeding (Frankham 1998;Keller and Waller 2002). These outcomes are in accordance with earlier results by Grill et al. (2007), who also found no evidence for lower genetic variation in Sardinian populations of Meadow Brown butterflies. DNA barcoding based on COI sequences cannot be reliably used to split Maniola further into a number of species, in contrast to many cases reported in the recent Figure 1. Frequency distribution of pairwise intra-and interlineage Kimura-2-distances of the COI sequences. To ensure the usability of DNA barcoding for species delimitation, a "barcoding gap" should exist between these two data series. In Maniola, however, genetic distances within lineages versus between lineages intermix. barcoding literature (Hausmann et al. 2011;Ratnasingham and Hebert 2013). The tree constructed from the barcoding region of the COI gene shows a clear split into only two lineages, which both range across boundaries of all species taxonomically described for this genus. Addition of further sequence markers did not enhance resolution of clades. All individuals (with the exception of  M. megala) invariably clustered according to a different pattern than expected from their current taxonomic affiliation.
The results of this study suggest that, instead of accepting seven distinct "good" species, it is more parsimonious to assign all Meadow Browns to one single genetic group. One may distinguish the two lineages A and B as different genetic entities, with some weak further substructure according to haplotypes.
It is understandable that a taxon morphologically as variable as Maniola tempted taxonomists of the 19th and 20th century to describe new species; but according to our genetic data, these descriptions delineate variation, not speciation. For example, the species M. chia and M. halicarnassus (which were proposed as distinct taxa as late as 1987 and 1990, respectively) probably should have been investigated more critically before describing them as new species. If there is endemicity on an island, it is usually observed in a number of unrelated groups of organisms (cf. the high number of endemic species on the Tyrrhenian islands or the Azores). On Chios, as a matter of fact, besides M. chia, only the terrestrial isopod Trachelipus buddelundi is supposed to be endemic ("Only known from its original description. A doubtful species."; Alexiou and Sfenthourakis 2013). Within the genus Maniola, wing patterns or genitalia morphology are not only determined through the genotype, but also by environmental conditions influencing larval development (Thomson 1973;Goulson 1993;Grill et al. 2004). Hence, even if significant and substantial, differences in phenotypes alone are insufficient to delineate taxonomic entities in that group of butterflies.
Interestingly, a dichotomylike the two lineages we foundwithin the "super species" (or species complex) M. jurtina has been postulated multiple times. Starting from Tauber (1970), also Hesselbarth et al. (1995) recognized groupings that were assigned to the more western M. jurtina sensu stricto on the one hand and the eastern M. telmessia on the other hand. Later on, Schmitt et al. (2005) reported the existence of an eastern and a western genetic lineage of Maniola jurtina based on allozyme data, wing patterns, and genitalia morphology as did Dapporto et al. (2011) based on genitalia morphology. Dapporto et al. (2011) rightly called M. jurtina "enigmatic", as they found contradictory patterns between allozyme and morphological data. Specimens from continental Italy, Sicily, and North Africa shared the same allozyme set, but their genitalia shape was more similar to specimens from the Balkans and some individuals from Eastern and Central Europe. Their explanation for this discrepancy is "recent gene flow in the wake of postglacial range expansions and shifts" (Dapporto et al. 2011). We hypothesize that these contrasting patterns could have come about through by chance-biased sampling from the genetic lineage A for the allozyme study and from lineage B for the genitalia analyses. At any rate, genetic lineages connected to the main Pleistocene glacial refugia, as postulated by Hewitt (1996Hewitt ( , 1999, for example, structured according to the three large Mediterranean peninsulas, could not be detected within Maniola. Hypothetical evolutionary scenario of the genus Maniola and possible expansion routes concerning the origin and expansion routes of the whole genus was much earlier suggested by Tauber (1970), who postulated that the western M. jurtina lineage and the eastern M. telmessia species complex diverged from North African ancestors. According to Tauber, the M. jurtina lineage spreads westward over Gibraltar and the Iberian Peninsula, while the M. telmessia complex expanded through Palestine, Lebanon, and Syria to the east. This hypothesis coincides with the fact that M. telmessia is the only Maniola "species" occurring in Israel, Jordan, Lebanon, and Syria (Tshikolovets 2011). According to Tauber's hypothesis, the western M. jurtina complex spreads across Europe and would have met the eastern M. telmessia complex in Asia Minor, giving rise to a sympatric occurrence of both lineages in the eastern Mediterranean nowadays.
Although over 40 years old, Tauber's (1970) "out-of-Africa" hypothesis seems to be a likely scenario for the evolution of the genetic diversity we find in Maniola today. A number of other recent studies revealed similar scenarios for other butterfly species Kodandaramaiah and Wahlberg 2009;Habel et al. 2010;Husemann et al. 2014). Another fact fitting to an African origin of Maniola is their current distribution: while closely related Satyrinae taxa such as Hyponephele sp. and Aphantopus sp. show trans-Eurasian distributions, the range of Maniola jurtina ends roughly at the Ural Mountains (Tshikolovets 2011).
If the geographic origin of Maniola lay in northern Africa (as it is for other nymphalid butterflies: Aduse-Poku et al. 2009;Kodandaramaiah and Wahlberg 2007), the butterflies could have spread northwestward over Gibraltar and the Iberian Peninsula and then eastward all over Europe and into western Asia. A convincing scenario is exactly what Tauber (1970) hypothesized based on morphological characters and paleo-ecological considerations: a split in the North African stem population, with one part dispersing to the west and the other one to the east. Considering our data, Lineage A (M. cypricola, M. telmessia and M. halicarnassus) could have expanded along an eastern pathway, and Lineage B (M. jurtina, M. nurag, M. chia, and M. megala) could have colonized western and Central Europe through a western migration route (Fig. 5). The few individuals that cluster with the "wrong" lineage might indicate a zone of rather recent intermixture in the Aegean Sea and Central Europe. Of course, to further investigate this hypothesis, more haplotypes from this hybrid zone as well as from northern Africa must be examined.

Conclusions
During the last decade, phylogenetic studies using sequence data revealed numerous examples for the unexpected discovery of many cryptic species. Our study presents the not so common opposite case: various distinct nominal species melting together to a single species complex. We did not find any proof for the existence of the seven morphologically defined species with the genetic methods used; instead, we found two but moderately distinct genetic lineages. As to the formation of the two lineages, we hypothesize an origin in Africa and two different expansion routes emerging from there, as postulated by Tauber (1970). Thus, we suggest addressing the whole genus as one "super species" M. jurtina. Numerous studies in the past few years have strikingly uncovered similar "out-of-Africa" examples, lending further support to this idea for the satyrine genus Maniola. Future studies, with a larger set of samples from the African continent, will be needed to evaluate our hypothesis. Most importantly, our results raise the controversial question whether oversplitting of species, despite all the contrasting evidence in recent bar code studies, might be more common than expected until now.

Supporting Information
Additional Supporting Information may be found in the online version of this article: Table S1. Sequences from GenBank with accession numbers used in this study. All sequences from Peña et al. (2006), except CytB sequences. Table S2. List of specimens used in this study and their collection localities and dates, when available.