Sibling species delimitation and nomenclature of the Merodon avidus complex ( diptera : Syrphidae )

A study of the relationships between 21 southern European, Moroccan and Turkish populations of the Merodon avidus species complex was carried out. Based on a parallel study of type material from several museums, documented diagnostic morphological characters, season of adult activity and geographical distribution, we justify the use of the following names for three closely related taxa in this complex: M. avidus (Rossi, 1790), M. moenium (Wiedemann, 1822), and M. ibericus Vujić nom. n. (new name for M. bicolor Gil Collado, 1930, preoccupied by M. bicolor Walker, 1852). Neotypes of Merodon avidus, M. spinipes (Fabricius, 1794) and M. quadrilineatus Lioy, 1864 are designated. Lectotype of Merodon moenium (Wiedemann in Meigen, 1822) is designated here. A cluster analysis of DNA barcoding sequences clearly separated M. ibericus, but not M. avidus and M. moenium, even though the lack of shared haplotypes, analysis of molecular variance (AMOVA), pairwise Φst values together with allozyme and ecological niche analyses revealed statistically significant percentage of variation among all three species in the Merodon avidus complex. Analysis of 5 diagnostic enzyme loci revealed the presence of genetic differentiation among the M. avidus/moenium complex populations investigated (Fst = 0.654) and species-specific alleles were found at the AAT locus. The presence of two separate related taxa within the M. avidus/ moenium complex was further supported by an UPGMA tree based on Nei’s (1978) genetic distances. The value of Nei’s measure of genetic identity (I = 0.520) between two large (meta) populations of M. avidus and M. moenium suggest that these taxa are sibling species. Populations from Djerdap (Serbia) confirmed the presence of temporal divergence between these species at a locality where they occur sympatrically, while spring and autumn populations from Umag (Croatia) provide an example of morphological plasticity within the species M. avidus. Ecological niche analysis contributed to the species delimitation. Review of the available genetic and ecological data confirmed our hypothesis that the M. avidus species complex, in addition to M. ibericus Vujić nom. n. from the Iberian Peninsula, consists of two sibling species in the rest of Europe and indicated their recent speciation. * Corresponding author.


IntroductIon
Merodon Meigen, 1803 is the largest European genus of the family Syrphidae (Vujić et al., 2012).Larvae of all the species of Merodon develop in bulbaceous plants, which are most diverse in the Mediterranean area, as is the genus Merodon (Radenković et al., 2011).This high diversity, together with the presence of cryptic species, taxonomic ambiguities and unclear relationships between taxa (Milankov et al., 2008;Marcos-García et al., 2011;Radenković et al., 2011;Vujić et al., 2012) makes the study of phylogenetic relationships within the genus a challenging undertaking.
As recognized until very recently Merodon avidus (Rossi, 1790) is a widespread European taxon, only absent in northern Europe (with the exception of southern Sweden), but also present in North Africa, the Middle East and Asia Minor (Dirickx, 1994).Great variation in the coloration of the antennae, thorax, legs and abdomen has made it difficult to correctly identify specimens (Hurkmans, 1993).M. avidus was described six times under different names: spin ipes (Fabricius, 1794), moenium (Wiedemann in Meigen, 1822), rufitibius (Rondani, 1845), graecus (Walker, 1852), quadrilineatus (Lioy, 1864).Hurkmans (1993) pointed out that there are no distinguishable geographic trends in morphology or colour of the various local forms.The intrapopulation variability in size and colour is explained as a consequence of the quantity and quality of food ingested during the larval stages.Almost a decade later, morphological analysis of tergites II and III, tibiae and mesoscutum, together with an allozyme variability analysis indicated that M. avidus is actually a geographically and genetically structured taxon, within which (Milankov et al., 2001) two cryptic species could be recognized.These were designated as Merodon avidus A (Mediterranean region and central parts of the Balkan Peninsula) and M. avidus B (mountainous regions on the Balkan Peninsula).
Subsequently, the molecular and phenotypic diversity of Merodon avidus from the Balkan Peninsula based on sequences of the 3´ end of the mitochondrial cytochrome c oxidase subunit I gene (COI) and wing parameters were analyzed (Milankov et al., 2009).Wing shape and the mtDNA COI marker did not reveal a clear delineation between M. avidus A and M. avidus B, unlike the diagnostic enzyme loci (Milankov et al., 2001).In spite of the high haplotype diversity within the M. avidus species complex and presence of unique haplotypes within each genetic The aim of the present study was to obtain a new insight into the relationships between populations within the Merodon avidus complex, by increasing the number of techniques used in the analysis, by considerably increasing the number of geographical areas sampled and by sampling generations of M. avidus A in different seasons at the same locality.The occurrence of more than one generation per year is not common among Merodon taxa.Usually they have only one generation per annum, as in M. avidus B.
Our principal objectives were to: (a) quantify the genetic variability of diverse natural populations of the putative sibling species in the Merodon avidus complex, (b) reveal boundaries of potential cryptic species, (c) examine genetic differentiation between these cryptic taxa and (d) evaluate and resolve their current taxonomic status.In addition, we examined a series of environmental factors that could drive the perceived genetic differentiation.

MAterIAl And MethodS nomenclatural study
Examination of type material was used to establish the nomenclatural status of the following European taxa: Merodon avidus (Rossi, 1790) (M.avidus A sensu Milankov et al., 2001), Mero don bicolor Gil Collado, 1930 and Merodon moenium (Wiedemann, 1822) (M.avidus B sensu Milankov et al., 2001), and their synonyms.The main depositories for types belonging to the Merodon avidus complex are: Natural History Museum, London, UK (BMNH); University of Novi Sad, Department of Biology and Ecology, Serbia (FSUNS); Museo Zoologico della Specola, Florence, Italy (LFC); Museo Nacional de Ciencias Naturales, Madrid, Spain (MNCN); Museum für Naturkunde, Berlin, Germany (ZMB) and Natural History Museum of Denmark, Copenhagen, Denmark (ZMUC).During the study of the collections in European museums, a thorough search was carried out to locate all true syntype specimens and designated lectotypes of taxa of the Merodon avidus complex.In the case of taxa for which no type material could be located, a neotype was designated as appropriate.

dnA barcoding analysis
Specimens for molecular analyses were collected in the field, from different parts of Central and Southern Europe, Morocco and Turkey.Adults of the Merodon avidus complex from Central and Southern Europe and Turkey were initially identified as M. avidus and M. moenium.In total this material consisted of 118 specimens from 19 populations.Six of the localities from which material was collected were in Serbia, five in Italy, one in Croatia and two each in Montenegro, Greece and Turkey (Fig. 1).The collected material is detailed in the Supporting Information Table S1.Specimens were frozen immediately on collection and stored at -20°C until the analysis.
For the DNA barcoding analyses, we selected one specimen of each determined species, M. avidus and M. moenium, from each locality.In addition, we selected specimens of M. ibericus from different localities and included them in the analysis (Supplemental Material Table S2).
Total DNA was extracted from legs and abdomens of the specimens using the SDS Extraction Protocol (Chen et al., 2010) with slight modifications.The cytochrome c oxidase subunit (COI) of mitochondrial DNA was amplified with forward primer LCO (5'-GCTCAACAAATCATAAAGATATTGG-3') and reverse primer HCO (5'-TAAACTTCAGGGTGACCAAAAAATCA-3') (Folmer et al., 1994).PCR reactions were carried out in 20 μl re-unit, M. avidus A and M. avidus B shared four haplotypes.Also, relatively high intraspecific divergence and phenotypic substructuring, detected within M. avidus A, raised the question of possible genetic subdivision between these populations (Milankov et al., 2004(Milankov et al., , 2009)).On the other hand, COI-3' fragment revealed differences between some morphologically indistinguishable taxa in the M. avidus complex, separating the Iberian Peninsula variety initially identified as M. avidus B sensu (Milankov et al., 2001) from the rest of the complex, which was then redefined as a valid species -M.bicolor Gil Collado, 1930(Milankov et al., 2009).The nomenclatural problem with using the name Merodon bicolor is discussed here.
Over the last decade integrative taxonomy was used for resolving complex taxonomical questions (Dayrat, 2005).Combining and evaluating different kinds of data (comparative morphology, population genetics, ecology, etc.) is the most reliable method for delimiting taxa.Integrative taxonomy has proved a useful approach for delimitating various taxa in Syrphidae (e.g.Ståhls et al., 2008Ståhls et al., , 2009;;Vujić et al., 2012;Nedeljković et al., 2013).In the case of the M. avidus complex, we previously highlighted the importance of using many characters (allozyme variability analyses, mitochondrial sequence data and quantitative traits) when delineating cryptic taxa (Milankov et al., 2009).Since in previous research a limited number of specimens was available, we focused on the M. avidus complex in a wide geographic area and analyzed a larger number of specimens.The analysis of the 3' fragment of COI gene in this study failed to discriminate between M. avidus A and M. avidus B, but the allozyme loci revealed the presence of two cryptic taxa (Popović et al., 2014).
Over several consecutive years we noticed that Merodon avidus A is morphologically variable, with the spring generation morphologically more similar to M. avidus B than the autumn generation.This led us to the assumption that environmental conditions could contribute to the seasonal predominance of the different phenotypes and could also be connected with the genetic separation within the M. avi dus complex.Analysis of the environmental preferences of taxa is widely used to investigate the environmental limits of taxa, the environmental conditions required by protected or invasive species, mechanisms of speciation and species delimitation (Stockman & Bond, 2007;Bond & Stockman, 2008;Nosil et al., 2009;Medley, 2010;Wooten & Gibbs, 2012;Zhou et al., 2012;Zhu et al., 2013;Petersen, 2013).Various authors have pointed out that differences in environmental preferences can be used as additional evidence for delimiting species (e.g.Rissler & Apodaca, 2007;Sánchez-Fernández et al., 2011).
The immature stages and host plants of the majority of species of Merodon are still unknown.The larva of M. avi dus A was recently found in the bulbs of Ornithogalum L. (Hyacinthaceae) and described by Andrić et al. (2014).The genus Ornithogalum is well distributed in the Mediterranean basin, almost all of Europe (except the North) and South and Central Africa, with a much greater range than the M. avidus complex (Popov, 2001).action volumes and the reaction mix consisted of 1 × PCR buffer, 2.5 mM MgCl 2 , 0.1 mM of each nucleotide, 1 U Taq polymerase, 2 pmol of each primer and 50 ng DNA.Thermocycler 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.The PCR products were purified using the Exo-Sap purification method, according to the manufacturer's protocol (Thermo Scientific).The fragments were sequenced using the same forward primer and all sequences were submitted to GenBank as COI barcodes (for accession numbers see Table S2).
The COI barcode sequences were clustered using Neighbourjoining (NJ), Maximum-likelihood (ML) and Maximum Parsimony (MP) analyses.Alignment of COI barcode sequences was performed using Clustal W algorithm (Thompson et al., 1994) as implemented in BioEdit version 7.2.5.(Hall, 1999).As an outgroup we used Cheilosia scutellata (GenBank Accession number AY533354).The final length of the sequences in the dataset after alignment and trimming was 631 nucleotides.DNA polymorphisms (number of polymorphic sites, h -haplotype diversity, π -nucleotide diversity, k -mean number of pairwise differences), parameters of overall genetic variability and haplotype frequencies were calculated using DnaSP v5 (Librado & Rozas, 2009).The NJ and ML trees were calculated using MEGA version 6 (Tamura et al., 2013), using the Tamura 3-parameter model with a gamma distribution of 0.05 (as it was determined using same software) and 1000 replicates.The bootstrap values were calculated for each tree using 1000 replicates.All base positions were treated as equally weighted characters.A maximum parsimony (MP) tree was constructed using NONA (Goloboff, 1999) and spawn with the aid of Winclada (Nixon, 2002), using heuristic search algorithm with 1000 random addition replicates (mult × 1000); holding 100 trees per round (hold / 100), max trees set to 100 000; and applying TBR branch swapping.Analysis of molecular variance (AMOVA) among and within the species analysed and pairwise Φ st values between three species were calculated using Arlequin 3.5.1.2 (Excoffier & Lischer, 2010).

Allozyme analysis
Since the initial DNA barcoding analysis clearly separated M. ibericus from the M. avidus / moenium pair, the latter was subjected to an allozyme variability analysis.In total 118 specimens from 19 populations (Fig. 1; Supporting Information Table S1) were included in this analysis.
In order to compare results for different populations and taxa, a mix of samples from different populations was loaded onto the same gel before the electrophoresis.Depending on metabolic function and regional distribution of enzymes, different body regions were used for analysis of enzyme variability (head: AAT, GPI, ME; thorax: GPD, IDH).Alleles were marked from "a" to "c" depending on the distance they moved, where "a" was the allele that moved the shortest.This was not in accordance with a previous study of these enzymes (Milankov et al., 2001), where different kinds of labels were used.Genotype and allele frequencies were calculated based on the banding patterns and genetic interpretation of zymograms.Statistical analysis of allozyme data was performed using the computer programme POPGENE version 1.32 (Yeh et al., 1999).The following parameters were The programme package Statistica 10 for Windows (StatSoft, 2011) was used to test significance of differences in allele frequencies between the morphologically defined species M. avidus and M. moenium.
The correlation between genetic distances and geographical distances was performed using Mantel's non-parametric test on pairwise distance matrices (Mantel, 1967), using the MANTEL procedure implemented in IBD version 3.23 (Jensen et al., 2005).Only accurate distributional data was used.Localities with geographic coordinates were used without modification.For records with only locality names coordinates were assigned using Google Earth (Google Inc., 2013), which was also used to confirm the validity of the localities with coordinates.All localities were represented in DivaGis (v7.5) (Hijmans et al., 2012).Nineteen bioclim variables and altitudes (2.5 arc-minutes resolution) were generated for each locality on the basis of the WorldClim data set (Hijmans et al., 2005).Principal component analysis (PCA) was used to summarize multidimensional environmental variables and compare the ecological niches of the species investigated.PCA was carried out using a normal varimax rotation of factor loadings.Only factors with an eigenvalue greater than 1 and that contribute more than 5% of total variance were considered significant.Variables with a factor loading greater than 0.8 were interpreted as significantly correlated with the correspondent factor.ANOVA and a Tukey HSD post hoc test were carried out to determine if the derived factor scores between the two spe-cies are significantly different.A scatter plot of PCA score values was used to graphically display the position of these species in environmental space.Statistical analyses were performed using Statistica 10 for Windows (StatSoft, 2011).

taxonomy of the Merodon avidus complex
Diagnostic morphological characters for separating Merodon avidus (M. avidus A in Milankov et al., 2001) and M. moenium (M. avidus B in Milankov et al., 2001) are defined and summarized in Table 1 of Milankov et al. (2001: page 17).They separate M. moenium from typical specimens of mostly the late summer / autumn generation of M. avidus.
M. avidus has a pair of whitish, microtrichose spots on tergite 2 and broad microtrichose bands on tergites 3 and 4 (Figs 2A, B, 3A); usually pale tibiae (Figs 5A, 6A) and slightly shorter body hairs, especially on the apical tergites (Fig. 4A).In M. moenium tergite 2 is undusted and shiny; there are narrow microtrichose stripes on tergites 3 and 4 (Figs 2C, 3C); the tibiae are always partly dark (Figs 5C, 6C) and the body hairs longer (Fig. 4B).Both species have orange lateral triangular spots on tergite 3, but in M. avidus the anterior part of tergite 3 is also predominately orangered, except medially, where it is narrowly black (in darker specimens at least small orange areas are present anterosublaterally) (Figs 2A, B, 3A), while in M. moenium tergite 3 is black (in the female it is anterolaterally orange-red, but with a black posterior margin, in contrast with M. avidus where the posterior margins of tergites 2-4 in both sexes are paler) (Figs 2C, 3C).
These morphological characters are very stable in populations of Merodon moenium, but can vary in M. avidus populations, with adults emerging in spring differing in appearance from adults of the late summer / autumn generations.Some early spring specimens of M. avidus (usually the first half of May) can hardly be distinguished from slightly paler specimens of M. moenium, in having darkened tibiae (Figs 5B, 6B), darkened tergite 3 (Fig. 2B) and less microtrichose tergites (Fig. 2B).In these cases we used geographical data (M.moenium was never found in the Eumediterrranean) and season (adults of M. moenium usually start to fly at the beginning of June) to identify the specimens.The study of the types of taxa synonymized with Merodon avidus was based on the above-mentioned diagnostic morphological characters, geographical distribution and when the adults were active (early spring IV-V, late spring-summer VI-VII, autumn VIII-IX).(Hurkmans, 1993).Hurkmans (1993) re-established the name, but without designation of a neotype.Our results for populations in the area where the type was collected indicate the existence of only one taxon of the M. avidus complex.In order to preserve the stability of nomenclature, we have designated a neotype from the material collected in Toscana province: ♂, Italy, Pisa, Parco Regionale Migliarino San Rossore, 10.v.2012 leg.Vujić (FSUNS).

Merodon ibericus
distribution.Merodon avidus occurs in the Mediterranean basin (Spain (coastal), southern parts of France, Italy, the Adriatic part of Croatia, southern Bosnia and Herzegovina, Montenegro, Albania, FYR Macedonia, Greece, Turkey, Israel and Cyprus), Serbia, eastern Bulgaria and the Crimea (Fig. 7A).
Merodon ibericus is endemic to the Iberian Peninsula, where it is widespread, mostly in the mountainous parts (Cantabrian mountain range, Pyrenees, Sistema Central and Portugal) (Marcos-García et al., 2007).Based on our results this species occurs also in the Atlas Mountains in Morocco (specimens included in analysis, see Supporting Information Table S2).
Merodon moenium mainly occurs in the continental part of Europe (France, Netherlands, Germany, Poland, Czech Republic, Denmark, southern Sweden, Switzerland, Slovakia, Slove-  Merodon avidus and M. moenium are not geographically isolated.They overlap in a narrow zone at the edge of their range.

dnA barcoding
A total of 27 COI mtDNA haplotypes, defined by 36 variable sites (singleton variable sites: 7 and parsimony informative sites: 29), were recorded in the 41 specimens of M. ibericus, M. avidus and M. moenium analyzed.Haplotype diversity for total dataset was h = 0.971 ± 0.013, nucleotide diversity, π = 0.0206 and the average number of nucleotide differences, k = 12.998.
There were fifteen haplotypes recorded for M. ibericus, seven for M. avidus and five for M. moenium, none of which were shared (Supporting Information Table S3).
In all trees, COI haplotypes formed two main clades: the first corresponds to M. ibericus and the second to M. avidus + M. moenium (Fig. 8, Supporting Information Figs S1 and  S2).
The high number of unique haplotypes recorded in both clades indicates marked intraspecific polymorphism.Pairwise haplotype nucleotide divergences ranged from 0.002 to 0.147 (Supporting Information Table S4).The pairwise nucleotide diversity values of the different haplotypes recorded in M. ibericus ranged between 0.002 and 0.015, while in M. avidus / M. moenium clade these values ranged between 0.002 and 0.009.
COI barcode failed to separate M. avidus and M. moe nium in cluster analyses, but the lack of shared haplotypes encouraged us to analyze molecular variance (AMOVA).The data obtained revealed statistically significant percentage of variation among species (Table 1).

Allozyme variability and population structure
Four of the five enzyme loci analyzed were polymorphic in all the populations of the Merodon avidus / moenium pair (GPD, IDH, GPI and ME) studied.At the AAT locus, however, there was one allele (a) in all populations that were initially identified as M. avidus, and another allele (b) in all populations that were identified as M. moenium (Table 3, Fig. 9).This indicates that the AAT locus can be considered diagnostic for delineating these cryptic taxa in the Mero don avidus complex, as previously suggested by Milankov et al. (2001).
In the 19 populations analyzed the mean number of alleles per locus ranged from 1.40 (Orjen, Montenegro, M. moenium) to 2.33 (Piemonte Cuneo, Italy, M. avidus).The mean number of alleles considering all loci and all populations of M. avidus and M. moenium was 2.4.With the exception of the AAT locus, which was monomorphic and species specific for M. avidus and M. moenium, no rare (frequency < 0.1) or unique alleles were found in any of the populations studied.However, IDH locus with "a" allele mainly present in M. avidus and "b" allele mainly in M. moenium populations, differed statistically significantly in allele frequencies in these two taxa (p ≤ 0.01) (Table 1, Fig. 9).In 6 of 11 populations of M. avidus, "a" was the only allele present, while it was a major allele (frequency > 0.5) in the population at Piemonte Cuneo, Italy (0.800).The lowest frequency of this allele was recorded in the population in the Bozdağ Mountains, Turkey (0.200).The major allele in the populations of M. moenium from Djerdap, Serbia, Trento Bleggio Superiore, Italy, Kopaonik, Serbia and Durmitor, Montenegro was "b", and it was a fixed allele in populations from Orjen, Montenegro and Stara Planina, Serbia.In populations from Fruška Gora, Serbia and Castiglione dei Pepoli Monte Baducco, Italy "b" was the rarest allele compared to other M. moenium populations and equal in frequency (0.5) with the "a" allele.There were also statistically significant differences between M. avidus and M. moenium at the ME locus in frequencies of the "a" (p ≤ 0.01) and "c" alleles (p ≤ 0.05).
The hierarchical structure of the M. avidus / moenium pair was analyzed using main parameters of F statistic for the total complex (divided into 19 populations), as well as for the same total population divided into 2 large (meta) populations (M.avidus and M. moenium) (Table 4).F statistics were also separately calculated for M. avidus (11 populations) and M. moenium (8 populations).
Mean F st value (measure of genetic differentiation between populations) was 0.571 for M. avidus and 0.434 for M. moenium.Differentiation among M. avidus populations  was mainly caused by differences in allele frequencies at GPI (F st = 0.678), IDH (F st = 0.664) and ME (F st = 0.585), followed by the GPD locus (F st = 0.405).The value of F is was negative for most of the loci analyzed (GPD, IDH and AAT).Positive value of F is for the IDH locus indicates an excess of homozygosity.However, the F st > F is relation recorded for each locus, including IDH, suggest that genetic drift and not local inbreeding affected this result.Differentiation among M. moenium populations was mainly associated with differences in allele frequencies at GPI (F st = 0.701) and ME loci (F st = 0.581).Negative values of F is were recorded for all loci, except IDH.Also, IDH was the only locus in M. moenium for which the results indicated that local inbreeding seemed to be the more important factor (F st < F is ).
Analysis of all the M. avidus / moenium material revealed considerable genetic differentiation among the 19 populations analyzed (F st = 0.654).The value of F is was only positive for the IDH locus and the F st value was higher than that of F is .On the other hand, when the sample was divided into the two presumed cryptic taxa (M.avidus and M. moenium), the F st < F is relation recorded for the IDH locus indicated inbreeding in subpopulations.This was also the case for the ME locus in the same sample.Mean F st value (0.238) implied a noticeable degree of genetic differentiation between M. avidus and M. moenium (meta) populations (Table 4).
The genetic structure of M. avidus and M. moenium populations was also estimated using values of expected (He) and observed (Ho) heterozygosity.Average values of Ho recorded for M. avidus (0.452) and M. moenium (0.433) were similar, whereas the He values were 0.388 and 0.435, respectively.However, the χ 2 test revealed that the observed genotype frequencies deviated significantly from Hardy-Weinberg equilibrium for IDH (p ≤ 0.01), ME (p ≤ 0.01) and GPI loci (p ≤ 0.05) in M. avidus, and IDH locus (p ≤ 0.05) in M. moenium (Table 3).
Analysis of allele distribution was quantified using unbiased measures of genetic identity and genetic distance (Nei, 1978).Average value of the genetic identity (0.520) between M. avidus and M. moenium was lower than the average value for the genetic distance (0.646) between these populations.In addition, a UPGMA tree based on Nei's (1978) genetic distance indicates two clusters: a M. avidus cluster and M. moenium cluster (Fig. 10).
The correlation between F st values (Weir & Cockerham, 1984) and geographical distances was calculated for each population pair from 19 populations.Mantel's test in the IBD program did not indicate either a significant correlation between genetic distance and geographic distance (r = -0.023)(Fig. 11) or a significant isolation by distance (p > 0.05 for 10,000 randomizations).

ecological niche analysis
To identify ecological preferences, PCA reduced the environmental variables to 4 PC axes that best described the geographical ranges of the species investigated in the M. avidus / moenium pair.The first 4 PCs explain 92% of total variation (Table 5).Separation of the species was highly significant (p ≤ 0.01) in terms of factor scores for all 4 PCs using ANOVA and the Tukey test (Table 6) and all PCs were used in the interpretation of results.
Precipitation accounted for 70% of the environmental variance (PC1 and PC2).PC1 explained most of the variation (43%) between M. avidus and M. moenium.It was negatively correlated with BIO18, BIO14 and BIO17, and related with precipitation in the warmest quarter, as well as in the driest month and quarter.PC2 was correlated with BIO16, BIO13, BIO19 and BIO12, and related to precipitation levels in the wettest quarter and month, and in the coldest quarter.PC3 was correlated with mean diurnal and annual temperature range (BIO2, BIO7, and BIO3) and explained 13% of variation, while PC4 was negatively correlated with altitude (alt) (Table 5).Altitude explained only a small part of the differences in ecological niches (9%) and M. avidus appeared to have a slightly greater altitudinal range than M. moenium.The scatter plot of PCs depicted the positions of M. avidus and M. moenium in environmental space (Fig. 12).Although ecological prefer-Fig.11.A Mantel's test for correlation between F st (Weir & Cockerham, 1984), genetic distance and geographical distance (km).In previous studies, the phenotypic diversity of the specimens and populations in the M. avidus complex led to the first hypothesis that this taxon is a complex of at least two cryptic taxa.Morphological differences were supported by an analysis of allozyme variability and fixed differences at two diagnostic allozyme loci (AAT and IDH) and confirmed the separate identity of M. avidus A and M. avidus B (Milankov et al., 2001).Genetic differences were also recorded between spring Mediterranean populations and montane populations that were morphologically almost inseparable.Phenotypic differences between the spring-early summer (May-June) and late summer-autumn generations (July-September) of Mediterranean samples pointed to  possible genetic differentiation between conspecific M. avidus A populations.However, Milankov et al. (2001) were limited by the small extent of the geographical area they sampled and concluded the variations were intraspecific.
Earlier results based on DNA sequences of COI, in line with a barcoding approach, brought additional ambiguities to the taxonomic issues associated with the M. avidus complex.Sequencing of the mtDNA COI-3' region revealed extensive haplotype variation in the M. avidus complex, but the selected marker failed to discriminate the evolutionarily independent genetic units previously identified as M. avidus A and M. avidus B (Milankov et al., 2009, Popović et al., 2014).On the other hand, there are three unique haplotypes in the morphologically indistinguishable species M. ibericus (described as Merodon bicolor, Milankov et al., 2009) from the Iberian Peninsula, which form a separate clade in the strict consensus tree based on a parsimony analysis (Milankov et al., 2009).However, this was based on few specimens (only three specimens included).
In order to include all the available genetic markers in our integrative taxonomy approach for delimitating sibling species in the Merodon avidus complex, we used the 5' fragment of the COI mtDNA gene (DNA barcode).This sequence is widely recognized as a "barcoding" sequence and is one of the few molecular markers used in genetic studies of Syrphidae.Different authors encourage the use of DNA barcoding for identifying species (Hebert et al., 2003;Tautz et al., 2003;Gaston & O'Neill, 2004).However, in some species in the family Syrphidae, the sequence divergence in COI is relatively low and even invariant.Mengual et al. (2008) report cases where morphological features and DNA barcoding data differ and others where there are differences in morphology but the COI sequences are identical.The use of DNA barcoding for identifying the species in the Merodon avidus complex led to interesting conclusions.Molecular phylogenetic analyses and NJ, ML and MP trees confirmed the presence of two clades, one corresponding to M. ibericus and the other to M. avidus / M. moenium.Thus, DNA barcoding supported the previous conclusion that M. ibericus is a separate species in the M. avidus complex characterized by a longer coalescent time of divergence and isolated diversification.On the other hand, cluster analyses of DNA barcoding sequences did not delimit M. avidus and M. moenium, whereas allozyme variability and ecological niche analyses did.Despite this, more detailed analyses of sequence variability revealed that these two sibling species do not share haplotypes.Detailed analysis of variable sites and AMOVA revealed statistically significant genetic divergence between M. avidus and M. moenium.Calculated Φ st values for these sibling species are statistically significant (p < 0.01) and thus support the results of the allozyme variability and ecological niche analyses.DNA barcoding sequences for these sibling species have the same polymorphic sites, which indicates that these species separated recently and a longer period of time is needed for the establishment of fixed differences.It is documented that recently separated species are more difficult to distinguish using barcoding as these species may share polymorphic sites that were polymorphic in the ancestral species (Austerlitz et al., 2009).In other words, the number of mutations at the COI locus separating two individuals increases with their coalescence time (Austerlitz et al., 2009).If M. avidus and M. moenium are sibling species that recently originated from a common ancestor, this would account for the shared COI-3' haplotypes in these taxa reported by Milankov et al. (2009) andPopović et al. (2014).Since the DNA barcoding sequence is less variable than the 3' fragment of COI (Herbert et al., 2003a, b) this might account for us not finding any shared haplotypes.
Furthermore, it is very likely that the high endemicity on the Iberian Peninsula (Marcos-García et al., 2007), probably due to geographical isolation during the Pleistocene glaciations (Milankov et al., 2008), contributed to the separation of M. ibericus from M. avidus and M. moenium.This would mean that a longer coalescence time for M. iberi cus and M. avidus / M. moenium than for M. avidus and M. moenium resulted in the better initial discriminating power of the COI-DNA barcoding approach.
Our enzyme variability analysis revealed species-specific alleles at the AAT locus and statistically very significant differences between the frequencies of alleles at the IDH locus (each of the alleles was dominantly present in each putative species).This is partially in accordance with the previous enzyme study of these species, which indicates fixed alleles in each cryptic taxon at the IDH locus, but the AAT alleles are not completely fixed (common allele has a frequency of 1.00 in M. avidus A population and 0.083 and 0.09 in two M. avidus B populations) (Milankov et al., 2001).Based on the present study we can confirm that the IDH and AAT loci are highly diagnostic and can be used to separate these two cryptic taxa in the Merodon avidus / moenium pair.We interpret the recorded differences in the presence of alleles / frequency data as a consequence of including a geographically broader sample in our study (in addition to specimens from the Balkan Peninsula, we included specimens from the Apennine Peninsula and Turkey), which also increased the reliability of our results.The value of Wright's F st recorded in this study (F st = 0.631) confirmed the presence of two taxa, which is in accordance with the previous result of 0.511 (Milankov et al., 2001).Moreover, in order to obtain more precise information about the genetic differentiation of the taxa studied, we also analyzed two large (meta) populations, of M. avidus and M. moenium, separately.The value of Nei's measure of the genetic identity (I = 0.520) of these two populations was important for our next conclusion.Based on the value of this variable for sibling species in the Drosophila wil listoni group (Diptera: Drosophilidae) Ī = 0.517 ± 0.024 (Ayala et al., 1975), which is similar to our result, we conclude that M. avidus and M. moenium are probably sibling species.This is in accordance with the fact that sibling species are reproductively isolated, but often morphologically indistinguishable (Steyskal, 1972).Their reproductive isolation is a consequence of the fact that they can interbreed, but their offspring cannot reproduce.However, their ability to interbreed could explain the occasional presence of other alleles at species-specific enzyme loci, as recorded here and in previous allozyme studies.All the above indicate that the M. avidus complex, apart from the Iberian endemic M. ibericus, consists of two distinct sibling species and justifies the initial decision of introducing different names for these taxa.
For the first time, populations from the isolated mountain Fruška Gora (Pannonian Plain, Serbia) were included in our allozyme analysis of the M. avidus species complex.Based on our results, these populations belong to the M. moenium cluster in the genetic distance tree, which is the opposite of what Milankov et al. ( 2009) conclude.Based on morphology and geography, populations from Fruška Gora fit well with others from the Pannonian Plain, all of which belong to the species M. moenium.
A central feature of the present study was to determine potential connections between morphology, genetic data and seasonal preferences of certain populations of the M. avidus species complex.Two key localities in resolving this issue were the Djerdap gorge (continental locality) and Umag (Adriatic coast-Mediterranean locality).Based on the presumed seasonal pattern of distribution of these species at localities where they both occur, the spring generation from Djerdap was initially identified as M. moenium and the autumn generation as M. avidus.The UPGMA genetic distance tree based on allozyme data grouped the spring population in the M. moenium cluster and the autumn population in M. avidus cluster, which corresponds with the prior determination of specimens (Fig. 10).This confirmed a temporal divergence of populations of M. avi dus and M. moenium at the same locality and justifies using this information to correctly identify morphologically similar specimens.It also indicates that difference in the seasonal occurrences of certain populations might lead to genetic differentiation of these sibling species.
Due to the typical Mediterranean climate at Umag, both of the populations sampled (spring and autumn generation) at this locality were identified as M. avidus although they differed slightly morphologically.Sampling two generations of presumably the same taxon at one locality was intended to clarify whether there was a genetic subdivision between M. avidus A populations, previously established because of the relatively high intraspecific divergence and phenotypic substructuring within this taxon (Milankov et al., 2004(Milankov et al., , 2009)).In a trial sample, we marked these two generations as two separate populations and included them in the genetic distance analysis.The resulting UPGMA tree grouped both populations in a cluster within the large M. avidus cluster.Nei's unbiased value of the genetic identity between them was 0.99.This confirmed the assumption that there is a continuous annual occupancy by M. avidus populations of localities with a Mediterranean climate, which is of great importance for the correct identification of specimens belonging to the M. avidus species complex.The fact that both populations belong to the same species indicated that the morphological differences between the Umag populations collected in different seasons probably resulted from variations in environmental conditions experienced by their developmental stages.Furthermore, their grouping into the nearest cluster and its very high genetic identity value indicated that there was no genetic substructuring within M. avidus at this locality.Consequently, we recorded both generations (June and September) of M. avidus from Umag (Croatia) as a single population in our final analysis (Fig. 10).Significant ecological divergence between M. avidus and M. moenium may also contribute to allozyme results in species determination.Since it is known that ecological preferences of closely related and cryptic species are unique (Penman et al., 2005), and that a highly consistent and significant positive association between ecological divergence and reproductive isolation exists (Funk et al., 2006), evidence that some taxa have different niches may be a useful criterion in species delimitation.The ecological niches of M. avidus and M. moe nium differed in four dimensions (PCs) that are correlated with precipitation, temperature and altitude.Precipitation explained the main differences, with 70% of total variation (PC1 and PC2), while temperature and altitude appeared to account for a smaller amount of variation.These ecological shifts in four independent axes and great % of variance accounted for by precipitation are important evidence of speciation (Nosil et al., 2009).
Precipitation during the warmest and driest quarter of the year is an important limiting environmental factor and is responsible for the narrower ecological niche of Mero don moenium.In contrast, Merodon avidus occurs across a wide temperature and precipitation range, mostly within the Mediterranean climate zone, in contrast to M. moenium that prefers colder areas with higher humidity.M. avidus also has a slightly greater altitude range than M. moenium.Adults of M. avidus are active from spring to autumn (in a Mediterranean climate), while M. moenium cannot be found after the middle of August.It is important to emphasize that overlap in environmental space occurs because these species are not geographically isolated: they share a narrow part of their range.Occasionally both species were recorded at the same locality but in those circumstances they do not appear simultaneously or in the same microhabitat.The seasonal preferences of these species were of great importance, especially in initial identification of specimens with doubtful morphological characteristics.
It can be assumed that environmental divergence of the species investigated is a result of the specific geological history of the area studied.A large part of the European continent was shaped during the last glaciations, which ended approximately 10,000 years ago (Savić, 2008).The present overlap in environmental space is indicative of recent speciation, which we previously hypothesized on the basis of genetic data.concluSIon From a parallel study involving re-examination of type material, documented diagnostic morphological characters, season of adult activity and geographical distribution, genetic variability data and ecological niche analyses, we justify the use of the following names for three closely related taxa in the Merodon avidus complex: M. avidus (Rossi, 1790) [M. avidus A in Milankov et al. (2001, 2009) and Marcos-Garcia et al. (2007)], M. moenium (Wiedemann, 1822) [M. avidus B in Milankov et al. (2001, 2009)] and M. ibericus Vujić nom. n. (new name for M. bicolor Gil Collado, 1930, preoccupied by M. bicolor Walker, 1852).Based on the genetic variability data presented in this study we confirm the presence of three cryptic taxa within the Merodon avidus complex and conclude that they can be recognized as sibling species.For the sibling pair of species M. avidus and M. moenium we further conclude that an ecological niche analysis, even when only physical environmental parameters are used, as in this instance, can aid in the differentiation of populations of putative cryptic species.The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree.The evolutionary distances were computed using the Tamura 3-parameter method and the rate variation among sites was modeled with a gamma distribution (γ = 0.05).

Fig. 10 .
Fig. 10.Dendrogram of populations of the Merodon avidus species complex based on Nei's (1978) genetic distance (Method = UPGMA, Modified from NEIGHBOR procedure of PHYLIP Version 3.5).Branch length values are indicated above the nodes.a -autumn flight, s -spring flight.

Fig. 12 .
Fig. 12. Scatter plot of factor loadings of tree PC scores, showing position of M. avidus and M. moenium in environmental space.A -PC1 against PC2; B -PC1 against PC3.
Fig. S2.Maximum Likelihood tree based on 5' fragment COI mtDNA sequences in the Merodon avidus complex.The percentage of trees in which the associated taxa clustered together is shown next to the branches.The tree is drawn to scale, with branch lengths measured in the number of substitutions per site.

Fig
Fig.S1.Neighbour-joining tree based on 5' fragment COI mtDNA sequences in the Merodon avidus complex.The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) are shown next to the branches.The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree.The evolutionary distances were computed using the Tamura 3-parameter method and the rate variation among sites was modeled with a gamma distribution (γ = 0.05).

table 3 .
Allelic frequency and values of heterozygosity for all loci in the Merodon avidus/moenium pair.

table 4 .
F statistics for populations of the M. avidus/moenium pair.Values in bold indicate greater importance of local inbreeding at particular loci, as well as the mean F st values averaged for all loci.
ences were significantly different, with high significance in all four PCs, they shared part of environmental space.Overlapping was related to localities where these species occur sympatrically.At these localities (Serbia: Malinik-Dubašnica, Djerdap; Greece: Pindos; Montenegro: Kanjon Komarnice, Orijen; Italy: Castiglione dei Pepoli) clear temporal succession or microhabitat separation was recorded.dIScuSSIon

table 5 .
Principal component analysis of 19 Bioclim variables and altitude associated with the occurrence of M. avidus and M. moenium.Significant factor loadings are printed in bold.

table 6 .
ANOVA and Tukey HSD post hoc test results of comparison factor scores between M. avidus and M. moenium.

table S2 .
List of specimens selected for DNA barcoding analysis.Accession number belongs to haplotype deposited in GenBank, specimens which share haplotype have same accession number.