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.