Molecular Phylogenetics and Evolution 178 (2023) 107635 Contents lists available at ScienceDirect Molecular Phylogenetics and Evolution journal homepage: www.elsevier.com/locate/ympev Slow and steady wins the race: Diversification rate is independent from body size and lifestyle in Malagasy skinks (Squamata: Scincidae: Scincinae) Francesco Belluardo a, b, c, *, A. Jesus Muñoz-Pajares a, c, d, Aurélien Miralles e, Daniele Silvestro f, g, h, Walter Cocca a, b, c, Fanomezana Mihaja Ratsoavina i, Andrea Villa j, Sam Hyde Roberts k, l, m, Marcello Mezzasalma a, c, n, Alexander Zizka o, Alexandre Antonelli g, h, p, q, Angelica Crottini a, b, c a CIBIO, Centro de Investigação em Biodiversidade e Recursos Genéticos, InBIO Laboratório Associado, Universidade do Porto, Campus de Vairão, Rua Padre Armando Quintas, 4485-661 Vairão, Portugal b Departamento de Biologia, Faculdade de Ciências, Universidade do Porto, Praça Gomes Teixeira, 4099-002 Porto, Portugal c BIOPOLIS Program in Genomics, Biodiversity and Land Planning, CIBIO, Campus de Vairão, Rua Padre Armando Quintas, 4485-661 Vairão, Portugal d Departamento de Genética, Universidad de Granada, Avenida de la Fuente Nueva S/N, 18071 Granada, Spain e Institut de Systématique, Evolution, Biodiversité (ISYEB), Muséum National d’Histoire Naturelle, CNRS, Sorbonne Université, EPHE, Université des Antilles, Paris, France f Department of Biology, University of Fribourg and Swiss Institute of Bioinformatics, Ch. du Musée 10, 1700 Fribourg, Switzerland g Department of Biological and Environmental Sciences, University of Gothenburg, Box 461, 405 30 Göteborg, Sweden h Gothenburg Global Biodiversity Centre, Box 461, 405 30 Göteborg, Sweden i Mention Zoologie et Biodiversité Animale, Domaine Sciences et Technologies, Université d’Antananarivo, B.P. 906, 101 Antananarivo, Madagascar j Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Edifici ICTA-ICP, 08193 Cerdanyola del Vallès, Spain k SEED Madagascar, Unit 7, Beethoven Street 1A, W10 4LG London, UK l Oxford Brookes University, Headington Campus, 0X3 0BP Oxford, UK m Operation Wallacea, Wallace House, Old Bolingbroke, PE23 4EX Spilsby, UK n Department of Biology, Ecology and Earth Sciences, University of Calabria, Via P. Bucci 4/B, 87036 Rende, Italy o Department of Biology, Philipps-University Marburg, Karl-von-Frisch-Straße 8, 35043 Marburg, Germany p Royal Botanic Gardens, Kew, TW9 3AE Richmond, UK q Department of Biology, University of Oxford, South Parks Road, OX1 3RB Oxford, UK A R T I C L E I N F O A B S T R A C T Keywords: Most of the unique and diverse vertebrate fauna that inhabits Madagascar derives from in situ diversification from Biogeography colonisers that reached this continental island through overseas dispersal. The endemic Malagasy Scincinae lizards are Evolutionary radiations amongst the most species-rich squamate groups on the island. They colonised all bioclimatic zones and display many Fossoriality ecomorphological adaptations to a fossorial (burrowing) lifestyle. Here we propose a new phylogenetic hypothesis for In situ diversification their diversification based on the largest taxon sampling so far compiled for this group. We estimated divergence Islands Madagascar times and investigated several aspects of their diversification (diversification rate, body size and fossorial lifestyle evolution, and biogeography). We found that diversification rate was constant throughout most of the evolutionary history of the group, but decreased over the last 6–4 million years and independently from body size and fossorial lifestyle evolution. Fossoriality has evolved from fully quadrupedal ancestors at least five times independently, which demonstrates that even complex morphological syndromes – in this case involving traits such as limb regression, body elongation, modification of cephalic scalation, depigmentation, and eyes and ear-opening regression – can evolve repeatedly and independently given enough time and eco-evolutionary advantages. Initial diversification of the group likely occurred in forests, and the divergence of sand-swimmer genera around 20 Ma appears linked to a period of aridification. Our results show that the large phenotypic variability of Malagasy Scincinae has not influenced diversification rate and that their rich species diversity results from a constant accumulation of lineages through time. By compiling large geographic and trait-related datasets together with the computation of a new time tree for the group, our study contributes important insights on the diversification of Malagasy vertebrates. * Corresponding author at: CIBIO, Centro de Investigação em Biodiversidade e Recursos Genéticos, InBIO Laboratório Associado, Universidade do Porto, Campus de Vairão, Rua Padre Armando Quintas, 4485-661 Vairão, Portugal. E-mail address: france89belluardo@gmail.com (F. Belluardo). https://doi.org/10.1016/j.ympev.2022.107635 Received 23 March 2022; Received in revised form 18 September 2022; Accepted 29 September 2022 Available online 5 October 2022 1055-7903/© 2022 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by- nc-nd/4.0/). F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 1. Introduction vertebrate assemblage at the end of the Cretaceous, when Madagascar was already isolated, has likely provided new colonisers with a depleted Island biotas have played a significant role in the formulation of and much less-competitive environment where new lineages could evolutionary theory (e.g., Darwin and Murray, 1859; Simpson, 1953; widely diversify (Vences et al., 2003; Yoder and Nowak, 2006; Vences Gavrilets and Losos, 2009; Valente et al., 2020). Limited area and long- et al., 2009; Crottini et al., 2012). The combination of these circum­ standing isolation are key features that make islands ideal settings to stances makes Madagascar an ideal setting to study in situ diversification study evolutionary radiations and diversification (Losos and Ricklefs, (Ali and Aitchison, 2008; Vences et al., 2009; Crottini et al., 2012; 2009; Warren et al., 2015). The role of in situ diversification for island Samonds et al., 2012), and several Malagasy endemic lineages are biodiversity was initially considered negligible for most islands systems salient examples of impressively species-rich evolutionary radiations (e.g., MacArthur and Wilson, 1963), but has gained increasing recog­ (Ganzhorn et al., 2014; Goodman and Benstead, 2003). However, only a nition in recent times (see e.g., Presgraves and Glor, 2010; Warren et al., limited number of Malagasy vertebrate radiations have been investi­ 2015). In situ diversification is positively correlated with island size, and gated to unveil their processes of in situ diversification (e.g., vangid there appears to exist a minimum threshold in surface area above which birds, Jønsson et al., 2012; lemurs, Herrera, 2017; pseudoxyrhophiine local diversification becomes the primary source of diversity, with its snakes, Burbrink et al., 2019). contribution to local diversity increasing linearly with area (Heaney, Here we focus on the diversification of the species-rich Malagasy 2000; Losos and Schluter, 2000; Kisel and Barraclough, 2010; Pre­ scincine lizards (subfamily Scincinae). This radiation endemic to sgraves and Glor, 2010; Warren et al., 2015). Madagascar and neighbouring islands (see 2.1 Samples and dataset) has The geographic, geological, and environmental features that provide recently been the object of intensive systematic and taxonomic research, suitable conditions for in situ diversification differ amongst oceanic and with several species described over the last three decades (e.g., Rax­ continental islands. Yet, the study of island evolutionary radiations has worthy and Nussbaum, 1993; Andreone and Greer, 2002; Köhler et al., been predominantly focused on oceanic archipelagos. The Galapagos 2009; Miralles et al., 2011a, 2011b, 2011c, 2012, 2016a, 2016b), and and Hawaii provide amongst the best-known examples (e.g., Darwin and the knowledge of their phylogenetic relationships that has greatly Murray, 1859; Lack, 1947; Warren et al., 2015; Meiri, 2017). The fact improved (e.g., Crottini et al., 2009; Erens et al., 2017; Köhler et al., that oceanic islands emerge from the sea surface as blank slates can 2010; Miralles et al., 2015; Miralles and Vences, 2013; Schmitz et al., facilitate lineage diversification (i.e., providing new land and reduced 2005; Whiting et al., 2004). Their richness and variety emphasise competition for early colonisers; Toussaint et al., 2015; Ali, 2017). Madagascar’s status as a global hotspot for skink diversity (Chapple However, small sizes and young ages are features that can limit the et al., 2021). Nonetheless, part of the diversity of this group likely re­ extent of in situ diversification (McPeek and Brown, 2007; Kisel and mains scientifically undescribed, as demonstrated by the large number Barraclough, 2010; Toussaint et al., 2015). In contrast, the rich biocli­ of candidate species identified in recent years (Miralles et al., 2016b; matic and environmental diversity that can be available on large con­ Miralles and Vences, 2013; Nagy et al., 2012). Molecular phylogenies tinental islands can provide ideal conditions for a multitude of suggest that the most recent ancestor of all Malagasy Scincinae reached diversification processes (Vences et al., 2009; Ali, 2017; Meiri, 2017; Ali Madagascar from Africa 65–45 Ma (Whiting et al., 2004; Crottini et al., and Meiri, 2019), and, if the islands are old enough, time may favour 2012; Miralles et al., 2015) and that in situ diversification is responsible prolonged diversification. On the other hand, the presence of rich and for the large variety of species and forms observed today. Malagasy structured biotas on continental islands may in turn limit the opportu­ Scincinae have colonised most geographic and bioclimatic regions of the nities for the diversification of new colonisers (McPeek and Brown, island, although species are differently distributed across Madagascar’s 2007; Toussaint et al., 2015; Ali, 2017). regions (Brown et al., 2016). The eastern rainforest hosts the majority of Madagascar is the fourth largest island in the world (ca. 587,000 the species of the group and the humid forest in northern Madagascar is km2), initially part of the Gondwana supercontinent until it separated the most relevant hotspot for scincines diversity (Brown et al., 2016). from Africa ca. 160–130 million years ago (Ma), it remained isolated This distributional pattern, shared with other amphibian and reptile from all other landmasses since ca. 80–65 Ma (Ali and Aitchison, 2008; groups (Brown et al., 2016), would suggest that Madagascar’s large Ali and Hedges, 2022). Madagascar’s Cretaceous vertebrate biota was environmental diversity and bioclimatic zonation played a role in the likely almost entirely depleted during the mass extinction at the diversification of Malagasy Scincinae. boundary between the Cretaceous and Paleogene (K-T boundary) Malagasy Scincinae display an exceptional array of ecomorpho­ (Samonds et al., 2013) and time-calibrated phylogenies and biogeo­ logical adaptations to the environments they inhabit (Fig. 1; Whiting graphic reconstructions suggest that most present-day biodiversity re­ et al., 2004; Glaw and Vences, 2007; Miralles et al., 2015, 2016a). Fully sults from subsequent colonisation events rather than vicariance quadrupedal species (genera Amphiglossus, Brachyseps, Flexiseps, and (Vences et al., 2003; Yoder and Nowak, 2006; Vences et al., 2009; Madascincus) mostly inhabit forested regions, with at least one genus Crottini et al., 2012; Samonds et al., 2012; Welt and Raxworthy, 2022). (Amphiglossus) specialised to live in close association with freshwater. Biogeographic studies revealed a peak in successful colonisations Limbless or limb-reduced species often exhibit several additional ad­ through multiple overseas dispersal events around the K-T boundary aptations to a fossorial (burrowing) lifestyle, such as body elongation, that continued through most of the Cenozoic (Crottini et al., 2012; modification of cephalic scalation, depigmentation, and regressed eyes Samonds et al., 2012, 2013). The subsequent in situ diversification of and ear-openings (Glaw and Vences, 2007; Miralles et al., 2011b, 2012, these colonisers is responsible for most of present-day Malagasy verte­ 2015, 2016a). Fossorial species mostly inhabit either the forest litter brate biota (Vences et al., 2003; Yoder and Nowak, 2006; Vences et al., (here termed ‘litter-burrowers’; genera Pseudoacontias and Paracontias 2009; Crottini et al., 2012; Samonds et al., 2012, 2013), which is except for Paracontias minimus that is found in sandy soils) or sandy composed of several endemic lineages, some of which represent high environments (here termed ‘sand-swimmers’; genera Grandidierina, taxonomic ranks (e.g., genera, subfamilies, families). Pygomeles, and Voeltzkowia). The most recent investigation of the evo­ The old age of Madagascar as an island (ca. 80–65 Ma), its rich lution of the body plan in Malagasy Scincinae suggested that the fully environmental heterogeneity, and large size have favoured intense and quadrupedal form (non-fossorial lifestyle) represents the plesiomorphic prolonged diversification (Vences et al., 2009; Kisel and Barraclough, condition from which different morphological adaptations to a fossorial 2010; Crottini et al., 2012). In addition, the extinction of the Malagasy lifestyle have convergently evolved multiple times (Miralles et al., 2 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 2015). In addition, the almost concurrent divergence of fossorial sand- (2a) Can an early burst in diversification rate, often observed in swimmer genera in different regions of the island during the Miocene evolutionary radiations (Schluter, 2000; Gavrilets and Losos, 2009; was suggested as evidence for the establishment of dry and sandy en­ Simões et al., 2016) as, for instance, in Madagascar’s vangid birds vironments in Madagascar during that period, although no biogeo­ (Jønsson et al., 2012), explain the pattern of diversity observed in this graphic reconstruction was performed to support this hypothesis group? Alternatively, did diversification rate in Malagasy scincines in­ (Miralles et al., 2015). crease towards the present (e.g., as previously identified in lemurs; Despite recent progress in the study of Malagasy Scincinae, we are Herrera, 2017), or followed a constant trend and a subsequent slow­ still lacking a comprehensive characterisation of the in situ diversifica­ down in the last million years as in pseudoxyrhophiine snakes from tion of the group. Following the recent systematic and taxonomic ad­ Madagascar (Burbrink et al., 2019)? vances outlined above and the availability of a large number of new (2b) Is there an association between small body size and increased samples, it is now possible to infer a new phylogenetic hypothesis of the rate of diversification in Malagasy Scincinae, as reported for other taxa group based on a more complete dataset (both in terms of number of in Madagascar (e.g., Wollenberg et al., 2011)? included taxa and loci). The availability of this dataset makes the esti­ (2c) Is there an association between the several parallel transitions mation of divergence times more reliable, which is, in turn, fundamental from quadrupedal to limb-regressed lineages (Miralles et al., 2015) and to performing informative macroevolutionary analyses within a robust increased diversification rates? phylogenetic and statistical framework. Using new macroevolutionary (2d) Did the bioclimatic zonation of Madagascar play a role in the tools that have become available recently (e.g., Beaulieu and O’Meara, diversification of Malagasy Scincinae, and is it possible that past large- 2016; Harvey and Rabosky, 2018; Höhna et al., 2022; May et al., 2016; scale climatic changes favoured diversification, as in the case of the Rabosky, 2014; Rabosky and Huang, 2016), we can now explore the divergence of the fossorial sand-swimmer genera (Miralles et al., 2015)? roles and contributions of intrinsic and extrinsic determinants of species diversification in the group by testing a variety of hypotheses on their in 2. Material and methods situ diversification. This contributes to the understanding of the evolu­ tion of Malagasy fauna and provides an accessible and replicable 2.1. Samples and dataset methodological framework that can be applied to study other radiations to identify possible shared patterns and improve our knowledge on the in Malagasy Scincinae comprise 65 described species in nine genera: situ diversification processes within Madagascar and other island sys­ Amphiglossus, Brachyseps, Flexiseps, Grandidierina, Madascincus, Para­ tems. To achieve a comprehensive characterisation of the diversification contias, Pseudoacontias, Pygomeles, and Voeltzkowia (Fig. 1; Supplemen­ of the Malagasy Scincinae we focussed on: tary Table S1; Glaw and Vences, 2007; Uetz et al., 2022). The genus (1) Obtaining a robust phylogenetic framework. The phylogenetic Mesomycterus Cope 1892 has recently been proposed to be an older relationships amongst major clades remain poorly resolved, and several synonym of Brachyseps Erens et al., 2017 (Shea, 2021), but since it has undescribed taxa have not yet been included in molecular phylogenetic not been used in the literature since its description, it qualifies as nomen analyses (Erens et al., 2017). A new phylogenetic hypothesis encom­ oblitum (Uetz et al., 2022), and it is not considered in this study. Most passing most described and known candidate species and an estimation species belonging to the Malagasy Scincinae are restricted to of their divergence times are needed to investigate the diversification of Madagascar, although two species of Flexiseps are endemic to the the group; neighbouring islands of Mayotte and Glorioso (Glaw and Vences, 2007; (2) Hypothesis testing: Uetz et al., 2022). Fig. 1. An overview of the Malagasy Scincinae diversity. (A) Pseudoacontias menamainty, Berara. (B) Brachyseps macrocercus, Ankaratra. (C) Madascincus CaNew Betampona (of the Madascincus “nanus” group), Betampona. (D) Madascincus igneocaudatus, Ifaty. (E) Flexiseps ornaticeps, Ifaty. (F) Amphiglossus astrolabi, Sainte Luce. (G) Voeltzkowia yamagishii, Ankarafantsika. (H) Voeltzkowia mira, Matsedroy. (I) Paracontias vermisaurus, Makira. (J, K) Pygomeles braconnieri, near Ifaty. (L) Gran­ didierina sp. aff. rubrocaudata Andranomaitso. (Credits: A - Franco Andreone; B - Frank Glaw and Miguel Vences; C, D, E, L - Gonçalo M. Rosa; F, H - Sam Hyde Roberts; G - Falk Eckhard; I - Miguel Vences; J, K - Aurélien Miralles). 3 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 We compiled a dataset comprising all individuals of Malagasy Scin­ posterior distributions with Tracer 1.7.1 (Rambaut et al., 2018) and cinae whose tissue samples were available to us (Supplementary considered a value of 200 Effective Sampling Size (ESS) as a minimum Table S2). These samples were sequenced for at least one of the three threshold to evaluate the fitness of each prior. We used FigTree 1.4.4 mitochondrial markers that are generally used for species identification (Rambaut, 2009) to visualise the 50%-majority rule consensus tree. of Malagasy Scincinae, which allowed us to identify and assign them Following the procedure described above, we inferred distinct to­ unequivocally to a mitochondrial interspecific lineage: 16S rRNA 3′ pologies for the concatenated alignments of mitochondrial and nuclear terminus (16S), cytochrome oxidase I gene (COI), and NADH dehydro­ markers, respectively. We ran two independent runs of 70-million gen­ genase subunit 1 (ND1) (with adjacent tRNAs) genes (see Supplemen­ erations each with MrBayes for both phylogenies and assessed topo­ tary Text S1 for details on laboratory work). To do this, we applied the logical congruence between the mitochondrial, nuclear and full skinks’ standard interspecific threshold of 6.1–6.5% genetic distance at concatenated (mitochondrial and nuclear) topologies with the Icong the COI marker, as proposed in Nagy et al. (2012), and used 16S and index (de Vienne et al., 2007; R script available on GitHub at https://gi ND1 to refine species identification when necessary. thub.com/damiendevienne/icong). We compiled a phylogenetic multi-locus dataset using the matrix of To explore the influence of alternative inference methods on tree Erens et al. (2017) as a starting point. This matrix contained 45 inter­ topology, we inferred a maximum likelihood (ML) phylogeny with specific lineages (either nominal species or candidate species as defined RAxML-HPC2 8.2.12 (Stamatakis, 2014) on the full multi-locus matrix by Vieites et al., 2009) and 13 loci. We augmented this matrix by adding (Supplementary Table S3). We used the same partition scheme esti­ one locus (COI) and sequences from 26 additional lineages (seven of mated with PartitionFinder2 for Bayesian analyses (see above). We ran which were identified for the first time; Supplementary Tables S2–S3) RAxML on CIPRES setting 1,000 replicates of rapid bootstrap and and removing three taxa corresponding to intraspecific lineages. Lineage GTRCAT model. In addition, we inferred a species tree with ASTRAL identification in our phylogenetic multi-locus dataset was performed 5.7.8 (Zhang et al., 2018) starting from gene trees of the 14 markers following the procedure described above. Our final molecular dataset available in our dataset. We inferred the gene trees with MrBayes included 68 ingroup taxa (47 of the 65 nominal species of Malagasy following the same procedure described above for concatenated data­ Scincinae and 21 candidate species; Supplementary Table S1) and 14 sets, with the only exception that we used jModeltest 2.1.10 (Darriba loci, which comprised four mitochondrial (12S rRNA, 16S, COI, and ND1 et al., 2012) on CIPRES to find the best model of sequence evolution for with adjacent tRNAs) and 10 nuclear markers (brain-derived neuro­ each gene. We performed two separate runs of 50-million generations. trophic factor (BDNF), oocyte maturation factor (CMOS), phosducin To estimate the species tree we ran ASTRAL with default parameters (PDC), recombination activating gene 1 and gene 2 (RAG1 and RAG2), using the 50%-majority rule consensus trees for each gene. alpha-enolase (ENO1), two non-overlapping segments of sacsin (SACS-A and SACS-B), leucine-rich repeat and WD repeat-containing protein 2.3. Dated phylogeny (KIAA1239), and titin (TTN)) (Supplementary Table S3; Supplementary Text S1). To maximise the number of available loci, we used sequences Given the absence of fossils in the ingroup, we used fossil-based of multiple congeneric species (chimera sequences) of Cordylus, Eume­ calibrations on the outgroups. We compiled a reduced dataset ces, Tiliqua, and Scincella as outgroups, as in Crottini et al. (2009) and including four Malagasy Scincinae and seven outgroups (Supplementary Erens et al. (2017) (Supplementary Table S3). Newly generated se­ Table S5). We selected one species for each of the outgroups (Cordylus, quences have been deposited in GenBank (OP256430–OP256524, Eumeces, Tiliqua, and Scincella) of the complete species-level dataset, and OP263988–OP264059, OP271309–OP271456; Supplementary the remaining outgroups were selected based on the fossil calibrations Table S3). For details on laboratory work and sample amplification see that were used (Supplementary Table S5); phylogenetic placement and Supplementary Text S1. age of fossil calibrations were expert-verified (by A. Villa) (see Supple­ mentary Table S6 for information on the selected fossils and justification 2.2. Phylogenetic analyses of their use to calibrate the tree). We aligned the dataset with MAFFT 7.310 and inferred the best-fit partitioning scheme and models of evo­ We performed phylogenetic analyses on a concatenated dataset lution with PartitionFinder2 2.1.1, under BIC, on the CIPRES Science composed of all mitochondrial and nuclear markers available in our Gateway. We estimated divergence times with BEAST 2.6.0 (Bouckaert dataset (Supplementary Table S3). We used PipeLogeny in R 4.0.2 et al., 2019) (see Supplementary Text S1 for all analysis details). The (Muñoz-Pajares et al., 2019; R Core Team, 2021), which is a recently height average and 95% Higher Posterior Density Interval (HPDI) of the developed pipeline that automatises the preparation of the input files for crown node of the Malagasy Scincinae was used to calibrate the corre­ several commonly used phylogenetic tools, which are then run sepa­ sponding node on the complete species-level dataset, and we estimated rately in their standalone software applications. We aligned sequences divergence times on this dataset with BEAST (Supplementary Table S3; with MAFFT 7.310 (Katoh and Standley, 2013), detected misaligned and see Supplementary Text S1 for all analysis details). hypervariable regions of the alignment to be excluded from the phylo­ genetic analyses with Gblocks 0.91b (Castresana, 2000), identified the 2.4. Diversification rate analyses best-fit partitioning scheme and models of evolution using Partition­ Finder2 2.1.1 (Lanfear et al., 2017), and performed Bayesian analyses We inferred diversification rates from the dated phylogeny pruned using MrBayes 3.2.7 (Ronquist et al., 2012). Due to the limited number from all outgroups. Whenever possible we set a sampling fraction of of available gene fragments and the unavailability of recently collected 76%, corresponding to the percentage of taxa included in the phylogeny tissue samples, we excluded A. splendidus and P. menamainty from the relative to the overall number of nominal species and candidate species phylogenetic analyses (see Supplementary Table S3 for the phylogenetic recognised for the group. Since several models and methods that we dataset). After the preparation of the input files, we ran Gblocks on the employed in the following analyses assume rate homogeneity across GBlocks Server (Castresana, 2000) using the least stringent parameter lineages, we first checked for possible heterogeneity of diversification settings. We used PartitionFinder2 on the CIPRES Science Gateway rate within the tree with BAMM 2.5.0 (Rabosky, 2014). BAMM identifies (Miller et al., 2010) to find the best-fit partitioning scheme and models heterogeneity in the diversification process through time and across of sequence evolution, which were tested under Bayesian Information branches within a phylogeny using the reversible-jump MCMC Criterion (BIC). We ran MrBayes on CIPRES with two runs of 100-million (rjMCMC). The priors of the evolutionary rate parameters were set based generations each, setting four independent chains for each run with on the function ‘setBAMMpriors’ of the package BAMMtools 2.1.7 sampling at every 1,000th generation and 40% burn-in of the Markov (Rabosky et al., 2014b) available on R 4.0.2. We tested different priors Chain Monte Carlo (MCMC) chain. We assessed run convergence and for the values of the expected number of diversification rate shifts (0.1, 4 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 0.5, 1, 5, and 10) and for each value we ran two MCMC chains of 10- relationships between diversification rates variation through time and million steps to assess reciprocal convergence, while we ensured lineage diversity accumulation through the nine models from Morlon single-chain convergence through the ESS (minimum threshold of 200). et al. (2010) (Supplementary Table S9), whose maximum likelihood We set a non-random incomplete taxon sampling, which accounts for the expressions are based on the coalescent process. We used the ‘fit_­ different numbers of included taxa amongst different genera of the coal_cst’ and ‘fit_coal_var’ functions in the R package RPANDA 1.8. The group. The output was inspected with the R package BAMMtools 2.1.7. tested models imply both saturated diversity scenarios with turnover We evaluated the statistical support of the number of diversification rate rate (the rate at which an extinct lineage is substituted by an emerging shifts based on Posterior Probability (PP) and Bayes Factor after dis­ one) and scenarios in which diversification is not diversity-dependent carding 15% of the posterior samples as burn-in. Moore et al. (2016) and diversity is not saturated (Supplementary Table S9). Model fitting criticised BAMM highlighting possible issues concerning the likelihood was evaluated through the AICc weights (Burnham and Anderson, function and the prior models employed in the software. To account for 2002). these uncertainties in diversification rate analyses and the several other Recently Louca and Pennell (2020) demonstrated that under a analyses in which we used this software (see 2.5 Trait evolution), we deterministic diversification process, infinite alternative homogeneous decided to run BAMM along with other tools to highlight potential birth–death models could produce the same deterministic lineages- discordance in the results that might be indicative of technical issues. through-time curve and be equally probable of generating the same We further investigated the possible heterogeneity in diversification tree. This set of alternative congruent models is defined as a congruence rates by examining the spectral density profile of the tree computed from class, characterised by the same identifiable variables of pulled diver­ the Laplacian graph of its branch lengths (Lewitus and Morlon, 2016). sification rate, product ρλ0 (between the fraction of sampled extant We obtained the eigenvalues and spectral density profile of the tree with species and speciation rate at present) and pulled speciation rate, which the function ‘spectR’ implemented in the R package RPANDA 1.8 is defined as the relative slope of the deterministic lineages-through- (Morlon et al., 2016). To test the number of modalities that characterise time curve. This approach differs from the hypothesis-driven model- the tree, that can represent the number of heterogeneities in diversifi­ fitting methods applied in some of the previous analyses and tradition­ cation rates, we used the function ‘BICompare’. We compared BIC values ally used in macroevolutionary studies (Morlon et al., 2020). Starting associated with the empirical tree and randomly bifurcating trees par­ from diversification rate estimates from CoMET analyses, we used the R ametrised from this to provide a significance test for the number of package ACDC (Höhna et al., 2022; Louca and Pennell, 2020) to modalities that were previously highlighted. construct the congruence class of all extinction and speciation rate We investigated the occurrence and timing of mass extinction events functions with the same likelihood score and pulled diversification and and tree-wide diversification rate changes through time with the CoMET speciation rates to explore the diversity of rate functions belonging to model within the R package TESS 2.1.0 (Höhna, 2013; May et al., 2016). the class and uncover the array of potential diversification models un­ This method implements a piecewise-constant speciation and extinction derlying Malagasy Scincinae diversification. We tested a series of rate model with the number of rate shifts through time that is inferred speciation and extinction rate functions (separately from each other) of through rjMCMC. We set the function ‘tess.analysis’ to account for different models and temporal trends towards the present (increasing empirical estimates of the hyperpriors of diversification-rates parame­ and decreasing exponential, linear and episodically changing birth­ ters. We ran the rjMCMC for 50-million steps of which the initial 14.5% –death models). We computed speciation, extinction, relative extinction was discarded as burn-in. We considered a minimum ESS of 200 to and net-diversification rates for each category of tested rate function and ensure convergence of the chains. The occurrence of mass extinction compared trends through time with the reference model (CoMET events and shifts in speciation and extinction rates was evaluated diversification rates estimates) to discern between functions belonging through Bayes Factor. to the same congruence class and those that did not. To further evaluate the TESS results, we investigated the possible time-dependence of diversification rates through a maximum likelihood 2.5. Trait evolution model-fitting approach. We tested six models reported in Morlon et al. (2011) accounting for both time-independent and time-dependent 2.5.1. Body size evolution variation of speciation and extinction rates as well as periods of We investigated body size evolution using snout-vent length (SVL) as expanding and declining diversity (Supplementary Table S7). These a proxy. We compiled a list of individual SVL measurements of speci­ models were fitted with the function ‘fit_bd’ implemented in the R mens from museum collections and complemented it with published package RPANDA 1.8 (Supplementary Table S7). The best-fitting model data to obtain a matrix of a total of 60 taxa and 341 SVL measures was chosen based on the corrected Akaike Information Criterion (AICc) (Supplementary Table S10). SVL measures were not available for Pseu­ weights (Burnham and Anderson, 2002; Morlon et al., 2016). doacontias CaNew. Since Pseudoacontias spp. are the largest Malagasy After detecting evidence for time-dependence of diversification rates Scincinae (Glaw and Vences, 2007), and their inclusion or exclusion may (see 3.3 Diversification rate), we investigated the possibility that time- have an influence on the analyses, we ran a replicate of these analyses variation was influenced by an accumulation of lineage diversity including the SVL of P. menamainty as a proxy value for the genus. We through time. We used the R package DDD 4.4 (Etienne et al., 2012, averaged SVL measures over conspecific individuals and log- 2016), which implements a framework for adaptive radiations, model­ transformed (Supplementary Table S11). We applied a global sam­ ling different aspects of diversity-dependence dynamics of diversifica­ pling fraction of 69% (with Pseudoacontias exclusion) and 70% (with tion rate and accounting for possible time-dependent decreases in Pseudoacontias inclusion) accounting for missing taxa to the following diversification rate that can follow an increase in species diversity and analyses, whenever possible. the presence of a diversity carrying capacity. Following Scantlebury We performed model-testing across different models of continuous (2013), we estimated the parameters (speciation rate, extinction rate, trait evolution (Brownian motion, Ornstein-Uhlenbeck process, and and carrying capacity) of different diversity-independent and diversity- Early Burst) with the function ‘fitContinuous’ in the R package geiger dependent models with the function ‘dd_ML’ and evaluated the best- 2.0.7 (Pennell et al., 2014) (Supplementary Table S12). We performed fitting model based on the AICc weights (Burnham and Anderson, the analyses on 100 trees randomly sampled from the posterior distri­ 2002) (Supplementary Table S8). The estimated parameters of the best- bution of the dated phylogeny. We considered the number of trees in fitting model were used in the function ‘dd_LR’, which implements a which a model resulted as the best-fitting (based on AICc weights; bootstrap likelihood ratio test of a diversity-dependent model against a Burnham and Anderson, 2002) to evaluate the statistical support of each constant-rates birth–death model (Etienne et al., 2016). model over the entire posterior sample and over the phylogenetic We tested the robustness of DDD results by investigating the possible inference uncertainty. In addition, we averaged AICc weights over the 5 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 100 trees for each tested model. We performed ancestral state recon­ prone to type I error compared with BiSSE and other SSE methods struction based on the best-fitting model on the 50%-majority rule (Beaulieu and O’Meara, 2016; Caetano et al., 2018). We tested 29 consensus tree (see 2.2 Phylogenetic analyses) using the function ‘con­ models comprising BiSSE-like models (with no hidden trait; Maddison tMap’ in the R package phytools 0.7-70 (Revell, 2012). et al., 2007) and several additional models with the presence of one or We investigated possible rate shifts in SVL evolution with BAMM. We two possible states of a hidden trait and considered both character- used the function ‘setBAMMpriors’ from BAMMtools to find a suitable dependent and character-independent evolutionary models (Supple­ set of priors for our dataset and tested different values of the expected mentary Table S14). We performed the analyses on the same 100 pos­ number of rate shifts priors for the evolution of this trait (0.1, 0.5, 1, 5, terior trees used for testing body size evolution (see 2.5.1 Body size and 10). We ran two MCMC chains of 10-million iterations. The output evolution) and each model statistical support over the 100 trees was was processed with BAMMtools. We assessed multiple and single-chain evaluated as described therein. We fitted the most statistically supported convergence as described for BAMM diversification rate analyses (see model across the 100 trees to the 50%-majority rule consensus tree for 2.4 Diversification rate analyses). We selected the most supported ancestral state reconstruction. We tested the possible association be­ number of trait evolution rate shifts based on the Bayes Factor and we tween trait evolution and diversification rate with STRAPP in BAMM­ chose the best configuration shift based on the highest Posterior tools. The analysis was run on the consensus tree with identical settings Probability. as in previous analyses of body size evolution (see 2.5.1 Body size We tested the possible correlation between SVL and diversification evolution). Unlike those, we used the Mann-Whitney U statistic to test rates with two tip-rate correlation methods: ES-sim (Harvey and the association with speciation and extinction rates. Rabosky, 2018) and ‘STructured RAte Permutations on Phylogenies’ (STRAPP) (Rabosky and Huang, 2016). These methods are more robust 2.6. Historical biogeography compared to other approaches (e.g., QuaSSE; Fitzjohn, 2010) and in particular with small phylogenies (Harvey and Rabosky, 2018; Rabosky We compiled a dataset of 391 geographic records relating to all and Huang, 2016). ES-sim tests the correlation between trait variation species included in our phylogenetic analyses (Supplementary across the tips of a phylogeny and diversification rates. It applies the Table S15). We compiled these species distribution records from own Inverse Equal Splits as summary statistic of branching pattern and tests field records, GenBank, taxonomic literature or from Glaw and Vences for deviations of the empirical phylogeny from a null distribution (2007). We revised all records following current Malagasy Scincinae simulated with no evolutionary model (Harvey and Rabosky, 2018). We taxonomy (Supplementary Table S15). We performed historical bioge­ ran ES-sim with the R code provided in Harvey (2017) on the same 100 ography analyses using the 5-region bioclimatic classification (rain­ trees sampled from the posterior distribution used in previous SVL an­ forest, dry deciduous forest, sub-humid forest, sub-arid spiny thicket, alyses (see above) and set 1,000 simulations for the null distribution. We and montane thickets) modified from Yoder and Nowak (2006) by assessed the number of trees over the total 100 for which the test was merging the Sambirano ecoregion with ‘rainforest’. We categorised statistically significant (p-value < 0.05) or not (p-value > 0.05). species into the regions with the R package speciesgeocodeR 2.0–10 We ran STRAPP using the function ‘traitDependentBAMM’ in (Töpel et al., 2016) and performed ancestral range reconstruction on the BAMMtools. In this analysis, a test statistic measuring the association dated phylogeny pruned from outgroups with the R package Bio­ between the trait and diversification rates estimated on a phylogeny is GeoBEARS 1.1.2 (Matzke, 2018) using the R codes provided on the compared between the empirical values and a null distribution gener­ PhyloWiki webpage (https://phylo.wikidot.com/biogeobears) and from ated by structured permutations of rates across the phylogeny (Rabosky Antonelli et al. (2018). We tested three biogeographic models that are and Huang, 2016). We set 100,000 permutations on the posterior sam­ implemented in the package within a maximum likelihood framework to ples drawn from the output of previous BAMM analyses on diversifica­ detect biases amongst tested methods: Dispersal-Extinction- tion rate (see 2.4 Diversification rate analyses) for the expected number Cladogenesis (DEC) (Ree and Smith, 2008), a likelihood interpretation of diversification rate shifts 0.1, 0.5, 1, 5, and 10. We used Pearson’s and of parsimony DIVA model (‘DIVALIKE’; Ronquist, 1997), and a Spearman’s correlation coefficients as test statistics to test the correla­ maximum likelihood interpretation of the Bayesian-framework model tion between trait evolution and speciation and extinction rates. implemented in the software BayArea (‘BAYAREALIKE’; Landis et al., 2013). While all three models use two free parameters of anagenetic 2.5.2. Fossorial evolution change of range evolution (dispersal and extirpation rates), they differ in We coded Malagasy Scincinae lifestyle as a binary character (fosso­ how they model cladogenetic inheritance of the ancestral ranges rial or non-fossorial) based on the phenotype of each species. We (Matzke, 2018). We did not include the null range and set the maximum included in the fossorial category limb-reduced species (which are range size to match the largest number of areas occupied by extant characterised by a comprehensive body plan adaptation to this lifestyle species. including other apomorphic traits, e.g., body elongation, regressed eyes; Glaw and Vences, 2007; Miralles et al., 2011b, 2012, 2015, 2016a), 3. Results while fully quadrupedal species were coded as ‘non-fossorial’: 0-fosso­ rial (genera Grandidierina, Paracontias, Pseudoacontias, Pygomeles, and 3.1. Phylogenetic hypothesis Voeltzkowia); 1-non-fossorial (genera Amphiglossus, Brachyseps, Flexiseps, and Madascincus) (Supplementary Table S13). Whenever possible, we We generated 168 sequences (Supplementary Table S3). The final accounted for incomplete taxon sampling (77%, computed separately multi-locus concatenated alignment contained 10,416 bp (2,456 bp for both states). from mitochondrial markers and 7,960 bp from nuclear markers), We explored fossorial evolution with the ‘Hidden State Speciation reduced to 9,573 after excluding the 843 positions identified as mis­ and Extinction’ (HiSSE) method implemented in the R package hisse aligned or hypervariable by Gblocks. The best-fitting partitioning 1.9.8 (Beaulieu and O’Meara, 2016), which estimates the effect of bi­ scheme included 12 partitions (see Supplementary Table S16, which nary characters on diversification rate within a maximum-likelihood also provides the best model of sequence evolution for each partition). framework for model-testing. Unlike the ‘Binary-State Speciation and Relationships amongst the different genera and all major clades Extinction’ (BiSSE) and other SSE methods (Maddison et al., 2007; retrieved full support (PP = 1), including the clades corresponding to the Fitzjohn, 2010), HiSSE accounts for additional unmeasured (hidden) recently proposed genera Amphiglossus, Brachyseps, and Flexiseps (Fig. 2; states, which may potentially influence diversification and transitions see also Erens et al., 2017). Pseudoacontias sp. (clade C) was found to be rates amongst states, and implements more realistically complex sister to all other Malagasy Scincinae, which are organised in two character-independent null models. These aspects make HiSSE less additional main clades: clade A containing the genera Madascincus and 6 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 Fig. 2. 50%-majority rule consensus tree inferred with MrBayes on the species-level dataset (Supplementary Table S3). Posterior Probability values are reported before each node as follows: ‘**’ = 1, ‘*’ = 0.99, ‘(*)’ = 0.95–0.98. Values below 0.95 not shown. Taxa analysed for the first time in an evolutionary context are highlighted in bold. Letters ‘A’, ‘B’, and ‘C’ indicate the three main clades of the group. The colours indicate genera. Outgroups are not shown. Representatives of each genus are displayed in photographs (in descending order): Brachyseps frontoparietalis, Amphiglossus reticulatus, Pygomeles trivittatus, Voeltzkowia yamagishii, Flexiseps tanysoma, Grandidierina fierinensis, Paracontias ampijoroensis, Madascincus sp. aff. mouroundavae Marojejy, and Pseudoacontias menamainty. Paracontias, and clade B containing the genera Amphiglossus, Brachyseps, (PP = 1). The few exceptions were: the position of Flexiseps meva as sister Flexiseps, Grandidierina, Pygomeles, and Voeltzkowia. The genus Grandi­ to the clade containing F. crenni, F. tanysoma, F. mandokava, F. johannae, dierina was found to be sister to the clade containing Amphiglossus, and F. ardouini (PP = 0.96); the monophyly of the clade containing all Brachyseps, Flexiseps, Pygomeles, and Voeltzkowia. In this clade, the genus Madascincus with the exception of Madascincus sp. CaNew Betampona Flexiseps was sister to Amphiglossus, Brachyseps, Pygomeles, and Voeltz­ and Madascincus sp. “baeus” (PP = 0.99); the sister relationship between kowia. In this other clade, Voeltzkowia was sister to the clade composed the clade composed of Madascincus sp. aff. minutus Marojejy and of Amphiglossus and Pygomeles and all together were sister of Brachyseps. Madascincus minutus Montagne d’Ambre with the clade of Madascincus Within the Madascincus clade, Madascincus sp. CaNew Betampona was minutus Tsaratanana and Madascincus sp. aff. minutus An’Ala (PP = sister to Madascincus sp. “baeus”, and together they were sister to all 0.99); the sister relationship between Madascincus minutus Tsaratanana other Madascincus. and Madascincus sp. aff. minutus An’Ala (PP = 0.79); the sister rela­ Nearly all nodes at the intra-generic positions were fully supported tionship of Madascincus polleni with the clade containing Madascincus sp. 7 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 aff. stumpffi Forêt d’Ambre, Madascincus sp. aff. stumpffi Marojejy, and mitochondrial and nuclear markers showed some topological differ­ Madascincus stumpffi (PP = 0.79). Interestingly, Flexiseps crenni was ences between each other and the full multi-locus concatenated dataset found to be sister to the clade made of F. ardouini, F. johannae, (mitochondrial with nuclear) (see Supplementary Figs S1–S2). However, F. mandokava, and F. tanysoma (confirming the tentative position pro­ the Icong index indicated a statistically significant overall topological posed by Erens et al., 2017). Paracontias holomelas had been excluded by congruence between all phylogenies (p-value < 0.01; 2.86 mitochon­ Erens et al. (2017) and Miralles et al. (2016b) because of problems with drial vs nuclear; 3.35 mitochondrial vs full; 3.63 nuclear vs full), sug­ the assignment of deposited sequences. After the revision of published gesting that both phylogenetic signals of mitochondrial and nuclear sequences and the generation of new sequences from recently collected markers contributed to the topology of the full dataset. The ML phy­ material for this taxon, we included P. holomelas in our phylogenetic logeny inferred with RAxML on the full concatenated dataset (Supple­ analyses and found this taxon to be sister of P. hildebrandti. mentary Fig. S3) showed differences from the MrBayes tree in some The phylogenies inferred from the concatenated datasets of inter-generic relationships (i.e., the position of the genus Flexiseps and Fig. 3. Time-tree and diversification rates through time of Malagasy Scincinae. (A) Divergence times estimate performed on the species-level dataset (Supplementary Table S3). The red circle indicates the node that was calibrated based on divergence times estimated with primary fossils information on the reduced dataset (Supplementary Fig. S5; Supplementary Table S5). The analysis was performed in BEAST. Outgroups are not shown. Geological times were abbreviated as follows: Cr. (Cretaceous), Pa. (Paleocene), Eo. (Eocene), Ol. (Oligocene), Mi. (Miocene), Pli. (Pliocene), and Ple. (Pleistocene). (B) Net-diversification rate. (C) Bayes Factor of speciation shifts. B–C were inferred with the CoMET model implemented in the R package TESS. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 8 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 Grandidierina in clade B), although most relationships were congruent bifurcating trees (15,052,056), confirming the single modality of the between the two phylogenies. However, differently from the Bayesian spectral density profile, which indicates the presence of homogeneous inference, this tree had overall lower statistical (bootstrap) support (e.g., diversification rates across the tree. most inter-generic relationships in Clade B received no statistical sup­ We investigated diversification rate through time and the possible port). Similarly, the species tree inferred with ASTRAL showed some occurrence of mass extinction events with the CoMET model (see Sup­ topological differences with the MrBayes phylogeny inferred from the plementary Figs S11–S12 for chain convergence). There was no evidence full multi-locus concatenated dataset (Supplementary Fig. S4) and of mass-extinction events (Bayes Factor < 2.02; Supplementary recovered limited statistical support: most inter-generic relationships Fig. S13). Net-diversification rate was found to be constant (0.02–0.12 and the monophyly of some genera (e.g., Flexiseps, Madascincus, Pygo­ lineages/myr) until ca. 6–4 Ma (Late Miocene-Pliocene), when it meles) were not statistically supported. The normalised quartet score decreased to the present-day values of 0.00–0.06 lineages/myr (Fig. 3B). was 0.78 suggesting moderate concordance between gene trees. Bayes Factor values supported the slowdown (Bayes Factor = 3.43–5.72; Considering the overall congruence between the trees inferred from the average Bayes Factor = 5.13; Fig. 3C), corresponding to “substantial mitochondrial, nuclear, and full multi-locus datasets with MrBayes, the evidence” (Jeffreys, 1961). Extinction rate was constant throughout the moderate concordance between gene trees, and the low statistical sup­ evolution of the group and close to 0 (Supplementary Fig. S13), although port of RAxML tree and ASTRAL species tree, we used the highly sta­ this might reflect a general inability of phylogenetic models to infer tistically supported topology inferred from the Bayesian analysis on the extinction (Louca and Pennell, 2021; see below ACDC results). Specia­ full multi-locus dataset (Fig. 2) in the following analyses. tion rate followed the trend that was recovered for net-diversification rate (Supplementary Fig. S13). 3.2. Dated phylogeny Time variation in diversification rate was confirmed by the maximum likelihood model-fitting approaches (Supplementary The concatenated alignment of the reduced dataset contained Table S7). Amongst the six models that were tested from Morlon et al. 11,147 bp, and the best-fitting partitioning scheme included 10 parti­ (2011), model 3 (speciation rate exponentially declining through time tions (Supplementary Table S17). The 95% confidence intervals of and null extinction rate) was the best-fitting (AICc weight = 0.61; divergence time estimates were not elongated and overlapping (Sup­ speciation rate at present = 0.05 lineages/myr; exponential rate varia­ plementary Fig. S5; Fig. 3A). The Scincoidea crown node was dated in tion = 0.03) (Supplementary Table S7). The analyses performed with the the Late Jurassic (160.29 Ma) (Supplementary Fig. S5), while the DDD 4.4 package suggested diversity-dependent diversification in the diversification of the Scincidae family was inferred to be 111.65 Ma group. The model DDLS (linear dependence of speciation rate with di­ (142.83–79.47 Ma 95% HPDI) at the end of the Early Cretaceous. The versity carrying capacity and null extinction rate) fitted the data the best crown node of the Scincinae subfamily was dated at 85.42 Ma (AICc weight = 0.37) with parameters estimates of speciation and car­ (112.79–59.39 Ma 95% HPDI) in the Late Cretaceous, which also rep­ rying capacity of 0.13 lineages/myr and 107, respectively (Supple­ resents the stem node of the endemic Malagasy group relative to the mentary Table S8). The bootstrap likelihood ratio test (Etienne et al., considered outgroups. The crown node of Malagasy Scincinae was 2016) confirmed the prevalence of a diversity-dependence against a estimated at the boundary between Paleocene and Eocene epochs constant-rates model (Likelihood ratio = 3.69; p-value = 0.003; power (55.14 Ma; 77.82–35.82 Ma 95% HPDI). The same node was dated at of the test = 0.93). Amongst the tested models from Morlon et al. (2010), 49.8 Ma (72.67–26.16 Ma 95% HPDI), at the beginning of the Eocene, in model 2 (saturated diversity with the turnover rate declining through the maximum clade credibility tree inferred from the complete dataset time) provided the best-fit (AICc weight = 0.77; turnover rate at present calibrated with secondary information (Fig. 3A). This estimate also = 0.03; exponential rate variation = 0.12; Supplementary Table S9). represents the divergence of Pseudoacontias (clade C) from the rest of the This model was followed by model 6 (non-saturated diversity, with group. The two other main clades diverged slightly after (44.34 Ma; declining speciation rate and null extinction rate) (AICc weight = 0.13), 65.6–23.4 Ma 95% HPDI). The inter-generic diversification in clade A which is conceptually identical to model 3 from Morlon et al. (2011) (see was complete by the mid-Oligocene (ca. 29 Ma), whereas the genera of above) (Supplementary Table S9). clade B diverged more recently, between the mid-Oligocene and With the exception of the episodically changing models and expo­ throughout the Miocene (ca. 28–13 Ma) (Fig. 3A). The genera with the nential models with decreasing trends towards the present (Supple­ most recent crown nodes were Amphiglossus (14.7 Ma; 23.24–6.92 Ma mentary Figs 18–20), the alternative extinction rate functions that we 95% HPDI) and Voeltzkowia (13.79 Ma; 21.61–6.62 Ma 95% HPDI). The tested with ACDC produced speciation and net-diversification rates that most recent cladogenetic event (the only one that was estimated during were overall consistent with the reference model of CoMET diversifi­ the Pleistocene) was dated to 1.8 Ma (3.19–0.72 Ma 95% HPDI) and it cation rate estimates, indicating that these models might be part of the represents the divergence between Flexiseps ardouini and F. johannae. congruence class computed from these (or be close to it) (Supplementary Figs 14–17). Models with increasing extinction rates towards the present 3.3. Diversification rates seemed more congruent, suggesting that such temporal trends might have characterised Malagasy Scincinae diversification. These scenarios BAMM analyses consistently supported a constant diversification shared a slowdown in net-diversification rate that was parallel in time rate across Malagasy Scincinae clades regardless of the choice of the with the reference model (suggesting that the drop estimated with expected number of rate shifts prior (Supplementary Figs S7–S9; see CoMET is robust), although for a few functions within some of the tested Supplementary Fig. S6 for multiple chains convergence). The Bayes category models, such a decrease started slightly before, questioning factors comparing scenarios with number of rate shifts >0 with the null their full congruence with the class (Supplementary Figs 14–17). These model (represented by the 0-rate shift scenario) were never above 1, scenarios also showed speciation rate trends that were consistent with meaning that the null model with no rate shifts could not be rejected the reference, although in some models rate magnitudes were larger and, consequently, supporting a unique evolutionary rate dynamic than in the reference, which is probably determined by higher extinction throughout the phylogeny. These results were confirmed by the analysis rates given similar net-diversification rates. Differently from the of the Laplacian spectrum of the tree in RPANDA (Lewitus and Morlon, extinction rate, none of the speciation rate functions that we tested 2016; Supplementary Fig. S10). The spectral density profile (principal seemed to be part of the congruence class as they often produced eigenvalue = 11,960.49; asymmetry = − 0.18; peakedness = 2.28) incompatible net-diversification and speciation rate trends with the showed a single modality with eigengap between the first and second of reference model and extinction rates that were unrealistically negative the eigenvalues sorted in descending order (Supplementary Fig. S10). (Supplementary Figs 21–25). The BIC value of the tree (7,510,067) was lower than for randomly 9 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 3.4. Trait evolution of the dated phylogeny suggested that the common ancestor to Malagasy Scincinae was a quadrupedal and non-fossorial species, and that limb 3.4.1. Snout-vent length reduction and fossorial lifestyle evolved five independent times in the The inclusion or exclusion of an SVL proxy for the genus Pseudoa­ evolutionary history of the group (Fig. 4B). Character-independent contias did not affect the analyses (Fig. 4A; Supplementary Fig. S26). We evolution was also confirmed by STRAPP analyses, which showed no consequently report results for the analyses including Pseudoacontias association (p-value > 0.05) between the fossorial lifestyle trait and the (Fig. 4A), and the outputs of the analyses excluding the genus Pseu­ evolutionary rates estimated on the phylogeny for any of the tested rate doacontias are provided in Supplementary material (Supplementary Figs shifts priors. S26, S31–S32). The Brownian motion model produced the best-fit over all 100 tested posterior trees (Supplementary Table S12). The ancestral state reconstruction performed on the 50%-majority rule consensus tree 3.5. Historical biogeography of the dated phylogeny reflected model testing results and suggested considerable body size variation within genera (Fig. 4A; Supplementary The results from BioGeoBEARS analyses were overall consistent be­ Fig. S26). In BAMM analyses, a model with a constant rate of trait tween DEC (Fig. 5), DIVALIKE (Supplementary Fig. S33), and BAYAR­ evolution across clades was the most supported (based on both Bayes EALIKE (Supplementary Fig. S34) biogeographic models. Dry deciduous Factor and Posterior Probability values) regardless of the choice of the forest and rainforest biomes were reconstructed as ancestral ranges for number of expected rate shifts prior (Supplementary Figs S28–S30; see the deepest nodes, suggesting a significant role of forested biomes in the Supplementary Fig. S27 for multiple chains convergence). ES-sim ana­ early diversification of Malagasy Scincinae. The current distribution of lyses were not significant (p-value > 0.05) for all 100 analysed posterior most genera is predominantly restricted to individual biomes. For trees, suggesting the absence of correlation between SVL evolution and instance, species in the fossorial genera of Clade B (Voeltzkowia, Gran­ diversification rates. Similarly, the association between trait evolution didierina, and Pygomeles) are restricted to either the sub-arid spiny and speciation and extinction rates investigated with STRAPP was not thicket or dry deciduous forest, and quadrupedal non-fossorial genera of significant (p-value > 0.05) for any of the tested priors for the expected the same Clade B (Flexiseps, Brachyseps, and Amphiglossus) are mostly number of diversification rate shifts. found in the forested humid biomes (rainforest and sub-humid forest) (Fig. 5). The 65 species included in our biogeographic analyses are not 3.4.2. Fossorial evolution equally distributed across biomes. The rainforest hosts the largest There was no evidence of fossorial-dependent diversification from number of species (34), followed by sub-humid forest (24), dry decid­ HiSSE analyses (Supplementary Table S14). BiSSE-like models received uous forest (23) and sub-arid spiny thicket and montane thickets with generally higher statistical support than HiSSE models, suggesting no equal species numbers (10). Dry deciduous forest has the highest num­ evidence for the presence of hidden traits that may have influenced ber of endemic species (i.e., with distributions exclusively restricted to diversification. The model with no difference in speciation and extinc­ individual biomes) (13), followed by rainforest (11), sub-humid forest tion rates between the two states and null transition from fossorial to (6), sub-arid spiny thicket (5), and montane thickets (3). non-fossorial state resulted the best-fitting in 81/100 analysed posterior trees (Supplementary Table S14). The ancestral state reconstruction 4. Discussion performed on this model fitted on the 50%-majority rule consensus tree The present study aimed to explore several facets of the in situ Fig. 4. Ancestral state reconstructions of traits evolution performed on the dated 50%-majority rule consensus tree. (A) Body size (snout-vent length (SVL)) evolution reconstructed with the R package phytools on the dataset including SVL measurement of P. menamainty as a proxy for the genus Pseudoacontias (Supplementary Table S11). The bars on the right represent SVL values of each species. (B) Fossorial lifestyle evolution reconstructed with the R package hisse. States are coded as follows: fossorial (white), non-fossorial (black) (Supplementary Table S13). 10 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 Fig. 5. Ancestral range reconstruction performed with the DEC model implemented in the R package BioGeoBEARS. Species distributions (Supplementary Table S15) were coded with the R package speciesgeocodeR based on Madagascar’s bioclimatic classification modified from Yoder and Nowak (2006) (sub-arid spiny thicket (a, blue), sub-humid forest (b, light blue), dry deciduous forest (c, green), montane thickets (d, yellow), and rainforest (e, red)). Estimated dispersal and extirpation parameters were 0.0284 and 0.1275, respectively. Ancestral ranges with the highest marginal probabilities are letter-coded before each node. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) diversification (diversification rate, phenotypic evolution, and bioge­ following paragraphs we discuss such findings for each tested diversi­ ography) of Malagasy scincine lizards. Based on diversification trends fication hypothesis as well as the new phylogenetic hypothesis and time frequently observed in evolutionary radiations and especially in insular tree. systems (e.g., Schluter, 2000; Wollenberg et al., 2011; Jønsson et al., 2012; Lee et al., 2013; Simões et al., 2016), a temporal trend in diver­ 4.1. Systematics and time-tree sification rate displaying an early-burst pattern and increased diversi­ fication rates for small-sized species and fossorial lineages was expected. In contrast to previous published phylogenetic hypotheses (Whiting Our work also aimed to determine whether the rich bioclimatic diversity et al., 2004; Schmitz et al., 2005; Crottini et al., 2009; Miralles et al., and possible past large-scale climatic changes of Madagascar influenced 2015; Erens et al., 2017), the monophyly of major clades and their the diversification of the group, as expected from the asymmetric dis­ phylogenetic relationships received full support (Fig. 2). This was likely tribution of species richness across biomes (Brown et al., 2016) and the the result of a more complete data matrix, both in terms of a broader almost concurrent divergence of fossorial genera living in dry and sandy range of taxa and additional molecular loci included. However, the short environments (sand-swimmers) (Miralles et al., 2015). Our findings internodes between Clade B genera might indicate that even a larger revealed a contrasting scenario in which we recovered: (1) a constant multi-locus dataset (e.g., obtained by a phylogenomic approach) will be diversification rate through most of the evolutionary history of the needed to reconstruct evolutionary trajectories in the context of a rapid group, (2) a random pattern in the evolution of body size, which was inter-generic divergence. This seems to be confirmed by the difficulty in found not to influence diversification rate, and (3) multiple evolutionary retrieving well-supported phylogenetic hypotheses using alternative transitions from quadrupedal ancestors towards fossoriality that were inference methods (see 3.1 Phylogenetic hypothesis). not associated with shifts in diversification rate. Furthermore, we Compared to previously published phylogenies, most differences are confirmed (4) a possible evolution of limb-regression in sand-swimmer related to inter-generic relationships within clade B. Amongst these, the genera as an adaptation to a period of general aridification of the Mal­ most remarkable are the recovery of the genus Grandidierina as the sister agasy climate around 20 Ma (Wells, 2003), suggesting a potential in­ lineage to the clade containing Amphiglossus, Brachyseps, Flexiseps, fluence of Madagascar’s bioclimatic zonation on the diversification of Pygomeles, and Voeltzkowia (similarly to the phylogenetic hypothesis the group. Altogether, these results support an evolutionary scenario in proposed in Schmitz et al., 2005), and the recovery of the genus which the in situ diversification process of Malagasy Scincinae consisted Voeltzkowia as the sister lineage to the clade containing Amphiglossus and of a constant accumulation of lineages through time and was determined Pygomeles. This work also included 21 candidate species, some of which by a combination of intrinsic and extrinsic factors, although none of were molecularly characterised for the first time (Supplementary them emerged as the main diversification driver of the group. In the Tables S2–S3). Amongst these, Madascincus sp. CaNew from Betampona 11 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 and Madascincus sp. “beaus” from Andasibe are likely representative of ascribed to diversity-dependent dynamics (Morlon, 2014; Morlon et al., the Madascincus nanus group (Andreone and Greer, 2002), which might 2010). Limited ecological resources may cause speciation and extinction qualify to be considered as a separate genus based on their phylogenetic rates to become equal and stabilise diversity to the carrying capacity of positions and marked genetic and morphological differences (Andreone the system (Etienne et al., 2012, 2016). Similar to what has been and Greer, 2002; Miralles et al., 2016b; Miralles and Vences, 2013). observed in Malagasy pseudoxyrhophiine snakes, we consistently The Malagasy vertebrate fossil record of the Cenozoic is largely detected diversity-dependence in the evolutionary history of Malagasy fragmented (Krause et al., 1997, 1999; Godfrey et al., 2020), limiting the Scincinae across all three tested methods (Morlon et al., 2010; Etienne use of fossil calibrations in dated phylogenies. To circumvent this limi­ et al., 2012, 2016; Burbrink et al., 2019). However, considering that tation, either secondary calibrations from existing large and dated slowdowns in macroevolutionary rates due to diversity-dependent dy­ phylogenies (e.g., as from Crottini et al. 2012; Zheng and Wiens, 2016) namics can be mistaken with alternative mechanisms (see Condamine or fossil data from outside Madagascar to calibrate selected outgroups et al., 2019; Etienne et al., 2016; Etienne and Rosindell, 2012; Harmon could be used (as in Blair et al., 2015; Burbrink et al., 2019). We chose to and Harrison, 2015; Moen and Morlon, 2014; Morlon, 2014; Pannetier use primary fossil calibrations on the outgroups of a reduced dataset et al., 2020), future studies accounting for additional non-phylogenetic (Supplementary Table S5) and used the resulting divergence time esti­ data types (e.g., geographical and ecological data) may help clarify the mates of the ingroup crown age as secondary calibration on the full clade diversity dynamics (Moen and Morlon, 2014; Pannetier et al., dataset. This approach received high statistical support and produced 2020). low uncertainty in the estimated divergence times (Fig. 3A), as also Alternative factors explaining this slowdown may be related to demonstrated by the consistency of our time estimates with estimates environmental and climatic changes that Madagascar underwent over available from other dated phylogenies that used different calibration the last million years. For instance, it has been proposed that Malagasy points (e.g., we estimated the crown node of the group at 49.8 Ma vs 47 grasslands originated during the Late Miocene and Pliocene (Hackel Ma in Crottini et al., 2012 and 44.17 Ma in Zheng and Wiens, 2016). The et al., 2017). This suggests that a major change in Malagasy biota is 95% HPDI for the same node also overlaps with the range reported in the concordant in time with the decreasing diversification rate observed for only other species-level dated phylogeny of the Malagasy Scincinae and Malagasy Scincinae. Considering that most species of this group are that used a secondary calibration approach to infer the crown age from forest dwellers (Glaw and Vences, 2007), the hypothesised expansion of the divergence times estimated in Crottini et al. (2012) (Miralles et al., grasslands may have had a significant impact on their diversification. A 2015). The average value estimated in that study is more recent (37.4 recent decline in diversification rates might also be consistent with a Ma), probably due to the absence of any representative of the basal protracted speciation model (Rosindell et al., 2010) whereby recent genus Pseudoacontias. Similarly, the average values of divergence times speciation events might not have reached a sufficient level of differen­ at the inter-generic level reported in Miralles et al. (2015) are slightly tiation to be reflected in current taxonomy. However, lineages included younger, although in agreement with the timing of the cladogenetic in the phylogenetic dataset were identified based on the interspecific events reported in our dated phylogeny (Fig. 3A). genetic distance thresholds proposed for Malagasy skinks (Nagy et al., 2012; see 2.1 Samples and dataset; Supplementary Tables S2–S3). This 4.2. Diversification rates approach likely helped to capture several incipient and recent speciation events even for taxa that, after a thorough investigation of their taxo­ The diversification rate was found to be constant with the exception nomic status with the generation of additional data (e.g., morphological of a sharp decrease around 6–4 Ma (Fig. 3B), at the boundary between analyses and surveys of the genetic structure based on a denser in­ the Miocene and Pliocene, towards the present. The extinction rate dividuals’ sampling) will not qualify to be formally described as new remained stable and almost null throughout the evolutionary history of species (i.e., Miralles et al., 2016b). the group (Supplementary Fig. S13). In the absence of a fossil record, Finally, and differently from what has been observed in other groups molecular phylogenies including only extant species are known to pro­ (Raxworthy and Nussbaum, 1995; Wilmè et al., 2006; Yoder and duce unrealistically low estimates of extinction rate (see Quental and Heckman, 2006; Boumans et al., 2007; Wollenberg et al., 2008; Vences Marshall, 2010; Rabosky, 2010; Stadler, 2013; Silvestro et al., 2018; et al., 2009), the estimation of divergence times (Fig. 3A) suggested that Louca and Pennell, 2021). The analyses we performed with ACDC sug­ Pleistocene climatic oscillations did not play a strong role in the diver­ gested that alternative extinction rate scenarios of higher rate magni­ sification of Malagasy Scincinae. A possible explanation for this might be tudes and with mostly constant or increasing trends towards the present that a relevant part of the diversification of the group (i.e., at the inter- might be part of the set of congruent models of our diversification hy­ generic level) was influenced by ancient large-scale climatic changes pothesis given, however, the same temporal trend and magnitude in net- that occurred over the long 50-million evolutionary history of the group, diversification rate (Supplementary Figs 14–17). Such analyses as our biogeographic reconstruction and the observation that most demonstrate the importance of exploring alternative congruent models genera are currently restricted to single biomes would suggest (Fig. 5; following the theoretical concepts introduced by Louca and Pennell see 4.4 Fossorial evolution and biogeography). An ancient influence of (2020) to identify diversification dynamics that would be difficult to climate on diversification could have determined climatic niche uncover using other approaches. The stable diversification rate observed specialisation following adaptive mechanisms and smaller distributional across most of the evolutionary history of this group stands in contrast ranges of species as a consequence of non-adaptive allopatric mecha­ with empirical evidence and theoretical assumptions of an initial burst nisms. These conditions might have made these animals less susceptible in speciation that is expected in evolutionary radiations, especially in to Pleistocene climatic oscillations compared to other Malagasy radia­ island systems where the access to novel geographic areas and habitats tions for what concerns novel opportunities of speciation, although an could trigger diversification (Schluter, 2000; Gavrilets and Losos, 2009; influence on phylogeographic patterns within species cannot be Jønsson et al., 2012; Simões et al., 2016). Interestingly, a recent study excluded (Everson et al., 2020). focusing on the evolutionary history of the most diverse Malagasy squamate radiation (snakes of the endemic Pseudoxyrhophiinae sub­ 4.3. Body size evolution family; Burbrink et al., 2019), using a similar methodological approach, also found a constant trend in diversification rate from the origin of the Our analyses suggested that body size of Malagasy Scincinae evolved group, followed by a sharp decline. Tellingly, this decrease in diversi­ homogeneously across clades, and we found no support for an influence fication rate is concurrent in time to the slowdown we detected for of body size on diversification rate in this group (Fig. 4A; Supplementary Malagasy Scincinae. Patterns of declining speciation rate towards the Figs S26–S32). With the exception of some species-poor genera (e.g., present are often recovered in molecular phylogenies and are frequently Amphiglossus and Pygomeles), the other clades did not evolve 12 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 unidirectionally towards larger or smaller sizes and exhibit large vari­ adaptations to arid and sandy environments where limbless species ability ranging from small to relatively large body sizes (Fig. 4A; Sup­ show faster subterranean locomotion than quadrupedal species (Gri­ plementary Table S12). zante et al., 2012; Bergmann et al., 2020; Morinaga and Bergmann, Body size frequency distributions in amphibians and squamates show 2020). Interestingly the three sand-swimmer genera of clade B (Gran­ considerable variation both at continental and global scales (e.g., Meiri, didierina, Voeltzkowia, and Pygomeles; Fig. 1G–H; Fig. 1J–L) all evolved 2008; Feldman et al., 2016; Amado et al., 2021). Small-sized species are during the Early-Middle Miocene (between ca. 23 and 14 Ma; as previ­ usually more frequent than large-sized species (e.g., Meiri, 2008; Feld­ ously observed by Miralles et al., 2015), while the origin of the litter- man et al., 2016; Amado et al., 2021) and associated with higher species burrower genera Pseudoacontias and Paracontias, which are restricted richness in analysed clades (Zimkus et al., 2012). However, a relation­ to humid forests, is older (ca. 50–30 Ma). Different divergence times of ship between body size evolution and diversification rate has rarely been these two groups may reflect the age of the environments they colonised, found (e.g., Feldman et al., 2016), while in most cases such a relation­ as suggested by our ancestral range reconstruction (Fig. 5). Forested ship was not proven (Burbrink et al., 2012; Rabosky et al., 2014a; biomes (dry deciduous and rainforest) were reconstructed as the oldest Zimkus et al., 2012; Wollenberg Valero et al., 2019). This might be ancestral ranges, where the diversification of the group likely started. because in phylogenies with extensive taxonomic and geographic sam­ The oldest reconstruction of the sub-arid spiny thicket biome dated to pling the effect of body size evolution on diversification can be masked ca. 18 Ma (Fig. 5), which could therefore match a time of general ari­ by major ecomorphological differences that can arise at the beginning of dification of Malagasy climate. It is during this time that the south- the evolutionary history of a group, with a strong impact on their early western sub-arid spiny thicket biome (where most Grandidierina spe­ diversification but with limited variation added subsequently (e.g., cies occur) may have established, and with arid sandy environments Crouch and Ricklefs, 2019). appearing within the dry deciduous biome, where the genus Voeltzkowia In the large anuran family of mantellid frogs from Madagascar, and some species of the other two sand-swimmer genera (Grandidierina clades with small body sizes were the most species-rich, tended to have and Pygomeles) are now distributed. Due to the latitudinal position of higher mitochondrial substitution rates and smaller range sizes (Wol­ Madagascar in the high-pressure arid belt at the boundary between the lenberg et al., 2011), and the lower vagility of these smaller species Cretaceous and Paleogene, the Malagasy climate was suggested to have could have favoured inter-populational divergence and speciation been mostly arid at that time, with the sub-arid spiny thicket considered (Pabijan et al., 2012; Wollenberg Valero, 2015; Wollenberg Valero et al., the oldest contemporary biome of the island (Wells, 2003). As 2019). Given the correlation between body size and range size in Mal­ Madagascar drifted north, it has been proposed that the island climate agasy squamates (Brown et al., 2016), a similar pattern could also be became more humid, and this climatic change triggered the rise of expected for the Scincinae, although we did not identify a relationship forested biomes and the contraction of the spiny thicket biome to its between body size evolution and diversification rate. Similarly, we can current extension in the southern part of the country (Wells, 2003). exclude a correlation between SVL and a specific ecomorphological However, evidence from divergence time estimates in plants suggests an adaptation (i.e., fossorial lifestyle evolution), due to the occurrence of older origin of dry and humid forests (see Buerki et al., 2013 for a re­ body size extremes in both analysed ecomorphological adaptations (i.e., view). In line with this proposal, Ohba et al. (2016) suggested that mesic fossorial and non-fossorial lifestyle). habitats (similar to the climatic conditions of the extant subhumid forest of the central highlands) were spread over most of the island during the 4.4. Fossorial evolution and biogeography early Paleogene, and that dry deciduous forests were present only in parts of the west and north-west of Madagascar, while the sub-arid spiny Our analyses suggested multiple, independent origins for limb thicket and the rainforest were likely not yet existing (Ohba et al., 2016). regression (Fig. 4B) and similar diversification rates between fossorial Assuming that our study group is indicative of the environment in which and non-fossorial species (Supplementary Table S14). Our ancestral it currently occurs, our results support the hypothesis that the arid state reconstruction is consistent with previous findings in suggesting biome may be younger, and that the rainforest biome may be older than multiple parallel transitions from fully quadrupedal to limb-reduced traditionally thought (Wells, 2003; Ohba et al., 2016) forms (Miralles et al., 2015). The opposite scenario has been rarely Overall, our ancestral range reconstruction suggested that Mada­ hypothesised within squamates (Ebel et al., 2020) and normally consists gascar’s biomes and past large-scale climatic changes played a role in in re-evolution of digits rather than the entire limbs (e.g., Kohlsdorf and the diversification of Malagasy Scincinae (Fig. 5). Today’s asymmetric Wagner, 2006; Wagner et al., 2018). These findings are also consistent distribution of species across the island bioclimatic zonation supports with Wiens et al. (2006), who found that, at the global level, most limb- the hypothesis that different biomes contributed differently to the in situ reduced fossorial clades evolved from fully quadrupedal ancestors diversification of the group. Forested biomes were reconstructed as the rather than from species with a similar body plan and non-fossorial oldest ranges (Fig. 5), and based on current knowledge they host most of lifestyle (arboreal, aquatic, saxicolous), which would imply a lower present-day diversity and endemism. Rainforest habitats are currently number of evolutionary steps. In addition, Wiens et al. (2006) found that home to 34 species (52% of total analysed diversity), 11 of which are colonisations of new landmasses by limb-reduced taxa are extremely found only in this biome. Considering that a distributional pattern with rare, probably due to their lower dispersal abilities. most species richness concentrated in the eastern rainforest is common The low vagility of fossorial species might be related to the generally across several Malagasy vertebrates (Lees and Colwell, 2007; Brown small distribution ranges of these animals (Wiens et al., 2006; Lee et al., et al., 2016), different authors hypothesised a role of this biome as 2013), a condition that has been suggested to facilitate inter- species pump (Crottini et al., 2012; Brown et al., 2016), although such a populational genetic differentiation and increase diversification rates role has not been yet proved in a robust statistical framework. For within clades (Pabijan et al., 2012; Lee et al., 2013; Wollenberg et al., instance, pseudoxyrhophiine snakes display a similar distributional 2011; Wollenberg Valero, 2015; Wollenberg Valero et al., 2019). Other pattern to Malagasy scincines relative to the island bioclimatic zonation evidence from large-scale macroecological studies suggests that fossorial (Burbrink et al., 2019). However, phylogenetic diversity of species as­ evolution is linked to slowdowns in diversification rate because of the semblages was randomly dispersed across biomes, suggesting that all high ecomorphological specialisations associated with this lifestyle, biomes contributed with the same extent to the in situ diversification of which make fossorial species more susceptible to environmental oscil­ the group (Burbrink et al., 2019). Again, the limited knowledge on lations (Bars-Closel et al., 2017; Cyriac and Kodandaramaiah, 2018). Madagascar’s climate at the time it became isolated (ca. 80–65 Ma, Ali However, in our model system, the transition to a fossorial lifestyle was and Aitchison, 2008) and its evolution until the present hampers a better not found to be associated with changes in diversification rates. understanding of the influence of climatic changes and bioclimatic re­ Limb reduction and fossoriality are frequently associated with gions on the diversification of Malagasy Scincinae and other endemic 13 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 radiations of the island. Resources, Writing – review & editing. Andrea Villa: Formal analysis, Resources, Writing – review & editing. Sam Hyde Roberts: Resources, 4.5. Comparison with diversification patterns in other radiations Writing – review & editing. Marcello Mezzasalma: Resources, Writing – review & editing. Alexander Zizka: Methodology, Validation, Writing Examples of decoupled phenotypic evolution and diversification rate – review & editing. Alexandre Antonelli: Conceptualization, Re­ suggested from our exploration of the in situ diversification of Malagasy sources, Writing – review & editing. Angelica Crottini: Conceptuali­ Scincinae have been previously described for Malagasy pseudox­ zation, Methodology, Validation, Investigation, Resources, Writing – yrhophiine snakes and mantellid frogs (Burbrink et al., 2019; Moen original draft, Writing – review & editing, Supervision, Project admin­ et al., 2021). The constant diversification rate for most of the evolu­ istration, Funding acquisition, Visualization. tionary history of Malagasy pseudoxyrhophiine snakes was not influ­ enced by the evolution of their large variety of dietary and habitat preferences (Glaw and Vences, 2007; Burbrink et al., 2019). Similarly, a Declaration of Competing Interest recent investigation of the radiation of mantellid frogs found a decou­ pling between an accelerating body shape rate and similar diversifica­ The authors declare that they have no known competing financial tion rates relative to other analysed frog radiations from outside interests or personal relationships that could have appeared to influence Madagascar (Moen et al., 2021). At the same time, diversification of the work reported in this paper. Malagasy Scincinae is discordant with the patterns proposed for lemurs and vangid birds (Jønsson et al., 2012; Herrera, 2017). In lemurs, a Acknowledgements steady and slightly increasing diversification rate towards the present is coupled with an early burst pattern in body mass evolution and a par­ We thank the Ministère de l’Environment, de l’Ecologie et des For­ allel expansion in the ecological space, concordant with the expectations ests (now Ministère de l’Environnement et du Développement Durable) of adaptive evolution (Herrera, 2017). Vangid birds are a textbook for issuing research, collection and export permits that enabled the example of adaptive radiation (Jønsson et al., 2012). In this group, compilation of the dataset used in this study. The numerous samples diversification rate had an early burst paralleled by an increase in body analysed in this study have been assembled over many years with the size disparity, which was then followed by a diversification slowdown. help of several colleagues whom we would like to acknowledge and Diversification rate had a second spike in more recent times, related to especially Gonçalo M. Rosa, Jean Nöel, Jean Honoré Velo, Honoré Lava, an increase in beak morphology evolution as an adaptation to multiple Georges, D. James Harris, Daniele Salvi, Emanuele Scanarini, Miary foraging techniques (Jønsson et al., 2012). Randriamialisoa, Stefano Faravelli, Vincenzo Mercurio, Fabio Mattioli, The lack of similar studies on a larger number of Malagasy radia­ Devin Edmonds, Iker A. Irisarri, Alexandra Lima, Solohery Rasamison, tions, along with discrepancies on the employed methodological Andolalao Rakotoarison, Emile Rajeriarison, and Anicet Haza. We are frameworks, hinder the comparisons among different radiations (see, for furthermore grateful to Jan Hackel for support in the interpretation of example, the methods employed in Jønsson et al., 2012; Moen et al., biogeographic results. We thank Meike Kondermann and Gaby Keu­ 2021, relative to those applied in this study), and prevent us from necke for help with laboratory work and Oliver Rauhut for providing identifying an unambiguous common in situ diversification pattern for information about the age of one of the fossils used for calibrations. We the evolution of Madagascar’s vertebrate biota (if existing), although a are in debt to Miguel Vences, Massimo Delfino, and Franco Andreone for few common trends can be mentioned. For instance, by comparing providing us with resources for this contribution. This work was carried studies that applied similar methodologies (i.e., on Malagasy scincines, out in the framework of collaboration agreements amongst the authors’ pseudoxyrhophiine snakes, and lemurs; Herrera, 2017, Burbrink et al., institutions and the Department of Animal Biology of the University of 2019), it is possible that these radiations share in their evolutionary Antananarivo and the Ministry of the Environment of the Republic of histories a continuous increasing trend in species diversity with multiple Madagascar. We thank three anonymous reviewers and associate editor intrinsic and extrinsic diversification drivers that were likely responsible R. Alexander Pyron for their comments that helped to improve an earlier for the diversity observed today. However, at least in the Malagasy version of this manuscript. scincines and pseudoxyrhophiine snakes it seems that none of these drivers had a predominant effect over the others. This might be due to some unique characteristics of Madagascar’s natural history, amongst Funding which a biogeographic history consisting in a long-standing isolation starting before the mass extinction at the K-T boundary that depleted Portuguese National Funds through FCT (Fundação para a Ciência e a Malagasy Cretaceous biota and provided plenty of opportunities for new Tecnologia) supported this work with the projects UIDP/50027/2020 colonisers that reached the island during the Cenozoic (Vences et al., and PTDC/BIA-EVL/31254/2017. FCT also funded the contract grant to 2009; Crottini et al., 2012). In addition, the large dimension of this AC [2020.00823.CEECIND/CP1601/CT0003], the PhD fellowship of FB continental island, the presence of a large environmental diversity, [PD/BD/128493/2017] and WC [SFRH/BD/102495/2014], and the broad bioclimatic conditions and a complex and structured biota, along post-doctoral fellowship of JMP [SFRH/BPD/111015/2015]. AV is with the occurrence of large-scale climatic changes over the time likely funded by the Agencia Española de Investigación, with a Juan de la played an important role in the diversification of Malagasy biota (Vences Cierva-Formación grant [FCI2019-039443-I/AEI/10.13039/ et al., 2009; Crottini et al., 2012). 501100011033]. AA acknowledges financial support from the Swedish Research Council [2019-05191] and the Royal Botanic Gardens, Kew. CRediT authorship contribution statement DS received funding from the Swiss National Science Foundation [PCEFP3_187012] and from the Swedish Research Council [VR: 2019- Francesco Belluardo: Conceptualization, Methodology, Validation, 04739]. Formal analysis, Investigation, Writing – original draft, Writing – review & editing, Visualization. A. Jesus Muñoz-Pajares: Conceptualization, Methodology, Software, Writing – review & editing, Resources. Data accessibility Aurélien Miralles: Conceptualization, Resources, Writing – review & editing. Daniele Silvestro: Conceptualization, Methodology, Valida­ Research data are available on the Mendeley Data platform (https:// tion, Writing – review & editing. Walter Cocca: Methodology, Re­ doi.org/10.17632/xfx2yhw949.2), which include datasets, input files sources, Writing – review & editing. Fanomezana Mihaja Ratsoavina: and R scripts of the analyses. 14 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 Appendix A. Supplementary material Fenner, A.L., Ford, S., Ganesh, S.R., Gardner, M.G., Geissler, P., Gillespie, G.R., Glaw, F., Greenlees, M.J., Griffith, O.W., Grismer, L.L., Haines, M.L., Harris, D.J., Hedges, S.B., Hitchmough, R.A., Hoskin, C.J., Hutchinson, M.N., Ineich, I., Supplementary data to this article can be found online at https://doi. Janssen, J., Johnston, G.R., Karin, B.R., Keogh, J.S., Kraus, F., LeBreton, M., org/10.1016/j.ympev.2022.107635. Lymberakis, P., Masroor, R., McDonald, P.J., Mecke, S., Melville, J., Melzer, S., Michael, D.R., Miralles, A., Mitchell, N.J., Nelson, N.J., Nguyen, T.Q., de Campos Nogueira, C., Ota, H., Pafilis, P., Pauwels, O.S.G., Perera, A., Pincheira-Donoso, D., References Reed, R.N., Ribeiro-Júnior, M.A., Riley, J.L., Rocha, S., Rutherford, P.L., Sadlier, R. A., Shacham, B., Shea, G.M., Shine, R., Slavenko, A., Stow, A., Sumner, J., Ali, J.R., 2017. Islands as biological substrates: Classification of the biological Tallowin, O.J.S., Teale, R., Torres-Carvajal, O., Trape, J.F., Uetz, P., Ukuwela, K.D.B., assemblage components and the physical island types. J. Biogeogr. 44, 984–994. Valentine, L., Van Dyke, J.U., van Winkel, D., Vasconcelos, R., Vences, M., https://doi.org/10.1111/jbi.12872. Wagner, P., Wapstra, E., While, G.M., Whiting, M.J., Whittington, C.M., Wilson, S., Ali, J.R., Aitchison, J.C., 2008. Gondwana to Asia: Plate tectonics, paleogeography and Ziegler, T., Tingley, R., Meiri, S., 2021. Conservation status of the world’s skinks the biological connectivity of the Indian sub-continent from the Middle Jurassic (Scincidae): Taxonomic and geographic patterns in extinction risk. Biol. Conserv. through latest Eocene (166–35 Ma). Earth-Science Rev. 88, 145–166. https://doi. 257, 109101 https://doi.org/10.1016/j.biocon.2021.109101. org/10.1016/j.earscirev.2008.01.007. Condamine, F.L., Rolland, J., Morlon, H., 2019. Assessing the causes of diversification Ali, J.R., Meiri, S., 2019. Biodiversity growth on the volcanic ocean islands and the roles slowdowns: Temperature-dependent and diversity-dependent models receive of in situ cladogenesis and immigration: Case with the reptiles. Ecography 42, equivalent support. Ecol. Lett. 22 (11), 1900–1912. https://doi.org/10.1111/ 989–999. https://doi.org/10.1111/ecog.04024. ele.13382. Ali, J.R., Hedges, S.B., 2022. A review of geological evidence bearing on proposed Crottini, A., Dordel, J., Köhler, J., Glaw, F., Schmitz, A., Vences, M., 2009. A multilocus Cenozoic land connections between Madagascar and Africa and its relevance to phylogeny of Malagasy scincid lizards elucidates the relationships of the fossorial biogeography. Earth-Sci. Rev. 232, 104103 https://doi.org/10.1016/j. genera Androngo and Cryptoscincus. Mol. Phylogenet. Evol. 53, 345–350. https://doi. earscirev.2022.104103. org/10.1016/j.ympev.2009.05.024. Amado, T.F., Martinez, P.A., Pincheira-Donoso, D., Olalla-Tárraga, M.Á., Davies, J., Crottini, A., Madsen, O., Poux, C., Strauß, A., Vieites, D.R., Vences, M., 2012. Vertebrate 2021. Body size distributions of anurans are explained by diversification rates and time-tree elucidates the biogeographic pattern of a major biotic change around the the environment. Glob. Ecol. Biogeogr. 30 (1), 154–164. https://doi.org/10.1111/ K-T boundary in Madagascar. Proc. Natl. Acad. Sci. U. S. A. 109, 5358–5363. https:// geb.13206. doi.org/10.1073/pnas.1112487109. Andreone, F., Greer, A.E., 2002. Malagasy scincid lizards: Descriptions of nine new Crouch, N.M.A., Ricklefs, R.E., 2019. Speciation rate is independent of the rate of species, with notes on the morphology, reproduction and taxonomy of some evolution of morphological size, shape, and absolute morphological specialization in previously described species (Reptilia, Squamata: Scincidae). J. Zool. 258, 139–181. a large clade of birds. Am. Nat. 193 (4), E78–E91. https://doi.org/10.1086/701630. https://doi.org/10.1017/S0952836902001280. Cyriac, V.P., Kodandaramaiah, U., 2018. Digging their own macroevolutionary grave: Antonelli, A., Zizka, A., Carvalho, F.A., Scharn, R., Bacon, C.D., Silvestro, D., Condamine, Fossoriality as an evolutionary dead end in snakes. J. Evol. Biol. 31, 587–598. F., 2018. Data and analyses scripts for Antonelli Zizka et al. (0.9). Zenodo. https://d https://doi.org/10.1111/jeb.13248. oi.org/10.5281/zenodo.1149127. Darriba, D., Taboada, G.L., Doallo, R., Posada, D., 2012. jModelTest 2: More models, new Bars-Closel, M., Kohlsdorf, T., Moen, D.S., Wiens, J.J., 2017. Diversification rates are heuristics and parallel computing. Nat. Methods 9 (8), 772. https://doi.org/ more strongly related to microhabitat than climate in squamate reptiles (lizards and 10.1038/nmeth.2109. snakes). Evolution 71, 2243–2261. https://doi.org/10.1111/evo.13305. Darwin, C., Murray, J. (Eds.), 1859. On the Origin of Species by Means of Natural Beaulieu, J.M., O’Meara, B.C., 2016. Detecting hidden diversification shifts in models of Selection, Or, the Preservation of Favoured Races in the Struggle for Life. John trait-dependent speciation and extinction. Syst. Biol. 65, 583–601. https://doi.org/ Murray, Albemarle Street, London. 10.1093/sysbio/syw022. de Vienne, D.M., Giraud, T., Martin, O.C., 2007. A congruence index for testing Bergmann, P.J., Morinaga, G., Freitas, E.S., Irschick, D.J., Wagner, G.P., Siler, C.D., 2020. topological similarity between trees. Bioinformatics 23, 3119–3124. https://doi.org/ Locomotion and palaeoclimate explain the re-evolution of quadrupedal body form in 10.1093/bioinformatics/btm500. Brachymeles lizards. Proceedings. Biol. Sci. 287, 20201994. https://doi.org/ Ebel, R., Müller, J., Ramm, T., Hipsley, C., Amson, E., 2020. First evidence of convergent 10.1098/rspb.2020.1994. lifestyle signal in reptile skull roof microanatomy. BMC Biol. 18, 185. https://doi. Blair, C., Noonan, B.P., Brown, J.L., Raselimanana, A.P., Vences, M., Yoder, A.D., 2015. org/10.1186/s12915-020-00908-y. Multilocus phylogenetic and geospatial analyses illuminate diversification patterns Erens, J., Miralles, A., Glaw, F., Chatrou, L.W., Vences, M., 2017. Extended molecular and the biogeographic history of Malagasy endemic plated lizards (Gerrhosauridae: phylogenetics and revised systematics of Malagasy scincine lizards. Mol. Phylogenet. Zonosaurinae). J. Evol. Biol. 28, 481–492. https://doi.org/10.1111/jeb.12586. Evol. 107, 466–472. https://doi.org/10.1016/j.ympev.2016.12.008. Bouckaert, R., Vaughan, T.G., Barido-Sottani, J., Duchêne, S., Fourment, M., Etienne, R.S., Rosindell, J., 2012. Prolonging the past counteracts the pull of the present: Gavryushkina, A., Heled, J., Jones, G., Kühnert, D., De Maio, N., Matschiner, M., Protracted speciation can explain observed slowdowns in diversification. Syst. Biol. Mendes, F.K., Müller, N.F., Ogilvie, H.A., du Plessis, L., Popinga, A., Rambaut, A., 61, 204. https://doi.org/10.1093/sysbio/syr091. Rasmussen, D., Siveroni, I., Suchard, M.A., Wu, C.-H., Xie, D., Zhang, C., Stadler, T., Etienne, R.S., Haegeman, B., Stadler, T., Aze, T., Pearson, P.N., Purvis, A., Phillimore, A. Drummond, A.J., 2019. BEAST 2.5: An advanced software platform for Bayesian B., 2012. Diversity-dependence brings molecular phylogenies closer to agreement evolutionary analysis. PLoS Comput. Biol. 15 (4), e1006650. https://doi.org/ with the fossil record. Proc. R. Soc. B Biol. Sci. 279, 1300–1309. https://doi.org/ 10.1371/journal.pcbi.1006650. 10.1098/rspb.2011.1439. Boumans, L., Vieites, D.R., Glaw, F., Vences, M., 2007. Geographical patterns of deep Etienne, R.S., Pigot, A.L., Phillimore, A.B., 2016. How reliably can we infer diversity- mitochondrial differentiation in widespread Malagasy reptiles. Mol. Phylogenet. dependent diversification from phylogenies? Methods Ecol. Evol. 7 (9), 1092–1099. Evol. 45, 822–839. https://doi.org/10.1016/j.ympev.2007.05.028. https://doi.org/10.1111/2041-210X.12565. Brown, J.L., Sillero, N., Glaw, F., Bora, P., Vieites, D.R., Vences, M., 2016. Spatial Everson, K.M., Jansa, S.A., Goodman, S.M., Olson, L.E., 2020. Montane regions shape biodiversity patterns of Madagascar’s amphibians and reptiles. PLoS ONE 11 (1), patterns of diversification in small mammals and reptiles from Madagascar’s moist e0144076. https://doi.org/10.1371/journal.pone.0144076. evergreen forest. J. Biogeogr. 2020 (47), 2059–2072. https://doi.org/10.1111/ Buerki, S., Devey, D.S., Callmander, M.W., Phillipson, P.B., Forest, F., 2013. Spatio- jbi.13945. temporal history of the endemic genera of Madagascar. Bot. J. Linn. Soc. 171, Feldman, A., Sabath, N., Pyron, R.A., Mayrose, I., Meiri, S., 2016. Body sizes and 304–329. https://doi.org/10.1111/boj.12008. diversification rates of lizards, snakes, amphisbaenians and the tuatara. Glob. Ecol. Burbrink, F.T., Ruane, S., Pyron, R.A., 2012. When are adaptive radiations replicated in Biogeogr. 25, 187–197. https://doi.org/10.1111/geb.12398. areas? Ecological opportunity and unexceptional diversification in West Indian Fitzjohn, R.G., 2010. Quantitative traits and diversification. Syst. Biol. 59, 619–633. dipsadine snakes (Colubridae: Alsophiini). J. Biogeogr. 39, 465–475. https://doi. https://doi.org/10.1093/sysbio/syq053. org/10.1111/j.1365-2699.2011.02621.x. Ganzhorn, J.U., Wilmé, L., Mercier, J.L., 2014. Explaining Madagascar’s biodiversity. In: Burbrink, F.T., Ruane, S., Kuhn, A., Rabibisoa, N., Randriamahatantsoa, B., Scales, I.R. (Ed.), Conservation and Environmental Management in Madagascar. Raselimanana, A.P., Andrianarimalala, M.S.M., Cadle, J.E., Lemmon, A.R., Routledge, London, pp. 17–43. Lemmon, E.M., Nussbaum, R.A., Jones, L.N., Pearson, R., Raxworthy, C.J., 2019. The Gavrilets, S., Losos, J.B., 2009. Adaptive radiation: Contrasting theory with data. Science origins and diversification of the exceptionally rich gemsnakes (Colubroidea: 323, 732–737. https://doi.org/10.1126/science.1157966. Lamprophiidae: Pseudoxyrhophiinae) in Madagascar. Syst. Biol. 68, 918–936. Glaw, F., Vences, M., 2007. A Field Guide to the Amphibians and Reptiles of Madagascar, https://doi.org/10.1093/sysbio/syz026. third ed. Vences & Glaw Verlag, Cologne. Burnham, K.P., Anderson, D.R., 2002. Model Selection and Multimodel Inference: A Godfrey, L.R., Samonds, K.E., Baldwin, J.W., Sutherland, M.R., Kamilar, J.M., Practical Information-Theoretic Approach, second ed. Springer, New York. Allfisher, K.L., 2020. Mid-Cenozoic climate change, extinction, and faunal turnover Caetano, D.S., O’Meara, B.C., Beaulieu, J.M., 2018. Hidden state models improve state- in Madagascar, and their bearing on the evolution of lemurs. BMC Evol. Biol. 20, 97. dependent diversification approaches, including biogeographical models. Evolution https://doi.org/10.1186/s12862-020-01628-1. 72, 2308–2324. https://doi.org/10.1111/evo.13602. Goodman, S.M., Benstead, J.P. (Eds.), 2003. The Natural History of Madagascar. Castresana, J., 2000. Selection of conserved blocks from multiple alignments for their use University of Chicago Press. in phylogenetic analysis. Mol. Biol. Evol. 17, 540–552. https://doi.org/10.1093/ Grizante, M.B., Brandt, R., Kohlsdorf, T., 2012. Evolution of body elongation in oxfordjournals.molbev.a026334. gymnophthalmid lizards: Relationships with climate. PLoS One 7 (11), e49772. Chapple, D.G., Roll, U., Böhm, M., Aguilar, R., Amey, A.P., Austin, C.C., Baling, M., https://doi.org/10.1371/journal.pone.0049772. Barley, A.J., Bates, M.F., Bauer, A.M., Blackburn, D.G., Bowles, P., Brown, R.M., Hackel, J., Vorontsova, M.S., Nanjarisoa, O.P., Hall, R.C., Razanatsoa, J., Malakasi, P., Chandramouli, S.R., Chirio, L., Cogger, H., Colli, G.R., Conradie, W., Couper, P.J., Besnard, G., 2017. Grass diversification in Madagascar: In situ radiation of two large Cowan, M.A., Craig, M.D., Das, I., Datta-Roy, A., Dickman, C.R., Ellis, R.J., 15 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 C3 shade clades and support for a Miocene to Pliocene origin of C4 grassy biomes. May, M.R., Höhna, S., Moore, B.R., 2016. A bayesian approach for detecting the impact J. Biogeogr. 45, 750–761. https://doi.org/10.1111/jbi.13147. of mass-extinction events on molecular phylogenies when rates of lineage Harmon, L.J., Harrison, S., 2015. Species diversity is dynamic and unbounded at local diversification may vary. Methods Ecol. Evol. 7 (8), 947–959. https://doi.org/ and continental scales. Am. Nat. 185 (5), 584–593. https://doi.org/10.1086/ 10.1111/2041-210X.12563. 680859. McPeek, M.A., Brown, J.M., 2007. Clade age and not diversification rate explains species Harvey, M.G., 2017. mgharvey/ES-sim: Initial Release (1.1). Zenodo. https://doi.org/10. richness among animal taxa. Am. Nat. 169, E97–E106. https://doi.org/10.1086/ 5281/zenodo.1067144. 512135. Harvey, M.G., Rabosky, D.L., 2018. Continuous traits and speciation rates: Alternatives Meiri, S., 2008. Evolution and ecology of lizard body sizes. Glob. Ecol. Biogeogr. 17, to state-dependent diversification models. Methods Ecol. Evol. 9 (4), 984–993. 724–734. https://doi.org/10.1111/j.1466-8238.2008.00414.x. https://doi.org/10.1111/2041-210X.12949. Meiri, S., 2017. Oceanic island biogeography: Nomothetic science of the anecdotal. Heaney, L.R., 2000. Dynamic disequilibrium: A long-term, large-scale perspective on the Front. Biogeogr. 9 (1), e32081. https://doi.org/10.21425/F59132081. equilibrium model of island biogeography. Glob. Ecol. Biogeogr. 9, 59–74. https:// Miller, M.A., Pfeiffer, W., Schwartz, T., 2010. Creating the CIPRES science gateway for doi.org/10.1046/j.1365-2699.2000.00163.x. inference of large phylogenetic trees. Gateway Computing Environments Workshop Herrera, J.P., 2017. Testing the adaptive radiation hypothesis for the lemurs of (GCE), New Orleans, pp. 1–8. https://doi.org/10.1109/GCE.2010.5676129. Madagascar. Mol. Phylogenet. Evol. 4, 16104. https://doi.org/10.1098/ Miralles, A., Vences, M., 2013. New metrics for comparison of taxonomies reveal striking rsos.161014. discrepancies among species delimitation methods in Madascincus lizards. PLoS ONE Höhna, S., 2013. Fast simulation of reconstructed phylogenies under global time- 8 (7), e68242. https://doi.org/10.1371/journal.pone.0068242. dependent birth–death processes. Bioinformatics 29, 1367–1374. https://doi.org/ Miralles, A., Köhler, J., Glaw, F., Vences, M., 2011a. A molecular phylogeny of the 10.1093/bioinformatics/btt153. “Madascincus polleni species complex”, with description of a new species of scincid Höhna, S., Kopperud, B.T., Magee, A.F., 2022. ACDC: Analysis of congruent lizad from the coastal dune area of northern Madagascar. Zootaxa 2876, 1–16. diversification classes. bioRxiv. https://doi.org/10.1101/2022.01.12.476142. https://doi.org/10.11646/zootaxa.2876.1.1. Jeffreys, H., 1961. The Theory of Probability. Oxford University Press, Oxford. Miralles, A., Köhler, J., Vieites, D.R., Glaw, F., Vences, M., 2011b. Hypotheses on rostral Jønsson, K.A., Fabre, P.H., Fritz, S.A., Etienne, R.S., Ricklefs, R.E., Jørgensen, T.B., shield evolution in fossorial lizards derived from the phylogenetic position of a new Fjeldså, J., Rahbek, C., Ericson, P.G.P., Woog, F., Pasquet, E., Irestedt, M., 2012. species of Paracontias (Squamata, Scincidae). Org. Divers. Evol. 11, 135–150. Ecological and evolutionary determinants for the adaptive radiation of the https://doi.org/10.1007/s13127-011-0042-6. Madagascan vangas. Proc. Natl. Acad. Sci. U. S. A. 109, 6620–6625. https://doi.org/ Miralles, A., Raselimanana, A.P., Rakotomalala, D., Vences, M., Vieites, D.R., 2011c. 10.1073/pnas.1115835109. A new large and colorful skink of the genus Amphiglossus from Madagascar revealed Katoh, K., Standley, D.M., 2013. MAFFT Multiple Sequence Alignment Software version by morphology and multilocus molecular study. Zootaxa 2918, 47–67. https://doi. 7: Improvements in performance and usability. Mol. Biol. Evol. 30, 772–780. org/10.11646/zootaxa.2918.1.5. https://doi.org/10.1093/molbev/mst010. Miralles, A., Anjeriniaina, M., Hipsley, C.A., Müller, J., Glaw, F., Vences, M., 2012. Kisel, Y., Barraclough, T.G., 2010. Speciation has a spatial scale that depends on levels of Variations on a bauplan: Description of a new Malagasy “mermaid skink” with gene flow. Am. Nat. 175, 316–334. https://doi.org/10.1086/650369. flipper-like forelimbs only (Scincidae, Sirenoscincus Sakata & Hikida, 2003). Köhler, J., Vieites, D.R., Glaw, F., Kaffenberger, N., Vences, M., 2009. A further new Zoosystema 34, 701–719. https://doi.org/10.5252/z2012n4a3. species of limbless skink, genus Paracontias, from eastern Madagascar. African J. Miralles, A., Hipsley, C.A., Erens, J., Gehara, M., Rakotoarison, A., Glaw, F., Müller, J., Herpetol. 58, 98–105. https://doi.org/10.1080/21564574.2009.9650029. Vences, M., 2015. Distinct patterns of desynchronized limb regression in Malagasy Köhler, J., Vences, M., Erbacher, M., Glaw, F., 2010. Systematics of limbless scincid Scincine lizards (Squamata, Scincidae). PLoS One 10 (6), e0126074. https://doi.org/ lizards from northern Madagascar: Morphology, phylogenetic relationships and 10.1371/journal.pone.0126074. implications for classification (Squamata: Scincidae). Org. Divers. Evol. 10, Miralles, A., Jono, T., Mori, A., Gandola, R., Erens, J., Köhler, J., Glaw, F., Vences, M., 147–159. https://doi.org/10.1007/s13127-010-0011-5. 2016a. A new perspective on the reduction of cephalic scales in fossorial legless Kohlsdorf, T., Wagner, G.P., 2006. Evidence for the reversibility of digit loss: A skinks (Squamata, Scincidae). Zool. Scr. 45, 380–393. https://doi.org/10.1111/ phylogenetic study of limb evolution in Bachia (Gymnophthalmidae: Squamata). zsc.12164. Evolution 60, 1896–1912. https://doi.org/10.1554/06-056.1. Miralles, A., Köhler, J., Glaw, F., Vences, M., 2016b. Species delimitation methods put Krause, D.W., Hartman, J.H., Wells, N.A., 1997. Late Cretaceous vertebrates from into taxonomic practice: Two new Madascincus species formerly allocated to Madagascar, Implications for biotic changes in deep time. In: Goodman, S.M., historical species names (Squamata, Scincidae). Zoosystematics Evol. 92, 257–275. Patterson, B.D. (Eds.), Natural Change and Human Impact in Madagascar. https://doi.org/10.3897/zse.92.9945. Smithsonian Institution Press, Washington, pp. 3–43. Moen, D., Morlon, H., 2014. Why does diversification slow down? Trends Ecol. Evol. 29, Krause, D.W., Rogers, R.R., Forster, C.A., Hartman, J.H., Buckley, G.A., Sampson, S.D., 190–197. https://doi.org/10.1016/j.tree.2014.01.010. 1999. The late Cretaceous vertebrate fauna of Madagascar: Implications for Moen, D.S., Ravelojaona, R.N., Hutter, C.R., Wiens, J.J., 2021. Testing for adaptive Gondwanan paleobiogeography. GSA Today 9, 1–7. radiation: A new approach applied to Madagascar frogs. Evolution 75, 3008–3025. Lack, D., 1947. Darwin’s Finches. Cambridge University Press, Cambridge. https://doi.org/10.1111/evo.14328. Landis, M.J., Matzke, N.J., Moore, B.R., Huelsenbeck, J.P., 2013. Bayesian analysis of Moore, B.R., Höhna, S., May, M.R., Rannala, B., Huelsenbeck, J.P., 2016. Critically biogeography when the number of areas is large. Syst. Biol. 62, 789–804. https:// evaluating the theory and performance of Bayesian analysis of macroevolutionary doi.org/10.1093/sysbio/syt040. mixtures. Proc. Natl. Acad. Sci. U. S. A. 113, 9569–9574. https://doi.org/10.1073/ Lanfear, R., Frandsen, P.B., Wright, A.M., Senfeld, T., Calcott, B., 2017. PartitionFinder 2: pnas.1518659113. New methods for selecting partitioned models of evolution for molecular and Morinaga, G., Bergmann, P.J., 2020. Evolution of fossorial locomotion in the transition morphological phylogenetic analyses. Mol. Biol. Evol. 34, 772–773. https://doi.org/ from tetrapod to snake-like in lizards. Proc. R. Soc. B Biol. Sci. 287, 20200192. 10.1093/molbev/msw260. https://doi.org/10.1098/rspb.2020.0192. Lee, M.S.Y., Skinner, A., Camacho, A., 2013. The relationship between limb reduction, Morlon, H., 2014. Phylogenetic approaches for studying diversification. Ecol. Lett. 17 body elongation and geographical range in lizards (Lerista, Scincidae). J. Biogeogr. (4), 508–525. https://doi.org/10.1111/ele.12251. 40 (7), 1290–1297. https://doi.org/10.1111/jbi.12094. Morlon, H., Potts, M.D., Plotkin, J.B., 2010. Inferring the dynamics of diversification: A Lees, D.C., Colwell, R.K., 2007. A strong Madagascan rainforest MDE and no coalescent approach. PLoS Biol. 8 (9), e1000493. https://doi.org/10.1371/journal. equatorward increase in species richness: Re-analysis of “The missing Madagascan pbio.1000493. mid-domain effect”, by Kerr J.T., Perring M. & Currie D.J. (Ecology Letters 9: Morlon, H., Parsons, T.L., Plotkin, J.B., 2011. Reconciling molecular phylogenies with 149–159, 2006). Ecol. Lett. 10, E4–E8. https://doi.org/10.1111/j.1461- the fossil record. Proc. Natl. Acad. Sci. U. S. A. 108, 16327–16332. https://doi.org/ 0248.2007.01040.x. 10.1073/pnas.1102543108. Lewitus, E., Morlon, H., 2016. Characterizing and comparing phylogenies from their Morlon, H., Lewitus, E., Condamine, F.L., Manceau, M., Clavel, J., Drury, J., 2016. laplacian spectrum. Syst. Biol. 65, 495–507. https://doi.org/10.1093/sysbio/ RPANDA: An R package for macroevolutionary analyses on phylogenetic trees. syv116. Methods Ecol. Evol. 7, 589–597. https://doi.org/10.1111/2041-210X.12526. Losos, J.B., Ricklefs, R.E., 2009. Adaptation and diversification on islands. Nature 457, Morlon, H., Hartig, F., Robin, S., 2020. Prior hypotheses or regularization allow inference 830–836. https://doi.org/10.1038/nature07893. of diversification histories from extant timetrees. bioRxiv. https://doi.org/10.110 Losos, J., Schluter, D., 2000. Analysis of an evolutionary species–area relationship. 1/2020.07.03.185074. Nature 408, 847–850. https://doi.org/10.1038/35048558. Muñoz-Pajares, A.J., Belluardo, F., Cocca, W., Crottini, A, 2019. PipeLogeny: An Louca, S., Pennell, M.W., 2020. Extant timetrees are consistent with a myriad of automated pipeline for phylogenetic reconstruction. Accessed 1st April 2020. diversification histories. Nature 580, 502–505. https://doi.org/10.1038/s41586- https://sites.google.com/site/pipelogenydownload/. 020-2176-1. Nagy, Z.T., Sonet, G., Glaw, F., Vences, M., 2012. First large-scale DNA barcoding Louca, S., Pennell, M.W., 2021. Why extinction estimates from extant phylogenies are so assessment of reptiles in the biodiversity hotspot of Madagascar, based on newly often zero. Curr. Biol. 31, 3168–3173. https://doi.org/10.1016/j.cub.2021.04.066. designed COI primers. PLoS ONE 7 (3), e34506. https://doi.org/10.1371/journal. MacArthur, R.H., Wilson, E.O., 1963. An equilibrium theory of insular zoogeography. pone.0034506. Evolution 17, 373–387. https://doi.org/10.1111/j.1558-5646.1963.tb03295.x. Ohba, M., Samonds, K.E., LaFleur, M., Ali, J.R., Godfrey, L.R., 2016. Madagascar’s Maddison, W.P., Midford, P.E., Otto, S.P., 2007. Estimating a binary character’s effect on climate at the K/P boundary and its impact on the island’s biotic suite. Palaeogeogr. speciation and extinction. Syst. Biol. 56, 701–710. https://doi.org/10.1080/ Palaeoclimatol. Palaeoecol. 441, 688–695. https://doi.org/10.1016/j. 10635150701607033. palaeo.2015.10.028. Matzke, N.J., 2018. BioGeoBEARS: BioGeography with bayesian (and likelihood) Pabijan, M., Wollenberg, K.C., Vences, M., 2012. Small body size increases the regional evolutionary analysis with R scripts. Accessed 1st February 2021. https://doi. differentiation of populations of tropical mantellid frogs (Anura: Mantellidae). org/10.5281/zenodo.1478250. J. Evol. Biol. 25, 2310–2324. https://doi.org/10.1111/j.1420-9101.2012.02613.x. 16 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 Pannetier, T., Martinez, C., Bunnefeld, L., Etienne, R.S., 2020. Branching patterns in Simpson, G.G., 1953. The Major Features of Evolution. Columbia University Press, New phylogenies cannot distinguish diversity-dependent diversification from time- York. dependent diversification. Evolution 75, 25–38. https://doi.org/10.1111/ Stadler, T., 2013. Recovering speciation and extinction dynamics based on phylogenies. evo.14124. J. Evol. Biol. 26, 1203–1219. https://doi.org/10.1111/jeb.12139. Pennell, M.W., Eastman, J.M., Slater, G.J., Brown, J.W., Uyeda, J.C., Fitzjohn, R.G., Stamatakis, A., 2014. RAxML version 8: A tool for phylogenetic analysis and post- Alfaro, M.E., Harmon, L.J., 2014. Geiger v2.0: An expanded suite of methods for analysis of large phylogenies. Bioinformatics 30, 1312–1313. https://doi.org/ fitting macroevolutionary models to phylogenetic trees. Bioinformatics 30, 10.1093/bioinformatics/btu033. 2216–2218. https://doi.org/10.1093/bioinformatics/btu181. Töpel, M., Zizka, A., Calió, M.F., Scharn, R., Silvestro, D., Antonelli, A., 2016. Presgraves, D.C., Glor, R.E., 2010. Evolutionary biology: Speciation on islands. Curr. SpeciesGeoCoder: Fast categorization of species occurrences for analyses of Biol. 20, 440–442. https://doi.org/10.1016/j.cub.2010.03.032. biodiversity, biogeography, ecology, and evolution. Syst. Biol. 66, 145–151. https:// Quental, T.B., Marshall, C.R., 2010. Diversity dynamics: Molecular phylogenies need the doi.org/10.1093/sysbio/syw064. fossil record. Trends Ecol. Evol. 25, 434–441. https://doi.org/10.1016/j. Toussaint, E., Hendrich, L., Shaverdo, H., Balke, M., 2015. Mosaic patterns of tree.2010.05.002. diversification dynamics following the colonization of Melanesian islands. Sci. Rep. R Core Team, 2021. R: A language and environment for statistical computing. R 5, 16016. https://doi.org/10.1038/srep16016. Foundation for Statistical Computing, Vienna. Accessed 3rd January 2021. https:// Uetz, P., Freed, P., Hošek, J., 2022. The reptile database. Accessed 2nd March 2022. htt www.R-project.org/. p://www.reptile-database.org. Rabosky, D.L., 2010. Extinction rates should not be estimated from molecular Valente, L., Phillimore, A.B., Melo, M., Warren, B.H., Clegg, S.M., Havenstein, K., phylogenies. Evolution 64, 1816–1824. https://doi.org/10.1111/j.1558- Tiedemann, R., Illera, J.C., Thébaud, C., Aschenbach, T., Etienne, R.S., 2020. 5646.2009.00926.x. A simple dynamic model explains the diversity of island birds worldwide. Nature Rabosky, D.L., 2014. Automatic detection of key innovations, rate shifts, and diversity- 579, 92–96. https://doi.org/10.1038/s41586-020-2022-5. dependence on phylogenetic trees. PLoS ONE 9, e89543. https://doi.org/10.1371/j Vences, M., Vieites, D.R., Glaw, F., Brinkmann, H., Kosuch, J., Veith, M., Meyer, A., ournal.pone.0089543. 2003. Multiple overseas dispersal in amphibians. Proc. R. Soc. B Biol. Sci. 270, Rabosky, D.L., Huang, H., 2016. A robust semi-parametric test for detecting trait- 2435–2442. https://doi.org/10.1098/rspb.2003.2516. dependent diversification. Syst. Biol. 65, 181–193. https://doi.org/10.1093/sysbio/ Vences, M., Wollenberg, K.C., Vieites, D.R., Lees, D.C., 2009. Madagascar as a model syv066. region of species diversification. Trends Ecol. Evol. 24, 456–465. https://doi.org/ Rabosky, D.L., Donnellan, S.C., Grundler, M., Lovette, I.J., 2014a. Analysis and 10.1016/j.tree.2009.03.011. visualization of complex Macroevolutionary dynamics: An example from Australian Vieites, D.R., Wollenberg, K.C., Andreone, F., Köhler, J., Glaw, F., Vences, M., 2009. Vast Scincid lizards. Syst. Biol. 63, 610–627. https://doi.org/10.1093/sysbio/syu025. underestimation of Madagascar’s biodiversity evidenced by an integrative Rabosky, D.L., Grundler, M., Anderson, C., Title, P., Shi, J.J., Brown, J.W., Huang, H., amphibian inventory. Proc. Natl. Acad. Sci. U. S. A. 106, 8267–8272. https://doi. Larson, J.G., Kembel, S., 2014b. BAMMtools: An R package for the analysis of org/10.1073/pnas.0810821106. evolutionary dynamics on phylogenetic trees. Methods Ecol. Evol. 5 (7), 701–707. Wagner, G.P., Griffith, O.W., Bergmann, P.J., Bello-Hellegouarch, G., Kohlsdorf, T., Rambaut, A., 2009. FigTree 1.3.1. Institute of Evolutionary Biology, University of Bhullar, A., Siler, C.D., 2018. Are there general laws for digit evolution in Edinburgh, Edinburgh. Accessed 10th May 2018. http://tree.bio.ed.ac.uk/software squamates? The loss and re-evolution of digits in a clade of fossorial lizards /figtree/. (Brachymeles, Scincinae). J. Morphol. 279, 1104–1119. https://doi.org/10.1002/ Rambaut, A., Drummond, A.J., Xie, D., Baele, G., Suchard, M.A., 2018. Posterior jmor.20834. summarization in bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67, 901–904. Warren, B.H., Simberloff, D., Ricklefs, R.E., Aguilée, R., Condamine, F.L., Gravel, D., https://doi.org/10.1093/sysbio/syy032. Morlon, H., Mouquet, N., Rosindell, J., Casquet, J., Conti, E., Cornuault, J., Raxworthy, C.J., Nussbaum, R.A., 1993. Four new species of Amphiglossus from Fernández-Palacios, J.M., Hengl, T., Norder, S.J., Rijsdijk, K.F., Sanmartín, I., Madagascar (Squamata: Scincidae). Herpetologica 49, 326–341. Strasberg, D., Triantis, K.A., Valente, L.M., Whittaker, R.J., Gillespie, R.G., Raxworthy, C.J., Nussbaum, R.A., 1995. Systematics, speciation and biogeography of the Emerson, B.C., Thébaud, C., 2015. Islands as model systems in ecology and dwarf chameleons (Brookesia, Reptilia, Squamata, Chamaeleontidae) of northern evolution: Prospects fifty years after MacArthur-Wilson. Ecol. Lett. 18 (2), 200–217. Madagascar. J. Zool. 235, 525–558. https://doi.org/10.1111/j.1469-7998.1995. https://doi.org/10.1111/ele.12398. tb01767.x. Wells, N.A., 2003. Geology and soils. In: Goodman, S.M., Benstead, J.P. (Eds.), The Ree, R.H., Smith, S.A., 2008. Maximum likelihood inference of geographic range Natural History of Madagascar. University of Chicago Press, Chicago and London, evolution by dispersal, local extinction, and cladogenesis. Syst. Biol. 57, 4–14. pp. 16–74. https://doi.org/10.1080/10635150701883881. Welt, R.S., Raxworthy, C.J., 2022. Dispersal, not vicariance, explains the biogeographic Revell, L.J., 2012. phytools: An R package for phylogenetic comparative biology (and origin of iguanas on Madagascar. Mol. Phylogenet. Evol. 167, 107345 https://doi. other things). Methods Ecol. Evol. 3, 217–223. https://doi.org/10.1111/j.2041- org/10.1016/j.ympev.2021.107345. 210X.2011.00169.x. Whiting, A.S., Sites Jr., J.W., Bauer, A.M., 2004. Molecular phylogenetics of Malagasy Ronquist, F., 1997. Dispersal-vicariance analysis: A new approach to the quantification of skinks (Squamata: Scincidae). African J. Herpetol. 53, 135–146. https://doi.org/ historical biogeography. Syst. Biol. 46, 195–203. https://doi.org/10.1093/sysbio/ 10.1080/21564574.2004.9635506. 46.1.195. Wiens, J.J., Brandley, M.C., Reeder, T.W., 2006. Why does a trait evolve multiple times Ronquist, F., Teslenko, M., Van Der Mark, P., Ayres, D.L., Darling, A., Höhna, S., within a clade? Repeated evolution of snakelike body form in squamate reptiles. Larget, B., Liu, L., Suchard, M.A., Huelsenbeck, J.P., 2012. MrBayes 3.2: Efficient Evolution 60, 123–141. https://doi.org/10.1554/05-328.1. bayesian phylogenetic inference and model choice across a large model space. Syst. Wilmè, L., Goodman, S.M., Ganzhorn, J.U., 2006. Biogeographic evolution of Biol. 61, 539–542. https://doi.org/10.1093/sysbio/sys029. Madagascar’s microendemic biota. Science 312, 1063–1065. https://doi.org/ Rosindell, J., Cornell, S.J., Hubbell, S.P., Etienne, R.S., 2010. Protracted speciation 10.1126/science.1122806. revitalizes the neutral theory of biodiversity. Ecol. Lett. 13, 716–727. https://doi. Wollenberg Valero, K.C., 2015. Evidence for an intrinsic factor promoting landscape org/10.1111/j.1461-0248.2010.01463.x. genetic divergence in Madagascan leaf-litter frogs. Front. Genet. 6, 155. https://doi. Samonds, K.E., Godfrey, L.R., Ali, J.R., Goodman, S.M., Vences, M., Sutherland, M.R., org/10.3389/fgene.2015.00155. Irwin, M.T., Krause, D.W., 2012. Spatial and temporal arrival patterns of Wollenberg Valero, K.C., Marshall, J.C., Bastiaans, E., Caccone, A., Camargo, A., Madagascar’s vertebrate fauna explained by distance, ocean currents, and ancestor Morando, M., Niemiller, M.L., Pabijan, M., Russello, M.A., Sinervo, B., Werneck, F.P., type. Proc. Natl. Acad. Sci. U. S. A. 109, 5352–5357. https://doi.org/10.1073/ Sites Jr., J.W., Wiens, J.J., Steinfartz, S., 2019. Patterns, mechanisms and genetics of pnas.1113993109. speciation in reptiles and amphibians. Genes 10, 646. https://doi.org/10.3390/ Samonds, K.E., Godfrey, L.R., Ali, J.R., Goodman, S.M., Vences, M., Sutherland, M.R., genes10090646. Irwin, M.T., Krause, D.W., 2013. Imperfect isolation: Factors and filters shaping Wollenberg, K.C., Vieites, D.R., van der Meijden, A., Glaw, F., Cannatella, D.C., Madagascar’s extant vertebrate fauna. PLoS One 8 (4), e62086. https://doi.org/ Vences, M., 2008. Patterns of endemism and species richness in Malagasy cophyline 10.1371/journal.pone.0062086. frogs support a key role of mountainous areas for speciation. Evolution 62, Scantlebury, D.P., 2013. Diversification rates have declined in the Malagasy 1890–1907. https://doi.org/10.1111/j.1558-5646.2008.00420.x. herpetofauna. Proc. R. Soc. B Biol. Sci. 280, 20131109. https://doi.org/10.1098/ Wollenberg, K.C., Vieites, D.R., Glaw, F., Vences, M., 2011. Speciation in little: the role of rspb.2013.1109. range and body size in the diversification of Malagasy mantellid frogs. BMC Evol. Schluter, D., 2000. The Ecology of Adaptive Radiation. Oxford University Press, Oxford. Biol. 11, 217. https://doi.org/10.1186/1471-2148-11-217. Schmitz, A., Brandley, M.C., Mausfeld, P., Vences, M., Glaw, F., Nussbaum, R.A., Yoder, A.D., Heckman, K.L., 2006. Mouse lemur phylogeography revises a model of Reeder, T.W., 2005. Opening the black box: Phylogenetics and morphological ecogeographic constraint in Madagascar. In: Fleagle, J., Lehman, S.M. (Eds.), evolution of the Malagasy fossorial lizards of the subfamily “Scincinae”. Mol. Primate Biogeography: Progress and Prospects. Kluwer Press, New York, Phylogenet. Evol. 34, 118–133. https://doi.org/10.1016/j.ympev.2004.08.016. pp. 255–268. Shea, G.M., 2021. The identity of twelve scincid generic names proposed by Cope in Yoder, A.D., Nowak, M.D., 2006. Has vicariance or dispersal been the predominant 1892 (Squamata: Scincidae). Zootaxa 5057, 577–589. https://doi.org/10.11646/ biogeographic force in Madagascar? Only time will tell. Annu. Rev. Ecol. Evol. Syst. zootaxa.5057.4.8. 37, 405–1301. https://doi.org/10.1146/annurev.ecolsys.37.091305.110239. Silvestro, D., Warnock, R.C.M., Gavryushkina, A., Stadler, T., 2018. Closing the gap between palaeontological and neontological speciation and extinction rate estimates. Nat. Commun. 9, 5237. https://doi.org/10.1038/s41467-018-07622-y. Simões, M., Breitkreuz, L., Alvarado, M., Baca, S., Cooper, J.C., Heins, L., Herzog, K., Lieberman, B.S., 2016. The evolving theory of evolutionary radiations. Trends Ecol. Evol. 31, 27–34. https://doi.org/10.1016/j.tree.2015.10.007. 17 F. Belluardo et al. Molecular Phylogenetics and Evolution 178 (2023) 107635 Zhang, C., Rabiee, M., Sayyari, E., Mirarab, S., 2018. ASTRAL-III: Polynomial time genes and 4162 species. Mol. Phylogenet. Evol. 94, 537–547. https://doi.org/ species tree reconstruction from partially resolved gene trees. BMC Bioinform. 19, 10.1016/j.ympev.2015.10.009. 153. https://doi.org/10.1186/s12859-018-2129-y. Zimkus, B.M., Lawson, L., Loader, S.P., Hanken, J., 2012. Terrestrialization, Zheng, Y., Wiens, J.J., 2016. Combining phylogenomic and supermatrix approaches, and miniaturization and rates of diversification in african puddle frogs (Anura: a time-calibrated phylogeny for squamate reptiles (lizards and snakes) based on 52 Phrynobatrachidae). PLoS One 7 (4), e35118. https://doi.org/10.1371/journal. pone.0035118. 18