Introduction
Ηoverflies (Diptera: Syrphidae), also known as flower flies, constitute a cosmopolitan and highly diverse insect group with more than 6100 taxa described globally (Thompson et al., Reference Thompson, Rotheray, Zimbado, Brown, Borkent, Cumming, Wood, Woodley, Norman and Zumbado2010). Hoverflies influence ecosystems in many ways, such as: (a) pollinating a wide range of flowering plants; (b) controlling plant pests of which they are effective predators; (c) having phytophagous larvae feeding on bulbs; and (d) effectively recycling nutrients from dead matter (Rotheray & Gilbert, Reference Rotheray and Gilbert2011). Because of their heterogeneous character and their wide distribution, hoverflies constitute an important and diverse group for ecological and biogeographical studies.
One of the most species-rich of the hoverfly genera is Eumerus Meigen, 1822, with 256 species recorded worldwide (Pape & Thompson, Reference Pape and Thompson2015), of which 89 occur in the western Palaearctic region (total species number of the entire Palaearctic is 140: Peck, Reference Peck, Soos and Papp1988; Speight, Reference Speight2014) and 72 in the Afrotropical region (Pape & Thompson, Reference Pape and Thompson2015). The only existing comprehensive key for Eumerus is that of Stackelberg (Reference Stackelberg1961) based on only some of the Palaearctic taxa, which was used by Speight (Reference Speight2014) to compile a list of European Eumerus. Due to the lack of an up-to-date European identification key, species delimitation is often not feasible, especially for poorly studied areas, e.g., the Mediterranean. In addition, the taxonomic status of a considerable number of taxa is uncertain, with confusing nomenclature and synonymies frequently present (Peck, Reference Peck, Soos and Papp1988). It is often unclear to which taxon many names refer, and a broad revision of the genus is needed, perhaps employing more sophisticated taxonomical tools, e.g., molecular systematics.
Over the last two decades, DNA barcoding has been introduced to taxonomy and has greatly expedited species identification, assisted in species delimitation and elucidated species evolution and biology (Rubinoff, Reference Rubinoff2006). DNA barcoding can be a fast and efficient way to identify species, to diagnose new species, and to provide molecular operational units for ecological and biodiversity studies (DeWalt, Reference DeWalt2011). However, opposing opinions exist regarding the application of DNA barcoding as a primary means of species delimitation. Rubinoff (Reference Rubinoff2006) claimed that mitochondrial DNA (mtDNA) is not adequate as a sole source of species-defining data because of reduced effective population size and introgression, maternal inheritance, recombination, inconsistent mutation rate, heteroplasmy and compounding evolutionary processes. As a consequence, there are cases where some sister taxa cannot be identified because they have identical or nearly identical DNA barcodes, giving a false negative signal of species differentiation (Moritz & Cicero, Reference Moritz and Cicero2004). In contrast, taxa with a wide geographical distribution may exhibit relatively large genetic divergence and, thus, might present a false positive signal (i.e., incorrectly indicating the occurrence of differentiated taxa) (Avise, Reference Avise2000). Notwithstanding the controversy about its effectiveness, DNA barcoding has been highly effective in species identification, especially for resolving insect taxonomy (Ball et al., Reference Ball, Hebert, Burian and Webb2005; Smith et al., Reference Smith, Fisher and Hebert2005). To overcome any false phylogenetic signal, it is recommended to use multiple and complementary tools to better delimit biodiversity (‘integrative taxonomy’: Dayrat, Reference Dayrat2005). This implies that morphological features, and molecular (e.g., DNA sequences) and biochemical (e.g., alloenzymes) data, as well as morphometric data should be utilized for species delimitation (Mengual et al., Reference Mengual, Ståhls, Vujić and Marcos-García2006).
The mitochondrial cytochrome-c oxidase subunit I gene (COI, cox1) is widely used in species identification due to its phylogenetic signal, which can discriminate species (Hebert et al., Reference Hebert, Cywinska, Ball and DeWaard2003a , Reference Hebert, deWaard and Landry2010; Hebert & Gregory, Reference Hebert and Gregory2005). A fragment of approximately 650 bp of the 5′ end of the COI, i.e., the ‘Folmer’ fragment (Folmer et al., Reference Folmer, Black, Hoeh, Lutz and Vrijenhoek1994), is the most frequently used gene fragment in DNA barcoding of animals, and this approach has been successfully applied to hoverflies (Pérez-Bañón et al., Reference Pérez-Bañón, Rojo, Ståhls and Marcos-García2003; Ståhls et al., Reference Ståhls, Vujić, Pérez-Bañón, Radenković, Rojo and Petanidou2009; Gibson et al., Reference Gibson, Kelso, Jackson, Kits, Miranda and Skevington2011; Marcos-García et al., Reference Marcos-García, Vujić, Ricarte and Ståhls2011; Radenković et al., Reference Radenković, Vujić, Ståhls, Pérez-Bañón, Rojo, Petanidou and Simić2011). Furthermore, primers amplifying a fragment of the 3′ end of the COI have been effectively applied to obtain longer sequences, i.e., >800 bp. According to Hebert et al. (Reference Hebert, Ratnasingham and deWaard2003b ) COI-3′ regions present similar sequence divergence profiles as COI-5′ fragments, which imply that COI-3′ can be used as an alternative DNA barcode. Indeed, the COI-3′ region has proven valuable for solving taxonomical uncertainties in various hoverfly genera (Milankov et al., Reference Milankov, Stamenković, Ludoški, Ståhls and Vujić2005, Reference Milankov, Ståhls, Stamenković and Vujić2008a , Reference Milankov, Ludoški, Ståhls, Stamenković and Vujić2009, Reference Milankov, Francuski, Ludoški, Ståhls and Vujić2010a , Reference Milankov, Francuski, Ludoški, Ståhls and Vujić b , Reference Milankov, Ludoški, Ståhls and Vujić2013; Mengual et al., Reference Mengual, Ståhls, Vujić and Marcos-García2006; Vujić et al., Reference Vujić, Radenković, Ståhls, Ačanski, Stefanović, Veselić, Andrić and Hayat2012; Francuski et al., Reference Francuski, Djurakic, Ludoški and Milankov2013).
Molecular and morphological studies have been carried out on many hoverfly genera such as Merodon (Milankov et al., Reference Milankov, Ståhls, Stamenković and Vujić2008a , Reference Milankov, Ståhls and Vujić b ; Ståhls et al., Reference Ståhls, Vujić, Pérez-Bañón, Radenković, Rojo and Petanidou2009), Cheilosia (Ståhls et al., Reference Ståhls, Stuke, Vujić, Doczkal and Muona2004; Milankov et al., Reference Milankov, Francuski, Ludoški, Ståhls and Vujić2010a , Reference Milankov, Francuski, Ludoški, Ståhls and Vujić b ), Chrysotoxum (Nedeljković et al., Reference Nedeljković, Ačanski, Vujić, Obreht, Dan, Ståhls and Radenković2013) and tribe Pipizini (Vujić et al., Reference Vujić, Ståhls, Ačanski, Bartsch, Bygebjerg and Stefanović2013). In these studies, an integrative taxonomic approach was applied, i.e., complementary use of molecular markers (mtDNA and/or nuclear gene fragments) and morphological characters, which often provided rather good taxonomic resolution (Milankov et al., Reference Milankov, Ståhls, Stamenković and Vujić2008a , Reference Milankov, Ståhls and Vujić b ). However, in certain cases, different taxa shared the same haplotypes (Milankov et al., Reference Milankov, Ståhls and Vujić2008b , Reference Milankov, Francuski, Ludoški, Ståhls and Vujić2010a , Reference Milankov, Francuski, Ludoški, Ståhls and Vujić b ; Ståhls et al., Reference Ståhls, Vujić, Pérez-Bañón, Radenković, Rojo and Petanidou2009) or possessed character differences of solely one nucleotide (Ståhls et al., Reference Ståhls, Vujić, Pérez-Bañón, Radenković, Rojo and Petanidou2009; Milankov et al., Reference Milankov, Francuski, Ludoški, Ståhls and Vujić2010b ). In the latter cases, this outcome indicated the insufficiency of a COI-based identification system alone, to delimit species or species complexes within the genus Merodon (Milankov et al., Reference Milankov, Ståhls and Vujić2008b ; Ståhls et al., Reference Ståhls, Vujić, Pérez-Bañón, Radenković, Rojo and Petanidou2009), and underlined the importance of integrative taxonomic inference.
Here, our aim is to validate usage of COI-3′ region for species delimitation of the genus Eumerus using different sequence lengths in order to identify species and explore intra- and interspecific variability. We address two specific questions: (1) Can COI-3′ barcoding reveal intra- and interspecific genetic variation in Eumerus? (2) Is a short sequence of COI-3′ (647 bp) sufficient to accurately resolve species taxonomy compared to a more elongated one (746 bp)?
Materials and methods
Taxon sampling
We used dry pinned specimens deposited in two entomological collections: the Melissotheque of the Aegean located at the University of the Aegean, Greece (M-UAegean; Petanidou et al., Reference Petanidou, Ståhls, Vujić, Olesen, Rojo, Thrasyvoulou, Sgardelis, Kallimanis, Kokkini and Tscheulin2013); and the collection of the Faculty of Sciences at the Department of Biology and Ecology of the University of Novi Sad, Serbia (FSUNS). The specimens used in this study were collected from 2009 to 2014. In total, 75 specimens were used from two geographical areas: 69 derived from the western Palaearctic and six from Afrotropical region (RSA) (table S1). Initial species identification was based on morphology. Specimens were assigned to 28 taxa (table S1) of which three had been previously undescribed (Eumerus aff. barbarus Coquebert, 1804; E. aff. rubiginosus Lyneborg, in litt. and E. aff. tarsatus Lyneborg, in litt.). The selected taxa are mainly of Palaearctic origin (22 taxa), while six taxa originated from the RSA. Identifications were carried out by Dieter Doczkal (Afrotropical taxa) and Ante Vujić (Palaearctic taxa). Taxa identification for specimens derived from the RSA was based on Lyneborg's revision (Lyneborg, in litt.). A list of the specimens used in the analysis, together with their GenBank accession numbers and collection data, is given in table S1.
To detect intraspecific variation, Eumerus amoenus Loew, 1848 was analyzed, since its broad distribution in the Mediterranean appears to be the widest within the genus. Two populations from the Aegean islands of Lesvos (three specimens) and Samos (seven specimens) were assayed. Measurements of interspecific variation were determined through the taxa, E. pulchellus Loew, 1848 and E. pusillus Loew, 1848 by comparing the genetic distances between adjacent and distant populations. For E. pulchellus, molecular diversity indices were calculated by using one specimen for each geographical area (i.e., Chios, Dadia, Lesvos, Limnos, Rhodes, Samos and Sardinia). Molecular diversity indices were also calculated for E. pusillus, but in this case, two or more samples per geographical area (Chios, Crete, Karpathos and Naxos) were selected.
DNA extraction and PCR amplification
Total genomic DNA was extracted using the head and/or two to three legs from each specimen. We used the protocol by Chen et al. (Reference Chen, Rangasamy, Tan, Wang and Siegfried2010) for Sodium Dodecyl Sulfate (SDS) extraction with the following modifications: (a) RNase A solution was not added; (b) the concentration of proteinase K solution was 40 mg ml−1; and (c) an additional step of chloroform/isoamyl alcohol (24:1) was applied. Samples were re-suspended in 30 µl of TAE buffer.
PCR amplifications of the COI-3′ were performed in a total volume of 25 µl, containing 25 ng µl−1 template DNA, 5 pmol µl−1 of each primer, 0.08 mM of dNTPs, 1× Reaction Buffer (Thermo Scientific, USA) and 1.25 units of Polymerase (Taq poly or Dream Taq poly, Thermo Scientific, USA). Amplifications were performed in an Authorized PCR Thermal Cycler (Mastercycler® personal, Eppendorf, Germany). Thermocycling conditions consisted of initial denaturation at 95°C for 2 min, 29 cycles of 30 s denaturing at 94°C, 30 s annealing at 49°C, 2 min extension at 72°C, followed by a final extension of 8 min at 72°C (Vujić et al., Reference Vujić, Ståhls, Ačanski, Bartsch, Bygebjerg and Stefanović2013). We employed universally conserved primers to amplify and sequence the COI-3′: forward primer C1-J-2183 (5′-CAACATTTATTTTGATTTTTTGG-3′) (alias JERRY) and reverse primer TL2-N-3014 (5′-TCCAATGCACTAATCTGCCATATTA-3′) (alias PAT) (Simon et al., Reference Simon, Frati, Beckenbach, Crespi, Liu and Flook1994). Amplified products were run on 1.5% agarose gels for visual inspection. Purification of the PCR products was done with the ExoSap-IT kit (USB, Cleveland, OH, USA) and clean products were thereafter Sanger sequenced in both directions on an ABI 3730 DNA analyzer (Applied Biosystems, USA) at the Sequencing Service laboratory of the Finnish Institute for Molecular Medicine (http://www.fimm.fi).
Sequence alignment
Two datasets of different sequence length were produced. Dataset A comprised sequences of 647 nucleotides in total length, obtained by forward sequencing of the COI-3′ region, from 75 Eumerus specimens and four outgroups. We chose Platynochaetus setosus Fabricius, 1794 (Accession No. KM224512), Merodon erivanicus Paramonov, 1925 (Accession No. KT157919) and two species of the genus Megatrigon Johnson, 1898 (Accession No. KT157920 and KT157921) as outgroups. Dataset B consisted of full length of the COI-3′ (746 nucleotides), acquired through bidirectional sequencing, for the aforementioned 75 Eumerus taxa and outgroups. All trees were rooted based on P. setosus sequence.
As required, sequences were edited by eye using BioEdit 7.2.5 (Hall, Reference Hall1999). For multiple sequence alignment we employed the L-INS-i algorithm, which is considered to be more accurate as an iterative refinement method incorporating local pairwise alignment information (Katoh et al., Reference Katoh, Kuma, Toh and Miyata2005). Alignments were implemented using MAFFT version 7 (http://mafft.cbrc.jp/alignment/server/index.html). Both datasets were trimmed to their final lengths using BioEdit 7.2.5 (Hall, Reference Hall1999). Polymorphic sites, DNA polymorphism and basic molecular diversity indices were calculated using DnaSP 5.10.01 (Librado & Rozas, Reference Librado and Rozas2009), which also generated nexus files. Sequences were checked for possible presence of stop codons (Buhay, Reference Buhay2009) using the Mesquite 2.75 system for phylogenetic computing (Maddison & Maddison, Reference Maddison and Maddison2011). All sequences were translated using the invertebrate mitochondrial code. The evolutionary models used for maximum-likelihood (ML) and Bayesian inference (BI) analyses were implemented in the HIV sequence database (http://www.hiv.lanl.gov/content/sequence/findmodel/findmodel.html) (Posada & Crandall, Reference Posada and Crandall2001).
Molecular data analysis
PTP models
Several studies have highlighted the ability of Poisson tree processes (PTP) models to reveal and resolve taxonomic issues (Leasi & Norenburg, Reference Leasi and Norenburg2014; Soldati et al., Reference Soldati, Kergoat, Clamens, Jourdan, Jabbour-Zahab and Condamine2014; Tang et al., Reference Tang, Humphreys, Fontaneto and Barraclough2014). Because this approach does not require ultrametrization of trees (and its associated biases), it constitutes a reasonable alternative to other species delineation models such as the General mixed Yule coalescent model (Pons et al., Reference Pons, Barraclough, Gomez-Zurita, Cardoso, Duran, Hazell, Kamoun, Sumlin and Vogler2006). In PTP models, the numbers of substitutions (branch lengths) represent speciations or branching events and, therefore they only require a phylogenetic input tree. PTP models have previously been implemented to reveal putative molecular species clusters (Zhang et al., Reference Zhang, Kapli, Pavlidis and Stamatakis2013). PTP analyses were conducted on the web server for PTP (available at http://species.h-its.org/ptp/) using the best ML tree resulting from the RA × ML analysis (see below).
Maximum parsimony (MP)
MP analyses were performed in the program NONA (Goloboff, Reference Goloboff1999), spawned in WINCLADA version 1.00.08 (Nixon, Reference Nixon2002). A heuristic search algorithm with 1000 random addition replicates (mult × 1000) was performed, holding 100 trees per round (hold/100), max trees set to 100,000, and applying TBR branch swapping.
Maximum likelihood
ML trees were generated using RAxML 8.0.9 (Stamatakis, Reference Stamatakis2006; Stamatakis et al., Reference Stamatakis, Hoover and Rougemont2008) in the Cipres Science Gateway (Miller et al., Reference Miller, Pfeiffer and Schwartz2010) under the general time-reversible (GTR) evolutionary model with a gamma distribution (GTR + G) (Rodriguez et al., Reference Rodriguez, Oliver, Marin and Medina1990) and 1000 bootstrap replicates.
Bayesian inference
BI topologies were assessed using MrBayes 3.2.6 (Huelsenbeck & Ronquist, Reference Huelsenbeck and Ronquist2001) in the Cipres Science Gateway (Miller et al., Reference Miller, Pfeiffer and Schwartz2010) with the GTR + G nucleotide substitution model (Rodriguez et al., Reference Rodriguez, Oliver, Marin and Medina1990) proposed by the Akaike information criterion (AIC). For BI, two analyses were run based on codon partition models: (1) partitioned (MBP), i.e., each codon position was treated separately, as they are subject to different evolutionary rates, and (2) non-partitioned (MBUP). The settings for the Bayesian Markov chain Monte Carlo (MCMC) process for the non-partitioned dataset A and partitioned datasets A and B included two runs of 10 × 106 MCMC generations (×4 chains) with a sampling frequency of 1000 generations. For the non-partitioned dataset B, the same settings were applied except that the number of MCMC generations, i.e., 15 × 106, was increased in order to diminish the autocorrelation. The relative burn-in was 10%. MCMC results were checked with the program Tracer 1.6.0 and all trees were displayed using FigTree 1.4.0 (both available at http://tree.bio.ed.ac.uk/software/).
Median-joining (MJ) network
The software NETWORK 4.6.1.2 (Bandelt et al., Reference Bandelt, Forster and Rohl1999) (http://www.fluxus-engineering.com/sharenet.htm) was applied to construct the MJ networks aiming to multistate the characters. For MJ reconstructions the mitochondrial haplotypes of outgroups were excluded.
Results
The proportion of gaps and completely undetermined characters in the alignment as generated in RA × ML was 0.09 and 0.06% for datasets A and B, respectively; the relevant distinct alignment patterns were 250 and 274. The genetic polymorphism of datasets A and B is shown in table 1.
Results generated for dataset A (i.e., forward sequencing of the COI-3′ region) and dataset B (i.e., bidirectional sequencing) in DNaSP 5.10.01, after excluding the outgroups.
Molecular analyses vs. morphological delimitation
Initial assignment to species level was based on morphology, with the 75 study specimens classified into 28 taxa. PTP models included the Eumerus taxa plus the four outgroups and predicted 32–49 taxa for dataset A and 31–46 taxa for dataset B. MP (figs 1 and 2), ML (figs S1 and S2) and BI analyses (figs 3, 4, S3 and S4) yielded similar tree topologies for both A and B datasets, with two main clusters and the nodes of the putative taxa strongly supported. Bootstrap values in MP and ML trees were generally low for both datasets, whereas posterior probability values were much higher for BI trees.
Maximum parsimony analysis for dataset A produced 72 equally parsimonious trees; the strict consensus tree is illustrated here. Length 989 steps, Consistency index (CI) = 33, Retention index (RI) = 71; filled circles denote unique changes, open circles non-unique. Bootstrap support values (>50) are illustrated above the branches.
Maximum parsimony analysis for dataset B produced 30 equally parsimonious trees; the strict consensus tree is illustrated here. Length 1153 steps, Consistency index (CI) = 32, Retention index (RI) = 71; filled circles denote unique changes, open circles non-unique. Bootstrap support values (>50) are illustrated above the branches.
Bayesian analysis of the dataset A (partitioned data). Values indicate Bayesian probability.
Bayesian analysis of the dataset B (partitioned data). Values indicate Bayesian probability.
In all trees, i.e., those generated by MP, ML and BI, the two major clades comprised the same groups of taxa in each dataset (table 2; figs 1–4 and S1–S4). One clade consisted of taxa belonging to the tricolor group (six taxa) with the second comprising the remaining 22 taxa. Our molecular-derived topologies revealed seven different groups of taxa, with all but eight of the 28 taxa assignable to one of these groups. We named these seven groups basalis, minotaurus, ornatus, pulchellus, strigatus, sulcitibius and tricolor (for more details see table 2). The eight ‘ungrouped’ taxa were E. torsicus Grković et Vujić, 2015, E. aff. rubiginosus, E. aff. tarsatus, E. alpinus Rondani, 1857, E. argenticornis Lyneborg, in litt., E. clavatus Becker, 1923, E. hungaricus Szilady, 1940, and E. tarsatus Lyneborg, in litt. (table 2). Slight differences were apparent in the topologies generated by datasets A and B. In dataset B, E. argenticornis, E. tarsatus, E. aff. rubiginosus and E. aff. tarsatus clustered together, close to the ornatus group (see figs 2, 4, S2 and S4), but in dataset A they appeared altogether, dispersed or clustered differently (see figs 1, 3, S1 and S3). Other observed discrepancies were as follows: (a) E. clavatus was in the sulcitibius group in all but the MP analysis of dataset B (fig. 2); (b) the position of E. alpinus differed in each analysis, and (c) branching topology for the tricolor group was only similar in MP, MBP and MBUP analysis of dataset A (figs 1, 3 and S3) and ML and MBUP analyses of dataset B (figs S2 and S4). Branching topology for the ornatus group was consistent across all analyses.
Taxa grouping as formed after the implementation of the MP, ML and BI analyses and the MJ network reconstructions.
The groups are in concordance with the taxa morphology. Taxa groups are separated with ≥19 mt steps.
Among the 28 taxa, three were previously undescribed, but were closely related to known taxa and, thus, were named E. aff. barbarus (collected in Morocco), E. aff. rubiginosus and E. aff. tarsatus (both collected in South Africa). Both molecular (high nodal support) and morphological (clear diagnostic features) aspects strongly supported the species delimitation of these three taxa. The E. aff. barbarus was included in the sulcitibius group, whereas the other two taxa were not collapsed to any of the seven supported groups (table 2).
We noted artifacts of Long Branch Attraction (LBA) artifacts in the ML and BI analyses of both datasets, but not in MP. LBA was due to sequences EU106, EU108, EU109, EU111, EU115 and EU117 (figs 1–4 and S1–S4). MJ network reconstructions supported our other analyses with sequences from both datasets grouping similarly to the clusters present in MP, ML and BI trees. Furthermore, the number of mutational steps between haplotypes was consistent with our phylogenetic reconstructions (figs 5 and 6).
The median-joining network of the haplotypes of the dataset A, as it was constructed by Network software ver. 4.6.1.2 (http://www.fluxus-engineering.com). Circle sizes are proportional to haplotype frequencies. The number of mutational steps is the one between each pair of haplotypes. When not stated, one mutational step interferes between the nodes/OTUs. Taxa and molecular taxa groups (supported by morphology) are depicted (apart of basalis and sulcitibius group).
The median-joining network of the haplotypes of the dataset B, as it was constructed by Network software ver. 4.6.1.2 (http://www.fluxus-engineering.com). Circle sizes are proportional to haplotype frequencies. The number of mutational steps is the one between each pair of haplotypes. When not stated, one mutational step interferes between the nodes/OTUs. Taxa and molecular taxa groups (supported by morphology) are depicted (apart of basalis group).
Intra- and interspecific variability
In both datasets A and B, we recorded a high number of mtDNA COI haplotypes and very rich haplotype (Hd > 0.95) and nucleotide diversity (Pi > 0.005) (table 3). No shared haplotypes between delimitated species were obtained (table 3).
Results generated for dataset A (i.e., forward sequencing of the COI-3′ region) and dataset B (i.e., bidirectional sequencing) in DNaSP 5.10.01 for E. amoenus (18 sequences), E. pulchellus (seven sequences) and E. pusillus (11 sequences).
The basic molecular diversity indices for E. amoenus, E. pulchellus and E. pusillus were calculated and are shown in table 3. For E. amoenus, the MJ network for dataset A showed one to nine mutational steps among haplotypes detected in the Lesvos population, and one to three mutational steps among haplotypes in the Samos population (fig. 5). The MJ network constructed using dataset B revealed one to eight mutational steps among E. amoenus haplotypes from Lesvos and one to four among haplotypes found in the Samos population (fig. 6). When selecting one specimen of E. pulchellus from each of seven geographical origins the MJ network analysis of dataset A showed one to six mutational steps, and dataset B one to seven mutational steps (figs 5 and 6). For E. pusillus, two or more samples from four different geographic origins were analyzed and MJ network reconstructions showed one mutational step for both datasets among different E. pusillus haplotypes (figs 5 and 6).
Discussion
A COI gene-based system has been successfully employed for species delimitation in various hoverfly genera such as the ruficornis group of the genus Merodon (Milankov et al., Reference Milankov, Ståhls and Vujić2008b ; Vujić et al., Reference Vujić, Radenković, Ståhls, Ačanski, Stefanović, Veselić, Andrić and Hayat2012), the Cheilosia vernalis complex (Ståhls et al., Reference Ståhls, Vujić and Milankov2008), the genus Chrysotoxum (Suk & Han, Reference Suk and Han2013; Nedeljković et al., Reference Nedeljković, Ačanski, Dan, Obreht-Vidaković, Ricarte and Vujić2015), the Afrotropical hoverflies (Jordaens et al., Reference Jordaens, Goergen, Virgilio, Backeljau, Vokaer and De Meyer2015) and the genus Platycheirus (Young et al., Reference Young, Marshall and Skevington2016). Our study is the first implementation of molecular tools to infer species delimitation in the genus Eumerus. We assessed the feasibility of using COI-3′ fragment for Eumerus taxonomic inference on 75 specimens assigned to 28 putative ‘taxonomical’ taxa clusters. Various tree inference approaches on genetic data conformed to morphological species assignment. Since species delimitation (through conventional classical taxonomy) of Eumerus has proven challenging in the past, generation of barcode sequences to diagnose species within this genus can prove very useful. A DNA barcode library for Eumerus is currently under construction, with more than 130 sequences having been generated in the last 6 months. The present study contributes to enriching accessible barcode records; an assessment of GenBank records on 11 March 2016 revealed that this study has provided more than 50% of the available Eumerus sequences and species to date. In addition, our study extends representation of both the number and geographic distribution of Eumerus species.
We generated two datasets that differed in terms of sequence length in order to test whether longer sequences improved taxonomic resolution by support values for nodes. Dataset B was 99 nucleotides longer and possessed 28 more parsimony-informative sites compared to dataset A. However, dataset A still provided high taxonomic resolution, confirming the efficacy of a COI system based on approx. 650 bp for species delimitation. Due to sample unavailability arising from fieldwork limitations (and the absence of available sequences in barcoding databases), few sequences – in some cases even only one – were obtained for some taxa. Ahrens et al. (Reference Ahrens, Fujisawa, Krammer, Eberle, Fabrizi and Vogler2016) discussed the issue of the singletons’ (the only representative sequence of a species) issue in DNA-based species delimitation studies and ascertained that ‘a high proportion of singletons has little impact on the accuracy of inferred species limits, and thus rarity (and singletons) should not be conflated with the much more pertinent population genetics parameters’.
The outcomes of our analyses were congruent for both datasets and indicated that the genus Eumerus is divided into two main lineages: the tricolor group and a lineage consisting of all the other taxa (both grouped and the ‘ungrouped’). Tree topologies differed slightly within and between dataset(s). Although longer sequences improved phylogenetic resolution, they did not fully resolve the position of some taxa, e.g., E. argenticornis, E. tarsatus, E. aff. rubiginosus and E. aff. tarsatus. As singletons were used in a few cases, employment of replicate reference specimens could be beneficial in determining the phylogenetic positions of unresolved Eumerus taxa. Our is only the second study to date presenting a hypothetical Eumerus phylogeny; Doczkal & Pape (Reference Doczkal and Pape2009) found indications of the genus being paraphyletic based on morphological characters (although this was not corroborated nor further investigated by them). Even if it was not the purpose of this study to deliberate the genus’ phylogeny, based on our results Eumerus could be monophyletic with two main lineages. In addition, the formation of several groups within the genus reveals certain affinities among the species. Such affinities have never been discussed before, except in Speight (Reference Speight2014) who only commented on the strigatus group that included the taxa E. consimilis Šimić et Vujić, 1996, E. narcissi Smith, 1928, E. ruficornis Meigen, 1822, E. sogdianus Stackelberg, 1952, and E. strigatus (Fallen), 1817. Based on our data, strigatus group is composed of at least two taxa (E. amoenus and E. consimilis) but we did not have sequences for the other taxa considered by Speight (Reference Speight2014) to incorporate them into our analyses. However, our findings do suggest that E. amoenus has a probable (genetic) affinity with the other taxa of Speight's (Reference Speight2014) strigatus group, but this needs further verification.
The large number of mitochondrial haplotypes we generated underlines the high genetic diversity of Eumerus, despite the employment of only one genomic region. Species did not exhibit shared haplotypes and were well separated as shown by the number of mutational steps in the MJ networks. This is not always the case for hoverflies (e.g., the genus Melanostoma Schiner, 1860, Haarto & Stahls, Reference Haarto and Ståhls2014). Therefore, in Eumerus, a COI gene-based system can yield unique, distinguishing mtDNA haplotypes between species.
We also studied intraspecific variation for E. amoenus, E. pusillus and E. pulchellus for which we had more abundant sequences and geographic spread. These taxa are considered ‘common’, i.e., they are widely distributed in all the Mediterranean peninsulas (Anatolian, Apennine, Balkan and Iberian). Eumerus amoenus, E. pulchellus and E. pusillus were analyzed further, for intraspecific genetic diversity. Eumerus pulchellus exhibited the highest haplotype diversity, with each single specimen per geographic origin having a different haplotype. Genetic distances among haplotypes inferred from the number of mutational steps and geographic distances between sampling locations (Greece, Italy and Montenegro), showed no clear pattern. Eumerus amoenus also presented very high haplotypic diversity, which was to be expected because sampling was performed over a wider geographical area. MJ network analysis of haplotypes generated a star-like pattern for this latter taxon, suggestive of past expansion (Bandelt et al., Reference Bandelt, Forster, Sykes and Richards1995). Nucleotide diversity values were similar for E. pulchellus and E. amoenus for dataset A, whereas they appeared to be a bit higher for E. amoenus and E. pusillus based on dataset B. Eumerus pulchellus and E. pusillus were represented by more than two specimens from neighboring and remote geographical areas, and in general showed the lowest values for molecular diversity indices among the three analyzed taxa. The star-like patterns observed for MJ haplotype networks for both datasets for E. pusillus may indicate a recent expansion of this taxon, resulting in the observed lower genetic diversity. Eumerus pusillus specimens originating from Crete, Karpathos and Naxos (one sequence, EU66) shared the same haplotypes in datasets A and B. Specimens from Chios and others from Naxos presented different haplotypes in both datasets. Further conclusions from these discrepancies are limited due to the low sample sizes for the population genetic analyses and, in order not to be speculative, we encourage further intraspecific analyses.
LBA is a sensitive issue in phylogenetic analyses. We used maximum parsimony and likelihood based methods (including BI), with these latter having been proven to be less sensitive to LBA-artifacts compared with maximum parsimony (Bergsten, Reference Bergsten2005). Nevertheless, long branches appeared in both ML- and BI-derived phylogenetic trees, but not for MP, with datasets A and B sharing the same LBAs. LBA-artifacts can arise due to several factors, such as poor taxon sampling and selection of highly divergent outgroups (for more details, see the review by Bergsten, Reference Bergsten2005). Here, our MJ network reconstructions for each dataset explain the LBAs as those taxa having the highest numbers of mutational steps. In addition, the geographic origins of samples should be taken into account. It is important to clarify that the ornatus clade included taxa from very geographically distant areas, i.e., Greece (Dadia) and South Africa (KwaZulu-Natal), so LBAs were to be expected. Given that long branched taxa occurred in all trees, we claim that our topologies remain robust. Adding more samples for LBA taxa could lessen the impact of the LBAs. However, taxon sampling can remain an issue for molecular and other data. We felt that inclusion of as many taxa as possible in our analyses was paramount and so we chose not to exclude any taxa.
To conclude, the present study reveals the adequacy of the COI gene fragment to delimit species in the genus Eumerus. Forward and bidirectional sequencing datasets led to similar results; the forward sequencing dataset appeared to be sufficient to identify taxa within Eumerus and the more enriched sequence dataset, i.e., bidirectional, provided slightly more information and mitigated certain (though not all) analytical problems. Moreover, we reveal high intraspecific diversity and a high number of mitochondrial haplotypes. Our findings confirm the potential of an integrative approach – combined usage of a COI barcoding system and morphological characters – to diagnose and delimit species within the genus Eumerus. For a complete revision of the genus, including phylogenetic inferences, we endorse the usage of additional molecular markers and/or longer mitochondrial sequence (>800 bp). More taxa and more specimens per taxa should also be sought in order to overcome the drawbacks faced in the present study.
Supplementary material
The supplementary material for this article can be found at http://dx.doi.org/10.1017/S0007485316000729.
Acknowledgements
This research was co-financed by the European Union (European Social Fund – ESF) and Greek national funds through the Operational Program ‘Education and Lifelong Learning’ of the National Strategic Reference Framework (NSRF) – Research Funding Program: THALES. Investing in knowledge society through the European Social Fund. We thank John O'Brien for proof reading and making constructive suggestions to the text.
Conflict of interest statement
We certify that there is no conflict of interest regarding the publication of this manuscript.
