We collected 85 adult specimens of 16 dipteran species of forensic interest from 24 localities in Belgium and three localities in France (Table 1). All Belgian specimens came from forensic cases. Three specimens from three species (Neomyia cornicina, Polietes lardarius and Eudasyphora cyanella) were collected on corpses but are currently not used for the calculation of the PMI. The French specimens of Chrysomya albiceps and Lucilia sericata were not collected on corpses, but were added because of their forensic interest. Morphological species identification was done by two taxonomic experts of Diptera (YB and ED), using five identification keys (D’Assis Fonseca 1968, Beĭ-Bienko 1988, Rozkošný et al. 1997, Gregor et al. 2002, Szpila 2012). Three Fannia specimens (Fanniidae) could not be identified to the species level and were considered as three putative different species. Specimens were deposited as vouchers at the National Institute of Criminalistics and Criminology in Brussels, Belgium (Table 1).

We extracted genomic DNA from one or two legs per specimen using the NucleoSpin Tissue Kit (Macherey-Nagel) and a final elution volume of 70 µl. Fragments of the COI marker were amplified using two primer pairs TY-J-1460/C1-N-2191 and C1-J-2183/TL2-N-3014 (Sperling et al. 1994, Wells and Sperling 1999). The fragment obtained with the first primer pair encompasses the barcode region of ca. 650 bp used for animals (Hebert et al. 2003). The assembly of the fragments obtained with both primer pairs generated a sequence of 1534 bp corresponding to the complete COI gene. Each 25 µl PCR reaction contained final concentrations of 0.2 mM dNTPs, 0.4 μM of each primer, 2.0 mM MgCl₂, 0.5 U of Taq DNA polymerase (Platinum, Invitrogen), 1 × PCR buffer and 2-4 µl DNA template. The thermal cycler program consisted of an initial denaturation step of 4 min at 94 °C, followed by 40 cycles of 30 s at 94 °C, 30 s at 45 °C and 90 s at 72 °C; with a final extension of 7 min at 72 °C. We cleaned PCR products using the NucleoFast96 PCR Kit (Macherey-Nagel) and sequenced them bidirectionally on an ABI 3130 Genetic Analyzer (Applied Biosystems) using the BigDye Terminator Cycle Sequencing Kit v3.1.

We assembled and aligned sequences in SeqScape v2.5 (Applied Biosystems) and confirmed the absence of stop codons using MEGA5 (Tamura et al. 2011). Sequences were deposited in BOLD (BOLD process ID’s are given in Table 1) and GenBank. All different haplotypes were extracted from the aligned sequences using the R package PEGAS (Paradis 2010). We calculated pairwise p-distances (i.e. the proportion of sites at which two sequences differ) and searched for haplotypes that were shared among species.

Haplotypes were then used as queries to search for most similar sequences in two public databases: GenBank (NCBI, National Centre for Biotechnology) and BOLD (the Barcode of Life Data Systems). These most similar sequences will be called “best matches” sensu Meier et al. (2006) in the following. In GenBank, searches were done using MegaBLAST, the Basic Local Alignment Search Tool (BLAST) optimised for highly similar sequences (Zhang et al. 2000, Morgulis et al. 2008). In BOLD, the in-built Identification System (IDS) was applied (Ratnasingham and Hebert 2007) on two different databases: the Public Record Barcode Database (341 580 sequences; 45 368 nominal species and 11 732 interim species, or candidate species that have not been described yet on 24 May 2013) and the Species Level Barcode Records (1 367 662 sequences; 127 679 species and 53 394 interim species on 24 May 2013). The first database comprises the same records as GenBank because both libraries regularly synchronize their published records. In BOLD, this database of public records is a collection of COI records of minimum 500 bp from the published projects of BOLD. The Species Level Barcode Records of BOLD is used by default in IDS. It contains, in addition to the published COI records, early data release of COI records with a species level identification and a minimum sequence length of 500 bp. These early releases contain all information necessary for barcodes (locality and date of sample collection, trace files and sequence information as well as voucher specimen and database identifiers), have passed computerized quality checks of BOLD but might include provisional taxonomic assignments (Hebert et al. 2010).

In total, we applied five search strategies by submitting the barcode sequences to 1) GenBank, 2) the Public Records of BOLD, 3) the Species Level Records of BOLD including early releases, as well as 4) by using the barcode sequences as queries in combination with a keyword, “barcode”, in GenBank, and 5) by submitting COI sequences longer than the barcode fragment (1412–1534 bp) to GenBank. The use of the keyword “barcode” allowed us to filter the GenBank reference sequences and obtain only best matches that are tagged as barcodes, not only in the field “keyword” but also in any field of GenBank records. Longer COI sequences have not been submitted to BOLD because BOLD was developed to accept sequences from the strict barcode region only. In BOLD, IDS returns a list of maximum 99 best matches and provides a species-level identification for best close matches showing less than 1% divergence (Ratnasingham and Hebert 2007). Since BLAST searches are based on approximate alignments (regions of local similarity between sequences), species assignments are usually preferably performed on the basis of local alignments. Hence we verified that the best hits and their percentages of sequence identity obtained from the MegaBLAST searches in GenBank were identical to those (= 1–p-distance) calculated with MEGA5 (Tamura et al. 2011) using local databases downloaded from GenBank and aligned with CLUSTAL W (Thompson et al. 1994). Identifications were made on the basis of the highly similar best matches (> 99% similarity), according to the “best close match” method of Meier et al. (2006). We qualified each best match with a similarity of > 99% as correct if it had the same species name as the query or as incorrect if it had a different species name than the query. In addition, the identification of a query was considered as unambiguous if all best matches with a similarity of > 99% had the same species name. If this was not the case, then the identification was ambiguous. For each identification, we made sure that best close matches included only records properly identified to the species level by excluding the few records with provisional identifications (a code instead of a nominal species name). We also verified whether the alignment of the query with each best match comprised at least 600 bp. When no best match of > 99% similarity was retrieved for a given query, the presence of conspecific and congeneric barcode sequences of > 500 bp was investigated in both public libraries. If present, their divergences (p-distances) with the queries were calculated using MEGA5 (Tamura et al. 2011).

Using the 42 haplotypes of the barcode region as a query yielded the same results in GenBank and in the Public Record Barcode Database of BOLD. Best matches of > 99% similarity were retrieved for 36/42 haplotypes, representing 11 out of 16 species (Table 1). These best matches were either identical (17/36) or differed from the query in less than three substitutions (19/36). We obtained at least one correct best match for each query. However, species identifications were either unambiguous (18 queries, 8 species) or ambiguous (18 queries, 3 species). For two queries, best matches included species of another genus: Musca domestica Linnaeus, 1758was found for Calliphora vicina and Chrysomya megacephala (Fabricius, 1794) for Lucilia ampullacea. In all other cases of ambiguous identification, best matches involved congenerics: Calliphora croceipalpis Jaennicke, 1867 was found for Calliphora vicina, Lucilia cuprina (Wiedemann, 1830) for Lucilia sericata and Lucilia porphyrina (Walker, 1856) for Lucilia ampullacea. Finally, the number of best matches with > 99% similarity varied from one to more than 99 per query (the number of best matches displayed by BOLD is limited to 99). For five species, less than five sequences with a similarity of > 99% were retrieved (Figure 1).

For six queries, the best matching similarities were < 93.5%. These included the haplotypes of Fannia sp1, sp2 and sp3, Muscina prolapsa and Neomyia cornicina. There were no COI sequences of Muscina prolapsa or of Neomyia cornicina in GenBank. For Fannia, fragments of the barcode region of > 500 bp were available for 14 specimens representing four species, viz. Fannia canicularis (Linnaeus, 1761), Fannia scalaris (Fabricius, 1794), Fannia brevicauda Chillcott, 1961 and Fannia serena (Fallen, 1825) but their p-distances with our three Fannia haplotypes ranged from 6.6% to 16.2%.

Using the Species Level Barcode Records dataset of BOLD (Table 1), highly similar best matches (> 99%) were retrieved for 40/42 queries (14/16 species). Correct best matches were retrieved for all specimens identified at the species level, but identifications were often ambiguous (25 queries, 6 species). This method yielded a higher proportion of best matches of > 99% similarity than when the search was restricted to public records (95% of the queries instead of 86%). However, the proportion of unambiguous identifications was smaller (38% instead of 50% of the queries; Table 2). Yet, in contrast to all the other searches, early-released sequences provided two correct matches for Muscina prolapsa, one match for Fannia sp2 and three matches for Fannia sp1 (correct at the genus level). The latter identification was ambiguous since two best matches showed 100% similarity with Fannia lustrator (Harris, 1780) and one showed 99.85% similarity with Fannia pallitibia (Rondani, 1866). The two queries for which best matches were of < 99% similarity were from Fannia sp3 and Neomyia cornicina. No barcodes were available for Neomyia cornicina in BOLD.

When both the barcode sequences and the keyword “barcode” were used as queries in GenBank, we retrieved best matches of > 99% similarity for Calliphora vicina, Calliphora vomitoria, Lucilia sericata, Protophormia terraenovae and Muscina levida (Tables 1 and 2). All best matches of > 99% similarity were correct and provided unambiguous identifications except for Lucilia sericata, which matched with both correct and incorrect species names (Lucilia cuprina and Lucilia sericata).

Haplotypes of longer COI fragments (1412–1534 bp) were also submitted to a MegaBLAST search on GenBank. Best matches of > 99% similarity were obtained for all haplotypes of Calliphora vicina, Calliphora vomitoria, Chrysomya albiceps, Lucilia ampullacea and Lucilia sericata. Like in the previous analysis, all best matches were correct and provided unambiguous identifications except for Lucilia sericata (best matches included Lucilia sericata and Lucilia cuprina).

With this study we contributed to the establishment of a local COI reference library for fly species of forensic importance in Belgium and France. As such, we provide the first barcodes for Muscina prolapsa and Neomyia cornicina. We also extended the geographic coverage of barcodes of species which hitherto were only sampled from a limited number of localities, e.g. Cynomya mortuorum and Polietes lardarius were each represented by only one barcode sequence from the UK (Kutty et al. 2008). Similarly, barcodes of Muscina levida were until now only available for samples from Canada, Germany (Renaud et al. 2012) and the USA (Nakano and Honda, unpublished). Conversely, barcodes of the other species sampled here were obtained from no more than five European countries. Ideally, a reliable reference library should comprise a large sampling of sequences, not only representing the European dipteran species that are currently used in forensics (whose development times have been studied under different temperature conditions), but also those of potential forensic interest (occurring on carcasses but whose biology has been less studied) and all their close relatives. Currently, 13 species belonging to 10 genera are being used in forensic investigations (Marchenko 2001, Grassberger et al. 2002, Richards et al. 2009, Velásquez et al. 2013). Hence, the geographic coverage of GenBank and BOLD is still far from comprehensive. Yet, we did not observe intraspecific COI divergences of > 1% at COI, neither among specimens sequenced in this study nor between them and their conspecific best matches in the public libraries (intraspecific distances among GenBank sequences were not calculated here). This indicates that geographic coverage does not always have to be complete to allow correct species identification. Nonetheless, a more comprehensive reference library may comprise more haplotypes, allowing a better assessment of the risk of incorrect identifications (Meier et al. 2006). Indeed, an increased sampling can result in a more difficult distinction between some closely related species (Bergsten et al. 2012) and this has considerable importance for courts.

For 86% of the barcode fragments used as queries, we retrieved highly similar conspecific sequences (> 99% similarity) from GenBank and BOLD. The more divergent best matches (< 99% similarity) obtained for the remaining 14% of the queries would have produced either incorrect (Muscina prolapsa and Neomyia cornicina) or doubtful identifications (Fannia) if all best matches were taken into account for identification. The better performance of the best close match method compared to the simple best match method has already been reported (e.g. Meier et al. 2006, Virgilio et al. 2010). However, even with the best close match method, our results revealed three issues that can hamper the DNA-based identification of forensically important flies in Belgium and France using GenBank or BOLD: These databases 1) do not include some fly species of forensic interest, 2) include sequences from misidentified specimens and 3) cannot always discriminate between closely related species. Below, we discuss these three issues in more detail.

Our results showed that some fly species collected at Belgian crime scenes are not represented by COI records in GenBank and BOLD. Muscina prolapsa, for which no barcode sequence is present in GenBank, colonises carrion and buried remains (Gunn and Bird 2011, Prado e Castro et al. 2012). Also, the identification of Fannia species of forensic interest (Prado e Castro et al. 2012) is hampered by their limited representation in GenBank and BOLD. Neomyia cornicina is currently not used for PMI estimation but the availability of reference sequences of such species collected on crime scenes can decrease the risk of incorrect identification and help to characterize the entomofauna surrounding the crime scene (Amendt et al. 2007).

Identifications based on the barcode fragment were ambiguous for 50% of the queries and for 27% of the species. Some ambiguous identifications can result from misidentified sequences in the libraries and could be corrected after re-examining the voucher specimens (Collins and Cruickshank 2013). In our study, the best matches with sequences from different genera could be the result of misidentifications: records of Musca domestica (GenBank accession number JQ350716) and Chrysomya megacephala (KC135926) matched our sequences of Calliphora vicina and Lucilia ampullacea, respectively.

Still, most ambiguous identifications involved closely related species that are not necessarily incorrectly identified (Stevens et al. 2002, Sonet et al. 2012). For example, Wells et al. (2007) and Wells and Stevens (2008) showed that the barcodes of several specimens of Lucilia cuprina (from Hawaii and Asia) are more similar to those of Lucilia sericata than to those of other Lucilia cuprina specimens. This explains the ambiguous identification obtained here for Lucilia sericata. In some cases, the arbitrary similarity threshold, below which matches cannot be used for identification, is too low. Consequently, best close matches with conspecific and allospecific sequences are considered for identification, even if all conspecific best matches are closer to the query than any of the allospecific ones. To solve this problem, the similarity threshold can be adapted according to the gap between intra- and interspecific distances observed in this particular group of species (Lefébure et al. 2006, Collins and Cruickshank 2013, Puillandre et al. 2012, Virgilio et al. 2012). Here, we only used an arbitrary threshold of 99% similarity. A stricter similarity threshold (e.g. 99.5%) wouldresolve ambiguous identifications obtained for Lucilia ampullacea, for Lucilia sericata (but not when early releases of BOLD are used) and for Cynomya mortuorum (Tables 1 and 2).

Similarity values between the query and its best matches can be calculated using several methods. Here, similarities with GenBank records were determined as 1 - p-distances but no explicit information was found on the exact method used by the IDS of BOLD to determine the similarity values. Even if the IDS of BOLD applied a different method than ours, – distances are standardly corrected using the Kimura 2-parameter model (Kimura 1980) in DNA barcoding (Hebert et al. 2003) – the two searches (1 and 2) using the same queries against the same public records resulted in an identical list of highly similar best matches. Several studies have indeed observed that biases due to different distance calculation methods are less severe with similar sequences than with divergent ones (Collins et al. 2012, Fregin et al. 2012).

It is striking that identifications provided by GenBank and BOLD for the barcode fragment were either ambiguous or involved a rather limited number of very similar reference sequences (Figure 1). Therefore, we tested alternative search strategies to optimise the number of best matches and minimise the number of ambiguous identifications. For this, we used different options offered by GenBank and BOLD by 1) including early releases from BOLD in the reference library, 2) adding the keyword “barcode” as a query in GenBank and 3) using longer COI sequences as queries in GenBank.

Including early releases as reference sequences in BOLD increased the number of best matches of > 99% similarity but also increased the proportion of ambiguous identifications (Table 1). Early releases might not have passed all controls that authors and reviewers make in the process of publication (e.g. Schindel et al. 2011). They are therefore more prone to errors. However, their early release allows the detection of errors and inconsistencies before publication, which is an efficient way to improve the quality of the reference libraries. In addition, they largely outnumber the published sequences and may include precious additional information such as rare haplotypes.

In order to improve the search for sequences that have been produced for DNA barcoding purposes, we added the word “barcode” to each query in GenBank. With this procedure, the number of best matches of > 99% similarity and the proportion of ambiguous matches drastically decreased. The same tendency was observed when longer COI sequences (1412–1534 bp) were used as queries. This is due to the smaller number of reference sequences that are tagged as barcodes or are longer than the standard barcode fragment. Therefore, this kind of search is currently only relevant for the identification of fly species of forensic interest that are well represented by longer COI reference sequences or that are tagged as barcodes. Moreover, longer DNA fragments are not always easy to sequence from degraded forensic samples (Mazzanti et al. 2010). Due to the limited number of best matches of > 99% similarity retrieved by these two options, it was not possible to assess their benefit when trying to minimise the proportion of ambiguous identifications.