Contrasting patterns of Andean diversification among three diverse clades of Neotropical clearwing butterflies

Abstract The Neotropical region is the most biodiverse on Earth, in a large part due to the highly diverse tropical Andean biota. The Andes are a potentially important driver of diversification within the mountains and for neighboring regions. We compared the role of the Andes in diversification among three subtribes of Ithomiini butterflies endemic to the Neotropics, Dircennina, Oleriina, and Godyridina. The diversification patterns of Godyridina have been studied previously. Here, we generate the first time‐calibrated phylogeny for the largest ithomiine subtribe, Dircennina, and we reanalyze a published phylogeny of Oleriina to test different biogeographic scenarios involving the Andes within an identical framework. We found common diversification patterns across the three subtribes, as well as major differences. In Dircennina and Oleriina, our results reveal a congruent pattern of diversification related to the Andes with an Andean origin, which contrasts with the Amazonian origin and multiple Andean colonizations of Godyridina. In each of the three subtribes, a clade diversified in the Northern Andes at a faster rate. Diversification within Amazonia occurred in Oleriina and Godyridina, while virtually no speciation occurred in Dircennina in this region. Dircennina was therefore characterized by higher diversification rates within the Andes compared to non‐Andean regions, while in Oleriina and Godyridina, we found no difference between these regions. Our results and discussion highlight the importance of comparative approaches in biogeographic studies.


| INTRODUCTION
The formation of mountains is a major geological event that results in profound changes in the topography, climatic conditions, and water drainage that are likely to influence the timing and geography of diversification. Mountains may act as a barrier that isolates populations on both sides or forms an island-like archipelago (e.g., Hughes & Eastwood, 2006), thereby driving vicariant speciation events.
Climatic turnover and complex topography along the slopes allow the establishment of a large variety of habitats, vegetation, predator, and pathogen communities and may in turn affect diversification (Badgley, 2010). A diversity of environmental and ecological conditions provides multiple opportunities for adaptation and ecological speciation. The distribution of poikilotherm phytophagous insects for example will be directly determined by temperature, rainfall, and solar radiation (e.g., Menéndez et al., 2007), as well as by the plant community that hosts their larval stages and that is known to also act as an important driver of diversification (Ehrlich & Raven, 1964;Janz, Nylin, & Wahlberg, 2006). From a biogeographic point of view, mountain ranges not only generate local diversification along their slopes, but they can also feed adjacent areas through dispersal events, potentially enhancing diversification in neighboring regions. Assessing the timing of diversification and dispersal events with respect to mountain uplift is therefore of primary importance in understanding the origins of many modern biotas.
The formation of the Andean cordillera that extends from northern Venezuela to southern Chile has been proposed as the main driver of diversification in the Neotropical region (Hoorn et al., 2010).
Several scenarios have been proposed to explain the role of the Andes in Neotropical diversification, but there has been confusion surrounding these scenarios and the actual processes underlying each of them. Chazot et al. (2016) proposed a clarified framework of four nonmutually exclusive diversification scenarios with respect to the Andean mountains, based on the assumption that a species pool of a biogeographic region results from the processes of speciation, extinction, dispersal, and the amount of time the region has been occupied. These scenarios and their predictions are as follows. (1) Cradle scenario. The Andes promote vicariant speciation and ecological speciation across and along the slopes, as supported for instance by the extremely high rates of speciation inferred in some Andean groups of plants (Madriñán et al. 2014). Under this scenario, Andean diversity is the result of increased speciation rates in Andean lineages compared to other regions.
(2) Museum scenario. The Andes may have provided more stable environments during periods of climate change and hence may have saved lineages from extinction. Under such a scenario, Andean diversity is the result of lower extinction rates of Andean lineages compared to other regions (Stebbins, 1974). (3) Species-attractor scenario. Lineages in areas adjacent to the Andes may have taken advantage of newly available Andean niches to colonize the slopes of the Andes multiple times (Brumfield & Edwards, 2007;Chazot et al., 2016). In this scenario, the colonization rate into the Andes is higher than the colonization rate out of the Andes. (4) Time-for-speciation scenario. If Neotropical clades historically originated in the Andes before spreading into the rest of the Neotropical region, they will have accumulated species over longer periods of time, regardless of differences in diversification and dispersal. Under such a scenario, the first colonization time of the Andes is higher than the first colonization time of non-Andean regions (time-for-speciation hypothesis, Stephens & Wiens, 2003).
Here, we investigated spatial and temporal patterns of diversification and assessed support for these four biogeographic scenarios in Andes, biogeography, Dircennina, Ithomiini, Lepidoptera, Neotropics, Oleriina, trait-dependent diversification two Neotropical butterfly clades, the ithomiine subtribes Dircennina and Oleriina (Nymphalidae: Danainae), and we compare them to published diversification patterns of a third ithomiine subtribe, Godyridina (Chazot et al., 2016). Both clades belong to the nymphalid tribe Ithomiini, one of the best-studied groups of Neotropical butterflies (Mallarino, Bermingham, Willmott, Whinnett, & Jiggins, 2005;Brower et al., 2006;Elias et al., 2009;Brower, Willmott, Silva-Brandão, Garzón-Orduña, & Freitas, 2014;Garzón-Orduña, Silva-Brandão, Willmott, Freitas, & Brower, 2015;Chazot et al., 2014Chazot et al., , 2016De-Silva et al., 2010, 2017. These three clades (the three largest ithomiine subtribes with 101 known species of Dircennina, 77 Godyridina, and 64 Oleriina, representing over 60% of ithomiine diversity) are endemic to the Neotropical region and occupy forest habitats from Central America to the Atlantic Forest, from the lowlands to high altitudes in the Andes. A recent study of spatial and temporal patterns of diversification in Godyridina, where nearly 60% of the species are found in the Andes, revealed that this subtribe originated probably at the interface between the upper Amazon region and the Central Andes about 17 million years ago (Chazot et al., 2016). The Godyridina diversification pattern conforms to the species-attractor scenario, with repeated colonization of the Andes, and the subtribe also underwent local radiations in the Northern Andes, in the Central Andes, and in the Upper Amazon (Chazot et al., 2016).
Concerning Dircennina and Oleriina, although both subtribes are diverse in Andean regions , 2017, they show contrasting pattern of species distribution among Andean and non-Andean regions. Dircennina are unusually species-rich in the Andes compared to the rest of the Neotropics (63% of their diversity is in the Andes) whereas Oleriina have similar species richness in Andean and non-Andean regions (49% of their diversity is in the Andes, with relatively high diversity also in Amazonia), which suggests the possibility of different scenarios of diversification. The three ithomiine subtribes therefore represent an excellent system for investigating how closely related clades have been affected by the Andes during their evolution and for improving our understanding of the mechanisms involved in diversification.
A phylogeny of Oleriina has been published (De-Silva et al., 2010, and the phylogeny of the richest dircennine genus, Pteronymia (53 species, i.e., about half of the subtribe), has just been generated (De-Silva et al., 2017). Here, we compile new and published sequences to confirm and revise, when needed, the taxonomy in the remaining dircennine genera, and to generate the first time-calibrated molecular phylogeny of the entire subtribe Dircennina. For the Oleriina, we use the time-calibrated phylogeny published by De-Silva et al. (2016).
We investigate the spatial pattern of species diversification of the two subtribes using the framework developed for the Godyridina (Chazot et al., 2016) as follows. (1) Using an explicit biogeographic model, we perform a biogeographic ancestral area reconstruction. (2) Using models of trait-dependent diversification, we assess the support of the four biogeographic scenarios outlined above for Andean diversification: Andean lineages had a higher speciation rate (cradle hypothesis), Andean lineages had a lower extinction rate (museum hypothesis), the Andes were colonized at a higher rate (species-attractor hypothesis) or earlier (time-for-speciation hypothesis) than non-Andean regions.

| MATERIAL AND
As the taxonomy of the genus Pteronymia has already been revised by De-Silva et al. (2017), here we only included one representative per species from that genus, corresponding to the consensus sequences of all individuals available for each species (see De-Silva et al., 2017 for more details). PCR conditions followed Elias et al. (2009) for COI-COII, EF1, and tektin and Wahlberg and Wheat (2008) for CAD, GAPDH, MDH, and RPS2.

| Dircennina individual-level phylogeny
We aligned sequences with CodonCode Aligner v6.0.2. The molecular dataset was partitioned by gene and codon positions. We performed a maximum-likelihood analysis including all individuals using IQ-TREE software as implemented in the W-ID-TREE server (Nguyen, Schmidt, von Haeseler, & Minh, 2015;Trifinopoulos, Nguyen, von Haeseler, & Minh, 2016). IQ-TREE automatically selected the best partition scheme, and we performed 1000 ultrafast bootstrap analyses (Minh, Nguyen, & von Haeseler, 2013). The tree can be found in Appendix S2.

| Dircennina species-level phylogeny and dating analyses
We computed a consensus sequence for each species and gene region, resulting in a dataset containing 86 species and 44 outgroups.
We used the "greedy" algorithm and linked rates implemented in PartitionFinder 1.1.1 (Lanfear, Calcott, Ho, & Guindon, 2012) to select the best models of substitution for optimized sets of nucleotides over all models implemented in BEAST. A time-calibrated phylogeny was generated using BEAST 1.8.2 (Drummond, Suchard, Xie, & Rambaut, 2012), using the best partition scheme and uncorrelated lognormal relaxed clock for each partition. We applied eight secondary calibration points based on Nymphalidae node ages obtained from Wahlberg 100,000 generations, resulting in 1,000 trees. The maximum clade credibility tree using the median of posterior distribution for node ages was extracted using TreeAnnotator (Drummond et al., 2012), applying a 10% burn-in (Appendix S4-S5). We also ran two independent analyses to assess the effect on node ages of using a Birth-Death tree prior or a Yule prior. Results were extremely similar (Appendix S5), and we used the Birth-Death tree for all analyses. Following Chazot et al.'s (2016) analyses of Godyridina diversification, we conducted two types of biogeographic analyses for each subtribe to assess the relative support for the four scenarios of diversification (cradle, museum, species-attractor, time-for-speciation).

| Biogeographic analyses
First, we divided the Neotropical region into nine areas to infer past distributions of ancestral lineages. We followed Morrone's (2014) classification of Neotropical biogeographic to define these nine areas ( Figure 1): Central America, lowlands adjacent to the western slopes of the Northern Andes, Central Andes, western and central cordilleras of the Northern Andes, eastern cordillera of the Northern Andes, Guiana Shield, upper Amazon, lower Amazon, and the Atlantic Forest. We inferred ancestral biogeographic areas under six different models of biogeographic reconstruction. Each model accounts for different range-changing processes, controlled by specific parameters. We used the models DIVALIKE, DEC, and BAYAREALIKE as implemented in BioGeoBEARS v.0.2.1 (Matzke, 2014), and each model was fitted with and without the founder-event speciation parameter (j). The differences among these three models are that DIVALIKE accounts for vicariant speciation events occurring for widespread ranges, DEC accounts for sympatric events of speciation where one descendant only inherits a subset of the ancestral range and BAYAREALIKE accounts for sympatric events of speciation where the two descendants inherit the entire ancestral range.
Finally, the founder-event speciation parameter (j) accounts for speciation events where one of the descendants occupies an area that was not occupied by the ancestor. These six models for each phylogeny were compared using AIC scores. De-Silva et al. (2016)  ing to the aforementioned procedure.
As a second step, we performed an analysis that aimed to statistically test the role of the Andes in the biogeography of the Neotropics and understand the origin of high Andean species richness. We assigned species to Andean or non-Andean regions and fit models of character state-dependent speciation and extinction to test whether the Andes (1) had a higher speciation rate than non-Andean lineages, as expected under the cradle hypothesis, (2) had a lower extinction rate as expected under the museum hypothesis, (3) had a higher rate of colonization into the Andes than out of the Andes as expected under the species-attractor hypothesis, or (4) were colonized before non-Andean regions as expected under the time-for-speciation hypothesis.
This follows the framework proposed by Chazot et al. (2016). Using locality and associated elevation records, species were assigned to either belonging to non-Andean (region 1) or Andean (region 2) regions.
Amazonian species whose distribution slightly overlaps with the low altitude Andean foothills were considered non-Andean (Chazot et al., 2016). We applied the ClaSSE model of diversification, which can incorporate up to 10 parameters (Goldberg & Igic, 2012). However, we reduced the number of models following Chazot et al. (2016). We constrained to 0 the speciation parameters involving a state transition in both descending lineages (λ 122 and λ 211 ) and the anagenetic state transition rates (q 12 and q 21 ) because these parameters correspond to unrealistic biogeographic events. Therefore, the most complex model included four speciation parameters that are biogeographically meaningful (within region speciation rates: λ 111 λ 222 , transition rates between regions: λ 121 , λ 212 ) and two extinction parameters (μ 1 , μ 2 ). The models fitted included all combinations of one or more parameters free to vary. The resulting 10 models were compared using the AIC scores. Differences of two units of AIC between models were considered as a significant improvement of the fit. In addition, we performed MCMC analyses on the best fitting models as a further exploration of parameter estimates (Fitzjohn 2012). We used exponential priors and a 20,000-steps chain. Parameters converged rapidly and we discarded the first 10% steps as burn-in before visualizing the distributions parameter estimates.
Finally, we used trait-dependent diversification models to infer ancestral states. However, ancestral state reconstruction is not implemented in ClaSSE. Hence, we had to use the BiSSE model for ancestral state reconstruction, which only allows anagenetic state transitions.
We fit the BiSSE models equivalent to those of the best fitting ClaSSE model and checked that parameter estimates were consistent with the ClaSSE analysis. Then, we used this model to make an ancestral state reconstruction that was compared to the biogeographic reconstruction.

| Time-dependent diversification
The ability of character state-dependent diversification models to confidently identify correlations between diversification rates and character state has been recently challenged (Maddison & FitzJohn, 2015;Rabosky & Goldberg, 2015). In particular, Rabosky and Goldberg (2015) have pointed at the effect of diversification rate heterogeneity in the tree on the probability of detecting false positives. We performed an additional test, which aimed at shedding light on the potential effect of rate heterogeneity on the results of ClaSSE. We tested for variation of speciation/extinction rate according to major radiations having occurring in one biogeographic region. These radiations, if associated with increased diversification rate, may strongly affect the results of ClaSSE because they bear a strong phylogenetic signal. For this purpose, we chose the approach developed by Morlon, Parsons, and Plotkin (2011) to estimate diversification rates through time. This is a maximum-likelihood method, which accommodates time-dependent birth-death processes that can vary across a tree and the positions of the rate shifts have to be specified a priori. Then, AIC comparisons can be used to determine whether a shift significantly improves the fit of the model or not.
For each subtribe, we defined a priori positions of shifts based on the time-tree configuration and the biogeographic reconstruction.
As such, we tested a single shift in Dircennina, a subclade within the genus Pteronymia (hereafter, Pteronymia-group) that rapidly diversi- In each case, we started by fitting models on the whole tree, but refining the sampling fraction to the different subclades for which we tested diversification shifts instead of a global sampling fraction. This provided a null hypothesis to which we compared the fit of models with an increasing number of additional shifts. When we fit models on the different subclades, the stem branch was included in the subclade (as the method was originally designed, Morlon et al., 2011), but the speciation event that led to the subclade was added to the remaining background tree to keep track of this cladogenetic event. The stem branch of the background tree was not included in the analyses. For each tree (whole tree, subclades, and corresponding background trees), we fit six models: constant speciation without extinction, time-dependent speciation without extinction, constant speciation with constant extinction, time-dependent speciation and constant extinction, constant speciation and time-dependent extinction, time-dependent speciation and time-dependent extinction, allowing to take into consideration all possible cases of constant, time-varying or null rates. Time dependency was modeled using an exponential function. We considered the constant speciation/no extinction model as the null model.
Models were compared using AIC scores. Only models with an AIC score greater than two units were considered significantly better than the null model.

| Dircennina phylogenetic tree and taxonomic changes
We implemented the taxonomic changes within the genus Pteronymia Other sympatric species of Episcada show weak or no morphological differences aside from wing pattern, so the lack of strong morphological differences between E. striposis and E. clausina is not unusual within the genus.

| Historical biogeography
We first estimated ancestral biogeographic areas using BioGeoBEARS.
Among the six models fitted on Dircennina, the DEC model performed better than DIVALIKE and BAYAREALIKE and the addition of the parameter j did not improve the fit (Table 1) For Oleriina, the best model was a DEC+j (Table 1)    λ112/λ212: cladogenetic transition rates; μ: extinction rates; df: degree of freedom (number of parameters); logL: log-likelihood; AIC: Akaike information criterion score; ∆AIC: difference between the model and the best fitting model.
Bold values indicate the model with the lowest AIC score and all the models falling into 2-unit AIC interval.

| Trait-dependent diversification
In both Dircennina and Oleriina, multiple models were within a 2-unit AIC interval (Table 2a). For Dircennina, the best fitting model involved two different speciation rates among regions. The Andean speciation rate (λ 222 = 0.230) was nearly twice as high as non-Andean speciation rates (λ 111 = 0.118). The second best fitting model recovered the same pattern but had in addition a twofold higher colonization rate from the Andes toward non-Andean regions (λ 212 = 0.094) than the other way round (λ 112 = 0.047).
MCMC analyses confirmed that the pattern of different speciation rates is strong (Figure 4). The root of Dircennina was clearly inferred to be Andean, unlike in the BioGeoBEARS reconstruction.
In Oleriina, the model with the lowest AIC had different colonization rates but three other models were found within an AIC interval of 2, involving either different speciation rates, different colonization rates or extinction rates (Table 2b). Among these four best models was also the "null" model, in which speciation rates within regions are equal, colonization rates between regions are equal, and speciation rates are different from colonization rates. We performed MCMC analyses for the first three "best" models ( Figure 4). In all cases, the results of the MCMC confirmed that differences between parameters were very small with posterior distributions largely overlapping ( Figure 4) contrasting with the strong difference in speciation rates found in Dircennina. In addition, ancestral state reconstructions performed using the parameters estimated by ClaSSE greatly diverged from the ancestral state reconstruction obtained with BiogeoBEARS (see also the discussion about BiSSE ancestral state reconstruction in Appendix S7). Only the ancestral state reconstruction performed with the "null" model was congruent among BiSSE, ClaSSE, and BiogeoBEARS (Appendix S7). Therefore, this model, which had the lowest number of parameters, was chosen for the ancestral state reconstruction in Figure 3 and discussed in the results. Given the results from the maximum-likelihood analyses, the MCMC, and the ancestral state reconstructions, we conclude that neither the cradle hypothesis, nor the species-attractor hypothesis nor the museum hypothesis is clearly supported in Oleriina and we discuss this result below.

| Time-dependent diversification
For Dircennina, the diversification rate shift in the Pteronymia-group (subclade of the genus Pteronymia) was significant (∆ AIC = 3.6 compared with null model, Table 3a and Appendix S8). For this subclade, the lowest AIC corresponded to a model of time-dependent speciation rate, but it was not significantly different from the null model (constant speciation without extinction, Appendix S8). The speciation rate inferred by the constant speciation rate model was 0.384 ( Figure 5).
In the background, a model of time-dependent speciation and extinction had the lowest AIC score, but it was not significantly different from the null model (Table 3a and Appendix S8). The speciation rate estimated for the background process by the constant speciation rate model was 0.225, that is, lower than that of the Pteronymia-group  Table 3b and Appendix S8). This model shows an initial speciation rate of 1.337, which rapidly decreases to 0.053 at present ( Figure 5). The best model for the onega-group was also a time-dependent speciation rate with no extinction, with an initial rate of 0.821 at the root of the clade followed by a decrease toward 0.082 at present ( Figure 5). The remaining background tree was characterized by a constant speciation rate of 0.155 without extinction, which means that both the makrena-group and the onega-group initially shifted toward much higher speciation rates (Table 3b, Figure 5 and Appendix S8). The accumulation of species reconstructed from the best fitting models showed that the makrena-group and the onega-group had a relatively similar accumulation of species during the last 6.7 million years and 9.6 million years, respectively, with a fast initial phase followed by a slowdown during the last 1-2 million years ( Figure 5).
F I G U R E 4 Posterior probability distribution of parameters obtained from the MCMC analyses performed on the best maximum-likelihood fitting ClaSSE models for Dircennina and Oleriina. As two or three models had a similar explanatory power in the maximum-likelihood analysis, we performed MCMC analyses for each model. The parameter constraints of each model are indicated on the left. First column: speciation rate, second column: extinction rate, third column: colonization rate. When the parameters were allowed to vary among regions, character state is indicated using colors: red: speciation or colonization into the Andes, blue: speciation or colonization out of the Andes  investigated using the same framework as here (Chazot et al., 2016), we highlight the similarities and differences between the diversification patterns of the three subtribes (Table 4).

| Diversification of Dircennina
Dircennina shows a strong asymmetrical spatial distribution of current species diversity, with a large fraction of species occurring in the Andes, mostly the northern Andes. In agreement with this pattern, we found that the diversification history of Dircennina is tightly associated with the Andes. Ancestral state estimation for deep nodes from BioGeoBEARS was unhelpful for deciphering the early history of Dircennina. However, the ancestral state estimation based on the character state speciation and extinction model clearly infers an Andean origin. Dircennina have likely occupied the Andean slopes and diversified in this area since their origin, about 15 (95% HPD: 13.39:17.09) million years ago, thereby supporting the time-for-speciation hypothesis (Figure 2). This pattern will need further investigation within a wider phylogenetic framework to infer more reliably the history of the earlier lineages.
A pattern of strong Andean diversification of Dircennina was supported by the trait-dependent diversification analysis. Speciation rate was at least twice as high in the Andes as in non-Andean regions, which supports the cradle hypothesis. Species richness of Dircennina is higher in the Northern Andes and in the Central Andes. We did not find any significant evidence of declining diversification through time, but we estimated a higher speciation rate for the Pteronymia-group, which largely diversified within the Northern Andes (Figures 2 and 5).
This shift may explain in a large part the difference in speciation rate between Andean and non-Andean regions identified by ClaSSE, but also suggests a difference in diversification rates between the Central and Northern Andes. Extinction rates in both Andean or non-Andean regions were always estimated to 0. This suggests that extinction explains very little of the current pattern of diversity and the museum hypothesis is not supported by our results.
Colonization events were strongly biased from Andean toward non-Andean regions. The ancestral state reconstructions based on trait-dependent diversification models recovered 22 out-of-the-Andes dispersal events but only 1 into-the-Andes colonization event and the second best ClaSSE model inferred a higher rate of colonizations toward non-Andean regions than conversely (although this was weakly supported by the MCMC analysis, Figures 2 and 4). This pattern is consistent with an out-of-the-Andes scenario, whereby Andean lineages feed adjacent areas (not only Amazonia) through dispersal events. Such a pattern has been observed in other ithomiine clades, such as the genera Ithomia and Napeogenes (Elias et al., 2009). In the case of Dircennina, these colonization events, as well as range expansions, reached far into Central America. This is likely due to these butterflies' ability to follow mountain chains in Central America, which allowed Andean lineages adapted to mid-and high altitudes to colonize this region (De-Silva et al., 2017). Conversely, colonizations of lowland regions such as the Amazon basin probably required a larger number of adaptations due to changes in climatic conditions (temperatures, moisture), host-plants, or predators, which may have hindered colonization of such regions. In support to this hypothesis, Chazot et al. (2014) showed that the altitudinal niche of ithomiine butterflies tends to be phylogenetically conserved.
The timing of the emergence of a connection between North and South America and the timing of biotic interchanges are highly controversial (see, e.g., Bacon et al., 2015a,b;Lessios, 2015). For some time, the dominant hypothesis was a very recent emergence of land masses  (Coates & Stallard, 2013;Coates et al., 1992).
However, recent publications, including fossil evidence (Bacon et al., 2015a;Bloch et al., 2016), are challenging this hypothesis and instead support a scenario in which a connection between North and South America emerged at least during the early Miocene (Farris et al., 2011;Montes et al., 2012Montes et al., , 2015. Our biogeographic estimation inferred one colonization that could potentially have happened as early as 9 million years ago but this was very poorly supported. All other dispersal events toward Central America occurred during the last 4 million years. Dircennina is relatively species-rich in the Atlantic Forest region compared to other ithomiine subtribes, but this diversity probably did not originate from local speciation events, except in a very limited number of cases. Similar to Central America, BioGeoBEARS inferred an early colonization of Atlantic Forest at the root of Episcada + Ceratinia, but this is very poorly supported (and rejected by the BiSSE reconstruction). Most Atlantic Forest species are the result of independent colonization events (many of them from Andean ancestral lineages).
This pattern appears to be relatively common, at least in butterflies, as exemplified by the ithomiine subtribe Godyridina (Chazot et al., 2016), the ithomiine genera Ithomia and Napeogenes (Elias et al., 2009)

| Comparing diversification patterns among Dircennina, Oleriina, and Godyridina
The comparison of the three subtribes revealed both common features and differences in their patterns of diversification (  Figure 2) and to a much lesser extent Oleriina (5 out-of-the-Andes dispersal events, 3 into-the-Andes dispersal events, Figure 3) reveal a pattern of Andean diversity contributing to that of adjacent areas.
Dircennina, Oleriina, and Godyridina (Chazot et al., 2016) may have all originated in the Andes ~17-15 million years ago. During most of the early and mid-Miocene, the upper Amazon region was covered by the large Pebas wetland (Hoorn et al., 2010;Wesselingh & Salo, 2006 to the Caribbean Sea (Wesselingh & Salo, 2006), likely received episodic marine incursions. As previously suggested (e.g., Chazot et al., 2016;Hughes et al. 2013, Wesselingh & Salo, 2006, this particular ecosystem, which was drained during the late Miocene (~10-8   The cradle hypothesis has been strongly supported in some groups of plants, and several studies have revealed extremely high speciation rates, such as in the Andean genus Lupinus (Fabaceae) (Hughes & Eastwood, 2006) and the bellflowers (Campanulaceae) (Lagomarsino et al., 2016). The high Andean Páramo ecosystem (2800-4700 m) has been reported to host the fastest speciation rates among Earth's biodiversity hot spots (Madriñán et al. 2014). However, this later study cannot be directly compared to our results supporting the cradle hypothesis in Dircennina and to a lesser extent in Godyridina (Chazot et al., 2016). The fast speciation rates reported in plants mostly occur in the Páramo, a geographically restricted ecosystem of high altitude in which ithomiine butterflies do not occur and which might be highly sensitive to recent Pleistocene climatic fluctuations and the timing of Andean uplift. Although the Andes affect speciation rates in many kinds of organisms, the patterns of diversification differ among lineages. In groups other than plants, the cradle hypothesis has found mixed support. For example, Hutter et al. (2013) did not find higher speciation rates in Andean glassfrog lineages (Centrolenidae), while Beckman and Witt (2015) found higher speciation rates in Andean goldfinches and siskins (Fringillidae). In Heliconius butterflies, Rosser et al. (2012) found that species richness peaked in the eastern slopes of the Andes and was characterized by very "young" species and they interpreted this pattern as a support for the cradle hypothesis. We did not find support for any pattern driven by extinction in the three ithomiine subtribes, and in fact, we found no signal of extinction, which argues against the museum hypothesis. The role of extinction has been poorly addressed perhaps because it has been rarely proposed as an explicit biogeographic scenario, but probably also because of the controversy surrounding the ability of current methods to reliably estimate extinction from molecular phylogenies of extant species (e.g., Rabosky, 2010; but see Morlon et al., 2011); hence examples are rare. For example, Hutter et al. (2013) did not find support for different extinction rates among Andean versus non-Andean regions in glassfrogs, but Antonelli and Sanmartín (2011) reported a lower extinction rate (combined with higher speciation rate) in the species-rich Andean subgenus Tafalla compared to the remaining non-Andean Chloranthaceae. Despite being rarely proposed as an explicit scenario of Andean diversification, the species-attractor hypothesis has been relatively well supported.

| Comparison beyond the Ithomiini
For example, Hall (2005) found repeated speciation events across altitudes as well as colonization events into the Andes in Ithomiola butterflies (Riodinidae). In plants, a large number of independent colonization events of the Andean slopes have been reported, for example, in Begonia (Begoniaceae) (Moonlight et al., 2015) and Bromeliaceae (Givnish et al., 2014). Finally, the time-for-speciation hypothesis, which is supported by our analyses on Dircennina and Oleriina, has also been supported in other groups such as the glassfrogs (Hutter et al., 2013).
In conclusion, we show that a strict evaluation of biogeographic scenarios within a common framework allows a "comparative" biogeographic approach. This approach helps to clearly decipher which processes are shared across different groups as well as why some groups differ from others. Here, we show that, at least for Oleriina and Godyridina, the Central Andean fauna appears to be old with slow diversification rates compared to the Northern Andean fauna, which is more recent and is diversifying at a faster rate. The three subtribes also show major dispersal and diversification events associated with the demise of the Pebas system, with a remarkable convergence in the timing of events. However, major differences also appear between these groups, especially when considering Amazonian diversification. Notably, repeated independent events of rapid diversification in Amazonia question the hypothesis that the Andes acted as a cradle (i.e., driving higher speciation rate) at the global scale and instead call for further investigation on the ecological or genetic characteristics explaining why some groups radiated in Amazonia and others not.

ACKNOWLEDGMENTS
This project was funded by an ATIP (CNRS, France) grant awarded to ME and LDS's postdoc was funded by this grant. ME also acknowledges support from the French National Research Agency grant SPECREP (ANR-14-CE02-0011-01). NC was funded by a doctoral fellowship from Ecole Doctorale 227 (France) and Lund University (Sweden). We thank authorities and counterparts at institutions in Peru, Ecuador, and Brazil (ICMBio) for delivering research and collection permits, as well as many assistants for their help in the field. Molecular work was performed at the GenePool (University of Edinburgh, UK), CBMEG-Unicamp (Brazil), and the Service de Systématique Moléculaire UMS2700 of the MNHN (France). AVLF thanks Luiza M. Magaldi for helping with molecular work and also ac-