Multilocus Species Trees Show the Recent Adaptive Radiation of the Mimetic Heliconius Butterflies

1Visual mimicry provides an excellent system in which to study the origins of biodiversity, as the targets of selection are clearly identifiable and the role of natural selection in promoting adaptation and ultimately speciation can be directly observed. Studies of mimetic assemblages have been instrumental in explaining many biological phenomena, but testing hypotheses regarding the evolution of mimicry depends heavily on our knowledge of systematic relationships between the participating taxa, particularly where mimics are closely related. Unfortunately, both strong selection on adaptive loci, which may facilitate adaptive introgression, and rapid radiation leading to incomplete lineage sorting (ILS), are likely to be common in many systems, but are especially prevalent in mimetic butterflies. These processes can significantly interfere with the estimation of phylogenies.

2A large body of recent work has been devoted to the issue of incongruence between the species tree (the true speciation history) and the gene trees evolving within. The traditional approach of concatenating the total genetic evidence into a supermatrix to obtain a global estimate of the predominant phylogenetic signal and hidden support, without consideration for the heterogeneity of individual partitions, has been to some extent superseded by multispecies coalescent (MSC) techniques (reviewed in:;;;; ). The majority of the new MSC algorithms are intended to model at least some of the sources of heterogeneity between different markers, most frequently focusing on the problem of ILS (e.g.; ), sometimes addressing hybridization (e.g.,; ), and in at least one case modeling discordance without specifying its potential source. Although the supermatrix approach remains popular and serves as an effective approximation of the species diversification history in most cases, its ability to properly assess the degree of statistical support for phylogenies has been brought into question and contrasted with the potential of the MSC techniques to assess confidence more realistically. Heliconiini are an especially interesting subject for a systematic study, where the purported robustness of MSC tools to gene flow and other sources of incongruence can be tested with real data.

3We estimate the phylogeny and investigate the link between the dynamics of speciation and macroevolutionary factors in the Neotropical butterfly tribe Heliconiini (Nymphalidae: Heliconiinae). Heliconiini are arguably the most thoroughly researched example of microevolution in the Neotropics, the most biologically diverse region of the world. They comprise the genus Heliconius and nine smaller genera, providing a spectacular example of a radiation where speciation is promoted by divergence in mimicry of aposematic wing patterns. Heliconiini, and especially Heliconius, form Müllerian mimicry rings of distantly related toxic species, in which mimetic species share the cost of educating avian predators by evolving similar wing patterns. Thus, Heliconiini are an excellent system for the study of convergence from both genomic and organismal perspectives (e.g.,;;;;;;;;; ). In addition to studies of wing pattern evolution, there has been a recent proliferation of comparative molecular and genomic studies of other traits including vision, chemosensation and cyanogenesis. A robust and stable molecular phylogeny is therefore especially desirable for this clade.

4Heliconius, the most species-rich genus of Heliconiini, has the key features of an adaptive radiation. All Heliconius are found in broadly similar Neotropical habitats and most of the species appear to originate on the eastern slope of the Central and North Andes. Heliconius is far more diverse than the other nine genera of the tribe, and individual species show high disparity in their dietary adaptations, and wing pattern diversity. Numerous field and laboratory experiments have demonstrated the adaptive value of wing patterns in protection from predators. Importantly, divergence in wing pattern also leads to both pre- and post-mating reproductive isolation and hence contributes to speciation. There is also evidence that ecological divergence in mating behavior, preference for microhabitat, and host plant use permit sympatric co-occurrence of closely related species. There is also evidence for putative key innovations, most notably pollen feeding, which enables unusually long lifespans and underlies much of the unusual ecology of the group. An expansion of olfactory receptor gene families plays a likely role in specific host and food plant preferences. In summary, Heliconius have undergone a burst of diversification associated with adaptive changes, some of which directly cause reproductive isolation. Nonetheless, one characteristic of adaptive radiations that has not yet been demonstrated in Heliconius is a temporal burst of species diversification.

5An unusual feature of Heliconius is the prevalence and importance of gene flow and hybridization, leading to a controversy over the validity of traditional species concepts in the clade. At least 26% of all species of Heliconiini occasionally produce interspecific hybrids in the wild, and at least one species, Heliconius heurippa, has resulted from homoploid hybrid speciation from parental forms diverged millions of years ago. Genome sequencing has shown that mimetic diversity in the Heliconius melpomene and silvaniform (Heliconius numata) clades is also explained by adaptive introgression of genes regulating the aposematic wing patterns. In addition, neutral gene flow seems to be widespread, influencing as much as 40% of the genome between H. melpomene and Heliconius cydno.

6Heliconiini systematics has a long history, starting with early morphological work, through allozymes and a combination of morphological and ribosomal DNA-restriction data, to studies based on the sequences of mitochondrial and nuclear markers. The most comprehensive study to date by attempted to address some of the difficulties by incorporating many taxa (38 of 46 Heliconius, 59 of 77 Heliconiini), considering two individuals of most species, and sequencing two mitochondrial (CoI/II and 16S) and four nuclear markers (EF1α, Wg, Ap, and Dpp). However, this data set is still potentially inadequate to address the challenges posed by Heliconiini systematics, as 3 of the loci (16S, Ap, and Dpp) were only sequenced for 12 representative species. Among the other three markers, 65% of the variable sites resided in the fast-evolving mitochondrial partition CoI/II, raising the possibility that the inferred relations are largely driven by the historical signal of the matriline. The relationships between Heliconius, Eueides and the other eight genera were not resolved with good support, as might be expected if most of the phylogenetically informative variation comes from a fast-evolving partition. Importantly, the individual data partitions were analyzed in concatenation. Despite these shortcomings, Beltrán and colleagues confirmed that the morphological and behavioral characteristics of the major clades do not correspond directly to their evolutionary history and that the traditionally recognized genera Laparus and Neruda are most likely nested within the crown genus Heliconius.

7The importance of Müllerian mimicry as a driver of individual speciation events has been well-established whereas the macroevolutionary processes governing the evolution of the group have been largely neglected in empirical studies (but see; ). Precise understanding of the evolution of the mimicry rings, as well as associated processes such as hybridization, requires knowledge of the relative timing of the divergence events and motivated the most widely cited study of the molecular clock in Arthropoda. created a chronogram from a partially unresolved tree, using a relaxed clock procedure and the CoI/II alignment. Importantly, the first dated phylogeny of Heliconiini is not calibrated to an absolute standard, making it impossible to make inferences about the relation of the diversification process and the contemporaneous geological and climatic events. A recent comparative study suggests that the majority of Heliconius lineages originated in the northeastern Andes and spread to other parts of the continent. The cumulative results of over 200 systematic studies demonstrate that most South American tropical clades have experienced periods of significantly elevated net diversification rate in response to Andean orogenesis, alterations in the hydrology and sediment dynamics of the present-day Amazon Basin, as well as local and global climatic changes. These processes can result in allopatry of incipient lineages, or in creation of new ecological niches and change to the species-level carrying capacity of the environment leading to ecological speciation. In particular, we can hypothesize that the diversification rate of Heliconiini increased during the periods of especially rapid Andean uplift around 23, 12 and 4.5 Ma.

8Here we aim to resolve the species tree of Heliconiini radiation and generate a data set including nearly all of the currently valid species in the tribe, sampling intraspecific diversity across the range of many species, combining whole-genome sequencing with Sanger sequencing. We apply a wide range of phylogenetic methods to reconstruct the species tree, including supermatrix, coalescent and network approaches, which allow us to assess the strength of the underlying signal of speciation. The power of our combined approach is harnessed to elucidate the importance of marker heterogeneity for the final assessment of systematic relationships, while realistically estimating the support values for our chosen topology. We date the time of individual divergence events with confidence, permitting an analysis of diversification dynamics. We thus present a comprehensive study of macroevolutionary dynamics in a mimetic system that has been studied intensively at the microevolutionary level.

9We sampled 180 individuals, including 71 of the 77 species in all genera of Heliconiini and 11 outgroup species. The specimens came primarily from our collection at the University of Cambridge, with additional specimens shared by museums and private collectors (online Appendix S1, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). We included five outgroup species from the sister tribe Acraeini and three from the related genus Cethosia. The diverse analyses used in this article require different sampling designs and the demands of all the techniques cannot be easily accommodated in a single data set. For example, the network analysis based on nucleotide distance produced much better supported and resolved trees when the 95% or more incomplete data from historical specimens were not used, whereas the various MSC techniques required the use of at least two individuals per species and had to be based only on taxa with intraspecific sampling. Thus, we distinguish four data sets. The complete data matrix includes all the data. The core data set excludes 14 individuals represented solely by short DNA fragments from historical specimens. The single-individual data set includes both modern and historical specimens, but with only the single best-sequenced individual per taxon. Finally, the *BEAST matrix contains only the 17 species of Eueides and Heliconius with extensive sampling of multiple representatives of each species.

10We used 20 nuclear and 2 mitochondrial loci as markers (Supplementary Tables S1 and S2; online Appendix S1, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). The selection includes the three classic molecular markers for Lepidoptera (CoI/II, EF1α, and Wg), two markers proposed by (16S and Dpp), eight new universal markers proposed by (ArgK, Cad, Cmdh, Ddc, Idh, Gapdh, Rps2, and Rps5), and nine highly variable loci identified by (Aact, Cat, GlyRS, Hcl, Hsp40, Lm, Tada3, Trh, and Vas). Additional Heliconius-specific primers were designed for Cmdh, Gapdh, and Idh. Details of the primers and PCR cycles are listed in Supplementary Table S1, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j. For most species, sequences of the three basic markers for multiple individuals were already published, and data for 26 individuals came exclusively from GenBank (online Appendix S1, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j).

11We generated Sanger sequences for 103 specimens (online Appendix S1, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). DNA was isolated from approximately 50 μg of thorax tissue using the DNeasy Blood & Tissue kit (Qiagen, Manchester, UK). PCR was carried out in a total volume of 20 μl, containing 1× Qiagen Taq buffer (Manchester, UK), 2.5 mM MgCl2, 0.5 μM of each primer, 0.2 mM dNTPs, 1 unit bovine serum albumin, 0.5 unit Qiagen Taq-Polymerase and 1 μl of the DNA extract. The following program was executed on a G-Storm cycler (Somerton, UK): denaturation 5 min at 94°C; 35 cycles of 30 s at 94°C, 30 s at the annealing temperature and 90 s at 72°C; final extension for 10 min at 72°C. The results were visualized by electrophoresis in 1.5% agarose gel stained with 1% ethidium bromide. PCR products were cleaned using the ExoSAP-IT system (USB, Cleveland, OH): 60 min at 37°C; 20 min at 80°C. We used gel purification with the Nucleo Spin Extract II kit (Macherey-Nagel, Düren, Germany) as needed. Sanger sequencing reaction was carried out with the BigDye Terminator v. 3.1 (AB, Foster City, CA): 2 min at 94°C; 25 cycles of 10 s at 94°C, 5 s at 50°C and 4 min at 60°C. The products were sequenced with the ABI 3730xl DNA Analyzer at the Sequencing Facility, Department of Biochemistry, University of Cambridge. We manually inspected the traces in CodonCode v. 4.0.4 using PHRED for quality assessment.

12At the time of our Sanger sequencing effort, whole-genome data generated for other studies became available from 57 individuals in 27 common species (;;; Dasmahapatra K.K., Mallet J., unpublished data) (online Appendix S1, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). 100 bp reads were generated using the Illumina Genome Analyzer II platform with insert size of 300–400 bp. We performed de novo assembly of the short reads in the program Abyss v. 1.3. Based on previous studies and preliminary results (Baxter S., personal communication), we chose k-mer length of 31, minimum number of pairs and minimum mean coverage as optimal settings. The 20 nuclear markers were mined from the assemblies by megaBLAST in Geneious v. 5.5.1 using reference sequences from the H. melpomene genome. The quality of the recovered sequences was assessed by alignment to previously generated amplicon sequences of the same loci from the same individuals.

13Mitochondrial sequences could not be recovered by de novo methods, presumably because the large number of reads from the highly abundant mitochondrial DNA contained a large enough number of erroneous sequences to interfere with the assembly. We reconstructed whole mitochondrial genomes of 27 species by mapping to the H. melpomene reference (The ), using the default settings in the Genomics Workbench v. 5.5.1. These data were analyzed separately from the 21 locus mixed nuclear-mitochondrial alignment.

14Short fragments of CoI/II and EF1α were sequenced from historical specimens up to 150 years of age, obtained from museum and private collections (online Appendix S1, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j) and processed in a vertebrate genetics laboratory to reduce the risk of contamination. Instruments and surfaces were cleaned with 5% bleach and irradiated with UV for 30 min prior to use. One to two legs were washed in water, immersed in liquid nitrogen in a test tube for 30 s and ground up, followed by an extraction into 20 μl of buffer using the QIAmp DNA Micro Kit (Qiagen, Manchester, UK). We treated every fifth extraction and every fifth PCR as a negative control with no tissue or DNA extract. PCRs were carried out in a 20-μl volume using 1 unit of Platinum HiFi Taq Polymerase (Invitrogen, London, UK) and 1× buffer, 2.5 mM MgCl2, 0.5 μM of each primer, 0.2 mM dNTPs, 1 unit bovine serum albumin, sterilized DNAse-free water and 1–5 μl of the DNA extract depending on concentration. To accommodate shearing of DNA with time, we designed and applied PCR primers spanning short fragments of 200–300 bp (online Appendix S2, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). We carried out amplification, product clean up and sequencing as above, partially accounting for possible cross-contamination by blasting the results against GenBank.

15Alignments for each locus were generated in CodonCode v. 3 to account for inverted and complemented sequences, and improved using MUSCLE v. 3.8. We visualized the alignments of the coding loci (all except the mitochondrial 16S and the tRNA-Leu fragment in the CoI/II sequence) in Mesquite v. 2.75 and checked translated sequences for stop codons indicating errors. The whole mitochondrial sequences were aligned to the Acraea issoria and H. melpomene references using the G-INS-i algorithm in MAFFT. The number of variable and parsimony informative sites was estimated for each locus in PAUP* v. 4. Models of sequence evolution implemented in MrBayes were selected in MrModelTest v. 2.3 based on the AIC. Xia's test in DAMBE v. 4.0 demonstrated saturation in the third codon position of CoI/II, prompting us to treat the third codon position of the fast-evolving CoI/II locus as a separate partition when estimating the gene tree. The Leucine tRNA (tRNA-Leu) fragment occurring in the middle of CoI/II displays very low variability and thus was included in one partition with the slower evolving first and second codon positions. Individual gene trees were estimated in MrBayes v. 3.1, using four runs of one chain, 10 million Markov Chain Monte Carlo (MCMC) cycles sampled every 1000, and 2.5 million cycles discarded as burnin based on the average standard deviation of split frequencies becoming less than 0.01 and a plateau in the log-likelihood values. The mitochondrial genes were concatenated due to their shared history, but treated as separate partitions with distinct models. All gene trees were visualized with FigTree v. 1.4.

16We investigated the cyto-nuclear discordance and other conflicts in the phylogenetic signal with several methods. To illustrate the global reticulate signal in the data, a NeighborNet network was built with the pairwise distances calculated under the F84 correction, the most complex model that could be fitted to the data in the program SplitsTree v. 4. We reduced the data set to a single best-sequenced individual per species to exclude the reticulations resulting from the expected recombination within species. Next, the topological disparity among individual loci was illustrated using Multi-Dimensional Scaling of pairwise Robinson–Foulds distance between the gene trees, as estimated by TreeSetViz v. 1.0 in Mesquite. Calculation of the RF required trimming the trees to the minimal set of 54 shared taxa from 27 species, using the R package APE. Finally, we investigated whether topologies and branch lengths of the individual loci are consistent enough to justify concatenation of the markers by means of a hierarchical likelihood-ratio test in Concaterpillar v. 1.5.

17We created a supermatrix of the 20 nuclear and 2 mitochondrial markers in Mesquite. To increase the efficiency of tree searches, optimal partitioning schemes for the complete and single-individual data sets were identified in PartitionFinder v. 1.1, using the Bayesian Information Criterion and a greedy search, followed by selection of nucleotide substitution models. Lists of optimal partitions for ML and Bayesian analyses can be found in online Appendix S2, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j. The maximum-likelihood (ML) phylogeny was searched for under the GTRGAMMA model in RAxML v. 8.1 with 350 bootstrap replicates under the GTRCAT approximation, where the number of repetitions was determined by the bootstopping criterion. To explicitly test the likelihood of various hypotheses for Heliconiini phylogeny, several alternative topologies were created in Mesquite, representing previously identified groupings, as well as possible placements of the enigmatic genera Cethosia, Laparus and Neruda (Supplementary Table S4, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j) (Brower 1994;;;;; ). We then re-estimated the ML tree using each topology as a constraint. The likelihood scores of the original and alternative trees were compared using the Shimodaira–Hasegawa (SH) test and the expected likelihood weights (ELW) based on 1000 bootstrap replicates. We note that due to computational constraints, the assumptions of the SH test were violated by comparing only a subset of all the possible topologies. A separate phylogeny was generated for the unpartitioned whole mitochondrial alignment in RAxML under the GTRGAMMA model with 1000 bootstrap replicates.

18A Bayesian chronogram was estimated using the program BEAST v. 1.8. To avoid incorrect estimates of the substitution rate parameters resulting from the inclusion of multiple samples per species, this analysis was based on a pruned alignment with one individual per species. The only exception is the inclusion of three races of H. melpomene and two races of Heliconius erato, where deep geographical divergences are found. The Vagrantini sequences were not used, as BEAST can estimate the placement of the root without an outgroup. Thus the analysis included 77 taxa and the optimal partitioning scheme was re-estimated appropriately (Table 1). We linked the topology, but modeled an uncorrelated lognormal clock and the substitution rate separately for each partition. Substitution rates were drawn from the overdispersed gamma distribution prior with shape parameter, scale parameter theta = 1 and starting value 0.001 for nuclear genes, and for the faster-evolving mitochondrial loci. We used a Birth–Death tree prior and empirical base frequencies to limit the computation time for the heavily parameterized model.

19As no fossils of Heliconiini or closely related tribes are known, we used secondary calibration points from the dated phylogeny of Nymphalidae. Prompted by the findings of, who demonstrate the potential for biases when secondary dating information is used, we compared having a single calibration point at the root with using all of the eight known split times. We also considered the impact of modeling each prior divergence time as normally distributed with the mean found by and the standard deviation matching the 95% credible intervals, or as uniformly distributed within the same bounds. For each of the four model combinations, four independent instances of the MCMC chain were run for 100 million generations each, sampling the posterior every 10,000 generations and discarding 10 million generations as burnin after convergence assessment in Tracer v. 1.6. Based on the marginal likelihood estimated under the AICM criterion, we chose the scheme with multiple calibrations modeled as uniformly distributed (Supplementary Table S3, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j), although we found that the ages of most of the nodes do not differ by more than 10% between the four schemes. While we recognize that the models should ideally be chosen based on the Path Sampling or Stepping Stone procedures, in practice we found these to be too computationally demanding for our relatively large data set. The input.xml file generated in BEAUti v. 1.8 (Drummond et al. 2010) can be found in the online Appendix S3, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j. To ensure that the results are driven by the data and not the priors, we executed an empty prior run. The Maximum Clade Credibility tree with mean age of the nodes was generated using LogCombiner v. 1.8 and TreeAnalyser v. 1.8. Parallel to the SH and ELW tests, plausible alternatives to the disputed and poorly supported nodes in the Bayesian chronogram (Supplementary Table S5, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j) were tested by calculating the posterior model ratios, a simple and robust alternative to the computationally expensive Bayes factor calculations. A posterior sample of 1000 phylogenies was filtered in PAUP* using constraint trees designed in Mesquite and the model odds were calculated as the frequency of the alternative grouping divided by the frequency of the clade observed in the MCC chronogram.

20To account for the heterogeneous phylogenetic signal resulting from gene flow, hybridization and ILS, we applied a variety of MSC analyses that take as input both the raw alignment and individual gene trees. We first used the established method of Minimizing Deep Coalescences (MDC), taking advantage of a dynamic programming implementation in the package PhyloNet. One hundred samples of 100 trees were drawn randomly without replacement from the distribution of Bayesian gene trees for the 21 loci, an MDC phylogeny was estimated for each sample and a 50% majority rule consensus was taken.

21Bayesian concordance analysis (BCA) is an MSC method that attempts to reconcile the genealogies of individual loci based on posterior distributions, regardless of the sources of conflict. BCA generates concordance factors (CFs), which show what proportion of loci contain a particular clade, and estimates the primary phylogenetic hypothesis from the best-supported clades. CFs offer a powerful alternative to traditional measures of support and can be conveniently estimated in the program BUCKy. We executed two runs of 1 million MCMC cycles in BUCKy based on the 21 posterior distributions of gene trees from MrBayes.

22Another approach to the MSC is to estimate the gene trees and the species tree simultaneously, explicitly modeling the sources of incongruence. We applied this technique using *BEAST, a program harnessing the power of BEAST and simultaneously implementing a powerful MSC algorithm that estimates the species tree and the embedded gene trees, as well as the population sizes of the lineages. Effective calculation of the population size parameters requires a thorough multilocus sampling of each species in the analysis, which forced us to reduce the data set to species with a minimum of three individuals from at least two distinct populations (the *BEAST data set). The final alignment included 87 terminal taxa in 17 species of Eueides and Heliconius. We re-estimated the individual substitution models for each partition (Supplementary Table S2, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j), used a constant population size coalescent tree model, and implemented other priors as described above for BEAST. We carried out four independent runs of 500 million cycles each, sampling every 10,000 cycles, generated a maximum clade credibility species tree and visually summarized the 21 gene trees by plotting in DensiTree. The online Appendix S4, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j, contains the.xml file for this analysis.

23New species of Heliconiini continue to be described based on morphological, genetic, or karyotypic evidence, yet the validity of the taxonomic status of some of these new lineages has been disputed (Mallet J., Brower A., personal communication). To explore the potential for phylogenetic determination of the number of species in the tribe Heliconiini, species were delimited with the novel Bayesian Poisson Tree Process (bPTP) algorithm. Using the ML tree with multiple individuals per species, two runs of 500,000 generations were executed with the first 10% discarded as burnin based on the convergence of log-likelihood values. The results were often poorly supported and included clusters and splits not compatible with the extensive knowledge of the biology of this group (Supplementary Fig. S7, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j), prompting us instead to use the accepted taxonomy of in our study of diversification dynamics in this clade.

24The formal analyses of diversification dynamics were based on the output of the Bayesian supermatrix analysis in BEAST. Initially we investigated the changes in the number of lineages by generating semi-log lineage through time (LTT) plots for Heliconiini and Heliconius in the R package paleotree based on a posterior sample of 1000 chronograms, plotting the 95% credible interval of the divergence time around the median curve. We also carried out preliminary analyses in the R package DDD, using an MCC tree from an analysis with a single root calibration. Comparisons of all possible models, including diversity-dependence and changes in rate regimes, suggested that the diversification rate of Heliconiini increased around 12 Ma, shortly before the appearance of Heliconius, which in turn started to diversify faster around 4.5 Ma.

25To achieve a more nuanced understanding of the diversification dynamic in this clade, we chose the reversible-jump MCMC (rjMCMC) approach implemented in the R package BAMM, which models both the variation in rates between lineages and the change of rates in time. BAMM is distinct from previously published approaches, as it also treats the location of diversification regimes on the tree as a variable to model. In effect, the algorithm does not give false confidence in single, point estimates of rate shift times. To account for the uncertainty surrounding the influence of dating priors and the divergence times, we carried out analyses with four chronograms: the MCC phylogeny, the two trees corresponding to the upper and lower boundary of the 95% credible interval, as well as the MCC tree from an analysis with a single, normally distributed calibration at the root. All the phylogenies were truncated to 1 Ma before the present, to account for the uncertainty in the estimates of recent speciation (protracted speciation) and extinction rates (pull of the present; ). Appropriate speciation and extinction prior values were chosen for each tree in the BAMMtools suite. We experimented with three alternative values of the Poisson process rate prior suitable for relatively small trees (1.0, 2.0, 4.0) and discovered that the results are very similar, as evidenced by insignificant differences in Bayes factors.

26Following the failure of automated species delimitation, incomplete sampling of the taxa was accounted for by coding in the sampling proportion of each genus and each major Heliconius clade according to the taxonomy of, with addition of new species of Philaethria and Neruda, and excluding Heliconius tristero (Mérot et al. 2013) (species sampling fraction = 0.904). A test with the most extreme of the plausible taxonomies (not recognizing H. heurippa, Eueides emsleyi, the undescribed Agraulis or the three recently proposed Philaethria species; sampling fraction = 0.956) demonstrated that the results are robust to the assumptions on the proportion of missing taxa. We assumed that rate changes may occur in increments of 0.5 Myr. For each of the four trees, two independent runs of four Metropolis-Coupled MCMC chains were executed for 2,200,000 generations and the first 200,000 generations were discarded as burnin. Convergence was assessed with the R package coda by plotting the log likelihood of all chains and calculating the effective sample size for log likelihood and the number of rate shifts. BAMMtools were used to analyze the posterior set of 2000 models; compare the frequently observed rate regimes with Bayes factors; identify the branches where rate shifts occurred; calculate the mean speciation, extinction and diversification rates for major clades; and plot the change of these parameters in time.

27We successfully combined three approaches to sequencing, which resulted in the high taxonomic coverage and intraspecific sampling necessary for the MSC methods, and a sufficient sampling of loci for each individual (online Appendix S1, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). Most of the data set consists of Sanger sequences from 108 individuals, 26 of which were already included in GenBank. We also obtained two classic lepidopteran markers CoI/II and EF1α from respectively 13 and 11 out of 14 historical specimens, thus adding eight species of Heliconiini that have not been sequenced previously.

28We capitalized on the availability of Illumina data by generating de novo assembly contigs for an additional 57 individuals. The N50 of the Abyss assemblies ranged from 552 to 1921 bp (average 1206) and all the nuclear markers were successfully recovered by megaBLAST from every assembly. Whole mitochondrial sequences of the same individuals were recovered by read mapping, with about a 400-bp stretch of the hypervariable control region incomplete in some sequences. We obtained a depth of coverage over 100× and high confidence in the base calls due to the high copy number of mtDNA in the tissue. Finally, we included the sequences extracted from the H. melpomene genome, as well as the previously published mitochondrial sequence of A. issoria.

29The sequence data for 20 nuclear and 2 mitochondrial genes encompass 70 out of 76 (92%) of the officially recognized species of Heliconiini, including 44 out of 46 species from the focal genus Heliconius. Although the taxonomic validity of some species is contested, we found that the diversification analysis is robust to altering the number of missing species. Some recognized taxa diverged very recently, as shown by the BEAST chronogram (Fig. 1a) in case of the Philaethria diatonica/Philaethria neildi/Philaethria ostara complex and the H. heurippa/H. tristero pair, the latter of which was recently reclassified as a race of Heliconius timareta (Mérot et al. 2013). However, the exact relationship between genetic differentiation and taxonomic species identity in the highly variable, mimetic Heliconiini remains unclear (e.g.,; ). Importantly, 36 species are represented by multiple individuals, usually from distant populations, allowing for more accurate estimation of sequence evolution rates, and detection of species paraphyly. The number of individuals represented by each marker ranges from 40% for Hcl to 98% for CoI/II, and only four specimens are represented exclusively by mitochondrial DNA (online Appendix S1, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j).

30New DNA sequences were deposited in GenBank under accession numbers KP072800-KP074896 and KP113715-KP114069 (online Appendix S1, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). Sequence alignments and phylogenetic trees were deposited in TreeBase (accession number 15531).

31Our nuclear markers span 11 out of 21 chromosomes (Supplementary Table S2, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j; ) and have both autosomal and sex chromosome Z-linked inheritance. We examined the conflict between individual markers in the entire tribe and in the genus Heliconius alone, using both gene tree summary methods and approaches utilizing the raw sequence alignments. The ML analysis of the core matrix in Concaterpillar rejected concatenation of any of the loci due to significant differences in both topology and substitution rate of individual partitions, but the exact nature of the discordance is unclear. A Multi-Dimensional Scaling ordination of pairwise RF distances between the gene trees does not reveal clustering by chromosome and the separation between many nuclear loci appears much greater than between nuclear and mitochondrial trees (Supplementary Fig S1, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). Consistent with this is the fact that the whole mitochondrial phylogeny of select taxa (Fig. 3) shows few differences from the tree based on the mixed marker supermatrix (Fig. 1a), highlighting that cytonuclear discordance is not the primary source of incongruence in the data set.

32The coalescent approaches reveal the high extent of marker conflict in the Heliconius data. Gene tree topologies from the explicit Bayesian modeling of ILS in *BEAST are highly varied, with a particularly high degree of reticulation in the H. melpomene/cydno and the Heliconius hecale (silvaniform) clades, where extensive horizontal gene flow has been observed previously (Supplementary Fig. S2, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j;; ). Another Bayesian method, BUCKy, infers the species tree in the presence of marker incongruence without modeling specific reasons for the observed discordance, and calculates the Bayesian CFs that illustrate the proportion of partitions in the data set that support a particular grouping. The concordance of the loci for Heliconius is strikingly low (Fig. 3), although the topology is consistent with the results of other analyses (Figs. 1a and 3; Supplementary Figs. S3–S7, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j).

33Further strong evidence of widespread incongruence comes from the NeighborNet network characterized by a high delta score of 0.276, which shows that the structure of the data is not entirely tree-like (Supplementary Fig. S4, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). This can be partially attributed to the effect of missing data, yet even a fit based on the 30 species with full or nearly complete sequences produced a delta score of 0.11, proving a substantial amount of non-bifurcating signal across the tribe. The most noticeable reticulations are found between nodes linking genera and the major clades of Heliconius, possibly due to pervasive gene flow during the diversification of the main extant lineages (Supplementary Fig. S4, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j).

34Although no trees are identical, the results from our Bayesian, ML, and distance-based network analyses of the supermatrix are very similar (Fig. 1a, Supplementary Figs. S4 and S5, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). Two basal nodes stand out as unstable. Cethosia is variably placed as a sister taxon to either Acraeini (MP, ML, NeighborNet) or Heliconiini (Bayesian), despite the reasonably extensive sampling of 11 loci for Cethosia cyane, while the position of Podotricha in relation to Dryas and Dryadula varies among all analyses. Most problematic are relations among the species within Eueides, where the position of 4 out of 10 taxa cannot be resolved with good support. The poor resolution for both Eueides and Podotricha can probably be attributed to insufficient site coverage, which produces high uncertainty due to patchily distributed missing data. We re-estimated the relations of Eueides based solely on a core set of 11 genes with coverage for at least 7/10 species and recovered a much better supported tree (Supplementary Fig. S6, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). We recommend studies focusing on Eueides use this specific phylogeny, but the exact relations of Eueides procula, Eueides lineate, and Eueides heliconioides remain unclear.

35Although our Bayesian maximum clade credibility tree is similar to the topology estimated by, significantly increasing the data set from 113 to 180 individuals and 5 to 21 loci allowed us to resolve many critical nodes. The major differences are in the relations of the genera other than Heliconius, which we infer to form a grade, with high support for Eueides as the sister genus of Heliconius (Fig. 1a). Importantly, we confirm that the enigmatic genera Laparus and Neruda are nested within Heliconius, as further supported by ELW and SH tests (Table 2), although Neruda is closer to the base of the genus than previously estimated (Fig. 1a). We find that the other so called “primitive” clade of Heliconius consists of two separate groups nested among other subgenera raising questions about the apparently unequal rate of morphological evolution in the genus.

36The ML tree is similar to the Bayesian phylogeny, although it puts less confidence in the deeper nodes linking the genera and subgenera (Supplementary Fig. S5, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j), as is typically the case. The nodes differing between the two trees are also the nodes that cannot be unequivocally confirmed by either the SH and the ELW test (Supplementary Table S4, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). For instance, the highest posterior probability in the ELW test (0.302) is given to the tree that does not lump Neruda and the “primitive” Heliconius, consistent with the BEAST chronogram in Figure 1a. Thus, we suggest that the Bayesian tree should be preferred as a more accurate picture of the systematic relationships, although the ML phylogeny based on the complete data set is still useful to uncover multiple polyphyletic species. Notably, Heliconius luciana is nested within Heliconius elevatus, and Heliconius wallacei is polyphyletic with respect to Heliconius astraea. These results must be interpreted with caution, as the inference relies on poorly covered museum specimens and may be sensitive to long branch attraction.

37Our whole mitochondrial phylogeny is largely consistent with the results of the multilocus supermatrix analysis and well-supported for 46/57 nodes (Fig. 3), despite the relatively limited taxonomic coverage of only 29 Heliconiini for which short-read data were available. In contrast to the multilocus data set, mitochondrial genomes are not very useful for resolving relationships between major clades within Heliconius, as the positions of Neruda, Heliconius xanthocles, and H. wallacei clades are poorly supported. An important deviation from the predominantly nuclear multilocus phylogeny is the placement of Heliconius pachinus as a sister taxon of H. timareta, rather than H. cydno. This may reflect the overall instability resulting from a high extent of reticulation in the melpomene/cydno/timareta assemblage. Furthermore, we find a surprising positioning of Heliconius hermathena within H. erato, which is also not supported by any of the analyses of the 21 locus matrix. The mitochondrial tree confirms previous observations of deep biogeographical splits in the widespread, highly diversified H. erato and its co-mimic H. melpomene. Within H. melpomene we find a well-supported distinction between races found to the west and east of the Andes, although our data place the individuals from French Guiana with the specimens from the Western clade, in contrast to a whole-genome phylogeny (Fig. 1a). Heliconius erato shows the opposite pattern in both mitochondrial and nuclear data, whereby the Guianian races form a fully supported clade in the monophyletic group of taxa from East of the Andes.

38The variety of MSC methods that we applied brings a new perspective to the phylogenetic signals in Heliconiini. We analyzed the small data set of the 17 best-sampled species with the Bayesian MSC algorithm implemented in *BEAST and recovered a tree which agrees with the supermatrix analyses (Supplementary Fig. S3, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j), except for the position of Neruda aoede, which was placed as a sister taxon to the H. xanthocles/L. doris clade with relatively low posterior probability. Furthermore, the mean ages of nodes are very close to those proposed in a supermatrix analysis (Table 1), with similar 95% credible intervals. Although the species tree is largely as predicted, we observe high levels of incongruence in the underlying distribution of gene trees (Supplementary Fig. S2, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). Differences in the depth of coalescence are clear throughout the tree and reticulation is again especially apparent in the H. melpomene/cydno clade. The estimated population size values are also consistent with a previous comparison based on two nuclear loci, showing a higher population size of H. erato ( individuals) when compared with H. melpomene.

39Similarly, we find the phylogeny derived by the gene tree summary approach BUCKy to be entirely consistent with the Bayesian analysis of concatenated sequence, although the recovered CFs are much lower than any other measure of support applied to our data. Importantly, most of the nodes connecting the major clades in the tree have CFs below 0.5, with the notable exception of the silvaniform/melpomene split (Fig. 3). The same nodes correspond to the reticulations in the NeighborNet analysis (Supplementary Fig. S4, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j), cases of low support in the MDC tree and its disagreement with the supermatrix analysis (Supplementary Fig. S7, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j), nodes that cannot be rejected in the ML and Bayesian tests of topologies (Supplementary Tables S4 and S5, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j), and the uncertain nodes in the whole mitochondrial tree (Fig. 3).

40Another summary analysis by MDC, although taking only a few minutes with 100 bootstrap replicates of the complete data set, returns a very different topology from the other techniques, showing a number of unexpected and poorly supported groupings. The lumping of all non-Heliconius genera, and the monophyletic Neruda/xanthocles/wallacei clade stand out in contrast to other proposed trees (Supplementary Fig. S7, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). Interestingly, many of the relations that are poorly supported in the supermatrix phylogenies are also not resolved in the consensus MDC tree, showing that MDC is highly conservative with regard to the placement of taxa unstable in individual gene trees.

41The phylogeny estimated under a relaxed clock model in BEAST shows diversification dynamics that differ from previous estimates, with the deeper splits between the genera of Heliconiini estimated as substantially older than previously inferred with mitochondrial data, but younger than estimated with a small sample of nuclear genes (Table 1). The most species-rich genera Heliconius and Eueides separated 18.5 (95% highest posterior density: 16.5–20.4) Ma and both started to diversify, respectively, 11.8 (10.5–13.4) Ma and 10.2 (8.9–11.7) Ma. The six major clades of Heliconius (corresponding to H. erato, Heliconius sara, H. xanthocles, H. wallacei, H. melpomene/silvaniforms, and Neruda) all started to diversify around 5 Ma (Fig. 1b, Supplementary Fig. S8b, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). The extinction rate has remained relatively constant throughout the history of the tribe, while the speciation rate (and thus net diversification) increased substantially in Eueides and Heliconius (Table 2). The LTT plot (Supplementary Fig. S8a, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j) and the plot of diversification rate for Heliconiini (Supplementary Fig. S9, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j) suggest rapid early emergence of Podotricha, Dryadula and Dryas, which later speciated at a very low rate (Table 2). This is followed by a period of stasis 18–12 Ma, roughly corresponding to the mid-Miocene, and a sudden but steady increase in the rate of diversification after 11 Ma. In the case of the 45 Heliconius, a shorter plateau is found between 9 and 6 Ma, and the number of extant lineages rises from 5 Ma onwards (Supplementary Fig. S8b, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). As expected, the LTT plots drop off sharply in the last million years, reflecting protracted speciation.

42ML modeling in DDD strongly supports the Birth–Death model with an increase in diversification rates 11 Ma as the best fit for Heliconiini (Akaike weight of 0.92; Supplementary Table S6, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). However, we find that the confidence in the point estimate of shift time is largely a function of forcing the model over the data. The Bayesian analysis with BAMM demonstrates that although there has likely been a large change in rate regimes across the phylogeny, it cannot be reduced to a single event at a strictly specified time, and this result is consistent for all four chronograms that we tested (values are reported for the MCC chronogram based on multiple uniform calibrations). Models with one or two shifts in diversification rate have the highest frequency in the posterior (0.430 and 0.281, respectively) and Bayes factors provide weak evidence for models with one shift (BF = 3.495) or two shifts (BF = 4.398) over models with no shift. Painting the mean rate per branch from 2000 individual models onto the phylogeny (Fig. 1b) and explicit testing of rate variation between clades (Supplementary Fig. S10, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j) show clearly that the rate changed on the branches leading to (i) the Dryas and Podotricha clade (BF >40); (ii) the Heliconius and Eueides clade (BF >10); and (iii) Heliconius. But even though there is evidence for rate increase on the Heliconius branch, no particular model with a specific shift time occurs at high frequency. On the contrary, the most frequent model in the 95% credible set is the one with a constant rate.

43We approached the problem of phylogeny reconstruction in a difficult mimetic assemblage through extensive intraspecific sampling of 22 markers from nearly all species in the clade, and compared the results between multiple philosophically distinct analytical approaches. As next generation sequencing technologies become widely accessible and the average number of loci used in systematic studies increases rapidly, MSC methods gain in importance as means of detecting and accounting for incongruence in multilocus data (e.g.,;; ). However, their relative merits and utility at different levels remain contested. The systematic relations of the tribe Heliconiini, which diverged from its extant outgroup about 47 Ma, can be effectively resolved with both MSC and supermatrix approaches, yielding highly similar topologies across a range of different subsampling schemes that correspond to the requirements of individual algorithms (Figs. 1a and 3; Supplementary Figs. S3–S7, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). Nonetheless, consistent with the recent radiation of the group, large effective population sizes and known hybridization between many species, we observed high heterogeneity among sampled fragments of the genome that differ markedly in both topology (Fig. 2, Supplementary Fig. S2, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j) and rates of evolution (Concaterpillar analysis). Such heterogeneity might have been expected to pose a significant challenge for the concatenation methods. However, the only method producing an obviously different phylogeny is MDC (Supplementary Fig. S7, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j), which fails to resolve 34 out of 62 nodes in the tree with bootstrap support above 0.9, and is the only method to suggest monophyly of the non-Heliconius genera. MDC derives a species tree from point estimates of gene trees and can be expected to perform poorly with a relatively limited number of gene trees that are not always fully resolved, leading to a complete polytomy in some of the clades. However, the MDC result is an indicator of instability, as well-resolved and consistent gene trees should produce a good quality MDC tree.

44Recent studies have proposed that likelihood-based MSC techniques should be preferred to integration over individual gene trees, due to their potential to capture synergistic effects between partitions, mirroring the phenomenon of hidden support in concatenation. Our results support this observation and further show that high degrees of conflict between many partitions can be reconciled by both supermatrix and MSC approaches to extract the predominant signal of speciation. The BCA in BUCKy assigns insignificant CFs to most of the nodes separating major subgenera of Heliconius. Two of these nodes (H. wallacei and N. aoede; Fig. 3) are also only weakly supported by the *BEAST coalescent model, and correspond to the areas of high reticulation in the NeighborNet network, reflecting conflicting signals (Supplementary Fig. S4, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). The same nodes are all assigned a posterior probability of one in the Bayesian supermatrix analysis (Fig. 1a), potentially leading to the erroneous conclusion that all the data point unequivocally to the inferred relations. The superiority of the MSC methods lies in the effective demonstration of incongruences, represented by lower support values or CFs assigned to the more difficult nodes.

45Our study highlights an important practical consideration in choosing the optimal analytical approach, where the requirements of the selected algorithm have to be reconciled with a realistic sampling of taxa. Our final goal of testing hypotheses regarding the diversification dynamics of Heliconiini can be met only if the number of included taxa is maximized. Despite a substantial effort we only managed to secure single samples of the rare or geographically restricted species, and some of them are represented solely by historical specimens with limited potential to generate extensive multilocus data. Six other species remain missing, as they are known only from a small number of specimens (Neruda metis, E. emsleyi) or found at low density in infrequently visited areas (Philaethria andrei, Philaethria browni, Philaethria romeroi, and Heliconius lalitae). Specimens of four of these were available to us, but poor preservation of material by desiccation made it impossible to extract reasonable quality DNA. Considering that our focal group is intensively studied and exceptionally well represented in research, museum, and private collections due to its aesthetic appeal, it would be considerably more challenging to obtain a complete sampling of many other groups. Another difficulty stems from the fact that the advanced coalescent techniques like BUCKy and *BEAST perform best with multiple samples per species, which should capture intraspecific diversity, and may require complete taxon coverage. Much of the uncertainty in the estimates can be attributed to missing data, which can negatively affect the estimation of both individual gene trees and the encompassing species tree. When fitting a NeighborNet network, we found that although the percentage of data missing from the matrix does not explain all of the observed reticulation and the high delta score, it causes these parameters to increase, thus suggesting that the data completeness at each alignment position must be maximized to identify genuine incongruence. We observe that in many cases of biological interest it will be a formidable challenge to generate the ideal data set that (i) has little missing data, (ii) comprises a large, genome-wide sample of loci, (iii) includes all taxa, and (iv) captures intraspecific variability. In case of Heliconiini, the supermatrix approach based on a limited number of markers (22) helped us to maximize taxonomic inclusiveness without compromising our ability to reconstruct a phylogeny in the light of conflicting biological signals.

46Heliconiini have been at the center of the debate about species concepts and designation criteria, providing empirical evidence for the permeability of species barriers. Species delimitation with the PTP based on relative branch lengths in the ML tree produced very surprising, although poorly supported results, strongly at odds with biological knowledge of Heliconiini. All the lineages in the H. melpomene/H. cydno group were lumped into a single Operational Taxonomic Unit (Supplementary Fig. S7, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). This is despite extensive evidence for at least three species (H. melpomene, H. cydno, and H. timareta), if not the currently recognized five (also including H. pachinus and H. heurippa): premating behavioral isolation, host plant usage and habitat preference, morphology, and genome-wide genetic structure. The small number of nucleotide differences between the biological species is consistent with interspecific gene flow, which led to a reduction in differentiation across the genome. Although automated species delineation is certainly a useful criterion for initial assessment of poorly understood taxa or for rapid inventory of entire communities, we argue strongly that traditional taxonomies based on accumulation of extensive biological knowledge should take precedence in case of intensively researched taxa like Heliconiini. This is an important result, considering that the mitochondrial barcode region, although frequently used to delimit arthropod species at low cost, has also been shown to be insufficient to distinguish species in Ithomiini, a butterfly tribe co-mimicking Heliconiini.

47Application of novel algorithms to Heliconiini data reveals a relatively stable topology and branch lengths, despite limited support (Figs. 1a and 3; Supplementary Figs. S2–S7, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). We uphold the previous reclassification of the species in the genera Neruda (Hübner 1813) and Laparus (Linnaeus 1771) as Heliconius based on their nested position. This placement is at odds with morphological evidence from adult and larval characters, which puts Neruda and Laparus together with Eueides, suggesting convergent evolution of homoplasious morphological characters. Nevertheless, the molecular evidence is decisive and we thus synonymize Laparus syn. nov. and Neruda syn. nov. with Heliconius. We also conclude that the traditional naming of some Heliconius as “primitive” (e.g., ) is phylogenetically unjustifiable and misleading. Instead, we propose to call these groups “the wallacei clade” (including H. astraea, Heliconius burneyi, Heliconius egeria, and H. wallacei) and “the doris clade” (Heliconius doris, Heliconius hecuba, Heliconius hierax, and H. xanthocles).

48The position of the enigmatic genus Cethosia remains unresolved, as it currently depends on the chosen method of analysis, and is unlikely to be established without a broad sampling of species using multiple markers. Cethosia has been variably considered to be either the only Old World representative of Heliconiini, a genus of Acraeini, or possibly a distinct tribe. Establishing the systematic relations between Acraeini, Cethosia, and Heliconiini is important for the study of Heliconiini macroevolution, as it will shed light on the dispersal route of Heliconiini into the Neotropics.

49Our phylogeny of Heliconiini brings novel insights into the diversification dynamics of the clade. Although most age estimates for Heliconius agree with other studies, the deeper nodes are older than previously suggested (Table 1). There is little agreement on the dates above the species level, and the studies to date either suffer from insufficient taxon sampling (;;; Table 3), or use markers unlikely to be informative above a relatively low level of divergence. For instance, the mean age of the split between Heliconius and Agraulis is estimated as 32 Ma in the study of, which includes only three species of Heliconiini; 26.5 Ma in, including one species per genus; or 23.8 Ma in the present study (Fig. 1a; Table 1). Conversely, find the divergences between the genera of Heliconiini to be much younger than we propose, likely due to an effect of using a fast-evolving mitochondrial locus without partitioning. estimate broadly similar dates, but do not formally test alternative hypotheses of diversification and their conclusions are limited by incomplete sampling of the most diverse genera.

50Our own ability to correctly falsify the hypotheses regarding the diversification of Heliconiini hinges on having a nearly complete phylogeny, yet none of our MSC analyses consider as many species of Heliconiini as the Bayesian supermatrix estimate. We are confident that the values proposed by the supermatrix method can be trusted, as both the topology and the branch lengths are consistent with the results inferred by *BEAST based on a smaller data set (Table 1; Supplementary Fig. S3, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). We find that the ages of the deeper nodes and the length of terminal branches are not inflated by the supermatrix method in comparison to *BEAST, contrary to the predictions from simulations, and both methods infer similar mean age for the observed splits, for instance 11.8 Ma for the basal divergence of Heliconius into the H. erato and H. melpomene lineages, or around 3.5 Ma for H. doris and H. xanthocles. Hence, we offer a new perspective on the dating of the Heliconiini radiation with a nearly complete set of divergence time estimates based on a well-resolved and supported tree.

51Despite evolving many properties unique among the Lepidoptera (reviewed by ), Heliconius has not undergone a rapid adaptive radiation involving an explosive growth in the number of lineages, followed by a dramatic slowdown, as might be expected due to a sudden availability and filling of ecological niches. In fact, the rate of diversification has been increasing steadily for the last 11 myr (Supplementary Fig. S10, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). This result is consistent with a recent critique of the traditional prediction that large radiations should display a pattern of the initially high net diversification rate decreasing as the ecological niches fill up. Such an expectation is reasonable for island radiations that constitute a large proportion of study cases to-date, but it does not necessarily apply to continental radiations in the tropics, where the scale and complexity of the ecosystems are likely to generate a number of suitable niches greatly exceeding even the cladogenetic potential of large radiations. To date, steady diversification of a widely distributed taxon has been demonstrated in the 129 African Synodontis catfish and the 293 Neotropical Furnariidae ovenbirds, but the generality of this pattern remains unknown and requires further verification in other well-sampled large groups like Heliconiini.

52Both ML (DDD) and Bayesian (BAMM) approaches provide evidence for substantial increases in speciation without changes in extinction rate, although we cannot completely exclude the possibility that extinction has not been modeled accurately and may have affected our findings. The first surge in speciation occurred early on the branch leading to the presently depauperate clade of Podotricha telesiphe, Podotricha judith, and the monotypic genera Dryas and Dryadula, and stands out against the low overall background rate of diversification (Fig. 1b, Supplementary Fig. S10, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). More substantial increases are observed in the Heliconius and Eueides clade (Heliconius mean net diversification: 0.196 new lineages per lineage per myr; non-Heliconius rate: 0.101 (Fig. 1b, Supplementary Fig. S11, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). Bayesian modeling shows that the gains have been gradual and cannot be correlated with a single environmental factor, although the speciation of Heliconius and Eueides may have been stimulated by the second rapid stage of Andean orogeny approximately 12 Ma, which strongly changed the elevation gradients in the Central and Eastern sectors. This was contemporaneous with climatic changes that affected the distribution of the rainforest, and followed shortly by the entrenchment of a major barrier, the Amazon, in its modern course 10 Ma. Nonetheless, our results do not strongly support such a correlation and highlight the dangers of the commonly used approach of selecting a single rate shift configuration (reviewed in ). For example, analysis of the Heliconiini data with DDD provided false confidence in a single distinct shift to a much higher rate of diversification. The lack of strong response by Heliconiini to the environmental perturbations is less surprising if we consider that intense changes on the continental scale have occurred nearly continuously from the start of Miocene (and thus Heliconiini) until the present. This has created a dynamic arena for constant evolution of new species and colonization of novel ecological niches, across a variety of environmental regimes, ecosystems, and geological formations.

53More puzzling is the gradual rise in the diversification rate of Heliconiini, clearly driven by the increased speciation rates of Heliconius (Fig. 1b, Supplementary Fig. S10, available on Dryad at http://dx.doi.org/10.5061/dryad.44b4j). One distinct possibility is that the extraordinary diversity of hundreds of mimetic patterns in the genus (see Fig. 1a for some examples) has facilitated speciation and stabilized extinction, leading to a positive feedback of higher number of species on the diversification rate. However, testing any ideas on the macroevolutionary impact of the mimetic phenotypes will be difficult. For instance, direct optimization of the evolution of color patterns on the phylogeny proved fruitless, as the signal of ancestral morphological states is lost behind the existing intraspecific diversity of widely varied races and the confounding effects of pattern elements introgressing between species.

54Some uncertainty surrounds the last few million years of evolution of Heliconiini. Early speculation regarding the drivers of speciation suggested diversification in allopatry, as the rainforest habitat occupied by most species in the group has undergone cycles of contraction and expansion in response to recent climatic variation. The hypothesis of vicariant cladogenesis has been subsequently criticized due to a lack of evidence for forest fragmentation in pollen core data, lack of temporal concordance of divergence of distinct lineages, and the likelihood of parapatric speciation. The decrease in observed diversification over the last million years may be due to protracted speciation, that is, our limited ability to delineate species in the assemblages of highly variable taxa such as Heliconius.

55In conclusion, we present a taxonomically comprehensive phylogeny of a large continental radiation characterized by extensive introgression and varying rate of diversification. Reticulate signals in phylogenetic data will increasingly have to be accounted for in study designs and could be common in other recent adaptive radiations, as recently reported for swordfish, mosquitoes, or broomrape plants. Ignoring the possibility of reticulation no longer seems a viable assumption in phylogenetic analysis of recent adaptive radiations.

56Conceived and designed the study: K.M.K., N.W., and C.D.J. Carried out the labwork, analyzed the data and drafted the manuscript: K.M.K. Provided samples and next-gen sequence data: C.D.J, J.M, K.K.D., N.W., and A.N. All authors have read and approved the final version of the manuscript.

57Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.44b4j.

58This work was supported by the Harvard Herchel Smith Trust, Emmanuel College and the Balfour Studentship from the Department of Zoology, Cambridge University [to K.M.K.]; the Leverhulme Trust Leadership Grant and the BBSRC grant H01439X/1 [to C.D.J.]; Kone Foundation [to N.W.]; the BBSRC [grant BB/G006903/1 to J.M. and K.K.D.].