North Andean origin and diversification of the largest ithomiine butterfly genus.

The Neotropics harbour the most diverse flora and fauna on Earth. The Andes are a major centre of diversification and source of diversity for adjacent areas in plants and vertebrates, but studies on insects remain scarce, even though they constitute the largest fraction of terrestrial biodiversity. Here, we combine molecular and morphological characters to generate a dated phylogeny of the butterfly genus Pteronymia (Nymphalidae: Danainae), which we use to infer spatial, elevational and temporal diversification patterns. We first propose six taxonomic changes that raise the generic species total to 53, making Pteronymia the most diverse genus of the tribe Ithomiini. Our biogeographic reconstruction shows that Pteronymia originated in the Northern Andes, where it diversified extensively. Some lineages colonized lowlands and adjacent montane areas, but diversification in those areas remained scarce. The recent colonization of lowland areas was reflected by an increase in the rate of evolution of species' elevational ranges towards present. By contrast, speciation rate decelerated with time, with no extinction. The geological history of the Andes and adjacent regions have likely contributed to Pteronymia diversification by providing compartmentalized habitats and an array of biotic and abiotic conditions, and by limiting dispersal between some areas while promoting interchange across others.

The Neotropical region is the most biologically diverse area on Earth for most organisms and numerous studies have identified the world’s longest mountain range, the Andes, as a major centre of biodiversity1 and source for adjacent areas in groups as diverse as birds2, reptiles3, insects456 and plants78. The Andes have been proposed as a major driver of diversification9. For instance, studies of some Andean plants have found some of the fastest diversification events reported, such as in the Andean Bellflowers, whose 550 species arose in the last 5 million years8. The Andes may have affected diversification rates in different ways, by offering scope for vicariant speciation due to the complex and intricate topography of the mountains, as well as by providing a large array of new environmental niches, thereby promoting adaptive speciation.

The timing of diversification within the Neotropical region is keenly debated and is linked to competing hypotheses about which biogeographic events have primarily driven speciation, extinction and dispersal events during the Cenozoic. However, the rate and geographic extent of surface uplift in the Andes is contentious, having varied through time and among different cordilleras1011, and assessing the role of the Andes and the different phases of uplift on the timing of diversification is complex. The closure of the Panamanian Isthmus during the last 5 million years may have allowed substantial biotic interchange between Central and South America (but see ref. 12) or the Pebas system, a large network of shallow lakes and wetlands, which occupied the upper Amazon region, may have constrained dispersal and promoted local diversification until its drainage around 10-7 million years ago1314. It is also suggested that climatic instability during the Pleistocene drove the diversification of extant Neotropical species, but the importance of this mechanism remains controversial9151617181920.

Insects represent the largest fraction of terrestrial biodiversity but analyses of insect diversification remain scarce compared to other groups such as vertebrates and plants. Butterflies are one of the best studied insect groups and around 7800 Neotropical species have so far been documented21. Over the last decade, an increasing number of studies have used molecular phylogenetic trees to investigate the timing and mode of diversification in a variety of butterfly groups. Such works provide scope for comparative approaches that can help decipher whether common drivers or distinct processes have shaped patterns of diversification. They also allow us to assess the extent to which inferred patterns of diversification match those found in other organisms, particularly well-studied plants222324 and birds225. Many butterfly studies indicate an important role for the Andes in the origin and diversification of new species. In some groups, diversification occurs mostly within the same elevational range, consistent with adaptation onto new resources (e.g., Hypanartia26, Lymanopoda27, some ithomiine butterflies628), whereas others show speciation across the elevational gradient (e.g., Ithomiola29). Even within some predominantly lowland groups the Andes have apparently played an important role, causing diversification contemporaneous with the Andean uplift and consequent major changes in climate and geography in the Neotropical region (e.g., Taygetis30, Dione31, Heliconius32). By contrast, the Andes seem to have had a limited impact on the diversification of some other groups, such as neotropical Troidini33.

The Neotropical butterfly tribe Ithomiini (Nymphalidae: Danainae) is a diverse group with ca. 48 genera and 390 species, which are ubiquitous in humid forest throughout the Neotropical region from sea level to around 3000 m elevation. Commonly known as the clearwing butterflies because of the transparent wings in many species, they are well-known because of their involvement in Müllerian mimicry rings, whereby co-occurring unpalatable species converge in wing colour pattern amongst themselves, other Lepidoptera and some other insects34. Their diversity and broad distribution makes them a relevant study system to investigate patterns of spatial and temporal diversification in the Neotropics and they have previously been the focus of a number of diversification studies463536. Most ithomiine clades have been found to show a peak of species richness in the Andes37, but the biogeographic histories that have led to this pattern are surprisingly diverse. The genus Napeogenes originated in the Andes and subsequently dispersed out of the mountains into the Amazon Basin6, whereas rates of colonization into the Andes from adjacent areas were found to be higher in the subtribe Godyridina35. The genus Oleria contains two main subclades, one of which diversified mainly in lowland Amazonian forests while the other diversified almost exclusively in high elevation Andean cloud forests4.

The genus Pteronymia Butler & Druce, 1872, which belongs to the largest ithomiine subtribe, the Dircennina, is one of the most species-rich ithomiine genera. The genus contains some 50 species (Lamas, 2004), although the species taxonomy has been highly confused in the past and is undergoing revision. Pteronymia butterflies occur throughout the Neotropics, with the most diverse communities found in east Andean cloud forests28. As part of our on-going effort to document patterns of spatial and temporal diversification of butterflies in the Neotropics, here we reconstruct a comprehensive, time calibrated phylogeny for Pteronymia using multi-locus molecular data and morphological characters. We first assess whether Pteronymia is monophyletic and whether species boundaries require re-definition. We then time-calibrate the phylogeny using a combination of larval host plant ages (Solanaceae) as maximal age constraints38 and estimates of divergence time between nymphalid butterfly genera from Wahlberg et al.39 to investigate biogeographical patterns of diversification of Pteronymia in relation with the Andean uplift. Specifically, we addressed the following questions: (1) Have the Andes acted as a centre of origin and a centre of diversification for the genus Pteronymia? (2) Have there been many interchanges between the different regions of the Neotropics, particularly between the Andes and other regions, and between the central and northern Andes? (3) How has the elevational niche of Pteronymia species evolved through time, and is there evidence for adaptive radiation across elevation ranges? (4) How has Pteronymia diversified through time?

Results

We obtained sequence data for a total of 166 Pteronymia specimens, representing 41 of the species recognized prior to our revision, and 47 of the species recognized after our revision (Table 1; see Supplementary Table S1). Species with no molecular data were P. alcmena, P. alicia, P. calgiria, P. fumida, P. glauca and P. peteri. A total of 87 morphological characters were coded for 52 species (after taxonomic revision, see below).

Table 1

List of the species in the genus Pteronymia, including revised status.

Pteronymia species	Previous name	Molecular data	Morphological data	Biogeographical area	Mean elevation [95% range] (m)
Pteronymia alcmena			x	A	1100 [1050–1200]
Pteronymia aletta		x	x	ABCE	693 [50–1450]
Pteronymia alicia			x	I	1235 [1100–1490]
Pteronymia alida		x	x	CDE	1933 [800–3000]
Pteronymia andreas (rev. stat.)	alissa	x	x	CDI	1337 [750–1860]
Pteronymia alissa		x	x	CE	1227 [300–2000]
Pteronymia artena		x	x	ACDE	1365 [900–2060]
Pteronymia calgiria			x	D	1719 [1300–2970]
Pteronymia carlia		x	x	H	910 [100–1700]
Pteronymia cotytto		x	x	AE	776 [80–1455]
Pteronymia donella		x	x	AB	386 [20–825]
Pteronymia euritea		x	x	H	553 [30–1250]
Pteronymia forsteri		x	x	F	632 [182–1350]
Pteronymia fulvimargo		x	x	A	1455 [1110–2060]
Pteronymia fumida			x	AC	870 [450–1000]
Pteronymia gertschi		x	x	C	950 [500–1000]
Pteronymia glauca			x	C	1230 [800–1400]
Pteronymia granica		x	x	C	1929 [1500–2250]
Pteronymia hara		x	x	ACDE	1474 [600–2400]
Pteronymia inania		x	x	CE	1659 [1200–2100]
Pteronymia latilla		x	x	ACE	1012 [170–1925]
Pteronymia laura		x	x	CE	1088 [800–1400]
Pteronymia lonera		x	x	A	1198 [500–1575]
Pteronymia medellina		x	x	C	1707 [1380–2000]
Pteronymia obscuratus		x	x	AB	257 [30–700]
Pteronymia olimba		x	x	D	1403 [1000–1525]
Pteronymia oneida		x	x	CE	1992 [1250–2650]
Pteronymia asopo (rev. stat.)	oneida	x	x	E	900 [900–900]
Pteronymia ozia		x	x	CDE	1446 [600–2500]
Pteronymia parva		x	x	A	842 [200–2060]
Pteronymia peteri			x	I	1040 [1040–1040]
Pteronymia picta		x	x	ACE	1104 [350–2060]
Pteronymia dispar (rev. stat.)	picta	x		C	1500 [1300–1800]
Pteronymia primula		x	x	BFG	389 [65–1050]
Pteronymia rufocincta		x	x	A	998 [210–1500]
Pteronymia sao		x	x	FG	332 [100–800]
Pteronymia serrata		x	x	CD	2091 [1700–2300]
Pteronymia sexpunctata		x	x	D	1265 [635–1500]
Pteronymia simplex		x	x	A	1544 [1100–2300]
Pteronymia tamina		x	x	CD	1667 [800–2075]
Pteronymia thabena (rev. stat.)	teresita	x	x	CDE	1670 [790–2400]
Pteronymia teresita		x	x	C	1184 [1000–1350]
Pteronymia ticida		x	x	CDE	1784 [1200–2400]
Pteronymia tucuna		x	x	F	777 [100–1350]
Pteronymia veia_EAST	veia	x	x	CE	1662 [1000–2400]
Pteronymia veia_WEST	veia	x	x	C	1729 [1300–2150]
Pteronymia vestilla		x	x	F	398 [100–900]
Pteronymia alina (rev. stat.)	zerlina	x	x	CD	1784 [700–2500]
Pteronymia zerlina		x	x	CE	1473 [635–2000]
Pteronymia spnov 1		x	x	E	1730 [1730–1730]
Pteronymia spnov 2		x	x	C	1765 [850–2530]
Pteronymia spnov 3		x	x	C	1898 [1700–2030]
Pteronymia spnov 4		x	x	A	656 [500–700]

Availability of molecular and morphological data, distribution area and elevation mean and 95% range are reported for each species.

Taxonomy

Molecular phylogenies of all Pteronymia specimens generated by maximum likelihood and Bayesian inference were largely congruent (see Supplementary Figs S1–S2), and showed that the genus Pteronymia is monophyletic. Several internal nodes had moderate to low support.

The molecular data suggested that several changes to the species-level taxonomy were warranted. Four cases concern species occurring on both slopes of the Andes, where the molecular data suggest that east and west Andean subspecies are not sister taxa, instead grouping with other related species. In all cases there are no genitalic characters that exclusively support the former classification, which was instead based on attempts to group similar allopatric phenotypes as subspecies of more widespread species21. We therefore split each of the original four species into east and west Andean species.

Pteronymia zerlina

(Hewitson, 1856)

This species formerly included taxa from both east and west of the Andes2140, ranging from Venezuela to western Ecuador and Bolivia. Sequenced material comes from eastern and western Ecuador, representing the taxa formerly known as P. zerlina pronuba (Hewitson, 1870) (west) and P. zerlina machay T. & L. Racheli, 2003 (east). Pteronymia zerlina zerlina (Hewitson, 1856) was described from ‘New Granada’ and the original illustrations and type material in the BMNH match specimens from the Cauca valley in west Colombia. Their smaller size and broad white forewing translucent band are similar to P. zerlina pronuba and it thus seems likely that these two taxa are conspecific. In the eastern Andes, populations throughout Peru apparently link east Ecuadorian P. zerlina machay with the Bolivian P. zerlina alina Haensch, 1909, and thus we treat P. alina as a distinct species (rev. stat.) and transfer to it the following Peruvian and Ecuadorian taxa: P. a. machay, P. a. mielkei Lamas, 2003 (rev. stat.). The status of Venezuelan and central Colombian populations is currently unknown, so for the moment we retain them in P. zerlina. The molecular results are consistent with significant differences in the immature stages of P. zerlina zerlina and P. zerlina machay (now P. alina machay) described by Bolaños et al.41, who also suggested the likelihood that these taxa were not conspecific.

Pteronymia veia

(Hewitson, 1853)

This species formerly included taxa from both east and west of the Andes2140, ranging from Venezuela to western Ecuador and northeastern Peru. Sequenced taxa include an undescribed taxon from western Ecuador and P. veia linzera Herrich-Schäffer, 1865 from eastern Ecuador. As with P. zerlina, the status of Venezuelan and Colombian populations is unknown, but unfortunately there are no described west Colombian taxa which can be reliably associated with the undescribed west Ecuadorian taxon. We treated east and west Ecuadorian taxa as distinct species, but for the moment do not make any nomenclatural changes.

Pteronymia alissa

(Hewitson, 1869)

This species formerly included taxa from both east and west of the Andes2140, ranging from Venezuela to western Ecuador and Bolivia. Sequenced taxa include the nominate subspecies from western Ecuador, and P. alissa andreas Weeks, 1901 from eastern Ecuador. The status of Venezuelan taxa is currently uncertain, so for the moment we retain them in P. alissa and just separate the east Andean P. andreas (rev. stat.) as a separate species, except for the subspecies dorothyae, which clusters with P. andreas in the phylogenetic trees (see Supplementary Figs S1–S2).

Pteronymia teresita

(Hewitson, 1863)

This species formerly included taxa from both east and west of the Andes2140, occurring in western Ecuador and from eastern Colombia to Bolivia. Sequenced taxa include west Ecuadorian P. teresita teresita and east Ecuadorian P. teresita thabena (Hewitson, 1869). We here separate out east Andean populations as a separate species, which both share a distinctive female with colorless translucent hindwing and yellow translucent forewing, including the following: P. thabena thabena, P. thabena denticulata Haensch, 1905 (rev. stat.).

Pteronymia oneida

(Hewitson, 1855)

This species formerly included taxa from western Colombia to Venezuela and along the eastern Andes to northern Peru2140. Pteronymia oneida asopo (C. & R. Felder, 1865) occurs in the Cordillera de la Costa in northern Venezuela, and appeared distantly related to east Ecuadorian P. oneida oneida in the molecular tree. As with related species (P. zerlina, P. veia), there are no genitalic characters that group asopo with oneida, and we therefore treat it as a distinct species, P. asopo (rev. stat.). The status of several other recently described40 Venezuelan taxa of P. oneida remains to be determined, and for the moment they are retained in P. oneida.

Pteronymia picta

(Salvin, 1869)

This species formerly included taxa ranging from Costa Rica to central and western Colombia. Specimen locality data from Colombia are too imprecise and unreliable to confirm whether the Colombian P. picta dispar Haensch, 1905, apparently restricted to the northern Cordilleras Occidental and Central, is sympatric or not with P. picta picta, the range of which appears to broadly encompass that of dispar in Colombia. The molecular data suggest that P. dispar (rev. stat.) is not closely related to P. picta picta and Central American P. picta notilla Butler & H. Druce, 1872.

After our taxonomic changes the genus Pteronymia now comprises 53 species, making it the most diverse ithomiine genus.

Morphological phylogeny

We performed a cladistic analysis of 87 adult and larval morphological traits, which resulted in a relatively poorly resolved tree, within which only several clades received moderate to strong support (see Supplementary Fig. S3). Notable clades include one containing six typically rare Andean species (alida clade: P. alida, P. sp. nov. 3, P. inania, P. lonera, P. teresita and P. thabena), supported by a number of genitalic characters. The immature stages (i. e., larvae and pupae) of this clade are also remarkably different in coloration and morphology from those of other Pteronymia species, and the distinctiveness of the genitalia and immature stages previously led to the description of a new genus, Talamancana, to include P. lonera42. A second, large and apparently well-defined clade contains P. zerlina and relatives, all of which have a very distinctive synapomorphy, a keel-like spine on the dorsal side of the aedeagus near the posterior tip (character 1:1). Genitalia barely differ among any of the eleven species in this clade, although the adult mimetic wing patterns and the immature stage morphology and biology show striking differences.

New calibration of the Solanaceae phylogeny

To calibrate the phylogeny of Pteronymia, we used a combination of secondary calibrations of Nymphalidae ages39 and maximum age constraints based on host-plant lineage ages (Solanaceae). To extract those ages, we used the molecular matrix of a previous Solanaceae phylogeny43 to generate a new phylogeny, which was calibrated with the stem age of the family extracted from a dated phylogeny of Angiosperms44. Our Solanaceae phylogeny shows similar node support and topology to those of the latest published phylogeny of Solanaceae43. Median lineage ages are on average about 25% older and have wider 95% credibility intervals (Table 2, see Supplementary Fig. S4). In most lineages the 95% credibility intervals presented here span the median age of the published phylogeny of Solanaceae43, but the median ages themselves fall almost always outside the 95% credibility intervals of the previous phylogeny (Table 2). The wider credibility intervals in our study are due to the use of the more conservative uniform prior on the calibration point as compared to the previous phylogeny of Solanaceae, which implemented a log-normal prior that tends to drive ages towards the mode of the prior distribution43. Consequently, the ages of the lineages used for calibrating the phylogeny of Pteronymia are older than those inferred previously43 (Table 2).

Table 2

Comparison of our estimates of ages of Solanaceae lineages with those of Särkinen et al. 43.

Solanaceae clades	Age estimates This study	Age estimates Särkinen et al. (2013)
Stem of Solanaceae*	64.0 my [48–83]	49my [46–54]
Crown of Solanaceae	42.2 my [26.8–55.2]	30.3 my [26.3–34.0]
Crown of “x = 12” clade (MRCA Solanum and Nicotiana)	39.5 my [21.9–44.4]	23.7 my [23.0–25.7]
Crown of Solanum	20.9 my [14.5–29.5]	15.6 my [13.1–17.5]
Stem of Capsicum*	19.0 my [11.5–25]	13 my [10–16]
Stem of Schultesianthus*	17.0 my [7–24]	13 my [9–19]
Stem of Brugmansia*	10.5 my [5–14]	7 my [4–10]

*Ages used as calibration in the Pteronymia phylogeny.

Dated combined Pteronymia phylogeny

We combined morphological and molecular data to generate a species-level phylogeny that was calibrated using secondary calibrations from a published Nymphalidae phylogeny39, and from Solanaceae lineage ages estimated in this study, using BEAST 1.7.545. This resulted in 901 trees (after burnin), from which we extracted the maximum clade credibility tree with median branch length (hereafter MCC tree). The combination of morphological and molecular data generated a phylogeny comprising all known extant species (Fig. 1). Many nodes were poorly supported, and this was mostly caused by the species represented only by morphological characters, whose placement was uncertain. Our calibration strategy that combined butterfly- and host-plant-derived secondary calibrations resulted in ages that were about 30 to 50% older than those inferred in a recent time-calibrated phylogeny of Ithomiini genera that relied on previous minimum age estimate of Solanaceae lineages38. By contrast, our ages were often slightly younger, but well within the credibility interval of the ages estimated in the Nymphalidae phylogeny39. Notably, the stem age of Pteronymia inferred in our study is 14.4 million years (my) [12.3–16.3], while it was 7.5 my [6.0–9.0] in the higher level Ithomiini phylogeny38 and 15.7 my [11.5–18.5] in the Nymphalidae phylogeny39. We repeated this analysis without the Solanaceae calibrations and this had little impact on most nodes of the phylogeny. The greatest difference was found for the divergence between the outgroup genera Athesis and Patricia (22.9 million years ago (mya) [20.0–24.0] under the combined calibration scheme, 25.4 mya [21.8–27.5] with Nymphalidae calibrations only). For the genus Pteronymia, ages under the two calibration schemes were extremely similar (regression ages Nymphalidae calibration only (N) versus combined Nymphalidae and Solanaceae (S): N = 0.989*S, r2 = 0.997). Trees generated under the combined calibration scheme were used in all subsequent analyses.

Figure 1

BEAST dated species-level phylogeny of the genus Pteronymia and outgroups, based on molecular and morphological characters.

Main clades and secondary calibration points based on butterfly (red circles) and host-plant ages (green circles) are indicated and corresponding age priors are shown in the table inserted in the figure. The figure was generated with FigTree (http://tree.bio.ed.ac.uk/software/figtree/) and edited with Adobe Illustrator 4 (http://www.adobe.com/uk/products/illustrator.html). Butterfly pictures were taken by Keith Willmott and edited in Adobe Photoshop CS4 (www.adobe.com/products/photoshop/) and correspond to butterfly names in red in adjacent phylogeny.

After splitting from its sister lineage (the clade composed of the genera Episcada, Ceratinia, and Haenschia) 14.4 mya [12.3–16.3], the genus Pteronymia started diversifying about 10.6 mya [9.0–12.2], when it split into two major, but poorly supported, clades: the P. sao clade, 17 species, and the P. oneida clade, 36 species. The two main clades appearing in the morphological phylogeny were also recovered in the combined analysis, with the exception that P. latilla and P. tucuna grouped with two species in the P. zerlina clade. The male genitalia of both of the former species are rather different to remaining members of the P. zerlina clade, lacking the distinctive male genitalic synapomorphy of that clade in addition to also lacking a gnathos, the loss of which occurs within the genus only in P. latilla, P. tucuna, P. sao and P. obscuratus. Otherwise, none of the relationships implied in the combined phylogeny seem to contradict any strong morphological evidence.

Because of topological uncertainty of the Pteronymia phylogeny, we performed subsequent analyses both on the MCC tree and on a random subset of 100 trees from the posterior distribution.

Spatial patterns of diversification

We used georeferenced records (Supplementary Fig. S5) to analyse the spatial patterns of diversification of Pteronymia across nine biogeographic areas (Fig. 2). We performed biogeographical analyses using the software RASP 2.146. The analyses on the MCC tree and on the 100 trees yielded very similar results (Fig. 3, see Supplementary Fig. S6), and only the analyses on the MCC tree are presented here. Our biogeographic reconstruction suggests that the most likely ancestral area for the genus Pteronymia is the Western/Central Northern Andes (hereafter, Northern Andes), i.e., the area comprising the slopes of the Western and Central cordillera of Colombian (and Ecuadorian) Andes (Fig. 3), although there is uncertainty as to whether the origin of the genus was limited to this area, or also spanned neighbouring regions (see Supplementary Fig. S7). The two main clades (P. oneida and P. sao clades) also originated and started diversifying in the same area, with some uncertainty as to whether ancestral lineages spanned larger regions for the P. oneida clade (see Supplementary Fig. S6). A large proportion (55%) of speciation events occurred within the Northern Andes. In particular, the most likely ancestral area for the young and diverse P. zerlina clade was the eastern slopes of the Northern Andes. Rapid diversification subsequently occurred within the last 3.6 my [3.0–4.3] in this clade, which is coincident with the final uplift of the Eastern Cordillera of Colombia and the Venezuelan Cordilleras c. 5-2 mya9. The Central Andes appear relatively species poor compared to the Northern Andes; only 13 species occur in the Central Andes. These are the result of multiple independent colonization events (10 events), and to a much lesser extent local diversification (e. g., the clade encompassing P. hara and its sister clade, three species in this region). The oldest colonizations of the Central Andes were recovered in the P. sao clade (in the last 5.8 my [4.4–7.0]), where such colonizations happened at least five times.

Figure 2

Biogeographical regions used in the DEC and SDEC models for reconstruction of ancestral areas and the dispersal probability matrices for the different time slices.

Species richness in each areas are shown in the circles. The map was generated using ArcGIS 9.3 (http://www.edit.com/software/arcgis/), and edited with Adobe Illustrator 4 (http://www.adobe.com/uk/products/illustrator.html).

Figure 3

RASP historical biogeography inference (best maximum likelihood estimates on the MCC tree).

Major paleoenvironmental events are indicated by large coloured rectangles (light pink: drainage of the Pebas system; light yellow: hypothesized closure of the Isthmus of Panama). Colours of the little squares at the node and tips of the tree correspond to colours of the biogeographical areas, as indicated in the map inserted (taken from Fig. 2). The figure was generated with R (https://cran.r-project.org/) and edited with Adobe Illustrator 4 (http://www.adobe.com/uk/products/illustrator.html).

A small number of lineages colonized the Upper and Lower Amazon regions. Only one dispersal event was followed by local diversification in Amazonia, in the vestilla clade (the clade encompassing P. dispar and its sister clade). Three other colonizations of Amazonia occurred independently throughout the phylogeny, resulting in an overall very low Amazonian diversity, particularly in the lower Amazon. All these events happened within the last 5.6 mya [4.8–6.7]. The two species that colonized the Atlantic Forest arose from Amazonian lineages, and did not result in local diversification. Conversely, a high number of independent colonization events (15, according to the maximum likelihood estimates) occurred from the Andes to Central America. Most of those colonizations occurred within the last 5.0 my [4.0–5.9]. Colonization time of Central America was highly uncertain for P. fumida and P. spnov 4, because these species split from their sister lineages 10.6 mya [9.0–12.2] and 7.0 mya [5.6–8.3], respectively, and may have therefore colonized Central America any time during those periods. Colonization of Central America was sometimes, but rarely given the high number of colonizations, followed by local speciation (e.g. P. lonera and P. teresita, P. alcmena and P. gertschi).

Evolution of elevational range

We investigated the evolution of the elevational range and mean elevation of Pteronymia species on the MCC tree and on 100 trees from the posterior distribution. The phylogenetic signal and the tempo of evolution of the elevational range and mean elevation were assessed by estimating simultaneously the maximum likelihood values of the λ and δ branch scaling parameters47, respectively. A λ value of one indicates that the phylogeny correctly represents the trait covariance among species (Brownian motion model of evolution), while λ < 1 indicates that the phylogeny overestimates the trait covariance among species. A δ value of one means that the trait evolves at a constant pace along branches of the tree; δ < 1 indicates early changes in the character values followed by a slowing down of the evolution rate; while δ > 1 indicates accelerated evolution rate and species-specific adaptation. For the mean elevation and the elevational range, estimates of λ across all trees from the posterior distribution and for the MCC tree were not significantly different from one, meaning that trait evolution does not differ from a Brownian motion model. Estimates of δ for the mean elevation were also not significantly different from one (Table 3). For the lower and upper boundaries of the elevational range, the estimates of δ were significantly higher than one in >50% of the trees of the posterior distribution (Table 3).The MCC tree showed significantly higher estimates of δ for the lower boundary of the elevational range, but only marginally significant for the upper boundary (Table 3). Values of δ higher than one indicate an acceleration of the rate of evolution of elevation range and species elevational specialisation. Results differ between the mean elevation and elevational range perhaps because the mean is less variable than the range (e.g., species with different elevational ranges may have similar elevation mean). Reconstructions of ancestral mean elevation and range boundaries accounting for the inferred δ values are depicted on Fig. 4.

Table 3

Maximum likelihood estimates of δ for the mean and boundaries of the 95% elevational range, for the 100-tree posterior distribution (average values ± standard deviation), and for the MCC tree.

Trees		Mean Elevation	95% lower boundary	95% upper boundary
	δ	1.85 ± 0.60	2.96 ± 0.17	2.89 ± 0.23
Posterior	LRT	1.36 ± 1.63	5.72 ± 1.89	5.26 ± 2.62
distribution	p-value	0.377 ± 0.250	0.028 ± 0.044	0.045 ± 0.048
	% trees p-value < 0.05	4	87	70
	δ	1.55	2.99	2.86
MCC	LRT	0.625	5.16	3.10
	P-value	0.429	0.023	0.078

Likelihood Ratio Tests values when compared with the null model (δ = 1), corresponding p-values and fraction of trees for which δ is significantly different from 1 are reported.

Figure 4

Ancestral reconstruction of the mean and boundaries of the 95% elevational range (left: lower boundary, middle: upper boundary, right: mean).

For the lower and upper boundaries of the elevational range, trees were rescaled according to the δ value inferred (Table 2). The figure was generated with R (https://cran.r-project.org/).

Temporal patterns of diversification

We first investigated heterogeneity among clades for speciation and extinction rates using MEDUSA48 on the MCC tree and on 100 trees from the posterior distribution. No significant shift of diversification rates was found on the MCC tree. Five trees of the posterior distribution (5%) had at least one significant shift of diversification rates. Each of those shifts was found in less than 5% of the trees, and therefore considered as non-significant. We then investigated whether diversification rates varied through time by fitting time-dependent models of speciation and extinction rates49. The best fitting model was an exponential time-dependent speciation rate without extinction (Table 4). According to this model speciation rate decreased from 0.646 event per lineage per million year for the MCC tree and 0.538 ± 0.024 for the 100 trees at the origin of the genus (crown) to 0.148 for the MCC tree and 0.159 ± 0.002 for the 100 trees at present (Fig. 5). All other models, including the constant speciation rate model, were rejected at the threshold of ΔAIC > 2, strengthening the support for a decreasing speciation rate through time.

Table 4

Models of time-dependent diversification fitted on 100 trees from the posterior distribution, ranked by increasing AIC score.

Model	Par	logL	AIC	ΔAIC	λ	α	μ	ß
BVAR	2	−122.93 (0.286)	249.86 (0.571)	0.00	0.161 (0.001)	0.110 (0.003)
BVAR - DCST	3	−122.93 (0.286)	251.86 (0.571)	2.00	0.161 (0.001)	0.110 (0.003)	5.00E-08 (5.00E-08)
BCST	1	−125.39 (0.258)	252.79 (0.516)	2.93	0.233 (0.001)
BVAR - DVAR	4	−122.93 (0.286)	253.86 (0.571)	4.00	0.161 (0.001)	0.110 (0.003)	5.00E-08 (5.00E-08)	−0.004 (0.003)
BCST - DCST	2	−125.39 (0.258)	254.79 (0.516)	4.93	0.233 (0.001)		<1.00E-08 (<1.00E-08)
BCST - DVAR	3	−125.39 (0.258)	256.79 (0.516)	6.93	0.233 (0.001)		<1.00E-08 (<1.00E-08)	0.005 (2.00E-04)

Mean values of parameters are indicated followed by the standard deviation in brackets. BCST = constant speciation, BVAR = time-dependent speciation, DCST = constant extinction, DVAR time-dependent extinction. logL = likelihood of the model, AIC = AIC score, ΔAIC = difference of AIC between the each model and the best fitting model, λ = speciation rate at present, α = coefficient of time variation of the speciation rate, μ = extinction rate at present, ß = coefficient of time variation of the extinction rate.

Figure 5

Speciation rate through time estimated by the best fitting model of diversification (Table 4).

The model was fitted on the MCC tree and on 100 trees. The dotted line corresponds to the speciation rate of the MCC tree. The plain line corresponds to the mean speciation rate from the 100 trees, and dashed lines correspond to the 95% confidence interval. The figure was generated with R (https://cran.r-project.org/).

Discussion

Our extensive molecular sampling of the genus Pteronymia encompassing multiple subspecies and combined with morphological and distributional data enabled us to redefine species boundaries in the genus, resulting in six additional recognized species in the genus. The resulting taxonomic changes now make Pteronymia the most species-rich ithomiine genus, with 53 species. Morphological characters that are typically useful in diagnosing ithomiine species, such as in the genitalia and wing androconia, were clearly unhelpful in these cases, although mimetic wing pattern and life histories, where known, often did show significant variation. As knowledge of the biology of different populations and more material for molecular study become available, additional changes to the species taxonomy may be needed in future. At deeper levels of the phylogeny, despite Pteronymia showing some of the greatest diversity within the Ithomiini in terms of the male and female genitalia, wing venation and androconia, immature stage morphology and biology, this variation proved remarkably unhelpful in resolving the phylogeny, perhaps due to high rates of morphological character evolution. The combination of morphological and molecular characters obviously increased the resolution of the phylogeny compared to morphology alone, but many clades were still surprisingly poorly supported, and in any case the phylogeny of Pteronymia was less resolved than those of other ithomiine genera6285051, even when considering only molecular characters. This may be due to incomplete lineage sorting, rapid diversification or hybridization. While hybridization and introgression seem common in other mimetic butterflies, such as Heliconius525354, little is known about such processes in ithomiine butterflies. Future genomic data may shed light on demographic processes and gene flow between species, which may have contributed to the poor support seen in some nodes in the Pteronymia phylogeny.

Our calibration strategy, based on a combination of secondary calibrations derived from host-plant and Nymphalidae phylogenies, required the recalibration of a published Solanaceae phylogeny43 (host-plants of most Ithomiini) because original ages were biased toward present, which precludes using those data as maximum calibrations. Our new calibration scheme, based on a secondary calibration extracted from a fossil-dated phylogeny of Angiosperms44, inferred ages of Solanaceae lineages that were about 25% older than previous estimates43, but mostly within the 95% confidence ranges of those estimates. The host-plant ages used in this study as maximum constraints are based on one of the ‘youngest’ hypotheses for Angiosperm diversification and recent studies suggest older ages of Solanaceae. A phylogenetic inference using genomic data found a much older origin of Angiosperms55, but this had a moderate impact on the stem age of the order containing Solanaceae (Solanales: ca. 92 mya55 versus 85.9 mya44, with overlapping 95% confidence ranges55) and presumably on that of Solanaceae (not inferred in that study). The recent description of a 52.2 my old Physalis fossil56, a solanaceous genus that is inferred to be 9.1 [5.9–12.9] my old in our study (stem age), is much more challenging for the ages of Solanaceae and Angiosperms as a whole. It is possible that this recently described Physalis fossil represents an earlier lineage of Solanaceae. The inflated calyx is found in many genera throughout the ‘berry’ clade of Solanaceae (i. e., the subfamily Solanoideae, where the stem is the MRCA of Nicotiana and Solanum and the crown the MRCA of Latua and Solanum57) and transcription factors governing this character appear to be plesiomorphic in the family585960. Placement of the fossil at an earlier diverging node, such as that of the berry clade, within Solanaceae is less contradictory in terms of Angiosperm evolution as a whole, and would not have a major effect in pushing back the stem and crown node age of Solanaceae. In terms of our findings here, older Solanaceae age estimates do not affect our time-calibrated Pteronymia tree given that butterfly clades appear to be much younger than their corresponding host-plant groups even when host-plant ages were inferred from one of the ‘youngest’ Angiosperm time-calibrated phylogenies44.

Our combined calibration strategy therefore resulted in ages that were consistent with those inferred in the Nymphalidae phylogeny39, but older than those inferred in the higher-level Ithomiini phylogeny38. This difference stems from several factors. The new ages of the Solanaceae lineages (this paper) used as maximum calibration points were 25% older than the previous estimates used in the higher-level Ithomiini phylogeny38, and we used the oldest boundary of the 95% credibility interval as the (hard) maximum age of corresponding ithomiine lineages, instead of the mean as in the higher-level Ithomiini phylogeny38. The rationale for this is that minimum (such as fossil-based) and maximum (such as host-plant-based) calibrations ought to be conservative and therefore account for uncertainty in calibration age61. In our case, where we implemented conservative uniform priors between maximum ages and present, node ages in the Pteronymia phylogeny were not necessarily attracted towards maximum ages, they were just allowed to go as far as those ages. Since host-plant ages provide maximum calibrations, they need to be combined with minimum calibrations, in a way similar to fossil-based minimum calibrations that need to be combined with at least one maximum calibration point61. Given the absence of fossils for Danainae, here we supplemented the host-plant-derived calibrations with secondary calibrations extracted from the Nymphalidae phylogeny39, which was calibrated with a combination of fossils and host-plant constraints. These calibrations provided both minimum and maximum ages, and therefore imposed a stronger prior on ithomiine ages than the host-plant-derived calibrations, but to be as conservative as possible we used a uniform prior spanning the 95% credibility interval of the Nymphalidae ages that were used for calibration.

Our results based on the biogeographic reconstruction of the genus Pteronymia and evolution of elevational range were consistent across trees of the posterior distribution despite the relatively poor resolution of the MCC phylogeny, and revealed the fundamental roles played by the Northern Andes in the diversification of the group, in multiple ways. The Northern Andes have probably been (1) the area of origin of the group (although the inference for the root is not well resolved, this region appears in all the potential ancestral areas); (2) the centre of early and sustained local diversification; and (3) a source of recent colonizations to lowland areas and to Central America, that led to an accelerated evolution of the elevational range.

The biogeographical pattern of diversification, where the Andes play a central role, resembles those found in other Andean ithomiine genera, such as Hypomenitis35, and to a lesser extent Napeogenes6 and a clade of Oleria4. These groups originated and diversified in the Northern Andes during the last 10 million years, a period during which the Northern Andes experienced different phases of intense uplift9. This period of orogeny involved great landscape transformations that may have affected the dynamics of Andean lineages, by isolating populations within deep valleys or on both sides of the Andes, but also in modifying the climatic conditions. Also, the slopes of the Andes offer a great number of opportunities for ecological speciation due to significant environmental turnover resulting in high habitat, host-plant and predator diversity. In the genus Pteronymia, the observed decrease of speciation rate with time and the low support for basal nodes may indicate a rapid diversification driven by adaptive factors. Elevations up to 2000m probably already existed by 10 million years ago in the Northern Andes62. Given the timing of diversification of the genera Hypomenitis, Napeogenes and Pteronymia in the last 15 to 5 million years, the majority of speciation events have likely not coincided with the appearance of newly available elevations. Instead, speciation was probably facilitated by an already well-established and diverse ecosystem. The slopes of the Andes harbour not only a diversity of habitats and host-plant communities6364, but also a diversity of mimicry rings28. Pteronymia is one of the most diverse ithomiine genera in terms of mimetic wing colour patterns. Shifts in colour pattern are known to drive speciation in other mimetic butterflies such as Heliconius6566, where colour pattern is considered as a ‘magic’ trait, i.e., a trait that is both under disruptive selection and associated with assortative mating656768. Mimetic butterflies often harbour multiple geographic races with different colour patterns, and it has been shown in Heliconius butterflies that interracial hybrids that display an intermediate, non-mimetic colour pattern suffer higher predation66. Shifts in colour pattern may therefore cause postzygotic reproductive isolation. In addition, Heliconius butterflies tend to prefer mates with their own colour pattern over conspecifics with a different colour pattern656768, thereby driving prezygotic reproductive isolation. In this case, loci involved in mate preference and in colour pattern are tightly linked69. In Ithomiini, experimental evidence for the role of colour pattern as a mating cue is absent due to the difficulty of maintaining and rearing Ithomiini in captivity, but observation suggests that this may be the case70. Moreover, shifts in colour patterns have been shown to be statistically associated with cladogenesis in phylogenies3671, consistent with a role of colour pattern in reproductive isolation. Shifts in colour pattern may therefore have contributed to Andean diversification in the genus Pteronymia.

In contrast with other ithomiine clades4635, little local diversification in Pteronymia occurred outside of the Northern Andes. An exception to this is the P. vestilla clade, which expanded into and partly diversified in lowland areas. According to our reconstructions, this lineage dispersed into the Upper Amazon from the Northern Andes and started diversifying locally around 5.6 mya [4.8–6.7]. Local diversification was then followed by expansions or dispersal into other areas, including the lower Amazon and Guiana shield, the Atlantic forest, and Western lowlands of Colombia and Ecuador. Only two extant lineages of this clade occur in Amazonia. The Upper Amazon region has experienced important environmental changes since the Miocene. During most of the Miocene, the region was covered by the large lake and shallow swamps of the Pebas system, which was connected northward to the Caribbean Sea and potentially westward to the Pacific Ocean9. This particular ecosystem may have had a major influence on the diversification of clades restricted to forest habitats1423 by limiting occurrence and therefore speciation to the margins of this region, but probably also by preventing dispersal between Amazonia and the Andes and between the Northern and the Central Andes (via the West Andean Portal, a low elevation gap between the Northern and Central Andes72). The demise of the Pebas system started around 10-8 mya and it was rapidly drained eastward, leading to the establishment of the modern Amazon basin, probably around 7 mya9. Our results suggest that dispersal and diversification in Amazonia did not occur simultaneously with the drainage of Pebas, but later on. The diversification of the P. vestilla clade in the Upper Amazon occurred fairly recently (5.6 mya [4.8–6.7]), as did the independent colonizations of the Upper Amazon by the ancestor of P. veia_WEST and P. tucuna (1.5 mya [0.9–2.2]), and the ancestor of P. sao and P. obscuratus (2.7 mya [1.8–3.6]). The timing of the colonization of the Upper Amazon by P. forsteri, which split from its sister clade 8.3 mya [6.8–10], is much less precise, since it could have happened any time during this period.

Because little local diversification in Pteronymia occurred outside of the Northern Andes, most of the diversity in non-Andean regions is due to colonization out of the Andes. Species diversity in the Central Andes and Central America built up through the accumulation of independent dispersal events, rarely followed by speciation events. Many of these events occurred without strong elevational shifts suggesting that dispersal may have been facilitated by the existence of similar ecological conditions in montane areas. By contrast, dispersal toward lowland areas, such as the Upper Amazon, may have entailed more adaptations to fit different bioclimatic conditions or host-plants, which likely explains the rare occurrence of such events.

There is a debate surrounding the timing of the closure of the Panama Isthmus. The hypothesis that it occurred very recently (5-3 mya) has been widely adopted in the literature (see12). However, both geological and paleontological findings suggest a possible much earlier appearance of land masses, possibly as early as the early or middle Miocene (e.g. refs 12, 73, 74, 75). In our biogeographical analysis, although interchanges between Central America and South America were allowed earlier, most of the colonization events of Pteronymia lineages toward Central America have occurred during the last 5.0 [4.0–5.9] million years, in agreement with the first hypothesis. However, the colonization time of Central America by P. fumida and P. sp. nov. 4, two species with long branches, is uncertain, and may have happened much earlier.

The pattern of diversification inferred for Pteronymia is very similar to biogeographic patterns of other taxa described in the literature. For example, in vertebrates, the Northern Andes were a major centre of diversification for glassfrogs (Allocentroleniae), which subsequently fed the adjacent areas through dispersal, including the Central Andes and the non-Andean regions3. A similar conclusion was reached in the Thraupini tanagers, which also diversified during the last 10 million years, with higher rates of colonization out of the Northern Andes (mostly toward the Central Andes) rather than into that region76. The bat genus Sturnira also diversified during the last 10 million years, from the Northern Andes toward the rest of the Neotropical region77. Similarly, in plants, many studies report the Northern Andes as a centre of diversification and a source for adjacent areas24. The Rubiaceae subfamily Cinchonoideae originated and mostly diversified in the Northern Andes, with subsequent colonization of both lowland and highland adjacent areas72. Extremely high rates of diversification were reported in the Andes for the genus Lupinus7. The Campanulaceae experienced higher rates of diversification in the Andes that were correlated with paleoelevations of the Andes8, suggesting that the Andes have directly affected diversification in this family. Taken together, these findings are consistent with the idea that, at least during the last 10 million years, the Northern Andes have been an important source of biodiversity in the Neotropics probably due to geological, climatic and biotic factors, with local diversification followed by lineage dispersal throughout other Neotropical regions.

Methods

Morphological characters and phylogeny

Prior to our study, 47 species were listed in the genus Pteronymia, but this figure increased to 53 after taxonomic revision based on new data (Table 1). Forty-six adult and 41 immature (i.e., larval and pupal) morphological characters were examined in 52 Pteronymia species (after our revision, see Supplementary Methods S1, Supplementary Table S2, Supplementary Figs S8–S9; P. dispar rev. stat. was not coded) and two outgroup species (Episcada apuleia and Dircenna adina).

A Maximum Parsimony topology was estimated in TNT v. 1.5-beta78 using the New Technology Search, with all four search methods–ratchet, tree-fusing, tree-drifting and sectorial, with default parameters, 100 random additional sequences and random seed equal 0, with all characters equally weighted. The majority-rule consensus tree, consistency index (CI) and the retention index (RI) were computed in Winclada79. The stability of each branch was estimated using the non-parametric bootstrap test, with 1,000 replicates and 100 random taxon additions, and Bremer support, using the script Bremer.run in TNT.

Molecular characters and phylogeny

We used a total of 166 Pteronymia specimens for molecular analyses, representing 41 of the species recognized prior to our revision, and 47 of the species recognized after our revision (Table 1; see Supplementary Table S1). Species with no molecular data were P. alcmena, P. alicia, P. calgiria, P. fumida, P. glauca and P. peteri.

We used de novo (ca. 85%) and published63550518081 (ca. 15%) sequences from five gene regions to infer a molecular phylogeny (see Supplementary Table S1): the mitochondrial region spanning the mitochrondrial genes cytochrome oxidase c subunit 1, leucine transfer RNA and cytochrome oxidase c subunit 2 (CO1, tRNAleu, CO2, 2356 bp), and the nuclear genes tektin (735 bp) and Elongation Factor 1 alpha (EF1A, 1259 bp). Species coverage was 98% for the mitochondrial fragment, 73% for tektin and 63% for EF1A. Primers and PCR conditions followed previously described conditions6. In addition, 52 ithomiine and danaine outgroup species were selected35 (see Supplementary Table S1). The dataset was then partitioned by gene and codon positions and the best models of substitution for optimized sets of nucleotides were selected over all models implemented in (1) RAxML82 and (2) MrBayes83, using the ‘greedy’ algorithm and linked rates implemented in PartitionFinder 1.1.184 (see Supplementary Table S3).

We performed a maximum likelihood phylogenetic inference using RAxML82 on the Cipres server85. In addition, we performed a Bayesian inference of the phylogeny using MrBayes 3.2.283 on the Cipres server85. Substitution models of each partition were re-estimated in MrBayes 3.2.2 using the reversible-jump MCMC86. Two independent analyses were run for 10 million generations, with four Monte Carlo Markov chains each and a sampling frequency of one out of 10,000 generations (resulting in 1,000 posterior trees). After checking for convergence, the posterior distributions of the two runs were combined, with a burnin of 10%. The maximum clade credibility tree with median node ages was computed using TreeAnnotator 1.6.245. The resulting tree was used to investigate topology and species boundaries.

Molecular dating

In order to estimate a time-calibrated phylogeny of Pteronymia, we combined two types of time constraints: the age of larval host-plants and age estimates from higher-level phylogenies of butterflies. Host-plant ages can be used as maximum age constraints in phylogenies of mono- or oligophagous herbivores3839, assuming that such herbivorous taxa diversified only after the emergence of their host-plant lineages. Most Ithomiini feed on Solanaceae, which represents a host-plant shift from the ancestral host-plants of Danainae (Apocynaceae87). In a previous study aiming at dating a higher-level phylogeny of the tribe Ithomiini38, ages of several Solanaceae lineages inferred from a dated phylogeny of Solanaceae43 were used as maximum calibrations for Ithomiini clades feeding on specific Solanaceae lineages. The ages of clades in the Solanaceae phylogeny43 were minimum age estimates because Solanaceae fossils known at that time were placed conservatively at the stem ages of lineages with which they shared morphological synapomorphies, and because fossils in general can only provide minimum age estimates for clades. Therefore, using those minimum Solanaceae ages as maximum calibrations for Ithomiini lineages may strongly underestimate the ages of the butterfly lineages, especially when using mean or median age estimates instead of older bounds. Here, we took advantage of a recent calibration of the family-level phylogeny of the Angiosperms based on 151 fossils44 to recalibrate the Solanaceae phylogeny using the inferred stem age of Solanaceae, i.e., 66.6 mya, as a calibration point in a Bayesian framework (Supplementary Methods S1.2). Due to the recent description of a 52.2 my old Solanaceae fossil placed in the extant genus Physalis56 that could have dramatic consequences on the age estimates of Solanaceae and Angiosperms as a whole, we also ran an analysis without calibration based on host-plant ages. Indeed, much older ages of Solanaceae lineages as those implied by this discovery would have hardly any impact on the ages of Ithomiini, and removing host-plant derived calibrations is therefore a conservative way of testing the influence of older host-plant ages. Since the two calibration strategies yielded almost the same ages (see results) we performed the biogeographic and diversification analyses on the tree calibrated using the combined calibration strategy detailed below. For further discussion on the recently described Solanaceae fossil and its potential effect on host plant age estimates, see Discussion above.

We extracted the new maximum ages of several Solanaceae lineages from our newly calibrated phylogeny in order to provide maximum calibrations for Ithomiini lineages that feed on them87. The choice of the host-plant calibrations (Fig. 1) followed that of the genus-level Ithomiini phylogeny study38, with the following exceptions. We excluded the Solanum calibration, because the two lineages feeding exclusively on Solanum (the Mechanitina and the clade comprising the Oleriina, Dircennina, Godyridina, Ithomiina and Napeogenina) do not have a sister relationship in our and in other Ithomiini phylogenies2835. We also excluded the Cestrum calibration (ithomiine subtribe Godyridina, excluding the genera Veladyris and Velamysta) because this highly diverse genus had a low sampling fraction of 20% in the Solanaceae phylogeny and was not resolved as monophyletic43. Finally, for simplicity we excluded the calibrations based on Solanaceae genera Brunfelsia and Lycianthes, because they apply to young butterfly lineages (the subtribe Methonina and the clade comprising the genera Oleria and Ollantaya, respectively3839) and would have no effect on the age estimates. A uniform prior was used for host-plant-derived calibrations. The upper boundary of the prior was set to the upper boundary of the 95% credibility interval of the age of the host-plant lineage and the lower boundary of the prior was set to 0 (present) (Fig. 1).

As a second source of calibration points we used age estimates from Wahlberg et al.’s39 dated Nymphalidae phylogeny. We defined seven secondary calibrations points which were set to a uniform prior bounded by the upper and lower boundaries of the 95% HPD of the ages inferred in the Nymphalidae phylogeny39 (Fig. 1).

Sequences of all Pteronymia specimens were combined into a consensus sequence for each species to maximize sequence coverage for each species35. We used PartitionFinder 1.1.184 to select the best partition scheme applying to this new dataset, where only the models implemented in BEAST were tested (see Supplementary Table S4). A time-tree was generated in BEAST 1.7.545 under an uncorrelated lognormal relaxed clock using the same outgroups as previously, and the dating procedure described above. To select the tree prior (Yule versus Birth-Death), we ran analyses with each type of prior and used Stepping Stone sampling88 to estimate marginal likelihood (MLE). This method is not implemented in BEAST 1.7.5, thus we performed this analysis on BEAST 1.8.2 on the molecular dataset only (BEAST 1.8.2 cannot handle morphological characters). The MLE were then used to compute Bayes Factors (BF), which supported the Yule model (BF = 4.44). We also ran analyses on the total evidence on BEAST 1.7.5 under the Yule and the Birth-Death model to compare age estimates. Both types of analyses yielded virtually identical ages (regression of ages under Birth-Death (BD) on ages under Yule (Y): BD = 0.9985*Y, r2 = 0.999). Analyses with and without morphological data also produced very similar ages (regression of ages with molecular data only (M) on ages with the combined evidence (C): M = 0.9892*C, r2 = 0.997). Analyses were run for 100 million generations each on the Cipres server85 and on a desktop computer depending on the version of BEAST, and trees and parameters were sampled every 100,000 generation. After checking for parameter ESS, a 10% burnin was applied to the posterior distribution and the maximum clade credibility tree with median node ages (hereafter, MCC tree) was computed using TreeAnnotator 1.6.245. Given the results on the tree prior (Yule versus Birth-Death, see above), subsequent analyses were performed on the trees generated under a Yule prior. Because many nodes of the phylogeny had moderate or poor support, we conducted most of the analyses outlined below on the MCC tree and on a random subset of 100 trees extracted from the post burnin posterior distribution of trees obtained from the BEAST run.

Geographic distribution and elevational range for all extant species were obtained from our own records, museum collections and collaborators (Table 1, Supplementary Fig. S5). The biogeographic history of Pteronymia was estimated using the Maximum Likelihood Dispersal-Extinction-Cladogenesis (DEC) model89 on the MCC tree, and its extension for multiple trees (Statistical DEC, or S-DEC) on the 100 trees extracted from the posterior distribution, using the software RASP 2.146. We defined nine biogeographical areas (Fig. 2): Central America (A), Western lowlands (B), Slopes of the Western/Central Cordillera of Ecuador and Colombia (C), Central Andes (D), Slopes of the Eastern Cordillera of Colombia and Venezuela (E), Upper Amazon (F), Lower Amazon (G), Atlantic Forest (H), and Guiana Shield (I), and the maximum number of ancestral areas in our analyses was set to four reflecting the maximum number of areas occupied by extant species. We implemented three time slices (11-8, 8-5, 5-0 million years ago, Fig. 1) with different dispersal probabilities, which took into consideration major geological events through time9.

Ancestral state reconstruction and evolution of elevational range

Elevational range (elevation range, i.e., boundaries of the elevational interval containing 95% of the records and mean elevation) for all extant species were extracted from the distributional database computed above (Table 1). Evolution of elevational range was investigated on both the MCC tree and the subset of 100 trees as follows. The phylogenetic signal and the tempo of evolution of the mean elevation and the boundaries of the elevational range were assessed by estimating simultaneously the values of the λ and δ scaling parameters47 that maximized the likelihood of the data using BayesTraits v290. A λ value of one indicates that the phylogeny correctly represents the trait covariance among species, while a value of 0 indicates that the trait evolution is independent of the phylogeny. A value of λ smaller than one indicates that the phylogeny overestimates the trait covariance among species. A δ value of one means that the trait evolves at a constant pace along the branches of the tree; δ < 1 indicates early changes in the character values followed by a slowing down of the evolution rate, such as that entailed by an adaptive radiation; while δ > 1 indicates accelerated evolution rate and species-specific adaptation. The MCC tree was then rescaled with the corresponding δ and λ values inferred for the mean elevation and the upper and lower boundaries of elevational range, when those differed from 1, such that the evolutionary rate of the elevation attributes on the transformed tree was constant. The resulting trees were used to infer ancestral values of the attributes (assuming a constant evolution rate), using the function contMap of the R package phytools91.

To investigate the pattern of speciation and extinction rate variations through time and across lineages, we chose not to use BAMM 2.5.092 because of recent criticisms on uninformative priors and biased estimates of diversification rates93. Instead, we implemented a two-step procedure. We first used MEDUSA48 implemented in the R package geiger94 to detect shifts on the selection of 100 random trees of the posterior distribution, and on the MCC tree. MEDUSA has also been criticised because of its inflated false discovery rate and biased estimates of diversification rates95. To overcome these shortcomings and to implement time-dependent models of diversification, we used the method developed by Morlon et al.49, a maximum likelihood approach that accommodates time dependent birth-death processes and enables to test for rate shifts. We used as a rule that if a shift was present in at least 5% of the trees of the posterior distribution in the MEDUSA analysis, this shift was tested with Morlon et al.(2011)’s method49. Shifts detected in fewer than 5% of the posterior distribution were considered non-significant and not tested with Morlon et al. (2011)’s method49. We fitted 6 models of diversification49: 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. Time-dependency was modelled using an exponential function. Models were compared using AIC scores. The root of the tree was always excluded from the analyses. Our phylogeny included all known species so we set the sampling fraction to 1. All models were fitted on the MCC tree and on the 100 trees sampled from the posterior distribution.

Additional Information

How to cite this article: De-Silva, D. L. et al. North Andean origin and diversification of the largest ithomiine butterfly genus. Sci. Rep. 7, 45966; doi: 10.1038/srep45966 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.