dataesr/scientific-paragraphs-categorization huggingface dataset:
We compiled a dataset comprising all individuals of Malagasy Scincinae whose tissue samples were available to us (Supplementary Table  S2 ). These samples were sequenced for at least one of the three mitochondrial markers that are generally used for species identification of Malagasy Scincinae, which allowed us to identify and assign them unequivocally to a mitochondrial interspecific lineage: 16S rRNA 3 ′ terminus (16S), cytochrome oxidase I gene (COI), and NADH dehydrogenase subunit 1 (ND1) (with adjacent tRNAs) genes (see Supplementary Text S1 for details on laboratory work). To do this, we applied the skinks' standard interspecific threshold of 6. 1-6. 5% genetic distance at the COI marker, as proposed in  Nagy et al. (2012) , and used 16S and ND1 to refine species identification when necessary.
We compiled a phylogenetic multi-locus dataset using the matrix of  Erens et al. (2017)  as a starting point. This matrix contained 45 interspecific lineages (either nominal species or candidate species as defined by  Vieites et al. , 2009)  and 13 loci. We augmented this matrix by adding one locus (COI) and sequences from 26 additional lineages (seven of which were identified for the first time; Supplementary Tables  S2-S3 ) and removing three taxa corresponding to intraspecific lineages. Lineage identification in our phylogenetic multi-locus dataset was performed following the procedure described above. Our final molecular dataset included 68 ingroup taxa (47 of the 65 nominal species of Malagasy Scincinae and 21 candidate species; Supplementary Table  S1 ) and 14 loci, which comprised four mitochondrial (12S rRNA, 16S, COI, and ND1 with adjacent tRNAs) and 10 nuclear markers (brain-derived neurotrophic factor (BDNF), oocyte maturation factor (CMOS), phosducin (PDC), recombination activating gene 1 and gene 2 (RAG1 and RAG2), alpha-enolase (ENO1), two non-overlapping segments of sacsin (SACS-A and SACS-B), leucine-rich repeat and WD repeat-containing protein (KIAA1239), and titin (TTN)) (Supplementary Table  S3 ; Supplementary Text S1).
To maximise the number of available loci, we used sequences of multiple congeneric species (chimera sequences) of Cordylus, Eumeces, Tiliqua, and Scincella as outgroups, as in  Crottini et al. (2009)  and  Erens et al. (2017)  (Supplementary Table  S3 ). Newly generated sequences have been deposited in GenBank (OP256430-OP256524, OP263988-OP264059, OP271309-OP271456; Supplementary Table  S3 ). For details on laboratory work and sample amplification see Supplementary Text S1.
Phylogenetic analyses We performed phylogenetic analyses on a concatenated dataset composed of all mitochondrial and nuclear markers available in our dataset (Supplementary Table  S3 ). We used PipeLogeny in R 4. 0. 2 (Muñoz-Pajares et al. , 2019; R Core Team, 2021) , which is a recently developed pipeline that automatises the preparation of the input files for several commonly used phylogenetic tools, which are then run separately in their standalone software applications. We aligned sequences with MAFFT 7. 310  (Katoh and Standley, 2013) , detected misaligned and hypervariable regions of the alignment to be excluded from the phylogenetic analyses with Gblocks 0. 91b  (Castresana, 2000) , identified the best-fit partitioning scheme and models of evolution using Partition-Finder2 2. 1. 1  (Lanfear et al.
, 2017) , and performed Bayesian analyses using MrBayes 3. 2. 7  (Ronquist et al. , 2012). Due to the limited number of available gene fragments and the unavailability of recently collected tissue samples, we excluded A. splendidus and P. menamainty from the phylogenetic analyses (see Supplementary Table  S3  for the phylogenetic dataset). After the preparation of the input files, we ran Gblocks on the GBlocks Server  (Castresana, 2000)  using the least stringent parameter settings. We used PartitionFinder2 on the CIPRES Science Gateway  (Miller et al. , 2010)  to find the best-fit partitioning scheme and models of sequence evolution, which were tested under Bayesian Information Criterion (BIC). We ran MrBayes on CIPRES with two runs of 100-million generations each, setting four independent chains for each run with sampling at every 1,000 th generation and 40% burn-in of the Markov Chain Monte Carlo (MCMC) chain.
We assessed run convergence and posterior distributions with Tracer 1. 7. 1  (Rambaut et al. , 2018)  and considered a value of 200 Effective Sampling Size (ESS) as a minimum threshold to evaluate the fitness of each prior. We used FigTree 1. 4. 4  (Rambaut, 2009)  to visualise the 50%-majority rule consensus tree. Following the procedure described above, we inferred distinct topologies for the concatenated alignments of mitochondrial and nuclear markers, respectively. We ran two independent runs of 70-million generations each with MrBayes for both phylogenies and assessed topological congruence between the mitochondrial, nuclear and full concatenated (mitochondrial and nuclear) topologies with the Icong index  (de Vienne et al. , 2007 ; R script available on GitHub at https://gi thub. com/damiendevienne/icong).
To explore the influence of alternative inference methods on tree topology, we inferred a maximum likelihood (ML) phylogeny with RAxML-HPC2 8. 2. 12  (Stamatakis, 2014)  on the full multi-locus matrix (Supplementary Table  S3 ). We used the same partition scheme estimated with PartitionFinder2 for Bayesian analyses (see above). We ran RAxML on CIPRES setting 1,000 replicates of rapid bootstrap and GTRCAT model. In addition, we inferred a species tree with ASTRAL 5. 7. 8  (Zhang et al. , 2018)  starting from gene trees of the 14 markers available in our dataset. We inferred the gene trees with MrBayes following the same procedure described above for concatenated datasets, with the only exception that we used jModeltest 2. 1. 10 (Darriba et al. , 2012)  on CIPRES to find the best model of sequence evolution for each gene.
Dated phylogeny Given the absence of fossils in the ingroup, we used fossil-based calibrations on the outgroups. We compiled a reduced dataset including four Malagasy Scincinae and seven outgroups (Supplementary Table  S5 ). We selected one species for each of the outgroups (Cordylus, Eumeces, Tiliqua, and Scincella) of the complete species-level dataset, and the remaining outgroups were selected based on the fossil calibrations that were used (Supplementary Table  S5 ); phylogenetic placement and age of fossil calibrations were expert-verified (by A. Villa) (see Supplementary Table  S6  for information on the selected fossils and justification 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 evolution with PartitionFinder2 2.
1. 1, under BIC, on the CIPRES Science Gateway. We estimated divergence times with BEAST 2. 6. 0  (Bouckaert et al. , 2019)  (see Supplementary Text S1 for all analysis details). The height average and 95% Higher Posterior Density Interval (HPDI) of the crown node of the Malagasy Scincinae was used to calibrate the corresponding node on the complete species-level dataset, and we estimated divergence times on this dataset with BEAST (Supplementary Table  S3 ; see Supplementary Text S1 for all analysis details).
Diversification rate analyses We inferred diversification rates from the dated phylogeny pruned from all outgroups. Whenever possible we set a sampling fraction of 76%, corresponding to the percentage of taxa included in the phylogeny relative to the overall number of nominal species and candidate species recognised for the group. Since several models and methods that we employed in the following analyses assume rate homogeneity across lineages, we first checked for possible heterogeneity of diversification rate within the tree with BAMM 2. 5. 0  (Rabosky, 2014). BAMM identifies heterogeneity in the diversification process through time and across branches within a phylogeny using the reversible-jump MCMC (rjMCMC). The priors of the evolutionary rate parameters were set based on the function 'setBAMMpriors' of the package BAMMtools 2.
1. 7 (Rabosky et al. , 2014b)  available on R 4. 0. 2. We tested different priors for the values of the expected number of diversification rate shifts (0. 1, 0. 5, 1, 5, and 10) and for each value we ran two MCMC chains of 10million steps to assess reciprocal convergence, while we ensured single-chain convergence through the ESS (minimum threshold of 200). We set a non-random incomplete taxon sampling, which accounts for the different numbers of included taxa amongst different genera of the group. The output was inspected with the R package BAMMtools 2. 1. 7. We evaluated the statistical support of the number of diversification rate shifts based on Posterior Probability (PP) and Bayes Factor after discarding 15% of the posterior samples as burn-in. Moore et al. (2016)  criticised BAMM highlighting possible issues concerning the likelihood function and the prior models employed in the software.
To account for these uncertainties in diversification rate analyses and the several other analyses in which we used this software (see 2. 5 Trait evolution), we decided to run BAMM along with other tools to highlight potential discordance in the results that might be indicative of technical issues. We further investigated the possible heterogeneity in diversification rates by examining the spectral density profile of the tree computed from the Laplacian graph of its branch lengths  (Lewitus and Morlon, 2016). We obtained the eigenvalues and spectral density profile of the tree with the function 'spectR' implemented in the R package RPANDA 1. 8  (Morlon et al. , 2016). To test the number of modalities that characterise the tree, that can represent the number of heterogeneities in diversification rates, we used the function 'BICompare'.
We compared BIC values associated with the empirical tree and randomly bifurcating trees parametrised from this to provide a significance test for the number of modalities that were previously highlighted. We investigated the occurrence and timing of mass extinction events and tree-wide diversification rate changes through time with the CoMET model within the R package TESS 2. 1. 0 (Höhna, 2013; May et al. , 2016). This method implements a piecewise-constant speciation and extinction rate model with the number of rate shifts through time that is inferred through rjMCMC. We set the function 'tess. analysis' to account for empirical estimates of the hyperpriors of diversification-rates parameters. We ran the rjMCMC for 50-million steps of which the initial 14. 5% was discarded as burn-in. We considered a minimum ESS of 200 to ensure convergence of the chains.
The occurrence of mass extinction events and shifts in speciation and extinction rates was evaluated through Bayes Factor. To further evaluate the TESS results, we investigated the possible time-dependence of diversification rates through a maximum likelihood model-fitting approach. We tested six models reported in  Morlon et al. (2011)  accounting for both time-independent and time-dependent variation of speciation and extinction rates as well as periods of expanding and declining diversity (Supplementary Table  S7 ). These models were fitted with the function 'fit_bd' implemented in the R package RPANDA 1. 8 (Supplementary Table  S7 ). The best-fitting model was chosen based on the corrected Akaike Information Criterion (AICc) weights  (Burnham and Anderson, 2002; Morlon et al. , 2016). After detecting evidence for time-dependence of diversification rates (see 3.
The estimated parameters of the bestfitting model were used in the function 'dd_LR', which implements a bootstrap likelihood ratio test of a diversity-dependent model against a constant-rates birth-death model  (Etienne et al. , 2016). We tested the robustness of DDD results by investigating the possible relationships between diversification rates variation through time and lineage diversity accumulation through the nine models from  Morlon et al. (2010)  (Supplementary Table  S9 ), whose maximum likelihood expressions are based on the coalescent process. We used the 'fit_-coal_cst' and 'fit_coal_var' functions in the R package RPANDA 1. 8. The tested models imply both saturated diversity scenarios with turnover rate (the rate at which an extinct lineage is substituted by an emerging one) and scenarios in which diversification is not diversity-dependent and diversity is not saturated (Supplementary Table  S9 ).
, 2020). Starting from diversification rate estimates from CoMET analyses, we used the R package ACDC  (Höhna et al. , 2022; Louca and Pennell, 2020)  to construct the congruence class of all extinction and speciation rate functions with the same likelihood score and pulled diversification and speciation rates to explore the diversity of rate functions belonging to the class and uncover the array of potential diversification models underlying Malagasy Scincinae diversification. We tested a series of speciation and extinction rate functions (separately from each other) of different models and temporal trends towards the present (increasing and decreasing exponential, linear and episodically changing birth--death models). We computed speciation, extinction, relative extinction and net-diversification rates for each category of tested rate function and compared trends through time with the reference model (CoMET diversification rates estimates) to discern between functions belonging to the same congruence class and those that did not.
Trait evolution 2. 5. 1. Body size evolution We investigated body size evolution using snout-vent length (SVL) as a proxy. We compiled a list of individual SVL measurements of specimens from museum collections and complemented it with published data to obtain a matrix of a total of 60 taxa and 341 SVL measures (Supplementary Table  S10 ). SVL measures were not available for Pseudoacontias CaNew. Since Pseudoacontias spp. are the largest Malagasy Scincinae  (Glaw and Vences, 2007) , and their inclusion or exclusion may have an influence on the analyses, we ran a replicate of these analyses including the SVL of P. menamainty as a proxy value for the genus. We averaged SVL measures over conspecific individuals and logtransformed (Supplementary Table  S11 ). We applied a global sampling fraction of 69% (with Pseudoacontias exclusion) and 70% (with Pseudoacontias inclusion) accounting for missing taxa to the following analyses, whenever possible.
We performed model-testing across different models of continuous trait evolution (Brownian motion, Ornstein-Uhlenbeck process, and Early Burst) with the function 'fitContinuous' in the R package geiger 2. 0. 7  (Pennell et al. , 2014)  (Supplementary Table  S12 ). We performed the analyses on 100 trees randomly sampled from the posterior distribution of the dated phylogeny. We considered the number of trees in which a model resulted as the best-fitting (based on AICc weights;  Burnham and Anderson, 2002)  to evaluate the statistical support of each model over the entire posterior sample and over the phylogenetic inference uncertainty. In addition, we averaged AICc weights over the 100 trees for each tested model. We performed ancestral state reconstruction based on the best-fitting model on the 50%-majority rule consensus tree (see 2.
2 Phylogenetic analyses) using the function 'con-tMap' in the R package phytools 0. 7-70  (Revell, 2012). We investigated possible rate shifts in SVL evolution with BAMM. We used the function 'setBAMMpriors' from BAMMtools to find a suitable set of priors for our dataset and tested different values of the expected number of rate shifts priors for the evolution of this trait (0. 1, 0. 5, 1, 5, and 10). We ran two MCMC chains of 10-million iterations. The output was processed with BAMMtools. We assessed multiple and single-chain convergence as described for BAMM diversification rate analyses (see 2. 4 Diversification rate analyses). We selected the most supported number of trait evolution rate shifts based on the Bayes Factor and we chose the best configuration shift based on the highest Posterior Probability.
We tested the possible correlation between SVL and diversification rates with two tip-rate correlation methods: ES-sim  (Harvey and Rabosky, 2018)  and 'STructured RAte Permutations on Phylogenies' (STRAPP)  (Rabosky and Huang, 2016). These methods are more robust compared to other approaches (e. g. , QuaSSE;  Fitzjohn, 2010)  and in particular with small phylogenies  (Harvey and Rabosky, 2018; Rabosky and Huang, 2016). ES-sim tests the correlation between trait variation across the tips of a phylogeny and diversification rates. It applies the Inverse Equal Splits as summary statistic of branching pattern and tests for deviations of the empirical phylogeny from a null distribution simulated with no evolutionary model  (Harvey and Rabosky, 2018). We ran ES-sim with the R code provided in  Harvey (2017)  on the same 100 trees sampled from the posterior distribution used in previous SVL analyses (see above) and set 1,000 simulations for the null distribution.
We explored fossorial evolution with the 'Hidden State Speciation and Extinction' (HiSSE) method implemented in the R package hisse 1. 9. 8  (Beaulieu and O'Meara, 2016) , which estimates the effect of binary characters on diversification rate within a maximum-likelihood framework for model-testing. Unlike the 'Binary-State Speciation and Extinction' (BiSSE) and other SSE methods  (Maddison et al. , 2007; Fitzjohn, 2010) , HiSSE accounts for additional unmeasured (hidden) states, which may potentially influence diversification and transitions rates amongst states, and implements more realistically complex character-independent null models. These aspects make HiSSE less prone to type I error compared with BiSSE and other SSE methods  (Beaulieu and O'Meara, 2016; Caetano et al. , 2018). We tested 29 models comprising BiSSE-like models (with no hidden trait;  Maddison et al.
Historical biogeography We compiled a dataset of 391 geographic records relating to all species included in our phylogenetic analyses (Supplementary Table  S15 ). We compiled these species distribution records from own field records, GenBank, taxonomic literature or from  Glaw and Vences (2007). We revised all records following current Malagasy Scincinae taxonomy (Supplementary Table  S15 ). We performed historical biogeography analyses using the 5-region bioclimatic classification (rainforest, dry deciduous forest, sub-humid forest, sub-arid spiny thicket, and montane thickets) modified from  Yoder and Nowak (2006)  by merging the Sambirano ecoregion with 'rainforest'. We categorised species into the regions with the R package speciesgeocodeR 2. 0-10  (Töpel et al. , 2016)  and performed ancestral range reconstruction on the dated phylogeny pruned from outgroups with the R package Bio-GeoBEARS 1.
1. 2 (Matzke, 2018)  using the R codes provided on the PhyloWiki webpage (https://phylo. wikidot. com/biogeobears) and from  Antonelli et al. (2018). We tested three biogeographic models that are implemented in the package within a maximum likelihood framework to detect biases amongst tested methods: Dispersal-Extinction-Cladogenesis (DEC)  (Ree and Smith, 2008) , a likelihood interpretation of parsimony DIVA model ('DIVALIKE';  Ronquist, 1997) , and a maximum likelihood interpretation of the Bayesian-framework model implemented in the software BayArea ('BAYAREALIKE';  Landis et al. , 2013). While all three models use two free parameters of anagenetic change of range evolution (dispersal and extirpation rates), they differ in how they model cladogenetic inheritance of the ancestral ranges  (Matzke, 2018).
However, the Icong index indicated a statistically significant overall topological congruence between all phylogenies (p-value < 0. 01; 2. 86 mitochondrial vs nuclear; 3. 35 mitochondrial vs full; 3. 63 nuclear vs full), suggesting that both phylogenetic signals of mitochondrial and nuclear markers contributed to the topology of the full dataset. The ML phylogeny inferred with RAxML on the full concatenated dataset (Supplementary Fig. S3 ) showed differences from the MrBayes tree in some 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. Grandidierina in clade B), although most relationships were congruent between the two phylogenies. However, differently from the Bayesian inference, this tree had overall lower statistical (bootstrap) support (e. g. , most inter-generic relationships in Clade B received no statistical support). Similarly, the species tree inferred with ASTRAL showed some topological differences with the MrBayes phylogeny inferred from the full multi-locus concatenated dataset (Supplementary Fig. S4 ) and recovered limited statistical support: most inter-generic relationships and the monophyly of some genera (e. g. , Flexiseps, Madascincus, Pygomeles) were not statistically supported. The normalised quartet score was 0. 78 suggesting moderate concordance between gene trees.
Dated phylogeny The concatenated alignment of the reduced dataset contained 11,147 bp, and the best-fitting partitioning scheme included 10 partitions (Supplementary Table  S17 ). The 95% confidence intervals of divergence time estimates were not elongated and overlapping (Supplementary Fig. S5 ; Fig. 3A ). The Scincoidea crown node was dated in the Late Jurassic (160. 29 Ma) (Supplementary Fig. S5 ), while the diversification of the Scincidae family was inferred to be 111. 65 Ma (142. 83-79. 47 Ma 95% HPDI) at the end of the Early Cretaceous. The crown node of the Scincinae subfamily was dated at 85. 42 Ma (112. 79-59. 39 Ma 95% HPDI) in the Late Cretaceous, which also represents the stem node of the endemic Malagasy group relative to the considered outgroups. The crown node of Malagasy Scincinae was estimated at the boundary between Paleocene and Eocene epochs (55.
14 Ma; 77. 82-35. 82 Ma 95% HPDI). The same node was dated at 49. 8 Ma (72. 67-26. 16 Ma 95% HPDI), at the beginning of the Eocene, in the maximum clade credibility tree inferred from the complete dataset calibrated with secondary information (Fig. 3A ). This estimate also represents the divergence of Pseudoacontias (clade C) from the rest of the group. The two other main clades diverged slightly after  (44. 34 Ma;  Ma 95% HPDI). The inter-generic diversification in clade A was complete by the mid-Oligocene (ca. 29 Ma), whereas the genera of clade B diverged more recently, between the mid-Oligocene and throughout the Miocene (ca. 28-13 Ma) (Fig. 3A ). The genera with the most recent crown nodes were Amphiglossus (14. 7 Ma; 23. 24-6. 92 Ma 95% HPDI) and Voeltzkowia (13. 79 Ma; 21. 61-6. 62 Ma 95% HPDI).
61; speciation rate at present = 0. 05 lineages/myr; exponential rate variation = 0. 03) (Supplementary Table  S7 ). The analyses performed with the DDD 4. 4 package suggested diversity-dependent diversification in the group. The model DDLS (linear dependence of speciation rate with diversity carrying capacity and null extinction rate) fitted the data the best (AICc weight = 0. 37) with parameters estimates of speciation and carrying capacity of 0. 13 lineages/myr and 107, respectively (Supplementary Table  S8 ). The bootstrap likelihood ratio test  (Etienne et al. , 2016)  confirmed the prevalence of a diversity-dependence against a constant-rates model (Likelihood ratio = 3. 69; p-value = 0. 003; power of the test = 0. 93). Amongst the tested models from  Morlon et al. (2010) , model 2 (saturated diversity with the turnover rate declining through time) provided the best-fit (AICc weight = 0.