DNA sequence analysis suggests that cytb-nd1 PCR-RFLP may not be applicable to sandfly species identification throughout the Mediterranean region

Molecular methods are increasingly used for both species identification of sandflies and assessment of their population structure. In general, they are based on DNA sequence analysis of targets previously amplified by PCR. However, this approach requires access to DNA sequence facilities, and in some circumstances, it is time-consuming. Though DNA sequencing provides the most reliable information, other downstream PCR applications are explored to assist in species identification. Thus, it has been recently proposed that the amplification of a DNA region encompassing partially both the cytochrome-B (cytb) and the NADH dehydrogenase 1 (nd1) genes followed by RFLP analysis with the restriction enzyme Ase I allows the rapid identification of the most prevalent species of phlebotomine sandflies in the Mediterranean region. In order to confirm the suitability of this method, we collected, processed, and molecularly analyzed a total of 155 sandflies belonging to four species including Phlebotomus ariasi, P. papatasi, P. perniciosus, and Sergentomyia minuta from different regions in Spain. This data set was completed with DNA sequences available at the GenBank for species prevalent in the Mediterranean basin and the Middle East. Additionally, DNA sequences from 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), from 19 countries, were added to the data set. Overall, our molecular data revealed that this PCR-RFLP method does not provide a unique and specific profile for each phlebotomine species tested. Intraspecific variability and similar RFLP patterns were frequently observed among the species tested. Our data suggest that this method may not be applicable throughout the Mediterranean region as previously proposed. Other molecular approaches like DNA barcoding or phylogenetic analyses would allow a more precise molecular species identification. Electronic supplementary material The online version of this article (doi:10.1007/s00436-015-4865-5) contains supplementary material, which is available to authorized users.


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).
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.
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).
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 MgCl 2 , 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 t he 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.
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).
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.

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.