Wing pattern morphology of three closely related Melitaea (Lepidoptera, Nymphalidae) species reveals highly inaccurate external morphology-based species identification

Wing morphology of the three closely related species of Melitaea – M. athalia (Rottemburg, 1775), M. aurelia (Nickerl, 1850) and M. britomartis Assmann, 1847 – co-occurring in the Balkans (SE Europe) was investigated in detail through visual inspection, morphometric analysis and multivariate statistical analysis. Results are compared to recent phylogenetic studies, searching for concordant patterns and discrepancies between the two approaches. The morphology of the genitalic structures is also compared with the results of the other two approaches. The main conclusions are as follows: (1) small albeit significant differences in wing morphology exist among the three species and (2) while the structure of male genitalia and phylogenetic position of the three species are concordant, they are (3) in discordance with the wing morphology. The present study represents another example where identification based on external morphology would lead to highly unreliable determinations, hence identification based on phylogenetic studies and/or genitalia is strongly recommended not only for the three studied species but also more broadly within the genus. Furthermore, we show that some of the characters generally used in the identification of these three Melitaea species should be avoided in future.


Introduction
The genus Melitaea Fabricius, 1807 (Lepidoptera: Nymphalidae) as it is known today (Leneveu et al. 2009), comprises approximately 80 species, all restricted to the Palaearctic region. A recent phylogenetic analysis showed that the genus can be divided into two clades ('Melitaea clade' and 'Didymaeformia clade') with high branch support from both maximum likelihood and Bayesian analyses (Leneveu et al. 2009). Within the 'Melitaea clade' five subclades have been recognized (sensu Leneveu et al. 2009): 'cinxia', 'diamina', 'arcesia', 'minerva' and 'athalia' groups. The latter subclade comprises fourteen species, three of which are present in the area of north-western Balkans: Melitaea athalia (Rottemburg, 1775), Melitaea aurelia (Nickerl, 1850) and Melitaea britomartis Assmann, 1847. These three species belong to two different monophyletic groups within the subclade 'athalia', one comprising M. aurelia and five other species, and the other M. athalia, M. britomartis and six other species (Leneveu et al. 2009). The split between the latter two monophyletic groups has been estimated to have Nota Lepi. 37(1) 2014: 75-90 | DOI 10.3897/nl.37.7966 occurred about 11 million years ago (Serravallian period) while the subsequent speciation occurred during the Messinian Period (5.3-7.1 Mya) (Leneveu et al. 2009), when the continuous decrease of temperature coincided with the aridification of the climate in Eurasia and a subsequent expansion of grassland of C4 plants (e.g. Maki et al. 2003;Quade et al. 1995). Although the separate evolution of those species has been in progress for a long period of time, the clear genetic differentiation (Bátori et al. 2012a) is not reflected in their external morphology (Tolman & Lewington 2008).
The external morphology of the three species is highly similar, rendering the accurate identification based on external morphology questionable (Tolman & Lewington 2008). There are several morphological characters that have been proposed for species identification of M. athalia, M. aurelia and M. britomartis: the coloration of the marginal line of the hindwing underside; the spacing of the marginal, submarginal, postdiscal and discal line on the forewing upperside; and the colour of the submarginal spots on the hindwing underside (Tolman & Lewington 2008). The high variability of those characters, however, sometimes prevents proper identification (Tolman & Lewington 2008). The genitalia of these species are highly species specific, and have been known to be a good identification tool (Urbahn 1952;Paulavičiūtė & Tamutis 2009). Nevertheless, even when genitalia based identification is used, only males can be unambiguously recognized to a species rank, whereas this is not always possible in females (e.g. see Urbahn 1952).
Melitaea athalia is a trans-European species (except for south-western Europe, where Melitaea celadussa Fruhstorfer, 1910, a former subspecies of M. athalia occurs (Leneveu et al. 2009)). M. aurelia and M. britomartis are absent from the majority of SW Europe (including the Iberian and Apennine peninsulas). M. britomartis is absent also from most of Central Europe being present only in eastern parts of Europe. A part of its distribution, however, reaches the NW Balkans -Slovenia, Croatia and Bosnia and Herzegovina (Koren & Jugovic 2012) -where these three species live in sympatry. In this area the distribution of M. britomartis is separated from the center of its distribution (Tolman & Lewington 2008).
The ecology of the three species in the area of NW Balkans (Slovenia) is similar in terms of flight period, altitudinal distribution and habitat requirements including foodplant (Koren & Jugovic 2012;Verovnik et al. 2012). Since their habitat requirements and foodplants partly overlap, they often occur in the same habitat (Bátori et al. 2012a;Koren & Jugovic 2012).
Morphological traits can reflect either historical isolation and/or local adaptation despite recurrent gene flow (Alexandrino et al. 2005). Some taxa remain morphologically unidentifiable despite detailed morphometric analysis, hence a consideration of a large set of morphometric characters from different anatomical regions may greatly increase the chances of revealing taxonomic differences within seemingly cryptic or morphologically highly variable species (e.g. Jugovic et al. 2011).
We wanted to check the concordance in the amount of morphological and genetic differences among the three closely related species. For this purpose a morphometric approach with subsequent multivariate statistical analysis was used in the M. athalia complex for the first time. We evaluated the following hypotheses: (1) the genetic divergence of the three species is accompanied by corresponding morphological differences exceeding the intra-species variability and (2) the level of accompanying morphological differences is concordant with the level of molecular distance between the taxa, meaning that M. aurelia should be the most distant in morphological space from the other two species.

Samples and species identification
Samples of the three species of Melitaea were collected over the last 40 years all over the NW Balkans (see Koren & Jugovic 2012 for the list). Altogether, samples from 56 localities were collected, sample sizes ranging from 1 to 12 for males and 1 to 3 for females. The genitalia of each specimen were isolated using the standard procedure for genitalia isolation (see Koren & Jugovic 2012 for details). Specimens were attributed to a species according to the genitalia structure that is highly species-specific in males (Urbahn 1952;Paulavičiūtė & Tamutis 2009). Using these methods we were able to identify all of the collected males. A total of 42 males and 24 females of M. athalia, 29 males and 6 females of M. aurelia and 24 males and 4 females of M. britomartis were used in further morphometric analysis.
Separation of these three species is possible using male genitalia as follows (see also Urbahn 1952): 1. uncus absent -> M. aurelia; 2. uncus present, long and slender, its processus longer than wide -> M. athalia; 3. uncus present, short and robus, its processus almost as long as wide -> M. britomartis. In females there are slight differences among the species, but not all specimens can be reliably identified in this way (Urbahn 1952;Bátori et al. 2012a).

Morphometric methods
Specimens used in this study were mounted and photographed on a millimeter grid from the same angle (90°) using a DSLR camera (Canon 450D). Subsequently, 25 metric characters were measured from the photographs using freeware ImageJ (Abramoff et al. 2004), 14 on the forewing and 11 on the hindwing (Table 1, Fig. 1). The measurements included a wide range of morphometric characters (i.e. lengths, surfaces and angles between measured distances of different parts of both wings). Moreover, one categorical character, i.e. marginal line colour in contrast to lunular colour (for categories see Table 1 and Appendix A) was also recorded. In order to exclude the impact of the animal's size, all metric characters were transformed into 18 additional ratios ( Table 1). The photographs of both wings were used for comparison and investigated in detail in order to find some other descriptive morphological species-specific characteristics and to describe intraspecific variability. Although no body asymmetry was detected, all characters were measured on the right side of an animal in order to exclude any possible influence of this phenomenon. The left side, however, was used in rare occasions when the right side was damaged.

Statistical analysis
Since in our samples males were prevalent and not all females could be reliably identified to a species, females were excluded from subsequent statistical analyses. With this approach, we also avoided the influence of possible sexual dimorphism. Only for the analysis of the marginal line colour were a few unambiguously identified females added to the sample (no sexual dimorphism was noticed in this character). We also had to remove from the analysis three aberrant individuals which lacked measuring points for some characters, hence they are commented on separately.
For each species, the Kolmogorov-Smirnov test (at p = 0.05) and Normal Q-Q plots were used to examine the normality of the data distribution, and the homogeneity of variances was evaluated visually through scatterplots. The multivariate analysis of variance (MANOVA, at p = 0.001) was used to test for significant differences between species. One-way analysis of variance (ANOVA) was used to assess the variation within a species for each character, significant variation in a character being accepted if the between species variation was significant at p < 0.001. Pearson's correlation coefficients (r) were computed to evaluate the extent to which each character contributes unique information; only one character was chosen to represent a pair or a group of characters where |r| > 0.9. The Durbin-Watson test (at p = 0.05 and p = 0.01) was applied to test for possible spatial (latitude, longitude, altitude) and temporal (year of collection) autocorrelations of morphometric data (Savin & White 1977;Farebrother 1980) in each of the three species. In subsequent multivariate statistical analyses, only selected characters (according to the limitations listed above) were used.  Table 1. For categorical character, see Appendix A.
Multivariate Principal Component Analysis (PCA) was used to identify the structure of our data i.e. to detect the possible influence of the species specific characteristics. Also, Discriminant Function Analysis (DFA) was carried out to examine possible separation of the three species. In DFA, the contribution of each species was weighted according to its sample size (number of specimens). Post-hoc Games-Howell and Bonferonni tests were performed to assess the rates of morphological divergences between pairs of species in details. The analyses were performed using Microsoft Excell (2010), SPSS 14.0 for Windows (2005) (Norusis, 2005) and Past (PAlaeontological STatistics) software (Hammer et al. 2001).

Results
The multivariate analysis of variance (MANOVA) showed significant differences between species (p < 0.001). Out of 43 morphometric characters, 12 characters (without a single ratio) were selected after Kolmogorov-Smirnov (p < 0.05), ANOVA (p < 0.001) and Pearson correlation (|r| < 0.9) tests. No spatial or temporal autocorrelation was detected in these characters neither at p = 0.01 nor at p = 0.05 (Durbin-Watson test). Statistic description of selected 12 characters (mean, standard error, 12th percentile, 88th percentile and extreme values) is presented in Table 2. When the Pearson correlation test was set to |r| < 0.7, only two morphometric characters (forewing height [FWH] and distance A4) were left after stepwise exclusion of characters, meaning that correlations among most pairs of morphological characters were high.
The principal result of the PCA (Fig. 2) run on 95 males from the NW Balkans, using 12 morphometric characters (Table 2), revealed a significant overlap of the three species along the first two principal components (PCs). PC 1 explained almost 70% of the total variance and PC2 explained an additional 8%. All characters are positively correlated to PC1, starting with the smallest M. aurelia and finishing with the largest of three species, M. athalia. Nevertheless, the overlap between the three species in the middle section along PC 1 is significant.
In the DFA (Fig. 3) run using the same morphometric characters as in the PCA, 95 males were analyzed to (1) provide adequate species-grouping according to external morphological characters and (2) detect intra-species variability of selected characters. Two Discriminant Functions (DFs) explain total variance, DF 1 explaining over 75% of the total variance. Species differ significantly only along DF 1 (p < 0.001), and only M. athalia differs significantly from the other two species (Games-Howell Test, p < 0.001). The characters most correlated with DF 1 (indicated by the standardized discriminant function coefficients: DC > 0.75, Table 3) are: forewing height (FWH), forewing lunule surface (FWLS), forewing lunule height (FWLH) and hindwing surface (HWS). For the differences among the three species in these four characters, see Fig. 4 and Table 2.
Although the DFA aims to find the differences between the a priori defined groups, the misclassification rates were high, especially when cross-validation process was employed (Fig. 5). Considering both DFs, only 67.4% specimens were correctly classified, and in the cross-validation procedure, the percentage dropped down to only 52.6% of correctly classified specimens. We then repeated DFA with the same twelve morphometric characters, adding the coloration of the marginal line in contrast to the lunular colour on the underside of the hindwing, not considering that this character was deviating from normal distribution (Kolmogorov-Smirnov test, all data pooled: Z = 2.642; p < 0.001). However, this did not improve the classification significantly (original grouping: 70.5%; cross-validation procedure: 55.8%).
In addition to the characters used in the DFA, two characters frequently mentioned in the literature as diagnostic for the recognition of (some of) Melitaea species should be mentioned.  For the characters like the spacing of the marginal, submarginal, postdiscal and discal line on the forewing upperside, the colour of the submarginal spots on the hindwing underside, the coloration of hairs on the palps as well as that of the marginal line in contrast to the lunular colour on the underside of the hindwing, no consistency with species attribution was found in our examination. For the latter, slight differences were found between the three species; however, species could not be identified with certainty using the lunular coloration on the hindwing underside due to the high intraspecific variability of this character. The coloration of the marginal line is in general equal to the coloration of the lunules in M. athalia, much darker in M. britomartis and only slightly darker in M. aurelia. Nevertheless, all of these categories were noticed in each species (Fig. 6, Appendix B).
The wing pattern and coloration show high intra-and interspecies variability. While in some species dark melanistic forms are common (e.g. in Melitaea britomartis ssp. michielii from Slovenian Karst and surroundings (Carnelutti 1992)), other peculiar forms also had been found. In some of these, the pattern deviated from normal form so much so that in these animals some  (HWS). Statistically significant differences between pairs of species (Bonferroni test) are shown with capital (p < 0.001) and small (0.001 < p < 0.05) letters. When only a trend in differences between species is shown (0.05 < p < 0.1), small letters with apostrophe are used. metric characters could not be measured. We found three such M. athalia specimens (i.e. with extremely reduced markings) from Kamniški Vrh and Smrekovec (both in Slovenia) and Vugrovec (Croatia).

Discussion
The genitalia structures (Urbahn 1952) and the phylogeny based on three genes (mitochondrial gene cytochrome oxidase subunit I, and two nuclear genes, elongation factor-1a, and wingless; Leneveu et al. 2009) of the three Melitaea species are in agreement, clearly showing that M. aurelia is more distant from the other two species. In our analysis this is proved by a genitalia based species identification of males without any questionable cases. However, the phylogenetic relationships and the differences in the genitalic morphology are not reflected in the external morphology of the three species' wing patterns.
All analysis conducted in this study indicates the important differences among the three species in the size of most of the measured characters, however, only in average values. Due to high variability of these characters, the clinal variation among species has been noticed (see the results of multivariate analyses). According to literature (e.g. Tolman & Lewington 2008), M. athalia is indeed the largest of the three species, with forewing height (FWH, Table  2, Fig. 4) being the most obvious character reflecting the difference in size from our study. Without exception, all measured characters had the smallest average values in M. aurelia on one side of the clinal variation and the largest in M. athalia on the other. High correlations among most of the characters show the stability of their wing shape, which can therefore be represented by a subset of all the characters. The wing shape stability is also supported by the exclusion of all ratios (that partly exclude the size impact and describe the wing shape) in ANOVA tests. In order to minimize the number of highly correlated characters, only one character was chosen from a pair or group of highly correlated characters. This goal was only partially achieved as further stepwise removal of characters (with Pearson correlation coefficients 0.7 < |r| < 0.9) would result in only two weakly correlated remaining characters (|r| < 0.7) and prevent the implementation of the multivariate analyses. Although the differences are small, M. athalia is the most distant from the other two in a powerful discriminant function analysis. All except one character that were included in the multivariate statistical analyses significantly separate M. athalia from the other two, and only slight differences in just one (out of 43) metric character (forewing lunule height, FWLH; ANOVA: 0.05 < p < 0.1) between M. aurelia and M. britomartis were shown. Although the differences are small, this further supports the discrepancies between the external morphology and phylogenetic results (sensu Leneveu et al. 2009). Considering the phylogenetic position of the three species, M. athalia should resemble M. britomartis more than M. aurelia but in our analysis the latter two were clustered more closely together. According to the phylogenetic data, the first split within the 'athalia' group that separated the clade containing M. aurelia from the clade containing the two other species occurred in the Tortorian period (approximately 11 MYA) and the later split from which also M. athalia and M. britomartis emerged happened during Messinian period (approximately 7.1 -5.3 MYA; Leneveu et al. 2009). Hence, the three lineages that emerged -M. athalia, M. aurelia and M. britomartis -underwent at least 5 million years of separate evolution. Despite such a long time of separation, no obvious external morphological differences have evolved. This is somewhat surprising when compared to the situation in some of their close relatives. For example, clear differences in external morphology have emerged in some species of Melitaea, although they have become separate evolutionary lineages more recently than M. athalia, M. aurelia and M. britomartis (e.g. M. asteria-M. aurelia, also within the same 'athalia' group; Tolman & Lewington 2008;Leneveu et al. 2009).
The allozyme polymorphism of these three species was studied in samples from the Carpathian basin, and revealed that M. britomartis and M. aurelia are more closely related, while M. athalia appears to be a further relative (Bátori et al. 2012a). This result contradicts the phylogenetic analysis of Leneveu et al. (2009). The latter analysis, however, finds support also in the genitalia structures as follows: in males of M. athalia and M. britomartis, the spines on the uncus are well developed, whereas they are completely absent in M. aurelia (Urbahn 1952). It should not be a surprise that the level of differentiation in genitalic structures corresponds to the phylogenetic position of the species since sexual differentiation usually represents one of the most rapid and obvious taxon-specific events during speciation in many animal groups, invertebrates in particular (Mayr & Ashlock 1991). In addition, the genitalic structures are in comparison to the external morphology less subjected to environmental factors (Cesaroni et al. 1994;Dapporto et al. 2011). For example, M. celadussa, recently elevated to a species level by Leneveu et al. (2009) from a subspecies rank of M. athalia, has not yet been reliably recognized on the basis of wing pattern and coloration, whereas some differences from M. athalia were found in male genitalia alone (e.g. see Higgins & Riley 1978;Leneveu et al. 2009). As long as external differences are not found, this species pair could be treated as cryptic (sensu Hawksworth 2010: "populations which are phylogenetically distinct, but distinguished by molecular or other features that are either not evident macroscopically or generally overlooked").
In contrast to the allozyme study, the phylogenetic analysis used as a framework for the explanation of our results (Leneveu et al. 2009) (1) included a vast majority of known species of Melitaea (whereas in Bátori et al. 2012a only M. athalia, M. aurelia and M. britomartis were included) and (2) three different genes provided support for the same topology within the 'athalia' group discussed herein. Hence, no serious consideration was given to the study of Bátori et al. (2012a).
No external characters were proven to be reliable for species delimitation even though the coloration of the marginal line of the hindwing underside shows the trend towards the correct identification. Although the majority of specimens correspond to the character states given in literature for a particular species, all three presumably species-specific categories of this character (Tolman & Lewington 2008) used in our analysis were present in each of the three species. Moreover, the percentage of misidentifications when this is the only character used would be very high (see Fig. 6). Other external characters suggested in literature (see Introduction) show even higher rates of variability, and no species-specific correspondence was found. Even when the twelve (out of 43, see Tables 1 and 2) most powerful characters are simultaneously used (see results of the DFA, Fig. 2) the misidentification rate reaches almost 50% in the cross-validation procedure (Fig. 5), rendering the appropriate identification of these species highly unreliable. To further support the high morphological variability within a single Melitaea species, as an example we should mention a comparison of populations of M. athalia from the Caucasian basin (Bátori et al. 2012b), where small, albeit significant differences in wing morphology were found among some of them (all contributing to a high intraspecies variability). Although partly reflecting the molecular differences between the populations, these differences cannot be assigned to a genetic background alone. They rather may through selection also be a result of an adaptation to local conditions (Bátori et al. 2012b), such as food availability, climate and microhabitat selection. This phenomenon is not uncommon since high levels of phenotypic plasticity are important for a response that ensures the survival of a population exposed to unstable environmental factors (Shapiro 1976;Brakefield & French 1999). High species plasticity together with accommodation to local conditions can in practice result in description of many subspecies and morphological forms within M. athalia as well as within some other Melitaea species. Nonetheless, the differences among populations in reality might have resulted from locally specific selection pressure (e.g. see Tolman & Lewington 2008). On the other hand, the three species exploit similar (and partly overlapping) resources (e.g see Verovnik et al. 2012;Koren and Jugovic 2012), which is probably the reason for the high similarity among them. It is important to note that the regional pattern of differentiation in M. athalia from the Caucasian basin was less expressed in genitalia than in wing characters (Bátori et al. 2012b).
The high variability of these species is further shown in their qualitative characters, for example in wing pattern and coloration. In M. athalia three out of more than 120 sampled specimens had a very peculiar coloration, with no visible marginal lines. Albeit unrecognized, we believe that differences in (some of) local conditions clearly show the importance of environmental factors for the wing morphology. This further demonstrates the importance of genitalia or phylogeny based identification in the herein investigated Melitaea species. Moreover, the latter two approaches should have a significant advantage also for the identification of (some) other species of the genus that should undergo thorough revision in the future (cf. Leneveu et al. 2009).

Conclusions
The morphometric analysis of three closely related Melitaea species (M. athalia, M. aurelia and M. britomartis) revealed the following: (1) small albeit significant differences in wing morphology exist among the three species; (2) only characters describing the size of the species with no characters describing their wing shapes (ratios) were statistically important for the separation of the species; (3) whereas the structure of the male genitalia and the phylogenetic position of the three species are concordant, the sequence of phylogenetic splits is not reflected in the rate of external morphological differences among them; (4) our study represents another example where external morphology based identification would lead to highly unreliable determinations, hence a use of genitalia based identifications is strongly recommended.