Synergistic effects of combining morphological and molecular data in resolving the phylogeny of butterflies and skippers

1. Introduction

Butterflies are arguably the best loved group of invertebrates and have been a source of inspiration for generations of natural historians and scientists. As a result, their generic- and specific-level classification is reasonably stable and the majority of taxa have been named (Ackery et al. 1999). Also unique for a large group of invertebrates is the immense biological knowledge amassed for many species, allowing butterflies to be used as a model group of organisms for wide ranging studies in ecology, evolution, population genetics and developmental biology (Boggs et al. 2003). Despite this wealth of information, the phylogenetic relationships and higher classification of the major groups of butterflies have remained contentious and competing hypotheses lack strong empirical support (Ehrlich 1958; Scott 1985; de Jong et al. 1996; Vane-Wright 2003).

The classification of butterflies and skippers has been based almost exclusively on the morphology of adult specimens for close to 250 years (Ackery et al. 1999), despite the rampant homoplasy in these morphologically variable insects. The usefulness of characters from immature stages has long been acknowledged (Müller 1886) but they have only recently been automatically incorporated in a phylogenetic context (Harvey 1991; Penz & Peggie 2003; Freitas & Brown 2004). Study of these features remains severely hampered by the lack of detailed descriptions or preserved specimens of larvae and pupae from which characters may be discerned.

Traditionally, the butterflies and skippers have been placed in two super-families and five families: Hesperiidae (skippers) are usually placed in their own super-family Hesperioidea while all other butterflies (Papilionidae, Pieridae, Lycaenidae and Nymphalidae) are placed in Papilionoidea. There is little agreement on the rank and monophyly of various groups and relationships among and within the families are largely unresolved (Vane-Wright 2003). For example, Riodinidae are sometimes separated from Lycaenidae, and Nymphalidae have been divided into as many as nine families. Papilionoidea and Hesperioidea have traditionally been considered sister taxa but recently the moth-like family Hedylidae (represented by the single genus Macrosoma) has been suggested to be more closely related to Papilionoidea than Hesperioidea (Scoble 1986, 1992), though this placement has been questioned (Weintraub & Miller 1987). Within Papilionoidea, Pieridae are either the sister group to the Papilionidae (Ehrlich 1958; Scott 1985) or to the (Nymphalidae+(Lycaenidae+Riodinidae)) (Kristensen 1976; de Jong et al. 1996; Weller et al. 1996; Ackery et al. 1999). Lycaenidae+Riodinidae are usually, but not always, considered to be the sister group of Nymphalidae (Scott 1985; Robbins 1988; de Jong et al. 1996).

DNA sequences have rarely been used to explicitly assess family level relationships of butterflies and skippers for exemplars of all families (Martin & Pashley 1992; Weller et al. 1996). These studies were based on short sequences of the nuclear 28S ribosomal DNA and the partial sequences of the mitochondrial gene ND1. The resulting trees were rooted with skippers, thus the monophyly of Papilionoidea was not tested. Molecular data have been utilized in a number of studies within lower taxa of butterflies, such as genera, tribes and subfamilies (Sperling 2003); however, there has been little coordination of effort among the studies cited in Sperling (2003). Thus, the many sequences available on public databases such as GenBank cannot be simply collated and analysed together because they do not represent homologous gene regions (Caterino et al. 2000). To alleviate this difficulty, we have agreed to sequence the same three genes for all studied taxa in our respective laboratories. The benefits of such cooperation are self-evident. The three genes we have chosen have been employed with great success in a variety of prior butterfly studies (Brower & Egan 1997; Brower 2000; Campbell et al. 2000; Caterino et al. 2001; Monteiro & Pierce 2001; Wahlberg et al. 2003; Megens et al. 2004), partly due to the availability of polymerase chain reaction primers that work well with all butterflies studied to date.

Recently, the utility of traditional morphological characters for inferring phylogenetic relationships has been heavily criticized and a purely molecular approach has been advocated (Hebert et al. 2003; Scotland et al. 2003). Such a stand is not supported by the fundamental characteristics of the two kinds of data (Miller et al. 1997; Baker & Gatesy 2002) nor by the actual use of morphological data in combination with molecular data (Jenner 2004). Here we investigate the relationships of higher taxa of butterflies and skippers using both traditional morphological data and new molecular data. We also demonstrate the synergistic impact of morphological data on the results of the simultaneous analysis of the combined dataset.

2. Material and Methods

Selection of taxa for sequencing was based on a published morphological study (de Jong et al. 1996) and covered all the major lineages in each butterfly and skipper family (see appendix A). In a few cases the genus coded for morphological characters was not available for molecular work. In such cases, a closely related genus was sequenced instead (see appendix A). Sequences of Cytochrome Oxidase subunit I (COI, 1531 bp), Elongation Factor-1α (EF-1α, 1225 bp) and wingless (403 bp) were generated in the laboratories of Brower, Pierce, Sperling and Wahlberg according to local protocols (see Brower 2000; Caterino et al. 2001; Monteiro & Pierce 2001; Wahlberg et al. 2003). Species identifications, voucher codes and deposition sites and GenBank accession numbers are given in appendix A. The morphological data of de Jong et al. (1996) were revised, with several characters being recoded, to yield a matrix of 99 characters (see electronic appendix for details). These were mainly scored from adult butterflies, comprising 39 wing venation, 19 leg, 14 head, 21 thoracic and two abdominal characters. In addition, four characters were taken from immature stages.

Heuristic parsimony searches were performed with NONA 2.0 (Goloboff 1998) via Winclada (Nixon 2002; 1000 random additions of taxa with 10 trees kept during each replication) for each dataset separately, for the combined molecular dataset and for the entire combined dataset. All characters were given equal weight. Out-groups were not specified a priori in any of the analyses and trees were rooted with a pyralid.

Bayesian phylogenetic analyses (Huelsenbeck et al. 2001) were performed using the programme MrBayes 3.1 (Ronquist & Huelsenbeck 2003). The GTR+G+I model of substitution (chosen using the programme MrModeltest; Nylander 2002) was fitted to each molecular partition and a rate variable model was fitted to the morphological data. Parameters were estimated for each of the genes and the morphology simultaneously (four partitions). Six chains (one cold and five heated) were run for 10 000 000 generations with trees sampled every 1000 generations. To check when stationarity was reached, likelihood values were graphically inspected and the first 1000 sampled trees were discarded as ‘burn-in’. Similar analyses were run without the morphological data.

Nodal support for the cladogram was assessed using Bremer support (BS) (Bremer 1988, 1994) values and bootstrap analysis. BS values and related indices (see below) were calculated with the aid of TreeRot (Sorensen 1999) and PAUP* (Swofford 2001). Bootstrap values were calculated with NONA 2.0 using 1000 pseudo-replicates with 10 random addition replicates per pseudo-replicate.

There are several indices related to BS that describe the interactions of datasets in a combined analysis framework at a node-by-node and character-by-character basis (Baker & DeSalle 1997; Gatesy et al. 1999). Each data partition contributes additively to the total BS in the combined analysis framework, giving values known as partitioned Bremer support (PBS; Baker & DeSalle 1997; Gatesy et al. 1999). A positive PBS value indicates that a given data partition supports a given node while a negative value indicates conflict from a given data partition.

Patterns of homoplasy differ in different datasets (Barrett et al. 1991; Chippindale & Wiens 1994; Olmstead & Sweere 1994) and analysing the different datasets together can bring forth the underlying phylogenetic signal common to all datasets (Baker & DeSalle 1997; Baker et al. 1998). BS values can allow us to compare the effects of combined and separate analyses through an index known as hidden Bremer support (HBS; Gatesy et al. 1999). HBS is the difference between the BS of a given node in the combined analysis and the sum of the BS of each partition for the same node in the most parsimonious tree for that partition. This index allows us to evaluate whether there is increased character support for a clade in combined analysis of multiple datasets compared with the sum of support for that clade in the separate analyses of the different partitions. Positive values indicate increased character support in the combined analysis whereas negative values indicate increased character conflict in the combined analysis.

3. Results

Our analysis is based on sequence data from two nuclear gene regions (1225 bp of EF-1α and 403 bp of Wingless) and one mitochondrial gene region (1531 bp of COI) for a total of 3159 base pairs and 99 morphological characters revised from a published study (de Jong et al. 1996) comprising a matrix of 57 taxa and 3258 characters (of which 1290 are parsimony informative; table 1). Analyses of the data partitions on their own results in many unconventional relationships (see electronic appendix). Our results for the morphology partition differ from those of de Jong et al. (1996), mainly due to our sparse out-group sampling (de Jong et al. sampled 15 out-group taxa).

Analyses based on the combined molecular data also fail to recover some of the traditionally recognized higher taxa as monophyletic groups (figure 1). The parsimony analysis finds many unconventional relationships (figure 1a). The low support for most nodes is probably due to large amounts of homoplasy and sparse sampling (there are about 20 000 recognized species of butterflies; Robbins 1982). Attempting to take into account these problems by using a parameter rich model (each gene region having unique parameters; see electronic appendix) and Bayesian methodology to estimate the phylogenetic relationships results in a more conventional phylogeny, although there are still several highly unconventional clades (figure 1b). For instance, the Bayesian analysis places Papilionidae as sister to the rest of the butterflies, skippers and hedylids with a high posterior probability (100%) and also places the nymphalid satyrine clade (sensu Wahlberg et al. 2003) as sister to the Riodinidae+Lycaenidae clade with a posterior probability of 100%. Such relationships have never been suggested previously in the literature.

• Download figure
• Open in new tab
• Download PowerPoint

In contrast, the analyses of the combined molecular and morphology dataset provide strong support for the monophyly of all traditionally recognized higher taxa, except Nymphalidae which has moderate BS (6), high posterior probability (100%) but no bootstrap support (less than 50%; figure 2). Both analyses place Riodinidae as sister to Lycaenidae with a monophyletic Nymphalidae sister to the former two. Pieridae is placed as sister to (Nymphalidae+(Lycaenidae+Riodinidae)). Hedylidae is placed as sister to Hesperioidea+Papilionoidea. Support values for most of the nodes describing higher taxa are strong (BS ≥ 9, bootstrap ≥ 80%, posterior probability ≥ 95%; figure 2, table 2). The exceptions are the nodes describing Hesperioidea+Papilionoidea, Pieridae as sister to (Nymphalidae+(Lycaenidae+Riodinidae)) and monophyly of Nymphalidae which have BS values of 6 or less, although the same nodes have posterior probabilities higher than 95%.

• Download figure
• Open in new tab
• Download PowerPoint

PBS values show that in the combined analysis all the data partitions contribute positively to the support of the relationships in figure 2 (table 2). Since the magnitude of BS is contingent on the number of parsimony informative characters, dividing the total PBS for a partition by the number of parsimony informative characters of that partition gives us an index of the relative informativeness of each partition. The index clearly shows that the morphological (0.95) and EF-1α (0.56) partitions contribute most of the support for the relationships in figure 2a. These two partitions support almost all the major nodes and there is weak conflict derived from them at three nodes (major nodes shown in bold in table 2). The COI (0.29) and wingless (0.30) partitions show a higher level of conflict at many of the nodes yet both show strong support for the major nodes of interest.

Despite the strong support for most clades, we find that each data partition shows a great deal of homoplasy when mapped on to the tree found in combined parsimony analysis (table 1). However, combining the different data partitions reveals strong hidden support for almost all of the major nodes (table 2), indicating that combined analysis enhances the inherent phylogenetic signal of the data compared with analysing the partitions separately. Partitioning the HBS values indicates that all four data partitions show substantial hidden signals that are revealed when analysed in combination with the other datasets.

4. Discussion

In the discussion of their analysis of butterfly relationships based on morphology (the dataset employed here), de Jong et al. (1996) expressed a degree of dissatisfaction with the power of morphological features to resolve many clades. They questioned whether additional morphological data would improve the results and called for the pursuit of alternative sources of evidence. However, they qualified these controversial statements by arguing that the most productive avenue would be to combine those alternative data with morphological characters, as we have done here. Many current phylogenetic studies do not take morphological data into account for a variety of reasons (see Baker & Gatesy 2002; Scotland et al. 2003; Jenner 2004). As we have shown, ignoring the morphological data in this case resulted in several spurious clades, regardless of the method of analysis used (figure 1), including a paraphyletic Nymphalidae with regard to Lycaenidae+Riodinidae, a sister relationship between Hedylidae and Papilionidae and a paraphyletic Papilionoidea with regard to Hesperioidea and Hedyloidea.

So we have an interesting problem: morphological data on their own provide weak support for, or fail to resolve, many nodes and molecular data on their own yield trees that contain implausible clades, some of which are strongly supported. A common claim from molecular systematists is that a large number of independent characters are needed to be able to estimate phylogenetic relationships robustly and reliably (Rokas et al. 2003). Yet for many groups of organisms, it is still not feasible to generate large amounts of DNA sequence data from a variety of gene regions. In Lepidoptera, for example, it has been a challenge to discover primers that are universal enough to amplify gene regions that are of broad phylogenetic utility in resolving all hierarchic taxonomic levels, from species to super-families (Friedlander et al. 1994). We feel that the propensity of limited molecular datasets to imply wrong, yet well supported, topologies is a problem that is not fully appreciated. Researchers now generally try to avoid the mistakes of early molecular systematics when dramatic rearrangements of phylogeny were often proposed on the basis of small and sparsely sampled datasets. However, there is still a tendency to reach very general conclusions from very specific and limited data. In particular, we feel that it should be more appreciated that adding more data is superior to adding ever more thorough analyses of existing data when the phylogeny is uncertain.

Our solution to this problem has been to combine our sequences with available morphological data. Morphological features are often uniquely derived complex structures that unequivocally unite major taxa, such as the osmeterial glands of papilionid larvae and the tricarinate antennae of nymphalid adults, and we see them as an intrinsic component of a robust treatment of this important group. Even though all characters are given equal weight and the morphological characters are vastly outnumbered by the molecular data, the low intrinsic homoplasy of the morphological features allows them to establish a structural framework in the combined analysis to which the morphological data contributes synergistically. Although the morphological data on their own were unable to give unequivocal answers in our analyses, the combined analyses (regardless of method used) have given us a robust phylogenetic hypothesis for the butterflies and skippers and the whole is greater than the sum of its parts, as is shown by the positive HBS values.

Our analyses of the combined dataset provide evidence to settle a long-standing question on the position and rank of Riodinidae. Our data give strong support for a Lycaenidae+Riodinidae sister group relationship, rather than embedding Riodinidae within Lycaenidae (de Jong et al. 1996) or placing them as sister to Nymphalidae (Robbins 1989). This result is concordant with an earlier molecular study that included a larger sample of taxa for the three groups but only used evidence from wingless (Campbell et al. 2000). The high BS and full congruence among the different datasets (table 2) indicates that only a great deal of contradictory new evidence could overturn this sister relationship. We believe that the current analysis provides a long-awaited empirical rationale for ranking Lycaenidae and Riodinidae as separate families.

The monophyly of Hesperiidae, Papilionidae and Pieridae is also strongly supported (figure 2) and there is little or no conflict between the data partitions at these nodes (table 2). These three families are generally considered to form monophyletic groups, though occasionally some small groups have been split off into their own families (e.g. Megathymidae from Hesperiidae and Baroniidae from Papilionidae). The case of Nymphalidae, however, is more complicated. Morphology, COI and EF-1α sequence data support the monophyly of the family while wingless conflicts moderately. The hidden BS of the partitions are positive, with the molecular partitions exhibiting larger values than the morphological partition. This suggests that there is a large amount of homoplasy (noise) in the molecular partitions for Nymphalidae that is overcome by the combined analysis of all data. The basal branch leading to Nymphalidae appears to be much shorter than the branches leading to the other families (figure 2b). One interpretation of such a pattern might be that the major lineages in Nymphalidae diverged very rapidly from one another. This may explain the difficulty in determining homologous states of morphological characters among nymphalids (Freitas & Brown 2004) and the difficulty in finding uniquely derived features to define the family (de Jong et al. 1996).

Relationships of taxa within families generally agree with recent family level studies (Brower 2000; Campbell et al. 2000; Caterino et al. 2001; Hall & Harvey 2002; Wahlberg et al. 2003; Freitas & Brown 2004), though some incongruent relationships can be noted. These incongruities are probably the result of insufficient taxon sampling in this study. The studies cited above included many more taxa per group and were generally able to resolve the relationships of major lineages within the families with higher support.

This is the first study to provide strongly supported answers to many of the questions about butterfly relationships so eloquently raised by Vane-Wright and co-workers (de Jong et al. 1996; Vane-Wright 2003). Our data show that the six currently recognized skipper and butterfly families (Hesperiidae, Papilionidae, Pieridae, Lycaenidae, Riodinidae and Nymphalidae) are monophyletic entities, that Lycaenidae and Riodinidae are sister taxa, that Nymphalidae are sister to Lycaenidae+Riodinidae, that Pieridae are sister to (Nymphalidae+(Lycaenidae+Riodinidae)) and lend support to the hypothesis that Hedylidae are more closely related to the butterflies and skippers than the geometrids within which they were once classified. Many of these taxa have been intuitively grouped based on striking morphological novelties since the time of Linnaeus but this study is the first to provide robust empirical support for relationships based on rigorous analysis of characters, both morphological and molecular. These results once again emphasize the importance of combining evidence from different sources as a means to dilute the bias of homoplasy within any individual data partition (Farris 1983; Miller et al. 1997).

APPENDIX A

See table 3.