dataesr/scientific-paragraphs-categorization huggingface dataset:
Squamata Phylogenetic Reconstruction To reconstruct a squamate time-calibrated phylogeny we targeted up to 15 genes (six mitochondrial and nine nuclear) and included 285 species in total. Phylogenetic relationships and divergence times were estimated with Bayesian Inference (BI) using BEAST2 v. 2. 6. 4  (Bouckaert et al. 2019). We calibrated the phylogeny with 13 calibration points (Table  S2 ), avoiding fixing nodes that were of potential biogeographic interest for our analyses. We also constrained higher-level clades to match recent and supported squamate topology  (Streicher and Wiens 2017; Šmíd et al. 2021; Tejero-Cicuéndez et al. 2022). For calibration points and prior configurations for the calibration analyses as well as composition, length, partitions, models, and run-specifications for the dataset used for phylogenetic analysessee extended methods and Tables  S2-S3.
ddRADseq  Protocol, Sequencing, Data Processing, and Dataset Building  We generated double digest restriction site-associated DNA (ddRAD) libraries with the protocol implemented in  Peterson et al. (2012)  using a combination of rare and common restriction enzymes (Sbf1 and Msp1, respectively) and sequencing 75 base pairs (bp) single-end reads for a total of 661 individuals. Illumina raw reads were processed using iPyRAD v. 0. 9. 78  (Eaton and Overcast). We filtered out low quality sites (Phred score < 33), reads (≥ 3 missing sites) and sequences (<10 reads, > 3 undetermined or heterozygous sites, or > 2 haplotypes). Filtered reads were clustered and aligned using an 89% within-and between-sample clustering thresholds. The minimum number of samples per locus was left by default  (> 4)  to retrieve the maximum number of loci possible for post-processing filtering.
We further filtered our datasets with an R script provided in Burriel-  Carranza et al. (2022b)  to identify and remove individuals and loci with high levels of missingness (Table  S4 ). We generated three dataset types: i) Concatenated loci files (c_loci) were generated with ipyrad after removing all individuals that did not pass the post-processing filtering steps and retaining only loci that were present in at least 60% of all specimens; ii) Concatenated SNPs (c_SNPs), which were obtained by only selecting the SNPs within each locus; and iii) Unlinked SNPs (uSNPs), where we selected one SNP per locus. See extended methods for further details.
Species Discovery We generated species delimitation hypotheses by following a two-step workflow. We first reconstructed phylogenomic relationships within each colonization event through maximum likelihood (ML) and Bayesian Inference (BI) methods to determine well supported monophyletic groups (see extended methods for details on the parameters used, and Table  S5  datasets  [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44]. Afterwards, we conducted ADMIXTURE v. 1. 3. 0 (Alexander et al. 2009; Alexander and Lange 2011)  analyses to evaluate a range of possible ancestral populations from a minimum number of populations of K=1 to a maximum number that ranged between K=8-20 (Table  S5 , datasets. We assembled different species delimitation models (SDM) by splitting each dataset into nonadmixed clades and further lumping monophyletic groups.
Species Delimitation and Time-Calibrated Species Trees (SNAPP) We used the ML and BI phylogenies together with the population ancestry from ADMIXTURE to design and evaluate up to 49 different configurations of species hypotheses with Bayes Factor Delimitation (BFD* with genomic data;  Leaché et al. 2014; Table S5, datasets 46-57). For each SDM we estimated a species tree with SNAPP v. 1. 5. 2 (Bryant et al. 2012 ) and conducted a path sampling algorithm to calculate its marginal likelihood. We then ranked each marginal likelihood and used Bayes factors to determine the best SDM for each independent colonization event (Table  S6 ). SNAPP parameters and priors were specific to each dataset (see extended methods). Path sampling analyses were run for 20 steps with the following parameters: 500,000 MCMC generations sampling every 1,000, with an alpha of 0.
3, 10% burn-in and a preburn-in of 50,000 generations. After identifying the best SDM for each dataset, we reconstructed a time-calibrated species tree with SNAPP dating the deepest node in the phylogeny as suggested by  Stange et al. (2018). We used the 'snapp_prep. rb' ruby script from (https://github. com/mmatschiner/tutorials) to prepare the SNAPP input file. Calibration dates (means and standard deviations) were extracted from the Squamata phylogenetic reconstruction (Table  S7 ) and were set to a normal distribution.
Estimation of genetic connectivity and barriers to gene flow We implemented Fast and Flexible Estimation of Effective Migration Surfaces (FEEMS;  Marcus et al. 2021)  to determine connectivity corridors and barriers to gene flow. We generated a total of 14 datasets of uSNPs composed by non-sympatric monophyletic lineages (see Table  S5  datasets 58-71) and a grid spacing ≈10 Km between centroids. All non-gecko groups were excluded from the analysis due to low spatial coverage. We selected a specific smoothing regularization parameter  (λ)  for each dataset, selecting the lowest, leave-one-out cross-validation value out of a range of λ between 1e -6 -100 (see Table  S5 ), and exported and visualized the resulting grid in shapefile format. Afterwards, we integrated all 14 FEEMS analyses into a single result to visualize common barriers to gene flow of the endemic reptile taxa of the Hajar Mountains.
Biogeographic Reconstruction Ancestral presence in the Hajar Mountains. -We conducted an ancestral state reconstruction of mountain occupancy to infer first colonization events and the age of colonization using the function 'make. simmap' within the R package 'phytools'  (Revell 2012)  in the all Squamate multilocus phylogenetic tree. We assigned the state Hajars/No Hajars to each tip based on current species distribution, selected the most likely model under AIC criteria (all-rates-different; ARD) and ran 1,000 simulations. This also allowed us to differentiate between first colonization followed by within-mountain speciation and independent colonization events within each mountain genus. Colonization events were considered when a node was present in the Hajar Mountains (more than 50% probability of its state as Hajars), and its parental node was No Hajars.
We incorporated uncertainty to the age of colonization by considering a time range defined as the branch length between the colonization node and its parental No Hajars node, plus the 95% highest posterior density (HPD) interval of both nodes. We trimmed the ancestral state reconstruction results to the deepest node of each mountain genus with the exception of Hemidactylus, for which we only considered the 'arid clade', and Omanosaura, for which we added members of two closely related Eremiadini genera, Acanthodactylus and Mesalina, since the whole genus is endemic to the mountain range. Within-mountain biogeographic reconstructions. -We used the R package BioGeoBEARS  (Matzke 2013)  to reconstruct ancestral ranges for each independent colonization event using the SNAPP time-calibrated species trees (see above).
2020 ) and packages within tidyverse  (Wickham et al. 2019)  for data manipulation and visualization purposes. We also explored the biogeographic patterns and ancestral elevations at the specimen level by conducting ancestral state reconstructions for each independent colonization event with the function 'make. simmap' within the R package 'phytools'  (Revell 2012). This analysis was implemented in each BI reconstruction. Phylogeographic traits were established according to the three discrete topographic discontinuities of the Hajar Mountains described above, and elevation traits were discretized into three categories: lowland (<300 m), montane (300-1,500 m), and high mountain (>1,500 m). Then, we selected the most likely model (ER, SYM or ARD) under AIC criteria (Table  S8 ) and ran 1,000 simulations.
Paleoclimate We further investigated the dependence of all biogeographic events together, and vicariance events only, on global mean surface temperature estimates  (Hansen et al. 2013)  for the last 24 and 21 my respectively (ages are defined by the first biogeographic event and vicariant event recovered in our data). We fitted a GLS model over 10,000 permutations using a permutation-based method within the R package 'RRPP' v. 0. 4. 0. (Collyer and Adams, 2019). To account for temporal autocorrelation in the GLS, we built matrices of covariance incorporating the decay of temporal autocorrelation from an autoregressive model (AR1) within the R package 'nlme' v. 3. 1-162  (Pinheiro et al. 2023).
This resulted in an increase in the average reads per sample from 2. 4 x 10 6 to 2. 7 x 10 6. Reads retained per sample after sample trimming ranged between 10. 5 x 10 6 and 1. 9 x 10 5 with a median of 2. 4 x 10 6 reads (Table  S4 ). Datasets generated for lineage discovery, population structure, and lineage delimitation analyses, presented 13. 14% missing data on average, the number of loci ranged between 18,982 and 491, and contained between 173 and 7 samples. For detailed summary statistics of each dataset, refer to supplementary Table  S5.
First biogeographic events started accumulating around 23 mya, but the overall biogeographic incidence is scarce until 5 mya when most of the events occurred (Fig. 5 ). When comparing the observed incidence of all biogeographic events through time with the model expectations (which reflect the dynamics of species accumulation alone), we did not find any clear deviances (Fig. 5c ). This means that the reconstructed biogeographic history falls within the random distribution generated with the 1,000 simulated biogeographic histories and that a large part of the observed reconstruction might be explained by the shape of the phylogeny, the number of lineages through time, and the best-fit biogeographic model (DIVALIKE in all datasets). However, when looking into each type of event separately (Fig.
DISCUSSION The challenges generally associated to the study of mountain diversity are even more pronounced in arid regions, which, despite being commonly regarded as bare and homogeneous, actually harbor almost one quarter of all continental vertebrate diversity  (Brito et al. 2014). In this study, we have assembled an extensive and unprecedented genomic dataset for the complete endemic reptile fauna of a mountain range from one of the major arid areas in the world. This has allowed us to study the reptile diversity of the region at multiple scales of biodiversity organization, obtaining from new highly-detailed insights (i. e. , species boundaries) to a broad perspective understanding of biodiversity patterns in the region (i. e. , general biogeographic patterns across the entire mountain range).
Biogeographic assessment By integrating 14 different datasets within the analyses of estimation of effective migration surfaces, we detected common dispersal barriers across the mountain range. This revealed that diversity is highly structured (Fig. 3 ) and that the main topographic units of the Hajar Mountains are shaping its reptile communities (Fig. 3o ). The most elevated regions of the mountain range (Jabal Shams in Central Hajars and Saih Hatat in the Eastern Hajars; Figs. 1) are acting as geographic barriers to gene flow while lower regions may act as corridors to lowland dwellers (Fig. 3 ). These topographic discontinuities seem to be influencing diversification by increasing allopatric speciation after species range expansions. All putative species but three are currently endemic or mostly occurring in a specific mountain block, and dispersal events between blocks are scarce (Fig.
Taxon sampling To infer the number of colonization events in the Hajar Mountains we assembled a dataset following the most updated reptile taxonomy (Uetz et al. 2022) including candidate species in the process of being described and representing a total of 284 squamate species and one outgroup (Table  S1 ). This dataset contained all available species from the genera with species endemic to the Hajar Mountains, at least one representative of each Squamata family, and the only extant species of the order Rhynchocephalia (the tuatara, Sphenodon punctatus) as outgroup. We added several key species to set 13 calibration nodes in the phylogenetic tree (Table  S2 ). The eight genera that contain endemic species to the mountain range are divided in five gekkotan genera: Asaccus (seven endemic species), Hemidactylus (two endemic species), Pristurus (seven endemic species, four of them are undescribed), Ptyodactylus (two endemic species) and Trachydactylus (one endemic species); the lacertoid genus Omanosaura (two endemic species); the iguanian genus Pseudotrapelus (one endemic species); and the snake genus Echis (one endemic species).
We sampled all the available species (i. e. not only the endemics to the Hajar Mountains) from those genera except for Hemidactylus, for which we considered the "arid clade" only  (Carranza and Arnold 2006). To investigate the diversity within each of the endemic reptile species of the Hajar Mountains we selected specimens spanning throughout the whole mountain range and including, when possible, individuals from all known genetic lineages of each species (Table  S3 ;  Carranza and Arnold 2012; Badiane et al. 2014; Tamar et al. 2016; de Pous et al. 2016; Garcia-Porta et al. 2017; Simó-Riudalbas et al. 2017; Tamar et al. 2019a Tamar et al. , 2019b;; Fattahi et al. 2020). For each monospecific group, we complemented the dataset with its closest described relative to root the phylogenetic trees.
We retrieved 53 ddRADseq samples of Trachydactylus from Burriel-  Carranza et al. (2023)  and the total intraspecific dataset included 661 specimens from 27 described species (19 endemic species and 8 outgroups). All the specimens were collected between 2006 and 2017, were preserved in 99% ethanol and stored at -20ºC until library preparation (see below "ddRADseq protocol").
Squamata phylogenetic reconstruction The final dataset to explore the colonization of the Hajar Mountains included a total of 285 species (284 squamates and the tuatara, see above We used a supermatrix approach, concatenating all the genetic markers for the 285 samples. A total of 1,804 sequences were aligned using two different protocols for the ribosomal and for the coding genes, respectively. Mitochondrial non-transcribed 12S and 16S ribosomal genes were aligned by performing multiple sequence alignments using MUSCLE implemented in Geneious Prime 2020. 2. 5. Poorly aligned positions of these two genes were eliminated with G-blocks  (Castresana 2000)  using low stringency options  (Talavera and Castresana 2007). Coding   S1 ). We reconstructed the phylogenetic relationships and estimated the divergence times with a Bayesian inference analysis using BEAST2 v.
2. 6. 4  (Bouckaert et al. 2019). We calibrated the phylogeny with 13 calibration points (Table  S2 ), avoiding fixing nodes that were of potential biogeographic interest for our analyses. We also constrained higher-level clades to match recent and supported Squamate topology  (Streicher and Wiens 2017; Šmíd et al. 2021; Tejero-Cicuéndez et al. 2022). We conducted four individual runs of 5x10 7 generations, sampling every 10,000 generations with a Yule process tree prior. The dataset contained 10 different partition schemes obtained from Partition Finder, and a discretized Gamma distribution with four categories was set for all partitions while independent substitution models were individually selected for each of them (Table  S3 ). Convergence between runs and stationarity of each chain was checked with Tracer v.
Data processing We processed raw Illumina reads using iPyRAD v. 0. 9. 78  (Eaton and Overcast 2020). We demultiplexed samples using their unique barcoded adapter sequences. Sites with Phred score < 33 (scores with quality under 99%), and reads with ≥ 3 missing sites were discarded. Within the iPyRAD pipeline, filtered reads were clustered and aligned using an 89% within-sample clustering threshold. As an additional filtering step, we filtered out consensus sequences that had low coverage (< 10 reads), excessive undetermined or heterozygous sites (> 3) or too many haplotypes (> 2). Consensus sequences were clustered between samples using an 89% Afterward, we removed all non-biallelic SNPs with Plink2  (Chang et al. 2015) , applied a minor allele frequency (maf) filter to prevent sequencing errors (maf < 0.
05) and removed monomorphic sites using Radiator  (Gosselin et al. 2017). Finally, we selected one SNP per locus to create a dataset of putatively unlinked SNPs (uSNPs) to account for Linkage Disequilibrium. If there was more than one SNP per locus, the highest-depth SNP was selected. In case there was more than one SNP with the highest depth value, one of those was randomly chosen. For specific information on applied filters, number of samples used, loci and SNPs recovered for each dataset refer to Table  S4. All ddRADseq analyses were performed using one of the following dataset types. For Maximum Likelihood (ML) and Bayesian Inference (BI) reconstructions, we used concatenated loci files (c_loci) generated with ipyrad after removing all individuals that did not pass the previously explained filters and retaining only loci that were at 60% of all individuals (Table  S4 ).
Species Discovery Phylogenetic reconstruction. -Using the c_loci alignments (Table  S4 , datasets 25-34) of each independent colonization event, we performed phylogenetic reconstructions through a maximum likelihood (ML) inference with RAxML-ng v. 1. 0. 2  (Kozlov et al. 2019) , with a GTR+GAMMA model, a total of 100 starting trees (50 random and 50 parsimony) and 1,000 bootstrap replicates. We also estimated Bayesian Inference (BI) with BEAST2 v. 2. 6. 4  (Bouckaert et al. 2019 ). To reduce computational load, for each independent colonization we used two independent sets of 600 randomly picked loci present in at least 60% of samples (Table  S4 , datasets 35-44). We calibrated each phylogeny by extracting the date of its deepest node from the squamate phylogenetic tree previously inferred and applying a normal distribution (See Table  S7  for specifications on the calibration nodes).
We selected a GTR model with 4 gamma categories, base frequencies were estimated, and a relaxed clock LogNormal was used. We conducted two individual runs of 10 8 generations sampling every 10,000 generations. Convergence was checked with Tracer v. 1. 7  (Rambaut et al. 2018 ) and a burnin of 40% was applied. Convergence between sets of loci was achieved in all cases with the exception of Ptyodactylus and Trachydactylus for which we had to use all available loci (17,588 and 5,115 loci respectively) to achieve convergence between runs. Admixture. -Using the uSNPs datasets of each described species (Table  S4 ; datasets 1-24), we inferred the population ancestry of each individual with ADMIXTURE v. 1. 3. 0  (Alexander et al. 2009; Alexander and Lange 2011). This program models the probability of the observed genotypes estimating simultaneously population allele frequencies along with ancestry proportions assuming Hardy-Weinberg equilibrium within populations.
We analyzed the uSNP datasets using the script 'admixture-wrapper' (https://github. com/dportik/admixturewrapper) to evaluate a range of possible ancestral populations from a minimum number of populations of K=1 to a maximum number of populations that ranged between K=8-20 (depending on the maximum number of individuals in each group). For each K, we generated 15 replicates and 15 fold cross-validations to determine the most probable K. We visualized the best K in a geographic map and generated first species delimitation hypotheses assigning each individual to the ancestral population from which the largest proportion of its genome is derived.
Species validation and time-calibrated species trees (SNAPP) We used the ML and BI phylogenies together with the population ancestry from For all analyses, mutation rates (u & v) were fixed to 1. The Yule prior (λ) representing the speciation rate was set to a gamma distribution and while alpha was set to 2, beta was estimated by calculating the expected tree height on the c_loci datasets (maximum observed divergence between any pair of taxa divided by 2). We then used 'pyule' (https://github. com/joaks1/pyule) to determine the mean value of lambda and calculated beta accordingly (λ = α × β). Theta prior (θ) was also set to a gamma distribution and the mean value of θ was estimated by averaging all genetic distances within each population. Each BFD* analysis was also conducted with default λ & θ priors and the most conservative result was selected.
Path sampling analyses were run for 20 steps with the following parameters: 500,000 MCMC generations sampling every 1,000, with an alpha of 0. 3, 10% burnin and a preburnin of 50,000. Stationarity of all runs was checked and each step was run until ESS >= 200. After identifying the best species delimitation model for each dataset, we reconstructed timecalibrated species trees with SNAPP dating the deepest node in the phylogeny as suggested by  Stange et al. (2018)  with a normal distribution. Calibration dates were extracted from the Squamata phylogenetic reconstruction (See Table  S7 ). To prepare the SNAPP input file we used the 'snapp_prep. rb' ruby script from (https://github. com/mmatschiner/tutorials). This script allows for incorporating and calculating the node ages while reducing the running time required for SNAPP by linking all population sizes.
Fast and flexible estimation of effective migration surfaces (FEEMS) To visualize spatial genetic structure and try to detect common genetic barriers in our data, we implemented FEEMS  (Marcus et al. 2021). This program allows the user to visualize spatially heterogeneous isolation-by-distance on a geographic map detecting effective migration surfaces and regions with less gene flow than expected. This approach can be useful to unravel geographic features that enhance or diminish the levels of genetic differentiation between populations. Moreover, in this study we integrate for the first time several FEEMS analyses into a single result to visualize common barriers to gene flow of the endemic reptile taxa of the Hajar Mountains. We constructed a discrete global grid system with 'dgconstruct' in the R package dggridR (Barnes et al.
2017) with a spacing resolution of 9 (about 10 km spacing between grid centroids) and clipped it with a self-generated perimeter of the mountain range (when samples spanned throughout the whole mountain range), or with a block specific range (when specimens were only present in one or two mountain blocks). We generated a total of 14 datasets of uSNPs composed of non-sympatric monophyletic clades (see Table  S5 ; datasets 58-71). Four datasets contained individuals spanning throughout all the Hajar Mountains and the rest were constrained to the region where the analysed lineages resided. All non-gecko groups were excluded from the analyses due to small sample size. FEEMS uses a smoothing regularization parameter λ to tune each edge weight on the graph. Large values of λ (e. g. λ = 100) promote FEEMS to fit a model with most edges nearing the mean value, while small values of λ (e.
g. λ = 0. 00001) promote overfitting of the data  (Marcus et al. 2021). To select the proper λ for each dataset, we selected the lowest, leave-one-out crossvalidation value out of a range of lambdas between 1e -6 and 100. We exported and visualized the resulting grid in shapefile format. With a custom R script (https://github. com/BernatBurriel/Hajar-Mountains/), we merged all 14 analyses into a single file to identify general patterns across the entire range of the Hajar Mountains. Since FEEMS's average migration weight equals 0 in all analyses, we were able to average all overlapping edges between the 14 analyses without having to apply any data transformation. We interpolated the results into a 1 km 2 resolution raster for better visualization.