Introduction
Some phenotypic traits may serve as key innovations that facilitate rapid diversification, influencing the evolutionary origins and trajectories of lineages (Mayhew, 2007; Nicholson et al., 2014; Schluter, 2000). One of the most visually spectacular phenotypic traits are the bright color patterns in several insect groups, and a number of studies have focused on the evolution of these color patterns and the possible roles they might have played in the evolution of insects (Berthier, 2005; Jiggins et al., 2001; Kemp, 2007; Mallet & Gilbert, 1995). In butterflies, for example, wing color pattern is a trait that has attracted the attention of naturalists since Darwin’s and Wallace’s observations in the 19th century (Darwin, 1880, 1874; Wallace, 1889, 1877), and it has proven to be key in intra- and interspecific interactions and speciation processes in butterflies (e.g. Lycaenidae; Fordyce et al., 2002; Heliconius ; Jiggins et al., 2006; Heliconius ; Mavarez et al., 2006).
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 biogeographic 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 biogeographic 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., 2019, 2016; Condamine et al., 2012; De-Silva et al., 2017, 2016; Toussaint 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. Preponini 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 withinPrepona in particular (Lamas, 2004). Lamas (2004) recognized 22 species in the tribe, but subsequent additions and revisions has resulted 25 described species are now recognized, distributed in three genera Archaeoprepona, Mesoprepona and Prepona(Ortiz-Acevedo et al., 2017; Turrent Carriles and García Días, 2019). 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, diversification rates and estimating ancestral geographical 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.
Materials and Methods
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. laertesthat we currently consider as distinct species (Neild, 1996; Turrent Carriles and García Días, 2019, Llorente et al. in prep, Ortiz-Acevedo et al. in prep.). 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 co-estimated 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 i) two molecular clocks, one for all mitochondrial partitions and the other for all nuclear partitions ii) multiple clocks, one for all mitochondrial partitions and one for each nuclear partitions iii) 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 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 incorporating 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 (Rambaut et al., 2014). 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 1 million generations and other settings were kept as default, and 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/software/tracer/) 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 (MGCL-FLMNH) 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 collection storage conditions by selecting as recently collected specimens as possible (Figure S1).
Photographs in JPEG format were taken using a standardized light box setup. Color was measured at three independent locations on the forewing as well as over the forewing as a whole. The three forewing locations included the discal cell (Cell 1), and the regions delimited by the discal cell and the veins M1 and M2(Cell 2), and CuA1 and CuA2 (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 mean and mode of the red, green, and blue channels (RGB) as well as total RGB values using ImageJ (Abramoff, 2004; Rasband, 1997; Schneider et al., 2012). In total, we measured color as a set 32 traits, 16 based on the mean and 16 on the mode. We considered both, the mean and the mode, because the mode in a skewed distribution describes better 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 as well as 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 rates 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 both mean and mode and 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: 1) BM, 2) rBM1, 3) jump-BM and 4) jump-rBM. We did not perform model averaging in cases in which -2 < ΔBIC < 2, because it has been suggested to be misleading (Taper and Ponciano, 2016).
For each model in each cell and measurement, we ran two chains each for 1250000 generations sampling every 1000. 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.