This study includes material from the mountain regions of Kirgizistan, Tadzhikistan, northern Afghanistan, northern Pakistan and India (e.g. mountain ranges Tian Shan, Hindu Kush, Karakorum, Himalaya) and the mountain regions in the Chinese provinces Qinghai, Gansu, Sichuan, Yunnan and the autonomous regions Tibet and Xinjiang Uygur. The Colias fauna of these Central Asian regions comprises about 34 species (Verhulst 2000) while the species number for Central Asia in broad sense is over 40 species.

The taxon sampling aimed to cover as many of the Colias species from this area as possible. Additionally, a few Colias species occurring in adjacent territories (e.g. Buryatia) were also available for molecular study. Whenever possible, several individuals of each species were analysed to assess intraspecific variation. The available specimens used for molecular study consisted of a total of 56 adult specimens covering 27 species of Central Asian Colias and two Colias species from adjacent territories (Table 1). The specimens are preserved as DNA voucher specimens and labelled accordingly, to be deposited in the collections of the Zoological Museum of Finnish Museum of Natural History, Helsinki, Finland (MZH) (DNA voucher specimens MZH_JL1-JL71). Species identifications were verified by JL based on easily recognizable diagnostic characters using the monograph by Verhulst (2000), while the taxonomy is according to Grieshuber and Lamas (2007). Additionally, we used 35 COI barcode sequences (17 species) of Palaearctic Colias species obtained from GB, as listed in Table 2.

Total genomic DNA was extracted form 2-5 legs of dried, pinned butterfly specimens using NucleoSpin® Tissue Kit (Machery-Nagel), according to manufacturer’s protocols, and resuspended in 50 µl ultrapure water.

The primer pair LCO-1490 (5’-GGTCAACAAATCATAAAGATATTGG-3’) and HCO-2198 (5’-TAAACTTCAGGGTGACCAAAAAATCA-3’) (Folmer et al. 1994) was used to amplify a ca. 650 bp fragment of the mitochondrial COI gene. The polymerase chain reactions (PCR) were done under the following parameters: initial heating 95 °C for 2 min, following 30 cycles of 94 °C for 30 s, 49 °C for 30 s and 72 °C for 2 min, followed by a final extension of 72 °C for 7 min. The primer pair RpS2 nF (5’-ATCWCGYGGTGGYGATAGAG-3’) and RpS2 nR (5’-ATGRGGCTTKCCRATCTTGT-3’) (Wahlberg and Wheat 2008) was used to amplify a ca. 400 bp fragment of the nuclear RpS2 gene. The PCR were carried out following the PCR cycling profile described in Wahlberg and Wheat (2008): initial heating 95 °C for 7 min, 40 cycles of 95 °C for 30 s, 50 °C for 30 s, 72 °C for 2 min, and a final extension period of 72 °C for 10 min. Sequencing of the double-stranded PCR product was carried out on an ABI PRISM® 377 Automated Sequencer (Applied Biosystems) following manufacturer’s recommendations. All PCR primers were used for sequencing. Sequences were inspected and edited using Sequence Navigator® (Applied Biosystems).

We analysed and clustered our sequence data using parsimony and Neighbour-Joining (NJ) of K2P-distances. We used parsimony and NJ for our newly generated COI sequence dataset, NJ for RpS2 sequences, parsimony for the concatenated COI and RpS2 sequences, and, finally, NJ for the combined COI sequences generated in this study and those in GB. All trees were rooted using Papilio glaucus (family Papilionidae) and Aporia crategi (Pieridae, subfamily Pierinae) as outgroup taxa.

Parsimony analysis was performed using NONA (Goloboff 1999) and spawn with the aid of Winclada (Nixon 2002), using a heuristic search algorithm with 1000 random addition replicates (mult*1000), holding 10 trees per round (hold/10), max trees set to 10 000 and applying TBR branch swapping. All base positions were treated as equally weighted characters. Nodal support was assessed with bootstrap resampling (1000 replicates) using Winclada (Nixon 2002). MEGA5 (Tamura et al. 2011) was used for NJ clustering using 1000 bootstrap replicates. The Kimura 2-parameter model was used for NJ clustering of the COI sequences, while the Tamura-Nei model with gamma distributed rates was chosen for the RpS2 sequences.

We obtained a 643 bp COI barcode for 56 Colias specimens, and a 409 bp fragment of RpS2 was obtained for 49 specimens (Table 1). A+T content of the COI sequences was 69.22%, and of the RpS2 45.0%. There were 115 parsimony informative sites for COI and 39 for RpS2.

Uncorrected pairwise divergences between ingroup taxa ranged between 1.09 and 4.09% (mean 2.77%) for COI and 0.0–1.7% (mean 1.0%) for RpS2. GenBank accession numbers are given in Table 1. Intraspecific uncorrected distances were up to 1.09% (in Colias thisoa) for COI, with specimens of most species differing by less than 4 nucleotide changes.

The parsimony analysis of the new COI sequences yielded four equally parsimonious trees (CI = 0.59, RI = 0.75) the strict consensus tree of which is presented in Figure 1. The NJ tree is presented in Figure 2.

The majority of the species could be identified with COI alone, as no COI haplotypes were shared between species. Both parsimony and NJ trees recovered 25 (out of 28) species as monophyletic groups (Figures 1–2). Neither Colias cocandica, nor Colias nebulosa formed monophyletic entities, as their sequences were scattered over various parts of the trees. The two samples of Colias tyche were not recovered as sister taxa, for sample MZH_JL5 appeared as sister taxon of Colias heos. The overall topologies of the parsimony and NJ trees were identical, except for the placement of Colias thrasibulus. Parsimony placed the taxon as sister to a clade of five taxa (Figure 1), while NJ placed it as sister to Colias romanovi (Figure 2). The external morphology of Colias thrasibulus is rather different from that of Colias romanovi, while some similarities can be found between Colias thrasibulus and Colias nina, Colias ladakensis, Colias tibetana and Colias cocandica (Figure 1). Only 17 of the 39 parsimony informative sites of RpS2 were variable among the 49 ingroup members. NJ only recovered few species as separate lineages due to the shallow divergences (Figure 3). The information content of this gene region is best interpreted as a character-based diagnostic table, as suggested by DeSalle et al. (2005). This gene region yielded species specific (diagnostic) haplotypes for 11 species out of 33 (Table 3).

The parsimony analysis of COI + RpS2 yielded nine trees of length 560 steps (CI = 0.63, RI = 0.72), the strict consensus tree of which is shown in Figure 4. Colias cocandica, Colias nebulosa and Colias tyche were not monophyletic and Colias thrasibulus had the same position as in the COI cladogram (Figure 1).

The strict consensus cladogram for all the available COI data resolved the taxa in the same positions as in the tree of the new COI sequences only. For ten species of the present study sequences were also available from GB. Sequences of most species clustered together as monophyletic entities, except for Colias nebulosa, Colias cocandica, Colias tyche and Colias regia. For Colias regia the GB sequence (GB accession no FJ663427) did not cluster together with our sequences. The GB barcodes of Colias erate and Colias croceus were shared by these two taxa.

Neither the Himalayan and south Tibetan adjacent mountain Colias fauna (Colias berylla, Colias ladakensis, Colias nina, Colias stoliczkana, Colias thrasibulus, Colias tibetana), nor the east Tibetan, Qinghai, Gansu and Sichuan species aggregates (Colias adelaidae, Colias grumi, Colias lada, Colias montium, Colias nebulosa, Colias sifanica, Colias wanda) were resolved as species clusters similar to the Tian Shan, Pamir and Hindukush species.

Several COI haplotypes were noted for a few species, even among specimens obtained from the same locality (e.g. Colias staudingeri and Colias thisoa). Taxa not resolved as monophyletic clusters were the species Colias cocandica and Colias nebulosa. All the included subspecies of Colias cocandica (Colias cocandica cocandica, Colias cocandica pljutshtshi and Colias cocandica hinducucia) showed distinct COI sequences, with Colias cocandica cocandica as most different.

Lukhtanov et al. (2009) tested the utility of COI barcodes for Central Asian butterflies by sampling specimens from a considerable geographical range. They observed that this substantially increased intraspecific variation reducing the interspecific divergences (“barcoding gap”), but that this did not hamper species identification. The present study shows that most Colias taxa form monophyletic entities that can be identified with COI data alone. The RpS2 gene region showed identical sequences in Colias cocandica pljutshtshi and Colias cocandica hinducucia (Table 3, Figure 3), differing by only three nucleotides from Colias cocandica cocandica. Based on the molecular data the recognition of these subspecies is not or weakly supported.

The fact that the three Colias nebulosa samples were scattered over different parts of the COI tree might be the result of a laboratory contamination due to carry over between samples. The Colias nebulosa samples were collected on the same day and in the same place. Colias nebulosa is morphologically distinct from other Colias species, excluding possible misidentification. The RpS2 data, however, could point to two morphologically cryptic species in sympatry (samples MZH_JL24 vs. MZH_JL9 and MZH_JL26), so that the different COI barcodes might represent numts, despite no apparent ‘signs’ (no indels). This discrepancy between morphology and DNA sequence data emphasises the necessity to use multiple samples to detect this sort of challenging issues.

Even though Colias cocandica and Colias nebulosa did not form monophyletic groups our results show that COI barcodes are useful for (1) identifying Palaearctic and Central Asian Colias, (2) pointing to a possible cryptic species, and (3) highlighting the necessity to further investigate the question on the subspecific rank of Colias cocandica cocandica.

The utility of RpS2 as a species barcode for Colias spp. is clearly more limited, since e.g. Colias heos, Colias lada, Colias nina, Colias thisoa of the subgenus Eriocolias and Colias tyche (subgenus Eucolias) have identical sequences (Table 3, Figure 3). Still, RpS2 yielded species specific (diagnostic) haplotypes for 11 species of the subgenus Eriocolias and for Colias hyale (subgenus Colias s.str.).

The strict consensus tree was more resolved than either of the trees resulting from separate analyses of the gene regions (Figure 4).

Although the concatenated data did not resolve the phylogenetic relationships among all Colias species, some observations can be made. The majority of the species confined to the adjacent Tian Shan, Pamir and Hindukush mountain ranges form a well supported clade. This includes Colias eogene, Colias regia, Colias romanovi, Colias marcopolo, Colias staudingeri, Colias christophi, Colias alpherakii and Colias wiskotti. Yet, Colias sieversi, which also occurs in these mountain ranges (Peter I and Khozratishoh mountains), was not included in this clade. Colias sieversi is morphologically most similar to Colias alpherakii, thus showing another case of disagreement between morphological and DNA sequence data. Colias thisoa, too, lives in the aforementioned mountain ranges, but it has a wider distribution, stretching from Turkey to the Altai Mountains. A third taxon, Colias cocandica cocandica, is considered closely related to Colias tamerlana (e.g. Verhulst 2000), a species occurring in southern Siberia and Mongolia. Thus, the origin of Colias thisoa and Colias cocandica cocandica may differ from that of the species confined to the Tian Shan, Pamir and Hindukush mountain range. One sample of Colias cocandica (MZH_JL43) was placed within this “mountain” clade, while the other two samples appeared as sister taxa to the Himalayan species Colias ladakensis. As with Colias sieversi, our DNA data disagree with the morphological characters, but it should be noted that this clade is not well supported. Conversely, two morphologically similar Himalayan species, viz. Colias nina and Colias ladakensis, were assigned to different clades. In the COI + RpS2 tree they were placed in different, more encompassing species clusters (Figure 4), in the COI NJ tree they were joined with Colias cocandica pljutshtshi and Colias cocandica hinducucia (Figure 2), while the COI cladogram resolved these taxa together with Colias adelaidae, Colias tibetana, Colias cocandica pljutshtshi and Colias cocandica hinducucia (Figure 1).

The analyses did not support the monophyly of the subgenera Eucolias and Eriocolias sensu Berger (1986). The Eucolias species Colias tyche was not resolved as a separate monophyletic lineage, but was resolved into Eriocolias. This is congruent with the results of Pollock et al. (1998) and Brunton (1998). Only the the subgenus Colias, here represented by Colias hyale, is supported as a distinct lineage, placed as sister to all other Colias sp.

The parsimony (Figure 5) and NJ analyses (Figure 6) of the larger matrix of Palaearctic COI barcodes (total COI) recovered the same species clusters, but some of the species show different placements (e.g. Colias thisoa, Colias christophi). This is not surprising as all internal nodes are very shallow. The samples of Colias tyche and Colias hyperborea show very low sequence difference, morphologically these taxa are different, and they largely share the same distribution area. An example of species that share the same distribution and that exhibit clear morphological similarities, and which as such were resolved as sister species in both analyses, includes Colias wiskotti and Colias alpherakii. Identification of Palaearctic Colias based on COI barcodes is in most cases possible, since shared haplotypes were recorded only for Colias erate and Colias croceus.

Intraspecific variation is notable between some of the recognized subspecies, both among our own samples and those downloaded from GB. The intraspecific variation can partly be explained by morphologically clearly distinct subspecies, such as those of Colias wiskotti, or by specimens from widely different localities, such the different specimens of Colias hyale (sample FJ663418 from Russia, FJ663421 from Kazakhstan, HQ004297 from Romania and MZH_JL35 and MZH_JL44 from SW Transbaikalia). However, notable intraspecific variation also occurs within populations, such as Colias thisoa aeolides with all samples originating from the same locality and date, but the limited sampling prevents conclusions on the reasons for this. It is apparent that the understanding of intraspecific variability of the COI barcode for Colias is presently very limited.

The combined COI data of our sequences and sequences downloaded from GB include species belonging to one additional subgenus, Neocolias, represented by Colias myrmidone and Colias erate. Only the subgenus Colias, represented by Colias hyale, is well supported as distinct lineage. Yet, one specimen of Colias hyale (FJ663419) clustered together with Colias erate (Neocolias) and Colias croceus (Eriocolias). The other subgenera were not resolved as clades according to present classification, in agreement with our results for the combined analysis.

Our findings generally support COI as a species specific barcode for Colias, but we also highlight the necessity of including multiple individuals of species in molecular barcoding studies. Problematic ‘cases’ of widely divergent barcodes or conflicting morphological and molecular ‘signals’ are found in most if not all barcoding studies, and this study makes no exception.