The roles of wing color pattern and geography in the evolution of Neotropical Preponini butterflies

Abstract Diversification rates and evolutionary trajectories are known to be influenced by phenotypic traits and the geographic history of the landscapes that organisms inhabit. One of the most conspicuous traits in butterflies is their wing color pattern, which has been shown to be important in speciation. The evolution of many taxa in the Neotropics has also been influenced by major geological events. Using a dated, species‐level molecular phylogenetic hypothesis for Preponini, a colorful Neotropical butterfly tribe, we evaluated whether diversification rates were constant or varied through time, and how they were influenced by color pattern evolution and biogeographical events. We found that Preponini originated approximately 28 million years ago and that diversification has increased through time consistent with major periods of Andean uplift. Even though some clades show evolutionarily rapid transitions in coloration, contrary to our expectations, these shifts were not correlated with shifts in diversification. Involvement in mimicry with other butterfly groups might explain the rapid changes in dorsal color patterns in this tribe, but such changes have not increased species diversification in this group. However, we found evidence for an influence of major Miocene and Pliocene geological events on the tribe's evolution. Preponini apparently originated within South America, and range evolution has since been dynamic, congruent with Andean geologic activity, closure of the Panama Isthmus, and Miocene climate variability.

As an alternative to trait-dependent diversification, geographic events can also promote diversification and range evolution.
Understanding the biogeography of diverse Neotropical clades is therefore also key to obtain insights about the origins of extant biodiversity. There have been many changes in Neotropical landscape's configuration that have influenced the distribution and evolutionary trajectories of organisms (Antonelli et al., 2018;Hoorn et al., 2010). For example, the Andean uplift that started around 30 Ma has been one of the key drivers of speciation and range evolution in the Neotropical biota (Smith et al., 2014). This uplift created major rearrangements of internal wetlands that provided novel habitats and were important barriers to dispersal (Chazot et al., 2019;Rahbek et al., 2019). The formation of the Panama Isthmus was another major biogeographical event that promoted the great biotic American interchange shaping current Neotropical biodiversity (Bacon et al., 2015).
Phylogenetic methods now make possible the investigation of the relative contributions of phenotypic traits and biogeographical events on the evolutionary history of different groups of organisms (Pinto-Sanchez et al., 2012;Weir et al., 2009). Combining knowledge of phylogenetic relationships, the timing of diversification and trait and distribution data permits testing of competing hypotheses about the origin and evolutions of diverse biotas (Mullen et al., 2011), and identifies traits that might have had a crucial role in the evolutionary history of organisms (Losos et al., 1997).
Butterflies are an excellent model system for the study of color evolution and the influence of geographic events on diversification.
Color pattern alone can be used to distinguish most of the 18,000 described butterfly species (Nijhout, 1991). The most studied functions of color include mimicry, predator avoidance, mate recognition, and sexual selection, and these functions are thought to have influenced diversification (Chazot et al., 2014;Jiggins, 2008;Kemp, 2007;Obara & Majerus, 2000). Similarly, novel habitats and barriers to dispersal created by Neotropical landscape rearrangement have influenced the evolution of a number of butterfly groups (Chazot et al., 2016(Chazot et al., , 2019Condamine et al., 2012;De-Silva et al., 2016Toussaint et al., 2019). The increasing availability of comprehensive dated phylogenetic hypotheses for butterflies (Chazot et al., 2020;Espeland et al., 2018) allows a more rigorous study of how different, and potentially conflicting, functions of color have generated trait and species distributions and diversity (Finkbeiner et al., 2014). Phylogenies also enable tests of how shifts in color pattern, in concert with changes in geographic range and habitat, have influenced speciation and diversification (Chazot et al., 2014;Jiggins et al., 2006).
The Neotropical region contains approximately 45% of species in the family Nymphalidae, the most species-rich family of butterflies (Chazot et al., 2020). The Nymphalidae contains such morphologically diverse groups of species that many currently recognized subfamilies were once treated as families. The family ranges from the often small and drab Satyrinae to some of the largest and most spectacularly colored butterflies, of which the Neotropical tribe Preponini (Charaxinae) is renowned as one of the most outstanding examples.
Preponine butterflies inhabit the forest canopy and are characterized by robust bodies and erratic flight. They are distributed from Mexico to Argentina, with a species richness peak in the Amazon basin. Wing color patterns among Preponini species exhibit a dramatic transition from dorsally blue and ventrally brown to dorsally red/orange and ventrally multicolored (Figure 1). This color pattern variation explains why some Prepona species were long classified in a separate genus, Agrias (Ortiz-Acevedo et al., 2017). The bright color patterns of preponines have attracted the attention of collectors and naturalists for over two centuries, resulting in hundreds of names for species within Prepona in particular (Lamas, 2004). Lamas Nevertheless, no study to date has focused on understanding the role of color and biogeography on the diversification of this group.
In this study, we attempt to understand to what extent coloration and biogeographical events may have shaped the diversity of the tribe by quantifying color change and diversification rates, and estimating ancestral geographic ranges. We hypothesize that the main driver of Preponini evolution is color shifts, but biogeographical events in the last 30 Ma also likely influenced the origins of particular clades. Given that the members of Prepona show drastic changes in color patterns, we hypothesize that there might be a congruence between shifts in diversification and phenotypic evolutionary rates.
Similarly, we expected that changes in diversification rates might be associated with major landscape changes, such as the Andes uplift and final closure of the Panama Isthmus.

| Phylogeny and divergence time estimation
We used DNA sequence data from Ortiz-Acevedo et al. (2017). The final matrix comprises sequences from two mitochondrial and four nuclear gene fragments from 29 Preponini specimens including 24 of the 25 generally recognized species and five representatives of geographically and/or morphologically distinct populations of Prepona deiphile and P. laertes that we currently consider as distinct species (LlorenteBousquets, J. E. and E. Ortiz-Acevedo unpublished data; Ortiz-Acevedo, E. unpublished data; Neild, 1996;Turrent Carriles et al., 2019). We used as outgroup the sister tribe Anaeomorphini (Espeland et al., 2018). Out of Preponini and Anaeomorphini, our taxon sample excludes only the recently described Prepona silvana from western Mexico, apparently a close relative of what we refer to as P. deiphile (CA), and P. "sahlkei," treated by Neild (1996) as a distinct species from the Guianas but as synonym of P. claudina claudina by Lamas (2004), and whose status is currently uncertain (Table S1).
We partitioned the data by codon position a priori and then searched for the optimal partitioning scheme using the greedy algorithm in PartitionFinder v2.1.1 (Lanfear et al., 2012) (Table S2).
Phylogenetic relationships and divergence times were inferred simultaneously using BEAST 2.4.5 (Bouckaert et al., 2014). The clock prior was set to an uncorrelated relaxed lognormal distribution, and site models were coestimated with the phylogeny using reversible jump by applying the bModelTest 0.3.3 plugin (Bouckaert and Drummond, 2017). Site models were unlinked, and tree models were linked. We ran three different analyses that differed in the number of molecular clocks used. The models consisted of (a) two molecular clocks, one for all mitochondrial partitions and the other for all nuclear partitions; (b) multiple clocks, one for all mitochondrial partitions and one for each nuclear partitions; and (c) multiple molecular clocks, one for each partition.
The tree prior was set to a Yule model. The birth rate prior was set to uniform with the lower boundary of 0 and upper boundary of 1,000, and the prior for the mean rate under the uncorrelated lognormal relaxed molecular clock (ucldMean) was set to a gamma distribution with alpha of 0.01 and beta of 1,000. We used a secondary calibration point for the common ancestor node for Preponini + Anaeomorphini assuming a normal prior distribution with mean at 32.9 Ma (95% credibility interval CI = 22.0-42.3 Ma) and sigma of 4.0 (Espeland et al., 2018). The value of sigma was selected to include the error associated with the primary dating study and incorporate the credibility interval estimated for the node Preponini + Anaeomorphini in the calibration (Forest, 2009). Three individual runs of 100 million generations were performed, sampling every 10,000 generations. The runs were combined in LogCombiner 2.4.5, and convergence was assessed in Tracer v1.6 . From the combined runs, we subsampled a total of 10,000 trees using LogCombiner to infer a maximum credibility tree in TreeAnnotator 2.4.5.
We used path sampling and stepping-stone sampling to estimate the marginal likelihoods of the different molecular clock analyses (Baele et al., 2013) and used Bayes factors (Fan et al., 2011) to select the best analysis. We used 100 steps with a chain length of one million generations, and other settings were kept as default, and we tested convergence by examining that the average standard deviation of split frequencies fell below 0.01, visually inspecting the trace files in Tracer v1.6 (http://tree.bio.ed.ac.uk/softw are/trace r/), and checking that the effective sample size values for parameters were higher than 200.

| Wing color pattern evolution
To reconstruct ancestral color patterns and estimate their rate of evolution, we photographed museum specimens deposited at the McGuire Center for Lepidoptera and Biodiversity, Florida Museum of Natural History (FLMNH-MGCL) collection. We sampled on average 12 individuals per species, but some species, such as Prepona amydon (n = 61), had higher sampling due to their phenotypic diversity. We reduced possible color pattern variation caused by adverse F I G U R E 1 Example Preponini species. Butterflies shown (from left to right across each row): Archaeoprepona licomedes, A. demophon, A. demophoon, A. amphimachus, A. meander, A. camilla, A. chalciope, A. priene, Mesoprepona pheridamas, Prepona dexamenus, P. laertes, P. werneri, P. gnorima, P. deiphile, P. praeneste, P. narcissus, P. aedon, P. hewitsonius, P. claudina, P. amydon. Photos by P. S. Padrón and A. Warren collection storage conditions by selecting as recently collected specimens as possible ( Figure S1).
Photographs in JPEG format were taken using a light box with a set of four daylight fluorescent Sylvania light bulbs and a Nikon D5300 camera body coupled with a Nikon 60 mm f/2.8G ED Auto Focus-S Micro-Nikkor Lens. We used a Kodak color separation guide and grayscale with a ruler included to calibrate the camera and the images before analysis. Color was measured at three independent locations on the forewing and over the forewing as a whole. Locations were delimited using wing veins, which delineate three homologous regions across species. The same region in the wing was measured consistently, irrespective of wing size or shape. The three forewing locations included the discal cell (Cell 1), and the regions delimited by the discal cell and the veins M 1 and M 2 (Cell 2), and CuA 1 and CuA 2 (Cell 3) ( Figure S2). Those regions were selected to represent regions containing the most significant observed variation in color across the forewing and to control for differences in wing size.
On each of the wing regions and for the whole wing, we measured the mode of the red, green, and blue channels (RGB) and total RGB values using ImageJ (Abramoff, 2004;Rasband, 1997;Schneider et al., 2012). In total, we measured color as a set of 16 traits. We preferred the mode over the mean because in a skewed distribution, it better describes where the bulk of the density is concentrated, while the mean may deviate from the mode because of large numbers in the tails of the distribution. Since we wanted to test for differences in evolutionary rates and the signature of evolution in the three different cells and the overall wing, we did not attempt to reduce the variable set by means of principal component analysis.
Additionally, since Preponini consists of three genera in two clades, one mainly blue and the other blue and red, we wanted to evaluate the rates of evolution of the different color channels independently.
To test for changes in evolutionary rates in color across lineages, we used the function rjmcmc.bm in the package "geiger 2.0" (Harmon et al., 2008) in R, which implements a Bayesian approach using a reversible jump Markov chain Monte Carlo (rjMCMC) process to compare among four different models of changes in evolutionary rates (Eastman et al., 2011). The first was a single rate Brownian motion (BM) model in which phenotypic rate evolution was constant across the tree. Next, we fitted a relaxed BM model (rBM1) in which evolutionary rates were allowed to shift multiple times across trees. Third, we kept evolutionary rate constant (no shift), but traits were allowed to pulse rapidly, representing jumps in the mean of the trait (jump-BM). Last, we combined rBM1 and jump-BM to allow phenotypic rates to shift several times and the mean trait to jump across the tree (jump-rBM, Eastman et al., 2011). Before running the models, we calibrated the rjMCMC and estimated that the most reasonable proposal width to initiate sampling of the Markov chain was eight for every model evaluated (see Eastman et al., 2011).
Model selection was performed using the Akaike information criterion (AIC) for MCMC samples (Raftery et al., 2007). The difference in AIC (ΔAIC) between models was used to select the best model. We assumed a model to be significantly better than another if ΔAIC > 2. In the case in which ΔAIC < 2, we selected the simplest model in the following order: (a) BM, (b) rBM1, (c) jump-BM, and (d) jump-rBM. We did not perform model averaging in cases in which −2 < ΔBIC < 2, because it has been suggested to be misleading (Taper & Ponciano, 2016).
For each model in each cell and measurement, we ran two chains each for 1,250,000 generations sampling every 1,000. We assessed convergence of the MCMC by inspecting the trace of each parameter and estimating the effective sample size. Finally, we used maximum-likelihood ancestral state reconstruction to visualize trait evolution (Felsenstein, 2004) using the ace function in "ape" (Paradis et al., 2004) in R. The shifts in R and B rate on the same branch might be due, however, to more or less rapid fluctuations between red and blue, or alternatively, they might be due to more or less rapid fluctuations in lightness. Since we hypothesize that changes were mainly due to changes in hue from red to blue rather than lightness, to account for the above-mentioned issue, we computed a lightness-independent index by regressing the RGB values against lightness and repeated the analyses using the residuals of this regression ( Figure S3).

| Diversification dynamic estimation
To investigate whether diversification rates varied through time, we used a series of likelihood-based diversification models in the R package RPANDA . One of the models implemented in RPANDA estimates the likelihood that speciation and extinction are constant or variable through time. Accordingly, we tested a combination of models in which we allowed speciation and extinction to be either constant or variable following a linear or an exponential model (see Table 1 for model specification and results for details on linear models). We compared the models using AIC.
In addition to the RPANDA time-dependent diversification analyses, we used alternative methods to support the estimation of diversification rates for the tribe. Overall, we estimated diversification rates using the constant-rates test (CRT; Pybus & Harvey, 2000), Magallón and Sanderson's estimator (Magallon & Sanderson, 2001), and Bayesian analysis of macroevolutionary mixtures (BAMM; Rabosky, 2014) and MEDUSA (Alfaro et al., 2009), and visualized the results with a lineage-throughtime plot (LTT). Furthermore, we used a recently developed method, which estimates unbiased speciation and extinction rates based on LTT (Louca & Doebeli, 2017;Louca & Pennell, 2020).
Details about these diversification analyses are provided in the Appendix S1.

| Trait-dependent diversification analyses
To test whether Preponini lineage diversification has been dependent on phenotypic evolution (i.e., color pattern), we used the equal-splits with simulated null model (ES-sim) method (Harvey & Rabosky, 2018 traits. In addition, ES-sim is useful for small phylogenies with few and slight changes in diversification rate, which might be the case in our study (Harvey & Rabosky, 2018). As a result of the ES-sim test, we obtained (a) a Pearson's correlation coefficient (ρ), which quantifies the strength of the relationship between traits and each lineage's diversification rate; and (b) the simulation-based two tailed probability value (p-value), which allows comparison of the data to a null hypothesis.

| Biogeographical-dependent diversification analyses
To test whether the Andean uplift or closure of the Panama Isthmus influenced Preponini diversification rates, we used a similar likelihood approach as for time-dependent diversification analyses implemented in RPANDA. In this case, we evaluated the fit of six models allowing speciation and extinction rates to be dependent on the paleoelevation of the Andes (three models) or the degree of connectivity between Central and South America (three models; Morlon, 2014; Morlon et al., 2011; see Table 1).
We used Andean paleoelevation from present time until 32 Ma based on the database provided by Lagomarsino et al. (2016) and information from Garzione et al. (2008). The closure of the Panama Isthmus was estimated using reports of migration rates for several groups of organisms by Bacon et al. (2015). In the latter study, they reported that migration rates between the two continents increased sequentially at four points in time, 41.1, 23.7, 8.7, and 5.2 Ma. We used these migration rate estimates to compute the probability of observing at least one migration event in each time period (i.e., 50-41.1, 41.1-23.7, 23.7-8.7, 8.7-5.2, and 5.2-present). We then used the cumulative probability through time to observe at least one migration event per million years as a proxy of the degree of connectivity between land masses (Table S3).

| Estimation of ancestral ranges
We used a Microsoft Access database to compile locality data from 4,050 Preponini specimens deposited at five butterfly collec- We used the R package BioGeoBEARS 1.1.1 (Matzke, 2018) to estimate the ancestral range of Preponini under the dispersal-extinction-cladogenesis (DEC) model (Ree & Smith, 2008) and a maximum-likelihood implementation of the dispersal-vicariance analysis (DIVALIKE) model (Ronquist, 1997). We did not include models with the parameter J since they seem less relevant to continental settings and have been shown to be difficult to compare statistically to other non-nested models (Ree & Sanmartín, 2018). We tested eight different hypotheses to evaluate the influence of different biogeographical events and distance among areas in the evolutionary history of Preponini (see Table 2 for details). Condamine et al., 2012; Montes et al., 2015). We allowed the probability of movement across areas to change in time, accounting for geographic position, for distance and for barriers to dispersal, and penalized accordingly (Tables S4-S11). We evaluated the effect of distance by penalizing dispersal probabilities with a factor proportional to the distance among the centroids of the areas (Tables S5,   S8-S9, S11).
Both DEC and DIVA are biased toward estimating widespread ranges in deep nodes (Buerki et al., 2011;Clark et al., 2008;Matzke, 2014;Ree & Smith, 2008). In an attempt to avoid this potential bias, the maximum number of areas any ancestor may occupy was set to six, since this is the maximum number of areas observed to be currently occupied by any single Preponini species (Ronquist & Sanmartin, 2011;e.g., Archaeoprepona demophoon). We reduced the number of potential ranges in the reconstruction by including only geographic ranges with adjacent areas. This resulted in 113 potential ranges from the 248 possible combinations. To identify the hypotheses with strongest support, we compared AIC among the eight hypotheses and two possible models. The gradual closing of the Panama Isthmus is the only geographic event that influenced Preponini dispersal.

Ho4
Restriction to dispersal by Andes na na yes The gradual uplift of the Andes mountain range is the only geographic event that influenced Preponini dispersal. The Andes uplift generated a reconfiguration of the Amazonian lowland landscape by aiding the formation of the Pebas system and subsequently the Acre system.

Ho5
Restriction to dispersal by distance + Andes yes na yes Preponini dispersal was influenced by a combination of geographic distance and the uplift of the Andes mountain range.

Ho6
Restriction to dispersal by Distance + Panama yes yes na Preponini dispersal was influenced by a combination of geographic distance and the closing of the Panama Isthmus.

Ho7
Restriction to dispersal by Panama + Andes na yes yes Preponini dispersal was influenced by a combination of both the closing of the Panama Isthmus and the uplift of the Andes mountain range.

Ho8
Restriction to dispersal by distance + Andes + Panama yes yes yes Preponini dispersal was influenced by a combination of geographic distance, the closing of the Panama Isthmus, and the uplift of the Andes mountain range.

| Phylogeny and divergence time estimation
The three analyses with different molecular clocks were equally supported (BF = ~1; Tables S12 and S13). The phylogenetic reconstruction is mostly congruent with published topologies and has good support overall (Figures 2 and 3, Figure S4) Figures 2 and 3, Figure S4). The red Prepona clade originated in the early Pliocene, ~5 Ma.

| Wing color pattern evolution
Although the overall signature of evolution suggests constant phenotypic evolutionary rates, we identified two jumps in evolutionary rate in the red channel Cell 1 in the red Prepona clade: the first one at the base of the red Prepona clade and the second along the branch leading to Prepona hewitsonius (Table 3, Figure 2, Figure S5).
Similarly, our analyses identified one jump in the total RGB of Cell 1 at the base of the red Prepona clade. The blue channel in Cell 2 showed 3 jumps (posterior probability > 0.5), located at the base of the P. hewitsonius + P. amydon clade, at the branch leading to P. hewitsonius, and at the tip leading to P. philipponi (W) (Figure 2, Figure S5).
The analysis using the lightness-independent index yielded similar results (Table S14, Figure S6), so we only show and discuss the results using the raw RGB data.

| Diversification dynamic estimation
The RPANDA analyses identified the time-dependent model with constant speciation and no extinction, and with exponentially varying speciation, as the best models (Table 1). MEDUSA (Appendix S1). Moreover, estimates of unbiased speciation and extinction rates were similar to those obtained using the maximum-likelihood approach in RPANDA (Appendix S1). Furthermore, despite recent critiques of approaches for estimating speciation rates, all models that we tested resulted in similar parameter estimates of between 0.12 and 0.17 LogSpecies per Ma (i.e., 1.12 and 1.18 species per Ma).

| Trait-dependent diversification analyses
We found that diversification rate was not dependent on any of the color traits. We found that the absolute value of Pearson's rho  (Table 4).

F I G U R E 2 Ancestral reconstruction of
Furthermore, none of the Pearson's rho was significantly different from the null expectation.

| Biogeographical-dependent diversification analyses
Diversification models accounting for the influence of biogeographical events on the diversification rates of Preponini suggest that speciation and both the Andean uplift and the closure of the Panama Isthmus were exponentially and positively correlated through time (Table 1). In fact, the model with exponentially varying speciation rate depending on Andean uplift had a slightly lower AIC than the time-dependent models (Table 1).

| Estimation of ancestral ranges
BioGeoBEARS analyses identified DEC as the best model with restriction to dispersal and an influence of the closure of the Panama Isthmus and Andean uplift (Ho7) as the hypothesis with strongest support ( Figure 3, Table 5). The tribe most likely originated in South America (node I in Figure 3, Figure S7). The genus Archaeoprepona retained its ancestor's range of origin (node II in Figure 3); other ranges with high probability also suggest a continental America origin ( Figure S7).

Similarly, Mesoprepona + Prepona likely originated in South America
(node III; Figure 3); the other probable regions of origin comprise broad ranges as for its sister genus ( Figure S7). The origin of Prepona was most likely in mainland South America, excluding high-elevation habitats of the western Andes (node IV; Figure 3). Finally, the clade that corresponds to the transition of color patterns from mostly blue to red Prepona species was found to have most likely originated in the Amazon and Eastern Andes (node V; Figure 3, Figure S7); the other ranges with high probability are narrow, where Andes, Amazon, and Chocó are also likely regions of origin ( Figure S7).

| D ISCUSS I ON
In our study, we hypothesized that shifts in color patterns and/or biogeographical events have contributed to Preponini diversification.

F I G U R E 3 Biogeographical range estimates for the tribe Preponini based on the DEC model. Dotted gray lines delimit the time slices used
in the stratified analysis. Upward black triangles denote periods of increased Andean uplift (Hoorn et al., 2010), downward black triangles denote periods of increased biological migration between Central and South America, and downward white triangle denotes initial faunal exchange between Central and South America at ~41 Ma (Bacon et al., 2015). Roman numerals in nodes refer to the main text.  (Figure 2; node V in Figure 3), we found that diversification rates in Preponini were not dependent on color evolution.
However, major landscape reconfiguration, happening through most of the tribe's evolutionary history, did seem to influence speciation rates and distribution patterns (

| Wing color pattern evolution
During the course of their evolution, Prepona underwent drastic changes in the color pattern of both wing surfaces (Figures 2 and 3).
These color patterns were the key features used by previous taxonomists to classify Prepona species in two genera. The former genus Agrias contained yellow, orange, blue, and red species, and Prepona, as delimited previously, contained mostly blue species. The former group being nested within the latter is supported by molecular and morphological data (Ortiz-Acevedo et al., 2017), while color patterns apparently transitioned rapidly from blue to red (Figure 2), resulting in Prepona containing butterflies with such strikingly different coloration patterns (Figures 1-3).
Our results suggest that there is stronger change in phenotype for red and blue RGB channels (Table 3) (Nadeau et al., 2016;Reed et al., 2011). Differential genetic expression in particular regions of the wing is also responsible for localized changes (Brakefield et al., 1996;Nadeau et al., 2016;Oliver et al., 2012), which might be a plausible explanation for the results we found in Cell 1. These localized changes are not always visible under analysis of the entire wing; thus, isolating those regions that are likely to change quickly allowed us to detect the jumps in phenotype.
Jumps and shifts in the blue channel are restricted to the tips of the phylogeny and might be associated with recent speciation events.
For example, we identified two jumps in the blue channel of Cell 2 ( Figure 2): (a) in the branch leading to Prepona amydon + P. hewitsonius Note: α is defined as the estimated value of the character at the root and σ 2 as the evolutionary rate (or variance).

TA B L E 3
Evolutionary rates for the mode color of the different areas of the wing measured and (b) in P. hewitsonius. These species are phenotypically extremely diverse in forewing coloration, ranging from dark, deep blue to light, bright blue, even within subspecies ( Figure S8). We also identified a change in the Cell 2 blue channel at the branch leading to Prepona philipponi (W). Prepona philipponi (W) has a different light blue tone compared with P. philipponi (E) ( Figure S9), but with only a single studied individual of the former taxon, this result needs corroborating with further samples.
A potential explanation for the dramatic change (jump) in color patterns within Preponini is involvement in mimicry rings with genera such as Callicore and Asterope (Nymphalidae, Biblidinae), which have remarkably similar color patterns on both wing surfaces (Descimon, 1977;Jenkins, 1987; Figure S10). Preliminary geographic distribution data show a correspondence in coloration and distribution among Prepona, Callicore, and Asterope (Ortiz-Acevedo, E. unpublished data). The inferred origin of the red clade in Prepona partially supports this hypothesis, since both Callicore and Asterope are predominantly Andean and Amazonian groups. In particular, the red Prepona showed no change in the Amazonian ancestral range as it underwent speciation ( Figure 3). Classical Müllerian and Batesian mimicry is plausible since species in both groups feed on toxic plants, including the families Sapindaceae and Erythroxylaceae, but whether they sequester plant toxins as caterpillars and retain them after emerging as adults is still unknown, as is their palatability. Alternatively, escape mimicry is also a possible explanation, since all species involved are fast flyers (Pinheiro & Freitas, 2014). Fast and erratic flying has been suggested as an antipredator defense mechanism, where predators learn to avoid these species due to the high-cost-low-benefit trade-off (Ruxton et al., 2004;van Someren & Jackson, 1959).

| Diversification of preponines
In studies of other Neotropical butterfly groups, natural history traits and paleoclimatic events have been shown to shape diversification rates by influencing speciation and/or diversification (Chazot et al., 2020;De-Silva et al., 2016;Sahoo et al., 2017). This seems to be the case for Preponini as well. We found support for exponentially varying speciation rates, similar to results in Neotropical Nymphalidae (Chazot et al., 2020). Extinction in the three top candidate models was a priori fixed to zero, as estimated for larger Neotropical butterfly clades during the same geological time (Chazot et al., 2020).

TA B L E 4 Trait-dependent diversification analysis for the color traits
The current diversification rate estimated by our likelihood analysis (0.15) is within the range of diversification rates estimated for other butterfly groups (Chazot et al., 2020;Peña & Espeland, 2015). We acknowledge, however, that the models with constant speciation have slightly weaker support than exponentially increasing speciation rate models. We thus cautiously interpret our results, favoring the model with the lowest AIC.
Our analyses are most consistent with an exponential increase in speciation promoted by Andean uplift (Table 1). In fact, this bio- in which the Pebas system started to drain and the Amazon river reached its current configuration, aided by the formation of the Acre system (Hoorn et al., 2010).
A recent study showed that phylogenetic trees alone do not have enough information to tease apart different evolutionary scenarios since a particular topology can either result from high speciation and constant extinction or decreased extinction and constant speciation (Louca & Pennell, 2020). Consequently, we would expect similar support for these alternative scenarios. Here we show, however, that support for models with constant speciation and no extinction or exponentially varying speciation with no extinction is much stronger than for models with constant speciation and varying extinction under an evidential statistics framework (Taper & Ponciano, 2016).
In a recent comment to Louca and Pennell (2020), Morlon et al. (2020) advocated for the kind of hypothesis-driven approach that we adopt here, suggesting that the conclusion that speciation and extinction rates are not identifiable parameters does not apply to all cases. Morlon et al. (2020) point out that the fact that there are many congruent diversification models for the same tree topology does not pose a direct challenge to the hypothesis-driven approach used here.

| Origin and biogeographical patterns
Although the high Amazonian species richness of the tribe would suggest an out-of-the-Amazon biogeographical model, the tribe Preponini originated from a widespread South American ancestor at ~27.5 Ma in the early Oligocene. By this time, South America, dominated by forests (Strömberg et al., 2013), was detached from Antarctica and moving north toward North America (Axelrod et al., 1991). This dating is consistent with other studies that used independent datasets (Peña & Wahlberg, 2008;Wahlberg et al., 2009).

Members of Preponini likely dispersed, colonized, and diversified in
new niches in Central America as the bridge between land masses became more continuous at ~23 Ma (Bacon et al., 2015;Montes et al., 2015). This event has been demonstrated to have played a major role in the diversification of multiple groups of organisms (Bacon et al., 2015).
It has been suggested that some of methods for ancestral range reconstruction are biased toward estimating widespread ancestral ranges (Buerki et al., 2011;Clark et al., 2008;Matzke, 2014;Ree & Smith, 2008), as we found in the ancestor of Preponini. To reduce this bias, we restricted the ancestral estimate to contain only six of the eight areas considered, as observed in the most widespread extant Preponini species (Ronquist & Sanmartin, 2011, i.e., Archaeoprepona amphimachus, A. demophoon, and A. demophon). In addition, climatic changes accompanying the evolution of Preponini make biological sense in light of our inferred ancestral range estimate, as we discuss below.
Early speciation events in the tribe are characterized by range maintenance (~13 Ma; nodes II and III; Figure 3). Subsequently, the ancestor of Archaeoprepona contracted in range, becoming restricted to the Atlantic forest region and subsequently dispersing to the Chocó.
The transition from Atlantic Forest to Chocó happened through an anagenetic expansion, including all South American areas except Western Andes, and a subsequent range contraction (Table S15).
The restriction of the ancestor of Archaeoprepona to Atlantic forest is coincident with a period in which tropical forests in South America were likely reduced by a decrease in global temperature (Kürschner et al., 2008;Pound et al., 2011). The late Miocene was probably characterized by grassland habitats, as suggested by high diversification of hoofed animals and changes in their dentition (Kürschner et al., 2008).
Preponini species are not currently known to inhabit these grassland habitats. The subsequent dispersal to the Chocó follows an increase in global temperature and the reappearance of widespread tropical forests (Kürschner et al., 2008;Pound et al., 2011).

After the colonization of Chocó tropical forests,
Archaeopreopona dispersed to central America at ~8 Ma during periods of intensified biological migration between Central and South America (Bacon et al., 2015). ing South and Central America, around a time of increased migration between these two continents (Bacon et al., 2015). The early branching Prepona laertes clade was found to be initially restricted to eastern South America, with subsequent dynamism in range contraction/expansion as it underwent speciation events.
More recent clades were found to have contracted their ranges at ~7 Ma, becoming restricted to the Chocó region, while the Amazon region suffered landscape reconfiguration resulting from the transition of the Pebas system to the Acre system.   (Toussaint et al., 2019).
Unfortunately, knowledge of host plants in Preponini is still rather incomplete, although, researchers, collectors, and enthusiasts are working together to fill this gap in knowledge. with Callicore and Asterope, which also show a high richness in the Amazon region, with a number of species restricted to this area.

ACK N OWLED G M ENTS
We are deeply grateful to everyone that was involved in the field and who donated material to be included in this study. We thank P.
S. Padrón and A. Warren for providing some photographs. We thank the collection managers in the institutions we visited for allowing us to access the specimens. We thank Dr. Lei Xiao for her invaluable help in the molecular laboratory. We are also thankful to our funding

CO N FLI C T O F I NTE R E S T
We declare no conflict of interests.