Materials and Methods
Additional data and methodology are described in SI Appendix .
Data Compilation. We extracted Neotropical clades that included at least 80% of the species distributed in the Neotropics as defined in (Antonelli et al. 2018c). From large-scale phylogenies, we relied on the time-calibrated phylogenies of frogs and toads (Hutter et al. 2017), salamanders (Pyron et al. , 2013; Pyron, 2014), Squamata (lizards and snakes) (Pyron & Burbrink 2014), birds (Jetz et al. 2012) (including only species for which genetic data was available), mammals (Bininda-Emonds et al. , 2007; Kuhn et al. , 2011), and plants (Zanne et al. 2014). In addition, phylogenies of particular lineages not represented in the global trees (or with improved taxon sampling) were obtained from published studies (SI Appendix , Tables S1-5) or reconstructedde novo (for caviomorph rodents, including 199 species, and phyllostomid bats, 170 species; see SI Appendix ). However, whenever possible, we preferred to extract phylogenies from a single dated tree rather than performing a meta-analysis of individual trees coming from different sources (e.g., ref. (Hoorn et al.2010; Jansson et al. 2013)), such that divergence times would be comparable. In addition, this procedure allows reducing the bias of focusing on particular taxonomic levels (i.e., individual studies often focus on genera) and thus comparing lineages of similar ages (Wiens 2017). When the most inclusive Neotropical clade was over the family rank, or for extremely species-rich groups, like the family Solanaceae with more than 2,700 species, we split the phylogenies into smaller clades (e.g., tribes) to take into account the heterogeneity of evolutionary histories likely characterized by different morphologies and key adaptations.
Diversification Analyses. We designed and applied a series of 10 birth-death diversification models estimating speciation(λ) and extinction (μ) rates for each of the 150 phylogenies with the R-package RPANDA 1.3 (Morlon et al. 2016). We first fitted a constant-rate birth–death model (null model) and a set of three models in which speciation and/or extinction vary according to time (Morlon et al. 2011): λ(t) and μ(t) . For time-dependent models, we measured rate variation for speciation and extinction with α and β , respectively: α andβ > 0 reflect decreasing speciation and extinction toward the present, respectively, while α and β< 0 indicate the opposite, increasing speciation and extinction toward the present.
We further investigated the effect of environmental changes, here approximated by mean global temperatures and paleo-elevations of the Andes (the main factors hypothesized behind the origin of Neotropical diversity). Temperature variations during the Cenozoic were obtained from global compilations of deep-sea oxygen isotope (δ18O) from ref. (Prokoph et al. 2008; Zachoset al. 2008), but also on other different paleotemperature curves, to assess the impact of paleotemperature uncertainty on our results (SI Appendix ). For Andean paleo-elevations we retrieved a generalized model of the paleoelevation history of the tropical Andes, compiled from several references (ref. (Lagomarsino et al. 2016) and references therein)). We examined whether speciation and/or extinction correlate with one of these variables using paleoenvironment-dependent diversification models (Condamine et al. 2013), which extends time-dependent models to account for potential dependencies between diversification rates and measured environmental variables. The elevation of the Andes caused dramatic modifications in Neotropical landscapes and could have impacted indirectly the diversification of non-Andean groups. Therefore, we decided to apply uplift models to all clades in our study, independently on whether their distribution is centered in the Andes or not. We fitted three models in which speciation and/or extinction vary continuously with temperature changes (λ(T) and μ(T)) , and three others with the elevation of the Andes (λ(A) and μ(A)) . In this case,λ00) is the expected speciation (extinction) rate under a temperature of 0°C (or a paleoelevation of 0 m for the uplift models). For the groups in which diversification was temperature- or uplift-dependent, we analyzed whether the speciation (α) and extinction (β) dependency was positive or negative. For temperature models, α (β ) > 0 reflect increasing speciation (extinction) with increasing temperatures, and conversely. For the uplift models, α(β ) > 0 reflect increasing speciation (extinction) with increasing Andean elevations, and conversely. We accounted for missing species for each clade in the form of sampling fraction(ρ) (Morlon et al. 2011) and assessed the strength of support of the models by computing Akaike information criterion (AICc), ∆AICc, and Akaike weights (AICω) to select the best-fit model.
Finally, we tested the effect of clade age, size and sampling fraction on the preferred model using Kruskal-Wallis rank sum test and performed multiple pairwise-comparison between groups with corrections for multiple testing using Wilcoxon rank sum test. We also evaluate if the estimated diversification trends (i.e. , increasing, constant, decreasing) differ between lineages characterized by different geographic distributions, altitudinal ranges and habitat preferences (see SI Appendix ).
Our study focuses on the effects of environmental factors over diversity, and in this sense we did not directly assess the effect of biotic factors. Only few models of that type exist, and none explicitly incorporate species interactions within clades. The most relevant of such models are diversity-dependent models where the number of extant species in a clade affects diversification rates (Etienne et al.2012). We did not fit these models for three additional reasons:(i) it is not straightforward to compare diversity-dependent models with other diversification models on the basis of AICc (Etienneet al. 2016); (ii) the sampling scheme currently implemented in diversity-dependent models assumes that exactly nspecies are sampled (n -sampling), while other models assume that each species is sampled with a fixed probability (ρ-sampling), and likelihoods associated to these two sampling schemes are not directly comparable (Stadler 2009; Lambert 2017); and (iii) previous simulations have shown that these models tend to over-fit the data when speciation decreases toward the present regardless the cause of decline (Etienne et al. 2016; Condamine 2018). This latter artifact might affect more than 78% of clades in our study, since we found most phylogenies are characterized by constant or decreasing rates through time (see Results ). Nonetheless, we acknowledge that the role of biotic interactions in explaining Neotropical diversity remains to be tested.
Acknowledgments. We thank all researchers who shared their published data through databases or to us directly (Drs. Arevalo, Martins, Fortes Santos, Simon, Lohmann, Mendoza, Swenson, Erkens, van der Meijden, and Freitas), Dr. Sanmartín and Dr. Manzano for invaluable comments on the manuscript and analyses. This work was funded by an “Investissements d’Avenir ” grant managed by the Agence Nationale de la Recherche (CEBA, ref. ANR-10-LABX-25-01) and by the ANR GAARAnti project (ANR-17-CE31-0009). A.A. is supported by the Knut and Alice Wallenberg Foundation, the Swedish Research Council, the Swedish Foundation for Strategic Research, and the Royal Botanic Gardens, Kew.