also might be incorrectly classified, i.e., should be lumped into single species. If most of the problematic cases refer to introgression and incomplete lineage sorting, then nuclear markers need to be developed to reliably discern between the affected species [31]. If recent speciation and generally low genetic distances (but reciprocally monophyletic haplotype lineages) are involved, then character based DNA barcoding may be more appropriate and would allow to sidestep the problem to find appropriate threshold values by searching ‘‘barcoding gaps’’. In any case, where not only species identification but species discovery is concerned, it is clear that DNA barcoding should be used as only one (in many groups the first preliminary) step in the recognition, diagnosis and description of species in terms of integrative taxonomy (e.g. [32]). Materials and Methods Data sampling The study was carried out in compliance with the institutional guidelines on animal husbandry and experiments of the the Zoological Museum of the University of Amsterdam. In addition, the authorization for the experiments was given by the Iranian (permission number: 3–5360) and Moroccan (p.n.: 04666 DCRF/ CPB/PFF) authorities. We sequenced 210 individuals of 145 nominal species for DNA sequences of cox1 and 16S gene fragments, of which 31 and 46 species were parapatric species for cox1 and 16S respectively. Parapatric species are defined as species that have at least one other closely-related species which inhabits a continuous range, the two species excluding each other geographically [15,33]. The range boundary between the two taxa has no dispersal barrier, allowing parapatric species to hybridise and display intergradation in their contact zones, yet they maintain distinct outside of these zones [34,35]. DNA was extracted from tissue or blood samples using DNeasy Tissue Kits (QIAGEN) following the manufacturer’s protocol. Polymerase chain reactions (PCR) and sequencing reactions follow protocols described by Aliabadian et al. (2007) [36] which can be summarized as follows. A fragment of cox1, was sequenced using two primer combinations that amplify a region of 612 bp starting from the 59 terminus of the mitochondrial cox1 gene: BirdF1 (59 - TTC TCC AAC CAC AAA GAC ATT GGC AC -39), BirdR1 (59 -ACG TGG GAG ATA ATT CCA EET CCT G- 39), and BirdR2 (59 -ACT ACA TGT GAG ATG ATT CCG AAT CCA G - 39) [9]. A fragment of 16S was sequenced for the same individuals using 16SA-L (light chain; 59-CGC CTG TTT ATC AAA AAC AT-39) and 16SB-H (heavy chain; 59-CCG GTC TGA ACT CAG ATC ACG T-39) [37]. PCR products were cleaned using QIAquick PCR Purification Kit (Qiagen). Sequencing reactions were resolved on ABI 3100 or ABI 3730 automated DNA sequencers. Genbank accession numbers of newly deter- mined sequences are FJ465179–FJ465383 and are listed in detail in Table S2). Our data set was complemented by cox1 and 16S sequences from GenBank, as available on 1 July 2006). For cox1, additional sequences were included from the Barcode of Life Data Systems website (http://www.barcodinglife.org/, as ac- cessed on 1 July 2006). Sequences were included provided they had a length of .612 (cox1) and .538 bp (16S) homologous to our sequences, with no more than 50 ambiguous or missing nucleotides. Cob sequences with a length of .1000 bp and no more than 50 ambiguous or missing nucleotides were retrieved from GenBank as well. Because of a probably high prevalence of misidentification, erroneous sequences or NUMTs in the DNA Barcoding in Birds Table 3. Number of individuals and taxa employed in this study. Individuals species genera families orders all birds cox1 cob 16S 2776 4614 708 9721 756 2087 498 2161 329 890 270 244 75 114 91 25 20 24 25 doi:10.1371/journal.pone.0004119.t003 Genbank sequences, we submitted these to a rigurous quality control. All sequences per gene were aligned and a Neighbor- joining tree produced. We identified, in this tree, all sequences clustering far from their known taxonomic or phylogenetic position, or characterized by extreme branch lengths, and omitted these sequences from further analysis. For cob, the gene were the largest data set (over 10,000 sequences) was initially downloaded, less than half of these were of sufficient length and quality. All sequences were aligned using Muscle, a multiple alignment software for protein and nucleotide sequences which allows multiple sequence comparison by log-expectation [38]. Probably erroneous sequences (with highly unlikely positions or extreme branch lengths, based on a neighbour-joining tree calculated with all sequences) were identified by eye and omitted. A total of 2776 sequences (756 nominal species) were kept for cox1, 708 (498 species) for 16S, and 4614 (2087 species) for cob (Table 3), and altogether 2719 species were included for at least one of the genes (Table S2). Among conspecific sequences, we verified that for many species, samples from distant localities were included and our analysis is thus not based on including only samples from the same locality and population. Data analysis Genetic distances were calculated to quantify sequence divergences among individuals using Kimura’s (1980) [39] two- parameter (K2P) models, theta, as implemented in MEGA 3.1 [40]. The K2P distance is the most effective model when genetic distances are low [8]. K2P distances were calculated at all taxonomic levels, intraspecific, intrageneric, intrafamilial, intraor- dinal, and, separately, between species in parapatric species pairs, following a published taxonomy [41] and unpublished data by CS Roselaar. For calculation involving higher taxonomic levels, pairwise comparisons of the previous levels were excluded (e.g., in comparisons of intrageneric, pairwise distances of intraspecific were removed, and only pairwise distances of those samples to other species were used). Average K2P distances were computed based on pairwise comparisons of all sequences for each species, and each pair of parapatric species. To calculate intra- and interspecific pairwise distances, based on output matrix of MEGA 3.1, we wrote a converter program SPD 1.1) in C language (this program will be available from the authors upon request). Altogether 3,823,995, 9,952,500, and 250,257 pairwise distanc- es were compared in this study for, cox1, cob, and 16S respectively. NJ trees of K2P distances showing inter- and intra specific variations were constructed using MEGA 3.1 (not shown here). A regression analysis was employed to assess the effect of sample size on intraspecific divergences for each gene using SPSS for Windows, version 11. PLoS ONE | www.plosone.org 7 January 2009 | Volume 4 | Issue 1 | e4119