Although the taxonomy of the ca 18 000 species of butterflies and skippers is well known, the family-level relationships are still debated. Here, we present, to our knowledge, the most comprehensive phylogenetic analysis of the superfamilies Papilionoidea, Hesperioidea and Hedyloidea to date based on morphological and molecular data. We reconstructed their phylogenetic relationships using parsimony and Bayesian approaches. We estimated times and rates of diversification along lineages in order to reconstruct their evolutionary history. Our results suggest that the butterflies, as traditionally understood, are paraphyletic, with Papilionidae being the sister-group to Hesperioidea, Hedyloidea and all other butterflies. Hence, the families in the current three superfamilies should be placed in a single superfamily Papilionoidea. In addition, we find that Hedylidae is sister to Hesperiidae, and this novel relationship is supported by two morphological characters. The families diverged in the Early Cretaceous but diversified after the Cretaceous–Palaeogene event. The diversification of butterflies is characterized by a slow speciation rate in the lineage leading to Baronia brevicornis, a period of stasis by the skippers after divergence and a burst of diversification in the lineages leading to Nymphalidae, Riodinidae and Lycaenidae.
With 18 363 described species , butterflies and skippers are a highly diverse group with a species-level taxonomy that is the best known among invertebrates. However, the evolutionary history of the group is not well known, including when the group evolved, and how it diversified, although recent studies on subgroups within the butterflies are suggesting intriguing hypotheses about where, when and how the clades originated and diversified [2–6].
Despite being one of the most familiar and best-studied groups of insects, the higher phylogeny of the true butterflies (superfamily Papilionoidea) has been a tough nut to crack. The family divisions of this alluring superfamily are well established and traditionally agreed upon to include the families Papilionidae, Pieridae, Nymphalidae, Lycaenidae  and recently Riodinidae [8,9]. It is also rather apparent that the family Hesperiidae (skippers) and Hedylidae (butterfly moths) are affiliated with the butterflies [10,11]. However, the inter-relationships of these families have long been contradictory, and results of even recent studies have not made the matter any easier [9,12–15]. Solving this systematic puzzle would allow better understanding of the evolutionary history of these popular organisms that serve as biological models in many fields of research.
The inability of molecular and morphological data to resolve phylogenetic relationships is often attributed to ancient, rapid diversification of lineages [16–18]. Such patterns are important signals of evolutionary history: could they be linked with geological or climatic events in the Earth's history such as the Cretaceous–Palaeogene (K–Pg) mass extinction [3,6]? Or perhaps, are they linked to biological innovations, such as colonizing novel host plants ? Knowledge of when major lineages have diverged would help in understanding patterns of diversification and discovering the links to the events that were responsible for the current diversity of butterflies and skippers.
Through time-calibrated phylogenies, it is also possible to investigate whether the high number of extant groups of butterflies and skippers is the result of rapid diversifications over a short period of time, or whether the numbers we see now are the expected numbers for an old lineage with even diversification rates (since a constant rate of diversification over long periods will certainly produce a very species-rich lineage) . This can be achieved by estimating past rates of speciation (birth: b) and extinction (death: d) along a phylogeny and for the lineages of interest. Estimating a phylogenetic hypothesis for all the ca 18 000 species of butterflies and skippers is impractical, and thus methods that are able to correct for incomplete phylogenies are needed (taking into account, the species richness of lineages) . Such methods have been recently developed and allow us to understand the patterns of ancient diversifications that are inferred based on estimated dates of diversifications.
The aim of this paper is to perform, to our knowledge, the most comprehensive phylogenetic analysis of the superfamilies Papilionoidea, Hesperioidea and Hedyloidea to date, in order to assess the relationships among their constituent families. We also investigated whether any lineage underwent unusually rapid diversification compared with the net diversification speed (r = b – d) across the phylogeny of butterflies and skippers. We present a morphological data matrix, which includes extensive data: 191 characters derived from adults, larvae and pupae, as well as a molecular data matrix that comprised eight protein-coding gene regions (one from the mitochondrial genome and seven from the nuclear genome). We use this data to investigate the evolutionary history of the butterflies and skippers.
(a) Molecular data sampling
Taxa for this study were selected based on previous analyses of families [3,4,6,9,20] and covered all subfamilies in each family (except Pseudergolinae in Nymphalidae and Aphnaeinae in Lycaenidae; electronic supplementary material, table S1). In a few cases (especially for Lycaenidae), specimens were not available for sequencing and thus previously published sequences were used . Twenty-five outgroups were initially chosen (electronic supplementary material, table S1) based on a recent study on all Lepidoptera , but for most analyses, the closest outgroups, i.e. five species of Thyrididae and two species of Callidulidae, were used. We sequenced one mitochondrial (COI) and up to seven protein-coding nuclear gene regions (EF-1α, Wingless, RpS5, MDH, GAPDH, CAD and IDH) accounting for a total of 6165 bp with gaps. PCR conditions and sequencing have been described in detail previously . Sequencing was performed mainly with an ABI 3730 capillary sequencer (Oulu) and an ABI PRISMR 3130xl capillary sequencer (Turku).
The morphological dataset includes characters derived from the scale vestiture and sclerotized structures of Papilionoidea (electronic supplementary material, tables S2 and S3). Characters were sought from adults, larvae and pupae by studying actual specimens from various museum collections, and supplemented by re-examining characters used in previous studies . For the study of adult morphology, dissections prepared from dry mounted specimens were used. The study of the larval characters was restricted to the last instar larva, some preserved in ethanol and others dry-inflated. Pupal characters were examined, whenever possible, from both entire preserved pupae and from exuviae from which the adult had emerged. Detailed methods used in the study of morphology are given in electronic supplementary material, table S2.
In total, 54 ingroup species were examined, mostly corresponding to the species that were sequenced (electronic supplementary material, table S1). Both male and female specimens were available for most taxa, and either larva, pupa or both available for 36 of these species. Of the outgroup, larvae and pupae were obtained only for Scardia boletella and Pterodecta felderi. The character list comprises 191 characters (143 binary states, 48 multi-states; electronic supplementary material, table S2). Of these characters, 45 are larval, 32 pupal and 114 are adult characters. Reasoning for major modifications in traditional character coding is reported in the electronic supplementary material, table S2. The data of the character coding are presented in the matrix (electronic supplementary material, table S3).
(c) Phylogenetic analyses
The data were analysed using maximum parsimony (MP), maximum likelihood (ML) and Bayesian inference (BI). MP and ML analyses were run on individual gene regions, combined molecular data, all data, molecular data with third codon positions removed and on a reduced dataset containing only taxa with six or more gene regions sequenced. MP was performed using Tree analysis using New Technology (TNT) . Morphological characters were treated as unordered  gaps in molecular characters as a fifth state. Bootstrap and Bremer support values for nodes were obtained using a script written for TNT . ML was performed with RAxML , with 100–500 bootstrap replicates. Each gene was assigned the GTR + Γ model, based on Akaike information criterion (AIC) values calculated using FindModel .
BI was performed with MrBayes v. 3.1  on the combined dataset. Missing nucleotides were coded as ‘?’. Parameter values were estimated separately for each gene region and the morphological data using the ‘unlink’ command and the rate prior (ratepr) set to ‘variable’. As in the ML analyses, each gene was assigned the GTR + Γ model. The morphological data were modelled under the rate variable model .
Bayesian inference of phylogeny and times of divergence were carried out using the program BEAST v. 1.5.4  using only the ingroup taxa. The dataset was partitioned into morphology and the gene regions and analysed simultaneously under the GTR + Γ model for each gene partition separately, and with a relaxed clock allowing branch lengths to vary according to an uncorrelated lognormal distribution . The tree prior was set to the birth–death process. Three fossils were used as calibration points. The split between Hypanartia and Vanessa (Nymphalidae) was constrained to be 34 Ma (s.d. ± 3 Ma) based on two Florissant fossils Prodryas persephone and Vanessa amerindica  as well as results from previous studies . Similarly, two other Florissant fossils, Stolopsyche libytheoides and Oligodonta florissantensis [31,32], were used to constrain the age of the first divergence in Pierini (Pieridae) in a similar manner. The third calibration point was taken from a previous study , where the divergence between Parnassiinae and Papilioninae was constrained to be between 50 and 80 Ma (uniform distribution). Initial analyses showed that there were serious problems with inferring relationships (the topology was found to be inverted with the root going between Lycaenidae and Riodinidae); thus we placed two constraints on the topology based on the combined analyses using TNT and MrBayes, i.e. all except Papilionidae were constrained to form a monophyletic group and Pieridae, Nymphalidae, Lycaenidae and Riodinidae were constrained to be a monophyletic group. All other priors were left to the defaults in BEAST. Parameters were estimated using four independent runs of 10 million generations each, with parameters sampled every 1000 generations. Convergence was checked in Tracer v. 1.4.6 and summary trees were generated using TreeAnnotator v. 1.5.3, both part of the BEAST package.
(d) Rates of diversification
We analysed patterns of diversification along the tree of butterflies and skippers using algorithms for modelling evolutionary diversification using stepwise AIC (MEDUSA) as implemented in R [19,33] in the package Geiger . MEDUSA fits alternative birth–death models to a phylogenetic tree taking into account taxonomic data in order to estimate changes in net diversification rates along the branches. MEDUSA obtains values of b and d based on the best-fit of the models according to comparisons of AIC scores . First, the MEDUSA algorithm estimates likelihood and AIC scores for the simplest birth–death model (with two parameters, b and d). Next, MEDUSA compares the AIC scores of the two-parameter model with an incrementally more complex model, which includes five parameters (for two birth parameters and two death parameters owing to one break point and a shift-location parameter). MEDUSA continues this process until the addition of parameters does not improve the AIC scores.
We obtained information about species numbers for major lineages of butterflies and skippers from several sources [1,35]. This information was assigned to 28 lineages in our phylogeny after pruning terminals belonging to the same monophyletic groups (electronic supplementary material, table S4). We let MEDUSA estimate up to 10 turnover points in our tree but present the pattern of diversification obtained using a moderate corrected AIC threshold (cut-off value = 4) . MEDUSA produces corrected AIC scores (AICc) , which is a bias correction to the AIC when sample size is small. We selected the model with the lowest AICc score, which in this case was the same when using a moderate cut-off value 4 for normal AIC scores, which means rejecting models that do not improve the AIC scores by more than 4 .
(a) Systematic implications
Two major results are apparent from our analyses: the paraphyly of butterflies, with Papilionidae being sister to the rest of the butterflies, skippers and butterfly moths, and the sister-group relationship between Hedylidae and Hesperiidae (figure 1 and electronic supplementary material, figure S1). We recovered the position of Papilionidae as sister-group to Hedylidae + Hesperiidae + the rest of butterflies in all the analyses except the parsimony analysis of morphology alone, where Hedylidae, Hesperiidae and the rest of Papilionoidea form an unresolved trichotomy (electronic supplementary material, figure S1a).
The position of Pieridae was highly unstable, it being found to be the sister-group to Lycaenidae + Riodinidae (figure 1b and electronic supplementary material, figure S1c,e), to Nymphalidae + Lycaenidae + Riodinidae (figure 1a and electronic supplementary material, figure S1b,f), to Papilionidae (electronic supplementary material, figure S1a) and finally to the clade excluding Papilionidae (electronic supplementary material, figure S1d), depending on the data and the method of analysis. All the analyses found Riodinidae and Lycaenidae to be sister lineages.
(b) Ages of lineages and rates of diversification
The analysis of times of divergence suggests that the lineages leading to Papilionidae, Hesperiidae + Hedylidae and the rest of the butterflies diverged quite rapidly from each other in the Early Cretaceous, some 110 million years ago (Ma; figure 2). Lineages leading to extant families all diverged from each other by 90 Ma, with Pieridae diverging from the common ancestor to the other butterflies at about 105 Ma, nymphalids from lycaenids and riodinids about 102 Ma, hedylids diverged from hesperiids about 99 Ma and finally riodinids diverging from lycaenids about 88 Ma.
These initial divergences appear to be followed by a period of no family-level divergences of extant lineages (figure 2). The first family that began diversifying was Nymphalidae, with the initial divergence happening at about the same time as the divergence of the Riodinidae and Lycaenidae lineages, i.e. in the Late Cretaceous. Strikingly, most within-family divergences leading to extant subfamily lineages occurred after the K–Pg boundary (also known as the K–T or Cretaceous–Tertiary boundary), the exceptions being the aforementioned Nymphalidae  and possibly Pieridae .
The diversification analysis suggests that the current diversity of butterflies and skippers is best explained by three shifts in the rates of diversification during their evolutionary history (figure 2). The best-fit model had a log-likelihood value of −288.05 and 11 parameters (four birth rates, four death rates and three shift-location parameters), and it had the best-corrected AIC score (AICc = 604.24). The diversification analysis found the net rate of diversification (r) of butterflies and skippers to be 0.051 lineages per million years (table 1). The lineage leading to the subfamily Baroniinae (Papilionidae) has strikingly low rates of speciation and extinction, leading to a very low net rate (r = 1.528 × 10e−16). Another decrease in the rate of speciation is inferred to have happened in Hesperiidae (r = 0.0327). The clade including Nymphalidae, Lycaenidae and Riodinidae is inferred to have an elevated net rate of diversification (r = 0.081) compared with the background diversification rates in butterflies and skippers.
(a) Systematic implications
The idea that Papilionoidea do not form a natural group has never been entertained in the history of the classification of these well-studied insects. Only recently have two studies [13,14] called into question the monophyly of Papilionoidea, although both studies suffered from minimal taxon sampling of butterflies and close relatives. We increased taxon sampling of butterflies and skippers and inferred that Papilionoidea as currently understood is a paraphyletic entity, with Papilionidae being sister to skippers, hedylids and all other butterflies. In addition, we found strong support for a sister-group relationship between Hedylidae and Hesperiidae. In this study, the identical shape of the third axillary sclerite of the forewing (51 : 1 and 52 : 0; electronic supplementary material, figure S2a), formerly considered an apomorphy of Hesperiidae, and the flat ridges on the mesophragma (38 : 1; electronic supplementary material, figure S2b) are morphological characters that could support the sister-group relationship between Hedylidae and Hesperiidae.
These results differ radically from a similar previous study , although we note that the Bayesian analysis of three genes in the previous study finds Papilionidae to be sister to the rest, including Hedylidae and Hesperiidae (figure 1b in the study of Wahlberg et al. ). That hypothesis was discounted as it changed with the addition of new data . Our current study adds new data, both morphological and molecular, and evidently the new data show a stronger phylogenetic signal, as the positions of Papilionidae, Hedylidae and Hesperiidae appear robust to various perturbations of the data (figure 1 and electronic supplementary material, figure S1).
One of the major uncertainties in butterfly phylogeny has been the position of Pieridae . Two hypotheses have dominated: Pieridae and Papilionidae as sister-groups [37–39] or Pieridae sister to Lycaenidae (s.l.) + Nymphalidae [9,40,41]. The MP analysis of the morphological data alone suggests that Pieridae is the sister-group of Papilionidae, in accordance with one hypothesis [37–39]. However, this relationship is only supported by five homoplastic characters (electronic supplementary material, figure S2c) and no molecular studies have suggested the sister relationship of Pieridae and Papilionidae. Most other analyses in our study place Pieridae with Nymphalidae and Riodinidae + Lycaenidae.
The difficulty in resolving the relationships between the families of butterflies and skippers is evidently owing to the rapid divergences of the lineages leading to Papilionidae, Hedylidae + Hesperiidae, Pieridae, Nymphalidae and Riodinidae + Lycaenidae (figure 2). It appears that these divergences, happening within a ca 8 Myr period in the Early Cretaceous (110–102 Myr ago), may not have left any synapomorphies behind which we could have detected to be used to unite lineages. Given the error inherent in estimating divergences that have happened over 100 Ma, such a time span is essentially not different from a hard polytomy, despite the visual representation of dichotomously branching lineages (figure 2).
Our results suggest that the division of butterflies, skippers and butterfly moths into three different superfamilies is not tenable. We suggest that the families Papilionidae, Hedylidae, Hesperiidae, Pieridae, Nymphalidae, Riodinidae and Lycaenidae be placed in a single superfamily Papilionoidea. This clade appears to be stable within Lepidoptera, and the constituent families are strongly supported clades that are well-defined morphologically [12,14]. Any changes of the relationships of the families based on increased data would not affect the classification as suggested here.
(b) Rates of diversification and the evolutionary history of butterflies
The diversification analysis found the net rate of diversification of butterflies and skippers to be low (r = 0.0506 lineages per Myr). It is notable that the species-rich skipper family Hesperiidae (ca 4000 species ) has a slower rate of diversification (r = 0.0327) than the overall tempo of diversification across our phylogenetic tree. Our results indicate that the lineage leading to extant Hesperiidae diverged early in the evolution of butterflies (around 100 Ma) and had a relatively long delay before it began to diversify after 65 Myr. The lineage thus survived the mass extinction event at the K–Pg boundary. It appears that despite having a slower net rate of diversification, the species richness of skippers is relatively high owing to a steady rise in the number of lineages over a long span of time. Whether the skippers were diverse in the Cretaceous is not known, thus we are unable to say whether there was a long period of time during which the skipper lineage survived with low diversity, or whether simply one lineage survived the K–Pg event. A previous study on the subfamily Satyrinae (Nymphalidae) suggests that a radical change in the environment allowed the butterflies to exploit the new habitat that became available with the drying up and cooling down of the Earth in the Oligocene . Clearly, a more detailed study of the evolutionary history of Hesperiidae is needed to investigate the reasons behind the ‘long fuse’ before the diversification of the skippers.
A previous study proposed that 10 or 12 lineages of nymphalids survived the K–Pg event and radiated with elevated rates of diversification during the Tertiary diversification . In this paper, we find support for our previous findings, as our Nymphalidae + (Riodinidae + Lycaenidae) clade has a higher net diversification rate (r = 0.0814) than the background tempo for all the butterflies and skippers (r = 0.0506; figure 2). A more detailed diversification analysis including more terminals for the species-rich Nymphalidae + (Riodinidae + Lycaenidae) clades would shed light on the different rates of diversification for specific clades, as there are many lineages with very few species, such as the riodinid Stygini, that might have suffered a deceleration of their net diversification rate. Similar findings may be found by studying Satyrinae at a finer scale, which is thought to have undergone a rapid radiation .
The rapid diversification of the lineages leading to the families, and the apparent deceleration of diversification until after the K–Pg boundary can be explained in several ways. It may be that the butterflies and skippers were more diverse in the Cretaceous, and suffered from extinctions during the K–Pg event, as has been suggested for Nymphalidae . An alternative scenario would suggest that the ecological conditions for diversification of day-flying Lepidoptera only became favourable in the Palaeogene, e.g. in conjunction with the diversification of angiosperm plants , on which the majority of butterflies and skippers are specialized as larvae .
It is interesting to note that our findings of the tempo of diversification of butterflies and skippers are comparable with the rates found in the evolutionary history of vertebrates , which have high net rates of diversification for species-rich lineages that began diversifying simultaneously with the butterflies, around 100 Ma, and have comparable numbers of extant species. For example, Percomorpha fishes, which include around 20 000 extant species, have a net diversification rate of 0.082 lineages per million years. This is comparable with our estimated rate of 0.0814 lineages per million years for the clade Nymphalidae + (Riodinidae + Lycaenidae) that include ca 12 000 butterfly species.
In summary, we found Papilionoidea to be paraphyletic with regard to Hedyloidea and Hesperioidea, as well as strong support for a sister relationship between the butterfly–moths (Hedylidae) and the skippers (Hesperiidae). Therefore, as currently understood, Papilionoidea is a paraphyletic entity, and we suggest it should include Hedylidae and Hesperiidae. As in previous molecular studies, our study recovered the higher level relationships among the rest of butterflies to be (Pieridae + (Nymphalidae + (Riodinidae + Lycaenidae))). The butterfly families probably diverged from each other in the Early Cretaceous, between 110 and 100 Ma, but extant lineages began diversifying only after the K–Pg event at 65 Ma. We found that the pattern of the diversification of butterflies is characterized by three changes in the rates of speciation along their history. The species Baronia brevicornis should be regarded as a ‘living fossil’, as it is a descendant of an old lineage with an extremely low rate of diversification (r = 1.528 × 10e−16). The skippers had a period of stasis following their divergence from their sister-group Hedylidae, while an acceleration in the rate of diversification occurred in the surviving lineages of the K–Pg event that led to the clade Nymphalidae + (Riodinidae + Lycaenidae) (r = 0.0814).
We thank the following institutions and persons for providing material for the morphological study: Finnish Museum of Natural History (MZH), Australian National Insect Collection (ANIC), Texas A&M University, McGuire Centre for Lepidotera and Biodiversity (Florida Museum of Natural History, FLMNH), Coll. Kimmo Silvonen, Nico Elfferich; Kimmo Silvonen for the photos of the pierid, lycaenid and nymphalid pupae and of the riodinid larva, and Dan Janzen and Winnie Hallwachs for the photo of the hedylid larva (voucher code 03-SRNP-13034.1; for more information see http://janzen.sas.upenn.edu/caterpillars/database.lasso) all in the electronic supplementary material, figure S1. We also thank Marjatta Mikkonen, Ville Heimala and Pekka Malinen of MZH for their help. We are grateful to Dan Janzen and Winnie Hallwachs (through US NSF award nos DEB0072730 and DEB0515699), Felix Sperling, Chris Müller and Heikki Roininen for providing tissue samples crucial for this study. We thank Fabien Condamine, Thomas Simonsen and two anonymous referees for constructive comments on the manuscript. This study was funded by the Academy of Finland (projects 110906 to L.K. and 118369 to N.W.).
- Received July 8, 2011.
- Accepted August 23, 2011.
- This journal is © 2011 The Royal Society