Mesoamerica is a cradle and the Atlantic Forest is a museum of Neotropical butterfly diversity: insights from the evolution and biogeography of Brassolini (Lepidoptera: Nymphalidae)

Regional species diversity is explained ultimately by speciation, extinction and dispersal. Here, we estimate dispersal and speciation rates of Neotropical butterflies to propose an explanation for the distribution and diversity of extant species. We focused on the tribe Brassolini (owl butterflies and allies), a Neotropical group that comprises 17 genera and 108 species, most of them endemic to rainforest biomes. We inferred a robust species tree using the multispecies coalescent framework and a dataset including molecular and morphological characters. This formed the basis for three changes in Brassolini classification: (1) Naropina syn. nov. is subsumed within Brassolina; (2) Aponarope syn. nov. is subsumed within Narope ; and (3) Selenophanes orgetorix comb. nov. is reassigned from Catoblepia to Selenophanes . By applying biogeographical stochastic mapping, we found contrasting species diversification and dispersal dynamics across rainforest biomes, which might be explained, in part, by the geological and environmental history of each bioregion. Our results revealed a mosaic of biome-specific evolutionary histories within the Neotropics, where butterfly species have diversified rapidly (cradles: Mesoamerica), have accumulated gradually (museums: Atlantic Forest) or have diversified and accumulated alternately (Amazonia). Our study contributes evidence from a major butterfly lineage that the Neotropics are a museum and a cradle of species diversity. applyparastyle “fig//caption/p[1]” parastyle “FigCapt”


INTRODUCTION
Terrestrial biomes are heterogeneous in terms of species composition and distribution, and some areas have accumulated greater diversity and endemicity than others. A combination of three evolutionary mechanisms explains this biogeographical pattern: species origination, extinction and dispersal (Goldberg et al., 2005;Jablonski et al., 2006). For example, a species-rich region can have attained its present diversity by recent and exceptionally high rates of speciation, becoming a 'cradle of diversity', or by the persistence of old lineages via low species extinction rates, being a 'museum of diversity' (e.g. Stebbins, 1974;Gaston & Blackburn, 1996;Bermingham & Dick, 2001;Richardson et al., 2001). Furthermore, a species-rich region might be a 'source of diversity' when dispersal out of the region is high or a 'sink of diversity' when regional species diversity is accrued from external sources via dispersal into the region (Antonelli et al., 2018b). Characterizing the interplay among speciation, extinction and dispersal is thus highly relevant for understanding the origin and evolution of total regional diversity and endemism (Wiens & Donoghue, 2004;Roy & Goldberg, 2007).
Much of the biodiversity and endemism of the world are concentrated in Neotropical rainforests in Central America, the Amazon Basin and the Atlantic Forest in south-eastern Brazil. However, there is no consensus on whether these currently separate rainforests were connected once or multiple times, nor for how long (Jaramillo & Cárdenas, 2013). It is also unclear whether Neotropical rainforests represent centres of ecological stability that share a common evolutionary history (Moritz et al., 2000). Mammals and birds in Neotropical biodiversity hotspots, including the Brazilian Atlantic Forest and rainforests in Mesoamerica and the neighbouring Chocó on the north-western (NW) side of the Andes (Myers et al., 2000), apparently underwent rapid diversification (Igea & Tanentzap, 2019). In other groups, a larger number of lineages seem to have dispersed out of those areas compared with the rest of the Neotropics (e.g. Machado et al., 2018;Toussaint et al., 2019), suggesting that such rainforests are both cradles and sources of diversity (Igea & Tanentzap, 2019). Other Neotropical regions are, in contrast, museums of diversity, given the old and gradual accumulation of species through time, such as in Amazonia for Euptychiina (Nymphalidae) (Peña et al., 2010) and Troidini butterflies (Papilionidae) (Condamine et al., 2012). It thus seems that the checkered spatial assemblage of rainforest areas where speciation rates are high plus those where extinction rates are low might explain why the Neotropical region as a whole is considered both a cradle and a museum of diversity (McKenna & Farrell, 2006;Moreau & Bell, 2013;Antonelli et al., 2018a). However, the impacts of geological and palaeoenvironmental changes in the Neotropics together with the roles of speciation and dispersal have seldom been elucidated jointly among taxa occurring in individual rainforest areas.
Here, we examine biogeographical patterns for butterflies in the tribe Brassolini (Nymphalidae: Satyrinae), an exclusively Neotropical monophyletic group (Freitas & Brown, 2004;Wahlberg et al., 2009;Espeland et al., 2018;Chazot et al., 2019a) consisting of 17 genera and 108 species found mainly in the rainforests of Mesoamerica, the NW slope of the Andes (e.g. including Chocó), the Atlantic Forest and Amazonia (for an overview, see Penz, 2007). Studies that focused on individual brassoline genera have used adult morphology (Penz, 2008(Penz, , 2009aGarzón-Orduña & Penz, 2009) and DNA data (Penz et al., 2011a;Shirai et al., 2017) to infer species-level phylogenies, but several deep nodes had weak or even conflicting phylogenetic support in tribal-level analyses (Penz et al., 2013). We infer a time-calibrated species phylogeny of Brassolini by merging all previous datasets and generating new DNA and morphological data to infer the diversification and biogeographical histories of these butterflies. From a biogeographical perspective, if physical connection among Neotropical biomes existed during the Neogene and Quaternary periods (i.e. the past 23 Myr) (Hoorn et al., 2010), we expect that this will be reflected in the timing and magnitude of dispersal and speciation trends through time. To assess this expectation, we estimate dispersal rates and speciation for Brassolini butterflies.

Morphological and total-evidence datasets
Our morphological dataset included 255 discrete characters (201 from adults and 54 from early stages) for 64 species, 44 of them corresponding to the same species included in the molecular dataset. The morphological matrix included all 17 Brassolini genera and was concatenated from previous studies (Penz, 2007(Penz, , 2008(Penz, , 2009aGarzón-Orduña & Penz, 2009;Penz et al., 2013) by eliminating duplicate characters and scoring new data for species missing entries. Specimens used were the same as those listed in the original publications.
The combined molecular and morphological dataset comprised all genera and 84 of the 108 Brassolini species (missing data coded as '?'). We called this the total-evidence dataset, which we used to infer a timecalibrated phylogeny.

distribution dataset and bioregion deliMitation
We advocate a data-informed bioregionalization of Brassolini geographical distribution. Specifically, we expect that parametric bioregionalization will inform better on the best-fitting delimitation of areas for Brassolini compared with any arbitrary assignment of areas within the hierarchical classification by Morrone (2014a, b) (subregions, dominions or provinces).
We delineated bioregions (sensu Vilhena & Antonelli, 2015) using geo-referenced occurrence points of Brassolini identified with taxonomic certainty at the species level, of which 6378 points came from GBIF (https://gbif.org; retrieved on 21 January 2019) and 881 points from the ATLANTIC BUTTERFLIES dataset (Santos et al., 2018). Both databases follow the current taxonomic understanding of Brassolini; thus, we did not encounter dubious species names. We used the R v.3.5.3 (R Core Team, 2020) package coordinatecleaner v.2.0-11 (Zizka et al., 2019) for automated cleaning of species occurrences. In addition to four observations of Caligo eurilochus from tropical butterfly houses in Germany that were removed a priori, the software flagged 81 potential errors, such as occurrence records in the ocean, geographical coordinates identical to country centroids, and occurrences in close proximity to biodiversity institutions, such as natural history museums (Supporting Information, Fig. S1). We uploaded the 7174 vetted geo-referenced occurrences, all within the known natural geographical range of Brassolini, to the Web application infoMap bioregions (Edler et al., 2017) at http://bioregions.mapequation.org/. We set the maximum and minimum values for cell sizes to four and two, and for cell capacities to 40 and four. We used different values of cluster cost, from 0.1 to 3.0, in steps of 0.1. The number of inferred areas remained stable between values of 1.8 and 2.4 and, thus, received the strongest support.

phylogenetic inference and divergence tiMe calibration
We ran preliminary phylogenetic analyses using each gene and morphological dataset to check for tree topology conflicts (Figs S2-S4, Tables S2-S3). The six molecular datasets were partitioned by codon position, and the morphological dataset was partitioned by using homoplasy scores calculated through implied weighting parsimony (Rosa et al., 2019). We calculated homoplasy measurements (f) using TNT v.1.5 (Goloboff & Catalano, 2016) under the default concavity parameter k = 3, and we used these values to subdivide the morphological dataset into 11 partitions (Supporting Information, Table S3). All seven phylogenetic analyses (six gene loci and morphology) were run using Mrbayes v.3.2.6 (Ronquist et al., 2012) via the CIPRES Science Gateway v.3.3 (Miller et al., 2010). We performed model averaging over all substitution models within the GTR model family (Huelsenbeck et al., 2004), taking into account rate variation across sites (+I and +Γ models) for the molecular partitions, and we applied the Markov (MKv) substitution model (Lewis, 2001) for the morphological partitions. The Bayesian inferences were run two independent times for 50 million generations, sampling trees every 5000 th generation and discarding the first 25% of the sampled distribution as burn-in. After checking that the logarithmic probabilities reached a stationary distribution, the average standard deviations of split frequencies were < 0.005, PSRF values close to 1.000, and the estimated sample sizes (ESS) > 200, we summarized the post-burn-in sampled trees using the 50% majority-rule consensus method.
We inferred time-calibrated species trees using the starbeast2 v.0.15.5 package (Ogilvie et al., 2017) available in BEAST v.2.6.3 (Bouckaert et al., 2014. We removed the species Opsiphanes camena from the dataset because it had only 325 bp of the COI locus, and its phylogenetic position was not stable in the preliminary analyses. We used a gene-based partitioning strategy and unlinked tree models in the multispecies coalescent framework. We performed model averaging over 31 substitution models available in the bModeltest v.1.2.1 package (Bouckaert & Drummond, 2017). We applied uncorrelated relaxed-clock models for all partitions because strict molecular clock models were rejected based on the stepping-stone method (Supporting Information, Table S4). We estimated two molecular clock rates, one for the mitochondrial locus and one for the nuclear loci, by using the gamma distribution for the clock rate priors [shape parameter (α) = 0.01 for mitochondrial clock and 0.001 for nuclear clock; scale parameter (β) = 1.0]. We used the Yule tree model plus two molecular clocks because it received decisive support, based on path-sampling steps, over other models with different combinations of one, two or six clocks and over the Yule and birth-death tree models (Supporting Information, Table S5). We added the morphological dataset onto the species tree, and we conditioned the characters to the Markov substitution model (Lewis, 2001). BEAST 2 assigns morphological partitions with unlinked site models based on characters having the same number of states. All morphological characters were considered to evolve under the same species tree and relaxed clock model, with an uninformative rate prior (1/x distribution).
Inasmuch as there are no described fossils of Brassolini, we relied on secondary calibrations to time calibrate the species tree. Conservatively, we used uniform distributions encompassing the 95% highest posterior density (HPD) from a large fossil-calibrated, genus-level butterfly phylogeny (Chazot et al., 2019a) to constrain six outgroup nodes (Table 1; Supporting Information, Fig. S5).
The analyses were run four independent times for 300 million generations, via CIPRES. We sampled 20 000 trees from the posterior distribution and discarded the first 25% as burn-in. We checked that the ESS values were > 200 using tracer v.1.7.1 (Rambaut et al., 2018). We merged the independent runs using logcoMbiner (part of the BEAST v.2.6.3 package), and we summarized the post-burn-in species trees as a maximum clade credibility (MCC) tree in treeannotator (part of the BEAST v.2.6.3 package).
As supplementary analyses, we inferred totalevidence phylogenetic trees using the concatenation a p p r o a c h a s s u m i n g t h a t a l l g e n e l o c i a n d morphological characters have evolved under the same tree topology. We carried out a gene-tree discordance test to assess the convenience of the multispecies coalescent framework for our dataset. We performed likelihood-based tree topology tests to evaluate the agreement of the molecular data with three tree topologies: (1) the morphology-based systematics of Brassolini (Penz, 2007) suggesting that Narope is sister to remaining genera in the subtribe Brassolina; (2) the multispecies coalescent  Table S6).

Within-area cladogenesis and dispersal aMong bioregions
We estimated ancestral geographical ranges and colonization rates among bioregions using the dispersal-extinction-cladogenesis (DEC) model (Ree et al., 2005) as implemented in the R package biogeoBEARS v.1.1.2 (Matzke, 2013). We used the presence/absence biogeographical dataset and the four defined bioregions: (1) Mesoamerica plus NW Andes; (2) Atlantic Forest; (3) Amazonia; and the (4) dry diagonal. The dry diagonal was kept separate to estimate dispersal between the Atlantic Forest and Amazonia across a large area of less suitable environments for Brassolini. We did not constrain any a priori dispersal multiplier, nor did we stratify dispersal rates across the phylogeny. This allowed us to evaluate whether the data alone were in agreement with the main geological and palaeogeographical events in the Neotropics. To account for ancestral state uncertainty, we carried out 100 biogeographical stochastic mappings (BSMs; Dupin et al., 2017). To account for divergence time uncertainty, we used 100 random topologies from the BEAST 2 posterior distribution. To account for phylogenetic uncertainty of 25 missing species, we added such lineages randomly to their currently assigned genera in the 100 posterior phylogenies.
In this case, we assumed that phylogenetic diversity per genus was maximized, and we placed missing taxa in the crown genera given our comprehensive taxonomic sampling. All genera are monophyletic, but we joined Selenophanes and Catoblepia conservatively [posterior probability (PP) = 1.0; Fig. 2] when assigning missing species, because the monophyly of Catoblepia had only moderate support (PP = 0.75). We also carried out the biogeographical analyses using the taxonomically incomplete inferred species trees to rule out any bias in the approach. Variation in taxonomic opinion and effort across localities might be an important source of error in macroevolutionary analyses (Faurby et al., 2016). After analysing trends in the year of description of each Brassolini species and the year of the latest revision of their current taxonomic status (e.g. when former subspecies were raised to species), we did not observe any taxonomic biases among Neotropical bioregions in terms of taxon discovery, description and taxonomic opinion (i.e. oversplitting species vs. multiple infraspecific taxa; Supporting Information, Fig. S6). We therefore estimated the number of dispersal events and within-area cladogenesis events by simulating areas on nodes based on the 10 000 pseudoreplicated biogeographical histories (100 BSMs × 100 posterior species trees). We calculated colonization rates through time as c XtoY (t 1 ) = d XtoY (t 1 )/Br(t 1 ), where d XtoY (t 1 ) is the number of inferred dispersal events from area X to area Y in a 2 Myr interval, and Br(t 1 ) is the total length of all branches in the 2 Myr interval (Antonelli et al., 2018b). In addition, we calculated the relative number of speciation events in area X through time as λ X (t 1 ) = s X (t 1 )/L X (t 0 ), where s X (t 1 ) is the number of within-area cladogenesis events in a 2 Myr interval, and L X (t 0 ) is the number of inferred lineages in area X in the preceding 2 Myr interval (t 0 ), calculated as the cumulative sum of within-area cladogenesis events and dispersal events into area X minus local extinction (Xing & Ree, 2017).

speciation Within bioregions
For a second, independent estimate of within-area speciation, we calculated speciation rates for all tips and branches in the MCC species tree using BAMM v.2.5.0 (Rabosky, 2014). To account for the 25 missing species, we generated clade-specific sampling fractions at the genus level (Supporting Information, Table S7) under the assumption of random taxonomic sampling (FitzJohn et al., 2009). We used the R package baMMtools v.2.1.6  to estimate prior values for the evolutionary rate parameters, and we assigned the expected number of shifts to 1.0, which is a conservative value that minimizes type I errors and is recommended for phylogenies including ~100 species (Rabosky, 2014). We ran the analysis for 10 million generations, sampling event data every 1000 generations. After discarding the first 20% samples, we checked that the ESS value was > 200 and the likelihood estimates reached a stationary distribution.
We used the function 'dtRates' in baMMtools to compute mean speciation rates along branches in ~2 Myr intervals (tau parameter set to 0.05, because the root of the phylogeny was 37.9 Mya). We assigned the most probable inferred area from the 10 000 BSMs to the points where branches cross time intervals (Igea & Tanentzap, 2019). When the most probable area for a node and its parental one were the same, we assumed that the branch connecting both nodes occurred in such an area. When inferred areas for a node and its parental one were different, we assumed that the dispersal event occurred stochastically along the branch. We computed mean speciation rates only for cladogenesis events that most probably occurred in a single bioregion.  (Table 2). Our morphological dataset comprised 219 binary characters, 32 characters with three states, and four characters with four states. The multispecies coalescent model recovered four major clades with PP = 1.0, namely, clade A, Dynastor and Dasyophthalma; clade B, Opoptera and Narope; clade C, Caligopsis, Eryphanis and Caligo; and clade D, Penetes, Catoblepia, Selenophanes, Mielkella, Orobrassolis, Blepolenis and Opsiphanes.
The MCC species tree depicts the following phylogenetic relationships: Bia: (clade A: (clade B: (Brassolis: (clade C: clade D)))), with divergence times among these clades between 29 and 34 Mya (entire HPD, 20.24-43.64 Mya). However, the relationships among clades B, C, D and the genus Brassolis received low support (PP = 0.56-0.72; Fig. 2). The concatenation-based analyses produced similar divergence times among major Brassolini clades compared with the multispecies coalescent species tree (between 26 and 30 Mya; entire HDP, 20.38-39.22 Mya) but recovered higher posterior probabilities for Brassolini deep nodes when using the molecular dataset (0.97-1.00; Supporting Information, Fig. S2) and the total-evidence dataset (0.83-0.96; Supporting Information, Fig. S3). Nevertheless, the consensus and MCC tree topologies among all concatenation-based phylogenetic analyses remained similar, except for the genus Brassolis as sister to clade D, and clade C sister to these two. The likelihood-based tree topology tests could not resolve this main discrepancy between the multispecies coalescent and the concatenation-based species tree topologies (Supporting Information, Table  S6). The tree topology tests, however, clearly rejected the relationship of Narope as sister to the subtribe Brassolina suggested by previous cladistic analyses of morphological characters (P-value = 0.0003).
Importantly, the total-evidence dataset supported a posterior probability of 1.0 for the crown node of non-monotypic genera, regardless of the phylogenetic method, except for two cases: Narope paraphyletic with respect to the monotypic Aponarope, and Catoblepia paraphyletic with respect to Selenophanes (see Consequences for Brassolini systematics and classification in the Discussion section).

geographical range evolution and colonization rates
The inferred ranges derived from the 10 000 pseudoreplicated biogeographical histories and plotted against the MCC species tree in Figure 2 suggested that Brassolini most probably (relative probability, Pr = 0.50) originated in Amazonia 38 Mya (entire HPD, 27.68-46.97 Mya). The early rapid radiation of Brassolini between 29 and 34 Mya also probably took place in Amazonia (Pr = 0.44-0.67) ( Supporting Information, Fig. S7).
When taking into account phylogenetic and divergence time uncertainties by using the Bayesian posterior distribution of Brassolini species trees, dispersal out of Amazonia had the highest rates compared with other bioregions (Fig. 3). The colonization rate from Amazonia to Mesoamerica showed a sharp increase at ~2-4 Mya. Dispersal into Mesoamerica from other bioregions probably occurred throughout the Miocene, but increased significantly at ~6-8 Mya (Supporting Information, Fig. S8). Dispersal out of Mesoamerica and the NW Andes, however, was very low through time, and only since ~2 Mya has it increased exponentially towards Amazonia. Dispersal from Amazonia to the dry diagonal increased throughout the Plio-Pleistocene and since then has equalled dispersal to Atlantic Forest. Dispersal from Amazonia into the Atlantic Forest was higher in the Oligocene-Miocene transition at ~23-27 Mya than from/to any other Neotropical region. However, for most of the Miocene and Pliocene, dispersal into or out of the Atlantic Forest from/to any other bioregion was low overall through time, rendering the Atlantic Forest isolated; only since the Pleistocene ~2 Mya has dispersal into and from the dry diagonal increased until the present. The dispersal patterns remained similar in supplementary biogeographical analyses using the taxonomically incomplete inferred species trees (Supporting Information, Figs S7-S9).

regional speciation through tiMe
The relative rates of within-area cladogenesis estimated in b i o g e o b e a r s and speciation estimated in BAMM were congruent (Fig. 4). Speciation rates increased through time in all areas except in the dry diagonal, which has had near zero endemic cladogenesis. Diversification in the Atlantic   Forest might have begun during the Oligocene-Miocene transition at ~20-24 Mya. During that time, the cladogenesis rate was slightly > 0.5 (Fig.  4A) and speciation rate began to increase (Fig. 4B), probably driven by the split of Penetes and, to a lesser extent, the stem lineage of Dasyophthalma. However, there seemed to be a peak in cladogenesis a t ~8 -1 0 M y a , s p e c i f i c a l l y i n t h e l i n e a g e s comprising extant Dasyophthalma, Orobrassolis and Blepolenis. Afterwards, diversification in the Atlantic Forest remained constant throughout most of the Miocene until ~2-4 Mya, when it decreased slightly towards the present. In contrast, speciation in Mesoamerica and the NW Andes began later, at ~8-12 Mya, remained constant during most of the Miocene and Pliocene, and increased sharply only at ~2-4 Mya, driven by the origin of extant endemic lineages within Caligo, Opsiphanes and Eryphanis. Speciation in Amazonia was higher episodically than in the Atlantic Forest and Mesoamerica throughout most of the Miocene, but it also increased dramatically during the Pleistocene (~2 Mya) by the split of several extant species within Bia, Eryphanis and Catoblepia. At present, the speciation rate in Mesoamerica and the NW Andes is similar to the Amazonian lineages, with both apparently being higher than in the Atlantic Forest.

DISCUSSION
Here, we investigated the evolutionary history of the butterfly tribe Brassolini and studied the role of regional speciation and dispersal to understand the origin of extant species across Neotropical rainforest biomes (Fig. 5). By using molecular and morphological datasets and the multispecies coalescent framework, we propose a revised hypothesis of relationships among genera and make three changes to the classification of Brassolini at the tribal and genus levels.

consequences for brassolini classification
Previously, Brassolini has been subdivided into three subtribes based on shared morphological similarities: Biina (monotypic), Naropina (Narope and Aponarope) and Brassolina (remaining Brassolini genera) (Casagrande, 1995;Penz, 2007). Nonetheless, here Narope is sister to Opoptera (clade B in Fig. 2) with strong statistical support, and this clade is nested within the subtribe Brassolina (sensu Casagrande, 1995). These relationships have not been recovered in any cladistic analysis of Brassolini to date. At the genus level, the main discrepancy among the multispecies coalescent and the concatenation-based phylogenies was in the phylogenetic position of the genus Brassolis. Morphological character states of adult and early stages in this genus are highly polymorphic and divergent in comparison to other Brassolini genera, thus, obscuring the phylogenetic relationships within Brassolis and among this and other lineages (Garzón-Orduña & Penz, 2009;Penz et al., 2013). Nonetheless, the six molecular markers aided in resolution of the apparent long-branch attraction between Brassolis and Dynastor seen in morphological phylogenies (Penz, 2007) and grouped Brassolis with clades C and D (Caligo and Opsiphanes groups) in both the multispecies coalescent and concatenation-based phylogenies.
We propose the following changes to the Brassolini classification (registered in ZooBank, LSID: urn:lsid:zoobank.org:pub:8B7AB3B5-79A6-464F-A5D9-99960FCEEC46): 1. The subtribe Naropina Stichel, 1925   Coloured directional arrows represent dispersal from source areas, following colours in the key. The width and shape of arrows represent the estimated relative colonization rate between two areas in biogeobears, as depicted in Figure 3. Bioregions were coloured to represent the estimated relative cladogenesis/speciation rates qualitatively, as either the centre of speciation or the area with comparatively lower speciation rates per time window.
probability (0.99) in our analysis, which justifies the revised generic assignment. Using COI barcoding sequences, Shirai et al. (2017) also found that the genus Catoblepia Stichel, 1902 was paraphyletic with respect to Selenophanes owing to the sister relationship between S. orgetorix and Selenophanes species.

species diversification and dispersal
One-third of Brassolini species are endemic to rainforests in Mesoamerica, the NW slope of the Andes including the Chocó, and the Brazilian Atlantic Forest. Rapid butterfly speciation has probably occurred in Mesoamerica and the NW Andes since the Pleistocene, suggesting that it acted as a cradle of diversity. In contrast, low and constant diversification since at least the early Miocene has shaped Brassolini endemic diversity in the Atlantic Forest. This finding, coupled with the survival of old lineages found mainly in montane habitats (e.g. Penetes, Dasyophthalma and Orobrassolis), suggests that the Atlantic Forest acts as a museum of diversity for Brassolini. Our results also suggest that Brassolini diversified mainly in Amazonia through episodic increases in cladogenesis rate throughout the Miocene from ~20 Mya (Fig.  4A), and such diversification events might have been accompanied by increasing dispersal out of this bioregion. The diversification and biogeographical history of Brassolini are consistent with the source, museum and cradle of diversity scenarios previously proposed for the Amazon rainforest. Amazonia has been the main source of Neotropical vertebrate and plant species diversity by means of dispersal during the past 60 Myr (Antonelli et al., 2018b). Rapid species diversification of some butterfly groups has occurred there since at least the Oligocene (e.g. Matos-Maraví et al., 2013;Chazot et al., 2016;Toussaint et al., 2019), and old lineages might have accumulated gradually in Amazonia since that time (Peña et al., 2010;Condamine et al., 2012). However, more recent diversification episodes supporting Amazonia as an alternating museum and cradle of diversity could have been triggered by palaeoenvironmental changes during the Miocene (Antonelli et al., 2009;Antonelli & Sanmartín, 2011;Chazot et al., 2019b).
The establishment of the Pebas wetland system in western Amazonia during most of the Miocene (~8-24 Mya) (Wesselingh & Salo, 2006) has probably hindered dispersal from Amazonia towards the NW Andes and the Pacific in some butterfly groups (Elias et al., 2009;De-Silva et al., 2016Chazot et al., 2019b). In the case of Brassolini, the Pebas wetland system has also probably hindered dispersal out of Amazonia (Fig. 3), but it might not have limited diversification drastically (Fig. 4). This wetland system was occasionally affected by marine flooding events, which were probably short lived, and at least two events, dated at ~13.7 and ~17.8 Mya (Jaramillo et al., 2017b), were concurrent with the early and rapid cladogenesis events of major Brassolini clades in Amazonia. The Pebas wetland might have increased opportunities for dispersal across Amazonia and the diversification of some typical coastal plants (Bernal et al., 2019), including many of the hostplants of Brassolini in the region: Arecaceae (palms), Poaceae (bamboos) and Zingiberales, such as Marantaceae (Beccaloni et al., 2008;Janzen et al., 2009). It is plausible that the evolutionary consequences of the Miocene marine incursions in western Amazonia might have been mixed and rather ephemeral (Musher et al., 2019), restricting dispersal and diversification in some butterfly lineages, such as Ithomiini (Chazot et al., 2019b), while providing ecological opportunities for diversification of others, such as Brassolini.
The increase in diversification of Amazonian lineages in the Pleistocene might be related to younger palaeoenvironmental changes, such as fluctuating river dynamics mediated by Quaternary glacial cycles (Ribas et al., 2012). For example, biogeographical structure in the genus Bia seems to be associated with major Amazon Basin rivers (Penz et al., 2017), suggesting that riverine barriers might have been involved in allopatric differentiation in this genus. Indeed, in some areas surrounding the Amazon River, the distance between river banks is much larger than the estimated daily dispersal of species such as Bia actorion (Tufto et al., 2012;Penz et al., 2015). Other hypotheses for the diversification of Neotropical butterflies in the Pleistocene have been put forward, such as rainforests acting as refugia (Garzón-Orduña et al., 2014), but this has not yet been corroborated in a macroevolutionary framework (Matos-Maraví, 2016). A phylogeographical evaluation including more samples per species would be necessary to estimate the influence of river barriers or rainforest refugia in differentiating Brassolini populations and, potentially, to explain the apparent increase in diversification of Amazonian butterflies during the Pleistocene.
Miocene expansions of wet forests between the Atlantic Forest and Amazonia across the dry diagonal might have been important determinants of ancient Brassolini diversity. The onset of a period of global cooling at ~10-15 Mya probably drove the expansion of open environments (Simon et al., 2009;Azevedo et al., 2020), but intermittent forest connections between Amazonia and the Atlantic Forest have been documented based on palynological data and biogeographical studies (e.g. Costa, 2003;Werneck, 2011;Prates et al., 2016;Trujillo-Arias et al., 2017). Although evidence of repeated Neogene connection between the Atlantic Forest and Amazonia is still scarce and restricted to phylogeographical studies (e.g. Ledo & Colli, 2017;Capurucho et al., 2018), the origin of the modern transcontinental Amazon River by ~10 Mya might have facilitated dispersal of rainforest taxa (Albert et al., 2018). Indeed, we found that late Miocene dispersal rates between these wet forest biomes increased slightly through time. In contrast, the Atlantic Forest has recently become a source of diversity for the present dry diagonal, because Brassolini dispersal increased sharply during the Pleistocene, in agreement with previous phylogeographical evidence (Costa, 2003;Batalha-Filho et al., 2013;Thomé et al., 2016). We suggest that episodic Miocene wet forest expansions might have represented ecological opportunities for Brassolini to disperse across Central Brazil, either by increasing the sizes of wet forest galleries in the dry diagonal biome or by fully connecting Amazonia and Atlantic Forest.
By the Plio-Pleistocene, speciation in Mesoamerica and the NW Andes seems to have increased and surpassed the speciation rate in the Atlantic Forest. Biotic dispersal into Mesoamerica from South America (e.g. Mullen et al., 2011;Bacon et al., 2015) and butterfly dispersal from Central America to other Neotropical areas (Toussaint et al., 2019) is likely to have occurred throughout most of the Miocene, facilitated by land availability, with uranium-lead dating confirming that segments of the Panama arc had already emerged by 13-15 Mya (Montes et al., 2015). However, only at the end of the Pliocene ~3-4 Mya did shallow seawaters recede fully from the region (Coates et al., 2004;Montes et al., 2012;O'Dea et al., 2016;Jaramillo et al., 2017a), thus increasing opportunities for terrestrial species diversification in lowland areas that were sporadically inundated by seawater during the Miocene (Hosner et al., 2016). In contrast, speciation in the Atlantic Forest remained low, and dispersal out of this region decreased until the mid-Pleistocene. The contraction of the Atlantic Forest since the mid-Miocene (Simon et al., 2009;Edwards et al., 2010;Strömberg, 2011;Werneck, 2011) might have restricted Brassolini speciation to montane habitats in south-eastern Brazil, where most extant species occur, while old lineages having low speciation rates through time survived in remnants of rainforest.

future directions
The evolutionary patterns inferred here indicate that geography has been an important determinant of the diversity and distribution of extant Brassolini butterflies. Other unmeasured or unknown traits might have driven the observed disparate trends in speciation through time; for example, larval hostplants. However, many Mesoamerican Brassolini feed on the same plant families (e.g. Arecaceae) and even genera (e.g. Euterpe and Bactris) as do Atlantic Forest taxa (Beccaloni et al., 2008;Janzen et al., 2009;Robinson et al., 2010). Thus, despite the limited data at hand, it seems that larval hostplants might have contributed little to the disparate diversification dynamics of Brassolini.
We are aware that the estimation of changes in diversification rate through time is challenging (Moore et al., 2016). However, further evaluation via simulated data suggested that BAMM can infer speciation rates accurately from extant-taxon phylogenies (Rabosky et al., 2017). Although low extinction rates might explain the museum of diversity pattern, we have focused conservatively only on regional species origination rates to explain lineage accumulation within bioregions. Furthermore, Louca & Pennell (2020) stressed that the identifiability of diversification processes is problematic in a model-selection approach, but this remains to be evaluated in methods that use rate shift models, such as BAMM (Laudanno et al., 2021). Reassuringly, we also used a fundamentally distinct measure of within-area cladogenesis (Xing & Ree, 2017) that corroborated the pattern of recent and rapid speciation in Mesoamerica and the NW Andes, and an older lineage accumulation and lower diversification in the Atlantic Forest. However, our proposed checkered assemblage of Neotropical cradles and museums certainly needs to be evaluated with other larger phylogenies of rainforest lineages, particularly in light of the differences in natural history among various butterfly groups.
We defined bioregions that fitted the distribution and composition of Brassolini species by using georeferenced occurrences, rather than by assuming any subjective criteria to adopt areas from hierarchical bioregion classifications (Vilhena & Antonelli, 2015). This data-driven approach revealed a clustering of butterfly communities from southern North America to Mesoamerica and the NW side of the Andes, which conforms with the likely emergence of landmass in Central America, Chocó and northeastern Colombia after the collision of the Panama Block and northwestern South America by the late Oligocene (Coates & Stallard, 2013;Jaramillo, 2018). Furthermore, the inference of large bioregions for Brassolini agrees with its apparent poor spatial structuring and high dispersal capabilities. For instance, given a suitable environment, ancestral Bia had the capability to disperse across pan-Amazonia in only 1463-3115 years (Penz et al., 2015). The delimitation of Mesoamerica plus the NW Andes and the Atlantic Forest as separate bioregions has probably been determined by ecological and habitat suitability rather than by geographical barriers to Brassolini dispersal; this being a similar biogeographical pattern found in other Neotropical rainforest taxa (Dexter et al., 2017;Pérez-Escobar et al., 2019).

conclusions
We investigated whether the diversification and biogeographical history of butterflies classified in the tribe Brassolini conforms to the premise that the Neotropics consists of a regional network of museums and cradles of diversity (Fig. 5). We found that species endemic to Mesoamerica plus the NW side of the Andes have originated by old dispersal events into the region, and that speciation rates increased only at ~2 Mya (cradle of diversity). In contrast, species endemic to the Brazilian Atlantic Forest probably arose under low and continuous speciation rates since ≥ 12 Mya (museum of diversity) and were likely to be connected intermittently to Amazonia since the Miocene. The dynamic evolutionary history of Amazonian Brassolini lineages, alternately diversifying and accumulating old lineages, appears to reflect palaeoenvironmental changes, including landscape reconfigurations in western Amazonia and palaeoclimate fluctuations during the Neogene and Quaternary. By focusing on a representative lineage of a major insect order, our study paves the way for understanding why the Neotropics should be considered both a museum and a cradle of diversity.

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article at the publisher's website: Figure S1. Map of the Neotropics showing the 7174 cleaned occurrences of Brassolini species and the 81 flagged as potential errors, which were excluded from the bioregion delimitation analyses. Figure S2. Single-gene, concatenated multi-locus and morphological trees. Each tree represents the 50% majorityrule consensus of 7500 posterior trees inferred in Mrbayes v.3.2.6. Posterior probabilities are shown on every node. A, CAD gene tree. B, COI gene tree. C, EF1α gene tree. D, GAPDH gene tree. E, RpS5 gene tree. F, wingless gene tree. G, concatenated multi-locus molecular tree. H, morphology-based tree. Figure S3. Total-evidence consensus tree using the concatenation approach. The tree represents the 50% majority-rule consensus of 7500 posterior trees inferred in Mrbayes v.3.2.6. Posterior probabilities are shown on every node. Figure S4. Strict consensus tree of the 13 equally parsimonious trees estimated using maximum parsimony in TNT v.1.5. Numbers next to nodes represent the contribution of genes as measured by partitioned Bremer support. The scores correspond to the genes CAD, COI, EF1α, GAPDH, RpS5 and wingless, respectively. Figure S5. Maximum clade credibility species tree using the total-evidence dataset and the multispecies coalescent model in BEAST v.2.6.3. Posterior probabilities and the 95% highest posterior density (HPD) intervals are shown on every node. The species tree is calibrated in millions of years. Calibration points, based on the study by Chazot et al. (2019a), are indicated by red arrows: (1) Figure S6. Brassolini taxonomic resolution across Neotropical bioregions. We found similar trends across Mesoamerica, Amazonia and Atlantic Forest in collecting, describing and revising the species status of all valid Brassolini taxa. Thus, we rule out any taxonomic bias affecting our biogeographical and diversification analyses. A-C, first, we compiled the year of description of every species restricted to Amazonia (A), Mesoamerica and NW Andes (B) and the Atlantic Forest (C). D-F, second, we compiled the year of the last revision of the specific status of Brassolini taxa restricted to Amazonia (D), Mesoamerica and NW Andes (E) and the Atlantic Forest (F). Vigorous taxonomic activity has occurred simultaneously in the three main rainforest biomes in the Neotropics during the mid-19 th century to the first quarter of the 20 th century, and during the late 20 th century to the present. Figure S7. Ancestral range probability based on the dispersal-extinction-cladogenesis (DEC) model and 10 000 biogeographical stochastic mappings, plotted against the maximum clade credibility (MCC) species tree of Brassolini. A, the most probable state accounting for missing species in the posterior phylogenies is plotted on every node. B, the probabilities of ancestral ranges accounting for missing species in the posterior phylogenies are shown as pie charts on every node. C, the most probable state using the sampled species trees is plotted on every node. D, the probabilities of ancestral ranges using the sampled species trees are shown as pie charts on every node. Bioregions were coded as follows: C, South American dry diagonal; F, Brazilian Atlantic Forest; M, Mesoamerica and Chocó; S, Amazonia. Figure S8. Dispersal rates through time calculated with 10 000 biogeographical stochastic mappings in biogeobears. Bioregions were coded as follows: C, South American dry diagonal; F, Brazilian Atlantic Forest; M, Mesoamerica and Chocó; S, Amazonia. The x-axis in every chart is on a million years scale. The y-axis represents the estimated dispersal rates (as events per lineage per million years) using the formula given by Antonelli et al. (2018b). 'rate.StoM', for example, is dispersal from source area 'S' to target area 'M'. Continuous lines are the median values; dark green ribbons represent the lower and upper quartiles (0.25 and 0.75 quantiles), light green ribbons the 0.1 and 0.9 quantiles, and dashed lines the 0.05 and 0.95 quantiles. A, estimates accounting for missing species in the posterior phylogenies. B, estimates based on the sampled species trees. Figure S9. Within-area cladogenesis through time calculated with 10 000 biogeographical stochastic mappings in biogeobears. Bioregions were coded as follows: C, South American dry diagonal; F, Brazilian Atlantic Forest; M, Mesoamerica and Chocó; S, Amazonia. The x-axis in every chart is on a million years scale. The y-axis represents the estimated relative number of cladogenesis events per million years using a formula modified from Xing & Ree (2017). 'rate.StoS', for example, is relative in situ cladogenesis in area 'S'. Continuous lines are the median values; dark green ribbons represent the lower and upper quartiles (0.25 and 0.75 quantiles), light green ribbons the 0.1 and 0.9 quantiles, and dashed lines the 0.05 and 0.95 quantiles. A, estimates accounting for missing species in the posterior phylogenies. B, estimates based on the sampled species trees. Table S1. Voucher locality information and associated genetic data deposited in GenBank or BOLD (ASARD codes). Table S2. The best-fitting partition scheme for the molecular dataset estimated by the program partitionfinder v.2.1.1. Table S3. The best-fitting partition scheme for the morphological dataset estimated by homoplasy scores (f) in the program TNT v.1.5. Table S4. Bayes factor comparison between the strict and relaxed clock models. Abbreviation: SSML, steppingstone marginal likelihood. The Bayes factor (BF) was calculated as twice its natural logarithm (2 loge BF) and, to account for the number of parameters (NP), we summed to this value the following: (NPrelaxed − NPstrict) × loge 0.01. The relaxed clock model was strongly preferred over the strict clock for all loci (BF > 10). Table S5. Bayes factor (BF) comparisons among six models involving two tree models (Yule and birth-death) and molecular clock partitions (one, two, mitochondrial and nuclear, and six, for each gene partition). Marginal L (likelihood) estimates were based on 25 path-sampling steps under thermodynamic integration. The BF with respect to the highest-likelihood model (Yule + 2 clocks) was calculated as twice its natural logarithm (2 loge BF). The Yule + 2 molecular clocks model received decisive support (BF > 10). Table S6. Tree topology test comparing the branching order of early divergent Brassolini lineages. The total-evidence concatenation-based consensus tree (Conc) has the following constraint topology (Bia: ((Dasyophthalma:Dynastor): ((Narope:Opoptera): ((Caligogroup): ((Brassolis: Opsiphanes-group)))))), which is highly similar to the unconstrained molecular tree (Unconstr). The multispecies coalescent tree (MSC) has the following constraint topology (Bia: ((Dasyophthalma:Dynastor): ((Narope:Opoptera): (Brassolis: (Caligo-group: Opsiphanes-group))))). The morphology-based systematics of Brassolini (Naropina) has the following constraint topology (Bia: (Narope: remaining Brassolina genera)). Abbreviations: bp-RELL, bootstrap proportion using the RELL method; c-ELW, expected likelihood weight; DeltaL, likelihood difference from the maximum likelihood tree topology; p-AU, P-value of approximately unbiased test; p-KH, P-value of one-sided Kishino-Hasegawa test; p-SH, P-value of Shimodaira-Hasegawa test; p-WKH, P-value of weighted Kishino-Hasegawa test; p-WSH, P-value of weighted Shimodaira-Hasegawa test. All tests performed 10 000 resamplings using the RELL method. The plus signs denote the 95% confidence sets. The minus signs denote significant exclusion. The tested branching orders are undecidable given our molecular dataset, except for Naropina, which is significantly rejected. Table S7. Sampling fractions used for taking into account missing species in the species tree for calculation of dispersal and speciation rates. The revised genera are monophyletic, and, given our comprehensive taxonomic sampling, we assumed that missing taxa are within crown nodes.

SHARED DATA
All the datasets and scripts used in this study can be found in TreeBASE (ID 25040) and Zenodo (DOI: http://doi. org/10.5281/zenodo.4500393).