Pieridae use three major host plant groups: the Fabales (Legumes), the Brassicales (glucosinolate-containing plants exemplified by the cabbages and Arabidopsis), and mistletoes. Phylogenetic reconstruction of almost 90% of the Pieridae genera (74 recognized genera plus six subgenera, based on 1,066 bp of the EF-1α gene) was recently completed (13). These results indicate that Fabales feeding is the ancestral state of Pieridae (Fig. 1). The Fabales feeders are the Dismorphiinae and nearly all Coliadinae, whereas the sister to the Coliadinae, the Pierinae, primarily feed on Brassicales (Fig. 1) (12). Within Pierinae, there are two subsequent derived shifts away from glucosinolate feeding onto mistletoes and other species. Thus, the Pierinae represent a single origin of glucosinolate feeding (Fig. 1).

Fig. 1. Chronogram of the Pieridae, showing host plant relationships and the NSP activity of representative species. Circles represent presence (filled) or absence (empty) of NSP activity in larvae from relevant genera. Vertical labels indicate higher taxa and (more ...)

The evolutionary appearance of the plant order Brassicales (Eurosid II, Dicotyledons) presented a radical new chemical challenge for insect herbivores, known as the glucosinolate-myrosinase system (14). All Brassicales have this system, which is one of the best and most widely studied chemical plant defenses (15–17). Its effectiveness as an anti-herbivore defense becomes apparent upon tissue damage, such as insect feeding. Tissue damage brings the formerly compartmentalized myrosinase enzyme into contact with nontoxic glucosinolates, which it hydrolyzes into breakdown products such as isothiocyanates (18–20). Whereas Homo sapiens may find these breakdown products enjoyable condiments (e.g., mustard, wassabi), they are well known to be highly toxic to many insect herbivores (15, 21–23).

We have identified two independent lepidopteran detoxification mechanisms for the glucosinolate-myrosinase defense system at biochemical and molecular levels by means of functional genomics approaches, beginning with glucosinolate sulfatase (GSS) in the diamondback moth Plutella xylostella (Plutellidae) (24). GSS desulfates glucosinolates, producing metabolites that no longer act as substrates for myrosinases. The second, called nitrile-specifier protein (NSP), has recently been identified for the pierid butterfly Pieris rapae (20). NSP, expressed solely in the larval midgut, promotes the formation of nitrile breakdown products instead of toxic isothiocyanates upon myrosinase-catalyzed glucosinolate hydrolysis (20). The GSS and NSP detoxification mechanisms are distinctly different from each other, as well as the other identified host plant detoxification mechanisms of the Lepidoptera: the cytochrome P450 monooxygenase gene superfamily in Papilio butterflies used against furanocoumarins (25) and the flavin-dependent monooxygenase system of the arctiid moth Tyria jacobaeae used against pyrrolizidine alkaloids (26).

With these issues in focus, we document a key innovation and its role in butterfly diversification by (i) assessing the evolutionary origin of a host plant detoxification mechanism involved in an ancient host plant shift, (ii) generating a robust estimate of when this detoxification mechanism evolved relative to the age of the colonized host plant family, and (iii) testing for an effect of the host shift on butterfly diversification rate.

We begin by asking a mechanistic question, focusing on whether all Pierid glucosinolate feeders use NSP, not NSP, or the other identified glucosinolate detoxification mechanism of the Lepidoptera, GSS. Midgut assays of larvae from 13 Pieridae species from North America, Europe, Africa, and Australia, representing each of the major subfamilies or groups, reveals complete concordance between NSP activity and glucosinolate feeding (Fig. 1; and see supporting information (SI) Fig. 3 and SI Table 1). Consistent with initial studies of these mechanisms on Pieris rapae, none of the butterfly species assayed displays any detectable GSS activity, and P. xylostella and other moths have no NSP activity. Within the Pieridae, glucosinolate feeding exactly matches the presence of NSP activity in larvae, suggesting a single evolutionary origin with a subsequent loss after secondary host plant shifts to Santalales (Fig. 1, SI Fig. 3, and SI Table 1). These results support NSP as a key innovation within the Pieridae enabling glucosinolate detoxification (27).

Fossil and molecular data agree that the Brassicales appeared by 90 to 85 Mya, which is much earlier than the parallel evolution of glucosinolates in the Euphorbiaceae (58 Mya) (28). Pierinae fossils with modern relatives in the Brassicales feeding clade appear ≈34 Mya, ≈50 Myr after the first fossil Brassicales (13, 28). Fossils inherently provide only minimum temporal estimates. Thus, the fossil record is compatible with either a long delay of insect adaptation after the Brassicales radiation, a contemporaneous radiation of Brassicales and their Pierid herbivores, or anything in between. Among the major insect orders, Lepidoptera have the poorest fossil record, probably as a result of their wing structure (8, 75). In stark comparison, the fossil record of beetles is very rich, and beetles provide some of the best evidence for coevolution (29), although this evidence has recently been questioned by a recalibrated phylogeny of the leaf beetles (30). In this article, we directly address the timing of the appearance of the glucosinolate-feeding Pierinae, using several independent molecular datasets and various calibration methods to generate a robust estimate of when Pieridae evolved relative to their Brassicales host plants.

Our temporal analysis used two standard molecular phylogenetic genes, cytochrome oxidase I (COI) and elongation factor 1α (EF-1α), with two and three different temporal estimations, respectively (Fig. 2). First, we compared COI sequence from Colias eurytheme (Coliadinae) and Pieris napi (Pierinae) and estimated the time of their divergence at 89.12 and 96.32 Mya. These two estimates are respectively derived (i) from an Insecta-wide average substitution rate of 0.022% per million years for COI 2nd position data calibrated on ancient, nonlepidopteran fossils (31), and (ii) from COI 1st and 2nd position data calibrated on butterfly fossils from within the Nymphalinae subclade of Nymphalidae (32). Second, we compared the EF-1α sequence distance between our Coliadinae and Pierinae representatives to the mean EF-1α sequence distance between temporally calibrated representatives in the Nymphalidae and Papilionidae butterfly families (32, 33). Nymphalidae estimates were calibrated as stated above, and the Papilionidae samples were calibrated by using biogeographic vicariance events inferred during the breakup of Gondwana. Molecular evolution rates of the nuclear EF-1α gene are likely to be similar among these families, as Nymphalidae and Papilionidae are Pieridae's nearest family clades (34). Together these EF-1α distance comparisons suggest that Coliadinae diverged from Pierinae 62 and >82 Mya, respectively (Fig. 2).

Fig. 2. Current and expected diversity, and divergence estimates for the Coliadinae and Pierinae subfamilies of Pieridae. (a) The x axis gives the age of clades (Mya) whereas the y axis gives the species numbers (Log scale) for each clade. The upper and lower (more ...)

Pieridae-specific temporal reconstruction used a Bayesian relaxed molecular clock method on EF-1α data from across Pieridae. Divergence estimates used the four known pierid fossils as minimum node dates (35, 36)(SI Figs. 4 and 5). Temporal reconstruction can be affected by assumptions used in analyses, such as tree topology, node calibration, clock-like evolution, and other prior estimates (37). In our analysis, we explored the potential effects of alternative topologies, fossil node placements, clock assumptions, and other prior estimates on temporal inferences. These alternatives had little effects on results except to increase age estimates (SI Table 2). The most temporally conservative model predicts the appearance of Pieridae with a 95% confidence interval of 166 to 79 Mya, which agrees with related clock-based estimates (13). In sum, our five temporal estimates are consistent and provide a robust assessment. The glucosinolate-feeding Pierinae diverged from Coliadinae in close evolutionary time to the appearance of the Brassicales, most likely within 10 million years (Figs. 1 and 2).

Ehrlich and Raven hypothesized, but were unable to test, that the “main diversification of [Pieridae] occurred after it became associated with the [Brassicales]” (8). Others have posited a relationship between host plant use and species abundance, but this relationship has not been tested (11, 12). We calculate and compare the standing diversity of the sister clades Coliadinae and Pierinae (38), which are a priori expected to have similar levels of diversification because they are sister clades (39). There are both more genera and species of Pierinae than Coliadinae (57 vs. 18 genera, 830 vs. 220 species). To assess diversification rate and extinction effects giving rise to these differences in extant taxa, we use the diversification rate of Dismorphiinae plus Pseudopontiinae (sister clade to Coliadinae plus Pierinae) as the background rate for generating expected species diversity confidence intervals through time. Coliadinae species diversity is at or below the expected numbers of species. However, there are significantly more species of Pierinae than expected, even if non-glucosinolate feeders are excluded, suggesting that the increased diversification of the Pierinae resulted from the host plant shift onto Brassicales (Fig. 2).

The evolution of the NSP glucosinolate detoxification gene was a novel, key innovation facilitating the host plant shift of ancestral Pierinae from Fabales to Brassicales feeding with macroevolutionary consequences. Our analysis of major clade representatives across Pieridae demonstrates that only glucosinolate-feeding Pierinae show NSP activity, indicating a single mechanistic basis of glucosinolate feeding originating within the Pierinae. Multiple independent temporal reconstructions provide a robust estimate placing the Pierinae host shift soon after the appearance of the Brassicales plant order. This host shift also resulted in significantly increased speciation compared with related clades. Together, these results provide strong support for the central tenet in Ehrlich and Raven's coevolution theory, that key biochemical innovations foster increased speciation (8, 40).

Ancestral Pieridae butterflies did not provide the initial herbivory pressure for Brassicales chemical defense evolution, nor are they solely exerting sufficient pressure for the continued evolution of the glucosinolate-myrosinase system of the Brassicales, as lepidopteran host plant coevolution is expected to be diffuse (8, 40, 41). Yet, comparative functional genomic analysis of the glucosinolate-myrosinase system in the family Brassicaceae reveals ongoing evolution in response to herbivory pressure (42–44). Multiple enzymes, with allelic variation, are involved in the biosynthesis of diverse glucosinolate structures and relative concentrations (45–47). This resulting intra- and interspecific glucosinolate profile variation results in herbivory level variation (48–51). In addition, signatures of positive selection at the molecular level have been identified at these genes, which are undergoing gene duplication and neofunctionalization, documenting ongoing glucosinolate metabolic evolution (52). Thus, the initial innovation that likely facilitated Brassicales escape and radiation from herbivore pressure is continuing to evolve in derived lineages with herbivory consequences.

Studies of plant–insect coevolution have evolved significantly since the days of Ehrlich and Raven (e.g., refs. 53–64). Butterfly support for the parallel coevolutionary process, i.e., cospeciation, has proven difficult to find, as cospeciation-driven diversification seems to be rare among insects and butterflies in particular, with Ehrlich and Raven's escape and radiate coevolutionary mechanism likely more common (8, 10, 40, 53, 65, 66). The only strong case for butterfly coevolution is found in the genus Papilio (Papilionidae), where ≈75% of this genus feeds on plants containing furanocoumarin-based chemical defenses, synthesized by plants using members of the cytochrome P450 monooxygenases gene superfamily (11, 25, 67). Papilio butterfly species diversity increases with host plant furanocoumarin diversity (11). Interestingly, furanocoumarins both induce and are detoxified by Papilio species using other members of the cytochrome P450 superfamily (68). Cytochrome P450 participation on both sides of the plant–insect interaction highlights the diverse functional roles members of this gene superfamily play in environmental response (68). In fact, the generalist Noctuid moth Helicoverpa, which is >100 Myr divergent from Papilio sp., independently uses other P450s to detoxify furanocoumarins, highlighting the repeatability of this gene superfamily's role in detoxification. Cytochrome P450s are also part of the glucosinolate biosynthetic pathway in Brassicales plants. Interestingly, in Lepidoptera specialized on glucosinolates, only the novel molecular mechanisms discussed in this article (GSS and NSP) have been identified. Given the apparently rapid host shift of the Pierinae onto Brassicales, we are investigating the evolutionary history of the NSP gene to assess its ancestral, pre-glucosinolate detoxification function.

Over 40 years ago, the authors who introduced us to coevolution using butterflies claimed that their “predictions cannot be tested” (8). The work presented here documents a key coevolutionary innovation in the Pieridae completely independent of other known lepidopteran detoxification mechanisms. Moreover, this work indicates a much more dynamic interaction between butterfly herbivores and their host plants than currently envisioned. We present our work here in the hope that evolutionary and ecological functional genomics will encourage increased mechanistic study of plant–insect interactions.

NSP and GSS Assay. Last instar larvae actively feeding on preferred host plants (SI Table 1) were used to test for NSP activity using methods described in ref. 20. Briefly, benzylglucosinolate and myrosinase were incubated with larval midgut extracts. Dichloromethane extracts of these assay mixtures were then analyzed by GC-MS with internal standard, with total ion current traces identifying NSP activity (see SI Fig. 3 for representative traces and more details). Methods for detecting GSS activity in samples were as reported in ref. 24, although a much greater amount of protein (up to 100 μg) was used in the assays in this study to detect potentially low enzyme activity levels.

Phylogenetic Reconstruction. Bayesian phylogenetic reconstruction, using a GTR + I + G nucleotide substitution model for a total of 1 million generations as implemented in Mr. Bayes v.3.1.2 (69, 70), was performed by using two different outgroups and different sets of the EF-1α dataset of Braby et al. (13). Resulting tree 1 used Papilio machaon (Papilionidae) as an outgroup, without representatives Pseudopontia or Dismorphiniae (SI Fig. 4). Tree 2 used all available EF-1α Pieridae sequences with Pseudopontia paradoxa as an outgroup (SI Fig. 5). Trees 1 and 2 place Dismorphiniae sister to Coliadinae plus Pierinae, concordant with previous studies (13, 71). Tree topologies differed in their terminal branching patterns, whereas subfamily relationships were stable (see SI Materials and Methods for analysis details).

Ancestral Node Dating.

COI. Recent analysis of Nymphalidae COI data, using fossil data for temporal calibration, estimated origins of the subgroup Nymphalinae at or older than 65 Mya (32). Using representative species Melitaea and Vanessa, from distant Nymphalinae clades, for 1st and 2nd position rate change per million year provides an estimate of 0.0004372 (26.8 changes per 943 bp; 26.8/943 = 0.02841; 0.02841/65 = 0.0004372 substitutions per million years). There are 43 1st plus 2nd position changes between C. eurytheme and P. napi, giving a temporal estimate of the Coliadinae and Pierinae divergence at 96.32 Mya (43/1021 = 0.04211; 0.04211/.0004372 = 96.32).

EF1α. The estimated EF-1α mean genetic distance and standard deviation between Coliadinae (n = 15) and Pierinae (n = 60) taxa is 12.928% (SE = 0.0095), whereas that between Papilio species (n = 55) and Pachliopta neptunus (Papilionidae) is nearly 3% lower at 10.029% (SE = 0.1115). Genetic distance was calculated by using a Tamura 3-Parameter distance model as implemented in Mega 3.1, with standard error calculated with 500 bootstrap replicates (72). Rate variation among sites was approximated by using a gama shape parameter = 1.4 [as per the estimated value for the dataset (33)]. The two Papilionidae groups are calculated to have diverged 89.1 to 82.5 Mya, suggesting a comparable or even older divergence between Coliadinae and Pierinae. Similar EF-1α calibrated estimates are available from Nymphalidae. By using representatives from Vanessa and Melitaea for a 65 million year split, there is a mean distance of 7.64% (SE = 0.0078). This distance is 60% of the divergence found between Coliadinae and Pierinae representatives, suggesting a Coliadinae Pierinae divergence much greater than 65 Mya.

EF1α using a relaxed clock method. Temporal estimation using a Bayesian relaxed molecular clock method was implemented by using MULTIDIVTIME (35, 73) on the two Bayesian analysis EF-1α-generated trees (SI Figs. 4 and 5). Prior estimates (SI Table 2) explored model response across many different prior values following a similar approach by Wiegmann et al. (74). To perform the most conservative analysis possible, all analyses used a prior maximum tree age set to the oldest known fossil Pierinae at 34 Mya. For all combinations of prior values, the Markov Chain was then sampled every 100 cycles for a total of 10,000 samples. Posterior distributions were approximated based on these 10,000 samples. Fossil node placement initially followed Braby et al. (13), but here fossils are assigned as a minimum date estimate and the effect of more conservative node placement was explored by using the next most ancestral node (SI Figs. 4 and 5 and SI Table 2). For each run, across two trees, the prior conditions, resulting estimates, and confidence limits are listed in SI Table 2.

Absolute Diversification Rate Calculations. There are both more genera and species of Pierinae than Coliadinae (57 vs. 18 genera; 830 vs. 220 species). These numbers are significantly different from equal number expectations (χ² = 10.8, P = 0.001; χ² = 193.5, P < 0.001). However, this simple calculation does not assess diversification rate and extinction effects giving rise to these differences in extant taxa. To address these effects, we perform a statistical analysis of a stochastic birth-and-death processes, employing a method-of-moments estimator following Magallon and Sanderson (38), which takes the age of the clades in question, together with their extant species diversities, and allows for them to be compared with their expected levels of diversification given different extinction rate assumptions (see SI Materials and Methods for more details). Sister clades Coliadinae and Pieridinae are by definition of equal age, and, for an unbiased, phylogenetically relevant estimator of their expected species diversity, we chose the diversification rate of their sister clade, Dismorphiinae plus Pseudopontiinae. By using this rate, 95% confidence limits were calculated for expected numbers of species across different clade ages with either no extinction or a high (90%) extinction rate (see SI Materials and Methods for more details; SI Table 3). Groups with a greater than expected number of species are considered to be “excessively species rich” (38). The 95% confidence interval for high extinction is generally greater than under no extinction due to having a higher variance; a high number of speciation events are needed to offset the high number of extinctions to maintain the same net diversification rate as when there is no extinction.

We thank K. Bargum, M. Clauss, J. Marden, S. Nylin, R. Schilder, N. Wahlberg, and W. Watt for comments on various stages of this manuscript; R. Oyama and J. Thorne for analysis assistance; and A. Bonkewitz and D. Heckel for specimen collection help. M. Borsch, S. Coombes, and M. Rowlings kindly permitted our use of their butterfly photographs (A. crataegi, A. cardamines, and C. hyale and C. evagore, respectively). This work was funded by the Max Planck Society and National Science Foundation Grant IBN-0412651 (to C.W.W.).