Deciphering the evolution of birdwing butterflies 150 years after Alfred Russel Wallace

One hundred and fifty years after Alfred Wallace studied the geographical variation and species diversity of butterflies in the Indomalayan-Australasian Archipelago, the processes responsible for their biogeographical pattern remain equivocal. We analysed the macroevolutionary mechanisms accounting for the temporal and geographical diversification of the charismatic birdwing butterflies (Papilionidae), a major focus of Wallace’s pioneering work. Bayesian phylogenetics and dating analyses of the birdwings were conducted using mitochondrial and nuclear genes. The combination of maximum likelihood analyses to estimate biogeographical history and diversification rates reveals that diversity-dependence processes drove the radiation of birdwings, and that speciation was often associated with founder-events colonizing new islands, especially in Wallacea. Palaeo-environment diversification models also suggest that high extinction rates occurred during periods of elevated sea level and global warming. We demonstrated a pattern of spatio-temporal habitat dynamics that continuously created or erased habitats suitable for birdwing biodiversity. Since birdwings were extinction-prone during the Miocene (warmer temperatures and elevated sea levels), the cooling period after the mid-Miocene climatic optimum fostered birdwing diversification due to the release of extinction. This also suggests that current global changes may represent a serious conservation threat to this flagship group.

simplified likelihood interpretation of the BayArea model that replicates BayArea's cladogenesis assumption (which postulates that nothing in particular happens at cladogenesis). The purpose of having a BayArea-like model was to test the potential importance of the cladogenesis model on particular datasets.
The three pairs of models were compared to see whether the inclusion of founderevent speciation (parameter J) significantly improved likelihood using AICc. We also used AICc and LRT to select the best-fit model out of the six that were implemented. The corresponding best-fit model was then used to make the final biogeographical reconstruction shown in Fig. 2 and Supplementary Fig. S3.
To visualize the tempo and mode of diversification of the group, we first reconstructed lineages-through-time (LTT) plots for the whole tree and for the genera Ornithoptera and Troides ( Supplementary Fig. S4).

Trait-dependent diversification.
If there are a priori reasons to believe that features of a species (e.g. ecological or morphological traits, geographical locations) influence diversification, it is relevant to use trait-dependent diversification models that concurrently model trait evolution and its impact on diversification 30,31 . In those models, species are characterized by an evolving trait, and follow a birth-death process in which speciation and extinction rates depend on the value of the trait at that time. Several models derived from the Binary State Speciation and Extinction model (BiSSE,31 ) have been used in various biogeographic contexts, e.g. the Geographic State Speciation and Extinction (GeoSSE,32 ). However, those models have rarely been applied in island biogeography.
Here we used 500 trees randomly selected from the Bayesian posterior distribution of dating analyses and we combined them with distributional data categorized in two traits: species occurring extra-Wallacea (trait A), and species occurring within Wallacea (trait B) [species can be widespread, i.e. occur within the two regions, and are coded AB]. We then fitted 12 GeoSSE diversification models to estimate the speciation rate for trait A (sA) and B (sB), the allopatric speciation rate between the two regions (sAB), the extinction rate (xA / xB), and the transition rate between trait A/B and B/A (denoted qAB / qBA).
Analyses were performed using the R-package diversitree 0.7-6 33 using the functions make.geosse to construct the likelihood functions for each model based on the data, and the functions constrain and find.mle to apply different diversification scenarios as detailed above.
For each GeoSSE model, we computed the AICc. We checked support for the selected model against all models nested within it using the Likelihood Ratio Test (LRT, significant at P < 0.05). The scenario supported by LRT and with the lowest AICc was considered the best.
If the model with the lowest AICc was not supported by LRT, the model with the least parameters was considered the best. Net diversification rates (r=λ -µ) were computed from speciation and extinction rate estimates. Finally we used the maximum clade credibility (MCC) tree and a MCMC approach for the best model to examine the credibility intervals of the parameter estimates. We used an exponential prior 1/(2r) to be as conservative as possible 33 and started the chain with the parameters obtained by maximum likelihood. We ran 20,000 steps of MCMC and applied a burn-in of 2,000 steps.

Palaeo-environment-dependent diversification.
To test the effect that past environment had on the diversification of birdwings, we used a         Table S1). (c)