Communicated by Ebbe S. Nielsen, Commonwealth Scientific and Industrial Research Organization, Canberra, Australia (received for review October 30, 1998)
The obligate mutualism between yuccas and yucca moths is a major model system for the study of coevolving species interactions. Exploration of the processes that have generated current diversity and associations within this mutualism requires robust phylogenies and timelines for both moths and yuccas. Here we establish a molecular clock for the moths based on mtDNA and use it to estimate the time of major life history events within the yucca moths. Colonization of yuccas had occurred by 41.5 ± 9.8 million years ago (Mya), with rapid life history diversification and the emergence of pollinators within 0–6 My after yucca colonization. A subsequent burst of diversification 3.2 ± 1.8 Mya coincided with evolution of arid habitats in western North America. Derived nonpollinating cheater yucca moths evolved 1.26 ± 0.96 Mya. The estimated age of the moths far predates the host fossil record, but is consistent with suggested host age based on paleobotanical, climatological, biogeographical, and geological data, and a tentative estimation from an rbcL-based molecular clock for yuccas. The moth data are used to establish three alternative scenarios of how the moths and plants have coevolved. They yield specific predictions that can be tested once a robust plant phylogeny becomes available.
Obligate pollination mutualisms such as the yucca-yucca moth and fig-fig wasp associations provide some of the classically cited examples of coevolution (1, 2). In these interactions, the adult insects serve as the exclusive pollinators of their hosts. Their larvae subsequently feed on host seeds, but because many seeds are left intact the interaction carries a net positive effect for the plants. Both plants and insects show obvious coadapted traits, such as specific pollen collection and deposition behaviors in the insects, and structural adaptations that mediate pollinator specificity in the flowers.
Given the strong associations between these organisms, they serve as excellent model systems for exploring many aspects of evolutionary biology, such as the origins of mutualism (3–7), the evolution of virulence (8), evolution of mating systems (9), stability and reversal of obligate mutualism (5, 7, 10, 11), and consequences of specialization on population structure in species (12). Because of this utility for many branches of evolutionary biology, it would be highly desirable to develop a strong phylogenetic framework for the organisms and to establish a timeline for their diversification. This framework would allow for macroevolutionary analyses of the insect-host association, assessing for example the importance of codiversification (13), and the role of intrinsic and extrinsic factors in driving diversification (14).
Here we use molecular data in conjunction with biogeographic and fossil data to develop a phylogeny for the yucca moth family, Prodoxidae, to estimate minimum ages of divergence and determine patterns of diversification within the family. This work builds on an earlier study (15) based on limited molecular data that suggested variable rates of diversification among prodoxid genera. The present results support the following conclusions: (i) an explosive radiation of the yucca moths and evolution of pollination behavior had occurred no later than 40 million years ago (Mya), (ii) a second burst of yucca moth diversification is coincident with Pliocene desertification, (iii) derived nonpollinating cheater yucca moths have co-occurred with pollinators for at least 1.26 ± 0.96 My, and (iv) the radiation of yuccas occurred much earlier than suggested by the sparse fossil record. The emerging picture of early diversification in this obligate mutualism is consistent with a model of pollinator colonization of partly diversified host plant taxa, rather than strictly congruent diversification.
Phylogenetic reconstruction was based on a 2.1-kb stretch of mtDNA, including most of the cytochrome oxidase I and II subunits and an intervening tRNAlys region (positions 1495–3603 in the Drosophila yakuba genome; ref. 16). The region was PCR-amplified, and both strands were sequenced either manually with Sequenase T7 polymerase or with Amersham Dye-terminator chemistry on an Applied Biosystems 377 automated sequencer following the manufacturer’s protocols (Amersham Pharmacia). All samples used for the analyses are listed in Table 1. For Prodoxidae, all but one genus were included, with a sufficient number of taxa for each genus to provide a stable topology. The unavailable genus Prodoxides from South America, known only from the type series, is basal or near-basal in the family (17), and its omission does not affect the current analyses. For the pollinating yucca moths and the derived nonpollinating cheater yucca moths, all taxa used in ref. 11 were included, resulting in representatives of all major life history categories. Outgroup taxa included single species from the three incurvarioid families (Adelidae, Heliozelidae, Cecidosidae) and four species from the Incurvariidae.
Samples used for the analyses
Relationships of the incurvarioid families have been established from morphological studies (17) (Fig. 1). Initial parsimony, neighbor-joining, and likelihood analyses of the complete COI-COII data set supported this phylogeny, except that the single available cecidosid, Cecidoses eremita, clustered with the basal Heliozelidae in all analyses. This finding also held true in an analysis based on amino acids for the COI and COII regions of the data. We hypothesized that this result was caused by an accelerated substitution rate in these lineages, leading to long branch attraction. This hypothesis was tested by using a parametric bootstrap procedure (18). Parameters for a five-taxon tree were estimated by assuming the topology in Fig. 1 and used as input data for the simulation program seq-gen, version 1.1 (19). After removing all third codon sites, maximum-likelihood estimates of base frequencies, transition bias (Tr/Ti), the distribution of substitution rates (approximation of the gamma parameter assuming four rate categories), and branch lengths were calculated on a tree constrained to have the topology presented in Fig. 1. Average stem node-to-tip branch lengths were estimated for Incurvariidae and Prodoxidae, giving single length estimates for each branch in the constraint tree (Fig. 2). By using the estimated parameters, 100 five-taxon data sets with sequence lengths equal to the analyzed data (1,182 bp) were generated under the evolutionary model of Hasegawa et al. (20) and Yang (21) and analyzed by using the unweighted parsimony algorithm in paup* (22). When the resulting most-parsimonious (MP) trees were tallied by topology, it became clear that the long branch leading to C. eremita was in fact causing problems. The true input topology was recovered only from 16% of the simulated data sets (Fig. 3). For this reason, we removed the C. eremita sequence from further analyses.
Proposed phylogeny for the families of Incurvarioidea, based on morphological data (17). Dashed line indicates that taxon was unavailable for DNA analysis, and placement is based on morphology alone.
Input parameters (topology, branch lengths, transition bias, and gamma parameter) used to generate simulated sequence data sets for the parametric bootstrap analysis. Branch lengths leading to prodoxids and incurvariids are averages for all sampled (see text).
Frequency of parametric bootstrap replicates from which each topology was recovered as a single MP tree, one of two MP trees, or one of three MP trees. The true topology (lane A) was only recovered from 10 of the 100 simulated data sets. A = Adelidae, H = Heliozelidae, I = Incurvariidae, C = Cecidosidae, P = Prodoxidae.
After removing the Cecidoses sequence, topologies were estimated with paup* (22) by using unweighted parsimony, maximum likelihood, and neighbor joining. The maximum-likelihood and neighbor-joining trees were estimated assuming the nucleotide substitution model of Hasegawa et al. (20). Bootstrap support values for each node were estimated in the parsimony and neighbor-joining analyses. Given the low level of sequence divergence among members of the Tegeticula yuccasella complex, all codon positions were used for initial estimates of the topology.
Likelihood ratio tests were used to establish the best-supported model of molecular evolution for estimating minimum ages for the key life history traits. Rate heterogeneity across sites and among branches was tested. The substitution model of Hasegawa et al. (20) was assumed. As with the simulation experiment, third codon positions were dropped from the analysis to avoid the effects of substantial saturation. Standard likelihood ratio tests were applied in which the null hypothesis is a special case of the alternative hypothesis (23, 24). For example, the model that specifies a single substitution rate for all branches in the phylogeny is a special case of the model that estimates a separate rate parameter for each branch. In such cases, where the null hypothesis is nested within the alternative, two times the negative log of the likelihood ratio statistic: 
approximates the χ2 distribution with the degrees of freedom equal to the difference in the number of parameters estimated under the null and alternative hypotheses (24).
Finally, boundary ages were estimated for key events in the evolutionary history within the Prodoxidae. The topology shown in Fig. 4 was assumed. At each node, we used the two-cluster test of Takezaki et al. (26) to test the null hypothesis of equal substitution rates for the two clades above. When the null hypothesis was not rejected, the “height” of the node was estimated as half the average genetic distance between pairs of taxa sampled from the separate clades above the node (26). The relative age of the node then was estimated as quotient of the height of node of interest divided by the height of the reference node that defines the split between Incurvariidae and Prodoxidae+Cecidosidae (Fig. 1). Standard errors on the age estimates were estimated from the variance in this ratio (27). Covariances between node heights were estimated as described by Ayala et al. (28). Finally, the minimum age of each node was estimated as the relative age multiplied by 95 Mya, the minimum age of the split between the Incurvariidae and Prodoxidae+Cecidosidae clades.
Phylogeny for Prodoxidae and three outgroup families; shown is one of 10 MP trees, with neighbor-joining and maximum-likelihood trees having identical topology. Numbers are bootstrap values based on 100 iterations. Moth families are given on right; names for unnamed taxa within the yuccasella complex are those given in ref. 11 and are based on host yucca and state (all of which are described in ref. 25). The terms early and late cheaters reflect use in ref. 11 and refer to the host phenological stage attacked. Asterisks identify taxa that were excluded in age estimations caused by significant rate deviations in two-cluster tests. Arrows indicate latest origin of four named traits to be explored.
Phylogenetic analysis of the entire COI-COII data set resulted in 10 MP trees that differed only in the placement of taxa within the yuccasella complex. One of the 10 MP trees was recovered in the neighbor-joining analysis of the full data set, assuming no rate variation among sites (Fig. 4). The maximum-likelihood tree recovered under the single rate model differed in that Vespina was pulled out of the Incurvariidae clade and placed below it. Both the maximum likelihood and MP/neighbor joining trees are consistent with an earlier phylogeny (15), including far fewer taxa, except that the basal genera within Prodoxidae (Greya and Lampronia+Tetragma) have switched positions. This difference holds up in all analyses of the current data. More importantly, inferences concerning yucca moth diversification are unaffected by alternative topological resolutions of the basal prodoxid genera.
Under the assumption that the tree in Fig. 4 represents the true topology, the maximum-likelihood estimate of the gamma-distribution parameter was 0.18. The likelihood for the model including rate heterogeneity was significantly better than that of the single rate model (−2 log Λ = 3,406, P ≪ 0.01). The estimated value of gamma was very robust, remaining unchanged given alternative resolution of poorly supported nodes within the topology. The observed rate heterogeneity across sites was accounted for in the age estimations.
The likelihood ratio test for rate variation among branches was run, assuming the topology shown in Fig. 4 and rate variation among sites. The molecular clock assumption (no rate variation among branches) was rejected (−2 log Λ = 134, P = ≪ 0.01). Given this result, we turned to the distance-based tests of Takezaki et al. (26). The two-cluster test identified rate variation between clades at two nodes. The node-to-tip test identified taxa within these nodes as having significantly higher (two Prodoxus species) or lower (Lampronia) than average substitution rates. Therefore, these three taxa were not included in the age estimation analyses.
We proceeded to estimate ages for nodes of particular significance in prodoxid life history evolution. Both biogeographic and fossil data support a minimum age of 95 My for the prodoxid stemgroup (Fig. 5). Cecidosid moths combine a classical Gondwanan distribution with very low dispersal ability (29). Six genera within the family are gall-makers on Anacardiaceae, and they make up separate monophyletic groups in South America and southern Africa (30). Their shared host utilization suggests an origin before breakup of West Gondwana, which was definitive for all but very able dispersers 95–100 Mya (31). This date substantially predates the proposed 70-My age of the Anacardiaceae (32, 33), but is consistent with other biogeographic evidence (34). Further support for old age of the Cecidosidae is a recently discovered, unnamed genus (35) endemic to New Zealand, now regarded as basal within the family (R. Hoare and E.S. Nielsen, personal communication), indicating family origin before the isolation of New Zealand, approximately 82 Mya (36). Cecidosid diversification thus indicates a conservative age of >95 My for the prodoxid-cecidosid divergence. Because the Cecidosidae has evolved at a faster rate than the other families, however, it cannot be used for calibration. Instead we applied the 95-My estimate to the next lower node, between the Incurvariidae and Prodoxidae+ Cecidosidae clades.
Phylogeny for Incurvariidae, Cecidosidae, and Prodoxidae, indicating key geological and biological events, and biogeographic benchmark dates used to infer minimum age of stemgroup Prodoxidae.
The fossil record is scant, but consistent with the age calculations based on biogeographic data. There are no prodoxid fossils, but 97-My-old leaf fossils with highly characteristic mines of the more derived group Ditrysia provide indirect evidence that the incurvarioid stemgroup was present at that time (37, 38), whereas a proposed 110- to 120-My-old incurvarioid wing (39) cannot be exclusively tied to the superfamily (40).
After age estimation for the family, we calculated minimum ages for nodes associated with major evolutionary events within the Prodoxidae (Fig. 6). Basal genera are confined to humid and semiarid habitats and feed on a variety of dicotyledonous families. Coincident with the invasion of arid habitats by Prodoxidae, woody monocots were colonized no later than 44.1 Mya, with extant members of the basal genus Mesepiola feeding on seeds of Nolinaceae. A very short internode to the Prodoxus+pollinator node indicates rapid colonization of the yuccas by 41.5 Mya. The pollinating genera Parategeticula and Tegeticula originated very closely in time 35–41 Mya. The estimated age of the node separating these two genera is younger than the basal node within Tegeticula, but not significantly so. In fact, all three genera feeding on Yucca arose so quickly that their ages overlap in the analysis. The distinct feeding and life history characteristics exhibited by Prodoxus, Parategeticula, and Tegeticula seem to have arisen in rapid succession shortly after the colonization of woody monocots (Fig. 6). Parsimonious trait reconstruction on the phylogeny places the colonization of woody monocots on the internode between the Lampronia+Tetragma and Mesepiola branches. and is consistent with a late burst of lineages quickly diversifying into all extant life histories. This reconstruction is consistent with the hypothesis that a burst of diversification occurred as the woody monocots were colonized, yielding all of the extant life histories observed in the monocot feeding prodoxids. Alternative scenarios of earlier colonization of the woody monocots require an early period of limited cladogenesis followed by explosive radiation, or consistent rate of cladogenesis with subsequent extinction of basal monocot feeders.
Phylogeny for Prodoxidae as in Fig. 4, showing estimated minimum ages (± SE in parentheses) for select nodes. Taxa are given on the right, with nested brackets within the yuccasella complex indicating cheater yucca moth clades. Prodoxid hosts are given at far right.
After the early rapid diversification of prodoxids on Yucca, a subsequent explosive radiation occurred within the T. yuccasella complex 3.2 ± 1.8 Mya, giving rise to all pollinator lineages on capsular-fruited yuccas and the nonpollinating cheater yucca moths so quickly that the exact topology among major branches cannot be established from data presented here. However, additional data on nucleotide variation among host specific populations show definitively that there have been at least two independent losses of pollinating behavior in the T. yuccasella complex (unpublished work). By using maximum sequence divergence within the two independent cheater lineages, we find that two losses of pollination behavior in these two lineages is simultaneous given the coarse resolution of our age estimates. The early and late cheater lineages identified by Pellmyr et al. (11) arose 1.26 ± 0.96 and 1.26 ± 0.95 Mya, respectively.
The availability of a timeline for diversification of the yucca moths makes possible the erection of explicit predictions for patterns of host diversification under different models of how this plant-moth association has evolved (25, 41). Specific models that generate different predictions include strict cospeciation, synchronous diversification without cospeciation, and asynchronous diversification without cospeciation. In a strict cospeciation model, which is the most constrained model, parallel topologies are predicted between the yuccas and the moths, and diversification rates also should be synchronous. This is the predicted outcome if host shifts are very unlikely, and diversification strongly linked between the taxa. This model already can be rejected (4, 11). If moths can colonize new hosts, but still respond evolutionarily to host diversification, synchronous patterns of diversification are predicted for the moths and their hosts, but their topologies are not predicted to be congruent. In the least constrained model, where host switching is common and factors driving diversification of yuccas and yucca moths are unlinked, no correlations are expected in either relative rates or timing of plant and pollinator diversification.
Information about phylogeny and history of diversification of yuccas is still limited (42), preventing strong tests of the above predictions, but some inferences about the history of the association can be drawn. The estimated Eocene age of the association between yuccas and prodoxids provides indirect evidence of far older age of Yucca than previously documented. The macrofossil record (excluding Pleistocene subfossils; ref. 43) is limited to a 14-My-old vegetative fragment sharing one synapomorphy with extant Yucca (44). Applying a molecular clock to rbcL data, Eguiarte (45) proposed that Agavaceae and Nolinaceae diversified about 47 Mya, a date slightly older than that inferred here for basal diversification of prodoxids onto these families. Geologically, this period coincides with the onset of uplifting of western North America (46–48), which led to widespread development of semiarid habitats where extant yuccas are constituent taxa. Alternatively, Axelrod (49) speculated that Yucca may have evolved as early as late Cretaceous. Although a proposed Eocene origin would not exclude simultaneous diversification of moths and plants, a Cretaceous origin would indicate that early yuccas persisted without moths as pollinators, and that the yucca-yucca moth association originated as moths colonized extant hosts rather than through parallel diversification.
The second explosive radiation of yucca moths, within the yuccasella complex, occurred 3.2 ± 1.8 Mya. The moths derived from this radiation were primarily those of capsular-fruited yuccas, which occur in the northern portion of the Yucca geographic range. The onset of this radiation coincided with rapid aridification (50), leading to true deserts and treeless steppes, thus creating or extending the primary extant habitats for the capsular yuccas and their moth associates. There is currently neither fossil nor phylogenetic information for the history of Yucca within this Pliocene window. Rigorous tests of the predictions above must await a robust plant phylogeny.
Obligate mutualisms have attracted particular interest in terms of their stability and endurance (10, 51). The presence of pollinating yucca moths for at least 40 My can be compared against available data for another obligate pollination mutualism, between fig wasps and figs. A fossil fig wasp of Oligocene age is known (52), but the presence of fig fruits from the Eocene onward (53) suggest that this mutualism has persisted as long as the yucca-yucca moth association. Furthermore, many associations between specific yuccas and yucca moths have persisted since the evolution of derived cheater yucca moths, which dramatically increased the cost to the plants (11), suggesting that these mutualisms retain evolutionary stability over long time spans.
We thank J. Donahue, T. Friedlander, M. Gentili, J.-F. Landry, R. Leuschner, Z. Mazanec, N. Ryrholm, and D. Wagner for moth samples and E.A. Herre for references to fossil fig wasps. Collecting permits were provided by Big Bend National Park and the Texas Department of Parks and Wildlife. Work was funded by National Science Foundation Grants 95–09056 and 95–28072 and a grant from the National Geographic Society.
