An extreme case of plant-insect codiversification: figs and fig-pollinating wasps.

It is thought that speciation in phytophagous insects is often due to colonization of novel host plants, because radiations of plant and insect lineages are typically asynchronous. Recent phylogenetic comparisons have supported this model of diversification for both insect herbivores and specialized pollinators. An exceptional case where contemporaneous plant-insect diversification might be expected is the obligate mutualism between fig trees (Ficus species, Moraceae) and their pollinating wasps (Agaonidae, Hymenoptera). The ubiquity and ecological significance of this mutualism in tropical and subtropical ecosystems has long intrigued biologists, but the systematic challenge posed by >750 interacting species pairs has hindered progress toward understanding its evolutionary history. In particular, taxon sampling and analytical tools have been insufficient for large-scale cophylogenetic analyses. Here, we sampled nearly 200 interacting pairs of fig and wasp species from across the globe. Two supermatrices were assembled: on an average, wasps had sequences from 77% of 6 genes (5.6 kb), figs had sequences from 60% of 5 genes (5.5 kb), and overall 850 new DNA sequences were generated for this study. We also developed a new analytical tool, Jane 2, for event-based phylogenetic reconciliation analysis of very large data sets. Separate Bayesian phylogenetic analyses for figs and fig wasps under relaxed molecular clock assumptions indicate Cretaceous diversification of crown groups and contemporaneous divergence for nearly half of all fig and pollinator lineages. Event-based cophylogenetic analyses further support the codiversification hypothesis. Biogeographic analyses indicate that the present-day distribution of fig and pollinator lineages is consistent with a Eurasian origin and subsequent dispersal, rather than with Gondwanan vicariance. Overall, our findings indicate that the fig-pollinator mutualism represents an extreme case among plant-insect interactions of coordinated dispersal and long-term codiversification. [Biogeography; coevolution; cospeciation; host switching; long-branch attraction; phylogeny.].

Processes affecting the diversification of insects are crucial to understanding the origin of biodiversity, because most animals are either insect herbivores, or natural enemies (predators or parasitoids) of these phytophages (Novotny et al. 2002). As primary consumers, most insect herbivores are involved in antagonistic interactions with plants and, although herbivores often exhibit host-specific coevolutionary adaptations to plant defenses (Ehrlich and Raven 1964), recent empirical studies have suggested that host plant lineages are generally older than their associated herbivores (Percy et al. 2004; Tilmon 2008; McKenna et al. 2009). Such patterns of asynchronous plant–insect diversification are consistent with the general paradigm that insect speciation results from colonization of novel host plants and subsequent reproductive isolation (Percy et al. 2004; Tilmon 2008; McKenna et al. 2009; Fordyce 2010).

Phytophagous insects are often enemies of plants, but some engage in beneficial pollination mutualisms. A charismatic example involves the ca. 750 species of figs (Ficus, Moraceae) and their pollinating wasps (Hymenoptera, Chalcidoidea, Agaonidae) (Fig. 1). Agaonids are the only pollen vectors for fig trees and agaonid larvae feed exclusively on the flowers of their Ficus hosts. Each partner is thus entirely dependent on the other for reproduction. Figs are also a major resource for frugivores and most animal-dispersed tropical tree species interact with vertebrates that also consume figs (Howe and Smallwood 1982). The fig-pollinator mutualism is therefore ecologically important in most tropical ecosystems (Shanahan et al. 2001). Many fig species reproduce irregularly, are relatively inaccessible in the forest canopy, or today are found only in rainforest remnants, such that coordinated sampling of Ficus and pollinator species for systematic study is difficult. These sampling challenges, coupled with the limitations of analytical tools for large data sets, have hindered progress toward understanding the global evolutionary history of the mutualism, despite the fact that many details of this intricate symbiosis were described almost a century ago (Janzen 1979; Wiebes 1979; Weiblen 2002; Cook and Rasplus 2003; Herre et al. 2008).

Figure 1

Classification and worldwide distribution of Ficus. The numbers of species per subgenus is represented as a proportion of total Ficus species richness. Breeding systems are indicated as either monoecious (M) or dioecious (D) and modes of pollination are indicated as passive (P) or active (A). *Agaon, Alfonsiella, Allotriozoon, Courtella, Elisabethiella, Nigeriella and Paragaon. **Deilagaon and Waterstoniella.

Species-specificity in fig pollination appears to be extreme compared with most other insect pollination mutualisms. Most fig species are pollinated by only one or a few wasp species and most wasps are associated with just a single fig species (Cook and Rasplus 2003; Molbo et al. 2003; Cook and Segar 2010). Pollinators are specifically attracted to volatile compounds emitted by figs (Hossaert-McKey et al. 1994) and access to the specially modified inflorescences is by means of distinctive mandibular appendages and detachable antennae (van Noort and Compton 1996). Pollination is either active (two-thirds of the fig species) or passive (one-third, mostly within subgenera Pharmacosycea, Ficus, Synoecia, and Urostigma) (Kjellberg et al. 2001). Active agaonid wasps collect pollen from the anthers of their native figs and store it in thoracic pollen pockets (Galil and Eisikowitch 1968; Ramirez 1978). Once inside a receptive fig, they remove pollen from their pockets and deposit it on the flower stigma each time they lay an egg (Galil and Eisikowitch 1968; Kjellberg et al. 2001). Passively pollinated figs produce large quantities of pollen through anther dehiscence and wasps are covered with pollen (Galil and Neeman 1977) before flying away from their natal figs.

Closely matching fig and pollinator traits might be products of coadaptation (Ramirez 1974; Wiebes 1979, 1982a; Kjellberg et al. 2001; Weiblen 2004) but, regardless, trait-mediated interactions have the potential to simultaneously affect the evolution of reproductive isolation among pollinator and fig populations; this is because fig wasps breed exclusively in pollinated figs. This line of reasoning has underpinned the hypothesis that cospeciation might account for patterns of fig and pollinator diversity. However, this notion runs contrary to the paradigm that insect speciation generally involves host-switching (Tilmon 2008) and so it remains a controversial proposition that requires rigorous testing.

Under the cospeciation scenario, phylogenies of figs and pollinators are expected to show substantial congruence. There is some evidence for this pattern (Herre et al. 1996; Machado et al. 2005; Rønsted et al. 2005; Cook and Segar 2010; Cruaud et al. 2011a), but recent studies have countered the underlying case for cospeciation with evidence of cryptic wasp species and relaxed partner specificity. At least 50 fig species are now known to have multiple pollinator species (Michaloud et al. 1985, 1996; Rasplus 1996; Kerdelhué et al. 1997; Lopez-Vaamonde et al. 2002; Greeff et al. 2003; Molbo et al. 2003; Haine et al. 2006; Moe and Weiblen 2010; Chen et al. 2012) and as many as 4 different wasp species are known to pollinate a single fig species (Machado et al. 2005; Cook and Segar 2010). Such cases occur in a broad taxonomic and geographic spectrum, although cases of pollinator species sharing multiple fig species have been reported mostly from monoecious figs in the Neotropics (Molbo et al. 2003) and the Afrotropics (Erasmus et al. 2007; Cornille et al. 2012; McLeish and van Noort 2012). In any event, evidence of relaxed host specificity and some incongruent fig-pollinator phylogenies (Machado et al. 2005) suggest that host shifting is a viable alternative explanation for fig-pollinator diversification.

Cospeciation has been hypothesized for the vertically transmitted endosymbionts of insects (e.g., Moran 2001; Jousselin et al. 2009) but this is not a plausible general model for the evolution of plant–insect associations, which are horizontally transmitted and not so integrated metabolically. Further, if the plant traits that mediate insect associations happen to be phylogenetically conserved, then host shifting among close relatives could also result in topologically congruent phylogenies (Percy et al. 2004). In addition, historical biogeography has the potential to confound the explanation of such patterns if synchronous plant–insect dispersal to new environments, followed by geographic isolation, results in cospeciation.

Another useful approach is to investigate patterns of temporal congruence (Page and Charleston 1998). Divergence time estimates for fig and pollinator clades are expected to be approximately equal in the event of coradiation, whereas insect lineages are expected to be younger than hosts in the case of host shifting (Percy et al. 2004; Tilmon 2008; McKenna et al. 2009).

Previous comparisons of fig and pollinator phylogeny have yielded rather different insights on the relative importance of host shifting and codiversification depending on the taxonomic scope of sampling (Cook and Segar 2010). Molecular phylogenetic trees appear roughly parallel when based on exemplars of Ficus sections and wasp genera (Herre et al. 1996; Jackson 2004; Cruaud et al. 2011a), but such deep taxonomic sampling is unlikely to detect host shifts among close relatives (Machado et al. 2005). On the other hand, regional comparisons of particular fig and pollinator clades have tended to reject cospeciation in favor of host-switching (Machado et al. 2005; Marussich and Machado 2007; Jackson et al. 2008; Jousselin et al. 2008), although not always (Weiblen and Bush 2002; Silvieus et al. 2008). A global test for codiversification therefore requires dense sampling of many fig and pollinator lineages across the entire geographic range, but a problem of this magnitude poses a further methodological challenge.

Tests of cophylogenetic hypotheses often employ tree reconciliation methods that infer evolutionary processes such as cospeciation, host shifts, duplications, and losses to account for topological incongruence between host and associate phylogenies (Page 1994). This approach has the power to model the relative contributions of different evolutionary processes to a given phylogenetic pattern, but biologically realistic scenarios become computationally intractable for large numbers of taxa (Merkle and Middendorf 2005; Ovadia et al. 2011). Genetic algorithms that incorporate dynamic programming to efficiently locate and evaluate samples from an extremely large universe of event-based solutions hold promise in this regard (Conow et al. 2010).

Here, we extended the application of a genetic algorithm to event-based tree reconciliation analysis for cophylogenetic problems involving >100 taxon pairs and applied randomization tests involving null models to test the codivergence hypothesis on an unprecedented scale. Nearly 200 pairs of interacting fig and fig wasp species were sequenced at 5 fig loci (providing up to a total of 5.5 kb DNA sequence) and 6 wasp loci (up to a total of 5.6 kb). Two supermatrices were assembled. On an average, wasps had sequences from 77% of 6 genes, figs had sequences from 60% of 5 genes, and overall, we generated 850 new DNA sequences for the purpose of this study. Maximum likelihood (ML) analyses of fig and wasp data sets and Bayesian phylogenetic analyses under relaxed molecular clock assumptions enabled the comparison of distance, event-based, and temporal congruence. Inferences from historical biogeography based on our global sample of fig and pollinator clades provided additional insight on the relative roles of dispersal and vicariance with respect to alternative hypotheses of diversification.

Materials and Methods

Taxonomic Sampling and DNA Sequencing

Ficus—We sampled 200 fig species (>1/4 of the circa 750 described species) that represent all Ficus sections recognized by Berg and Corner (2005) (Appendix S1 in the Supplementary Material Online, doi: 10.5061/dryad.hr620). Four taxa belonging to the tribe Castilleae s.l., Antiaropsis decipiens, Castilla elastica, Poulsenia armata, and Sparattosyce dioica, were included as outgroups (Datwyler and Weiblen 2004; Rønsted et al. 2005; Zerega et al. 2005; Clement and Weiblen 2009). Total genomic DNA was extracted from 20–30 mg of dried leaf-fragments or herbarium material following Rønsted et al. (2008). Ficus phylogeny was reconstructed using 5 genes: ITS (891 bp), ETS (528 bp), glyceraldehyde 3-phosphate dehydrogenase (G3pdh, 769 bp), chloroplast expressed glutamine synthetase region (ncpGS, 1630 bp), and granule-bound starch synthase (waxy region, 1734 bp).

Amplification of ITS, ETS, and G3pdh was performed following Rønsted et al. (2008). The ncpGS region (Emshwiller and Doyle 1999) was amplified using Moraceae-specific primers 3F (5′ GTT GTG ATT WAC CAT GCT) and 4R (3′ AGA TTC AAA ATC GCC TTC) designed for this study. Amplification of ncpGS consisted of 4 min at 94°C followed by 36 cycles of: 1 min denaturation (94°C), 1 min annealing (50°C), and 2-min extension (72°C). After the last cycle, the temperature was kept at 72°C for a final 5-min extension and then lowered to 4°C. The GBSSI or waxy region (Mason-Gamer et al. 1999; Clement 2008) was amplified using Moraceae-specific primers 3F (5′ GAT CGY GTG TTT GTR GAT CAC C) and 10R (3′ GCA ACT GAA TGA GAC CAC A). Amplification of waxy consisted of 3 min at 94°C followed by 2 cycles of 94°C for 1 min, 58°C for 1 min, 72°C for 2 min, 2 cycles of 94°C for 1 min, 56°C for 1 min, 72°C for 2 min, 2 cycles of 94°C for 1 min, 54°C for 1 min, 72°C for 2 min, 2 cycles of 94°C for 1 min, 50°C for 1 min, 72°C for 2 min, and 24 cycles of 94°C for 1 min, 48°C for 1 min, 72°C for 2 min. After the last cycle, the temperature was kept at 72°C for a final 20-min extension and then lowered to 4°C. Amplified products were purified with the Qiagen PCR purification kit (Qiagen Inc.) following the manufacturer's protocols. ITS, ETS, G3pdh, and ncpGS were sequenced directly from PCR products whereas waxy was cloned using the TOPO-TA PCR cloning kit (Invitrogen, Carlsbad, CA). Nine clones were screened for inserts, and plasmids were isolated from 3 of these using the Qiagen plasmid prep kit. Multiple copies of waxy are known in the Rosales (Evans et al. 2000) and therefore it was necessary to ensure that phylogeny reconstruction was performed with orthologous copies. Two copies have been detected in Moraceae, GBSS1 and GBSS2, which were easily distinguished on the basis of size and intron alignment (Silvieus et al. 2008; Clement W., unpublished data). Analyses were based solely on GBSS1 because GBSS2 was encountered less commonly in figs.

Cycle sequencing reactions were carried out following Rønsted et al. (2008). For sequencing of the ncpGS region, internal primers 1F (5′ TCW TGW GCT GAA AAG CAT), 2F (TTT AAT CTC CAG ACT CSA), and 5F (5′ TAG TTC ACT CTA AAG GGT) were designed for this study in addition to the primers used for amplification. Some 50% of the sequences were obtained from de novo sequencing for the purpose of this study and have been deposited in GenBank (Appendix S1 in the Supplementary Material Online). Other sequences, mostly deposited by coauthors, were obtained from existing databases.

Agaonidae.—Ninety-three percent of the 200 wasps and figs used in this study are true associates; i.e., even if they were not collected together simultaneously, the agaonid species is the pollinator of the fig species. In the very few cases where the corresponding agaonid was not available in our collection, we used instead the pollinator of a closely related species of fig (Appendix S2 in the Supplementary Material Online). This was always a wasp species that was a close congener of the actual pollinator. As the phylogenetic position of Agaonidae within the large and complex superfamily Chalcidoidea is as yet unknown (Gibson et al. 1999; Munro et al. 2011), 4 divergent members of the superfamily served as outgroups: Sycophaga (Sycophaginae), Ficomila (Eurytomidae), Megastigmus (Torymidae), and Trichogramma (Trichogrammatidae).

All material was collected alive in the field and fixed in 95% ethanol. With very few exceptions, Agaonidae sequences were obtained from the nondestructive extraction of a single wasp specimen (corpse kept as voucher). DNA was extracted from a single individual that was incubated at 56°C overnight (with gentle “shaking” steps by inverting the tubes) and using the Qiagen DNeasy kit following the manufacturer's protocol. When destructive extraction was used, vouchers were selected among specimens sampled from the same tree and the same fig after careful identification by J.Y.R., Sv.N., and R.U. Vouchers are deposited at CBGP, Montferrier-sur-Lez, France. To infer phylogenetic relationships between agaonid species, we combined 2 nuclear protein-coding genes [F2 copy of elongation factor-1a (EF1a, 516 bp), Wingless (Wg, 403 bp)]; 2 mitochondrial protein-coding genes [cytochrome c oxidase subunit I (COI, 1536 bp), cytochrome b (Cyt b, 749 bp)], and 2 ribosomal genes [28S rRNA (D2–D3 and D4–D5 expansion regions, 1520 bp), 18S rRNA (variable regions V3–5, 787 bp)]. Extraction, amplification, and sequencing protocols follow Cruaud et al. (2010) for CytB, COI (barcode fragment), Wg, 28S, and 18S rRNA, Weiblen (2001) for COI [C1-J-2183 (Jerry)—TL2-N-3014 (Pat) fragment], and Cruaud et al. 2011b for EF1a. Both strands for each overlapping fragment were assembled using Geneious v5.4.2 (Drummond et al. 2007).

Sixty-seven percent of the sequences were obtained from de novo sequencing for the purpose of this study and have been deposited in GenBank (Appendix S2 in the Supplementary Material Online). Other sequences, mostly deposited by coauthors, were obtained from public databases.

Phylogeny Reconstruction

Protein-coding genes and hypervariable regions were aligned using ClustalW 1.81 with the default settings (Thompson et al. 1994). Alignments of protein-coding genes were translated to amino acids using Mega 4 (Tamura et al. 2007) to detect frameshift mutations and premature stop codons, which may indicate the presence of pseudogenes. Alignment of sequences encoding rRNA was based on secondary structure models (Gillespie et al. 2006), following Cruaud et al. (2010). Phylogenetic trees were estimated using both ML and Bayesian methods. We selected separate models of molecular evolution for different genomic regions including mitochondrial genes, rRNA stems, rRNA loops + regions of ambiguous alignment, and individual nuclear genes using the Akaike information criterion implemented in MrAIC.pl 1.4.3 (Nylander 2004).

For each data set, we performed ML analyses and associated bootstrapping (1000 replicates) using the MPI-parallelized RAxML 7.0.4 software (Stamatakis 2006b). GTRCAT approximation of models (Stamatakis 2006a) was used for ML boostrapping (1000 replicates). Bootstrap percentage (BP) >70% was considered as strong support (Felsenstein and Kishino 1993). Bayesian analyses (BA) were conducted using a parallel version of MrBayes 3.1.1 (Huelsenbeck and Ronquist 2001). We assumed across-partition heterogeneity in model parameters by considering the parameter m (Nylander et al. 2004; Marshall et al. 2006; McGuire et al. 2007). Parameter values for the model were initiated with default uniform priors and branch lengths were estimated using default exponential priors.

To improve mixing of the cold chain and avoid it converging on local optima, we used Metropolis-coupled Markov chain Monte Carlo (MCMCMC) simulation with each run including a cold chain and 3 incrementally heated chains. The heating parameter was set to 0.02 in order to allow swap frequencies from 20% to 70%. For both figs and pollinators, we ran 2 independent runs of 30 million generations. All the values were sampled every 3000 generations. For the initial determination of burn-in, we examined the plot of overall model likelihood against generation number to find the point where the likelihood started to fluctuate around a constant value. Convergence of the chains was evaluated using the online application AWTY (Nylander et al. 2008) and the results were based on the pooled samples from the stationary phases of the 2 independent runs. Given that posterior probabilities (PP) may overestimate clade support, for reasons discussed elsewhere (Suzuki et al. 2002; Cummings et al. 2003; Erixon et al. 2003; Simmons et al. 2004), only clades with PP > 0.95 were considered strongly supported. All analyses were conducted on a 150core Linux Cluster at CBGP, Montferrier-sur-Lez, France.

Test of alternative hypotheses.—To assess whether certain alternative relationships among recovered clades could be statistically rejected, we performed AU (Shimodaira 2002) and SH (Shimodaira and Hasegawa 1999) tests in the CONSEL package (Shimodaira and Hasegawa 2001). The program makermt was used to generate K = 10 sets of bootstrap replicates (r1 = 0.5, r2 = 0.6, r3 = 0.7, r4 = 0.8, r5 = 0.9, r6 = 1, r7 = 1.1, r8 = 1.2 r9 = 1.3, r10 = 1.4). Each set consisted of 100 000 replicates of the row sums (10 times the default number of replicates). RAxML was used to compute the per-site log likelihoods for all topologies tested. To assess the relative support for competing phylogenetic hypotheses, we also conducted AU and SH tests on recently published data sets (Lopez-Vaamonde et al. 2009; Cruaud et al. 2010), which placed Tetrapus as sister to all other agaonids with strong support (PP = 0.99 and BP = 55/59; PP = 1.00, respectively).

Effects of missing data.—There is debate in the literature as to the effect of missing data on the accuracy of phylogenetic analyses. Simulation results based on limited numbers of characters (Lemmon et al. 2009) indicated that nonrandom distributions of missing data can result in strong support for nodes that share no supporting characters. However, other empirical and simulation studies have concluded that taxa with extensive missing data can be accurately placed in phylogenetic analyses, and that adding characters with missing data is generally beneficial, if the overall number of characters is large and data are analyzed with appropriate methods (see Wiens and Morrill 2011). To assess the impact of missing data on our analyses, we performed 2 sets of additional analyses.

First, we built new (“complete species”) trees using only the more completely sequenced taxa (figs with more than 3 genes; wasps with more than 5 genes). Then, we used AU and SH tests to compare the full (all species and genes) tree, pruned of incompletely sequenced taxa, with the matching “complete species” tree. Second, we built new (“complete genes”) trees by removing gene fragments for which <60% of the taxa were available. We then used AU and SH tests to test if the full tree differed significantly from the “complete genes” tree. Taxa were pruned from the combined ML tree using the APE package (Paradis et al. 2004) in R 2.14.0 (http://www.R-project.org).

Bayesian relative rate tests and long-branch attraction artifact.—We tested constancy of evolutionary rates among agaonid species using both BEAST 1.5.3 (Drummond and Rambaut 2007) (coefficient of variation statistic and average rate for each branch of the chronogram, see “Molecular Dating” section) and a Bayesian relative rate (BRR) test (Wilcox et al. 2004). For the BBR test, the PP distributions of lengths for all branches from the most recent common ancestor (MRCA) of the ingroup to each of the terminal taxa were based on 1000 randomly chosen post-burn-in trees from the BA of the mitochondrial (mtDNA), nuclear (nuDNA) and combined data sets, respectively. Following Wilcox et al. (2004), we considered rates of evolution significantly different between 2 taxa if the 95% confidence interval of the PP distribution of the summed branch length did not overlap. Branch length estimates were compiled using Cadence v1.08b (Wilcox et al. 2004).

Long-branch attraction (LBA) artifacts can be difficult to detect, but methods have been proposed and we applied these to our data (Bergsten 2005). For computation time reasons, all additional analyses were conducted using RAxML only.

Removing first and third codon positions, which are fast evolving, may be a way to reduce LBA. However, this can also compromise tree resolution (Källersjö et al. 1999; Savolainen et al. 2002; Stefanovié et al. 2004). The RY-coding strategy (Woese et al. 1991), by discarding fast-evolving transitions and reducing compositional bias, constitutes a better approach (Phillips and Penny 2003; Philippe et al. 2005). We therefore compared the topologies obtained with or without RY-coding of: (i) the third (nt3) and (ii) first (nt1) and third mtDNA codon positions.

Long-branch extraction is another approach advocated for cases where LBA is suspected (Pol and Siddall 2001). Because LBA to the outgroups is the most frequent problem, analyses were conducted without the outgroups (Bergsten 2005).

Finally, the different sensitivity of parsimony and ML methods can help to detect if LBA is playing a major role (Brinkmann et al. 2005). We therefore performed parsimony analysis on our data set to detect potential shifts in position of agaonid groups. Parsimony analyses were conducted with TNT version 1.1 (Goloboff et al. 2008), using New Technology Search: 1000 replicates of random addition sequences, followed by random sectorial searches with default options, 100 cycles of ratchet and 3 rounds of tree-fusing. All substitutions were equally weighted and gaps treated as missing data. Robustness of topologies was assessed by bootstrap procedures using 1000 replicates.

Cophylogenetic Analyses

We tested the congruence between fig and wasp phylogenies using both distance and event/topology-based methods. The former generate patristic distance matrices between species in each phylogeny and then test for correlations between the 2 matrices. In contrast, event-based methods use evolutionary events [cospeciation, duplication, host-shifts, lineage sorting, and “failure to diverge” (Page and Charleston 1998; Charleston and Perkins 2006)] to map the associate phylogeny to the host one. A cost is assigned to each event type and we seek to find mappings that minimize the total cost. Statistical analyses can be performed by comparing the best costs found for the host–parasite data set against those of randomized instances.

We used the distance-based method, ParaFit, developed by Legendre et al. (2002) and implemented in the program CopyCat (Meier-Kolthoff et al. 2007). ParaFit evaluates the global hypothesis of host-associate cospeciation with a matrix permutation test of codivergence. This test combines 3 types of information: the associate phylogeny and the host phylogeny both described by their respective matrices of patristic distances, and the observed host-associate links. Each matrix representing associates and hosts is transformed into a matrix of principal coordinates. The association is then described by a new matrix, which includes both matrices of principal coordinates and the matrix of association. Patristic distances were computed from fig and wasp ML-phylogenetic trees. Tests of random association (null hypothesis) were performed using 9999 permutations globally across both phylogenetic trees. Although the distance-based approach is computationally simple, it only yields a measure of overall phylogenetic congruence and no information on the relative distribution of underlying evolutionary events that might have produced the pattern.

Event-based methods have the advantage of modeling evolutionary processes directly, but are computationally intensive (Charleston 2009). The problem of finding a mapping (event-based reconstruction) of minimum total cost has been shown to be computationally intractable (“NP-complete”) (Ovadia et al. 2011). Some existing software packages, e.g., TREEMAP 1.0; (Page 1994) and 2.02; (Charleston and Page 2002), use exhaustive searches, which are prohibitively slow and also permit only limited numbers of species. Other programs use heuristics (Merkle and Middendorf 2005), which are fast but may converge on suboptimal or invalid solutions (e.g., ancestral speciation inferred to have occurred after speciation of descendants nodes). For this reason, our analyses used a genetic algorithm to search a sample of the possible solution space with a dynamic programming step that efficiently evaluates the cost of each such sample. This approach, which finds solutions of near-optimal cost, was first implemented in the Jane software package (Conow et al. 2010) and was validated using a number of existing data sets in the literature (Libeskind-Hadas and Charleston 2009). However, the sheer size of our data sets put the analysis far beyond the computational limits of the original version of Jane, which also lacks support for randomization tests. We therefore substantially optimized and improved the existing Jane cophylogeny software package, resulting in a new system, Jane 2, which is capable of performing event-based analyses of very large data sets. Jane 2 and its tutorial are freely available for research and educational purposes at http://www.cs.hmc.edu/~hadas/jane/index.html.

We used Jane 2 with the following parameter values: the number of “generations” (iterations of the algorithm) was set to 40 and the “population” (number of samples per generation) was set to 1000. We explored 3 different cost models, each with 2 types of randomization test. The first model set costs per event as cospeciation = 0 and all other events = 1. This correponds to the TreeMap cost scheme so that a duplication event actually contributes two to the total cost because each of the 2 daughter lineages contributes one duplication event. The cost of the best solution was compared with the costs found in 100 randomizations in which the tip mappings were permuted at random, a method advocated by Aldous (2001). The second randomization involved 100 randomly generated pollinator trees, of the same size as the actual wasp pollinator tree, with random tip mappings. The random pollinator trees were constructed using the Yule model with beta parameter equal to −1.

In the second model, we used costs of 0 for cospeciation, 1 for each duplication, 1 for each host switch, and 2 for each loss event. In the third model, we set the cospeciation cost at -1 and all other costs to 0, where a negative cost maximizes the number of inferred cospeciations. For the second and third cost models, we used the same 2 randomization tests described for the first model.

All the analyses were performed at Harvey Mudd College (Claremont, CA) on a heterogeneous cluster of commodity computers comprising a total of 168 cores. On a single commodity computer (e.g., a dual core Macintosh), a single fig/wasp tree required ~3 h of computation and thus 100 randomized trials required several hours on our cluster.

Molecular Dating

We used the uncorrelated log-normal relaxed clock method implemented in BEAST 1.5.3 (Drummond and Rambaut 2007) and the same modeling strategies as for MrBayes and RAxML analyses. We assumed a Yule tree prior and we used default priors for all other parameters. We used 2 runs of 60 million generations with sampling every 6000 generations for figs, and 2 runs of 240 million generations with sampling every 24 000 generations for wasps. The 2 separate runs were then combined using LogCombiner 1.5.3. We ensured convergence using TRACER 1.5 (Drummond and Rambaut 2007). Following the removal of 10% burn-in, the sampled posterior trees were summarized using TreeAnnotator 1.5.3 to generate a maximum clade credibility tree and calculate the mean ages, 95% highest posterior density intervals (95% HPD) and PP. We used independent calibration points to estimate divergence ages of the main Ficus and agaonid clades. Following Rønsted et al. (2005), crown-group Ficus was assigned a uniform prior distribution with a minimum age of 60 Ma based on fossilized achenes (Collinson 1989) and a maximum age of 198 Ma based on converging molecular estimates for the origin of the angiosperms (Bell et al. 2005). Given uncertainties over the age of Dominican amber (Grimaldi 1994; Iturralde-Vinent and MacPhee 1996, 1999), crown-group Pegoscapus and Tetrapus were assigned uniform prior distributions with minimum ages of 15 Ma and maximum ages of 60 Ma based on Dominican amber fossil (Poinar 1993; Penalver et al. 2006). For both fig and wasp phylogenies, nodes including taxa endemic to La Réunion were modeled with a normal distribution with a mean of 8 Ma and SD of 1 Ma based on the proposed age for the Mascarene archipelago (McDougall and Chamalaun 1969; McDougall 1971).

Ancestral Area Reconstructions and Evolution of Pollination Mode

We inferred the evolution of pollination mode and the ancestral areas for figs and their pollinators using both ML and parsimony approaches implemented in Mesquite 2.73 (Maddison and Maddison 2008). Pollination modes and ancestral areas were inferred on the ML topologies. For ML optimization, we used a stochastic Markov model of evolution (Mk1). The Likelihood Decision Threshold was set to 2 log-likelihood units. Character data for Ficus and Agaonidae were obtained both from the literature (Kjellberg et al. 2001; Berg and Corner 2005) and from our examination of flowers, pollen pockets, and coxal combs. Following Lopez-Vaamonde et al. (2009), current species distributions were categorized into 4 character states: (0) Afrotropics, (1) Australasia, (2) Neotropics, (3) Eurasia. However, because several taxa occur in both Eurasian and Australasian regions and a couple of taxa occur in both Eurasian and Afrotropical regions, and Mesquite requires unique character states, we also defined 2 other states: (4) Australasia + Eurasia and (5) Afrotropics + Eurasia. We took into account all published geographic localities for Ficus and agaonids, museum specimens and about 3000 samples of fig wasp communities that we collected over the last 15 years. We also used the dispersal-extinction-cladogenesis model implemented in Lagrange (Ree and Smith 2008), using the same raw data and 4 character states. Dispersal rate between all areas was set to 1 during the whole period considered (data available upon request).

Results

DNA Sequence Data

The completeness of taxa in the combined data matrices is different for fig and wasps (Appendices 1–2 in the Supplementary Material Online and Supplementary Table S1). On an average, wasps have sequences from 77% of the 6 genes and 67% of the species were sequenced for at least 5 gene regions. On an average, figs have sequences from 60% of the 5 genes and 70% of the species were sequenced for at least 3 regions. Plastid regions provide little phylogenetic information within Ficus, enforcing the use of more informative single copy nuclear regions. These are known to be notoriously difficult to amplify from plants in general (Rønsted et al. 2007) and this was also the case for Ficus in the present study. Indeed, ncpGS and waxy matrices only include 24 and 23% of the taxa, respectively. Models chosen by MrAIC for each partition were as follows. Ficus data set: GTR + Γ (ETS, ITS, ncpGS, and waxy),GTR + I + Γ (G3pdh); Agaonid data set: GTR + Γ (mtDNA), GTR + I + Γ (EF1a, Wg, rRNA stems), HKY + I + Γ (rRNA loops). Given that α and the proportion of invariable sites can not be optimized independently from each other (Gu 1995) and following Stamatakis' personal recommendations (RAxML manual), we used GTR + Γ with 4 discrete rate categories for all partitions. As RAxML does not implement the HKY model, we used GTR instead.

Wasp Phylogeny

Our phylogenetic trees (Fig. 2 and Supplementary Fig. S1), reconstructed using ML and Bayesian approaches provide several new insights into the systematics of fig wasps.

Figure 2

BEAST chronograms of the evolutionary history of figs and fig wasps. Groups of figs and their associated genera of pollinators are represented using the same color. Ficus subgenera and Agaonidae subfamilies according to current classifications are delimited by colored rectangles (Pharma. for Ficus subgenus Pharmacosycea). Pie charts at main nodes show the likelihood of different geographic areas of origin as inferred by Mesquite (see “Materials and Methods” section). Gray rhombuses show clades of fig species from Continental Asia, whereas gray arrows indicate hypothesized southward migration of clades. Squares correspond to node supports: Black square: BP > 70% and PPMrBayes or PPBEAST> 0.95; white square: BP > 70% or PPMrBayes or PPBEAST> 0.95. Details in this figure can be viewed at greater magnification at Systematic Biology online.

Monophyly of the genera and intergeneric relationships.— Fifteen agaonid genera are recovered as monophyletic with strong support (Agaon, Alfonsiella, Allotriozoon Ceratosolen, Courtella, Deilagaon, Elisabethiella, Eupristina, Kradibia, Nigeriella, Pegoscapus, Pleistodontes, Tetrapus, Valisia, and Waterstoniella). In contrast, Platyscapa is polyphyletic, and Dolichoris is paraphyletic with respect to Blastophaga psenes (the pollinator of F. carica and type species of the genus Blastophaga), indicating the need for taxonomic rearrangements (Cruaud et al. 2012).

The relationships among the major clades are unclear (Fig. 2). BEAST analysis places Ceratosolen + Kradibia (subfamily Kradibiinae) as the sister group to the remaining Agaonidae with strong support (PPBEAST = 0.98), but this position is not strongly supported by MrBayes (PPMrBayes = 0.88) and ML analyses (BP = 43). BA place Tetrapus (monogeneric subfamily Tetrapusinae) nested within the Agaonidae with strong support (PPMrBayes = 1.00, PPBEAST = 1.00), although this position is only moderately supported by ML analyses (BP = 67).

Phylogenetic placement of the genus Tetrapus.—By not placing Tetrapus as sister to all other agaonids, our topology challenges all previous molecular studies by ourselves and others (Herre et al. 1996; Machado et al. 1996, 2001; Lopez-Vaamonde et al. 2009; Cruaud et al. 2010; Supplementary Table S3). This result deserves further examination, so we have conducted additional analyses on not only our current data set, but also previously published data sets. We provide here a summary of the main results (see the Appendix S3 in the Supplementary Material Online for further details):

Both AU and SH tests fail to reject alternative topologies in which either Tetrapus, or the clade of pollinators associated with subgenus Synoecia and subsection Frustescentiae (corresponding to Group 4 in Cruaud et al. 2010), is constrained to be the sister group to all other Agaonidae (Supplementary Table S2). Furthermore, AU and SH tests also fail to reject alternative positions of Tetrapus using 2 previously published data sets (Lopez-Vaamonde et al. 2009; Cruaud et al. 2010) that recover Tetrapus as sister to all other Agaonidae (Supplementary Table S2).

Tetrapus is recovered nested within the Agaonidae in all the analyses conducted to assess the impact of missing data on the accuracy of our phylogeny. AU and SH tests showed that phylogenetic trees pruned of incompletely sequenced taxa and trees built only on gene fragments for which at least 60% of the taxa were available, were not significantly different from the original trees including all available data (Supplementary Fig. S2B,C and Supplementary Table S2). Therefore, our analyses show that missing data are not responsible for the position of the genus Tetrapus.

Examination of branch lengths (mtDNA, nuDNA, and combined trees, Supplementary Figs. S1 and S3) indicates considerable variation in rates of molecular evolution among agaonid lineages This result is confirmed by the BRR tests (Supplementary Figs. S4 and S5) and the BEAST outputs (95% credible interval for the coefficient of variation of rates is not abutting against zero for each partition and covariance values span zero). Furthermore, a long branch leading to Tetrapus, is visible in both the nuDNA tree (Supplementary Fig. S3b) and ML and Bayesian combined trees (Supplementary Fig. S1). BRR tests and branch-specific rates inferred by BEAST reveal a lineage-specific increase in nucleotide substitution rates on this branch, and this is also the case for the branch leading to the outgroups (Supplementary Fig. S4).

RY-coding of first and third mtDNA codon positions does not result in significant topological changes, but increases support for Tetrapus nested within the Agaonidae (Supplementary Fig. S2D,E and Supplementary Table S2).

Unrooted and rooted topologies appeared congruent (Supplementary Fig. S2J), showing that rooting does not alter the ingroup topology. Furthermore, the unrooted topologies from Lopez-Vaamonde et al. (2009) and Cruaud et al. (2010) do not show conflicts with the topology presented here (Supplementary Fig. S7b,c).

Parsimony analysis of the combined data set recovers Tetrapus as sister to the remaining Agaonidae (BP = 64) (Supplementary Fig. S6).

We conclude that neither our study nor previous ones have a strong basis for inferring which group is sister to all other agaonids. Accordingly, the placement of Tetrapus remains unresolved. However, we suggest that the repeated recovery of Tetrapus as sister to all other agaonids in previous studies may be due to LBA to the outgroups and we await further studies.

Ficus Phylogeny

The Ficus phylogenetic trees (Fig. 2 and Supplementary Fig. S9) are globally congruent with previous hypotheses (Herre et al. 1996; Weiblen 2000; Jousselin et al. 2003; Rønsted et al. 2005, 2008; Cruaud et al. 2011a; Xu et al. 2011) (Supplementary Table S5).

Monophyly of the subgenera and infrageneric relationships.— Several moderately to strongly supported clades broadly correspond to currently recognized sections or subsections based on previous molecular phylogenetic studies (see Rønsted et al. 2008) and morphology (sections Pharmacosycea, Oreosycea, Americana, Galoglychia, Adenosperma s.l., Sycomorus s.l., Sycocarpus, Eriosycea, and subsections Malvanthera, Conosycea, Urostigma, Ficus, and Frutescentiae). Only 3 of the 6 Ficus subgenera currently recognized based on morphology (Berg and Corner 2005) are recovered as monophyletic with strong support. These are: Sycomorus (BP = 71, PPMrBayes = 0.75, PPBEAST = 1.00); Sycidium (BP = 100, PPMrBayes = 1.00, PPBEAST = 1.00); and Synoecia (BP = 100, PPMrBayes = 1.00, PPBEAST = 1.00). Relationships of deeper nodes are not strongly supported. The first split within Ficus is between section Pharmacosycea (BP = 100, PP = 1.00, PPBEAST = 1.00) and the remainder of Ficus (BP = 39, PPMrBayes = 0.88, PPBEAST = 0.85). The next split is between a clade with all members of subgenus Urostigma except subsection Urostigma (BP = 100, PPMrBayes = 1.00, PPBEAST = 1.00) and a clade with members of subsection Urostigma, subgenus Sycomorus and all other dioecious figs (BP = 66, PPMrBayes = 0.95, PPBEAST = 1.00).

Exploration of bias in the Ficus phylogenetic trees.—Previous molecular studies are similar in recovering section Pharmacosycea (pollinated by the genus Tetrapus) as sister to the other Ficus species. However, with the exception of the BEAST analysis by Xu et al. (2011), this relationship is supported by parsimony only (Supplementary Table S5). The difference in likelihood scores between our best ML tree and the trees from analyses constrained to place either subgenus Sycomorus or a clade of subgenera (Sycomorus, Sycidium, Ficus, and Synoecia) sister to the remaining Ficus were not significant (Supplementary Table S2). This confirms that relationships within Ficus are unstable along the backbone of the tree and should be regarded as uncertain.

Analyses conducted to assess the impact of missing data on the accuracy of our phylogeny resulted in topologies that were congruent with the topology estimated from the global data set (Supplementary Table S2). It is noteworthy that using only Ficus species for which at least 3 gene regions were available slightly increases node support, but deeper nodes remain unresolved (not shown).

Cophylogenetic Comparisons

All our analyses rejected the null hypothesis of no correlation between fig and wasp phylogenies. Using distance-based methods, the global test of cospeciation (Parafit) rejected a random association between host and pollinator taxa (ParaFitGlobal = 1.37866, P ⩽ 0.01). Further, 176 of the 200 tests of individual host-associate pairs resulted in significant associations between figs and their agaonid pollinators (P ⩽ 0.01) (Supplementary Table S6).

In event-based analyses, exact results depend on the weights assigned to different speciation events. Under the classic TreeMap cost-model of zero for cospeciation and one for other events (Charleston and Page 2002), Jane 2 inferred 198 cospeciation events, 204 duplications, 102 host shifts and 61 losses between fig and wasp phylogenies, accounting for an optimal cost of 367. Whatever the cost model used, the number of cospeciation events inferred by Jane 2 was always significantly greater than expected by chance (Supplementary Fig. S10).

This topological correlation suggests codiversification, but does not establish a time line, so we next used independent relaxed molecular-clock dating techniques to test for contemporaneous divergence (Percy et al. 2004). We found strong temporal congruence between both stem and crown mean ages of most partner clades and between the ages of inferred cospeciation events (Fig. 3). Codiversification test results were not sensitive to the order of deep branching in the phylogenies.

Figure 3

Temporal evidence for fig and fig wasp codivergence. a) Correlation between stem and crown mean ages of major fig and wasp groups (with 95% HPD). b) Temporal congruence of the 198 cospeciation events inferred by Jane 2.

Evolution of Pollination Modes

Parsimony and likelihood reconstruction on the wasp topology both inferred the ancestral pollination mode as ambiguous. Parsimony inferred passive and active pollination as equiprobable ancestral conditions. Similarly, the likelihood difference between the 2 states was not significant (proportional likelihoods of 0.53 and 0.47, respectively) (Supplementary Fig. S11). Using the Ficus topology, parsimony again inferred active and passive pollination as equibrobable. However, likelihood favors active pollination as the ancestral condition (proportional likelihood of 0.91 versus 0.09 for passive pollination). Overall, the reconstructions reveal that pollination modes are homoplastic with several independent shifts between states (passive/active) along both phylogenies (Supplementary Fig. S11).

Biogeographic Analyses

Our dating analyses indicate that the current pantropical distribution of the mutualism cannot have resulted simply from vicariance following the breakup of Gondwanaland. Instead, our ancestral area reconstructions suggest that figs and their pollinators arose simultaneously in Eurasia (Mesquite proportional likelihood = 0.72 for figs and 0.97 for wasps, Supplementary Fig. S12) during the Late Cretaceous ~75 Ma (74.9 Ma for figs and 75.1 Ma for wasps; Fig. 2, Supplementary Fig. S13, and Table 1). Mesquite and Lagrange results were similar, indicating that fig wasps most probably arose in Eurasia. However, Lagrange reconstructions for the fig phylogeny were equivocal due to a basal polytomy. The Eurasian region was proposed as the ancestral area of origin for Ficus in one of the alternative reconstructions that fall within 2 log-likelihood units of the optimal scenario (data not shown, but available upon request). Although the concordance in means crown ages is striking, the PP density around the mean estimate is quite wide (101.9–60.0 for figs and 94.9–56.2 for wasps; Table 1).

Table 1.

Comparison of mean age estimates (Ma) for selected nodes in the fig and wasp phylogenies

Nodes	Ficus mean age Ma (95% HPD)	Agaonidae mean age Ma (95% HPD)
Crown Ficus/Crown Agaonidae	74.9(101.9–60.0)	75.1(94.9–56.2)
Crown Pharmacosycea/Tetrapusa	16.2(25.7–8.2)	15.9(22.0–9.3)
Stem Pharmacosycea/Tetrapusa	74.9(101.9–60.0)	62.1(79.0–45.2)
Crown Sycomorus/Ceratosolen	49.1(67.4–34.0)	59.7(75.4–43.3)
Stem Sycomorus/Ceratosolen	59.4(83.9–43.5)	69.0(87.6–51.5)
Crown sect. Adenospermae/sg. Strepitus	35.8(51.9–23.1)	41.6(57.1–26.7)
Crown sect. Sycocarpus/sg. Rothropus	39.5(54.3–26.5)	43.4(58.2–30.6)
Crown sect. Sycomorus/sg. Ceratosolen	35.1(50.7–23.3)	52.7(67.8–38.8)
Stem sect. Adenospermae/sg. Strepitus	49.1(67.4–34.0)	56.5(72.0–41.5)
Stem sect. Sycocarpus/sg. Rothropus	45.0(62.8–31.9)	56.5(72.0–41.5)
Stem sect. Sycomorus/sg. Ceratosolen	45.0(62.8–31.9)	52.7(67.8–38.8)
Crown Sycidium/Kradibia	38.9(56.6–26.3)	53.6(70.0–39.7)
Stem Sycidium/Kradibia	49.2(69.3–33.0)	69.0(87.6–51.5)
Crown Synoecia/Wiebesia	25.4(37.0–15.3)	33.6(46.4–22.7)
Stem Synoecia/Wiebesia	39.2(56.1–26.7)	40.5(54.7–27.6)
Crown Frutescentiae/Blastophaga	31.8(46.0–19.4)	15.1(22.0–9.3)
Crown Eriosycea/Valisia	22.2(34.8–11.9)	34.8(46.5–24.2)
Stem Eriosycea/Valisia	41.6(59.5–28.4)	55.0(70.0–41.4)
Crown “F. pumila group”/Wiebesiab	20.4(30.9–11.3)	32.4(45.3–20.4)
Stem “F. pumila group”/Wiebesiab	24.4(36.4–14.6)	50.5(63.9–38.1)
Crown Malvanthera/Pleistodontes	28.5(41.7–17.8)	24.1(34.4–14.9)
Stem Malvanthera/Pleistodontes	43.4(61.0–29.4)	50.3(62.0–37.1)
Crown Oreosycea/Dolichoris	31.0(46.4–17.6)	41.0(52.8–30.2)
Stem Oreosycea/Dolichoris	46.2(64.6–31.0)	48.2(61.5–36.5)
Crown Urostigma/Platyscapa	26.6(40.2–14.8)	28.5(37.6–17.9)
Stem Urostigma/Platyscapa	61.2(84.7–43.9)	43.5(55.5–32.7)
Crown Americana/Pegoscapusa	20.5(29.3–13.1)	18.6(23.8–15.0)
Stem Americana/Pegoscapusa	32.3(46.1–22.1)	38.2(48.5–28.8)
Crown Galoglychia/Afrotropical pollinatorsc	28.3(40.3–18.6)	32.9(41.4–24.5)
Stem Galoglychia/Afrotropical pollinatorsc	32.3(46.1–22.1)	35.0(44.0–26.4)
Crown Conosycea/(Deilagaon, Eupristina, Waterstoniella)	36.3(51.3–23.1)	30.8(39.1–22.9)
Stem Conosycea/(Deilagaon, Eupristina, Waterstoniella)	43.4(61.0–29.4)	32.9(41.4–24.5)
Crown Cyathistipulae/Agaon	11.3(18.3–5.8)	21.6(28.8–15.0)
Stem Cyathistipulae/Agaon	20.0(31.4–11.0)	29.8(38.5–22.7)
Crown Caulocarpae/Courtella	17.7(25.8–10.4)	23.8(31.0–16.9)
Stem Caulocarpae/Courtella	24.9(35.3–17.2)	32.9(41.4–24.5)

Note: 95% lower and upper highest posterior distribution inferred by BEAST is reported between parentheses.

aCrown-group Pegoscapus and Tetrapus were assigned uniform prior distributions with minimum ages of 15 Ma and maximum ages of 60 Ma based on Dominican amber fossils.

b“F. pumila group” refers to the clade including F. pumila, F. oleifolia, and F. deltoidea.

c Allotriozoon excepted.

Overall, our analyses favor an Eurasian origin for both Ficus and their pollinators. Indeed, in most Eurasian clades, Sino-Himalayan figs and their associated pollinators appear sister to the rest of the species (Fig. 2, gray rhombus). The overall biogeographical signal was similar across the different methods used and showed instances of dispersal resulting in southward range expansion. The major lineages of figs and pollinators split during the Tertiary and it appears that they then spread southward from Eurasia (Fig. 4), as reflected by the branching order of several clades (Fig. 2, gray arrows). The major lineages subsequently diversified within the Paleotropics and Neotropics during the Miocene.

Figure 4

Hypothetical biogeographical scenario of mutualism diversification. Scenarios are presented on 4 different maps for clarity. Solid arrows: figs; dashed arrows: wasps. Nodes with BP < 70% or PP < 0.95 are collapsed.

Discussion

Codiversification

Our analyses provide both topological and temporal lines of evidence to indicate that figs and fig wasps may represent the first significant case of long-term (~75 myr) codiversification in an insect–plant association. The existence of mutualism per se appears insufficient for codiversification, because speciation in other intimate and sophisticated insect pollination mutualisms (e.g., Yuccas and Yucca moths, Glochidion and Epicephala moths) seems to be driven by host shifting and host tracking rather than cospeciation (Smith et al. 2008; Kawakita and Kato 2009). A plausible explanation for the significant pattern of cospeciation in the fig–fig wasp mutualism is the unusually strong phenotypic coadaptation of key traits, such as the specificity of the chemical mediation between partners (Grison-Pige et al. 2002), the lock-and-key shapes of fig ostioles and wasp heads (van Noort and Compton 1996; Kjellberg et al. 2001).

Despite a history dominated by codiversification, there are also some clear mismatches between fig and wasp phylogenies (Fig. 2). Our analyses support some ancient host-shifts (e.g., by the pollinators of Eriosycea, Conosycea, and F. carica), implying that coadapted pollinators are sometimes replaced by other wasp species without collapse of the mutualism. Finally, several host shifts occur at shallow nodes, such as between Ficus species in the section Americana (Supplementary Fig. S10).

Overall, our tree reconciliation analyses suggest that fig-pollinator history includes numerous species duplications and host shifts, as well as cospeciation events. However, most host shifts are inferred to have occurred between relatively closely related fig species, consistent with observations of extant wasp species occasionally sharing 2 closely related fig species (Molbo et al. 2003; Erasmus et al. 2007). If more distant host shifts were common, the congruence of fig and wasp phylogenies would be eroded rapidly, even if cospeciation remained common (Machado et al. 2005; Cook and Segar 2010). Considering the uncertainty of the sister to all other fig-pollinating wasps, it should be noted that an alternative topology with Tetrapus as sister to all other fig-pollinating wasps would mirror the position of Ficus section Pharmacosycea as sister to all other figs and should therefore increase cophylogenetic signal.

Biological observations and phylogenetic trees show that pollinators of figs are clustered into groups that are consistently associated to Ficus sections, subsections, and even to some species groups of figs. These inter-and intra-generic wasp clusters are highly diverged and relatively old and groups of wasps rarely experience shifts to other groups of figs. Considering only resolved nodes of both phylogenies (Fig. 2, white boxes), we observed only 4 shift events between fig subgenera [Blastophaga psenes, Wiebesia cf callida, 3 pollinators of Frutescentiae (Wiebesia pumilae, W. quadrupes and W. sp. ex F. oleifolia) and Valisia spp.] and 5 shift events between fig sections [Platyscapa bergi and P. sp. (ex F. glaberrima) to Conosycea, Ceratosolen vissali and Ceratosolen sp. (ex F. semivestita) to Sycocarpus and Adenosperma, respectively, and K. subulatae and K. sessilis to Palaeomorphe]. This could be explained by: (i) the allopatry of many fig and agaonid groups, (ii) the differences between habitats of their host figs (e.g., forest canopy versus savannah), and (iii) their host specificity due to intricate coadaptation of phenotypes, including key traits involved in their reproduction. Consequently, we hypothesize that these figs and associated wasp groups evolve largely independently as closed systems (see also Machado et al. 2005; Cook and Segar 2010) and rarely exchange genes or pollinator species, a kind of higher level of lineage sorting.

Recent analyses suggest that the stability of fig/pollinator associations can be erratic before complete lineage sorting has occurred, or before ecological/geographical isolation of fig groups (Machado et al. 2005; Jackson et al. 2008; Jousselin et al. 2008; Renoult et al. 2009; Cornille et al. 2012). During that period of time, pollinator duplication, extinction, and hosts shifts within local groups of related figs sharing similar phenotypic traits may occur frequently. However, these events should not disrupt long-term phylogenetic correlations, if lineages sort over evolutionary time (Cook and Segar 2010). Future work should also seek to understand the patterns and processes of cospeciation and other processes between closely related figs and wasps such as within fig sections. Previous studies at this level have focused on sections Americana (e.g., Machado et al. 2005; Jackson et al. 2008), Galoglychia (e.g., Jousselin et al. 2008), and Sycomorus s.l. (Weiblen and Bush 2002) and future studies should focus on adding also more dioecious and Eurasian clades and to explain why the degree of cospeciation appears to vary between clades.

Wasp and Fig Systematics

Our phylogenetic trees provide several new insights into the systematics of figs and fig wasps and a sound evolutionary framework for future studies in community and behavioral ecology. The question of which groups of wasps and figs are sister to the rest of agaonids and figs, respectively, remains open. Statistical support for the deeper nodes of the phylogeny is low and precludes us from drawing any definite conclusion. Our additional analyses show that there is little support for Tetrapus as sister to all other Agaonidae based on molecular data (see Appendix S3 in the Supplementary Material Online for details). Instead, it appears that Kradibiinae (Ceratosolen + Kradibia) or Group 4 (most Wiebesia species and pollinators of subsection Frustescentiae) are good candidates for the sister taxon to all other agaonids. We raise the possibility that an LBA artifact may have confounded all previous molecular analyses resulting in the inference of Tetrapus as sister to all other agaonids (Appendix S3 in the Supplementary Material Online).

One notable difference between this and previous studies concerns taxon sampling, which can be critical in phylogenetic analyses (e.g., Philippe et al. 2005). By including a larger number of species and using sampling that reflects the known diversity of each group, we may counteract potential problems with long branches (e.g., Bergsten 2005). Increasing taxonomic sampling to break up long branches has been applied repeatedly, often with the conclusion that earlier studies were misled (e.g., Soltis and Soltis 2004; Stefanovié et al. 2004; Leebens-Mack et al. 2005; Phillips et al. 2010; Philippe et al. 2011).

To further investigate Tetrapus placement, we analyzed several morphological characters used in the past to assess the relationships between agaonid genera (Ramirez 1978; Wiebes 1982b; see Appendix S3 in the Supplementary Material Online and Supplementary Fig. S14 for details). We show that there is no evidence from these morphological characters to support Tetrapus as sister to all other agaonids. On the contrary, several independent morphological characters support Tetrapus as nested within the family and closely related to Dolichoris (Appendix S3 in the Supplementary Material Online).

Further studies using more genes and increased taxonomic sampling of both ingroups and outgroups of figs and wasps are still needed and should contribute to resolving the higher taxonomic group relationships with more confidence and determine the earliest divergence among the figs and their pollinating wasps.

Pollination

A number of authors (including some coauthors of this study) have previously proposed passive pollination as the ancestral mode for the mutualism, followed by a single shift to active pollination and several independent reversions to passive (Machado et al. 2001; Jousselin et al. 2003; Herre et al. 2008; Jandér and Herre 2010). This hypothesis has intutive appeal as most other insects that pollinate do so passively, but it based primarily on the fact that Tetrapus wasps are passive pollinators and appeared as the sister of all other pollinators in previous phylogenetic analyses. Similarly, their host figs (Pharmacosycea) appeared as sister to all other figs. However, our new phylogenetic trees support a different phylogenetic position for Tetrapus.

Consequently, the issue of ancestral pollination mode must be revisited. The phylogenetic tree itself is in question, but we also highlight a key issue about the interpretation of a given phylogeny. The previous conclusion that passive pollination is ancestral relies on the assumption that “basal” branches of the trees are more informative about ancestral character states. However, there is no reason to assume that traits found in Tetrapus/Pharmacosycea are “more primitive” or represent traits of the common ancestors of both sister groups (Krell and Cranston 2004; Crisp and Cook 2005; Lamm and Redelings 2009). At the present time, the ancestral pollination mode should be considered as equivocal and our analyses imply that it remains so.

Of our 4 reconstructions, 3 find the ancestral state equivocal, whereas one (ML on fig phylogeny) favors active pollination. This indicates that further studies are needed to infer the ancestral pollination mode with more confidence. Importantly, these results were established on fully bifurcating trees, but in reality the backbones of both trees are not strongly supported and may change in future studies. Recent advances in our understanding of the morphological evolution of Moraceae, and in particular of an expanded tribe Castilleae, the figs closest relatives, may also shed new light on the evolution of the mutualism (Clement and Weiblen 2009).

Cobiogeography

Molecular divergence time estimates point to a Cretaceous origin for the mutualism, but differ with respect to biogeographic scenarios (Supplementary Tables S3 and S5). Previous biogeographic analyses of fig-wasps have argued in favor of Gondwanan vicariance (Machado et al. 2001). However, a previous study by Lopez-Vaamonde et al. (2009) reconstructed ancestral areas of fig-pollinating wasps using a phylogenetic tree with Tetrapus as sister to the remainder of the fig-pollinating wasps. These authors concluded that the MRCA of all extant fig-pollinating wasps was most likely Asian, although a southern Gondwanan origin could not be rejected. A Laurasian origin with subsequent dispersal has been proposed for figs and their nearest relatives (Zerega et al. 2005).

Our analyses indicate that the fig-wasp mutualism was already in existence ~75 Ma in Eurasia and our independently derived mean date estimates for figs (75.1 ± 19.4 Ma) and wasps (74.9 ± 21.0 Ma) crown groups are remarkably similar, although the size of the confidence intervals introduces a degree of uncertainty. Despite differences in sampling and dating algorithms, the dates obtained correspond well with most previous estimates (Supplementary Tables S3 and S5).

In addition, the hypothesis of an Eurasian origin of the mutualism is supported by several other lines of evidence: (i) the presence in Asia of 70% of the major Ficus clades; (ii) the early divergence of Sino-Himalayan fig and wasp species in most Eurasian clades (e.g., F. henryi, F. sarmentosa, F. tikoua, F. nervosa); (iii) the fact that pollinators of the subsection Frustescentiae are found only in Continental Asia (Fig. 2); (iv) the age estimates for Moraceae in general and Ficeae (Dorstenieae and Castilleae) in particular (Zerega et al. 2005); and (v) the fact that the oldest fig and wasp fossils are known only from the Northern Hemisphere (Collinson 1989; Compton et al. 2010) (see our review of the literature on fig fossils in Appendix S4 in the Supplementary Material Online). Finally, Burnham and Graham (1999), analyzing the origin of the tropical component in northern Latin American vegetation also suggest that Ficus arrived from the north. Accordingly, current data support the conclusion that the mutualism probably originated in the tropical forests of Eurasia (Otto-Bliesner and Upchurch 1997).

Pharmacosycea and Tetrapus divergence is dated to 74.9–62.1 Ma (mean stem figs-mean stem age wasp; Table 1), before South America split from Antarctica. However, rather than explaining the South American colonization of Pharmacosycea/Tetrapus by trans-antarctic routes, we propose that both lineages might have reached the New World across North Atlantic land bridges (Tiffney 1985), dispersing through the evergreen woodland and tropical forest belts of Eurasia (Fig. 4). South America may have been colonized later via “stepping-stone” volcanic islands. Indeed, most Pharmacosycea species inhabit the Northern Andes and there are none in Chile and Patagonia, which have vegetation similar to late Cretaceous Antarctica (Poole et al. 2003). Because figs are also absent from the exceptionally good fossil flora of Patagonia (Wilf et al. 2003), trans-Antarctic dispersal seems unlikely but cannot be completely ruled out. Although the Laguna del Hunco flora is not considered a tropical flora, this Patagonian flora hosts one Papuacedrus species very closely related to extent tropical Papuan species (Wilf et al. 2009). In Papua New Guinea and in Papua Barat (Indonesia), this conifer is found in the same mountainous habitats as Malvanthera fig tree species at altitudes >2000m (for example in the Arfak mountains, Kebar Valley, Bulolo-Wau, Mt Kerewa). Accordingly, despite not being considered a tropical flora, the Laguna del Hunco flora could have also have included Ficus and the fact that Ficus appears absent from this exceptionally good flora, supports the later arrival of Ficus from Eurasia.

Based on our biogeographic analyses, the major lineages of figs and pollinators split during the Tertiary and spread southward from Eurasia (Fig. 4), possibly in response to the cooling climate (Davis et al. 2002). Subsequent diversification occurred within continents during the warmer Miocene epoch (Zachos et al. 2001).

The general scenario of fig-wasp codiversification is illustrated by the charismatic hemi-epiphytic or “strangler” figs (the subgenus Urostigma clade, Fig. 4), which evolved ~52–50.3 Ma, during a period of global warming. A first clade, probably living west of the Turgai straits (Akhmetiev and Beniamovski 2009), dispersed southward into Africa to form section Galoglychia and into South America to form section Americana, some 32.3–38.2 Ma. Another clade, probably occurring in east Eurasia, spread to India and Sundaland to form section Conosycea and to Australasia to form section Malvanthera, ~50.3–43.4 Ma. This latter dispersal was probably via stepping-stones through the Ninety East Ridge (Carpenter et al. 2010), because direct dispersal from Sundaland to Australia was impossible before 25 Ma (Hall 2002). Today, each tropical continent has its own major endemic radiation of strangler figs, stemming from these ancient dispersal processes. Interestingly, pollinator biogeography shows a few discrepancies with this scenario for fig dispersal. Indeed, the genus Pleistodontes pollinating Malvanthera figs is sister to all other Urostigma pollinators. Therefore, we propose that Conosycea was colonized by a host shift of an ancestral Galoglychia/Americana pollinator in southern Eurasia before spreading to southern Sunda.

Conclusion

Based on multiple lines of evidence (fossils, Moraceae history, branching pattern, and ancestral area reconstructions), we infer an Eurasian origin for the fig/pollinator mutualism. We show that the mutualism arose ~75 Ma, confirming previous estimates (Supplementary Tables S3 and S5). Because that time, the insects and plants have diversified together leaving a strong long-term signal of phylogenetic congruence, confirming previous studies based on smaller data sets (Supplementary Tables S3 and S5). This is not due to strict cospeciation alone, but reflects a history with large amounts of cospeciation and insufficient host shifts to alter the marked phylogenetic matching. This is the only known example of long-term insect/plant codiversification and we are not aware of other candidates for such a pattern. Other insect/plant pollination mutualisms do not appear to be characterized by phylogenetic congruence, and we propose that strong codiversification of figs and their pollinators is driven by their unusually high level of phenotypic trait matching. Figs and their pollinators have spread across the globe to occupy all tropical continents, where they play important ecological roles in forests and savannahs. Their numerous interactions with other species, such as vertebrate frugivores, mean that the evolution of entire tropical ecosystems has been influenced strongly by this unique strong pattern of codiversification between fig trees and their pollinating insects.

Supplementary Material

Supplementary material, including data files and/or online-only appendices, can be found in the Dryad data repository at http://datadryad.org, doi: 10.5061/dryad.hr620. Matrices are also available in TreeBASE (No. TB2:S13315) at http://www.treebase.org/.

Funding

This work was supported by the ANR projects “NiceFigs” and “BioFigs/National Science Council, Taiwan, R. O. C. code: 98WFA0100291” led by M.H.M. and L.S.C.; by a grant from the Carlsberg Foundation and a Marie Curie grant “CoEvol” to N.R.; grants from the Leverhulme Trust, NERC, the Royal Society, and the European Research Council to V.S.; NSC [project 99-2923-B-002-001-MY3 to L.S.C]; Fapesp [grant 04/10229-4 to R.A.S.P.]; NERC to J.M.C.; NSF [grant 0753306 Harvey Mudd College and R.L.H., B.C., J.P., and T.S] and the Howard Hughes Medical Institute [grant 52006301 to Harvey Mudd College and R.L.H., B.C., J.P., and T.S.]. Funding to pay the Open Access publication charges for this article was provided by The Natural History Museum of Denmark.

Author Contributions

J.Y.R. and V.S. designed the research. As.C., Bh.C., W.C., J.M.C., R.D.H., E.J., C.K., F.K., M.H.M., Sv.N., R.A.S.P., P.Y.Q., J.Y.R., R.U., C.L.V., G.D.W., and D.R.Y. provided material or data. Sv.N., R.U., and J.Y.R. identified the fig wasps. F.K., J.Y.R., G.D.W., and N.R. identified the figs. G.G., R.J.Z., and A.C. performed and coordinated fig wasp DNA sequencing. N.R. and W.C. performed and coordinated fig DNA sequencing. As.C., N.R., and J.Y.R. performed the analyses. J.Y.R. examined fig wasp morphological characters and reviewed the literature about fig fossils. R.L.H., A.Y., T.S., J.P., and B.C. implemented statistical tests and performed analyses with Jane2. As.C., N.R., J.Y.R., J.M.C., G.D.W., and V.S. wrote the article, with major comments from R.H., C.K., and F.K. All authors commented on the article.