Identification of the species of the Cheilosia variabilis group (Diptera, Syrphidae) from the Balkan Peninsula using wing geometric morphometrics, with the revision of status of

The present study investigates phenotypic differentiation pat - terns among four species of the Cheilosia variabilis group (Dip-tera, Syrphidae) using a landmark-based geometric morphomet - ric approach. Herein, wing geometric morphometrics established species boundaries that confirm C. melanopa and C. redi stat. nov. as evolutionarily independent entities or separate species within the previously described C. melanopa species. Seventeen samples of C. lasiopa , C. variabilis , C. melanopa and C. redi from 11 localities on the Balkan Peninsula were analysed. Based on wing shape and size, both species boundaries and intraspecific divergent units were delineated indicating the possible presence of evolutionarily independent units within the analysed taxa. Sig nificant differences in wing size were obtained among the ana - lysed species and canonical variate analysis showed that wing shape was sufficiently different to allow the correct classification of 99.7% and 100% individuals of males and females, respec -tively. Considerable wing size variation was detected among males and females of conspecific populations of C. melanopa and males of C. variabilis . Contrary to partial discrimination based on the male wing shape among conspecific populations of C. redi and C. variabilis , clear separation of conspecific popula tions of C. melanopa (for both sexes) and females of C. redi and C. variabilis was observed. The UPGMA cluster analysis based on squared Mahalanobis distances revealed a close similarity be - tween C. melanopa and C. redi , whereas C. lasiopa was the most divergent species. Results presented in this study utilize wing shape and size as new taxonomic characters in delimitating the closely related species and populations of the Cheilosia variabi- lis group.


Introduction
The Balkan Peninsula harbours ecologically, morphologically and genetically diverse species of hoverflies (Diptera, Syrphidae), including a number of endemic and cryptic species (e.g., Milankov et al., 2008Milankov et al., , 2009. A pattern of genetic and phenotypic divergence of taxa reflects the complex biogeographical history of the region (Hewitt, 2000(Hewitt, , 2004Pincell et al., 2005) coupled with the evolutionary history of population bottlenecks followed by range expansion Ståhls et al., 2008). In addition, hidden phenotypic and genetic units and evolutionary diversification among populations of the hoverflies on the Balkan Peninsula have been observed Francuski et al., 2009). Therefore, study of hoverfly taxa subdivided into genetically and phenotypically distinct populations has implications in taxonomy, systematics and conservation.
The genus Cheilosia Meigen, 1822 (Diptera, Syrphidae) comprises 300 species distributed in the Palaearctic region, more than 80 Nearctic and at least 50 from the Oriental region (Ståhls et al., 2004). High intraspecific diversity, morphological similarity of the Cheilosia species, presence of cryptic species and numerous synonyms create various problems in taxonomy and systematics studies (e.g., Ståhls et al., 2008;Ludoški et al., 2008).
The subjects of the present study are C. lasiopa Kowarz, 1885 andC. variabilis (Panzer, 1798) and two Balkan subspecies of C. melanopa (Zetterstedt, 1843) (Vujić, 1996), all of which form the C. variabilis group of species (van Veen, 2004). Using a combination of morphological characters for subgeneric classification of the genus Cheilosia, Barkalov (2002) placed the four closely related taxa into subgenus Cheilosia s. str, while some authors kept these taxa in another Cheilosia group (Becker, 1894;Haarto and Kerppola, 2007). The parsimony analysis of the mitochondrial DNA (mt-DNA) sequences of the cytochrome c oxidase subunit I (COI) and nuclear 28S ribosomal RNA gene (Ståhls et al., 2004) revealed subdivision of the genus Cheilosia into 13 subgenera according to Barkalov (2002). The phylogenetic analysis classified the C. variabilis and C. melanopa species within the monophyletic clade of the Cheilosia s. str. subgenus (Ståhls et al., 2004). Van Veen (2004) provided a morphology-based key for accurate identification of the C. variabilis group taxa. The taxa belonging to the C. variabilis group share the following characters: eyes haired, face with long hairs on the sides of the facial knob, equal to or longer than hairs on eye rim, and black legs (van Veen, 2004). The C. variabilis group has been known as a challenging group of hoverflies for studies from taxonomical and biogeographical perspectives. For example, the analysis of the C. melanopa specimens from the Balkans and adjacent regions indicated the existence of two subspecies: C. melanopa melanopa described as a smaller subspecies predominantly covered with black hairs, and C. melanopa redi which is larger and covered with reddish hairs (Vujić, 1996). Cheilosia m. melanopa is distributed in areas that are more northern and at higher altitudes in the southern part of this range, while C. m. redi appears in Central and Southern Europe at low altitudes (Vujić, 1996). Both C. variabilis and C. lasiopa are distributed in Northern, Central and Western Europe, Siberia and Mongolia, while C. variabilis also occupies some areas in Africa (Morocco) (Speight, 2007). While the larval development of the subspecies is unknown, the larva of C. variabilis have been described in the roots of Scrophularia nodosa (Dušek, 1962). The immature stages of C. lasiopa initially live in the leaf and stem bases of Plantago lanceolata and later in the rootstock of that plant (Stuke and Carstensen, 2000).
Landmark-based geometric morphometric technique is considered to be the most rigorous morphometric method (Gilchrist et al., 2000;Debat et al., 2003) and a useful tool in accurate identification and delimitation of numerous insect species (e.g. Debat et al., 2006;Jirakanjanakit et al., 2008;Marsteller et al., 2009). Since it has been known that wing shape exhibits a high heritability in nature (Bitner-Mathé and Klaczko, 1999;Moraes et al., 2004), wing morphology is of a primary importance to entomologists interested in systematics. Wing size and shape have different genetic properties, size heritability being generally low, while wing shape is less sensitive to environmental changes and highly heritable (Bithner-Mathé and Klaczko, 1999;Birdsall et al., 2000). However, wing size is usually associated with body size, and thus potentially linked to a number of fitness components (Reeve et al., 2000). Assuming that wing shape is related to flight ability, wing shape might also influence fitness (Kölliker-Ott et al., 2003). Although the genetic determination of the final details of morphometric characters is still relatively poorly known, Birdsall et al. (2000) suggested that genes that regulate wing shape are more tightly connected than those for wing size, and the former are probably associated with gene determination of wing veins. Therefore, investigation of the wing measurements has both taxonomic value (cf. cryptic species of the Drosophila buzzatii Patterson & Wheeler, 1942 cluster, Moraes et al., 2004), and importance in study of adaptive evolution (Soto et al., 2008;Bischoff et al., 2009).
To date, wing geometric morphometrics have been successfully used to resolve taxonomic uncertainties of some species groups (Merodon avidus (Rossi, 1790), Milankov et al., 2009; Cheilosia canicularis (Panzer, 1801), Ludoški et al., 2008), and to define divergent phenotypic units in hoverflies (the Merodon ruficornis Meigen, 1822group, Francuski, et al., 2009. Moreover, the only morphological traits that might be used in delineating some cryptic species (a group of reproductively isolated but morphologically inseparable species) such as Merodon bicolor Gill Collado, 1930 vs M. avidus  Considering all of these, we examined the usefulness of wing landmarks for species delimitation and for quantification of intra-and interspecific variation of species of the C. variabilis group on the Balkan Peninsula. The goal of the present study was identification of taxa within the C. variabilis group using landmarkbased geometric morphometrics on wing traits. By introducing morphometric techniques, we also re-evaluated the taxonomic status of C. melanopa mela nopa and C. melanopa redi. Then, we used wing traits to test the differences among conspecific populations and recognize phenotypic units within the species. Thus, phenotypic variation across spatially and temporally fragmented populations was used to examine the utility of wing size and shape for detecting levels of intraspecific variation as well. Finally, phenetic relationships among species of the C. variabilis group based on our findings were discussed. The results reported herein provide a taxonomic framework and give us better insight into intra-and interspecific variation of the C. variabilis group.

Sample collection
During the period 1981-2000, 475 specimens of the C. variabilis group species were collected from 11 localities on the Balkan Peninsula ( Fig. 1, Table 1) for use in morphological studies. The study material considered 371 male and 104 female specimens from 17 populations. Due to the sexual wing shape dimorphism  these species exhibit, analyses were performed separately on males and females.
In the preliminary analyses of phenotypic variability, morphological diagnostic traits (see Results) allowed identification of C. melanopa melanopa and C. m. redi as two separate species. Morphological differences in wing shape and size were considered substantial for species level differentiation. Prior to geometric morphometric analysis, identification of specimens was performed using other morphological characters of adults (Vujić, 1996).

Morphometric and statistical analysis
Wing size and shape variation was observed from 475 specimens of the C. variabilis group, using the landmark-based geometric morphometric method (Bookstein, 1991;Rohlf and Marcus, 1993;Marcus et al., 1996;Adams et al., 2004). The right wing of each specimen was mounted separately in Hoyer's medium between microscope slides and digital images were captured using a stereomicroscope Leica MZ12.5 and attached video camera Leica DFC320 connected to a PC computer. Sixteen landmarks positioned at vein intersections or terminations were collected using TpsDig 1.40 and expressed as x,y coordinates in a Cartesian space (Rohlf, 2004) (Fig. 2).
Wing size was computed as centroid size (CS; the square root of the sum of squared distance between each landmark and the wing centroid), an isometric estimator of size. One-way analysis of variance (ANOVA with post hoc Tukey HSD test) was used to test differences in CS among species and populations. To examine wing shape variation, the digitized landmark data were subjected to Generalized Procrustean superimpositions to standardise size of the landmark configurations and remove differences due to translation and rotation (Rohlf and Slice, 1990), and then a thin plate spline analysis was done. The resulting matrix (w; 'weight matrix' of Rohlf et al., 1996) containing the non-uniform and uniform shape variables were calculated was used for a discriminant analysis combined with canonical variate analysis (CVA) to examine the pattern of among and within species/ population variation in total shape space. We used the percentages of correct classification to evaluate the discriminatoring power of wing shape. Both centroid size and w matrix were obtained utilizing TpsRelw software (version 1.44; Rohlf, 2006) and all statistical analyses were calculated using Statistica for Windows (version 8.0).
Shape changes associated with canonical axes were visualized as deformations obtained by regressing w matrix onto individual canonical axis scores (Rohlf et al., 1996). Multiple regressions and visualization of deformation grids were computed using TpsRegr 1.31 (Rohlf, 2005). The phenetic relationships among species were determined by UPGMA cluster analysis (Sneath and Sokal, 1973) based on squared Mahalanobis distances derived from CVA.

Cheilosia melanopa redi was described from Central
Europe and the Balkan Peninsula as a subspecies of C. melanopa based on clear differences in coloration of hairs on the body. Comparing to C. melanopa (8-10mm), C. redi is a larger species (10-12mm). In addition, diagnostic characters of C. redi are: males -body hairs reddish (from yellow-reddish to brown-reddish), black hairs can be found on tergites, sometimes on mesoscutum (on postalar calli) and scutellum margin, wings with brown veins; females -body covered with adpressed yellow-reddish hairs, legs paler than in C. melanopa, tibiae pale in basal third and on apex. Specimens of C. melanopa are characterized by the following features: males -body hairs black, except the combination of black and pale hairs on tergites, differently arranged, and pale hairs laterally and adpressed black hairs in the middle of tergites, tibiae brownish basally; females -face with pale hairs, mesoscutum with pale and black hairs (longer are black), abdomen with pale hairs laterally and adpressed black hairs in the middle of tergites; tibiae brownish basally (Vujić, 1996). Contrary to differences of body hairs, differential features in male genitalia structure, which was proved to be the main taxonomic character for recognizing the species of the genus Cheilosia (Vujić, 1996), appeared to be identical in both subspecies (Table 2; Fig. 3).
Cheilosia m. redi is distributed in the Pannonian plain, lowlands and hills of the central Balkan Peninsula. Cheilosia m. melanopa occurs on high mountains in Northern and Central Europe. A study conducted during 1998 discovered sympatric and sinhronic occurance of both, C. m. melanopa and C. m. redi 'subspecies', only on the Durmitor Mt in Montenegro (locality Kanjon Sušice). Results of wing geometric morphometrics obtained in this paper (see below) strongly support an independent taxonomic position at the species level of C. melanopa and C. redi. Based on the observed morphometric data, previously recognized morphological differences (Table 2) and co-existance of both taxa, we propose formal elevation of C. redi to the species level.

Diagnostic characters
The analysis of variance (ANOVA) of wing centroid size (CS) revealed significant differences between the   (Fig. 4, Table 3).
The first canonical axis accounting 77.8% of total shape variation clearly separated C. melanopa from C. redi with 100% classification success. The thin plate spline visualization showed that male shape differences were associated with displacement of landmarks 5, 6, 7 and 8 (Fig. 5c), while for female specimens, landmarks 1, 2, 6, 7 and 8 have the highest contributions for C. melanopa/C. redi wing shape differences (Fig. 5d).

Identification of species
The ANOVA on CS of the wings showed a highly significant result among species of the C. variabilis group when sexes were considered separately (male: F 3,367 = 131.41, P<0.001; female: F 3,100 = 72.47; P<0.001; Tukey HSD test: P<0.001, except for males of the C. lasiopa/C. melanopa species pair) (Fig. 4, Table 3).
CVA with species as grouping variable performed on the w matrix gave satisfactory results for overall species comparison. In males, the three axes extracted by CVA were highly significant and contribute to shape differentiation among species (Wilks' Λ = 0.059, F 84,10 = 19.04, P<0.001) (Fig. 6a). The percentage of individuals correctly classified was in total 99.7%; except for the C. melanopa species (94.7%) analysed species were classified with rate of 100%. Similarly, CVA showed a clear interspecific discrimination of females within the C. variabilis group (Wilks' Λ = 0.003, F 84,22 = 15.87, P<0.001) (Fig. 6b).
There was a substantial divergence in wing shape to classify correctly 100% of all female individuals.

Recognition of phenotypic units
The character CS varied significantly between two C. melanopa populations; wings of Durmitor Mt specimens were considerably larger than C. melanopa wings from Šar Mt, both in males (F 1,17 = 13.04, P<0.01) and females (F 1,12 = 93.69, P<0.001). To assess the wing shape variation between studied populations, we used CVA conducted on the w matrix. Although CVA evidenced no significant population differences (Wilks' Λ = 0.00003, F 84,6 = 2.59, P = 0.09), percentages of correct classification were 100% (Fig.  7a).The first two canonical axes were significant and explained 87.6% and 10.2% of the total shape variation. The second axis (CV2) contributed to shape differentiation between two populations. Deformation grids which represents deformations along the CV2 demonstrated that population differences primarily associated with displacement of landmarks 4, 7, 8 and 12 (Figs 7b, 7c).
The ANOVA showed no wing size differences between conspecific populations of C. redi (DUB, FG,   Females from three analysed populations clearly differed in their wing shape, while male specimens were partially overlapped (Fig. 8). Contrary to females (F 8,48 = 2.11, P = 0.53), significant wing size differences between male specimens of  the C. variabilis populations were found (F 10,31 = 12.51, P<0.001). Males from Kopaonik and Suva Mts had considerably larger wings than males from Eastern Alps, Dubašnica, Durmitor and Fruška Gora Mts (P<0.05, Tukey HSD test). The ANOVA also exhibit significant size variation (P<0.05, Tukey HSD test) between the populations of the Fruška Gora Mt/ Stara Mt and Vršačke Mt/Suva Mt pairs (not shown). CVA revealed population differences within both sexes. In males, six of ten axes extracted by CVA were highly significant, contributed to shape differentiation among populations (Wilks' Λ = 0.141, F 280,27 = 2.22, P<0.001) and led to correct classification rates between 20% and 100% (not shown). When female specimens were compared (Wilks' Λ = 0.0001, F 224,19 = 2.03, P<0.001), significant wing shape differences allowed 98.3% of individuals to classify correctly with 100% classification success for individuals of all analysed populations except those from Fruška Gora Mt (83.3%) (Fig. 9).

Phenotypic relationships
The UPGMA cluster analyses of the squared Mahalanobis distances computed from w matrix of female shape parameters (Table 3) clustered C. melanopa and C. redi in the same branch, while C. lasiopa was the most divergent species (Fig. 10). The UPGMA dendrogram obtained from male shape parameters provided the same phenetic relationships among analysed species (not shown).

Species delimitation
This study is the first to quantify subtle wing geometric morphometric traits among populations of two previously nominated subspecies, C. melanopa melanopa and C. m. redi from the Balkans (Vujić, 1996). To date, subspecies discrimination was solely based on morphological characters such as body hairs (black and reddish in C. melanopa melanopa and C. melanopa redi, respectively) and wing veins (dark and brown in C. m. melanopa and C. m. redi, respectively) (Vujić, 1996). Although the morphology of male genitalia has been widely used in identification of hoverfly species, the diagnostic value of it was not recognized to be important in delimitation of the previously nominated subspecies. Similarly, lack of taxonomic value of male genitalia has been observed for cryptic taxa within the morphologically defined spe-cies Cheilosia vernalis (Fallen, 1817)  and closely related species pair Cheilosia nigripes (Meigen, 1822) and C. vicina (Zetterstedt, 1849) (Barkalov and Ståhls, 1997). However, based on wing shape and size presented herein the status of the focal taxa was clarified as evolutionarily independent biological entities, the C. melanopa and C. redi species. Wing centroid size (C. melanopasmaller, and C. redi -larger wings) and wing shape were found to be significantly different between the analysed species (Table 3). Moreover, observed divergent wing shape phenotypes in sympatric populations from Durmitor Mt clearly demonstrated C. melanopa and C. redi as evolutionarily independent entities. Wing shape differences were associated with displacement of outer landmarks in females (landmarks 5, 6, 7, 8) and landmarks that located at proximal and posterior parts of wing in males (1,2,6,7,8), all which influence wing's width and length. The consistency of this pattern of interspecies variation is also evident in other Cheilosia species (J. Ludoški, unpublished data). Therefore, associated with previous findings, the current study has located the wing's region of taxonomic importance. The occurrence of similar patterns of wing variation in different species of the genus Cheilosia might reflect conservative genetic constraints, but only further analysis will offer a precise explanation.
We quantified morphological differences of the species within the C. variabilis group using wing geometric morphometrics as well. Considerable interspecific differences in wing size were estimated; C. melanopa had the smallest wings, followed by C. lasiopa, C. redi and C. variabilis. Documented wing  shape variation also revealed clear phenotypic differentiation among the species of the C. variabilis group. Such subtle wing morphometric distinctiveness among closely related species was also found in other species groups within the family Syrphidae (Ludoški et al., 2008;Francuski et al., 2009;Milankov et al., 2009), and demonstrated wing traits as a powerful taxonomic tool.

Intraspecies phenotypic diversity
Our findings revealed high phenotypic divergence among conspecific populations within the C. melanopa, C. redi and C. variabilis species. Contrary to wing size, wing shape was found to be significantly different between females of C. redi. The extent of wing variation between populations of C. melanopa originating from Durmitor and Šar Mts included discrete variation of wing shape primarily associated with displacement of landmarks 4, 7, 8 and 12. Deformation grids showed that most morphological variation occurred in the distal-posterior part of wing, a region that has been supposed to be important in aerodynamic performance and courtship song in dipteran species (Birdsall et al., 2000). We also found significant interpopulation wing size variation between males of C. variabilis. In comparison, the wing shape of females of C. variabilis revealed phenotypic interpopulation differentiation. Hence, we considered that allopatrically distributed populations of the C. melanopa, C. redi and C. variabilis species have likely maintained a separate and distinct evolutionary history that is probably the result of isolation in refugia such as Balkan mountains (Dubašnica, Fruška Gora, Durmitor and Šar Mts). Likewise, intraspecific morphological variability and distinctive phenotypic units were also observed in some species of the Cheilosia (J. Ludoški, unpublished data) and Merodon genera (Francuski et al., 2009;Milankov et al., 2009). Among hoverfly taxa, especially those originating from the Balkans, the presence of divergent units and cryptic diversity is likely to be influenced by an assortment of environmental factors and biogeographical processes (Francuski et al., 2009;Milankov et al., 2009).
In summary, in order to derive a conclusion about the status of divergent phenotypic units it would be necessary to include ecological traits and molecular markers, and more sampling. Bearing in mind that recognition and delimitation of distinctive phenotypic and genetic units are essential in detecting species and genetic diversity, our findings need to be taken into ac-count in taxonomic revisions and ongoing conservation plans. Indeed, determination of species borders as well as intraspecies evolutionarily significant units (ESU, according to Ryder, 1986) is proposed to be crucial for identifying areas of endemism, biodiversity hotspots, and consequently, areas of conservation priority (Myers et al., 2000;Moritz, 2002). Finally, considering that phenotypic diversity presented herein shed light on morphological evolution, we would like to address the use of subtle morphometric traits in the study of intraspecies diversity of other insect species as well.

Phenetic relationships
Analysis of wing geometric morphometric traits in populations from the Balkans has provided phenetic relationships among members of the C. variabilis group. The UPGMA cluster analysis presented a close relationship between the C. melanopa/C. redi species pair indicating that these are recently diverged taxa, while C. lasiopa was the most distant species. The wing differences may indicate phylogenetic patterns but this should be proven by further detailed phylogenetic classification of the C. variabilis group species and evaluation of the phylogenetic significance of wing shape. Results from this study are of interest in their own right, and would, also, be highly useful in understanding the biogeographic history of the region. Finally, observed phenotypic diversity along with ongoing molecular analysis would provide broader insights into the phylogeny of the genus Cheilosia. This is of special interest since phylogenetic relationships among members of the genus Cheilosia was inferred using only two members of the C. variabilis group. In that paper, Ståhls et al. (2004) using cytochrome c oxidase subunit I mitochondrial DNA (COI mtDNA) and nuclear 28S rRNA gene, placed C. variablis and C. melanopa redi as closely related taxa within the Cheilosia s. str. subgenus.