Introduction

Hematophagous females of certain species of phlebotomine sandflies (Diptera: Psychodidae: Phlebotominae) are vectors of human and different animal pathogens worldwide (Pigott et al. 2014). They are responsible for the transmission of diseases caused by protozoans (Leishmania), viruses (Phlebovirus, Vesiculovirus, and Orbivirus), and bacteria (Bartonella bacilliformis) (Antoniou et al. 2008; Depaquit et al. 2010; Ready 2013; Maroli et al. 2013). (i) In the Mediterranean basin, different species of sandflies are vectors of Leishmania tropica, L. major, L. infantum, and L. donovani, the causative agents of leishmaniasis; (ii) the main viral diseases transmitted by phlebotomine sandflies in this region are caused by Phlebovirus (family Bunyaviridae), like the sandfly fever (SF) caused by the SF Sicilian virus and the SF Naples virus, as well as the summer meningitis caused by Toscana virus.

Both leishmaniasis and sandfly-borne viral diseases represent a serious threat to public veterinary health and are considered (re-)emerging infections in the Mediterranean area. Therefore, vector surveillance should be among the pillars of the activities devoted to the design and implementation of effective control measures against these diseases (WHO/EMRO 2008; Depaquit et al. 2010; Charrel et al. 2012; Maroli et al. 2013; Antoniou et al. 2013; Ejov and Dagne 2014). Although fragile in nature, sandflies are able to adapt to a broad variety of environmental conditions including urban or peri-urban settings; the recent leishmaniasis outbreak in Madrid (Spain) is an example of this plasticity (Tarallo et al. 2010; Galvez et al. 2011; Carrillo et al. 2013; Ready 2013; Lisi et al. 2014). Their potential to spread to northern Europe cannot be omitted; the northward spread of leishmaniasis in Italy has been well documented, and the potential establishment of vector species in Germany is a cause of concern (Maroli et al. 2008; Biglino et al. 2010; Fischer et al. 2010; Melaun et al. 2014; Medlock et al. 2014). Due to their public health importance, surveillance of sandfly populations is critical to assess both their geographical distribution and that of the diseases they transmit, as well as routes of introduction to non-endemic areas.

As only some species of sandflies are confirmed vectors, an important aspect of surveillance is the appropriate identification of the obtained specimens. Traditionally, the identification of sandflies at species level has relied on the morphological analysis of anatomical structures such as the pharynx, spermathecae and cibarium of the females, and the terminal genitalia of the males (Killick-Kendrick et al. 1991). This approach requires a high degree of expertise and can be challenging when the specimens are not properly preserved. To contribute to this task, different molecular methods have been proposed as complementary approaches for specimen identification. It has been suggested that the amplification of a DNA target encompassing a fragment of the mitochondrial genes cytb and nd1 and the subsequent digestion of the amplicons with the endonuclease Ase I would be a useful tool for the rapid identification of the most common phlebotomine sandflies in the Mediterranean region (Latrofa et al. 2012). However, it is important to take into account that this study was based on specimens collected from a limited geographical area of Italy (Basilicata). In addition, in a recent report, Bounamous et al. (2014) failed to effectively distinguish some phlebotomine species in Algeria, even belonging to different genera (like P. ariasi and S. schwetzi).

In the present study, we aimed to assess whether the methodology proposed by Latrofa et al. (2012) was also useful for sandfly species identification in specimens from a wider geographical origin, including Spain and other endemic areas, where sandfly-borne diseases are endemic and entomological identification is required for surveillance.

Material and methods

We used two different approaches to test the PCR-RFLP protocol proposed by Latrofa et al. (2012). The first one consisted on a wet-lab PCR-Ase I RFLP analysis on selected representative sequences of each species collected in Spain (Assembly 1: phlebotomine sandflies from Spain (N = 155), see below); the second approach was an in silico analysis of the DNA sequences of the sandflies from assembly 1 plus a collection of 277 DNA sequences retrieved from the GenBank database belonging to 12 different species of sandflies that were collected by other authors in other Mediterranean and Middle East countries (Assembly 2: DNA sequences retrieved from the GenBank (N = 277), see below).

Assembly 1: phlebotomine sandflies from Spain (N = 155)

The specimens of this assembly were captured at different sites from mainland Spain and the Balearic Islands between June and October 2013, at different time-points within the sandfly activity season in our region (May–October). The specimens were collected using CDC miniature light traps (John W. Hock, Company, Gainesville, FL) and stored in 70° ethanol at −20 °C until further analyses. Sampling sites were selected according to the different bioclimatic zones in Spain (Fig. 1).

Fig. 1
figure 1

Geocoded surveyed sites and number of sandflies collected shown on a bioclimatic zone map of Spain. Numbers indicate sampling locations: Agés (Burgos) (1), Escarabajosa de Cabezas (Segovia) (2), Zaorejas (Guadalajara) (3), Majadahonda (Madrid) (4), Vicálvaro (Madrid) (5), Cañamero (Cáceres) (6), Buenache de Alarcón (Cuenca) (7), Castilblanco de los Arroyos (Sevilla) (8), Achivel (Murcia) (9), Moratalla (Murcia) (10), and San Joan de Labritja (Ibiza) (11). Colored filled circles indicate phlebotomine species. Circle size indicates the number of specimens according to the scale shown

Both male and female were morphologically identified based on their genitalia. The pharyngeal armature and cibarium were also studied in the female specimens. The tip of the abdomen (between segments VI and VII) and head of the females were cutoff and cleared in Mark André medium (Abbonnenc 1972). The specimens were mounted onto glass slides in Hoyer medium (Upton 1993) and identified at species level, according to the taxonomic keys proposed by Gil et al. (1989). Mounted parts of the specimens were kept for future reference. As morphological identification was based on the study of the male genitalia and female head and spermatechae, the remaining of the body was kept at −20 °C for DNA extraction and downstream molecular analyses. A detailed description of sampling locations, species (P. perniciosus, P. ariasi, P. papatasi, and S. minuta), and number of specimens included in assembly 1 is presented in Fig. 1 and Table 1.

Table 1 Description of sampling locations, species, and number of specimens included in assembly 1 (N = 155)

Assembly 2: DNA sequences retrieved from the GenBank (N = 277)

In order to increase our study sample and to complete the picture of representative sandflies in the Mediterranean region beyond the specimens captured in Spain (assembly 1), we retrieved from GenBank a total of 277 DNA sequences corresponding to 13 different phlebotomine species (P. ariasi, P. balcanicus, P. caucasicus, P. chabaudi, P. chadlii, P. longicuspis, P. neglectus, P. papatasi, P. perfiliewi, P. perniciosus, P. riouxi, P. sergenti, and S. minuta) containing the cytb-nd1 sequence encompassed by the primers PhleF and PhleR described by Latrofa et al. (2011). These sequences were identified by BLASTn search using as query different DNA sequences obtained from specimens in assembly 1. Alternatively, taxon-specific key (Phlebotomus, Sergentomyia) and gene identifier (cytb, cyt b, cyt-b, cytochrome b, cytochrome-b) were also used to search for candidate sequences. Selected DNA sequences included those from species endemic to the Mediterranean basin (including Italy and Algeria), the Middle East (Afghanistan and Iran), and other regions (Portugal) where the former species are also present. Details on the specimens included in assembly 2 are presented in Table 2. Further details related to these sequences are presented as supplementary material (Supplementary material 1).

Table 2 Description of the sampling location, species, and number of specimens included in assembly 2 (N = 277), for which cytb-nd1 sequences were retrieved from the GenBank

DNA extraction and cytb-nd1 PCR from assembly 1

Genomic DNA was independently extracted from each specimen of assembly 1. Specimens were grinded and processed with the QIAamp DNA Mini Kit (QIAgen, Germany) in accordance with the manufacturer’s instructions. DNA was eluted in 50 μL of PCR-grade water and stored at −20 °C.

The mitochondrial DNA target encompassing the cytb-nd1 region was amplified according to Latrofa et al. (2012) with minor modifications. Two microliters of DNA was used in a 25-μL final volume PCR reaction, including standard reaction buffer 2 mM MgCl2, 0.2 mM each dNTP, and 0.7 U of Thermus sp. DNA polymerase (Biotools, B&M Labs, Spain), and 15 pmol of each primer PhleF (5′-AAT AAA TTA GGA GGA GTA ATT GC-3′) and PhleR (5′-GCC TCG AWT TCG WTT ATG ATA AAT T-3′) (Sigma-Genosys).

Amplification was performed on a 9800 Fast Thermal Cycler (Applied Biosystems), standard ramp enabled, with the following conditions: initial denaturation of 5 min at 94 °C, 40 cycles at 94 °C for 1 min, 52 °C for 1 min, and 72 °C for 1 min, and final extension at 72 °C for 10 min. All the amplified products were resolved on 2 % agarose gels stained with Pronasafe Nucleic Acid Staining (Laboratorios CONDA, Spain) and visualized under UV light.

Ase I RFLP analysis

Selected representative amplicons obtained at the cytb-nd1-PCRs from assembly 1 were digested with the endonuclease Ase I as described by Latrofa et al. (2012). Digested products were resolved on 2 % agarose gels as described above and their sizes estimated by comparison to a 100-bp DNA ladder (Nippon Genetics, Europe GmbH).

DNA sequencing

Both strands of the cytb-nd1 PCR products from assembly 1 were sequenced in duplicate using the same PCR primer set. The Big-Dye Terminator Cycle Sequencing Ready Reaction Kit V3.1 and the automated ABI PRISM 377 DNA sequencer (Applied Biosystems, Foster City, CA) were used. Sequences obtained were analyzed and edited using the software Molecular Evolutionary Genetics Analysis (MEGA) version 5.2 (Tamura et al. 2011).

Alignment and trimming of the DNA sequences from assembly 1 and assembly 2

DNA sequences generated from the 155 sandfly specimens collected in Spain (assembly 1) and the 277 DNA sequences obtained from the GenBank (assembly 2) were loaded in MEGA 5.2 and aligned using the ClustalW algorithm (Thompson et al. 1994). The sequences were trimmed to the cytb-nd1 sequence flanked by the primers PhleF and PhleR. A P. perniciosus DNA sequence retrieved from the GenBank (Accession Number JF766956) was used as reference sequence to guide both alignment and trimming. Minor editing to complete the 3′-end was conducted in some cases, provided all other specimens of the same species had identical sequence. Primer sequences were added to the end of each sequence to simulate the PCR product.

In silico Ase I RFLP

A virtual Ase I RFLP was performed on all trimmed sequences flanked by the PhleF/PhleR primers from assemblies 1 and 2 using the program NEBcutter V2.0 from NEW ENGLAND BioLabs® Inc. (Vincze et al. 2003). which allows to produce theoretical digests with restriction enzymes. Settings for the enzyme Ase I included a virtual run of the digested products on a 2 % agarose gel and the full resolution of a 100-bp ladder for DNA sizing, L = 70 mm.

Results

The sequences corresponding to a total of 432 sandfly specimens (assemblies 1 and 2) were tested by the cytb-nd1-PCR Ase I-RFLP method in silico. In addition, a subset of sequences of assembly 1 was also experimentally assessed. Overall, this study covered 13 sandfly species from five different subgenera: 249 Larroussius specimens (152 P. perniciosus, 6 P. longicuspis, 57 P. ariasi, 26 P. perfiliewi, 2 P. chadlii, 6 P. neglectus), 33 Phlebotomus specimens (P. papatasi), 127 Paraphlebotomus specimens (57 P. sergenti, 13 P. riouxi, 32 P. chabaudi, 25 P. caucasicus), 2 Adlerius specimens (P. balcanicus), and 21 Sergentomyia specimens (S. minuta).

cytb-nd1-PCR and experimental Ase I RFLP on specimens from assembly 1

DNA from all specimens in assembly 1 was successfully amplified using the cytb-nd1-PCR method. Obtained amplicon sizes ranged from 472 to 482 bp depending on the phlebotomine species considered (Online Resource 2). Ase I RFLP analysis was conducted in a limited number of DNA isolates corresponding to representative haplotypes. In all cases, the digested PCR products showed identical restriction patterns to those predicted by the in silico analysis (Fig. 2a, b). Interestingly, Ase I RFLP profiles for P. ariasi-variant 1 and P. perniciosus-variant 2 produced a practically identical pattern comprising bands of 26 and c455 bp. On the other hand, P. perniciosus variant 1 returned a distinct pattern with two bands of approximately 100 and 350 bp. The P. papatasi specimen yielded a band of c375 bp, whereas S. minuta generated a c230-bp band (Fig. 2 and Table 3).

Fig. 2
figure 2

Ase I RFLP patterns of the cytb-nd1 PCR products on selected representative sequences of each species collected in Spain (assembly 1). A 100-bp ladder (lane M), P. perniciosus variant 1 (lane 1), P. perniciosus variant 2 (lane 2), P. ariasi variant 1 (lane 3), P. papatasi (lane 4), and S. minuta variant-2 (lane 5). a Wet-lab PCR-Ase I RFLP analysis and b in silico Ase I RFLP pattern based on the cytb-nd1 DNA sequence of the specimens assessed in a using NEBcutter V2.0 software

Table 3 Summary of patterns obtained by AseI RFLP analysis with the 432 DNA sequences considered in this study

Alignment and trimming of cytb-nd1 DNA sequences from specimens included in assemblies 1 and 2

A final alignment was built with the 155 DNA sequences generated from the assemblies 1 and 2. The alignment (primer sequences included) encompasses the nucleotides at positions 440 to 920 of the P. perniciosus sequence used as reference (GenBank Acc no. JF766956) and comprises cytochrome b gene (partial cds), tRNA-Ser (complete sequence), and NADH dehydrogenase subunit 1 gene (partial cds). DNA sequence analysis showed a different PCR product size for each species, ranging from 472 to 491 bp. This is mainly due to the sequence length variability of the intergenic mitochondrial DNA spacer-1 (Igs1), lying between cytb gene and tRNA-Ser (Supplementary material 2).

The 155 DNA sequences from assembly 1 were deposited in GenBank under the following accession numbers: P. perniciosus KP685413-KP685538, P. ariasi KP685539-KP685550, P. papatasi KP702248 and KP702249, and S. minuta KP702250-KP702264.

In silico assessment of Ase I RFLP profiles on cytb-nd1 DNA sequences obtained from assemblies 1 and 2

The theoretical in silico Ase I RFLP analysis performed on the 432 cytb-nd1 DNA sequences, comprising 13 different sandfly species, returned 16 (I–XVI) different RFLP profiles. A virtual image of the 16 RFLP patterns is given in Supplementary material 3. Specific details for each specimen are provided in online resource 1. These results were identical to those obtained by experimental Ase I RFLP on the selected subset of sandfly specimens from assembly 1 (Fig. 2a, b).

Single, but not necessarily distinctive, RFLP patterns were obtained in all specimens from the following species: P. papatasi (pattern VI), P. neglectus (IX), P. perfiliewi (X), P. riouxi (VII), P. chadlii (XIII), and P. balcanicus (VII). Among all species, only P. papatasi gave a unique species-specific pattern that could be unequivocally differentiated. In addition, P. riouxi and P. balcanicus shared a single RFLP pattern (VII) that was also identical to the variant 1 observed in 31 out of 32 specimens of P. chabaudi. The remaining species presented two distinctive patterns each: P. perniciosus (151 specimens assigned to pattern I, and the remaining one to pattern II), P. ariasi (56 pattern IV, one pattern V), P. longicuspis (5 pattern III, one pattern I), P. chabaudi (31 pattern VII, one pattern VIII), P. sergenti (52 pattern XI, 5 pattern XII), P. caucasicus (24 pattern XIV, one pattern XI), and S. minuta (16 pattern XV, 5 pattern XVI).

Seven species could not be accurately differentiated due to shared RFLP patterns. Thus, pattern I was common to 151 P. perniciosus and one P. longicuspis, pattern VII was common to all P. riouxi, P. balcanicus, and 31 out of 32 P. chabaudi, and pattern XI was common to one P. caucasicus and 52 out of 57 P. sergenti. In addition to this, a number of in silico RFLP patterns were very similar, but not identical. These included pattern III (26, 90, 104, 261 bp) obtained in 5/7 P. longicuspis and pattern X (26, 90, 104, 260 bp) from all P. perfiliewi analyzed, and also patterns II (26, 455 bp) in 1/152 P. perniciosus, IV (26, 454 bp) in 56/57 P. ariasi, VII (26, 446 bp) in 31/32 P. chabaudi, 13/13 P. riouxi, and 2/2 P. balcanicus, IX (26, 465 bp) in 6/6 P. neglectus, XI (26, 450 bp) in 52/57 P. sergenti and 1/25 P. caucasicus, and XIII (26, 451 bp) in 2/2 P. chadlii (Fig. 3).

Fig. 3
figure 3

In silico Ase I RFLP analysis obtained with NEBcutter V2.0 showing RFLP patterns that are similar or almost undistinguishable. Patterns are indicated in Roman numerals at the bottom of the image. The virtual run shows the digested products on a theoretical 2 % agarose gel and the full resolution of a 100-bp ladder. P. longicuspis_1 (lane 1), P. perfiliewi (lane 2), P. pernicious_2 (lane 3), P. ariasi_1 (lane 4), P. chabaudi_1 (lane 5), P. riouxi (lane 6), P. balcanicus (lane 7), P. neglectus (lane 8), P. sergenti_1 (lane 9), P. caucasicus_1 (lane 10), and P. chadlii (lane 11)

Discussion

Sequence analysis of a cytb gene region from 13 different sandfly species provides in silico and experimental evidence demonstrating that the cytb-nd1 PCR-RFLP assay may not be applicable to phlebotomine species identification throughout the Mediterranean region as initially proposed by other authors. Our data confirm that it is not possible to obtain a unique, distinctive RFLP pattern for each individual sandfly species tested excluding P. papatasi and S. minuta. We found a significant proportion of species returning identical or almost undistinguishable RFLP profiles. This fact would represent a pitfall for vector incrimination in regions where two coendemic species share the same or almost identical RFLP pattern. This is the case for P. perniciosus pattern II and P. ariasi pattern IV, both identified in Spanish specimens, being both vector of L. infantum. Importantly, this lack of discriminatory power has been previously reported using the same methodology by Bounamous et al. (2014). In that survey, the authors were unable to effectively discriminate among 17 sandfly species present in Algeria, even when a double digestion with Ase I and Mnl I restriction enzymes was used.

It is important to take into consideration that the study by Latrofa et al. (2012) was restricted to sandfly specimens from the south of Italy. In our subset of Italian specimens (n = 6) from the same area, we found similar findings except for S. minuta. Therefore, pattern XVI (the one showed by Latrofa et al. 2012) was identified in 5 S. minuta specimens, whereas the remaining one (Gen Bank Acc no. JF766981) presented pattern XV, as in all 15 S. minuta specimens from Spanish origin.

The similarity of certain RFLP patterns would be particularly problematic for phlebotomine species identification if the patterns correspond to sympatric species. This seems to be the case for Algeria and Tunisia, where P. longicuspis and P. perfiliewi are coendemic and thought to be potential vectors of L. infantum (Ready 2013). Discrimination between both species in these countries would be very difficult to achieve if we obtain patterns III (as in five out of six P. longicuspis in this work) or X (in all 26 P. perfiliewi studied here). In this study, one Moroccan P. longicuspis specimen showed PCR-RFLP pattern I, as did the majority of the P. perniciosus specimens analyzed. These two species coexist in several countries and are remarkably similar morphologically. In Spain, where only P. longicuspis male specimens have been described so far, the presence of P. longicuspis and P. perniciosus has been questioned (Guernaoui et al. 2005).

On the other hand, the similarity found among patterns II (in 1/152 P. perniciosus), IV (in 56/57 P. ariasi), VII (in 31/32 P. chabaudi, P. riouxi, and P. balcanicus), and IX (P. neglectus) would impair the definitive identification due to geographical distribution overlapping between P. ariasi and P. perniciosus; this would be particularly important in Italy, where P. ariasi, P. perniciosus, and P. neglectus coexist. Although P. ariasi and P. perniciosus are proven vectors of L. infantum, this is not the case of P. neglectus (Ready 2013). It is important to point out that P. riouxi and P. chabaudi are similar from the morphological and taxonomic point of view. Both species are known to occur sympatrically in some geographical regions including Algeria and Tunisia. According to our sequence analysis data, no significant differences were found at the nucleotide level to conclusively consider P. riouxi and P. chabaudi as independent species, as previously proposed by Bounamous et al. (2014). This notion is also supported by the phylogenetic data obtained with other molecular markers, including cytb and elongation factor 1-α (Tabbabi et al. 2014). Overall, this information provides molecular and evolutionary evidence suggesting that P. riouxi should be considered either as a synonym or as a subspecies of P. chabaudi.

Pattern VII was the only profile shared by two different subgenera (Paraphlebotomus and Adlerius). Thus, it would be useful to analyze the RFLP patterns in other sandfly species belonging to Adlerius subgenus endemic to the Mediterranean basin. Similarity between patterns XI (in 52/57 P sergenti and in 1/25 P. caucasicus) and XIII (P. chadlii) would be also a problem regarding sandfly identification because these three phlebotomine species coexist in some geographical areas (see www.sandflycatalog.org). Whereas P. sergenti is a proven vector of L. tropica (Ready 2013). this potential role has also been suggested for P. caucasicus (Artemiev 1978; Killick-Kendrick 1999). However, P. chadlii has not been described as suitable vector of the disease to date.

An important contribution of this work is the submission to the GenBank of 155 new DNA sequences for thecytb-nd1 region from Spanish specimens, including 126 P. perniciosus sequences to be added to the single sequence already available at GenBank from this country. Additionally, a total of 15 new S. minuta sequences add to the only six previously available sequences, all corresponding to Italian specimens.

Our results demonstrate that neither in Spain nor throughout the Mediterranean region is possible to accurately differentiate phlebotomine sand flies using the cytb-nd1 PCR-RFLP method. Because sequences from Middle East countries were underrepresented in this study, we cannot draw definitively conclusions on the performance of this method with samples from regions other than the Mediterranean basin. Nevertheless, further work at the subregional level is needed to confirm the extent of these findings. Taking into account the sequence variability observed within the cytb-nd1 region, we believe that a phylogenetic analysis or a DNA barcoding approach would allow more precise species identification, as it has been proposed by other authors (Esseghir et al. 1997; Bounamous et al. 2008; Kruger et al. 2011; Latrofa et al. 2011). Nevertheless, before these approaches are fully developed, they should be based on a solid morphological identification of the specimens.