Molecular systematics and evolution of the recently discovered “Parnassian” butterfly (Parnassius davydovi Churkin, 2006) and its allied species (Lepidoptera, Papilionidae)

Molecular systematics and evolution of the recently discovered “Parnassian” butterfly (Parnassius davydovi Churkin, 2006) and its allied species (Lepidoptera, Papilionidae)

Gene 441 (2009) 80–88 Contents lists available at ScienceDirect Gene j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / g ...

1MB Sizes 2 Downloads 110 Views

Gene 441 (2009) 80–88

Contents lists available at ScienceDirect

Gene j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / g e n e

Molecular systematics and evolution of the recently discovered “Parnassian” butterfly (Parnassius davydovi Churkin, 2006) and its allied species (Lepidoptera, Papilionidae) Keiichi Omoto a,⁎, Takahiro Yonezawa a,b, Tsutomu Shinkawa c a b c

Hayama Center for Advanced Studies, The Graduate University for Advanced Studies, Hayama, Japan School of Life Sciences, Fudan University, Shanghai, China The University of the Air, Chiba, Japan

a r t i c l e

i n f o

Article history: Received 23 June 2008 Received in revised form 27 October 2008 Accepted 28 October 2008 Available online 14 November 2008 Keywords: Swallowtail butterfly Mitochondrial DNA Apollos Molecular phylogenetic tree Divergence time Sphragis

a b s t r a c t The nucleotide sequence of 807 bp of the mtDNA-ND5 locus of Parnassius davydovi (Churkin, S. 2006. A new species of Parnassius Latreille, 1804, from Kyrgyzstan (Lepidoptera, Papilionidae). Helios (Moskow) 7,142158), was determined. This butterfly was unexpectedly discovered recently in Kyrgyzstan, and we wished to shed light on its molecular phylogenetic relationship to other Parnassian butterflies, as well as to the related taxa in the subfamily Parnassiinae of the family Papilionidae. Using the ML method with the GTR+I+Γ model, we inferred the phylogenetic tree for 60 Parnassius individuals together with materials of the related genera in the subfamily Parnassiinae (Hypermnestra, Archon, Luehdorfia, Bhutanitis, Allancastria, Zerynthia and Sericinus) with Papilio machaon as an out-group. It was found that P. davydovi is a distinct species most closely related to P. loxias in clade VI among the eight clades, or species groups of Parnassius. The morphological diversity in the form of sphragis, the attachment to the female abdomen formed by the male during copulation, is characteristic to this clade, and we inferred the order of emergence of the different sphragis forms during evolution. Attempts to estimate the divergence times between related taxa were also made. It was inferred that the relatively rapid radiation of Parnassian butterflies started at about 24 MYA BP, while P. davydovi diverged from P. loxias at about 10 MYA BP. © 2008 Elsevier B.V. All rights reserved.

1. Introduction In July 2005, a new species of the genus Parnassius (Lepidoptera, Papilionidae) was discovered unexpectedly in the Mold-Too Mts. in Kyrgyzstan, Central Asia, and the news spread rapidly with excitement among the butterfly collectors of the world (Churkin, 2006). The butterflies of this genus, often called “Apollos,” are so popular among collectors and taxonomists that the term Parnassiology has been in use since the 1920s. About fifty species are known for this butterfly group, occurring mostly in high altitude areas of Central Asia, the Himalayas, and Western China. Because of enthusiastic searching by “Apollo” hunters for over 100 years in remote mountain areas of Eurasia, it was considered that a “new” species would not be found any more. The first molecular systematic study of this genus was carried out by Professor Omoto and his colleagues who attempted to shed light on the evolutionary history of the subfamily Parnassiinae of the family Papilionidae (Omoto et al., 2004). They examined mitochondrial DNA of essentially all the fifty species of the genus Parnassius and most

Abbreviations: NJ, neighbor joining; ML, maximum likelihood; BP, bootstrap probability; MYA BP, million years before present. ⁎ Corresponding author. Center for Advanced Studies, The Graduate University for Advanced Studies, Hayama, 240-0193, Japan. Tel./fax: +81 46 858 1578. E-mail addresses: [email protected], [email protected] (K. Omoto). 0378-1119/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.gene.2008.10.030

species of the related genera (Hypermnestra and Archon) of the tribe Parnassiini of the subfamily Parnassiinae, together with the tribe Zerynthiini (Zerynthia, Allancastria, Sericinus and Luehdorfia). In the study, using the sequence data of about 800 bp of the ND5 locus, the NJ and MP molecular phylogenetic trees were inferred. It was shown that the genus Parnassius was a monophyletic group comprised of eight clades, or species groups, which were mostly found to correspond to the subgenera hitherto proposed by a number of “Parnassiologists”, with some exceptions. Judging from the branching pattern in the molecular phylogenetic tree, the genus Parnassius was assumed to have experienced a relatively rapid radiation at a certain geological time, probably during the late Tertiary period (20–30 MYA BP). The sister genus to Parnassius was found to be Hypermnestra, rather than Archon, the butterfly which superficially looks exceedingly similar to those of Parnassius. Contrary to the previously supported classification based on morphological characteristics (Hancock, 1983; Häuser, 1993), the genus Archon was found to belong to the tribe Zerynthiini, with very different looking butterflies such as of the genus Luehdorfia, rather than Parnassiini, at the DNA level (Omoto et al., 2004). Further study adding data of the 16S and ND1 sequences, as well as more detailed statistical analyses, confirmed the result mentioned above (Katoh et al., 2005). Nazari and his colleagues reported detailed molecular systematic studies using a large set of data of five mtDNA

K. Omoto et al. / Gene 441 (2009) 80–88

loci (COI, COII, ND5, ND1, 16S) together with two nuclear DNA loci (EF1α and wg). They confirmed the sister-group relationships between Parnassius and Hypermnestra, as well as between Archon and Luehdorfia (Nazari et al., 2007). More recently, Michel et al. (2008) reported molecular phylogenetic analyses of Parnassiinae using the sequence data of four mtDNA loci (16S, ND1, ND5, COI). Their results essentially confirm the previously mentioned results with more detailed discussion about the relationships among species of the genus Parnassius. The present study aims at determining the molecular systematic position of Parnassius davydovi, and by examining the molecular versus morphological diversities, shedding more light on the evolutionary history of some butterflies of the genus Parnassius. Moreover, based on the data given by Nazari et al. (2007), but using a different statistical method, we report some new estimates for divergence times among the taxa in the Parnassiinae.

81

RAxML ver. 7.0.4 (Stamatakis et al., 2008) with the GTR+I+Γ model (Rodriguez et al., 1990, Yang, 1996). The Bayesian analysis was carried out using Mr. Bayes ver. 3.1.2 (Ronquist and Huelsenbeck, 2003) with the GTR+I+Γ model under two sets of four simultaneous chains of 1,000,000 generations (Burnin was set to be 100,000 generations. Sample frequency was set to be per 100 generations.). The convergence of the log-likelihood was confirmed by checking that the PSRF (the potential scale reduction factor) values of all parameters were close to 1.0. NJ analysis was carried out by MEGA ver. 4.0 (Tamura et al., 2007) with the MCL (Maximum Composite Likelihood)+Γ model (Tamura et al., 2004, Yang, 1996). The shape parameter (α) of gamma distribution, which was estimated by RAxML, was applied to this analysis. The bootstrap analyses were carried out with 100 replications for the ML method and 1,000 replications for the NJ method to assess the confidence of internal nodes. 2.3. Estimation of divergence time

2. Materials and methods 2.1. Material, DNA extraction and sequencing One male Parnassius davydovi specimen, collected in July, 2006, by Sergei Churkin at altitude 2,500–2,600 m in the Moldo-Too Mts., Tien Shan Mountain Range, Kyrgyzstan, was obtained in a dried and mounted condition (Fig. 1). Two legs were taken from the specimen and sent to the laboratory of one of us (T. S.) for mtDNA sequencing. Extraction of the whole genomic DNA from the leg samples was carried out by means of the method commonly used for Parnassius and other Papilionid groups (Yagi et al., 1999, 2001, 2006). Part of the ND5 locus of mitochondrial DNA was amplified by PCR and the nucleotide sequence of 807 bp was determined using a 377-18 DNA sequencer with a Big-Dye Terminator Cycle Sequence Kit (Applied Biosystems). The detailed procedures of PCR and sequencing followed that of Omoto et al. (2004). 2.2. Phylogenetic analyses The nucleotide sequence data used in this study are summarized in Table 1. These data were aligned automatically by Clustal W (Thompson et al., 1994) and carefully checked by eye. All ambiguous sites were excluded from the analyses. The phylogenetic trees were inferred with the ML method (Felsenstein, 1981), the Bayesian method (Yang and Rannala, 1997), and the NJ method (Saitou and Nei, 1987). In all the trees, Papilio machaon hippocrates (Papilioninae) was used as an out-group. The ML analysis was carried out using

We estimated divergence times following two different methods. For the higher-taxa relationships (genus or tribe level), the divergence times were estimated using the method of Thorne et al. (1998), which does not assume the constant rate of the molecular clock. For this analysis, we used TTT (Thornian Time Traveler; http:// abacus.gene.ucl.ac.uk/). This software requires two stages of analyses. The first stage is the estimation of the branch lengths with the ML method and the second stage is the estimation of the divergence time with the Bayesian method. The nucleotide sequence data of ND1, ND5, COI, COII, 16S, LeutRNA, and EF1α (accession numbers are summarized in supplemental Table S1) were used in this analysis. The gaps and missing sites were completely excluded from the analyses. To take account of the different tempo and mode among different genes and among different nucleotide sites of codon, we separated gene data into several partitions. The branch lengths of each partition and their variance-covariance matrix were independently estimated using the ESTBNEW program of TTT with the HKY+Γ model (Hasegawa et al., 1985, Yang, 1996). The parameters of this model were estimated using the BASEML program of PAML ver. 4.0 (Yang, 2007) in advance. Then the divergence times were estimated using the MULTIDIVTIME program of TTT under the chains of 1,000,000 generations (Burnin was set to be 100,000 generations. Sample frequency was set to be per 100 generations.) The constraints of the divergence times were given as follows. The split of Papilionini (Papilio) and Troidini (Troides) was between 82.5 and 89.1 MYA BP (Gaunt and Miles, 2002). The initial split of genus Papilio (P. machaon and P. thoas) was between 35 and 65 MYA BP (Zakharov et al., 2004). The split of Allancastria cretica and A. cerisyi was between 3 and 11 MYA BP (Nazari et al., 2007). For the lower-taxa relationship (species level), we assumed a molecular clock. In this analysis, we used the BASEML program of PAML ver. 4.0 with GTR+Γ model. The parameters (and rate matrix) were estimated for each nucleotide codon site. 3. Results and discussions 3.1. The molecular systematic position of P. davydovi

Fig. 1. Parnassius davydovi Churkin, 2006. Male, collected in the Moldo-Too Mts, Tien Shan Mountain Range, Kyrgyzstan in July, 2006. Scale bar indicates 10 mm.

In the ML tree in Fig. 2, it is shown that P. davydovi is closest to P. loxias belonging to clade VI as mentioned below. The NJ and Bayesian trees showed essentially the same branching pattern (data not shown). Concerning the Bayesian tree, the PSRF value of the all parameters ranged from 1.00 to 1.03. In the original description of P. davydovi, S. Churkin already mentioned the morphological similarity between these two species, and comparing the wing patterns and male genitalia he concluded that the two species are “good” species, rather than subspecies (Churkin, 2006).

82

K. Omoto et al. / Gene 441 (2009) 80–88

Table 1 Nucleotide sequences of ND5 of butterflies used in this study No.

Species

Parnassiinae Genus Parnassius 1 P. acco 2 P. acco 3 P. acco 4 P. acdestis 5 P. actius 6 P. andreji 7 P. Apollo 8 P. Apollo 9 P. Apollonius 10 P. arcticus 11 P. ariadne 12 P. augustus 13 P. autocrator 14 P. boedromius 15 P. bremeri 16 P. cardinal 17 P. cephalus 18 P. charltonius 19 P. clodius 20 P. davydovi 21 P. delphius 22 P. delphius 23 P. epaphus 24 P. eversmanni 25 P. eversmanni 26 P. eversmanni 27 P. eversmanni 28 P. glacialis 29 P. glacialis 30 P. hardwickii 31 P. hide 32 P. honrathi 33 P. huberi 34 P. hunningtoni 35 P. imperator 36 P. inopinatus 37 P. jacquemontii 38 P. loxias 39 P. maharaja 40 P. maharaja 41 P. maximinus 42 P. mnemosyne 43 P. nomion 44 P. nordmanni 45 P. Orleans 46 P. patricius 47 P. phoebus 48 P. ruckbeili 49 P. schultei 50 P. simo 51 P. simonies 52 P. smintheus 53 P. staudingeri 54 P. stenosemus 55 P. stoliczkanus 56 P. stubbendorfii 57 P. stubbendorfii 58 P. szechenyii 59 P. tenedius 60 P. tianschanicus Genus Hypermnestra 61 H. helios Genus Zerynthia 62 Z. rumina 63 Z. polyxena Genus Archon 64 A. apollinus 65 A. apollinaris 66 A. apollinaris Genus Allancastria 67 A. cerisyi 68 A. caucasica

Subspecies

Locality

Accession number

Reference

gemmifer tulaishani prezwalskii lampidius actius

schultei zarraensis simonies smintheus Chitralica mulkilensis stoliczkanus hoenei kosterini szechenyii Britae insignis

Karo-la, S. Tibet, China Qilianshan, Gansu, China Tibet, China Karo-la, S. Tibet, China Tianshan, Xinjiang, China Sichuan, China Ural, South Russia Mt. Ventaux, S. France Kazakhstan N.E. Yakutia,Russia Altai, Russia Tibet, China Tajikistan Tianshan, Xinjiang, China Middle Amur, F.E. Russia Tajikistan Qamdo, E. Tibet, China Kyrgyzstan Montana, USA Kyrgyzstan Kyrgyzstan Tianshan, Xinjiang, China Qilianshan, Gansu, China Hokkaido, Japan Magadan, Russia Middle Amur, F.E. Russia Alaska, USA Hokkaido, Japan Sichuan, China E. Nepal Wenquan, Qinghai, China Uzbekistan Tanggulashan, C. Tibet, China Karo-la, S. Tibet, China Qinghai, China Afghanistan Alai, Kyrgyzstan Kyrgyzstan Qinghai, China Ladakh, India Tianshan, Xinjiang, China Kyrgyzstan Primorye, F.E. Russia Dombay, N. Caucasus, Russia Qinghai, China Tianshan, Xinjiang, China Magadan, F.E. Russia N. Hami, Xinjiang, China Karo-la, Tibet, China Ladakh, India Kyrgyzstan Colorado, USA Shandur, P., N.W. Pakistan Baralacha P., Ladakh, N.W. India Ladakh, N.W. India Hokkaido, Japan Magadan, Russia Qinghai, China Artyk, N.E. Yakutia, Russia Alai, Kyrgyzstan

AB095617 AB095618 AB095652 AB095621 AB095622 AB095643 AB095635 AB095636 AB095646 AB095639 AB094970 AB095645 AB095634 AB095629 AB095611 AB095644 AB095616 AB095630 AB095624 AB440273 AB095632 AB095655 AB095610 AB094971 AB095607 AB095608 AB096089 AB095627 AB095628 AB094969 AB095613 AB096091 AB095631 AB095633 AB095612 AB095641 AB095647 AB096090 AB095614 AB095615 AB095651 AB095626 AB095609 AB094968 AB095623 AB095620 AB095654 AB095638 AB095619 AB095640 AB095649 AB095653 AB095637 AB095656 AB095650 AB095625 AB095657 AB095642 AB095658 AB095648

Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. this study Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al. Omoto et al.

helios

Uzbekistan

AB095659

Omoto et al. (2004)

Rumina

Aix-en-Provence, France District of Voronezh, Russia

AB095660 DQ351054

Omoto et al. (2004) Nazari et al. (2007)

apollinus apollinaris bostanchii

Izmir, Turkey Rijab, Kermanshah, Iran Ploedokhtar, Lorestan, Iran

AB095661 DQ351051 DQ351052

Omoto et al. (2004) Nazari et al. (2007) Nazari et al. (2007)

cerisyi

Izmir, Turkey Bolu Daglari, Bolu Pro., Turkey

AB095662 DQ351057

Omoto et al. (2004) Nazari et al. (2007)

uralensis venaissinus Apollonius

augustus boedromius bremeri elwesi romanovi

albulus Juldussica tsaidamensis daisetsuzanus eversmanni felderi thor glacialis Tsingtaua hardwickii Aksobhya honrathi

musagetus variabilis tashkorensis labeyriei Maharaja giganteus mandschuricus groumi xinjiangensis ochotskensis

(2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004) (2004)

K. Omoto et al. / Gene 441 (2009) 80–88

83

Table 1 (continued) No.

Species

Parnassiinae Genus Allancastria 69 A. cretica 70 A. louristana Genus Luehdorfia 71 L. japonica 72 L. puziloi 73 L. chinensis 74 L. longicaudata Genus Bhutanitis 75 B. thaidiana 76 B. mansfieldi 77 B. lidderdali Genus Sericinus 78 S. montela Papilionidae Genus Papilio 79 P. machaon

Subspecies

yessoensis

Hippocrates

Locality

Accession number

Reference

Lassithi, Crete Island, Greece Malavi, Lorestan, Iran

DQ351056 DQ351055

Nazari et al. (2007) Nazari et al. (2007)

Kyoto, Japan Hokkaido, Japan

AB095663 AB095664 AB016826 AB016828

Omoto et al. (2004) Omoto et al. (2004) Makita et al. (2000) Makita et al. (2000)

Dongchuan, Yunnan, China

AB026728 AB026727 DQ351053

Makita et al. (2000)⁎ Makita et al. (2000)⁎ Nazari et al. (2007)

Kyoto, Japan

AB095665

Omoto et al. (2004)

Kyoto, Japan

AB095666

Omoto et al. (2004)

⁎ Published only in database.

Since P. loxias also occurs in Kyrgyzstan, but not sympatric to P. davydovi, one may wonder that the latter is only a subspecies of the former species in view of the fact that there are ample examples of morphological characters, such as wing pattern, leading to an erroneous classification. The present study, however, supports that they are in fact “good” species from each other, judging from the relatively long branch lengths from the common node to the two OTUs in the molecular phylogenetic tree shown in Fig. 2. As already proposed by Omoto et al. (2004) and confirmed by Michel et al. (2008), the butterflies of the genus Parnassius are made up of eight cluster groups (clades), or species groups, using molecular data. In general, they correspond to the subgenera hitherto proposed by various “Parnassiologists” as listed in Table 2. While Michel et al. (2008) use the subgeneric names to denote the cluster groups, we use the clade number system since only a part of them is discussed in the present study. 3.2. Diversity of the sphragis form in clade VI One of the important key characters for “Parnassiologists” in classifying the subgenus is the form of sphragis, the attachment to the end of the female abdomen which has been formed by the male secretion during copulation. As shown in Table 2, a nearly unique type of sphragis is found to occur in all clades except clade VI, to which the present target species P. davydovi belongs. A total of eight species are involved in this clade: P. acdestis, P. imperator, P. augustus, P. loxias, P. davydovi, P. autocrator, P. charltonius and P. inopinatus. Morphologically, all of these are relatively big, strongly built and brilliantly marked butterflies, except for P. acdestis, which is rather small, resembling in the wing pattern, as well as in the form of sphragis, the butterflies of clade V (see Table 2 for a list of different sphragis forms). Four types of sphragis are distinguished among the female butterflies of clade VI as seen in Fig. 3. The first type has two small brownish horns, possessed by P. acdestis. This sphragis type is characteristic of the female butterflies of clade V. The second sphragis type is possessed by P. imperator and P. augustus. It is a big, brownish fan-like form, with a rather irregular contour. Because of this peculiar sphragis, together with the wing pattern characteristics, these two species were given the subgeneric name Eukoramius. The third type of sphragis is big, light-brown, and rolled like a snail's shell. It is possessed by P. charltonius and, in a reduced and simplified form, by P. inopinatus. Another subgeneric name, Kailasius, was proposed for these two species, in addition to P. loxias and P. autocrator, although the female of the latter two species are characterized by the fourth type of sphragis, a big, yellowish square form and laterally flat. The

female P. davydovi was found to carry a sphragis of this type. The similarity in the form of sphragis found among P. loxias, P. davydovi and P. autocrator is in agreement with the molecular proximity between these three species as seen in Fig. 3. In Fig. 3, we indicate the time of emergence of the different sphragis types during the evolution of the eight species in clade VI. P. acdestis has been classified under subgenus “Koramius” together with the species of clade V. This “subgenus” is a paraphyletic group (Fig. 2). Therefore, morphological characters seen in this subgenus can be recognized as primitive characters (symplesiomorphy), and those seen in clade VI (except for P. acdestis) can be recognized as derived characters (apomorphy). It is evident that P. acdesis has inherited the ancestral morphological characters from “Koramius”, judging from the similar wing pattern and sphragis (Table 2). Later, three groups, P. imperator/ augustus, P. loxias/davydovi/autocrator, and P. charltonius/inopinatus, with markedly different sphragis types, emerged independently. The biological function of the sphragis is regarded as preventing secondary insemination after copulation, but the significance of the marked variation in the sphragis form, particularly among the butterflies of cluster VI, is unknown. 3.3. Statistical tests for some monophyly hypotheses Omoto et al. (2004) firstly indicated that the tribe Parnassiini (Parnassius, Hypermnestra and Archon) is not a monophyletic group. This relation was confirmed by Katoh et al. (2005), Nazari et al. (2007) and Michel et al. (2008). In this study, we also confirmed this relation (Fig. 2). This relation for a polyphyletic hypothesis was statistically supported by the KH test (Kishino and Hasegawa, 1989) and the AU test (Shimodaira, 2002) (Table 3). Moreover, the monophyletic relationship between Bhutanitis and Luehdorfia was suggested from morphological studies (Hiura 1980, Hancock 1983, Nazari et al. 2007). However, this relation was not supported from recent molecular phylogenetic studies (Nazari et al., 2007, Michel et al., 2008). Our study also failed to support Bhutanitis and Luehdorfia monophyly (Fig. 2). The ML tree (Bhutanitis and Luehdorfia polyphyly) was significantly supported by the KH test and the AU test (Table 3). The phylogenetic position of P. davydovi among genus Parnassius was also enigmatic before our molecular systematic analyses. Based on the similarity in external morphological characters between P. davydovi and P. charltonius, one may wonder that a monophyletic relationship exists between the two species. However, the ML tree in this study clearly shows that P. davydovi and P. loxias are monophyletic (Table 3).

84

K. Omoto et al. / Gene 441 (2009) 80–88

Fig. 2. The ML tree as inferred from the nucleotide sequence data of ND5 (807 bp) with GTR+I+Γ model. The nodal numbers indicate the BPs and those larger than 50% are shown. The clades, denoted with roman numbers (I–VIII), correspond to those presented by Omoto et al. (2004). The branch lengths are proportional to the estimated numbers of the nucleotide substitution (see scale bar).

3.4. The estimation of divergence times The estimated divergence times among the subfamily Parnassiinae are shown in Fig. 4. Although we used nearly the same data set and the same calibration points as Nazari et al. (2007), we obtained much younger divergence times than their estimates. Interestingly, when we estimated divergence times with concatenated data instead of partitioning among genes and codon positions, they became much older. For example, the divergence time between Parnassius and Hypermnestra was estimated to be 35.0+/-5.8 MYA BP from the partitioned data, and it was 50.3 ± 5.5 MYA BP based on the concatenated data. The latter value is similar to the result of Nazari et al. (2007, Table 5 and Fig. 14). Divergence time estimations directly affected the estimated branch lengths and it must be taken into consideration that, if

multiple substitutions occur, branch lengths will be under-estimated. The concatenated model does not take into account the difference of the mode and the ratio of the substitutions among each gene. Alternatively, the partitioned model seems to be more realistic, since it takes into account these differences, and is relatively more effective in evaluating the correct numbers of the multiple substitutions. Actually, by using the AIC (Akaike Information Criterion; Akaike, 1973), it is indicated that the partitioned model fit to the data better than the concatenated model, the smaller the AIC, the better the model. In this case, the AIC of the partitioned model was (44,041.0) and the AIC of the concatenated model was (47,116.4). Igawa et al. (2008) demonstrated that if branch lengths are underestimated and the calibration point set in the deep node, the divergence time may be overestimated. When genes with a high evolutionary rate are used for the estimation, this tendency

K. Omoto et al. / Gene 441 (2009) 80–88

85

Table 2 Classification of the eight Parnassian clades and the different sphragis forms Clade

Species group

Sphragis form

Subgenus

I II III IV V VI

The The The The The The The The The The

1 2 3 4 5 6 7 8 9 10

Parnassius s.str. Latreille 1804 Lingamius Bryk 1935 Tadumia Moore 1902 Sachaia Korshunov 1988 Koramius Moore 1902 Eukoramius Bryk 1934 Kailasius Moore 1902 Ditto Kreizbergia Korshunov 1990 Driopa Korshunov 1988

VII VIII

Apollo group Hardwickei group Acco group Tenedius group Delphius group Imperator group Charltonius group Loxias group Simo group Mnemosyne group

The eight clades (species groups) of Parnassius proposed by the molecular phylogenetic tree (Omoto et al., 2004) and the corresponding subgenera proposed previously by various “Parnassiologists” based on morphological characteristics, including the form of sphragis. Note that only in clade VI multiple forms of sphragis are found. Explanation of the sphragis forms. 1. Small, dark brown, and sharply pointed. 2. Large, whitish pouch laterally flattened. 3. Large, whitish pouch laterally flattened, completely surrounding the abdomen. 4. Small, white and formed by two or three bent strips. 5. Small, brownish, with two horns. 6. Large, brownish, flattened dorso-ventrally into a fan shape, with the outermost extremities bent upwards. 7. In P. charltonius, it is large, creamy, and rolled back on itself like a snail's shell. In P. inopinatus, it is diminished and flattened dorso-ventrally, indented at the back. 8. In P. loxias and P. davydovi, it is medium sized, light grey, and flattened laterally. In P. autocrator, it is large, lemon yellow, and laterally flattened. 9. Mostly absent. 10. Large, white or creamy or brownish pouch with a hole caudally located.

Fig. 3. Phylogenetic relationships among clade VI (magnified view of Fig. 2). The nodal numbers indicate the BPs of the ML analysis (above branches) and the posterior probabilities of the Bayesian analysis (below branches). The branch lengths are proportional to the numbers of the nucleotide substitution estimated by ML (see scale bar).

86

K. Omoto et al. / Gene 441 (2009) 80–88

Table 3 Comparison of the various phylogenetic hypotheses of Parnassinae

ML tree tribe Parnassiini monophyly Luehdorfia +Bhutanitis monophyly P.davydovi +P.charltonius monophyly

Log-likelihood score

KH test

AU test

b− 7791.126N − 15.17 ± 8.77 − 33.38 ± 11.67 − 13.84 ± 8.38

– 0.042 0.002 0.049

– 0.026 0.002 0.034

The results of the comparison of the various phylogenetic hypotheses of Parnassinae based on the nucleotide sequence data of ND5 (807 bp) with the GTR+Γ model. The log-likelihood of the ML tree is given in angled brackets, and the differences in the log− likelihoods of alternative trees from that of the ML tree are given with 1SE(after±). The numbers of the KH test (Kishino and Hasegawa, 1989) and the AU test (Shimodaira, 2002) indicate p values. The calculations of the log-likelihood scores and the KH test were performed by using PAML ver 4.0 (Yang, 2007) and the AU tests were performed by using CONSEL (Shimodaira and Hasegawa, 2001).

will be remarkable. Michel et al. (2008) used only the mitochondrial genes and estimated divergence times with a concatenated model. They set the calibration point on the root of the tree (they assumed the split of Parnassiinae and Papilioninae as 100 MYA BP). Their estimated divergence times were much older than ours. They seem to be overestimated, due to the reason mentioned above. The time of the relatively rapid radiation of different species groups in the genus Parnassius was estimated at 24.3 MYA BP (Fig. 4). This is regarded approximately as the period that the uplift of the Tibetan Plateau accelerated (Zachos et al., 2001). It may be postulated that the rapid geological change associated with climatic

changes may have been a driving factor for the creation of new lineages of Parnassius, which is essentially adapted to the cold climate. 3.5. Divergence times within the genus Parnassius We tested the molecular clock hypothesis among species of the genus Parnassius with the likelihood ratio test, and it was not rejected (p = 0.49). Therefore, divergence times among Parnassian species are regarded as following the molecular clock (Fig. 5). When the date of the rapid dispersal event of the eight major Parnassian≐clades (node 1 in Fig. 5) was calibrated, 24.3 MYA BP, the split of P. davydovi and P. loxias was estimated to be 10.8 MYA BP. Taking into account the 95% confidence interval of the divergence time estimation, the range was estimated to be 7.5–14.7 MYA BP. When the split date between clade V and clade VI (node 2 in Fig. 5) was calibrated at 19.5 MYA BP, the split of P. davydovi and P. loxias was estimated to be 10.0 MYA BP. Taking into account the 95% confidence interval of the divergence time estimation, the range was estimated to be 6.5–14.2 MYA BP. From these findings, the split between P. davydovi and P. loxias was estimated to be approximately 10 (6.5–14.7) MYA BP. This age is much older than the ages of the splits between subspecies among the genus Parnassius (e.g., P. apollo, P. acco, P. maharaja, P. delphius, P. stubbendorfii, P. glacialis and P. eversmanni) as indicated in Fig. 5. Although only one individual of P. davydovi was used in our analyses, it appears to be clear that P. davydovi is not a subspecies of P. loxias but a distinct, “good” species.

Fig. 4. The estimated divergence times of the subfamily Parnassinae as inferred from the nucleotide sequence data of ND1, ND5, COI, COII, 16SrRNA, Leu-tRNA, and EF1a without assuming the constant rate of the molecular clock. The nodal numbers indicate the estimated divergence times with 1SE (after±) in MYA BP. For conducting MCMC, we set sample numbers to be 10,000, sample frequency to be 100, and burnin to be 100,000. We assumed the monophyly of three major clades of Parnassinae: (1) Parnassius+Hypermnestra, (2) Archon+Leuhdorfia, and (3) Allancastria+Zerynthia+Bhutanitis+Sericinus, respectively (Omoto et al., 2004, Katoh et al., 2005, Nazari et al., 2007, this study) and also assumed monophyly of Allancastria+Zerynthia (Nazari et al., 2007, this study). Abbreviations: Pa, Paleocene; E, Eocene; O, Oligocene; M, Miocene; Pli, Pliocene; Ple, Pleistocene.

K. Omoto et al. / Gene 441 (2009) 80–88

87

Fig. 5. The estimated divergence times among the genus Parnassius as inferred from the nucleotide sequence data of ND5 assuming the molecular clock. The phylogenetic relationships among the eight clades were treated as multifurcation. In this analysis, we treated clade VIII as a monophtyletic group. Although the monophyly of clade VIII was not supported from the ND5 data (Fig. 2), it was supported from the ND1, ND5, COI and 16SrRNA data of Michel et al. (2008). Hypermnestra helios was used as an out-group. The nodes with the circled digits were used as calibration points for the estimation of the divergence time among the genus Parnassius.

Further studies including the large data set with respect to number of genes, both mitochondrial and nuclear, are in progress for elucidating the detailed phylogenetic relationships among the species of the genus Parnassius, as well as among higher taxa of the subfamily Parnassiinae. Acknowledgments We wish to express our special thanks to Dr. Masami Hasegawa and two anonymous referees for the valuable comments. We thank Mr. Yuichi Kawasaki who played an intermediary role for one of us (K. O.) in obtaining the specimen of P. davydovi from Mr. Sergei Churkin. Thanks are also to Mr. Akio Masui who helped us in preparing Fig. 1 of this paper. Invaluable support by Drs. Kazuo Umetu, Tadao

Matsumoto, and Masaru Nonaka are also gratefully appreciated. We thank Mr. Michael Kryshak for reading the manuscript. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.gene.2008.10.030.

References Akaike, H., 1973. A new look at the statistical model identification. IEEE Trans. Autom. Contr. AC 19, 716–723. Churkin, S., 2006. A new species of Parnassius Latreille, 1804, from Kyrgyzstan (Lepidoptera, Papilionidae). Helios (Moskow) 7, 142–158.

88

K. Omoto et al. / Gene 441 (2009) 80–88

Felsenstein, J., 1981. Evolutionary trees from DNA sequences — a maximum likelihood approach. J. Mol. Evol. 17, 368–376. Gaunt, M.W., Miles, M.A., 2002. An insect molecular clock dates the origin of the insects and accords with palaeontological and biogeographic landmarks. Mol. Biol. Evol. 19, 748–761. Hancock, D.L., 1983. Classification of the Papilionidae (Lepidoptera): a phylogenetic approach. Smithersia 2, 1–48. Hasegawa, M., Kishino, H., Yano, T., 1985. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22, 160–174. Häuser, C.L., 1993. Critical comments on the phylogenetic relationships within the family Papilionidae. Nota Lepidopterol. 16, 34–43. Hiura, I., 1980. A phylogeny of the genera of Parnassiinae based on analysis of wing pattern, with description of a new genus (Lepidoptera: Papilionidae). Bull. Osaka Mus. Nat. Hist. 33, 71–85. Igawa, T., Kurabayashi, A., Usuki, C., Fujii, T., Sumida, M., 2008. Complete mitochondrial genomes of three neobatrachian anurans: a case study of divergence time estimation using different data and calibration settings. Gene 407, 116–129. Katoh, T., Chichvarkhin, A., Yagi, T., Omoto, K., 2005. Phylogeny and evolution of butterflies of the genus Parnassius: inferences from mitochondrial 16S and ND1 sequences. Zool. Sci. 22, 343–351. Kishino, H., Hasegawa, M., 1989. Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in Hominoidea. J. Mol. Evol. 29, 170–179. Makita, H., Shinkawa, T., Kazumasa, O., Kondo, A., Nakazawa, T., 2000. Phylogeny of Luehdorfia butterflies inferred from mitochondrial ND5 gene sequences. ENT Sci. 3, 321–329. Michel, F., Rebourg, C., Cosson, E., Descimon, H., 2008. Molecular phylogeny of Parnassiinae butterflies (Lepidoptera: Papilionidae) based on the sequences of four mitochondrial DNA segments. Ann. Soc. Entomol. Fr. (n. s.) 44, 1–36. Nazari, V., Zakharov, E.V., Sperling, F.A.H., 2007. Phylogeny, historical biogeogeography, and taxonomic ranking of Parnassiinae (Lepidoptera, Papilionidae) based on morphology and seven genes. Mol. Phylogenet. Evol. 42, 131–156. Omoto, K., Katoh, T., Chichvarkhin, A., Yagi, T., 2004. Molecular systematics and evolution of the “Apollo” butterflies of the genus Parnassius (Lepidoptera: Papilionidae) based on mitochondrial DNA sequence data. Gene 326, 141–147. Rodriguez, F.J., Oliver, J.L., Marin, A., Medina, J.R., 1990. The general stochastic model of nucleotide substitution. J. Theor. Biol. 142, 485–501. Ronquist, F., Huelsenbeck, J.P., 2003. MRBAYES 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574.

Saitou, N., Nei, M., 1987. The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol. Biol. Evol. 4, 406–425. Shimodaira, H., Hasegawa, M., 2001. CONSEL: for assessing the confidence of phylogenetic tree selection. Bioinformatics 17, 1246–1247. Shimodaira, H., 2002. An approximately unbiased test of phylogenetic tree selection. Syst. Biol. 51, 492–508. Stamatakis, A., Hoover, P., Rougemont, J., 2008. A Rapid Bootstrap Algorithm for the RAxML Web-Servers. Syst. Biol. 57, 758–771. Tamura, K., Dudley, J., Nei, M., Kumar, S., 2007. MEGA4: Molecular Evolutionary Genetics Analysis (MEGA) software version 4.0. Mol. Biol. Evol. 24, 1596–1599. Tamura, K., Nei, M., Kumar, S., 2004. Prospects for inferring very large phylogenies by using the neighbor-joining method. Proc. Natl. Acad. Sci. 101, 11030–11035. Thompson, J.D., Higgins, D.G., Gibson, T.J., 1994. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignments through sequence weighting, position specific gap penalties and weight matrix choice. Nucl. Acids Res. 22, 4673–4680. Thorne, J., Kishino, H., Painter, I., 1998. Estimating the rate of evolution of the rate of molecular evolution. Mol. Biol. Evol. 15, 1647–1657. Yagi, T., Sasaki, G., Takebe, H., 1999. Phylogeny of Japanese Papilionid butterflies inferred from nucleotide sequences of the mitochondrial ND5 gene. J. Mol. Evol. 48, 42–48. Yagi, T., Katoh, T., Chichvarkhin, A., Shinkawa, T., Omoto, K., 2001. Molecular phylogeny of butterflies Parnassius glacialis and P. stubbendorfii at various localities in East Asia. Genes & Genet. Syst. 76, 229–234. Yagi, T., Sasaki, G., Omoto, K., 2006. Phylogeny, biogeography and wing pattern convergence of the subgenus Achillides (Lepidoptera, Papilionidae, Papilio) revealed by a mitochondrial DNA sequence analysis. Trans. Lepid Soc. Japan 57, 137–147. Yang, Z., 1996. Among-site rate variation and its impact on phylogenetic analyses. TREE 11, 367–372. Yang, Z., Rannala, B., 1997. Bayesian phylogenetic inference using DNA sequences: a Markov chain Monte Carlo method. Mol. Biol. Evol. 14, 717–724. Yang, Z., 2007. PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586–1591. Zachos, J., Pagani, M., Sloan, L., Thomas, E., Billups, K., 2001. Trends, rhythms, and aberrations in global climate 65 Ma to present. Science 292, 686–693. Zakharov, E.V., Caterino, M.S., Sperling, F.A., 2004. Molecular phylogeny, historical biogeography, and divergence time estimates for swallowtail butterflies of the genus Papilio (Lepidoptera: Papilionidae). Syst. Biol. 53, 193–215.