Ehrlich and Raven proposed a model of coevolution where major host plant shifts of butterflies facilitate a burst of diversification driven by their arrival to a new adaptive zone. One prediction of this model is that reconstructions of historical diversification of butterflies should indicate an increase in diversification rate following major host shifts. Using reconstructed histories of 15 butterfly groups, I tested this prediction and found general agreement with Ehrlich and Raven's model. Butterfly lineages with an inferred major historical host shift showed evidence of diversification rate variation, with a significant acceleration following the host shift. Lineages without an inferred major host shift generally agreed with a constant-rate model of diversification. These results are consistent with the view that host plant associations have played a profound role in the evolutionary history of butterflies, and show that major shifts to chemically distinct plant groups leave a historical footprint that remains detectable today.
Ehrlich & Raven's (1964) seminal paper on phytochemically based host plant affiliations of butterflies spurred great interest in the evolutionary potential of reciprocal selection in generating the vast diversity of herbivorous insects. The coevolutionary model posited by Ehrlich and Raven is that the evolution of phytochemical novelty releases plants from herbivore pressure, thereby placing them in a ‘new adaptive zone’ that facilitates evolutionary radiation. Similarly, evolution of a trait in a herbivore that allows it to tolerate or ameliorate plant defences might lead to a burst of diversification. This model is often described as escape and radiate, and predicts that plant and herbivore lineages have undergone repeated bursts of diversification as a consequence of their antagonistic interactions (Thompson 1999).
Owing to the great abundance of herbivore species, and diverse guilds they occupy, overall selection pressure imposed by herbivores on plants might be considered diffuse (Rausher 1996, 2001). Thus, there might be no reason to expect that escape from any particular herbivore will facilitate plant diversification. Some have even questioned the premise that herbivores are important enough as selective agents to have had profound effects on plant diversification (Jermy 1984). There is, however, little debate over the importance of phytochemical variation as an important determinant of herbivore niche breadth (Janz & Nylin 1998; Cornell & Hawkins 2003; Futuyma & Agrawal 2009). One observation by Ehrlich and Raven was that related groups of butterflies tend to feed on related groups of plants (Brues 1920; Dethier 1954; Ackery 1988; Janz & Nylin 1998; Nishida 2002). Thus, the reciprocal selection pressure hypothesized by Ehrlich and Raven might be extremely asymmetrical, with escape from herbivores through chemical novelty having little impact on plant diversification, whereas shifts to chemically novel hosts create new adaptive zones that facilitate bursts in herbivore diversification rates (Farrell 1998; Wheat et al. 2007).
Evidence for bursts in diversification might be detected in the gene trees used for phylogenetic hypotheses (Nee et al. 1992, 1994a,b; Pybus & Harvey 2000; Nee 2001; Rabosky 2006b), and this expectation has been used to examine radiations in various groups (Lovette & Bermingham 1999; Barraclough & Vogler 2002; Harmon et al. 2003; Kozak et al. 2006; McKenna & Farrell 2006; Phillimore & Price 2008; Steeman et al. 2009; Winkler et al. 2009). In the case of butterflies, if host shifts facilitate rapid diversification, there is the expectation that branches on the tree following a host shift will be relatively short, and the accumulation of lineages across the tree will not follow a constant-rate model (Schluter 2000). In this study, I examine patterns of variation in diversification rate across multiple butterfly groups to assess whether major host shifts are associated with bursts of diversification. Here, I identify a major host shift as an inferred shift to a chemically distinct host plant group, usually inferred at the plant family level, to the exclusion of the major plant groups used by its sister lineage and ancestor. This definition excludes cases of host breadth expansion or contraction, and is consistent with Ehrlich and Raven's implied definition that host shifts place herbivores in a new adaptive zone.
2. Material and methods
(a) DNA data and phylogenetic reconstruction
In total, patterns of diversification were examined for 15 butterfly groups, representing members of Papilionidae, Nymphalidae and Pieridae. Of these datasets, eight include situations where a chemically based host shift is inferred to have occurred, and a shift in diversification rate might be expected. The remaining seven datasets do not include an inferred host shift, and therefore serve as a historical check, or negative control, in the analyses.
Sequence data for the various butterfly groups examined were obtained from GenBank and aligned by eye using MacClade v. 4.07 (Maddison & Maddison 2005). Regions of each data matrix with missing data were excluded from alignments, and taxa with incomplete sequences were removed. Only coding regions were used for analyses. Gene trees were estimated using a partitioned Bayesian phylogenetic reconstruction implemented in MrBayes v. 3.1.2 (Huelsenbeck & Ronquist 2001; Ronquist & Huelsenbeck 2003). For each model, two independent Markov Chain Monte Carlos were performed with one cold chain and three hot chains run for 2 × 106 generations with a burn-in of 1 × 105 generations. Partitions were defined by codon position for each gene. Nuclear and mitochondrial data were analysed separately (Degnan & Rosenberg 2006; Kubatko & Degnan 2007). Models of sequence evolution for each partition were selected using Modeltest 3.7 based on Akaike information criterion (AIC) (Posada & Crandall 1998). Ultrametric trees were estimated using penalized likelihood (Sanderson 2002), with the smoothing parameter estimated by cross-validation implemented in the program r8s v. 1.7 (Sanderson 2003). Where chronograms were provided in the primary literature, branching times were estimated from the original published figure using ImageJ (Abramoff et al. 2004).
(b) Host plant information
Data on host plant use were obtained from the primary literature, usually including the original study providing sequence data. Additional data were also obtained from field guides, books and the online database, HOSTS (http://www.nhm.ac.uk/research-curation/research/projects/hostplants/index.html).
(c) Diversification rate variation
Inferences of rate variation for each ultrametric tree were based on the vector of branching times. All analyses of rate variation were implemented in the statistical programming language R (R Development Core Team 2005) and included functions from the APE (Paradis et al. 2004) and LASER (Rabosky 2006a) packages. Two approaches were used to explore rate variation: (i) the constant-rate (CR) test, which examines the accumulation of lineages through time using the test statistic γ (Pybus & Harvey 2000), and (ii) fitting candidate rate-constant and rate-variable models to branching times (Rabosky 2006b). Both of these approaches test for a deviation in the vector of branching times against a null expectation of a constant rate.
(i) CR test
Under a pure-birth model of diversification, γ follows a standard normal distribution. Values of γ < −1.645 indicate a pattern of recent decrease in diversification rate, and might indicate whether lineage accumulation is consistent with a history of early, rapid diversification (Pybus & Harvey 2000). The CR test assumes complete taxon sampling. Incomplete taxon sampling can increase type I error because it tends to prune nodes at the tip of a tree, inflating what Nee et al. (1994a) called ‘the push of the past’, or the apparent high rate of early cladogenesis. To address type I error owing to incomplete taxon sampling, Pybus & Harvey proposed the Monte Carlo CR (MCCR) test, which calculates the critical value of γ at α = 0.05 for a given dataset by generating a distribution of pure-birth trees with incomplete taxon sampling. The 0.05 quantile of calculated γ statistics for 10 000 trees simulated in the program Phylogen (Rambaut 2002) was used to calculate the critical γ. The value for the possible number of taxa simulated for the analyses herein is based either on the opinion of authors of the study that was the source of sequence data, or by the number of taxa suggested by the websites http://www.funet.fi/pub/sci/bio/life/intro.html and http://tolweb.org/Lepidoptera. Discrepancies in total taxa among these sources did not affect statistical inferences presented here. Taxonomic inflation can present a problem for interpreting the pattern of cladogenic events near the present (Isaac et al. 2004). One approach to this problem is to truncate an arbitrary amount of time from the tree before analysis to reduce the influence of closely related taxa (e.g. Phillimore & Price 2008). A similar method was employed here, where values of γ were explored as nodes are serially truncated from the tip of the tree. The stability of γ values below the critical threshold as nodes are truncated provides an additional line of evidence regarding rate variation and reduces the influence of over-split taxa.
(ii) Model fitting
Likelihoods and parameter estimates of rate-constant models (pure birth and birth death) and rate-variable models (density-dependent exponential (DDX), density-dependent logistic (DDL) and two pure-birth rate (2-rate)) for each vector of branching times were calculated using the R package LASER (Rabosky 2006a). The AIC score for the best-fit rate-constant (AICRC) and rate-variable (AICRV) models was used to arrive at the best model given the data. Type I error rates were addressed by calculating a critical value for rejecting a constant-rate model by examining the distribution of differences between AICRC and AICRV of 10 000 pure-birth trees with taxon sampling simulated in Phylogen (Rabosky 2006b).
3. Results and discussion
Within the Papilionidae, Papilio and Parnassius exhibit a major inferred historical host shift from the Aristolochiaceae. Papilio is characterized by a shift primarily to Rutaceae and Apiaceae (Berenbaum 1983; Zakharov et al. 2004), and Parnassius has an inferred shift primarily to Crassulaceae and, shortly thereafter, to Papaveraceae and Fumariaceae (Omoto et al. 2004; Michel et al. 2008). Analyses of branching times largely indicated an acceleration followed by a decrease in diversification rate following these major host shifts (table 1 and figure 1a,c). A similar result was detected for Papilio and Parnassius using branching times of previously proposed chronograms (Zakharov et al. 2004; Michel et al. 2008). These results are consistent with the hypothesis that following a host shift, inferred to occur near the base of each tree, there is a burst followed by decreases in diversification rate.
Two other Papilionidae genera, Parides and Graphium, were examined (Makita et al. 2003; Silva-Brandao et al. 2005). Parides is a member of the tribe Troidini, which feed exclusively on Aristolochiaceae (Silva-Brandao et al. 2005; Silva-Brandao & Solferini 2007), and no major host shifts are inferred to occur throughout its history. Consistent with predictions of Ehrlich and Raven, Parides showed no evidence of rate variation (table 1 and figure 1b). Graphium is not involved with a major host shift, feeding primarily on Annonaceae and Magnoliaceae as does the closely related Eurytides, and similarly showed no evidence of rate variation (table 1).
Pieridae is characterized by at least one major host shift by the Pierinae from Fabaceae to Brassicaceae (Braby & Trueman 2006; Wheat et al. 2007). The Pieridae was examined using a genus-level phylogeny (Braby et al. 2006). The MCCR test indicated a decrease in diversification rate (table 1); however, this result should be interpreted with caution because, as it is a genus-level phylogeny, there is nowhere near-complete taxon sampling (75 terminal taxa in this tree compared with approx. 1100 possible species). Furthermore, the MCCR test assumes that incomplete sampling is random and, as this is a genus-level tree, this is not the case. However, γ remained below the critical value when more than 50 per cent of the nodes were truncated from the tree, where taxon sampling is near complete. Model fitting indicated rate variation across the tree, with an increase in diversification associated with the emergence of Pierinae and their inferred shift to Brassicaceae (table 1 and figure 1d). An analysis of Pierinae separately from the mostly Fabaceae-feeding Coliadinae showed that the detection of rate variation for Pieridae is largely owing to a burst of diversification following the shift to Brassicaceae. Models describing rate variation were consistently favoured over constant-rate models as nodes were truncated from the tree until a time prior to 40 Myr before present (figure 2). Wheat et al. (2007) described biochemical innovations that facilitated this shift by Pierinae. Using methodologies different from those employed here, they concluded that this major shift resulted in increased diversification rate. Analyses of their published chronogram gave similar results (table 1). The increase in diversification rates following the major shift to Brassicaceae is consistent with the escape and radiate hypothesis.
Nymphalidae is the most species-rich group of butterflies and has an extraordinarily rich host plant repertoire. Within the Nymphalinae clade, there exists a positive correlation between host plant diversity and species diversity among clades, and it has been argued that expansion of host plant repertoire (not shifts, per se) has been the important driver of diversification for this group, particularly for Nymphalini (Janz et al. 2006; Weingartner et al. 2006). Nylin & Wahlberg (2008) argue against the escape and radiate scenario for Nymphalini, as most species feed on ancestral urticalean rosids host, and are not characterized by major host shifts (Janz et al. 2001). Examinations of reconstructed gene trees support their argument (table 1).
Within Nymphalinae, an inferred colonization of Lamiales occurs at the origin of the Kallimoid clade (sensu Wahlberg 2006). An examination of branching times proposed by Wahlberg (2006) suggests that colonization of Lamiales was followed by a period of increased diversification rate, and although not densely sampled, the significance of the CR test was robust to 30 per cent node truncation. Nested within Kallimini (sensu Nylin & Wahlberg 2008), Junonia shares much of the same host plant repertoire as closely related groups and is not predicted to show a pattern of rate variation. Examination of Junonia indicated no evidence for rate variation (table 1). The tribe Melitaeini feed primarily on Asteridae, most of which contain iridoids, with the exception of Asteraceae and Acanthaceae. There have been host shifts within Melitaeini to non-iridoid containing plant families, notably in the Phyciodes and Chlosyne groups (sensu Wahlberg 2001). An examination of Melitaeini indicated acceleration in diversification rate following this shift (table 1 and figure 1e).
The Ithomiini (Danainae) are characterized by a presumed host shift early in their history from Apocynaceae to Solonaceae (Brower et al. 2006; Willmott & Freitas 2006). Analyses were largely consistent with the escape and radiate hypotheses, indicating a burst of diversification following the shift (table 1 and figure 1f). An examination of Ithomia, which itself is not involved with a shift, failed to detect rate variation (table 1).
Three tribes included in the subfamily Heliconiinae were examined: Argynnini, Acraeini, and Heliconiini. Argynnini feed mostly on Violaceae, though the most basal group, Euptoieta, also include Passifloraceae among their host plants (Simonsen 2006). No evidence for rate variation was detected in Argynnini (table 1). Excluding Euptoieta, thereby focusing on Violaceae-feeding groups, similarly failed to detect rate variation. This result is inconsistent with the escape and radiate hypothesis. Acraeini feed on various host plant families known to contain cyanogenic glycosides in the Old World, including Passifloraceae, although some species are highly polyphagous (Ackery 1988; Silva-Brandao et al. 2008). A shift to the non-cyanogenic glycoside containing Asteraceae is associated with their invasion of the New World (Silva-Brandao et al. 2008). Examination of diversification rates of Acraeini provided no evidence that arrival to the New World and the shift to Asteraceae had an impact on diversification rate (table 1). The Heliconii feed almost exclusively on Passifloraceae, usually on Passiflora (Brown 1981). There was evidence for acceleration in diversification rate at the origin of this group (table 1 and figure 1g). If the high degree of specialization on Passifloraceae is considered a host shift, this pattern is consistent with the escape and radiate hypothesis; however, other members of Heliconiinae are recorded to use Passifloraceae, and considering it a discrete shift in the spirit of Ehrlich and Raven's model might not be appropriate. Further examination of Heliconii showed that the ‘advanced genera’ (sensu Brown 1981), Heliconius and Eueides, characterized as slow-flying and unpalatable, are largely responsible for the pattern of decreased diversification rate (table 1 and figure 1h). Similarly, an examination of Heliconius showed evidence for a decrease in diversification rate (table 1). This apparent rate variation is not explicitly predicted by the escape and radiate hypothesis, as it is not associated with a host shift; however, this hypothesis does not exclude other possible mechanisms for rate variation. In particular, unpalatability might provide a new adaptive zone. Similarly, pollen feeding, which permits rapid larval growth and a long lifespan (Gilbert 1972; Beltran et al. 2007) in Heliconius, might too be an innovation that facilitates a burst in diversification rate.
One member of the Satyrinae, the genus Bicyclus, was examined. Bicyclus is mostly associated with Poaceae and is nested within other grass-feeding clades (Pena & Wahlberg 2008) and, thus, is not expected to exhibit rate variation owing to the process proposed by Ehrlich and Raven. Analyses failed to reject a constant-rate model of diversification based on the estimation of the gene tree constructed for this study. However, the chronogram proposed by Monteiro & Pierce (2001) indicates a reduction in diversification rates through time (table 1), raising the possibility that rate variation has occurred in Bicyclus.
(d) Summary of analyses
Ehrlich and Raven's model predicts that following a major host shift, a new adaptive zone is entered, resulting in a rapid accumulation of new, distinct lineages. The results of this study show that, overall, the pattern of lineage accumulation (i.e. diversification) is consistent with this prediction. Six of the eight datasets where historical host shifts are inferred to occur showed evidence of rate variation with an increase in diversification rate followed by a subsequent decrease. Across all datasets examined here, there was a higher likelihood of detecting a rate shift for clades associated with a host shift compared with clades where no host shift is inferred (Fisher's exact test, p = 0.03). Combined p-values show greater strength of evidence for rate variation in clades with a host shift (Stouffer's method: p ≪ 0.001 for MCCR and model fitting) compared with clades with no inferred host shift (p = 0.07 and p = 0.20 for MCCR and model fitting, respectively; Whitlock 2005). These results suggest major host shifts result in a burst of diversification that leaves a macroevolutionary footprint that can be detected in phylogenetic hypotheses.
The datasets examined here vary in completeness of taxon sampling and clade age, both factors that might affect the analyses. However, there was no systematic bias in completeness of taxon sampling between groups with or without inferred host shifts (t-test: t = 0.08, d.f. = 12, p = 0.93). Although many of the datasets lacked dense taxon sampling, there was no relationship between completeness of taxon sampling and significance of γ (r = −0.29, p = 0.28), or in confidence of rejecting constant-rate models (r = −0.28, p = 0.30). Similarly, using maximum pairwise distance (Kimura 2-parameter) as a surrogate for clade age, there was no age difference between groups with or without an inferred host shift (t-test: t = 0.30, d.f. = 12, p = 0.77), nor was the maximum distance correlated with the significance of γ (r = 0.06, p = 0.83) or confidence of rejecting constant-rate models (r = 0.05, p = 0.87).
(e) Host shifts and butterfly diversification
An important challenge in this study was to determine an unambiguous definition for a ‘major host shift’. The taxonomic grain at which host plants are considered is an arbitrary decision, whether it is at the level of shifts among congeners, or shifts across higher taxonomic ranks (Futuyma 1976; Futuyma et al. 1995; Farrell 1998; Janz et al. 2001; Wahlberg 2001; Silva-Brandao & Solferini 2007; Silva-Brandao et al. 2008; Forister et al. 2009). Most of the analyses presented here were consistent with the spirit of Ehrlich and Raven and others (Brues 1920; Dethier 1954), where plant family is the rank considered. Plant family level is often associated with characteristic defensive chemistry; however, many defensive compounds associated with particular plant families are not restricted to those families (e.g. furanocoumarins, glucosinolates and cardenolides). Furthermore, apparent host shifts or specialization can occur among lineages within a family. For example, there are some genera of Ithomini that are largely restricted to genera within the Solanaceae. The biological activity and suite of defensive compounds can even vary among conspecifics (Bowers & Stamp 1993; Zangerl & Berenbaum 1993). Where information was available, such as the shift in Melitaeini in regards to iridoids (Wahlberg 2001), that information was used to characterize an important host shift. Similarly, when chemical similarities among plant families were known, such as furanocoumarins in Apiaceae and Rutaceae host plants of some Papilio (Zakharov et al. 2004), it was considered a single major shift to a chemically distinct group. However, if plant chemistry is the important determinant of host use, the use of any taxonomic rank to determine ‘major host shifts’ provides an inherently coarse approximation (Janzen 1979). Finally, the salient aspects of a plant lineage's chemical phenotype existing at the time of colonization are not necessarily a good predictor of the chemotype today. If we accept that herbivore pressure imposes selection on plant chemistry (Berenbaum & Zangerl 1998), it is unlikely that the chemistry of today is identical to the phenotype present at the time of colonization.
Although the analyses presented here largely support Ehrlich and Raven's hypothesis that major host shifts result in a burst of diversification, other hypotheses have been proposed that link host plant use patterns and herbivore diversification. One hypothesis, particularly germane to butterflies, is the ‘oscillation hypothesis’ (Janz et al. 2006; Wahlberg 2006; Weingartner et al. 2006; Janz & Nylin 2008). This hypothesis was developed, in part, to explain the richness of the Nymphalidae and the staggering diversity of plant families used by many lineages. This hypothesis posits a repeated cycling between periods of host plant specialization and expansion, and predicts a sustained increase in speciation rate. It also predicts that newly incorporated host plants might be suboptimal, but used because of retained plasticity in host use, suggesting that complete host shifts to chemically distinct plant lineages might take an extended time (Nylin & Wahlberg 2008; Nylin & Janz 2009). Whereas the escape and radiate hypothesis emphasizes the importance of a new adaptive landscape as an engine of diversification, the oscillation hypothesis champions the importance of host breadth, facilitating the emergence of new lineages. A comparative analysis of Nymphalidae by Janz et al. (2006) showed a positive relationship between diversity of host use and species richness among clades, supporting the oscillation hypothesis.
It is possible that the role host plant shifts and diet breadth play in the diversification of butterflies varies across groups. The predicted pattern of accelerated diversification following a major host shift was supported in Papilionidae and Pieridae; however, this pattern was less pronounced in Nymphalidae. Interpreting patterns consistent with host shifts was particularly challenging for some groups of Nymphalidae. For example, Argynnini are largely characterized by a shift and subsequent specialization on Violaceae; however, the basal lineage, Euptoieta, uses both Violaceae and the inferred ancestral host, Passiflora (Ackery 1988; Simonsen 2006). Other species in the tribe use Violaceae in addition to other plant families (Ackery 1988), consistent with the oscillation hypothesis. The dual use of ancestral and derived host plants was less a problem in the Papilionidae and Pieridae, making inferences about host shifts clearer because they appear as discrete shifts. It is possible that the difficulty in determining whether major host shifts have occurred in some groups of Nymphalidae is owing to taxonomic rank for which data were available. For example, in reference to tribes examined within Heliconiinae, the level of tribe might not be appropriate to ascertain the importance of major host shifts if major shifts occurred prior to the origination of the currently defined groups. For the Heliconiinae, this major shift might have been the colonization of Passifloraceae. Subsequent specialization on other plant families, such as Violaceae, might have been facilitated by colonization of Passifloraceae, inviting further caution for using plant taxonomic rank to infer major host shifts (Nyman 2010).
(f) Methodological considerations
There are important methodological considerations that might affect conclusions based on analyses presented here. First, there are inherent limitations for reconstructing ancestral character states, and thus for inferring when a host shift occurred (Maddison 1994, 2006; Cunningham et al. 1998; Cunningham 1999; Paradis 2008). This might be less a problem here, as most datasets did not include host shifts within the analysed tree; rather, the available data are from studies examining relationships within groups characterized by major host shifts at the base of the tree. As taxonomic sampling becomes more extensive in phylogenetic studies of butterflies, for example, as has recently been done for the Nymphalidae (Wahlberg et al. 2009), there will be opportunities for more thorough comparisons of diversification rate variation (i.e. subtrees in a larger phylogenetic context and comparisons with sister clades). Additionally, recent studies have suggested that extinction events can lead to patterns consistent with a burst of diversification that are in fact owing to changes in extinction rates (Rabosky & Lovette 2008; Crisp & Cook 2009; Quental & Marshall 2009). This might be of less concern for this study because of its comparative nature, and there was an a priori expectation for the hypothesis being tested. Whether host shifts increase diversification rate, or increase the likelihood of lineage persistence through periods of increased extinction rate, might be of more semantic concern than evolutionary relevance in this particular context.
Another consideration is the role method of phylogeny reconstruction and chronogram development plays on inferred diversification patterns. Although nearly all the chronograms constructed for this study were in general agreement with chronograms proposed by previous studies, the agreement was not perfect. For example, the analyses of Ithomia failed to detect evidence of rate variation; however, Jiggins et al. (2006) concluded that Ithomia showed a decrease in diversification rate, though later analyses failed to detect this pattern (Elias et al. 2009). There is no reason to assume a systematic bias among trees examined here, as each step in the process followed the same criteria. Future comparative studies of diversification rate should similarly use consistent approaches across datasets, and caution should be exercised when generalizing results based on a single tree.
Ehrlich and Raven proposed a model of coevolution that explicitly addressed how evolutionary novelty enabling major shifts in host plant use might drive patterns of diversification. As phylogenetic data continue to accumulate, there will be more opportunities to examine variation in diversification rates as it relates to host plant associations. The evidence presented here suggests that, for butterflies, the interactions with their host plants, and occasional shifts among chemically distinct plant groups, have played an important role in their diversification history.
I thank the following people for discussion as this project developed and/or advice on improving this manuscript: C. C. Nice, B. M. Fitzpatrick, C. D. Hulsey, T. J. Near, Z. Gompert, M. L. Forister, B. C. O'Meara and A. M. Shapiro. This work was supported by the University of Tennessee and the US National Science Foundation (DEB-0614223).
- Received May 7, 2010.
- Accepted June 16, 2010.
- © 2010 The Royal Society