The haploid chromosome number n=28 was found in prometaphase I, MI and MII cells of seven studied individuals (Table 1, Fig. 1). All chromosome elements formed a gradient size row. The karyotype contained no exceptionally large or small chromosomes.

Figure 1.

Karyotypes in male meiosis of Melitaea gina from Iran. a sample Q183, prometaphase I, n = 28 b sample Q153, late prometaphase I, n = 28 c sample Q183, MI, n = 28 d sample Q155, M I, n = 28. Scale bar corresponds to 10µ in all figures.

Bayesian analysis of the barcode region recovered the M. didyma complex as a monophyletic clade (Fig. 2), which agrees with Leneveu et al. (2009). Despite low karyotype variability, the clade was found to have unexpectedly high mitochondrial haplotype diversity. These haplotypes were clustered in 23 highly diverged haplogroups called chitralensis, deserticola, didyma, didymoides, gina, gina2, interrupta, latonigena, liliputana, mauretanica, mixta, neera, neera2, occidentalis, protaeoccidentis, saxatilis, sutschana, sutschana2, sutschana3, turkestanica, turkestanica2, turkestanica3 and turkestanica4 (Figs 2–6, Suppl. material 1). These haplogroups had high support (Bayesian posterior probability from 0.95 to 1) and were associated with particular geographical areas (Fig. 7).

Figure 2.

The Bayesian tree of Melitaea based on analysis of the cytochrome oxidase subunit I (COI) gene. Numbers at nodes indicate Bayesian posterior probability.

Figure 3.

Fragment of the Bayesian tree of Melitaea didyma complex (haplogroups neera and liliputana) based on analysis of COI gene. Numbers at nodes indicate Bayesian posterior probability.

Figure 4.

Fragment of the Bayesian tree of Melitaea didyma complex (haplogroups interrupta, occidentalis, saxatilis, lathonigena, didymoides, sutschana, sutschana 2, sutschana 3) based on analysis of COI gene. Numbers at nodes indicate Bayesian posterior probability.

Figure 5.

Fragment of the Bayesian tree of Melitaea didyma complex (haplogroups turkestanica, turkestanica 2, didyma) based on analysis of COI gene. Numbers at nodes indicate Bayesian posterior probability.

Figure 6.

Fragment of the Bayesian tree of Melitaea didyma complex (haplogroups mixta, chitralensis, mauretanica, protaeoccidentis, neera2, gina and deserticola) based on analysis of the COI gene. Numbers at nodes indicate Bayesian posterior probability.

Figure 7.

Distribution ranges of western COI haplogroups of Melitaea didyma complex.

The uncorrected mean p-distances between the haplogroups were high (up to 9.1% between turkestanica4 and deserticola) (Table 2). The majority of them were much higher than the ‘standard’ 2.7–3.0% DNA barcode threshold usually used for allopatric taxa as an indicator for their species distinctness (Lambert et al. 2005, Lukhtanov et al. 2015).

Most of the haplogroups were found to be allopatric. However, in some cases barcodes’ clusters did not correspond to the simple allopatric geographical distribution. The sample Melitaea gina M22 (haplogroup gina2) was found in sympatry with the haplogroup gina in north-west Iran. The distance between gina and gina2 was 6.5%. Haplogroups turkestanica4, turkestanica3 and turkestanica2 were highly diverged (up to 7.4%) as compared with the haplogroup turkestanica and were found in sympatry with the haplogroup turkestanica (Fig. 8). In Slovenia, the specimen BPAL3090-15 (haplogroup neera2) was found together with the haplogroup neera. The distance between neera and neera2 was 1.7%.

Figure 8.

Localization of neera and turkestanica haplogroups (yellow circles – neera, black – turkestanica, green – turkestanica2, red – turkestanica3, blue – turkestanica4)

Two samples with the turkestanica haplotypes (haplogroup turkestanica), one from Aktobe (Kazakhstan) and one from Samara (Russia) were found in sympatry with M. dimyma neera haplotypes (haplogroup neera). In Karabiryuk (Kazakhstan), two samples with the neera haplotypes (haplogroup neera) were found in sympatry with M. didyma turkestanica haplotypes (haplogroup turkestanica and turkestanica4).

The genus Melitaea is known to be characterized by relatively low interspecific chromosome number variation. The representatives of basal clades (see phylogeny in Leneveu et al. 2009), the taxa of M. cinxia (Linnaeus, 1758), M. diamina (Lang, 1989), M. athalia (Rottemburg, 1775), M. trivia ([Denis & Schiffermüller], 1775) and M. phoebe ([Denis & Schiffermüller], 1775) species groups demonstrate n=30-31 (Federley 1938, de Lesse 1960, Robinson 1971, Larsen 1975, Hesselbarth et al. 1995). These haploid numbers are modal ones not only for Melitaea, but also for the family Nymphalidae and for the order Lepidoptera in whole (Robinson 1971, Lukhtanov 2000, 2014). Most likely, one of them (probably, n=31, see Lukhtanov 2014) represents an ancestral lepidopteran condition preserved in the basal lineages of Melitaea.

The younger lineages, the M. fergana Staudinger, 1882 and M. didyma species groups, were found to possess lower chromosome numbers varying from n=27 to n=29-30. Within the M. fergana species group, M. athene Staudinger, 1881, the only karyologically studied species, was found to have n=29 (with n=30 as a rare intra-individual variation) (Lukhtanov and Kuznetsova 1989). The species-rich M. didyma group consists of three complexes: a complex of taxa close to M. ala, a complex of taxa close to M. persea and a complex of taxa close to M. didyma. Within these complexes the following chromosome numbers were found: n=29 in M. ala (Lukhtanov and Kuznetsova 1989), n=27 in M. persea (de Lesse 1960) and different numbers from n=27 to n=29-30 in species of the M. didyma complex (Table 3).

Together with M. deserticola (n=29, Larsen 1975), M. gina occupies a basal position within the M. didyma complex (Fig. 6). Therefore analysis of M. gina was crucially important for understanding chromosome number evolution in this complex. Our study revealed M. gina to have n=28, a number previously observed in M. didyma from Italy (de Lesse 1960) and M. didyma neera from the Kazakh Altai (Lukhtanov and Kuznetsova 1989). Taking into account absence or relatively low level of interspecific chromosome number variation in the M. didyma complex and presence of intraspecific variation (Table 3), we conclude that in this group chromosome numbers have relatively low value as taxonomic markers (but see: Lukhtanov and Kuznetsova 1989).

Despite low level of chromosome number variability, the M. didyma complex was found to have unexpectedly high level of mitochondrial haplotype diversity. These haplotypes were clustered in 23 highly diverged haplogroups (Fig. 2). 12 of these haplogroups are associated with nine traditionally recognized and morphologically distinct species M. deserticola, M. gina, M. didymoides, M. saxatilis, M. sutschana (this species was devided recently in M. sutschana and M. yagakuana Matsumura, 1927, see Oorschot and Coutsis 2014), M. latonigena (this species was devided recently in M. latonigena and M. latonigenides Oorschot and Coutsis, 2014, see Oorschot and Coutsis 2014), M. interrupta, M. mixta and M. chitralensis.

The rest of the haplogroups belong to the well-known west-palearctic species M. didyma. Despite intrapopulation and seasonal variability, this species is very homogenous with respect to morphology, including the structure of genitalia, a character which is most useful for species separation in Melitaea (Suschkin 1913, Higgins 1941, Oorschot and Coutsis 2014). In accordance with this homogeinity, in the recent revision (Oorschot and Coutsis 2014) all populations of this species, except for Central Asian populations, were considered as members of the same subspecies M. didyma didyma. The populations from Central Asia were treated by Oorschot and Coutsis (2014) as a separate subspecies M. didyma turkestanica.

If we follow the opinion of experts in Melitaea taxonomy (Kolesnichenko et al. 2011, Oorschot and Coutsis 2014) and accept the traditional taxonomic treatment of the species M. didyma, we should acknowledge that this species is particularly unusual in the haplotypes we obtained. First, it is clearly polyphyletic with respect to COI gene, and the lineages of M. didyma are intermixed with other well recognized species on the tree (Figs 2–6). Second, the number of distinct COI lineages within M. didyma is unusually high (11 lineages) and their genetic differentiation is extreme. The majority of these haplogroups are allopatric, but some of them have sympatric (neera/neera1, turkestanica/turkestanica2, turkestanica/turkestanica3, turkestanica/turkestanica4) or partially sympatric (neera/turkestanica, occidentalis/didyma) distribution. The mean uncorrected pairwise distances between the lineages is up to 7.4% if the lineages turkestanica3 and turkestanica4 are considered (Table 2). The lineages turkestanica3 and turkestanica4 are the most diverged lineages of M. didyma. Together with gina2, on the tree (Fig. 2) they have an intermediate position between the lineage (M. didyma + M. deserticola + M. gina) and the lineage (M. persea + M. casta). It even appears as a sister group to (M. persea + M. casta), but with a very low support (0.54). However, even if the lineages gina2, turkestanica3 and turkestanica4 are not considered, the distances between M. didyma groupings remains high, up to 4.1% between turkestanica2 and liliputana, i.e. much deeper than the “standard” DNA barcode species threshold (2.7-3%) (Hebert et al. 2003, Lukhtanov et al. 2016).

There are two theoretically possible explanations for this pattern. First, M. didyma sensu auctorum can be a mix of multiple species that mostly have allopatric distribution ranges, but some of them are sympatric. Second, the recovered haplogroups (at least the allopatric ones) can represent highly diverged intraspecific lineages. Of course, a combination of the first and the second hypotheses is possible, and a part of the haplogroups could represent different species, and another part of the haplogroups could represent intraspecific variations.

In our opinion, the second hypothesis seems to be more plausible. There are the following arguments for the second scenario. First, no morphological differences between the bearers of these haplogroups are known (except for lighter, more yellowish wing colour in the three M. didyma turkestanica haplogroups as compared with other haplogroups). The second (and the most convincing) argument is based on our field obseravtion of butterfly habitats and ecological preferences. In ecology the competitive exclusion principle, also known as Gause’s law is one of the most important rule (Gause 1934, Hardin 1960). In complete accordance with this rule, in case of sympatry the most closely related species pairs, such as M. didyma/M. interrupta, M. didyma/M. latonigena and M. gina/M. saxatilis demonstrate clear niche differentiation (M. didyma and M. gina are more xerophilous, whereas M. interrupta, M. latonigena and M. saxatilis are more mesophilous taxa). This was not a case for sympatric haplogroups neera/neera2, turkestanica/turkestanica, turkestanica/turkestanica3 and turkestanica/turkestanica4 (Fig. 8). The bearers of these haplogroups were not only morphologically identical, but also were found to fly exactly syntopically and synchronously. This pattern is hardly compatible with non-conspecifity of these haplogroups.

M. didyma neera and M. didyma turkestanica are differentiated ecologically (Pazhenkova et al. 2015), however, there was no ecological separation between bearers of the neera and turkestanica haplogroups in cases of their sympatry. In Samara and Aktobe, where the haplogroup neera was predominant, both haplogroups were found in M. didyma neera biotope (steppe), and in Karabiryuk where the haplogroup turkestanica was predominant, both haplogroups were found in M. didyma turkestanica biotope (desert) (Fig. 8). This pattern corresponds more to a result of haplotype introgression than to co-habitation of two ecologically differentiated species.

Interestingly, the haplogroup turkestanica2 is not related to the haplogroup turkestanica and is a derivative from West-European haplogroup didyma. This pattern can be treated as a result of ancient introgression. Generally, footprints of ancient and more recent introgression are both an evidence for transparency of boundaries between M. didyma populations.

The mega-analysis of species-level para- and polyphyly in DNA barcode gene trees was recently conducted by using a huge data set (4977 species and 41,583 specimens of European Lepidoptera) (Mutanen et al. 2016), however without in-depth-analyses of particular cases. This study resulted in conclusion that cases of species’ polyphyly in COI trees arising as a result of deep intraspecific divergence were negligible, and the detected cases reflected misidentifications or/and methodological errors. Despite this, our analysis demonstrates that species-level polyphyly in DNA barcode based on deep intraspecific divergence may be a real phenomenon.

The M. didyma complex consists of at least 23 COI haplogroups, the majority of which demonstrated a strict attachment to particular geographic ranges: chitralensis (north Pakistan); deserticola (north Africa, Israel, Jordan, Lebanon, Syria); didyma (west Europe); didymoides (Asian Russia, Mongolia, North China); gina (W Iran, Azerbaijan); interrupta (Caucasus, NE Turkey); latonigena (Asian Russia, north-east Kazakhstan, Mongolia, north-west China); liliputana (Armenia, Turkey, Syria, Lebanon, Israel); mauretanica (south Spain); mixta (Tajikistan, Kyrgyzstan, Uzbekistan, Pakistan, Afghanistan); neera (east Europe, north Caucasus, west Siberia, north Kazakhstan); occidentalis (Spain); protaeoccidentis (north Africa); saxatilis (north Iran); sutschana (Russian Far East, Korea, north-east China) and turkestanica (Kazakhstan, Kyrgyzstan, Uzbekistan, Tajikistan, west China). With few exceptions (e.g. deserticola/protaeoccidentis, deserticola/liliputana), the ranges of these haplogroups do not overlap substantially (Fig. 7), and we hypothesize that mitochondrial diversity was formed in allopatry. Given the deep level of genetic differentiation between the lineages, we assume that there was a long period of allopatric differentiation when the lineages were separated by geographic or/and ecological barriers. Under generally accepted maximum 2.3% (Brower 1994) and minimum 1.5% uncorrected pairwise distance per million years (Quek et al. 2004) for COI sequence of various arthropod taxa, this period can be estimated to be as long as 0.5–5.0 My. In our opinion, this is an evidence that each of these haplogroups evolved in one of the main west-palaearctic late Pliocene and Pleistocene refugia in north Africa (protaeoccidentis, deserticola), the Iberian Peninsula (occidentalis, mauretanica), the Balkan Peninsula (neera), the Middle East (liliputana, saxatilis, gina) and Central Asia (turkestanica, mixta, chitralensis). The presence of additional diverged minor haplogroups neera2, turkestanica2, turkestanica3, turkestanica4, gina2, which could originate allopatrically in small isolated spots, but currently exist in secondary sympatry with major haplogroups neera, turkestanica and gina, agrees well with the refugia-within-refugia concept (Gòmez and Lunt 2007, Karaiskou et al. 2014). Interestingly, the area of the most diverged haplogroup turkestanica3 is close to the area of the recently described subspecies M. didyma carminea (Kolesnichenko et al. 2011).

We tentatively suggest interpreting the main clusters discovered within M. didyma sensu stricto (M. didyma didyma, M. didyma mauretanica, M. didyma occidentalis, M. didyma protaeoccidentis, M. didyma liliputana, M. didyma neera and M. didyma turkestanica) as subspecies because each of them has its own distribution range and is distinct with respect to mtDNA (i.e. represents by a monophyletic lineage or a combination of two or three monophyletic lineages). As a result we propose the following classification:

M. didyma (Esper, [1779])

M. didyma didyma (Esper, [1779])

M. didyma mauretanica Oberthür, 1909

M. didyma occidentalis Staudinger, 1961

M. didyma protaeoccidentis Verity, 1929

M. didyma liliputana Oberthür, 1909

M. didyma neera Fischer de Waldheim, 1840

M. didyma turkestanica Sheljuzhko, 1929

M. didymoides Eversmann, 1847

M. sutschana Staudinger, 1892

M. latonigena Eversmann, 1847

M. interrupta Colenati, 1846

M. mixta Evans, 1912

M. chitralensis Moore, 1901

M. deserticola Oberthür, 1909

M. saxatilis Christoph, 1873

M. gina Higgins, 1941

We provide the first record of M. gina in Azerbaijan (sample BPAL1697-12, Azerbaijan, Shamkir, 27 June 2011, collector V. Tikhonov).

We also record M. didyma turkestanica as a new taxon for Russia and Europe (samples BPAL3168-16, BPAL3169-16, BPAL3170-16, BPAL3173-16 Russia, Astrakhanskaya oblast, Bogdinsko-Baskunchaksky zapovednik, 24 May 2008, collector S. Nedoshivina).