The enigmatic Leiosaurae clade: Phylogeography, species delimitation, phylogeny and historical biogeography of its southernmost species

Journal Pre-proofs The enigmatic Leiosaurae clade: Phylogeography, species delimitation, phylogeny and historical biogeography of its southernmost species Martin M. Femenias, Luciano J. Avila, Jack W. Sites Jr., Mariana Morando PII: DOI: Reference:

S1055-7903(19)30350-1 https://doi.org/10.1016/j.ympev.2019.106725 YMPEV 106725

To appear in:

Molecular Phylogenetics and Evolution

Received Date: Revised Date: Accepted Date:

7 June 2019 24 December 2019 24 December 2019

Please cite this article as: Femenias, M.M., Avila, L.J., Sites, J.W. Jr., Morando, M., The enigmatic Leiosaurae clade: Phylogeography, species delimitation, phylogeny and historical biogeography of its southernmost species, Molecular Phylogenetics and Evolution (2019), doi: https://doi.org/10.1016/j.ympev.2019.106725

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier Inc.

The enigmatic Leiosaurae clade: phylogeography, species delimitation, phylogeny and historical biogeography of its southernmost species. Martin M. Femenias a, Luciano J. Avila a, Jack W. Sites, Jr. b, Mariana Morando a, c aInstituto

Patagónico para el Estudio de los Ecosistemas Continentales, Consejo Nacional de

Investigaciones Científicas y Técnicas (IPEEC-CONICET), Boulevard Almirante Brown 2915, U9120ACD, Puerto Madryn, Chubut, Argentina. bDepartment

of Biology and M.L. Bean Life Science Museum, Brigham Young University (BYU), Provo, UT

84602, USA; current address: Department of Biology, Austin Peay State University, Clarksville, Tennessee 37044 cUniversidad

Nacional de la Patagonia San Juan Bosco, Sede Puerto Madryn, Boulevard Almirante Brown

3700, U9120ACD, Puerto Madryn, Chubut, Argentina. Keywords: Matuasto; Leiosauridae; Plio-Pleistocene; Speciation; Conservation; Patagonia

Abstract The clade Leiosaurae is composed of poorly-known species endemic to the southern region of South America. The difficulties of finding these lizards in the field, and their highly conserved morphology, have limited our taxonomic knowledge and understanding of their evolutionary histories. Here, we use data collected over 9 years to study the phylogenetic history, genetic diversity, and biogeographic history of almost all the southernmost species of Leiosaurae (except P. nigroigulus), including: Leiosaurus bellii, Diplolaemus darwinii, D. bibronii, D. sexcinctus and D. leopardinus. We use a fragment of the mitochondrial cytochrome-b gene to resolve general phylogeographic patterns, and add another mitochondrial gene and eight nuclear genes to perform species delimitation and phylogenetic analyses associated with divergence times. We found evidence for three putative new species-level taxa within L. bellii and five within Diplolaemus species, indicating high levels of geographic structure. We used a time-calibrated phylogeny to estimate ranges of ancestral distributions and to generate new hypotheses about their historical biogeography.

1. Introduction Knowledge of biological diversity is essential for ecological studies of regional communities and ecosystems, species richness, and for conservation of natural communities and ecosystems (Bickford et al., 2007; Camargo and Sites, 2013; De Queiroz, 2007; Sites and Marshall, 2003). The southernmost arid region of South America is recognized for its high degree of endemism in different vertebrate groups (Lamoreux et al., 2006). The geological and climatic fluctuations that occurred since the Miocene and especially during the Quaternary period greatly contributed to the origin and evolution of the region's biological diversity. Several examples of how geological events such as glaciations, the formation of paleobasins, and marine transgressions, have acted on ancestral populations of current species are well documented for several clades of animals and plants in Patagonia (Breitman et al., 2012; Cosacov et al., 2013; Fontanella et al., 2012; Sérsic et al., 2011). Although these works mark a strong baseline for ecological and biogeographical studies, the diversity of other clades that are underrepresented in the scientific literature, remains to be explored with a phylogeographic approach. This lack of knowledge not only makes it impossible to apply appropriate conservation actions, but also gives us a biased understanding of the processes that underlie the evolution of species. Leiosaurae is among the southernmost distributed lizard clades in the world and its species are commonly known as “matuastos”. It includes 18 species, almost all of which are difficult to find in the field; thus, they are underrepresented in herpetological collections, and one of the least-studied Patagonian lizard clades. Furthermore, the highly conserved morphologies of these lizards limit distinguishable variation among species, which in turn limits our understanding of the evolutionary history of the clade. Some phylogenetic hypotheses have been proposed based on molecular and morphological data (Abdala et al., 2009; Cei et al., 2003; Frost et al., 2001; Pyron et al., 2013; Reeder et al., 2015; Schulte et al., 2003; Townsend et al., 2011; Victoriano et al., 2010; Wiens et al., 2012; Zheng and Wiens, 2016), but the limited number of taxa sampled in these studies compromises phylogenetic inference. A molecular phylogeny based on 10 markers for 14 species, including all described species of the genera Diplolaemus (n = 4) and Leiosaurus (n = 4),

and six of ten species of Pristidactylus, was presented by Morando et al. (2015). This study advanced understanding of the evolutionary history of this clade; however, some relationships are poorly supported, species limits are poorly resolved for some taxa, and diversification patterns remain to be evaluated in a biogeographical context. The southernmost species of the Leiosaurae clade, including Leiosaurus bellii, Pristidactylus nigroiugulus, and all Diplolaemus species are widely distributed in regions characterized by high endemism in other species groups (Andrade-Díaz et al., 2017; Corbalán et al., 2011; Morando et al., 2007; Roig et al., 2009). However, the species diversity within some of these groups remains unclear due to between-species conserved morphologies. For example, morphometric and meristic analyses in Diplolaemus reveal that morphology is in some cases so conserved among species that other data are required to discriminate species (Cei, 1986; Scolaro et al., 2003; Victoriano et al., 2010). Thus, exploring the diversity of these groups with molecular approaches is essential to assess the underlying species diversity. In this work, we employed molecular approaches to delineate population identities and to infer evolutionary histories of the southernmost distributed Leiosaurae species. Under the hypothesis that species diversity in this clade is greater than that inferred so far with morphological data, the objectives of this study are to: (1) infer phylogeographic histories for the southernmost distributed species of the Leiosaurae clade (L. bellii and D. bibronii, D. darwinii, D. sexcinctus, D. leopardinus; except P. nigroigulus; (2) clarify the identity of geographically isolated populations; (3) provide a time-calibrated multilocus phylogenetic hypothesis of this clade; and (4) elucidate its biogeographic history. Our molecular phylogeny reveals geographically concordant clades that, based on species delimitation methods, can be considered to represent candidate species. We suggest hypotheses of how the geologic and climatic fluctuations that occurred in PlioPleistocene could have modeled the biogeographic histories of current species. We provide a useful baseline for further field work to fill distributional gaps, as well as future species-delimitation studies, and further extension/refinement of historical biogeographic scenarios. This and further studies will provide reliable information for designing appropriate conservation plans.

2. Materials and Methods 2.1. Taxon sampling and laboratory procedures We collected 318 individuals representing the clade Leiosaurae, from 232 localities between the years 2006 and 2015. The geographic sampling includes the entire known distributional range of Diplolaemus and Leiosaurus species, as well as the six species of Pristidactylus (P. achalensis, P. araucanus, P. fasciatus, P. nigroiugulus, P. scapulatus and P. torquatus) (Tables S1 and S2). As outgroups, we used two samples of Urostrophus vautierii from Brazil and a 12S sequence of Urostrophus gallardoi from GenBank (AF338325.1, Table S2) that are part of the Anisolepae clade (=Enyaliinae subfamily), which is strongly supported as the sister clade of Leiosaurae (Pyron et al., 2013; Reeder et al., 2015; Wiens et al., 2012). We extracted DNA from tissues with the Qiagen DNeasy Tissue Kit following the protocol provided by the manufacturer. We amplified 10 fragments, including two mitochondrial: cytochrome-b (primers of Palumbi 1996 and Whiting et al. 2004) and 12S (primer of Wiens et al. 1999); five nuclear: DmX-like protein 1 gene (DMXL-1, primers of Werneck et al. 2012), C-Moloney oncogene sarcoma (CMOS, primers of Wiens et al. 1999), Exophilin 5 gene (EXPH5, primers of Townsend et al., 2008), intron 8 and flanking exon regions of RNA binding motif protein X-linked (RBMX, primers of Gamble et al. 2011), NK-tumor recognition gene (NKTR, primers of Townsend et al. 2011); and three Anonymous Nuclear Loci (ANL; Diplolaemus: ANLD90, ANL-D92, ANL-D77; primers of Morando et al. 2015). All nuclear and mitochondrial genes were amplified following PCR protocols of Morando et al. (2003) and Morando et al. (2004), and ANL protocols followed those of Morando et al. (2015). Double-stranded PCR amplified products were checked by electrophoresis on a 2% agarose gel, purified using a Multi-Screen PCR (l) 96 (Millipore Corp., Billerica, MA), and directly sequenced using the BigDye Terminator v 3.1 Cycle Sequencing Ready Reaction (Applied Biosystems, Foster City, CA). All cycle sequencing reactions were purified using Sephadex G-50 Fine (GE Healthcare) and MultiScreen HV plates (Millipore Corp.). Samples were analyzed with an ABI3730xl DNA Analyzer at the BYU DNA

Sequencing Center. We used Sequencher v4.10. (Gene Codes Corporation Inc.™ 2007) to edit and to align all genes. We translated protein-coding regions to amino acids to check for stop codons, and alignments were checked by eye. Alignments are available at TreeBASE (http://purl.org/phylo/treebase/phylows/study/TB2:S20311) We used MAFFT v7.266 (Katoh and Standley, 2013) for alignment, following the G-INS-i algorithm as the alignment strategy. The partition schemes and the nucleotide substitution models were selected using PartitionFinder v1.1.1 (Lanfear et al., 2012), and we tested all schemes and 56 possible substitution models for all genes with Bayesian (BIC) and Akaike information criterion (AICc) methods. 2.2. Phylogeographic and phylogenetic analyses 2.2.1. Gene trees We used 26 previously published cytochrome-b sequences and added 294 new sequences, giving a total of 320 individuals for this study (Table S1). Gene trees were estimated by Bayesian (BI) and Maximum Likelihood (ML) methods using unique haplotypes resolved with DnaSP v5 (Librado and Rozas, 2009). We performed Bayesian inference with BEAST2 v2.4.3 (Bouckaert et al., 2014), assigning best-fit partition schemes and substitution models (Table 1). Four hot chains were run for 1x108 generations of Monte Carlo Markov Chains (MCMC), and sampled at intervals of 1000 generations. We discarded 25% as burn-in and used the equilibrium samples to generate the consensus tree. We used posterior probability values greater than 0.95 to represent nodes with high support (Huelsenbeck and Ronquist, 2001). We obtained the maximum likelihood tree with RAxML v8.2.7 (Stamatakis, 2014) using a GTR model with gamma distribution as suggested by the best-fit partition scheme. A bootstrap analysis was used with 1000 quick searches, and considered values greater than 70% as a strong nodal support (Hillis and Bull, 1993). 2.2.2. Genetic diversity We constructed statistical parsimony networks (Templeton et al., 1992) with TCS version 1.21 (Clement et al., 2000), and assigned connection limits of 95%. The number of haplotypes (h), polymorphic sites (S) and

haplotype (Hd) and nucleotide (π) diversities were calculated for the cytochrome-b gene with DnaSP v5 (Librado and Rozas, 2009). We calculated pairwise between-clade distances with Arlequin 3.5 (Excoffier and Lischer, 2010). 2.2.3. Demographic history To evaluate deviations from the neutral model we estimated Fu´s Fs (Fu, 1997), Tajima´s D (Tajima, 1983) and R2 statistics (Ramos-Onsins and Rozas, 2002) with DnaSP v5 (Librado and Rozas, 2009) based on cytochrome b sequences. The same program was used to obtain the significance of each statistic through coalescent simulations for segregating sites (Hudson, 1990), using 5000 replicates. All of these estimates were made at the species level, and in some cases for haploclades within species. Signs of demographic or range expansion were assessed through pairwise distances with DnaSP v5 (Librado and Rozas, 2009) and Extended Bayesian Skyline Plots (EBSP, Drummond et al., 2005) using BEAST 2 v2.4 (Bouckaert et al., 2014), and Tracer v1.6 (Rambaut et al., 2014). For the reconstruction of EBSP, we used two mitochondrial and eight nuclear loci for D. darwinii, D. bibronii, D. sexcinctus sensu stricto, and D. lineage 2; and used two mitochondrial and seven nuclear loci for L. bellii. We could not perform these analyses on any remaining Diplolaemus due to small sample sizes, which can affect the estimates (Drummond et al., 2005; Ho and Shapiro, 2011). The substitution models and partition schemes for each dataset were selected with PartitionFinder v1.1.1 (Lanfear et al., 2012), and we applied a strict molecular clock for each gene with substitution rates obtained from Morando et al. (2015). We ran 1x108 generations of MCMC sampled at intervals of 1x106 generations, and registered trees every 5000 generations, excluding the first 10% as burn-in. There is no published information on the generation time of Leiosaurae species; based on our field observations, we assume that the generation time should be between 5 and 10 years. Therefore, we used a generational time of 10 years to rescale our Y-axis to the effective population size. We used R version 3.4.1 (R Core Team, 2017) to graph all results. 2.3. Species delimitation

We used two analyses based on a multispecies coalescent model (Heled and Drummond, 2010) to delimit putative species, with the programs BP&P v4 (Rannala and Yang, 2013; Yang and Rannala, 2010) and STACEY (Jones, 2017; Jones et al., 2014). We used two mitochondrial loci (cytochrome b and 12S) and eight nuclear loci (DMXL-1, CMOS, EXPH5, RBMX, NKTR, ANL-D77, ANL-D90, ANL-D92) for a total of 114 individuals assigned to 24 groups (including Urostrophus species as a single outgroup, Table S2). The individuals used for these analyses, their locations and genes for each individual are summarized in Table S2. Table S3 indicates the percentage of missing data. In the BP&P analysis we implemented species delimitation and species-tree estimation algorithm (A11) with a prior of 1, which assigns equal probabilities to the rooted species trees (Yang, 2015). We used combinations of priors for ancestral population sizes and root ages by adjusting the prior gamma; these estimated the best model as inv gamma IGθ(3,0.002) and IGτ(3,0.0055). We implemented a second delimitation analysis using the STACEY package v1.2.2 (Jones, 2017; Jones et al., 2014), as implemented in BEAST2 v2.4 (Bouckaert et al., 2014), and with the same matrix as the BP&P analysis with 24 minimal clusters (SMC; Jones et al., 2014). For each locus we assigned the corresponding partition scheme and molecular evolution models (Table 1). Ploidy values were fit to 2.0 and 0.5 for nuclear and mitochondrial genes, respectively, and the species tree was estimated with a Yule model (Yule, 1924) with collapsed height nodes of ε = 1x10-6, and birth diff rate of 100 a priori. The collapsed weight ω was estimated with a uniform prior on [0,1]. The relative death rate was assigned with a uniform distribution between -0.5 and 0.5. We assigned lognormal priors to growth, popPriorScale and clock rates. The growth rate 'bdcGrowth' was assigned with 182.5 in a LogNormal distribution (mean Ln (182.5) = 5.2 and sd = 3). All settings were assigned following the recommendations described in the program manual. We ran the MCMC chain for 100 million generations, storing the state every 10 million generations and logging trees every 50,000 generations. We used Species Delimitation Analyser v1.8.0 (http://www.indriid.com/software.html), and assigned the same collapse height and simcutoff 1 in order to obtain posterior probabilities for grouping clusters.

2.4. Species tree Results of the species-delimitation analyses with new “putative species” (identified by a number or letter) were used for the reconstruction of the species tree with the StarBEAST2 package (Ogilvie et al., 2017), as implemented in BEAST2 v2.4 (Bouckaert et al., 2014). We used independent matrices for each locus for a total of 87 individuals (see percentages of missing data per locus in Table S3). All matrices were combined in a BEAST analysis assigning the corresponding partition schemes and the molecular substitution models (Table 1) and random local clocks (RLC), for each locus. We assigned 24 taxa a priori and fixed the ploidy of each locus with values of 2.0 and 0.5 for nuclear and mitochondrial genes, respectively. In order to obtain population sizes, we set a constant population model and selected a birth-death model that allows a nonzero rate of extinction, and a priori assigned the speciation rate to 182.5. Because substitution rates within species should not vary between biologically unlikely ranges, we assigned a range of 0.1-10 for substitution variation rates between specimens of the same species (i.e. the fastest rate cannot be 100 times larger than the slowest rate within the same species). We assigned molecular clocks a priori with means obtained in STACEY analyses, using a LogNormal distribution in real space with standard deviation equal to 1. We ran 1x108 generations of MCMC sampled at intervals of 1x106 generations, and registering trees every 5000 generations. ESS values greater than 200 and subsequent distributions of the estimated parameters and their highest posterior density intervals (HPD), were visualized in Tracer v.1.6 (Rambaut et al., 2014). The species tree was obtained with TreeAnnotator v2.4.2, again with the initial 10% of generations discarded as burn-in. The consensus tree was visualized with FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/), and DensiTree v2.0 (Bouckaert, 2010) was used to visualize the set of gene trees. 2.5. Calibration and divergence times After substitution rates were estimated, the tree was calibrated in order to estimate the divergence times of each lineage and mutation rates for each gene. A new multilocus analysis (with the same nuclear and mitochondrial dataset that was used as input in the species tree) was performed following a Yule model calibrated with an a priori birth rate of 1.0. The means of molecular clocks for each locus were assigned

from the results obtained in the previous analysis with a 1/X distribution, allowing calibration of the tree. Albino et al. (2016) recognized two new fossil specimens belonging to the genus Pristidactylus based on fragments of jaws found in the Santa Cruz Formation located in the Southern Geological Basin of Santa Cruz province. This Formation is geochronologically dated to 18-16.20 MYA based on tephrochronological (Perkins et al., 2012; Vizcaíno et al., 2012) and radiometric methods (Cuitiño et al., 2016). Pristidactylus sp. A described by Albino et al. (2016) is characterized by location of the anterior inferior alveolar foramen within the splenial bone, a characteristic shared only with P. torquatus, although it is located dorsally in the fossil. The species Pristidactylus sp. B shares characters with fossils assigned to Pristidactylus sp., from the Colhuapense Formation in Chubut province (Albino, 2008). This Formation has been dated 20-20.4 MYA (Ré et al., 2010). Two of the specimens analyzed in this study have well-preserved splenial bones, and also are characterized by location of the anterior inferior foramen within this bone. Given these characteristics, the fossil can be placed at the base of the clade Leiosaurae, before the divergence of P. torquatus, assuming that the plesiomorphic state is a jaw with more teeth (between 13 and 18), and a lower alveolar foramen located inside the splenial bone. Thus, we assigned an exponential distribution to the basal node of the Leiosaurae clade, using a displacement of 16.20 MYA as the lower bound. We assigned mean=0.9 for this distribution so that the a priori node's mean is 17.1 MYA, the midpoint of the dated range for these fossils. We ran 1x108 generations of MCMC sampled at intervals of 1x106 generations, and registering trees every 5000 generations. We discarded 10% of initial generations as burn-in. ESS values greater than 200 and estimated mutation rates were visualized in Tracer v.1.6 (Rambaut et al., 2014). The calibrated species tree was obtained with TreeAnnotator v2.4.2, discarding the initial 10% of generations as burn-in. The consensus tree was visualized with FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/), and DensiTree v2.0 (Bouckaert, 2010) was used to visualize the set of gene trees. 2.6. Historical biogeography

We performed a probabilistic inference of the ancestral species range using BioGeoBEARS (BioGeography with Bayesian - and likelihood - Evolutionary Analysis in R Scripts) (Matzke, 2014, 2013), as implemented in R version 3.4.1 (R Core Team, 2017). The models tested were the Dispersion/Extinction/Cladogenesis (DEC; Ree, 2005; Ree and Smith, 2008), Dispersion-Vicariance Analysis model (DIVA; Ronquist and Cannatella, 1997), and Bayesian Inference model (BayArea, Landis et al., 2013). In addition, these three models were further tested by including a founding effect parameter "j" (DECj, DIVAj, and BayAreaj) (Matzke, 2014). We compared the adjustment of the six models by Akaike Information Criterion, and compared significant differences between nested models by Natural Logarithm of Likelihood (LnL) in R version 3.4.1 (R Core Team, 2017). We used a calibrated ultrametric species tree obtained in StarBEAST2 (see section Phylogenetic Analysis). This tree was pruned with Mesquite v3.51 (Maddison, 2018), including all clades delimited for Diplolaemus. Biogeographically significant areas were defined to test the possible historical scenarios, following the results of previous phylogeographic works in the region, based on lizards, plants, and rodents (Sérsic et al., 2011). The delimited areas were: A, the region extending from the Deseado River (northern limit) to the south; B, the area between Deseado River (southern limit) and the Chico and Chubut rivers (northern limit); C, the area between the Chico River (in the Chubut province) and the Chubut River; D, the central plateau of Chubut and Río Negro; E, the central area of Somun Curá Plateau; F, the Limay River basin to southern Neuquén; G, the Neuquén River basin to northern Neuquén; H, the area between Colorado and Negro rivers in Neuquén province; I, the area between Colorado and Negro rivers in Rio Negro province; J, northwestern Santa Cruz province and southwestern Chubut province; K, the region between northwestern Chubut and southwestern Rio Negro; and L, the pre-cordillera region in Mendoza province. 3. Results 3.1. Taxon sampling and laboratory procedures We sequenced all the markers for all species of the genera Diplolaemus and Pristidactylus, but it was not possible to obtain sequences for the anonymous locus ANL-90 for Leiosaurus and Urostrophus. Table 1

describes all the markers, their lengths in base pairs and the best partition schemes and substitution models for each matrix. The new sequences included in this work were deposited in GenBank with accession numbers MN700660-MN700857, MN75506-MN705799, MN714928-MN715126, MN756546-MN756593 (Tables S1 and S2). 3.2. Gene tree We obtained 156 unique haplotypes of the cytochrome-b fragment used for the gene tree reconstruction shown in Figure 1 with their respective support values for Bayesian and maximum likelihood analysis. Leiosaurus bellii is inferred as a well-differentiated clade (PP = 1, BT = 93) structured into three haploclades (Lineage A, Lineage B, and Lineage C), each with high statistical support (PP = 1, BT > 70). The genus Diplolaemus is inferred as a clade with moderate support (PP = 0.96, BT < 70), within which D. darwinii is inferred with modest support (Bayesian but not ML), as the sister group to all other species. The D. darwinii clade includes two haploclades that may correspond to species lineages: a southern Lineage A: PP = 1, BT = 73, and Lineage B that is not statistically supported. Diplolaemus leopardinus was identified as a haploclade with high statistical support (PP = 1, BT = 99) in both Bayesian and ML analyses. The Bayesian analysis infers D. leopardinus as the sister taxon to (PP = 1) a paraphyletic Diplolaemus sexcinctus (sensu Victoriano et al. 2010), within which D. bibronii is nested. The ML analyses nested D. leopardinus within D. bibronii, and both of these nested within D. sexcinctus sensu Victoriano et al. (2010), but not statistically supported (Fig. S1). Six haploclades within D. sexcinctus are hypothesized to correspond to species-level lineages. We identified as D. sexcinctus sensu stricto the haploclade that includes individuals from the type locality. The other five haploclades (1-5) previously assigned to D. sexcinctus are inferred with high statistical support, but D. s. Lineage 1 and D. s. Lineage 2 are not strongly supported in the ML analyses (BT < 70). Diplolaemus bibronii was inferred as a haploclade with high statistical support in the Bayesian analysis (PP = 0.98); within it, D. bibronii Lineage A, which included five haplotypes distributed in the northernmost range of the species (northwestern of the Pampa del Castillo and Pampa de Salamanca mountain ranges), is

paraphyletic with respect to a well-supported, widely distributed southern clade, D. bibronii Lineage B (PP = 1, BT = 71; Fig.1). 3.3. Genetic diversity for cytochrome b haplotypes Haplotype diversity is high (Hd > 0.8, Table 2) for all species. At the intraspecific level, only Diplolaemus sexcinctus Lineage 4 and D. bibronii Lineage A have low haplotype diversities (Hd = 0.46 and Hd = 0.42 respectively). Indexes of nucleotide diversity are relatively low for the southernmost species D. darwinii and D. bibronii (0.01-0.02). Diplolaemus sexcinctus (sensu Victoriano et al. 2010) shows high nucleotide diversity indexes, while the six species lineages within it have lower indexes by one to two orders of magnitude (Table 2). The index of nucleotide diversity in L. bellii is almost double that of its lineages A and B, and an order of magnitude higher than Lineage C (Table 2). However, only three samples are included in Lineage C; therefore, these indexes may not represent good estimates of the real diversity. Breitman et al. (2012) estimated that uncorrected pairwise cytochrome b distances > 3% can be used for the diagnosis of candidate species within the Patagonian lizard genus Liolaemus. Although there is no similar analysis for Leiosaurae, we used this approximate uncorrected pairwise cytochrome b threshold distance as indicator of putative candidate species, here referred as lineages. All uncorrected pairwise cytochrome b distances among Diplolaemus species and hypothesized species lineages within D. sexcinctus, and Leiosaurus bellii are greater than 3%, except for the distance between D. sexcinctus sensu stricto and Lineage 5 (2.59%). The intra-group distances vary between 0.35% and 1.60% (Table 3). Maximum parsimony network analysis inferred Lineages 1, 3 and 4 of D. sexcinctus (sensu Victoriano et al. 2010) as separate networks (Fig. 2), while D. s. Lineage 2 includes a network plus four unconnected haplotypes. Haplotypes of D. sexcinctus sensu stricto were grouped into a single network that also includes five haplotypes of D. s. Lineage 5 (localities 97, 107, 111, 113, 114, located at the southwestern-most part of its distribution). The remaining D. s. Lineage 5 haplotypes (H6 and H9) were grouped into a separate network on the central Somun Curá plateau (localities 121-124, 126, 128, 143, 146). The southernmost D. leopardinus haplotypes in Mendoza (locs. 207-210) were grouped into a single network, and a singleton

haplotype (H41) is isolated in northern Mendoza province (Fig. 2, loc. 212). Most haplotypes from D. darwinii (Fig. 3, right panel) were grouped into a single network, D. darwinii Lineage A, while those corresponding to D. darwinii Lineage B are two separate haplotypes (H2, H8) from Pampa del Castillo (loc. 61). Diplolaemus bibronii Lineage A haplotypes are distributed north of Pampa del Castillo (Fig. 3, left panel, locs. 64, 71, 80-82), and resolved as two simple networks and a singleton (H16, loc. 64). Haplotypes south of Chubut province are part of a network that corresponds to D. bibronii Lineage B (Fig. 3). These haplotypes are distributed mainly in Santa Cruz province south of Pampa del Castillo, but haplotype 15 occurs in the northernmost locality of the distribution of this species, almost center of the Chubut province (loc. 86). Three networks inferred for Leiosaurus bellii correspond to the putative species lineages identified in the gene tree (Fig. 1). Leiosaurus bellii Lineages A and B grouped respectively all the haplotypes distributed to the south and north of the Negro and Neuquén rivers. In L. bellii Lineage A, most of the haplotypes are distributed mainly around the Somun Curá plateau, but four (H2, H5, H16, H36) are distributed north of the Limay River (locs. 165, 175, 177, 178). L. bellii Lineage C has three haplotypes, two distributed north of the Colorado River (locs. 211, 206) and one (Hap 35) distributed in central Neuquén province (loc. 178) between the Neuquén and Limay rivers (Fig. 4). 3.4. Demographic histories Diplolaemus bibronii Lineage B and D. darwinii showed significant and negative Tajima's D indicating a deviation from the neutral model by virtue of a non-random process such as demographic contraction/expansion based on cytochrome b locus. In addition, for D. bibronii Lineage B, the Fu's Fs was high and negative and R2 was significant, which would indicate that this deviation of the neutral model may be due to a recent population expansion. Similarly, D. sexcinctus Lineage 1 may have experienced a recent demographic expansion since R2 was significant and this statistic is particularly sensitive with small sample sizes (Table 2). Within L. bellii, Lineage B also showed signs of recent demographic expansion with significant Fs and R2 statistics (Fig. 4, Table 2). Multilocus Extended Bayesian Skyline Plots (EBSP) are

consistent with these results and identified recent changes in the effective population size (Figs. 4, 5). The EBSP for D. darwinii showed a drastic constriction of the effective population size during the last 0.5 MYA and a recent population expansion. Conversely, the EBSP of D. sexcinctus sensu stricto showed an expansion of the effective population size since the last 0.2 MYA. The analysis of D. bibronii and D. sexcinctus Lineage 2 showed little variation in their effective sizes. Diplolaemus bibronii declined steadily and it was a little more abrupt during the last 0.5 MYA with a slight increase since 0.2 MYA (Fig. 5). Similarly, L. bellii showed an expansion of the effective population size since the last 1 MYA (Fig. 4). 3.5. Species delimitation using two mitochondrial and eight nuclear genetic markers Based on the species-delimitation methods implemented here, we obtained two models with high posterior probability. The model estimated with BP&P (PP = 0.80) inferred D. bibronii as a single species lineage (PP = 0.89), suggesting 23 groups in total (including Urostrophus) with posterior probabilities > 0.80. The model estimated by STACEY (PP = 0.98) delimited 24 species lineages, including two from D. bibronii. In this second analysis, most of the lineages were inferred with posterior probabilities > 0.80 except for P. nigroiugulus (PP = 0.66), P. fasciatus (PP = 0.66), D. s. Lineage 1 (PP = 0.31), D. s. Lineage 3 (PP = 0.60), and D. s. Lineage 4 (PP = 0.60). We conservatively use the model estimated with BP&P to define independent lineages in the reconstruction of species trees and biogeographic history. The posterior probabilities estimated for each lineage are summarized in Figure S3. 3.6. Species tree and divergence times using two mitochondrial and eight nuclear genetic markers The species tree inferred most relationships with high statistical support (Fig. 6). We obtained the Leiosaurae as a monophyletic group (PP=0.98). Except Pristidactylus torquatus, which was inferred as the sister species to all other taxa, the three genera of Leiosaurae were inferred as monophyletic groups with high support (PP = 1-0.99), but relationships among them have low support. Within Leiosaurus, all species relationships are highly supported, and within L. bellii the (Lineage B + Lineage C) topology has moderate support (PP = 0.85). Within the rest of Pristidactylus species relationships had moderate support (PP =

0.82-0.94). Within Diplolaemus (PP=1), a strongly supported (PP=1) (D. darwinii Lineage A + D. darwinii Lineage B) clade was inferred as the sister taxon to a clade comprising all other terminals with high support (PP=1). Diplolaemus leopardinus is inferred as the sister taxon to the clade comprising D. bibronii, D. sexcinctus and D. s. lineages 1-5, but without statistical support. Within this clade, subclades (D. s. Lineage 5 + D. sexcinctus sensu stricto) and (D. s. Lineage 2 + D. bibronii) received high statistical support (PP = 0.99 and PP = 0.96 respectively). Finally, clade (D. s. Lineage 3 + D. s. Lineage 4) has moderate support (PP = 0.74), and the position of D. s. Lineage 1 is not resolved (Fig. 6). Calibration and divergence times indicate an origin of the Leiosaurae early in the Miocene. Pristidactylus torquatus separated from a common ancestor of the rest of the Leiosaurae clade ~ 16.8 MYA, and the crown ancestors of the other genera may have originated during the Miocene. Within the Diplolaemus clade the ancestor of D. darwinii diverged from a common ancestor of all other species perhaps during the Pliocene, while the origin of the remaining species is inferred during the Quaternary (Fig. 6). 3.7. Historical biogeography The best model in the BioGeoBEARS analysis for Diplolaemus was BAYAREALIKE + j (Relative probability = 0.70, Table 4A). When comparing the model with and without the 'J' parameter of founder-effect, we found significant differences, indicating that incorporation of the founder effect parameter provided better likelihood fit to the data. The best model for L. bellii was DIVALIKE (Table 4 B), and we found no significant differences in models including and excluding the J parameter, i.e. the two models confer the same likelihood on the data. Our biogeographic reconstruction (Fig. 7) suggests that ancestral populations of Diplolaemus were distributed over a wide region in southern Patagonia (ABCJ), across what are now the provinces of Chubut and Santa Cruz. The most recent common ancestor (MRCA) of D. darwinii could have been retracted in the southern region of Patagonia (AJ), and the ancestral populations of D. leopardinus possibly migrated towards the north early in the Pleistocene, colonizing the region of current Mendoza province. Similarly, the

ancestral populations of D. s. Lineage 3 and D. s. Lineage 4 probably migrated to the north occupying the current northern Neuquén province (region G in Fig. 7). Later, D. s. Lineage 3 colonized the southwest of current Chubut province (region C, Fig. 7). Ancestral populations of D. sexcinctus sensu stricto, D. s. Lineage 1 and D. s. Lineage 5 probably emerged from a refuge northwest of Chubut province (region K, Fig. 7). Diplolaemus sexcinctus and D. s. Lineage 5 then migrated eastward, remaining in the Somun Curá plateau for the last 500 thousand years (Fig. 7). The estimated range of distribution for the ancestral populations of L. bellii corresponds to the entire centralnorthern region of Patagonia (Fig. 7). Leiosaurus bellii Lineage B + L. bellii Lineage C diverged from L. bellii Lineage A ~ 0.48 MYA. Most probably L. bellii Lineage A remained in a reduced ancestral range in the south (B, C and D in Fig. 7) and L.bellii Lineage B and Lineage C occupied the area defined as L in the province of Mendoza. Probably the populations expanded after the divergence of L. bellii Lineage B and Lineage C, so the current ranges are more extensive (Fig. 7). 4. Discussion 4.1. Genetic structuring, species boundaries and phylogenetic relationships Our results provide evidence of intraspecific genetic structure that mostly corresponds to discrete geographic areas, revealing a previously unknown diversity for the Leiosaurae. This evidence is based on: 1) the cytochrome-b gene tree that inferred haploclades with high statistical support with both IB and ML analyses; 2) maximum parsimony networks; 3) pairwise mitochondrial genetic distances >3% among the different species-level lineages; and 4) multilocus methods that delimited putative species congruent with the identity of each cytochrome b clade as a genetically isolated population. Diplolaemus darwinii is well-differentiated (Figs. 1, 6), with mitochondrial genetic distances >16% compared to the rest of the Diplolaemus species/lineages (Table 3). The molecular differentiation found in this species with respect to others of the genus is consistent with the serological (Cei, 1975), meristic, and color-pattern analyses (Victoriano et al., 2010). Slopes of Pampa del Castillo and Pampa de Salamanca separate

geographically D. darwinii Lineages A and B within this species (Fig. 3). This region has previously been identified as a possible refugium for other lizards (Avila et al., 2006; Breitman et al., 2012) and plants (Cosacov et al., 2010). However, our sample of D. darwinii Lineage B includes two haplotypes north of these slopes (Fig. 3), highlighting the need for additional sampling to resolve its taxonomic status. Diplolaemus bibronii is separated by cytochorme b genetic distances > 5% from all other congeners (Table 3). Intraspecific variation is low (π ~0.011, Table 2) in this species, but the same geographic features identified for D. darwinii (slopes of Pampa del Castillo and Pampa de Salamanca) separate two groups (Lineage A and Lineage B, Figs. 1, 3) that are delimited with STACEY but not with BP&P (Fig. S3). Diplolaemus leopardinus was described by Werner (1898) and assigned to Santiago (Chile) as an in error type locality. Later the locality was corrected by Donoso-Barros (1966), suggesting a distribution in the region of Lonquimay, Mari Menú and Pino Hachado in Chile, and the region of Aluminé Lake in Neuquén, Argentina. However, Victoriano et al. (2010) based on morphological and meristic characters, concluded that individuals assigned as D. leopardinus by Donoso-Barros (1966) from Laguna del Laja (Chile), correspond to the species D. sexcinctus, and suggested that the distribution of D. leopardinus is limited to southern Mendoza province. Here we included for the first time populations from Mendoza in a phylogeny, and provide evidence corroborating Victoriano et al. (2010); the distribution of D. leopardinus is restricted to Mendoza province. The individuals from Paso Pichachén (near to Laguna del Laja National Park, Bío-Bío Region) (#14190 and #14191) included in our analyses are inferred as D. s. Lineage 2 (Fig. 1) in northern Neuquén province, independent from Mendoza samples in the cytochrome-b gene tree. Diplolaemus leopardinus and D. s. Lineage 2 are separated by mitochondrial genetic distances of 8.28% (Table 3), and were inferred as independent lineages in the species-delimitation analyses (Fig. S3). The species tree shows D. leopardinus as a well-differentiated lineage compared to D. sexcinctus (sensu Victoriano et al.2010) and D. bibronii, while D. s. Lineage 2 is inferred as closely related to D. bibronii in the species tree (Fig. 6). In summary, our analyses infer D. leopardinus as a well differentiated lineage distributed in the Mendoza province, and not grouped with the Bío-Bío populations.

The nomenclatural problems for D. leopardinus have also confused the taxonomic history of D. sexcinctus. Cei (1986) recognized two morphotypes in Diplolaemus populations from southern Chubut to Mendoza provinces: "sudmendocina" (populations in southern Mendoza) and "altopatagónica" (Andean and subAndean populations). Subsequently, Cei et al. (2003) formally described these morphotypes on the basis of comparisons among individuals from Santa Cruz, Chubut, Río Negro and Mendoza provinces, and described a new species, D. sexcinctus, widely distributed from northern Chico River (Chubut) to Mendoza province. However, Cei et al. (2003) did not include populations from Bío-Bío region (Chile), which at the time were assigned to D. leopardinus (Donoso-Barros, 1966). In addition, as discussed in the previous paragraph, populations from Mendoza, which included the holotype of D. sexcinctus by Cei et al. (2003), should be assigned to D. leopardinus. Our results clarify relationships among different populations throughout this wide distribution of D. sexcinctus. Based on the cytochrome-b gene, we inferred six well-differentiated clades mostly corresponding to defined geographic areas (Fig. 1), and almost all are separated by pairwise genic distances >3% (except D. s. Lineage 5). The multilocus species-delimitation analyses show support for recognition of three species with BP&P and STACEY (D. s. Lineage 2, D. s. Lineage 5, and D. sexcinctus sensu stricto), and three others with STACEY alone (D. s. Lineage 1, D. s. Lineage 3, D. s. Lineage 4). Diplolaemus Lineage 5 is phylogenetically grouped with D. sexcinctus, and the two are sympatrically distributed on the Somun Curá plateau. As mentioned before, D. s. Lineage 2 is grouped with D. bibronii, and phylogenetic relationships among the remaining lineages were not resolved. The genetic divergence and geographic structure of these lineages provide evidence for their recognition as “candidate species”, but an integrative taxonomic study based on additional molecular markers, multiple classes of morphological data, and ecological information, is necessary to assess their status (Sukumaran and Knowles, 2017). Within L. bellii three lineages were identified; L. bellii Lineage A and L. bellii Lineage B are allopatrically separated by the Negro River (Fig. 1), and by mitochondrial genetic distances > 5% (Table 3). L. belli Lineage C is distributed in northwestern Neuquén province and south of Mendoza province (Fig. 4).

Although its geographic distribution overlaps with the other two lineages of L. bellii, L. bellii Lineage C differs from both by mitochondrial genetic distances > 5%. The three lineages were resolved as independent networks separated by several mutational steps, and with the species-delimitation analyses all three were inferred as independent lineages (Fig. S3). This evidence suggests strong genetic structuring that may correspond to new species. Given that the type locality of L. bellii is unclear, we conservatively do not assign any of these lineages to the type species; further integrative taxonomic analyses are needed to clarify species boundaries in this clade. 4.2 Divergence, biogeographic history, and recent demographic changes in Diplolaemus. The BayArea+J model suggests that processes of extinction and dispersal explain the biogeographic history of Diplolaemus. Cladogenesis within the genus began during the Pliocene ~ 3.35 MYA, which is within the range suggested for the second Eulaemus lizard radiation (Olave et al., 2015). The MRCA of Diplolaemus was probably widely distributed throughout southern Patagonia, and we estimated that D. darwinii diverged early (Pliocene) in the history of the genus, a possible explanation for the presence of keeled caudal scales exclusively in this species (a plesiomorphic character for the genus; Etheridge and Williams, 1985; Victoriano et al., 2010). This divergence coincides with an increase in the aridity of the region, and the formation of steppe vegetation, which persisted from ~ 3.35 MYA (Jakob et al. 2009). This ecological shift may have driven the divergence of D. darwinii, which is restricted now to the steppe environment. The Great Patagonian Glaciation (GPG) occurred ~ 1-1.2 MYA (GPG, Mercer, 1976), and coincides with the divergence between the two D. darwinii lineages. One hypothesis is that ancestral populations separated during the regional uplift of the Chubut river basin, and remained isolated during the GPG (Lapido, 1981; Rabassa, 2008). The slopes of Pampa del Castillo and Pampa de Salamanca probably kept these populations isolated during the interglacial periods, and while we did not detect demographic changes during the GPG, our EBSP analysis showed an abrupt population size reduction after the Coldest Pleistocene Glaciation (CPG, ~ 0.7 MYA, Rabassa et al., 2008). We hypothesize that the CPG could have led to drastic population reduction, retracting some populations towards isolated refuges; subsequently,

when conditions were favorable, these populations could have experienced an expansion. This hypothesis is supported by neutrality deviation (Tajima’s D = -1.89, p=0.014), low nucleotide diversity indexes (π=0.011, θw=0.022; Table 2), but high haplotype diversity (Hd=0.94). These haplotypes are geographically structured in regions where glacial refuges have been proposed for other species (Breitman et al. 2012, Sérsic et al. 2011). Similarly, D. bibronii, the other southern distributed species, also shows a slight reduction of its Ne at ~0.5 MYA (post-CPG), but this species shows a subsequent increase in population size during the last 0.2 MYA in the EBSP analyses (Fig. 5), which coincides with its divergence with respect to D. s. Lineage 2 (Fig. 7). Geological evidence indicates that the ice sheet moved quickly in this region of Patagonia in the Middle Pleistocene (Hein et al., 2017; Rabassa et al., 2005), and these fluctuations may have influenced the historical demographic changes of animal and plant populations (Breitman et al., 2012; Cosacov et al., 2010; Sérsic et al., 2011). However, the responses to these post-GPC fluctuations seem to have been different for these two lizard species. While Diplolaemus bibronii seems that it was able to maintain its population size, D. darwinii seems to have suffered a 10-fold population size reduction (Fig. 5). The neutrality test shows signs of population expansion in the southernmost lineage of D. bibronii (D. bibronii Lineage B, Table 2). Therefore, we suggest that D. bibronii may have remained in a refuge in the slopes of Pampa del Castillo and Pampa de Salamanca during the CPG and during the post-CPG period, the river reactivation may have separated the D. bibronii lineages A and B. Breaks have been suggested in this region for plants (Cosacov et al., 2010) and for other species of lizards (Breitman et al., 2012). Subsequently, the fall of sea-level exposed large areas to the east that could have allowed the population expansion of D. bibronii Lineage B. On the other hand, we suggest that D. darwinii remained in the southern regions of Patagonia during the CPG and subsequently experienced a drastic population reduction. The divergence time between Diplolaemus sexcinctus Lineage 3 and D. s. Lineage 4 from populations distributed in central and southern Patagonia was estimated at ~ 1.17 MYA. This divergence is coincident with the Great Patagonian Glaciation (~ 1-1.2 MYA; Mercer, 1976). Immediately after the GPG, the increased erosion related to the tectonic ascent of the Patagonian Andes increased the deepening of the

valleys (Rabassa, 2008; Rabassa et al., 2011). This could have isolated the ancestral D. s. Lineage 4 in the center of Neuquén province, while uplift of the Chubut River basin could have isolated the ancestral D. s. Lineage 3 in southeastern Chubut province, during the interglacial period. We found exclusive haplotypes of these lineages in parsimony networks (Fig. 2) and pairwise genetic distances >3% (Table 3), congruent with a scenario of isolation and limited gene flow. From Pliocene to Holocene, volcanic eruptions were active in northwestern areas of the current Somún Curá Plateau (Kay et al., 2007, 2006; Stern et al., 1990; Fig. 2), favoring allopatric isolation of populations. This high basaltic plateau was not affected by marine transgressions or ice sheets. Thus, this region served as a refugium for various species of animals and plants (Nicola et al., 2014; Sérsic et al., 2011), and this isolation could have promoted cladogenesis in multiple groups of organisms. The lineages currently distributed in this region (Diplolaemus sexcinctus sensu stricto, D. s. Lineage 5 and D. s. Lineage 1) diverged ~ 0.98 MYA when the ice sheet began to retreat, probably occupying the central-western region of the Río Negro province and north-central Chubut province (Fig. 7). We suggest that the ancestral populations of these lineages could have remained in a glacial refuge during the GPG. Probably the postglacial fluvial reactivation of the Limay River favored the vicariant process that separated D. s. Lineage 1 from the ancestral populations of D. s. Lineage 5 + D. sexcinctus sensu stricto. We estimated that this divergence occurred ~ 0.88 MYA (Fig. 7); subsequently, D.s. Lineage 1 experienced a recent population expansion (Table 2). We further propose that this divergence could be linked to the isolation of ancestral populations in the center of the Somún Curá plateau during the Coldest Pleistocene Glaciation (~ 0.7 MYA, Rabassa et al., 2011). We estimate the divergence of D. sexcinctus sensu stricto and D. s. Lineage 5 to have occurred ~0.5 MYA. Subsequently, D. sexcinctus sensu stricto experienced a population increase in the last 0.2 MYA (Fig. 5). Probably after the CPG, when the conditions were more favorable, large areas were exposed to the north and center of the Somún Curá plateau. This could have favored the expansion of the populations that remained in the center of the plateau during the GPG. Diplolaemus leopardinus diverged from a common ancestor of D. bibronii and D. sexcinctus ~ 1.43 MYA. The volcanic activity and the

glaciofluvial cycles of the Colorado River basin that occurred in the Pleistocene could have promoted their dispersion, subsequent isolation and founder-effect-mediated cladogenesis. Southern Mendoza province has been recognized to harbor various endemic species of lizards (Corbalán et al., 2011; Medina et al., 2014) and other vertebrates (García-R et al., 2012; Ojeda, 2010), a different ensemble of biota from the southern Neuquén province, consistent with the differentiation that we found for D. leopardinus within this genus. The inference of the biogeographic history for L. bellii suggests a wide ancestral distribution range in the north-central Patagonia. We estimated that the L. bellii Lineage A divergence with respect to L. bellii Lineage B + L. bellii Lineage C began ~ 0.48MYA. Probably the deglaciation that caused the fluvial reactivation of the Colorado River caused the vicariant event. After this separation we estimated an immediate reduction of the ancestral ranges of each lineage with a subsequent expansion (Fig. 7). This process is also supported with the EBSP analyses, in which we estimated a growing Ne from the GPG. We also suggest that the expansion of L. bellii Lineage B, probably occupying the southern regions, was the most recent since we found deviations from the neutral model suggesting recent demographic changes (Table 2). 5. Conclusions Many works have explored the diversity of Patagonian biota with phylogeographic approaches; however, the clade Leiosaurae had not been studied with this framework given the difficulty of finding them in the field. In this work, we obtained samples collected over 9 years of continuous field work (2006-2015) and performed molecular analyses at the intra- and interspecific levels. Our results revealed an unknown diversity. We found eight lineages (5 within Diplolaemus and three within L. bellii), and we suggest that most of them could correspond to putative new species. For example, within Diplolaemus species we found strong evidence for D. s. lineages 1-4 based on geographical structuring, cytochrome b distances and multilocus species-delimitation methods. However, D. s. lineage 5 and D. sexcinctus sensu stricto are co-distributed in the central region of the Somun Curá plateau and show cytochrome b distances of 2.59%. Although

multilocus species delimitation methods suggest that they could be two candidate species, additional evidence is required to corroborate it. For L. bellii we found strong evidence for recognition of lineages A and B. These lineages are separated by high pairwise cytochrome b distances, and the Colorado River seems to be a barrier separating them. Leisosaurus bellii lineage C is separated by high cytochrome b distances from L. bellii lineages A and B, but is co-distributed with them. We recommend future in-depth integrative analyses incorporating morphological and ecological data to test their specific status. Our work also provides useful information for planning of conservation actions. First, we clarify the relationships and distributions that have brought confusion in the Diplolaemus species. Our results indicate that D. leopardinus is a geographically restricted species to Mendoza province, congruent with the results of Victoriano et al. (2010), and populations of Bío-Bío (Chile) and northern Neuquén province should be assigned to D. s. Lineage 2. In addition, we found that D. darwinii populations distributed in southern Patagonia have a very low nucleotide but high haplotype diversities, and probably experienced a population bottleneck in the recent past. Therefore, we warn that its conservation status should be revised. Finally, we generated hypotheses of the historical biogeography of each of the lineage/species we studied. We found that the geological and paleoclimatic events from the Pliocene to the present seemed to have produced the diversity of these species. These hypotheses can be contrasted with the biogeographic history of other species in the region, to help our understanding on the origin and evolution of Patagonia current diversity. 6. Acknowledgments We thank all past members of the Grupo de Herpetología Patagónica (GHP) for assistance in field collections, in the laboratory, and/or assistance in animal curation procedures. We also want to thank the AE Allan Larson for his positive comments on an early draft and three anonymous reviewers that made useful observations that helped to improve the manuscript. Financial support was provided by the following grants: ANPCYT-FONCYT (LJA) PICT 506/2006, ANPCYT-FONCYT 1397/2011; 1252/2015, PIP-CONICET 0336/13, UNPSJB 2015 (MM), and a doctoral fellowship (MF) from Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET). The NSF-PIRE award (OISE 0530267) supported collaborative research on Patagonian Biodiversity to the following institutions (listed alphabetically): Brigham Young University, Centro Nacional Patagónico (AR), Dalhousie University, Instituto Botánico Darwinion (AR), Universidad

Austral de Chile, Universidad de Concepción (CH), Universidad Nacional del Comahue (AR), Universidad Nacional de Córdoba (AR) and University of Nebraska. We thank the fauna authorities from Río Negro, Neuquén, Mendoza, Chubut and Santa Cruz Provinces for collection permits. 7. Declarations of interest: none.

8. Bibliography Abdala, V., Manzano, A., Nieto, L., Diogo, R., 2009. Comparative myology of Leiosauridae (Squamata) and its bearing on their phylogenetic relationships. Belgian J. Zool. 139, 109–123. Albino, A.M., 2008. Lagartos iguanios del Colhuehuapense (Mioceno Temprano) de Gaiman (provincia del Chubut, Argentina). Ameghiniana 45, 775–782. Albino, A.M., Brizuela, S., Vizcaíno, S., 2016. The southernmost fossil record of squamates. AmphibiaReptilia 1–16. https://doi.org/10.1163/15685381-00003078 Andrade-Díaz, M.S., Hibbard, T.N., Díaz-Gómez, J.M., 2017. Identifying endemism areas: an example using neotropical lizards. South Am. J. Herpetol. 12, 61–75. https://doi.org/10.2994/SAJH-D-16-00038.1 Avila, L.J., Morando, M., Sites Jr., J.W., 2006. Congeneric phylogeography: hypothesizing species limits and evolutionary processes in Patagonian lizards of the Liolaemus boulengeri group (Squamata: Liolaemini). Biol. J. Linn. Soc. 89, 241–275. https://doi.org/10.1111/j.1095-8312.2006.00666.x Bickford, D., Lohman, D.J., Sodhi, N.S., Ng, P.K.L., Meier, R., Winker, K., Ingram, K.K., Das, I., 2007. Cryptic species as a window on diversity and conservation. Trends Ecol. Evol. 22, 148–155. https://doi.org/10.1016/j.tree.2006.11.004 Bouckaert, R.R., 2010. DensiTree: making sense of sets of phylogenetic trees. Bioinformatics 26, 1372– 1373. https://doi.org/10.1093/bioinformatics/btq110 Bouckaert, R.R., Heled, J., Kühnert, D., Vaughan, T., Wu, C.-H., Xie, D., Suchard, M.A., Rambaut, A., Drummond, A.J., 2014. BEAST 2: a software platform for Bayesian evolutionary analysis. PLOS Comput. Biol. 10, e1003537. Breitman, M.F., Avila, L.J., Sites Jr., J.W., Morando, M., 2012. How lizards survived blizzards: phylogeography of the Liolaemus lineomaculatus group (Liolaemidae) reveals multiple breaks and refugia in southern Patagonia and their concordance with other codistributed taxa. Mol. Ecol. 21, 6068– 85. https://doi.org/10.1111/mec.12075 Camargo, A., Sites, J., 2013. Species delimitation: a decade after the renaissance, in: Pavlinov, I. (Ed.), The Species Problem - Ongoing Issues. InTech, pp. 225–247. https://doi.org/10.5772/52664 Cei, J.M., 1986. Reptiles del centro-oeste y sur de la Argentina. Herpetofauna de las zonas áridas y semiáridas. Monogr. IV. Muzeo Reg. di Sci. Nat. Torino, Italy. Cei, J.M., 1975. Diferenciación serológica de Diplolaemus darwinii y Diplolaemus bibronii en poblaciones alo-simpátridas. Physis 34, 209–210. Cei, J.M., Scolaro, J.A., Videla, F., 2003. A taxonomic revision of recognized Argentine species of the leiosaurid genus Diplolaemus (Reptilia, Squamata, Leiosauridae). Facena 19, 87–106. Clement, M., Posada, D., Crandall, K.A., 2000. TCS: a computer program to estimate gene genealogies.

Mol. Ecol. 9, 1657–1659. https://doi.org/10.1046/j.1365-294x.2000.01020.x Corbalán, V., Tognelli, M.F., Scolaro, J.A., Roig-Juñent, S.A., 2011. Lizards as conservation targets in Argentinean Patagonia. J. Nat. Conserv. 19, 60–67. https://doi.org/10.1016/j.jnc.2010.05.004 Cosacov, A., Johnson, L.A., Paiaro, V., Cocucci, A.A., Córdoba, F.E., Sérsic, A.N., 2013. Precipitation rather than temperature influenced the phylogeography of the endemic shrub Anarthrophyllum desideratum in the Patagonian steppe. J. Biogeogr. 40, 168–182. https://doi.org/10.1111/j.13652699.2012.02776.x Cosacov, A., Sérsic, A.N., Sosa, V., Johnson, L.A., Cocucci, A.A., 2010. Multiple periglacial refugia in the Patagonian steppe and post-glacial colonization of the Andes: the phylogeography of Calceolaria polyrhiza. J. Biogeogr. 37, 1463–1477. https://doi.org/10.1111/j.1365-2699.2010.02307.x Cuitiño, J.I., Fernicola, J.C., Kohn, M.J., Trayler, R., Naipauer, M., Bargo, M.S., Kay, R.F., Vizcaíno, S.F., 2016. U-Pb geochronology of the Santa Cruz Formation (early Miocene) at the Río Bote and Río Santa Cruz (southernmost Patagonia, Argentina): implications for the correlation of fossil vertebrate localities. J. South Am. Earth Sci. 70, 198–210. https://doi.org/10.1016/j.jsames.2016.05.007 De Queiroz, K., 2007. Species concepts and species delimitation. Syst. Biol. 56, 879–886. https://doi.org/10.1080/10635150701701083 Donoso-Barros, R., 1966. Reptiles de Chile. Ediciones de la Universidad de Chile, Santiago de Chile. Drummond, A.J., Rambaut, A., Shapiro, B., Pybus, G., 2005. Bayesian coalescent inference of past population dynamics from molecular sequences. Mol. Biol. Evol. 22, 1185–1192. https://doi.org/10.1093/molbev/msi103 Etheridge, R., Williams, E., 1985. Notes on Pristidactylus (Squamata: Iguanidae). Breviora 483, 1–18. Excoffier, L., Lischer, H.E.L., 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567. https://doi.org/10.1111/j.1755-0998.2010.02847.x Fontanella, F.M., Olave, M., Avila, L.J., Sites, J.W., Morando, M., 2012. Molecular dating and diversification of the South American lizard genus Liolaemus (subgenus Eulaemus) based on nuclear and mitochondrial DNA sequences. Zool. J. Linn. Soc. 164, 825–835. https://doi.org/10.1111/j.10963642.2011.00786.x Frost, D.R., Etheridge, R., Janies, D., Titus, T.A., 2001. Total evidence, sequence alignment, evolution of polychrotid lizards, and a reclassification of the Iguania (Squamata: Iguania). Am. Museum Novit. 3343, 1–39. https://doi.org/10.1206/0003-0082(2001)343<0001:TESAEO>2.0.CO;2 Fu, Y.-X., 1997. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147, 915–925. García-R, J.C., Crawford, A.J., Mendoza, Á.M., Ospina, O., Cardenas, H., Castro, F., 2012. Comparative phylogeography of direct-developing frogs (Anura: Craugastoridae: Pristimantis) in the southern Andes of Colombia. PLoS One 7, e46077. https://doi.org/10.1371/journal.pone.0046077 Hein, A.S., Cogez, A., Darvill, C.M., Mendelova, M., Kaplan, M.R., Herman, F., Dunai, T.J., Norton, K., Xu, S., Christl, M., Rodés, Á., 2017. Regional mid-Pleistocene glaciation in central Patagonia. Quat. Sci. Rev. 164, 77–94. https://doi.org/https://doi.org/10.1016/j.quascirev.2017.03.023 Heled, J., Drummond, A.J., 2010. Bayesian inference of species trees from multilocus data. Mol. Biol. Evol. 27, 570–580. https://doi.org/10.1093/molbev/msp274

Hillis, D.M., Bull, J.J., 1993. An empirical test of bootstrapping as a method for assessing confidence in phylogenetic analysis. Syst. Biol. 42, 182-192. https://doi.org/10.1093/sysbio/42.2.182 Ho, S.Y., Shapiro, B., 2011. Skyline-plot methods for estimating demographic history from nucleotide sequences. Mol. Ecol. Resour. 11, 423–34. https://doi.org/10.1111/j.1755-0998.2011.02988.x Hudson, R.., 1990. Gene genealogies and coalescent process, in: Futuyma, D.J., Antonovics, J. (Eds.), Oxford Surveys in Evolutionary Biology. Oxford: Oxford University Press, pp. 1–44. Huelsenbeck, J.P., Ronquist, F., 2001. MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics 17, 754–755. https://doi.org/10.1093/bioinformatics/17.8.754 Jakob, S.S., Martinez-Meyer, E., Blattner, F.R., 2009. Phylogeographic analyses and paleodistribution modeling indicate Pleistocene in situ survival of Hordeum species (Poaceae) in Southern Patagonia without genetic or spatial restriction. Mol. Biol. Evol. 26, 907–923. Jones, G., 2017. Algorithmic improvements to species delimitation and phylogeny estimation under the multispecies coalescent. J. Math. Biol. 74, 447–467. https://doi.org/10.1007/s00285-016-1034-0 Jones, G., Aydin, Z., Oxelman, B., 2014. DISSECT: an assignment-free Bayesian discovery method for species delimitation under the multispecies coalescent. Bioinformatics 31, 991–998. Katoh, K., Standley, D.M., 2013. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 30, 772–780. https://doi.org/10.1093/molbev/mst010 Kay, S.M., Ardolino, A.A., Gorring, M.L., Ramos, V.A., 2007. The somuncura large igneous province in Patagonia: interaction of a transient mantle thermal anomaly with a subducting slab. J. Petrol. 48, 43– 77. https://doi.org/10.1093/petrology/egl053 Kay, S.M., Burns, W.M., Copeland, P., Mancilla, O., 2006. Upper Cretaceous to Holocene magmatism and evidence for transient shallowing of the Andean subduction zone under the northern Neuquen basin, in: Kay, S.M., Ramos, V.A. (Eds.), Evolution of an Andean Margin: A Tectonic and Magmatic Perspective from the Andes to the Neuquen Basin (35–39S Lat.). Geological Society of America, p. 407. Lamoreux, J.F., Morrison, J.C., Ricketts, T.H., Olson, D.M., Dinerstein, E., McKnight, M.W., Shugart, H.H., 2006. Global tests of biodiversity concordance and the importance of endemism. Nature 440, 212–214. Landis, M.J., Matzke, N.J., Moore, B.R., Huelsenbeck, J.P., 2013. Bayesian analysis of biogeography when the number of areas is large. Syst. Biol. 62, 789–804. https://doi.org/10.1093/sysbio/syt040 Lanfear, R., Calcott, B., Ho, S.Y.W., Guindon, S., 2012. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol. Biol. Evol. 29, 1695. https://doi.org/10.1093/molbev/mss020 Lapido, O.R., 1981. Descripción Geológica de la Hoja 44g -Cañadon Iglesias-, Provincia de Chubut. Buenos Aires, Argentina. Librado, P., Rozas, J., 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. https://doi.org/10.1093/bioinformatics/btp187 Maddison, W.P. and D.R.M., 2018. Mesquite: a modular system for evolutionary analysis. http://mesquiteproject.org. Matzke, N.J., 2013. BioGeoBEARS: BioGeography with Bayesian (and likelihood) evolutionary analysis in R Scripts. https://cran.r-project.org/web/packages/BioGeoBEARS/index.html. Matzke, N.J., 2014. Model selection in historical biogeography reveals that founder-event speciation is a

crucial process in island clades. Syst. Biol. 63, 951–970. https://doi.org/10.1093/sysbio/syu056 Medina, C.D., Avila, L.J., Sites Jr., J.W., Morando, M., 2014. Multilocus phylogeography of the Patagonian lizard complex Liolaemus kriegi (Iguania: Liolaemini). Biol. J. Linn. Soc. 113, 256–269. https://doi.org/10.1111/bij.12285 Mercer, J.H., 1976. Glacial history of southernmost South America. Quat. Res. 6, 125–166. https://doi.org/10.1016/0033-5894(76)90047-8 Morando, M., Avila, L.J., Baker, J., Sites, J.W., 2004. Phylogeny and Phylogeography of the Liolaemus darwinii Complex (Squamata: Liolaemidae): evidence for Introgression and Incomplete Lineage Sorting. Evolution (N. Y). 58, 842–861. https://doi.org/10.1554/03-009 Morando, M., Avila, L.J., Sites Jr., J.W., Sullivan, J., 2003. Sampling strategies for delimiting species: genes, individuals, and populations in the Liolaemus elongatus-kriegi complex (Squamata: Liolaemidae) in Andean-Patagonian South America. Syst. Biol. 52, 159–185. https://doi.org/Doi 10.1080/10635150390192717 Morando, M., Avila, L.J., Turner, C.R., Sites Jr., J.W., 2007. Molecular evidence for a species complex in the patagonian lizard Liolaemus bibronii and phylogeography of the closely related Liolaemus gracilis (Squamata: Liolaemini). Mol. Phylogenet. Evol. 43, 952–973. https://doi.org/10.1016/j.ympev.2006.09.012 Morando, M., Olave, M., Avila, L.J., Baker, E., Sites, J.W., 2015. Molecular phylogeny of the lizard clade Leiosaurae endemic to southern South America. Herpetologica 71, 322–331. Nei, M.,1987. Molecular Evolutionary Genetics. New York: Columbia University Press. Nicola, M. V., Sede, S.M., Pozner, R., Johnson, L.A., 2014. Phylogeography and palaeodistribution modelling of Nassauvia subgenus Strongyloma (Asteraceae): exploring phylogeographical scenarios in the Patagonian steppe. Ecol. Evol. 4, 4270–4286. https://doi.org/10.1002/ece3.1268 Ogilvie, H.A., Bouckaert, R.R., Drummond, A.J., 2017. StarBEAST2 brings faster species tree inference and accurate estimates of substitution rates. Mol. Biol. Evol. 34, 2101–2114. https://doi.org/10.1093/molbev/msx126 Ojeda, A.A., 2010. Phylogeography and genetic variation in the South American rodent Tympanoctomys barrerae (Rodentia: Octodontidae). J. Mammal. 91, 302–313. Olave, M., Avila, L.J., Sites, J.W., Morando, M., 2015. Model-based approach to test hard polytomies in the Eulaemus clade of the most diverse South American lizard genus Liolaemus (Liolaemini, Squamata). Zool. J. Linn. Soc. 174, 169–184. https://doi.org/10.1111/zoj.12231 Palumbi, S.., 1996. PCR and molecular systematics, in: Hillis, D., Moritz, C., Mable, B. (Eds.), Molecular Systematics. Sinauer Press, USA. Perkins, M.E., Fleagle, J.G., Heizler, M.T., Nash, B., Bown, T.M., Tauber, A.A., Dozo, M.T., 2012. Tephrochronology of the Miocene Santa Cruz and Pinturas formations, Argentina, in: Vizcaíno, S.F., Kay, R.F., Bargo, M.S. (Eds.), Early Miocene Paleobiology in Patagonia: High-Latitude Paleocommunities of the Santa Cruz Formation. Cambridge University Press, USA, pp. 23–40. Pyron, R.A., Burbrink, F.T., Wiens, J.J., 2013. A phylogeny and revised classification of Squamata, including 4161 species of lizards and snakes. BMC Evol. Biol. 13. https://doi.org/doi:10.5061/dryad.82h0m R Core Team, 2017. R: A software language and environment for statistical computing. https://www.Rproject.org/

Rabassa, J., 2008. Late Cenozoic Glaciations in Patagonia and Tierra del Fuego, in: Sciences, J.R.B.T.-D. in Q. (Ed.), The Late Cenozoic of Patagonia and Tierra Del Fuego. Elsevier, pp. 151–204. https://doi.org/http://dx.doi.org/10.1016/S1571-0866(07)10008-7 Rabassa, J., Coronato, A., Martínez, O., 2011. Late Cenozoic glaciations in Patagonia and Tierra del Fuego: an updated review. Biol. J. Linn. Soc. 103, 316–335. https://doi.org/10.1111/j.10958312.2011.01681.x Rabassa, J., Coronato, A.M., Salemme, M., 2005. Chronology of the Late Cenozoic Patagonian glaciations and their correlation with biostratigraphic units of the Pampean region (Argentina). J. South Am. Earth Sci. 20, 81–103. https://doi.org/10.1016/j.jsames.2005.07.004 Rambaut, A., Suchard, M.A., Xie, D., Drummond, A.J., 2014. Tracer 1.6, http://beast.bio.ed.ac.uk. Ramos-Onsins, S.E., Rozas, J., 2002. Statistical properties of new neutrality tests against population growth. Mol. Biol. Evol. 19, 2092–2100. Rannala, B., Yang, Z., 2013. Improved reversible jump algorithms for Bayesian species delimitation. Genetics 194, 245–253. https://doi.org/10.1534/genetics.112.149039 Ré, G.H., Bellosi, E.S., Heizler, M., Vilas, J.F., Madden, R.H., Carlini, A.A., Kay, R.F., Vucetich, M.G., 2010. A geochronology for the Sarmiento Formation at Gran Barranca, in: Madden, R.H., Carlini, A.A., Vucetich, M.G., Kay, R.F. (Eds.), The Paleontology of Gran Barranca: Evolution and Environmental Change through the Middle Cenozoic of Patagonia. Cambridge University Press, UK, pp. 46–60. Ree, R.H., 2005. Detecting the historical signature of key innovations using stochastic models of character evolution and cladogenesis. Evolution 59, 257–265. Ree, R.H., Smith, S.A., 2008. Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis. Syst. Biol. 57, 4–14. https://doi.org/10.1080/10635150701883881 Reeder, T.W., Townsend, T.M., Mulcahy, D.G., Noonan, B.P., Wood Jr., P.L., Sites Jr., J.W., Wiens, J.J., 2015. Integrated analyses resolve conflicts over squamate reptile phylogeny and reveal unexpected placements for fossil taxa. PLoS One 10, e0118199. Roig, F.A., Roig-Juñent, S., Corbalan, V., 2009. Biogeography of the Monte Desert. J. Arid Environ. 73, 164–172. http://doi.org/10.1016/j.jaridenv.2008.07.016 Ronquist, F., Cannatella, D., 1997. Dispersal-vicariance analysis: a new approach to the quantification of historical biogeography. Syst. Biol. 46, 195–203. Schulte, J.A., Valladares, J.P., Larson, A., 2003. Phylogenetic relationships within Iguanidae inferred using molecular and morphological data and a phylogenetic taxonomy of Iguanian lizards. Herpetologica 59, 399–419. Scolaro, J.A., Videla, F., Cei, J.M., 2003. Algunos modelos de especiación geográfica que interpretan aspectos de la diversidad herpetológica andino-patagónica. Hist. Nat. 2, 73–83. Sérsic, A.N., Cosacov, A., Cocucci, A.A., Johnson, L.A., Pozner, R., Avila, L.J., Sites Jr., J.W., Morando, M., 2011. Emerging phylogeographical patterns of plants and terrestrial vertebrates from Patagonia. Biol. J. Linn. Soc. 103, 475–494. https://doi.org/10.1111/j.1095-8312.2011.01656.x Sites, J.W., Marshall, J.C., 2003. Delimiting species: a Renaissance issue in systematic biology. Trends Ecol. Evol. 18, 462–470. https://doi.org/https://doi.org/10.1016/S0169-5347(03)00184-8 Stamatakis, A., 2014. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313. https://doi.org/10.1093/bioinformatics/btu033

Stern, C.R., Frey, F.A., Futa, K., Zartman, R.E., Peng, Z., Kyser, T.K., 1990. Trace element and Sr, Nd, Pb, and O isotopic composition of Pliocene and Quaternary alkali basalts of the Patagonian plateau lavas of southernmost South America. Contrib. to Mineral. Petrol. 104, 294–308. https://doi.org/10.1007/BF00321486 Sukumaran, J., Knowles, L.L., 2017. Multispecies coalescent delimits structure, not species. Proc. Natl. Acad. Sci. 114, 1607 LP – 1612. Tajima, F., 1983. Evolutionary relationship of DNA sequences in finite populations. Genetics 105, 437–460. Templeton, A.R., Crandall, K.A., Sing, C.F., 1992. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics 132, 619–633. Townsend, T.M., Mulcahy, D.G., Noonan, B.P., Sites, J.W., Kuczynski, C.A., Wiens, J.J., Reeder, T.W., 2011. Phylogeny of iguanian lizards inferred from 29 nuclear loci, and a comparison of concatenated and species-tree approaches for an ancient, rapid radiation. Mol. Phylogenet. Evol. 61, 363–380. https://doi.org/10.1016/j.ympev.2011.07.008 Victoriano, P.F., Coronado, T.M., Ortiz, J.C., 2010. A multivariate analysis of taxonomic limits in Diplolaemus Bell 1843. Gayana (Concepción) 74, 23–36. https://doi.org/10.4067/S071765382010000100006 Vizcaíno, S.F., Kay, R.F., Bargo, M.S., 2012. Early Miocene Paleobiology in Patagonia: High-latitude Paleocommunities of the Santa Cruz Formation, Cambridge University Press. Werneck, F.P., Gamble, T., Colli, G.R., Rodrigues, M.T., Sites, J.W., 2012. Deep diversification and longterm persistence in the South American “dry diagonal”: integrating continent-wide phylogeography and distribution modeling of geckos. Evolution 66, 3014–34. https://doi.org/10.1111/j.15585646.2012.01682.x Werner, F., 1898. Die Reptilien und Batrachier der Sammlung Plate, in: Fischer, G. (Ed.), Zoologische Jahrbücher. Abteilung Für Systematik, Geographie Und Biologie Der Tiere. Zoologische Jahrbücher, Jena, Germany, pp. 244–278. Whiting, A.S., Sites, J.W., Bauer, A.M., 2004. Molecular phylogenetics of Malagasy skinks (Squamata: Scincidae). African J. Herpetol. 53, 135–146. https://doi.org/10.1080/21564574.2004.9635506 Wiens, J.J., Hutter, C.R., Mulcahy, D.G., Noonan, B.P., Townsend, T.M., Sites Jr., J.W., Reeder, T.W., 2012. Resolving the phylogeny of lizards and snakes (Squamata) with extensive sampling of genes and species. Biol. Lett. 8, 1043 LP – 1046. Wiens, J.J., Reeder, T.W., Montes de Oca, A.N., 1999. Molecular phylogenetics and evolution of sexual dichromatism among populations of the Yarrow’s Spiny Lizard (Sceloporus jarrovii). Evolution (N. Y). 53, 1884–1897. https://doi.org/10.2307/2640448 Yang, Z., 2015. The BPP program for species tree estimation and species delimitation. Curr. Zool. 61, 854– 865. Yang, Z., Rannala, B., 2010. Bayesian species delimitation using multilocus sequence data. Proc. Natl. Acad. Sci. 107, 9264–9269. Yule, G.U., 1924. A mathematical theory of evolution, based on the conclusions of Dr. J. C. Willis, F.R.S. Philos. Trans. R. Soc. London. Ser. B, Contain. Pap. a Biol. Character 213, 21 LP – 87. Zheng, Y., Wiens, J.J., 2016. Combining phylogenomic and supermatrix approaches, and a time-calibrated

phylogeny for squamate reptiles (lizards and snakes) based on 52 genes and 4162 species. Mol. Phylogenet. Evol. 94, 537–547. https://doi.org/10.1016/j.ympev.2015.10.009

Graphical abstract

Highlights: 1- The southernmost species of Leiosaurae (all species of the genus Diplolaemus and Leiosaurus bellii) are poorly known and have a wide range of distributions in regions where a high degree of endemism has been recognized. 2- Our phylogenetic/phylogeographic analyses resolve five putative new species-level taxa in Diplolaemus and three in L. bellii, and the genetic structure in both groups corresponds to independently defined geographical areas. 3- We present updated, time-calibrated phylogenetic hypotheses. 4- Different post-Pliocene geological and paleoclimatic events in Patagonia, such as the “Great Patagonian Glaciation” and the “Coldest Pleistocene Glaciation”, are the most likely vicariant “drivers” of cladogenesis in these lizards.

Figure 1. Cytochrome b tree for Leiosaurae species and distribution of Diplolaemus (A) and L. bellii (B) lineages. The numbers in each node indicate posterior probabilities/boostrap values >70. Low values of PP or BT are indicated with ‘-‘. The scale bar of cytochrome b tree is in units of substitutions per site. Figure 2. Statistical parsimony networks based on cytochrome b haplotypes of Diplolaemus sexcinctus and D. leopardinus lineages. Number and abundance of each haplotype are indicated in networks and distributions, and locality numbers are indicated in maps. Haplotypes are colored per lineage and networks are separated based on a 95% parsimony criterion. Figure 3. Statistical parsimony networks based on cytochrome b haplotypes of D. darwinii and D. bibronii. Number and abundance of haplotypes are indicated in networks, and distributions and locality numbers are indicated in maps. Haplotypes are colored per lineage and networks are separated based on a 95% parsimony criterion. Figure 4. Statistical parsimony networks based on cytochrome b haplotypes of L. bellii. Number and abundance of haplotypes are indicated in networks and distributions and locality numbers are indicated in maps. Networks are separated based on a 95% parsimony criterion. Bottom right: Extended Bayesian skyline plot estimated with mitochondrial and nuclear DNA sequence data of L. bellii. The Y-axis indicates effective population size (Ne) based on an expected generation time of 10 years, while the X-axis indicates mean time in millions of years before present. Dotted lines indicate median effective population size, whereas thin grey polygons show the boundary of the 95% central posterior density (CPD) intervals. Figure 5. The Y-axis indicates effective population size (Ne) calculated using an expected generation time of 10 years, while the X-axis indicates time in millions of years before present. Dotted lines indicate median effective population size, whereas thin grey polygons show the boundary of the 95% central posterior density (CPD) intervals. Figure 6. Time-calibrated species tree of Leiosaurae obtained with mitochondrial and nuclear DNA sequences using *BEAST. Posterior probabilities (left, in italic) and age (right) in millions of years before present are separated by slash symbol in each node. Bars in each node indicate the 95% highest posterior densities (HPD) of divergence dates. Figure 7. BioGeoBEARS estimation of ancestral areas for Diplolaemus and L. bellii species. Letters in nodes indicate most likely ancestor areas estimated by the BAYAREALIKE+j model in Diplolaemus, and the DIVALIKE model in L. bellii. The numbers indicate age in millions of years before present. The map on the right shows the biogeographic areas delimited for analysis based on main phylogeographical patterns (breaks and refuges) in the area of distribution of species/lineages. SC: Santa Cruz Province, CH: Chubut Province, SomC: Somún Curá Plateau, RN: Rio Negro Province, NQ: Neuquén Province, MZA: Mendoza Province.

Table 1. Partition schemes and molecular evolution models for each partition defined by PartitionFinder. The best molecular evolution models and partition schemes were selected with the AICc and BIC criteria. For each analysis, the locus/loci used and the total length of fragments in base pairs are indicated. Partition schemes indicate the position of the codon separated by "," or "+" to indicate separation or non-separation respectively. The molecular evolution model and the number of sites for each partition are indicated in the last columns. Analysis

Locus

Fragmen t Length (bp)

Gene tree

Cyt-b

809

Partitio n Scheme

1,2,3

TVMef+G+I

Sites per parti tion 270

HKY+I

270

K81

269

Model

12S

1185

1+2+3

SYM+G

1185

ANL77

768

1+2+3

HKY+G

768

ANL90

588

1+2+3

TrN

588

ANL92

575

1+2+3

TVMef+G

575

CMOS

522

1+3,2

TrN

348

HKY

174

TVMef+I+G

270

Cyt-b

809

1,2,3

TrN+I+G

270

GTR+I+G

269

JC+I

320

BP&P/STACEY DMXL

960

1,2,3

EXPH5

889

1+3,2

NKTR

1062

1,2,3

RBMX

660

1+2+3

HKY+I

320

HKY+G

320

TrN+I

593

K81UF

296

TIM+I+G

354

HKY+G

354

HKY

354

HKY+I

660

12S

1185

1+2+3

GTR+G

1185

ANL77

768

1+2+3

TVM+G

768

ANL90

588

1+2+3

TrN+I

588

ANL92

575

1+2+3

TVMef+G

575

CMOS

522

1+2,3

Species Tree Cyt-b

DMXL EXPH5

809

960 889

1,2,3

1,2,3 1+3,2

K80

348

TrN+I

174

TVMef+I+G

270

HKY+I+G

270

GTR+I+G

269

JC

320

HKY+I

320

HKY+G

320

TrN+I

593

NKTR RBMX

1062 660

1,2,3 1+2+3

TVM

296

TIM+I+G

354

HKY+G

354

HKY

354

HKY+I

660

Table 2. Diversity genetic indexes and neutrality tests for species/lineages of Leiosaurae based on cytochrome b data. Diversity indexes Species /Lineages

n

h

Hd

Lineage A

21

5

0.424

Lineage B

82

19

0.818

Total

103

24

Lineage A

21

Total

S

Neutrality tests

θw ±S. D

π ± S.D

D

p-value

Fs

p-value

R2

p-value

0.01594±0.00573

0.01225±0.00435

-0.9194

0.1798

8.3987

0.9958

0.1158

0.3846

71

0.01778±0.00487

0.01210±0.00148

-1.6250

0.0264

-4.4066

0.0726

0.0435

0.0106*

0.847

73

0.01748±0.00462

0.0108±0.00134

-1.2384

0.0852

-0.4120

0.5124

0.0591

0.0994

11

0.929

15

0.00515 ±0.00133

0.00429±0.00060

-0.6144

0.3082

-0.98843

0.0702

0.1081

0.2414

23

13

0.941

68

0.02277±0.00783

0.01186±0.00459

-1.8966

0.0140*

-0.2990

0.4518

0.0896

0.1052

6

4

0.800

19

0.01039±0.00530

0.00857±0.00344

-1.0896

0.1606

1.7979

0.7938

0.2478

0.8282

D. sexcintus sensu stricto

42

23

0.958

62

0.01797±0.00555

0.01323±0.00088

-0.9430

0.1844

-3.1630

0.1790

0.1238

0.9644

D. s. Lineage 1

6

6

1

29

0.01584±0.00786

0.01463±0.00253

0.4820

0.3838

0.7178

0.2044

0.1153

0.0350*

D. s. Lineage 2

22

9

0.861

56

0.01915±0.00673

0.01616±0.00324

-0.6195

0.2966

4.2396

0.9464

0.1084

0.2988

D. s. Lineage 3

11

6

0.855

18

0.00766±0.00343

0.00753±0.00114

-0.0798

0.5138

1.2072

0.7258

0.1517

0.4590

D. s. Lineage 4

8

3

0.464

10

0.00481±0.00248

0.00356±0.00208

-1.2830

0.1168

2.4141

0.8974

0.2680

0.9726

D. s. Lineage 5

15

7

0.819

30

0.01150±0.00458

0.01135±0.00186

-0.0554

0.5330

2.7039

0.8906

0.0881

0.4216

104

54

0.980

175

0.04175±0.01039

0.04163±0.00162

-0.0096

0.5746

-2.0426

0.3714

0.0931

0.4365

Lineage A

44

25

0.925

53

0.01506±0.00467

0.01603±0.00090

0.2272

0.6690

-2.8970

0.1934

0.1153

0.6404

Lineage B

14

13

0.989

36

0.01399±0.00557

0.01156±0.00082

-0.7529

0.2378

-4.5230

0.0208*

0.0946

0.0430*

Lineage C

3

3

1

4

0.00330 ±0.00237

0.03300 ±0.00123

-

-

-

-

0.3118

0.5780

47

28

0.934

89

0.02491±0.00736

0.02170±0.00298

-0.4593

0.3616

-2.3360

0.2464

0.9419

0.3784

D. bibronii 46

D. darwinii

D. leopardinus D. sexcinctus sensu Victoriano et al. 2010

Total L. bellii

Total

*p-values <0.05 and their respective statistic values are shown in black, and p-values <0.1 are shown in italic. n: number of sequences, h: haplotype number, Hd: haplotype diversity, S: polymorphic sites, θw: Watterson estimator, π: nucleotide diversity, Fs: Fu’s Fs, D: Tajima’s D, R2: Ramos-Onsins and Rozas’s R2, s.d.: standard deviation.

Table 3. Pairwise genetic distances of cytochrome-b for the species/lineages of Diplolaemus (A) and L. bellii (B). Intra-group distances are shown in the diagonal, under and above it, corrected (equation 10.21 of Nei, 1987) and uncorrected distances respectively. All distances are expressed in percentages. Values greater than 3% are shown in bold. A)

D. s. Lineage 1

D. s. Lineage 2

D. s. Lineage 3

D. s. Lineage 4

D. s. Lineage 5

1.60

3.97

6.18

6.45

6.25

6.20

6.77

2.44

1.45

6.13

6.95

6.24

6.51

6.17

5.00

5.03

0.75

4.57

3.91

4.12

5.40

5.47

6.05

4.02

0.35

4.31

3.94

5.67

4.89

4.95

2.97

3.58

1.13

2.59

5.38

4.75 5.43

5.13 4.91

3.09 4.49

3.11 4.96

1.37 4.28

1.31 4.20

5.39 1.07

6.72 14.88

7.13 14.71

6.17 15.80

7.93 16.82

6.96 16.42

6.69 15.88

5.99 16.27

D. s Lineage 1 D. s. Lineage 2 D. s. Lineage 3 D. s. Lineage 4 D. s. Lineage 5 D. sexcinctus sensu stricto D. bibronii D. leopardinus D. darwinii

D. sexcinctus sensu stricto D. bibronii

B)

L. bellii Lineage A L. bellii Lineage B L. bellii Lineage C

L. bellii Lineage A 1.60 3.96 5.31

L. bellii Lineage B 5.34 0.32 4.78

L. bellii Lineage C 6.27 5.53 1.15

D. le

Table 4. Comparison between biogeographic models for ancestral area reconstructions of Diplolaemus (A) and L. bellii (B). The models compared were Dispersal–Extinction– Cladogenesis (DEC), Dispersal–Vicariance Analysis (DIVA) and BayArea, as well as these three models allowing for founder-event speciations (+J). The LnL, AICc and AICc weights are shown for each model. A) Num. LnL Params. d e j AICc AICc_wt DEC -48.14 2 0.082 0.009 0.000 102.00 0.032 DEC+J -44.40 3 0.058 0.000 0.905 98.80 0.157 DIVALIKE -48.44 2 0.163 0.676 0.000 102.58 0.024 DIVALIKE+J -45.04 3 0.070 0.045 0.341 100.07 0.083 BAYAREALIKE -50.43 2 0.125 0.885 0.000 106.57 0.003 BAYAREALIKE +J -42.91 3 0.045 0.095 0.193 95.81 0.701 B) LnL -10.9 -10.9 -10.54 -10.54 -11.88 -11.88

Num. Params. 2 3 2 3 2 3

d 0.43 0.43 0.43 0.43 0.76 0.77

e j 6.90E-08 0 1.00E-12 1.00E-05 1.00E-12 0 1.00E-12 1.00E-05 2.85 0 2.9 1.00E-05

AICc 25.8 27.8 25.09 27.09 27.77 29.77

DEC DEC+J DIVALIKE DIVALIKE+J BAYAREALIKE BAYAREALIKE+ J LnL: log-likelihood, d: rate of dispersal, e: rate of extinction, j: relative probability of founder-event speciation at cladogenesis, AIC: Akaike's information criterion, AICc: Akaike's information criterion corrected for small sample size.

AICc_wt 0.26 0.096 0.37 0.14 0.098 0.036

Author statement: Martin Femenias: Conceptualization, Methodology, Software, Formal analysis, Investigation, Data Curation, Writing - Original Draft, Visualization. Jack W. Sites Jr.: Writing - Review & Editing, Luciano J. Avila: Investigation. Mariana Morando: Data Curation, Writing - Review & Editing, Supervision, Project administration.