DNA barcode analysis of butterfly species from Pakistan points towards regional endemism

DNA barcodes were obtained for 81 butterfly species belonging to 52 genera from sites in north-central Pakistan to test the utility of barcoding for their identification and to gain a better understanding of regional barcode variation. These species represent 25% of the butterfly fauna of Pakistan and belong to five families, although the Nymphalidae were dominant, comprising 38% of the total specimens. Barcode analysis showed that maximum conspecific divergence was 1.6%, while there was 1.7–14.3% divergence from the nearest neighbour species. Barcode records for 55 species showed <2% sequence divergence to records in the Barcode of Life Data Systems (BOLD), but only 26 of these cases involved specimens from neighbouring India and Central Asia. Analysis revealed that most species showed little incremental sequence variation when specimens from other regions were considered, but a threefold increase was noted in a few cases. There was a clear gap between maximum intraspecific and minimum nearest neighbour distance for all 81 species. Neighbour-joining cluster analysis showed that members of each species formed a monophyletic cluster with strong bootstrap support. The barcode results revealed two provisional species that could not be clearly linked to known taxa, while 24 other species gained their first coverage. Future work should extend the barcode reference library to include all butterfly species from Pakistan as well as neighbouring countries to gain a better understanding of regional variation in barcode sequences in this topographically and climatically complex region.

Introduction DNA barcoding has emerged as a useful tool for the identification and discovery of animal species. It employs sequence diversity in a 648 base pair fragment near the 5′ end of the mitochondrial cytochrome c oxidase subunit I (COI) gene as a tool for species discrimination (Hebert et al. 2003a). Barcoding has been shown to discriminate species across the animal kingdom (Tyagi et al. 2010;Virgilio et al. 2010) including fishes, mammals, birds, insects, crustaceans and many other groups (Hebert et al. 2004a;Foottit et al. 2008;Hastings et al. 2008;Hubert et al. 2008;Hou et al. 2009;Wong et al. 2009;Clare et al. 2011). Reflecting the rapid growth in barcode coverage (Jinbo et al. 2011), BOLD, the Barcode of Life Data System (Ratnasingham & Hebert 2007), now includes records for more than 261K animal species. The order Lepidoptera has received particular attention (Hajibabaei et al. 2006;Silva-Brandao et al. 2009;Hebert et al. 2010;Kim et al. 2010) with 691K barcode records on BOLD (Feb 3, 2013), including data for 9124 named butterfly (Papilionoidea, Hesperioidea) species from 194 countries.
The effectiveness of DNA barcoding has spurred efforts to construct DNA barcode reference libraries for various animal groups (Ekrem et al. 2007;Guralnick & Hill 2009;Janzen et al. 2009;Lee et al. 2011;Zhou et al. 2011;Webb et al. 2012). These libraries not only aid the documentation of biodiversity (Janzen et al. 2005;Naro-Maciel et al. 2010) including endangered species (Elmeer et al. 2012;Vanhaecke et al. 2012), but can disclose endemism (Bossuyt et al. 2004;Quilang et al. 2011;Sourakov & Zakharov 2011). Because Lepidoptera have been selected as a model group for intensive analysis, the order is well represented on BOLD, but some regions such as South-East Asia have seen little investigation. Barcode records are available for a significant fraction of the Central Asian butterfly fauna (Lukhtanov et al. 2009) and for a smaller number of species from Western India (Gaikwad et al. 2012). However, these studies fail to provide coverage for many species known from Pakistan (Roberts 2001). The current study had the primary goals of testing the effectiveness of DNA barcodes in the identification of butterfly species from Pakistan and comparing these records with those from other regions to gain a better sense of the extent of intraspecific variation.

Specimen sampling
Butterflies were collected at 107 locations across central and northern Pakistan ( Fig. 1) during 2009-2012. These sites included three different climatic zones: tropical, subtropical and temperate, with altitudes ranging from 127 to 2660 m, and both agricultural and forested environments. Each specimen was labelled, assigned a code number and deposited in the arthropod collection at the National Institute for Biotechnology and Genetic Engineering (NIBGE), Faisalabad, for subsequent morphological and molecular analysis. Using standard guides to the fauna (Malik 1973;Hasan 1994;Roberts 2001), the 407 specimens were assigned to 81 species belonging to 52 genera. Two species (Lasiommata sp. MA01 and Polycaena sp. MA01) could only be identified to a generic level, but were included in the analysis. Specimen data and images are available on BOLD (Ratnasingham & Hebert 2007) in the project MABUT (Barcoding Butterflies of Pakistan). Fifty-nine of the 81 species were represented by more than one specimen (range 2-20). All sequences generated in this study are available on BOLD (Process IDs: MABUT001-10 to MABUT312-12; MABUT326-13 to MABUT388; MAIMB133-09 to  and on GenBank under the following accession nos: KC158311-KC158471, HQ990321-HQ990449, HQ990705, HQ990728-HQ990729, GU681850-GU681851, GU681855-GU681856, GU681859, GU681870 and GU681872-GU681875.

DNA extractions and PCR amplifications
A single leg was removed from each specimen with a sterile forceps and transferred to a 96-well microplate preloaded with 30 lL of 95% ethanol in each well. DNA extraction, PCR amplification and sequencing were performed at the Canadian Centre for DNA Barcoding (CCDB) following standard protocols (Ivanova et al. 2006(Ivanova et al. , 2007Ivanova & Grainger 2007a,b,c). DNA extractions were performed by following the protocols developed for invertebrate barcoding (Ivanova et al. 2006). Amplification of the COI-5′ barcode region was performed with primer pair LepF1 (ATTCAACCAATCA TAAAGATATTGG)/LepR1 (TAAACTTCTGGATGTCC AAAAAATCA) (Hebert et al. 2004b) using the following PCR conditions: 94°C (1 min); 5 cycles of 94°C (30 s), 45°C (40 s), 72°C (1 min); 35 cycles of 94°C (30 s), 51°C (40 s), 72°C (1 min); and final extension of 72°C Fig. 1 Map of Pakistan and neighbouring nations showing collection localities for this study as well those for specimens examined in a prior study (Lukhtanov et al. 2009).
(10 min). PCRs were carried out in 12.5 lL reactions containing standard PCR ingredients and 2 lL of DNA template. PCR products were analysed on 2% agarose Egel â 96 system (Invitrogen Inc.). Amplicons were sequenced bidirectionally using BigDye Terminator Cycle Sequencing Kit (v3.1) on an ABI 3730XL DNA Analyzer. The forward and the reverse sequences were assembled and aligned using CodonCode Aligner (CodonCode Corporation, USA). Sequences were also inspected and translated in MEGA V5 (Tamura et al. 2011) to verify that they were free of stop codons and gaps.

Data analysis
The sequence from each specimen was compared with barcode sequences on GenBank using 'Blast' and with sequences on BOLD using the 'Identification Request' function. Prior studies have revealed that most different species of Lepidoptera show >2% sequence divergence at CO1 (Hebert et al. 2003b), and researchers have used a 2% pairwise distance threshold for species delimitation (Strutzenberger et al. 2011). For the barcode-based identity analysis, we also used a threshold of 2% divergence. DNA barcodes for 9124 butterfly species from 194 countries are currently available on BOLD, all readily available for sequence comparisons. In addition, the results were compared with those of prior studies in Central Asia (353 butterfly species) (Lukhtanov et al. 2009), Korea (83 species) (Kim et al. 2010) and India (40 species) (Gaikwad et al. 2012). ClustalW nucleotide sequence alignments (Thompson et al. 1994) and NJ clustering analysis were performed using MEGA V5 (Tamura et al. 2011). The Kimura-2-Parameter (K2P) (Kimura 1980) distance model was used, along with pairwise deletion of missing sites, with nodal support estimated using 500 bootstrap replicates. The online version of Automatic Barcode Gap Discovery (ABGD) (Puillandre et al. 2012) was used for both pairwise distance analyses and to generate distance histograms and distance ranks. The presence or absence of a 'barcode gap' (Meyer & Paulay 2005) was also determined for each species as a test of the reliability of its discrimination. Using the barcode gap criterion, a species is distinct from its nearest neighbour (NN) if its maximum intraspecific distance is less than the distance to its NN sequence. The 'Barcode Gap Analysis' (BGA) was performed using BOLD. Species identification success by 'Best Match' and cluster analysis was performed using TaxonDNA (Meier et al. 2006). The relationship between geographical distance and intraspecific genetic distance was analysed separately for each species (with at least three individuals and three locations) using the Mantel test (Mantel 1967) and by linear regression using XLSTAT (version 2013.3.02; Addinsoft, Inc., NY, USA).

Results
Barcode sequences greater than 500 base pairs (bp) were recovered from 374 of the 407 specimens (92%), providing at least one sequence for each of the 81 butterfly species. When these sequences were compared with those in the BOLD and NCBI databases, close sequence matches (<2% divergence) were detected for 55 of the species from Pakistan, while 26 lacked a match. The highest number of matches involved records from India (15), Central Asia (11) and Korea (10). Figure 2 presents results from the ABGD and BGA analyses. Distance values show a gap between the intraspecific and the interspecific distances ( Fig. 2A). As well, both the maximum and mean distances to NN are higher than the respective intraspecific distances for all species (Fig. 2B). Nearest neighbour distances were more than 3% for all but three species pairs: Tarucus balkanicus vs. T. rosaceus (1.70%), Junonia orithya vs. J. hierta (2.49%) and Celastrina huegelii vs. C. argiolus (2.64%). Intraspecific distances could not be determined for the 22 species with just a single representative, but NN distances were greater than 4% for 21 of them.
NJ clustering analysis showed that each of the 81 species formed a monophyletic cluster (Fig. 3). Species with two or more barcode sequences were analysed for species identification using TaxonDNA. When a 3% threshold was employed, 100% of the species were correctly identified using the 'Best Match or Best Close Match' criterion. Analysis of the 374 sequence records using TaxonDNA led to the recognition of 78 clusters at a 3% threshold and 80 clusters at a 2% threshold. At the 3% threshold, 75 of the 78 clusters were comprised of a single species, with the largest pairwise intraspecific distance being 2.88%, while 79 of the 80 clusters were a single species at the 2% threshold with the largest pairwise intraspecific distance being 1.67%.
Genetic divergences increased with taxonomic rank (Table 1; Fig 2) with little overlap between conspecific and congeneric distances. Intraspecific divergences ranged from 0.0 to 1.6% with a mean of 0.2%, while divergences for the species in a genus ranged from 1.7 to 14.3% with a mean of 8.0%. The distances within families ranged from 3.9 to 19.2% with a mean of 13.1%. Fifty-five species were represented by at least one conspecific from another country, but in most cases, there was little increase in intraspecific divergence when they were included in the analysis (Table 2). Seventeen species showed a three-fold or more increase in intraspecific distances (Table 2, bold-faced numbers), but their maximum intraspecific divergence remained <3%, and mean divergence was <1% in all cases except Colotis amata (max = 3.20%, mean = 1.17%) ( Table 2).
The relationship between geographical and genetic distances was quantified by plotting geographical distances against intraspecific variation (K2P). Table 2 provides species-wise Mantel correlation statistics, while Fig. 4 shows the overall trend between geographical distance and intraspecific genetic divergence. Some species showed a strong correlation between the two parameters, as genetic distances increased with geographical distance, but others did not show a significant relationship between the two variables (Table 2). Overall, this analysis showed a weak relationship (R 2 = 0.22; y = 8E-05x + 0.250) between the geographical extent of a species and its maximum intraspecific divergence (Fig. 4).

Identification success for the butterflies of Pakistan
This study has begun the construction of a DNA barcode reference library for the butterflies of Pakistan. Cluster analysis revealed that all 81 species examined in the study formed a monophyletic cluster which corresponded perfectly with the taxa recognized on morphological criteria. Although three species pairs showed limited divergence (<3%), maximum intraspecific divergence was always lower than the NN distance, enabling the separation of all species. Even the most closely related (1.70%) species pair, Tarucus balkanicus and T. rosaceus, was separated with strong bootstrap support in the NJ tree. Our results confirm the usefulness of DNA barcoding in identifying the butterflies of Pakistan, but the sample size was low for some species and 75% of the fauna awaits analysis.
When sequences for butterfly species from Central Asia (Lukhtanov et al. 2009) were included, eight species pairs formed paraphyletic clusters. Among these pairs, the NN distance between Aglais caschmirensis (from Pakistan) and A. nixa (from Uzbekistan) was 0.2%, while that between A. caschmirensis and A. urticae (from Kazakhstan) was 1.4%. Although NN distances for these sister species pairs were small, barcode-based identifications  were possible as reported by Tavares & Baker (2008) in their study on sister species of birds.
'Barcode Gap Analysis' showed that NN distance for all the species was higher than the maximum intraspecific distance. The Barcode Index Number (BIN) system (Ratnasingham & Hebert 2013) provided further evidence of the genetic distinctiveness of the species as it assigned the 81 species to 80 BINs with only T. balkanicus and T. rosaceus sharing a BIN. When identity analysis was performed using Best Match/Best Close Match at a 3% threshold, all the species were correctly identified. Other studies have generally reported similar results (Janzen et al. 2005;Lukhtanov et al. 2009;Gaikwad et al. 2012) with a few exceptions. For example, Gaikwad et al. (2012) found that intraspecific divergence was higher (7.8%) in the butterfly Lethe europa than the distance to its NN (7.4%). Such cases can, of course, arise through a failure to discriminate sibling taxa. Bortolus (2008) has emphasized the importance of detailed taxonomic study in cases where DNA barcode results are discordant with taxonomic assignments. Costa et al. (2012) have reinforced this conclusion, noting the need for a ranking system to register the certainty of identifications for specimens used to develop reference barcode libraries. These suggestions reinforce the importance of an integrative approach to species delimitation by considering morphological, genetic, ecological and geographical information, rather than considering taxonomic identifications as facts against which to 'test' DNA barcoding (e.g. Smith et al. 2008). Nevertheless, focusing on one region of the genome is useful to the community for generating a comparable set of sequences across a large number of diverse taxa and geographical regions.
Genetic divergence patterns with increasing geographical distance: a regional Asian perspective The within-species divergence values for most species in the study were under the 2%. In most cases, the addition of conspecific sequences from other countries increased the intraspecific distance, but the relationship between geographical distance and the level of intraspecific divergence was not strong. In a few cases, substantial intraspecific distances were observed between specimens from the same region. For example, Pelopidas mathias collected from sites in Pakistan <250 km apart showed 1.54% divergence. On the other hand, Deudorix epijarbas from Pakistan and Taiwan (4832 km) lacked barcode divergence. Other species showed regional variation that was not linked to distance. For example, specimens of Lampides boeticus from Pakistan and Queensland Australia were just 0.4% divergent, but specimens from Papua New Guinea were 1.9% divergent. These results reinforce previous conclusions that geographical distance is often associated with an increased genetic divergence, but that the increase is too small to impede the identification of species (Lukhtanov et al. 2009;Bergsten et al. 2012;Gaikwad et al. 2012).
Diversity hotspots and endemism in Asia underscores the need for regional barcode libraries Although Pakistan and neighbouring Central Asia are only 700 km apart, prior studies have indicated that there is little overlap in their butterfly faunas. In fact, just 42 species (14%) are shared among the 320 butterfly species from Pakistan (Roberts 2001) and the 353 species from Central Asia (Lukhtanov et al. 2009). Their distinctive faunas undoubtedly reflect the effectiveness of the Pamir mountain chain, which rises to more than 5000 m, as a dispersal barrier. This limited overlap suggests the presence of multiple regions of endemism in this segment of Asia, mirroring a pattern of low overlap between the biodiversity hotspots in the Western Ghats (India) and Sri Lanka (Bossuyt et al. 2004). Although India and Sri Lanka are on the same continental shelf, and the strait separating them does not exceed 70 m in depth, limited biotic interchanges have left the two areas with an unexpectedly large number of endemics. This fact highlights     The number of individuals of a species included in the analysis is indicated in brackets. A double dash indicates that a given species was presented by only one specimen, and thus, maximum intraspecific divergence is not presented, while bold highlighting is used to indicate those species that exhibit a three-fold or greater increase in intraspecific variation when records outside of Pakistan were included. *Insufficient data to run the Mantel test.
the need to expand barcode coverage for all animal groups from the various subregions in southern Asia. Certainly, barcode reference libraries based on species from other nations will only permit the identification of a fraction of Pakistan's biodiversity.

Supporting Information
Additional Supporting Information may be found in the online version of this article: