Defining species boundaries in the Merodon avidus complex ( Diptera , Syrphidae ) using integrative taxonomy , with the description of a new species

Several recent studies have detected and described complexes of cryptic and sibling species in the genus Merodon (Diptera, Syrphidae). One representative of these complexes is the Merodon avidus complex that contains four sibling species, which have proven difficult to distinguish using traditional morphological characters. In the present study, we use two geometric morphometric approaches, as well as molecular characters of the 5’-end of the mtDNA COI gene, to delimit sibling taxa. Analyses based on these data were used to strengthen species boundaries within the complex, and to validate the status of a previously-recognized cryptic taxon from Lesvos Island (Greece), here described as Merodon megavidus Vujić & Radenković sp. nov. Geometric morphometric results of both wing and surstylus shape confirm the present classification for three sibling species-M. avidus (Rossi, 1790), M. moenium Wiedemann in Meigen, 1822 and M. ibericus Vujić, 2015-and, importantly, clearly discriminate the newly-described taxon Merodon megavidus sp. nov. In addition to our geometric morphometric results, supporting characters were obtained from molecular analyses of mtDNA COI sequences, which clearly differentiated M. megavidus sp. nov. from the other members of the M. avidus complex. Molecular analyses revealed that the earliest divergence of M. ibericus occurred around 800 ky BP, while the most recent separation happened between M. avidus and M. moenium around 87 ky BP. European Journal of Taxonomy 237: 1–25 (2016) 2


Introduction
The genus Merodon Meigen, 1803 (Diptera: Syrphidae: Merodontini) has become the largest European hoverfly genus due to several recent studies describing many new taxa (Marcos-García et al. 2007;Popov 2010;Radenković et al. 2011;Vujić et al. 2007Vujić et al. , 2012Vujić et al. , 2013aVujić et al. , 2013bVujić et al. , 2015)).A considerable number (37 of 57 species in southeastern Europe) of the taxa are morphologically and/or genetically cryptic, with restricted distributional ranges in particular mountain ranges or islands (Vujić et al. 2016).Larvae of Merodon species are phytophagous on the underground storage organs of plants (Amaryllidaceae, Iridaceae and Hyacinthaceae) (Hurkmans 1993;Andrić et al. 2014).Larvae of Merodon avidus (Rossi, 1790) were recently found in the bulbs of Ornithogalum umbellatum L. (Hyacinthaceae) (Andrić et al. 2014).This host plant is an important foodstuff for adults, as well as for larval development and is widespread in Europe, Southwestern Asia and North Africa.Distribution of this host plant is wider than, but overlaps with, the range of the M. avidus complex (Andrić et al. 2015;Popović et al. 2015).
Taxa of the Merodon avidus complex have been the subject of many studies in the last decade due to perceived taxonomic difficulties (Milankov et al. 2001(Milankov et al. , 2009;;Ståhls et al. 2009).The M. avidus complex is characterized by a considerable morphological variability, especially in the coloration of the antennae, thorax, abdomen and legs (Popović et al. 2015).This colour variability has been explained by the differential availability of trophic resources during the larval stage (Hurkmans 1993).Popović et al. (2015) provided justifications for the identification of three species from the M. avidus complex based on new diagnostic morphological characters, records of the seasonal activity and geographical distribution of bi-voltine adults, nuclear allozyme and mtDNA COI sequence analyses, and descriptions of the ecological preferences of the taxa.However, difficulties in distinguishing the species of this complex based on morphological characters remain, despite the numerous studies on the subject.Spring generations of Merodon avidus are very similar to those of M. moenium Wiedemann in Meigen, 1822 based on external morphology, and, so, are easily confused using existing diagnostic features (e.g., Milankov et al. 2001).Furthermore, it is important to note that species of the M. avidus complex are not distinguishable by traditional visual identification of the structures of male genitalia under a stereo microscope.Additionally, in their analysis of COI barcodes of all Merodon species from Lesvos Island (Greece), Ståhls et al. (2009) revealed the presence of one cryptic taxon (Merodon sp.nova 2, herein described as M. megavidus sp.nov.) within the M. avidus complex, but did not present morphological diagnostic characters to enable characterization of the new taxon.
In the present study, two different geometric morphometric approaches were applied to quantify wing and surstylus shape variability among all hitherto described species of the Merodon avidus complex, including the new fourth taxon M. megavidus sp.nov.Insect wing shape is highly heritable and constitutes an important character for separating species (Birdsall et al. 2000).Geometric morphometric analysis of the wing shape has been successfully used in taxonomic studies of multiple hoverfly taxa (Francuski et al. 2009;Vujić et al. 2013b;Nedeljković et al. 2013Nedeljković et al. , 2015)).The structure and/or shape of parts of the male genitalia are very informative and thus useful for Syrphidae taxonomy and systematics (e.g., Hippa & Ståhls 2005).Differences in male genitalia structure (particularly the surstyli and aedeagal parts) detected in morphological taxonomic studies of hoverflies have revealed them to be an important mechanism of isolation between species (Rotheray & Gilbert 2011).However, differences in genitalia structure between closely-related species can be very small, as shown for many hoverfly genera (Dušek & Laska 1964;Hippa 1990;Nedeljković et al. 2013Nedeljković et al. , 2015;;Vujić et al. 2013bVujić et al. , 2015)).Geometric morphometric analysis of surstylus shape can reveal subtle shape differences that are not detectable or quantifiable by traditional visual examination (Mutanen & Pretorius 2007).This approach has only recently been applied to the taxonomy of the Syrphidae, having been successfully implemented in the genus Chrysotoxum Meigen, 1803 to distinguish four species (Nedeljković et al. 2013(Nedeljković et al. , 2015)).The male genitalia of Merodon species are sclerotized and rigid and thus stable structures, so they are well suited for geometric morphometric analysis.
Our study had three objectives: (1) to further clarify the species borders of all taxa within the M. avidus complex using integrative taxonomy (geometric morphometrics of wings and male surstylus shape and mtDNA COI sequences); (2) to estimate times of divergence among investigated taxa; and (3) to provide descriptions and diagnostic characters of the new species.A typological (or morphological) species concept is applied in this study, integrating all available data.

Studied material
A total of 444 specimens belonging to the Merodon avidus complex from Bulgaria, Croatia, Greece, Italy, Montenegro, Morocco, Serbia, Spain and Turkey was analysed (Fig. 1; Appendix 1).Specimens of M. avidus were sampled during spring, summer and autumn.All specimens were identified by Ante Vujić and Snežana Radenković, and labelled using unique codes.For the description of the new species, the terminology follows McAlpine (1981) for non-genitalic morphology and Marcos-García et al. (2007) for morphology of the male terminalia.

Wing morphometry
Wings of all 444 available specimens were analysed using a landmark-based geometric morphometric approach (Appendix 1).Population analysis was carried out for 418 specimens of 23 populations (Appendix 1, marked with *).Populations with small sample sizes were excluded from the analysis to avoid statistical errors.The left wing of each specimen was removed using micro-scissors and mounted in Hoyer's medium on a microscope slide.Eleven homologous landmarks were digitized at vein intersections or terminations that could be reliably identified and best represented wing shape using TpsDig 2.05 software (Rohlf 2006).Generalized least-squares Procrustes superimposition was first applied to the landmark data to remove non-shape variation in terms of location, scale and orientation, and also to superimpose the objects in a common coordinate system (Rohlf & Slice 1990;Zelditch et al. 2004).For the wing shape analysis, partial warp scores were calculated (Zelditch et al. 2004).Procrustes superimposition and partial warps were computed using the free IMP software CoordGen7.14and CVAgen 7.14a (Sheets 2012).MorphoJ v.2.0.software (Klingenberg 2011) was used for visualization of thin-plate spline deformation.Surstylus morphometry Shape analysis of the posterior part of the left surstylus was carried out on 125 specimens of the M. avidus complex using a semi-landmark geometric morphometric approach (Appendix 2).The posterior part of the left surstylus (Fig. 2A: psl) was removed using a scalpel and placed on its side in glycerol on a microscope slide, and a coverslip was placed on top of the surstylus to immobilize it.Due to a lack of homologous anatomical loci, 30 semi-landmarks were digitized using the 'resample curve by length' option in TpsDig 2.05 software (Rohlf 2006).CoordGen 7.14 with an integrated Semiland module was used for semi-landmark superimposition using a distance-minimizing protocol, which minimizes the shape differences due to the arbitrary nature of semi-landmark positions along the curve (Bookstein 1997;Zelditch et al. 2004).

Statistical analyses
Principal component analysis (PCA) was used to test variability of wing and surstylus shape without a priori defining groups.Multivariate analysis of variance (MANOVA) and Fisher LSD post-hoc tests were used to test the variability connected with shape differences between species.Additionally, canonical variates analysis (CVA) and discriminant function analysis (DA) were used to test significance in wing and surstylus shape differences, to produce distance matrices and to graphically present results.The phenetic relationships among taxa were determined by UPGMA analysis based on squared Mahalanobis distances computed from the DA applied to wing variables.All statistical analyses were conducted using Statistica for Windows (Statsoft Inc. 2015: version 12).

Molecular analysis
The 5'-end of the cytochrome c oxidase 1 (COI) gene ('barcode' region) was amplified for three specimens of M. megavidus sp.nov.from Lesvos Island, Greece.DNA was extracted from a leg of each specimen using a slightly modified SDS Extraction Protocol (Chen et al. 2010).Genomic DNA vouchers are conserved at the Faculty of Sciences, Department of Biology and Ecology, University of Novi Sad (AU827, AU885, AU886) and the Finnish Museum of Natural History, Helsinki (MZH S532).DNA barcodes were amplified with the forward primer LCO (5'-GCTCAACAAATCATAAAGATATTGG-3') and the reverse primer HCO (5'-TAAACTTCAGGGTGACCAAAAAATCA-3') (Folmer et al. 1994).PCR amplifications were carried out in 20 μl reaction volumes using the following mix of components: 1 × PCR buffer (Thermo Scientific), 2.5 mM MgCl 2 , 0.1 mM of each nucleotide, 1 U Taq polymerase (Thermo Scientific), 2 pmol of each primer and 50 ng template DNA.PCR conditions were: initial denaturation for 2 min at 95°C; 30 s denaturation at 94°, 30 s annealing at 49°C, 2 min extension at 72°C/30 cycles; and the final extension for 8 min at 72°C.Obtained amplicons were purified using the Exo-Sap purification protocol (Thermo Scientific).Forward sequencing using the PCR primer was conducted at the Sequencing Service Laboratory of the Finnish Institute for Molecular Medicine, Helsinki, Finland (FIMM).In order to confirm the taxonomic status of M. megavidus sp.nov.specimens, we created a dataset consisting of 22 DNA barcode sequences, 18 of which represented GenBank-accessioned sequences of M. avidus, M. moenium and M. ibericus Vujić, 2015 (6 sequences per species, accession numbers indicated in Fig. 6, and one sequence for M. megavidus sp.nov.(generated for the Ståhls et al. 2009 study), in addition to three newly-generated M. megavidus sp.nov.DNA barcode sequences in this study).Sequences of Eumerus sulcitibius Rondani, 1868 (Acc.No. KT157875) and Platynochaetus setosus (Fabricius, 1794) (Acc.No. KM224512) (both Merodontini) were used as outgroups.Sequence alignment was conducted using the Clustal W algorithm (Thompson et al. 1994), as implemented in BioEdit (Hall 1999), and final modifications were done by hand.The total length of the dataset after alignment and trimming was 634 nucleotides.We constructed a Neighbour-joining (NJ) tree, and phylogenetic trees using Maximum Likelihood (ML) and Maximum Parsimony (MP) analyses.The NJ and ML trees were constructed using MEGA version 6 (Tamura et al. 2013) using a Tamura-3 parameter model of nucleotide substitution with Gamma correction for among-site variation in substitution rates (estimated in the same software).MP analysis was done using NONA (Goloboff 1999), with the aid of Winclada (Nixon 2002), using the heuristic search algorithm and 1000 random addition replicates (mult x 1000); holding 100 trees per round (hold ⁄ 100), max trees set to 100 000; and applying TBR branch swapping.Bootstrap values were calculated for each tree using 1000 replicates.The genetic relationships among species were also tested using a Median-Joining (MJ) network (Bandelt et al. 1999) generated with Network v4.6.1.3(available from http://www.fluxus-engineering.com/sharenet.htm) by applying the default settings (ɛ = 0 and the variable sites weighted equally = 10), with additional postprocessing with the maximum parsimony (MP) option.
Pairwise Φst values among species of the Merodon avidus complex were calculated, together with an exact test of population differentiation, using ARLEQUIN 3.5.1.2(Excoffier & Lischer 2010).Furthermore, we calculated the genetic distance between defined species using MEGA version 6 (Tamura et al. 2013) in order to estimate times of divergence between genetic clusters based on uncorrected p-distances divided by the pairwise evolutionary rate/MYR as described in Pröhl et al. (2010).We used the pairwise sequence divergence of the COI gene, relative to the mutational rate of 2.3% per million years, as estimated for various arthropod taxa (Brower 1994).

Diagnosis
Medium-to large-sized species (13-18mm); black mesoscutum with four white microtrichose longitudinal stripes; tapering orange and black abdomen with white, transverse, microtrichose bands on tergites 2-4 (exceptionally without bands on tergite 2); tarsi reddish-orange dorsally; hind femur medium wide and slightly curved (Fig. 3C-D), with very short pile posteroventrally.Merodon megavidus Vujić & Radenković sp.nov.belongs to the avidus complex (male genitalia in all species identical in shape, as on Fig. 2).Merodon megavidus sp.nov.can be separated from the other members of the complex by larger size, golden body pile, bright orange colour of the pale parts of legs and extremely short pile on hind femur (Fig. 3C-D).These characteristics contrast with other species from the complex, which have yellow to grayish pale body pile and longer pile on the hind femur (Fig. 3A-B).

Etymology
The name megavidus refers to the large size (Greek word megas means "large") and great similarity with Merodon avidus.

Description
Male (Figs 2C, 3, 4A, 5A) Head (Fig. 4A).Antenna (Fig. 4A) orange, first flagellomere 1.8-2.0times as long as wide, 2.0 times longer than pedicel, concave, apex acute; arista: second, third and basal part of fourth flagellomeres pale, fourth flagellomere dark brown in apical ⅔ and thickened basally, 1.4 times longer than first flagellomere; with short, dense microtrichia.Face and frons black, covered with long golden pile and silver, dense microtrichia.Oral margin shiny black, except for the lateral microtrichose areas (Fig. 4A).Vertical triangle isosceles, shiny black except in front of the anterior ocellus that has pale microtrichia, covered with long orange pile except for black pile on the ocellar triangle.Ocellar triangle equilateral.Eye contiguity about 12 ommatidia long.Vertical triangle: eye contiguity: ocellar triangle = 1.5 : 0.7 : 1. Eye pile dense, white.Occiput with orange pile, along the eye margin with dense white microtrichia and posteriorly with metallic, bluish-greenish lustre.
THorax.Mesoscutum and scutellum black with bronze lustre, covered with relatively long, dense, erect golden pile.Side of mesoscutum above wing-base with a patch of black pile.Mesoscutum with two lateral and two submedian, longitudinal, white microtrichose stripes.Proepimeron, posterior anepisternum, anteroventral and posterodorsal part of katepisternum, anepimeron, metasternum and katatergite with long golden pile and grey-green microtrichia.Wing hyaline, with dense microtrichia; veins dark brown except for light brown C, Sc and R1.Calypter pale yellow.Haltere with light brown pedicel and yellow capitulum.Legs orange, except for the black basal ¾ of the front-and mid-femora.Pile on legs golden.Hind femur (Fig. 3C) moderately thickened and curved, about 3.6 times as long as deep.Pile on hind femur very short.abdomen (Fig. 5A).Dark with white microtrichose bands, tapering, 1.4 times longer than mesonotum (including scutellum).Tergites orange and reddish except for black tergite 1 and central parts of tergites 2-3 (and 4) (Fig. 5A); orange-reddish parts of variable size on tergite 3 and 4, laterally and along microtrichose bands.Tergites 2-4 each with a pair of white microtrichose marks (exceptionally absent only on tergite 2); tergites 3-4 with wide, oblique bands (Fig. 5A).Pile on tergites golden.Sternites translucent, orange to brown towards the tip of the abdomen, covered with long yellow pile.Male genitalia (Fig. 2).Similar to all species of the M. avidus complex.Anterior lobe of surstylus broad and hairy (Fig. 2A-B); posterior lobe of surstylus ellipsoidal at ventral margin (Fig. 2A-B); cercus rectangular, without prominences (Fig. 2A).Hypandrium elongate and sickle-shaped, without lateral projections (Fig. 2C); lingula long (Fig. 2C).2D, 4B, 5B) Similar to the male except for typical sexual dimorphism and for the following characteristics: first flagellomere broader and longer; frons with two wide (about 0.34 width of frons) lateral silver microtrichose longitudinal stripes; frons in the widest part about 0.25 width of head; white microtrichose longitudinal stripes on mesoscutum more visible; broad stripe of black pile between wing bases; tergites predominately red except for tergite 1 and darkened parts of tergites 2-4 in some specimens (Fig. 5B); white, microtrichose, transverse bands on tergites 3-4 (Fig. 5B); tergites 2-3 with black pile on dark parts; white microtrichose bands solely with pale pile.obtained from two specimens morphologically similar to M. avidus taken from Lesvos Island.In the key prepared as supporting material for Ståhls et al. (2009), M. avidus and M. sp.nova 2 key out together, without any morphological differences being described.

Distribution and habitat data
Lesvos Island (Greece).Maquis shrubland.

Molecular data
Merodon megavidus sp.nov.clearly differs from other members of M. avidus complex based on our barcoding fragment of COI.All conducted phylogenetic analyses resulted in similar tree topologies (Figs 6-7; Appendix 3).Sequences from the Merodon avidus complex formed three separate clusters: one cluster represented M. ibericus, a second comprised all M. megavidus sp.nov.sequences with highly-significant bootstrap values (100), the third grouped sequences of M. avidus and M. moenium together.All four M. megavidus sp.nov.sequences were identical, defining one haplotype unique to the species.The number of mutational steps between M. megavidus sp.nov.and M. avidus is 14, with 13 and 22 mutational steps between M. megavidus sp.nov.and M. moenium and M. ibericus, respectively (Fig. 8).
Our UPGMA tree based on genetic distances among species revealed genetic relationships among four taxa of the Merodon avidus complex that support and strengthen the branch positions in wing and surstylus phenograms described below (see Fig. 9).
Significant pairwise genetic divergence (ϕst value) was detected in each pairwise comparison between M. megavidus sp.nov.and each of the other species M. ibericus, M. avidus and M. moenium (0.420, 0.492 and 0.529, respectively).Sequence divergence (uncorrected p distance) of the COI gene was used to assess relative divergence times between the four Merodon taxa (Table 1), indicating an initial separation of M. ibericus from the rest of the complex around 800 ky BP.Divergence between M. megavidus sp.nov.and M. avidus/M.moenium occurred around 500 ky BP.The most recent separation happened between M. avidus and M. moenium, around 87 ky BP.
Additionally, the phenogram constructed based on squared Mahalanobis distances of wing shape showed that all 26 analysed populations grouped according to species (Fig. 13).

Discussion
Differentiation and taxonomic status of the cryptic species of the Merodon avidus complex has been controversial for many years.An initial molecular analysis showed that M. avidus is a genetically and geographically structured taxon with at least two cryptic species, which were designated as M. avidus A and M. avidus B (Milankov et al. 2001).Further investigations based on the mtDNA COI marker expanded on this by finding two additional cryptic species from the Iberian Peninsula (M.bicolor Gil Collado, 1930, now M. ibericus) and Lesvos Island (Greece) (M.sp.nova 2, now M. megavidus sp.nov.) (Milankov et al. 2009, Ståhls et al. 2009).In Milankov et al. (2009), wing shape analysis could not distinguish taxa that could be differentiated by allozyme loci, even though these taxa differed in wing size.Also, they found a great deal of similarity between allopatric metapopulation pairs of M. avidus Fig. 12. Thin-plate spline deformation grids showing wing shape differences between analysed species.Differences between the species have been exaggerated five-fold to make them more visible.
A and M. avidus B; populations of M. avidus A from Macedonia and the Pannonian plain overlapped with M. avidus B populations from Durmitor, Stara Planina and Kopaonik in terms of wing shape.Only a population of M. avidus A from Morinj (Montenegro) was clearly separated from other populations.Subsequent analysis showed that this substantial overlap in wing shape was a result of incorrect species identifications due to the morphological similarities between spring generations of M. avidus and M. moenium (Popović et al. 2015).Previously, it had been considered that M. moenium (M.avidus B) was a species distributed on mountainous regions of the Balkan Peninsula.However, it is now known that while M. moenium is predominantly distributed in the continental parts of Europe, it also occurs in some parts of the Mediterranean basin and the Black Sea coast, where it can occur sympatrically with M. avidus.Further, allozyme analysis of morphologically different spring and autumn generations from Umag (Croatia) clarified this misclassification issue by showing that both morphotypes belonged to M. avidus (Popović et al. 2015).
Our molecular analysis using the barcoding COI fragment for an additional three specimens from Lesvos Island supports the presence of a morphologically cryptic taxon (M.megavidus sp.nov.), previously revealed by Ståhls et al. (2009).All our phylogenetic analyses placed M. megavidus sp.nov.as a separate independent cluster from M. ibericus and the M. avidus/M.moenium branch (Figs 6-9).
According to estimated divergence times, all diversification occurred in the Pleistocene (2.6 to 0.0117 MYA).This geological epoch was marked by repeated (at least 20) glacial and interglacial periods, which influenced speciation and the distributions of many recent taxa in Europe (Julius & Kukla 1977;Martinson et al. 1987;Perissoratis & Conispoliatis 2003).Recent studies have revealed many examples of insect species that have altered their ranges and/or evolved as a response to repeated isolation during glacial-interglacial cycles across three large peninsulas of Southern Europe (Hewitt 2001;Konstantinov et al. 2009;Dapporto 2010;Nicholls et al. 2010;Zhu et al. 2013).
Fig. 13.UPGMA phenogram constructed using squared Mahalanobis distances of wing shape for populations of species of the M. avidus complex.
The first mitochondrial diversification in the M. avidus complex took place in the Calabrian stage of the Early Pleistocene (around 800 ky BP), when M. ibericus diverged from a common ancestor.The Günz-Mindel interglacial corresponds to the approximate period when separation of M. megavidus sp.nov.from the M. avidus/M.moenium lineage occurred.The most recent diversification, i.e., between M. avidus and M. moenium, most likely took place at the end of the Riss-Würm interglacial or the beginning of the Würm glaciation period.The low species resolution ability of the 3' and 5' fragments of COI to separate the M. avidus and M. moenium lineages (Milankov et al. 2009;Popović et al. 2014Popović et al. , 2015) ) can be explained by recent speciation and by the founder effect associated with postglacial recolonization of northern Europe (Shikano et al. 2010).Both species share an identical haplotype of the COI 3'-end (Popović et al. 2014), although the 5'-end fragment haplotypes of DNA barcodes are not shared, they give low support in cluster analyses (Popović et al. 2015).
The most probable scenario for postglacial colonization of Europe is expansion of the M. avidus and M. moenium lineages from the Balkan Peninsula into Central Europe (a "grasshopper" pattern of colonization according to Hewitt 1999).The Pyrenees acted as a significant geographical barrier preventing the dispersal of M. ibericus to other European areas.
Our geometric morphometric analyses show that all investigated species of the M. avidus complex can be successfully discriminated based on wing and surstylus shape.The high overall success rate of classification indicates that wings, and especially surstylus shape, have meaningful interspecific discriminatory power.It is important to emphasize that, contrary to previous geometric morphometric studies, we found highly significant differences in wing shape between M. avidus and M. moenium.The phenogram based on squared Mahalanobis distances for wings is congruent with the UPGMA tree based on genetic distances among species.Wing shape, as a highly heritable structure, has greater importance in insect taxonomy than wing size (Birdsall et al. 2000), which has been confirmed by earlier studies of hoverflies in which wing shape has been successfully used for identification and delimitation of species (Vujić et al. 2013b;Nedeljković et al. 2013Nedeljković et al. , 2015)).
Additional evidence for species distinctiveness is provided by the highly significant differences in surstylus shape between M. avidus, M. ibericus, M. megavidus sp.nov.and M. moenium.The main differences in surstylus shape are connected to the posterior margin of the posterior part of the surstylus lobe, which is involved in gripping the female during copulation.It could be hypothesized that these shape differences contribute to a mechanism of sexual isolation, especially since M. avidus has a distinct surstylus shape compared to partly sympatric M. moenium.Although the exact mechanism of sexual isolation in Syrphidae is not known, traditionally the shape of male genitalia is deemed a significant mechanism of isolation between species (Rotheray & Gilbert 2011).Given that the morphology of male genital structures is considered one of the fastest evolving traits in animal groups with internal fertilization (Soto et al. 2013), and that insect genitalia are conspicuously variable even in closely related taxa that are otherwise morphologically very similar (Hosken & Stockley 2004), we assert that significant differences in surstylus shape is strong evidence for species delimitation.Recent studies have found that natural and sexual selection (and their interaction) caused insect genital evolution (Hasson et al. 2009;House et al. 2013), and that the degree of morphological differences in the structures of genitalia is associated with geographic distance (Soto et al. 2013).Despite their recent speciation and geographic proximity, M. avidus does not exhibit a similar surstylus shape to M. moenium.Other studies have found that other factors, such as specific host plant relationships Soto (2012) or sexual selection (Hosken & Stockley 2004), may promote further evolution in genital morphology.
Our population-level morphometric analysis of wing shape is in concordance with species delimitation findings.All populations clustered according to species.Furthermore, in this analysis, M. megavidus sp.nov.from Lesvos Island is distinct from all other species and/or populations, and it is important to underline that this taxon was clearly differentiated from population(s) of M. avidus also from Lesvos Island.

Fig. 1 .
Fig. 1.Map of population sampling locations of the Merodon avidus complex from the Western Palaearctic.

Fig. 7 .
Fig. 7. Maximum Likelihood tree based on a 5' fragment of COI mtDNA sequences from the Merodon avidus complex.Bootstrap values (1000 replicates) are shown next to the branches.The tree is drawn to scale, with branch lengths proportional to the number of substitutions per site.

Fig. 8 .
Fig. 8. Median-joining network of the mtDNA 5'-end of the COI gene.Circle sizes are proportional to haplotype frequencies.Each branch represents one mutational step; if more than one mutational step is present, it is denoted by the given number.

Fig. 9 .
Fig. 9. UPGMA tree based on pairwise genetic distances for four species from the Merodon avidus complex.

Fig. 10 .
Fig. 10.Differences in wing shape among species of the M. avidus complex.A. Scatter plot of individual scores of CV1 vs CV2.B. Scatter plot of individual scores of CV2 vs CV3.

Fig. 11 .
Fig. 11.UPGMA phenogram constructed using the squared Mahalanobis distances of wing shape for species of the M. avidus complex.

Fig. 14 .
Fig. 14.Differences in the posterior part of the surstylus among species of the M. avidus complex.A. Scatter plot of individual scores of CV1 vs CV2.B. Scatter plot of individual scores of CV2 vs CV3.

Fig. 15 .
Fig. 15.Differences in the posterior part of the surstylus among species of the M. avidus complex.A. UPGMA phenogram constructed using squared Mahalanobis distances.B. Thin-plate spline deformation grids showing overall shape differences between analysed species.

Table 1 .
Below diagonal -pairwise p distances; above diagonal -estimated divergence times in years.
List of specimens used for wing geometric morphometric analysis, by geographical area and species.Appendix 2. List of specimens used for surstylus geometric morphometric analysis, by geographical area and species.