Interpretation of mitochondrial diversity in terms of taxonomy: a case study of Hyponephele lycaon species complex in Israel (Lepidoptera, Nymphalidae, Satyrinae)

Abstract It is difficult to interpret mitochondrial diversity in terms of taxonomy even in cases in which a concordance exists between mitochondrial, ecological and morphological markers. Here we demonstrate this difficulty through a study of Israeli Hyponephele butterflies. We show that samples commonly identified as Hyponephele lycaon are represented on Mount Hermon in Israel by two sympatric groups of individuals distinct both in mitochondrial DNA-barcodes (uncorrected p-distance = 3.5%) and hindwing underside pattern. These two groups were collected in different biotopes. They also tended to be different in length of brachia in male genitalia, although the latter character is variable. We reject the hypothesis that the discovered COI haplogroups are selectively neutral intraspecific characters. We hypothesize that they represent: either (1) two different biological species, or (2) a consequence of a strong positive selection acting at intraspecific level and resulting in two intraspecific clusters adapted to low and to high elevations. If we accept the first hypothesis, then provisionally these two haplogroups can be attributed to transpalearctic Hyponephele lycaon sensu stricto and to Hyponephele lycaonoides, previously known from Iran and East Turkey.

Within the genus, Hyponephele lycaon (Rottenburg, [1775]) is the best known and the most common species broadly distributed in the temperate zone of the Palearctic from Portugal in the west to Far East Russia in the east (Samodurov et al. 2001, Eck weiler andBozano 2011). In south Palearctic it is replaced by closely related allopat ric taxa H. maroccana Blachier, 1908 (North Africa), H. galtscha (GrumGrshimailo, 1893) (Tajikistan) and H. sifanica (GrumGrshimailo, 1891) (China) (Eckweiler and Bozano 2011). One more species, H. lycaonoides D. Weiss, 1978 was described from Zagros Mountains in Iran. Hyponephele lycaonoides was shown to be sympatric with H. lycaon in Iran (Weiss 1978, Eckweiler and Bozano 2011, Tshikolovets et al. 2014. Hyponephele lycaonoides was also reported for Turkey (Koçak 1989, Eckweiler andBo zano 2011), but the reports for Turkey were questioned in the comprehensive analysis of Turkish butterfly fauna made by Hesselbarth et al. (1995). Male genitalia struc tures are commonly used for distinguishing between H. lycaon and H. lycaonoides, and specimens with long brachia are attributed to H. lycaon, whereas specimens with short brachia are attributed to H. lycaonoides (Weiss, 1978). However, male genitalia are variable in both H. lycaon and H. lycaonoides, and intermediate forms are reported to be common (Eckweiler and Bozano 2011). Moreover, Hesselbarth et al. (1995) con sidered these traits (the long and short brachia) as intraspecific variations, rather than speciesspecific characters. Unfortunately, until now nobody used molecular markers to test the nonconspecifity of H. lycaon and H. lycaonoides.
In our study we analysed mitochondrial DNA barcodes and morphological and ecological markers to show that butterflies commonly identified as Hyponephele lycaon are represented in Israel by two sympatric groups of individuals. We further discuss different possible evolutionary and taxonomic interpretations of the pattern discovered.

Materials and methods
In the course of our DNA barcode survey of Israeli butterflies (2012-2015) we found butterflies similar to H. lycaon on Mount Hermon in northern Israel. They were col lected in a small area situated between 33°17'12"N, 35°45'49"E, at 1440 m and 33°18'38"N, 35°47'07"E, at 2050 m. The distance between these extreme points of the collecting was 3460 m (measured using Google Earth map). Some of the butterflies were collected in the forest zone at 1450-1600 m above sea level, other were collected in the subalpine zone with predominance of xerophytous vegetation at 1800-2050 m above sea level (Table 1).
DNA barcodes, 658 bp fragments within mitochondrial gene, cytochrome oxidase subunit I (COI), were sequenced at the Canadian Centre for DNA Barcoding (CCDB, Biodiversity Institute of Ontario, University of Guelph) using standard highthrough put protocol described in deWaard et al. (2008). DNA was extracted from a single leg removed from each voucher specimen employing a standard DNA barcode glass fibre protocol (Ivanova et al. 2006). All polymerase chain reactions (PCR) and DNA sequencing were carried out following standard DNA barcoding procedures for Lepi doptera as described previously (Hajibabaei et al. 2005). Photographs of specimens used in the analysis are available in the Barcode of Life Data System (BOLD) at http:// www.barcodinglife.org/. All voucher specimens are deposited at the Zoological Insti tute of the Russian Academy of Sciences and could be identified by the correspond ing unique BOLD Process IDs, that were automatically generated by BOLD, and by GenBank accession numbers ( Table 1).
The procedure of phylogenetic inference was described previously (Vershinina and Lukhtanov 2010, Talavera et al. 2013, 2015a. Briefly, the sequences were aligned using BioEdit version 7.1.7 software (Hall 1999) and ed ited manually. Phylogenetic relationships were inferred using Bayesian Inference and  (Ronquist et al. 2012). A GTR substitution model with gamma distributed rate variation across sites and with proportion of invariable sites was specified before running the program as suggested by jModelTest (Posada 2008). Two runs of 10,000,000 generations with four chains (one cold and three heated) were performed. Chains were sampled every 10,000 generations, and burnin was deter mined based on inspection of log likelihood over time plots using TRACER, version 1.4 (available at http://beast.bio.ed.ac.uk/Tracer). For comparison we used additional COI barcodes of Hyponephele downloaded from GenBank (Lukhtanov et al. 2009, Dinca et al. 2011.
Butterfly photographs were taken with Nikon D810 digital camera equipped with a Nikon AFS Micro Nikkor 105 mm lens. Genitalia photographs were taken with Leica M205C binocular microscope equipped with Leica DFC495 digital camera, and processed using the Leica Application Suite, version 4.5.0 software.

Results
During a 20122015 survey of Israeli fauna, H. lycaonsimilar butterflies were found only on Mount Hermon in northern Israel. We never observed H. lycaonsimilar but terflies in other parts of Israel, although the distantly related H. lupinus (Costa, 1836) was found not only in the northern, but also in central Israel. Thus, our observations support the finding that the geographic range of H. lycaon species complex is restricted in Israel to the northernmost part of the country (Benyamini 2002).
Molecular analysis of H. lycaonsimilar samples (Table 1, Fig. 1) revealed two dis tinct mitochondrial haplogroups (I and II) that were strongly differentiated with re spect to the COI gene. These two haplogroups differed from one another by 23 fixed nucleotide substitutions in the studied 658 bp fragment of the mitochondrial COI gene. When looking at the level of primary polypeptide structure, these differences translate to two fixed amino acid substitutions in the studied fragment. The minimal uncorrected COI pdistance between these two haplogroups was found to be as high as 3.5 %. Hyponephele lupinus from Israel was found to be closely related to H. interposita and distant from all the taxa of the H. lycaon complex.
With a single exception (female sample BPAL269514|CCDB17968_C11, Fig.  1, Table 1), the representatives of these two COI haplogroups were collected in differ ent biotopes (Fig. 2). The butterflies of haplogroup I were found on grassy slopes in the forest zone (14501600 m above see level) (Fig. 2a). The butterflies of haplogroup II were found in steppe lands of the subalpine zone (1800-2050 m alt.), where xerophy tous thorny cushion vegetation formed by Onobrychis cornuta and Astragalus species (Fabaceae) was predominant (Fig. 2b).
Standard χ2test was used to distinguish between random vs. nonrandom distribu tion haplogroups I and II in the low (forest) and high (subalpine) zones. Empirical and expected frequencies of COI haplogroups I and II in low and high altitude belts were compared ( Table 2). The calculated χ 2 was larger than the tabular value (9.558 vs. 6.635,  df = 1, 0.01 level of significance). Therefore, we reject the H 0 hypothe sis and conclude that haplogroup I butterflies are significantly more frequent in the lower zone, whereas haplogroup II butterflies are significantly more frequent in the higher zone. The representatives of these two clusters were also different in the pattern on the hindwing underside (Figs 3 and 4). In haplogroup I this pattern had more contrast with clearly visible medial band, whereas in haplogroup II the hindwing underside was paler and had less contrast.  A standard χ2test was used to distinguish between random vs. nonrandom as sociation between haplogroups I and II and hindwing underside pattern (Table 3). The calculated χ 2 of 12.860 was larger than the tabular value (12.860 vs. 10.83, df = 1, 0.001 level of significance). Therefore, we reject the H 0 hypothesis and conclude that COI haplogroup I is significantly associated with contrast pattern of the hindwing underside, whereas COI haplogroup II is significantly associated with pale pattern of the hindwing underside.
The representatives of these two COI haplogroups also tended to be different in the length of the brachia in male genitalia (Fig. 5), although the latter character had high variability. Males of haplogroup I often had long brachia (Fig. 5a, b), whereas males of haplogroup II were mainly characterized by reduced brachia (Fig. 5c, d).

Evolutionary interpretation of the discovered pattern
The COI genetic distance between haplogroups I and II (3.5 %) is higher than the 'stand ard' 2.7-3.0% DNAbarcoding threshold commonly used as a tentative indicator for spe cies distinctness of the taxa compared (Lambert et al. 2005, Lukhtanov et al. 2015a. It is known that COI barcodes alone are not sufficient for making any taxonomic decisions, since trees inferred from single markers sometimes display relationships that reflect the evolutionary histories of individual genes rather than the species being studied (Nichols 2001). Mitochondrial introgression ) and Wolbachia infection (Rit ter et al. 2013) can lead to additional bias when inferring taxonomic conclusions based on mitochondrial genes. Typically, multiple molecular markers or a combination of mor phological and molecular markers are required for inferring taxonomic hypotheses. In our research such additional information is represented by ecological characteristics (altitude belts). Less attention was attributed to the wing pattern, since we were not sure that it was an independent character. As the wing pattern strongly correlated with the ecology (low versus high elevation), one could hypothesize that the morphological difference is a consequence of phenotypic plasticity, i.e. ability of the same genotype to result in different phenotypes in response to changes in the environment (Price et al. 2003).
Three alternative explanations can account for bimodal sympatric distribution of mitochondrial markers. First, the diverged COI sequences may be selectively neutral intraspecific characters. Both preservation of a variety of ancestral haplotypes and mito chondrial introgression due to complex phylogeographic history could be responsible for such a neutral polymorphism (Avise 2000). Second, bimodal sympatric distribution of mitochondrial markers may be a result of a strong positive habitatrelated selection work ing at intraspecific level and resulting in two COI clusters associated with different alti tude belts (Cheviron and Brumfield 2009). Third, bearers of two diverged haplogroups may represent two different biological species (Avise 2000).
In our case the first hypothesis (neutral polymorphism) can be easily rejected. It predicts that the COI haplogroups I and II should be stochastically (i.e. randomly) distributed within high and low altitude belts. This prediction is not supported by χ2 test that demonstrated significantly nonrandom distribution of the COI haplogroups.
The second hypothesis (strong intraspecific positive selection) offers a more exotic, but not improbable, explanation. As COI sequence can be translated into a subunit of cytochrome c oxidase, a functional protein in mitochondria involved in energy metabolism (Kirk and Freedland 2011), this gene should be under natural selection (Castoe et al. 2008). Different haplotypes at this locus (or other linked mitochondrial genes) may be favoured in different environments. This could trigger a rapid sweep to fixation of a novel haplotype. This may result in sympatric clusters that differ in mi tochondrial genes while exchanging alleles freely throughout the rest of the genome. Interestingly, such groups maintained by habitatrelated selection could be considered species according to the genotypic cluster species concept (Coyne and Orr 2004, p. 448-449). The positive habitatrelated selection of mitochondrial genome, despite its theoretical plausibility, has so far relatively low empirical support, although there are some data confirming mitochondrial evolution along temperature and altitude gradi ents (Cheviron andBrumfield 2009, Quintela et al. 2014).
The third hypothesis (two different species) seems to be a more likely explanation in the case of haplogroups I and II, especially if one takes into account the high level of genetic divergences between the haplogroups and concordance between molecu lar (Fig. 1), ecological (Fig. 2, Table 2) and morphological (Figs 3 and 4, Table 3) characters. More samples, especially from the intermediate elevation (1600-1800 m), and analysis of additional nuclear molecular markers across altitudinal transect will be required in future research to support or to reject the second (positive selection) and the third (two species) hypotheses and to reveal potential nuclear gene flow between haplogroups I and II.

Taxonomic interpretation of the discovered pattern
The presence of two sympatric, ecologically differentiated groups within H. lycaon complex in the Middle East is not a completely novel issue. A similar situation is known to exist in Iran and East Turkey (Weiss 1978, Eckweiler and Bozano 2011, Tshikolovets et al. 2014. It is accepted by Hyponephele genus experts Bozano 2011, Tshikolovets et al. 2014) that in Iran and Turkey these two groups repre sent two different species: H. lycaon and H. lycaonoides (but see the alternative opinion: Hesselbarth et al. 1995). Although we understand that this point of view requires an additional justification, we may accept it as a working hypothesis until further investi gations and taxonomic revisions justify or falsify it.
If the species status of the discovered haplogroups will be confirmed in further stud ies, we suggest that, following Weiss (1978), the name H. lycaon (Rottenburg, [1775]) can be used for the Israeli taxon characterized by the contrast pattern on the hindwing underside and the predominance of longer brachia in male genitalia. Correspondingly, the name H. lycaonoides D. Weiss, 1978 can be used for the Israeli taxon characterized by the less contrasted pattern of the hindwing underside and the predominance of re duced brachia in male genitalia. However, this nomenclatural decision should be con sidered as a tentative one. First, despite recent revisions of the genus (Samodurov et al. 1995, 1996, 1997, 1999a, b, 2000, 2001, Eckweiler and Bozano 2011, no one studied typespecimens of numerous taxa that were described as subspecies and variations of H. lycaon. We cannot exclude that the name lycaonoides is a synonym of one of the pre viously described taxa, e.g. of libanotica (Staudinger, 1901). Second, molecular markers have never been used for analysis of taxonomic structure of H. lycaon species complex in its whole distribution range. Therefore, we will not be surprised if the true genetic and taxonomic structure of this group will be revealed as much more complex than a simple combination of two sympatric clusters as discovered in Iran, Turkey and Israel.