Integrative taxonomy and preliminary assessment of species limits in the Liolaemus walkeri complex (Squamata, Liolaemidae) with descriptions of three new species from Peru

Abstract Species delimitation studies based on integrative taxonomic approaches have received considerable attention in the last few years, and have provided the strongest hypotheses of species boundaries. We used three lines of evidence (molecular, morphological, and niche envelopes) to test for species boundaries in Peruvian populations of the Liolaemus walkeri complex. Our results show that different lines of evidence and analyses are congruent in different combinations, for unambiguous delimitation of three lineages that were “hidden” within known species, and now deserve species status. Our phylogenetic analysis shows that L. walkeri, L. tacnae and the three new species are strongly separated from other species assigned to the alticolor-bibronii group. Few conventional morphological characters distinguish the new species from closely related taxa and this highlights the need to integrate other sources of data to erect strong hypothesis of species limits. A taxonomic key for known Peruvian species of the subgenus Lioalemus is provided.


Introduction
The issue of species delimitation (building explicit hypotheses about species lineages and their geographic boundaries) has received considerable attention in the last decade due in part to an emerging consensus about species concepts and new approaches for testing species boundaries (Sites and Marshall 2003, de Queiroz 2007, Knowles and Carstens 2007, Wiens 2007, Padial and De la Riva 2010, Padial et al. 2010, Hart 2011, Zapata and Jiménez 2012, Camargo and Sites 2013. The ontological General Lineage Concept (GLC) defines a species as a group of separately evolving meta-population lineages, originally proposed by Mayden (1997Mayden ( , 2002 and de Queiroz (1998de Queiroz ( , 2005. This definition is generally supported by a consensus view in evolutionary biology (Padial and De la Riva 2010, Padial et al. 2010, Hart 2011, Zapata and Jiménez 2012, but see Hausdorf 2011). The GLC distinguishes the primary property (species are separately evolving meta-population lineages) that is shared by most previous competing species concepts (e.g., biological, phylogenetic, ecological species concept, etc.), from secondary properties (e.g., reproductive isolation, character fixation, niche differentiation, etc.) that arise at different times during the processes of speciation (de Queiroz 2007). These secondary properties are lines of evidence that are relevant to inferring the species boundaries (de Queiroz 2005(de Queiroz , 2007.
In addition to this agreement with respect to GLC, there is a growing number of new empirical methods of species delimitation (SDL; Pons et al. 2006, Knowles and Carstens 2007, Kubatko et al. 2009, Carstens and Dewey 2010, Flot et al. 2010, Hausdorf and Hennig 2010, Martínez-Gordillo et al. 2010, Gurgel-Gonçalves et al. 2011. These new methods for testing hypotheses of species boundaries have been accommodated under the new term "integrative taxonomy" (IT; Dayrat 2005, Padial and De la Riva 2010, Padial et al. 2010. Methods such as the multi-locus coalescent to infer species limits without monophyletic lineages, ecological niche modeling (ENM) to assess spatial distributions of closely related species, and multivariate tolerance regions to test for discontinuities or gaps in morphology, have all been used in new integrative taxonomic studies (Omland et al. 2006, Knowles and Carstens 2007, Raxworthy et al. 2007, Rissler and Apodaca 2007, Vasconcelos et al. 2012, Zapata and Jiménez 2012. Character fixation as well as discontinuities or gaps have been used as a SDL criterion to assess species limits based on genetic and morphological characters (Marshall et al. 2006, Zapata andJiménez 2012). Fixed differences and gaps in morphology suggest that some evolutionary force (e.g., absence of gene flow, natural selection) prevent two putative taxa from homogenizing (Wiens andServedio 2000, Zapata andJiménez 2012). Often analysis of variance or discriminant analysis have been used to evaluate morphological differentiation in SDL studies, but these statistics, even if significant, evaluate central tendencies and not gaps in morphology, and the latter may be more relevant for testing species boundaries (Zapata and Jiménez 2012). In addition to character fixation and gaps in morphology, niche envelopes can be used to assess the status of uncertain populations which are separated from closely related species by areas that are outside of the climatic niche envelope, and where gene flow between these species is unlikely because it would involve crossing unsuitable habitat (Wiens and Graham 2005). Ecological niche modeling (ENM) can summarize niche envelopes and this approach has also been used in SDL studies (e.g., Raxworthy et al. 2007, Rissler andApodaca 2007).
Well-supported hypotheses of species boundaries are essential because species are used as basic units of analysis in several areas of biogeography, ecology, and macroevolution, and from the broader perspective of evolutionary theory, delimiting species is important in the context of understanding many evolutionary mechanisms and processes (Sites and Marshall 2003, Wiens 2007. Among animal groups, lizards have been used extensively in evolutionary studies ranging from community ecology, behavioral ecology, multiple origins of body elongation coupled with limb reduction/loss, multiple origins of novel reproductive modes, including parthenogenesis and viviparity , as well as phylogeography and speciation studies (Camargo et al. 2010).
SDL studies in lizards have included molecular markers, morphological characters and/or models of species distributions (Camargo et al. 2010). In particular, several clades of the genus Liolaemus Wiegmann, 1834 have been studied intensively using molecular and morphological data to delimit species and infer phylogeographic histories (Morando et al. 2008, Victoriano et al. 2008, Breitman et al. 2011a, 2012, and for testing hypotheses about evolutionary processes (Olave et al. 2011) and performance (in accuracy and precision) of different SDL methods (Camargo et al. 2012). This South American genus includes ~ 230 species (Breitman et al. 2011b), and extends from central Peru to Tierra del Fuego, and from sea level on both Atlantic and Pacific coasts to almost 5000 m in elevation. Species diversity is highest in the Andes and adjacent arid regions, and new species descriptions are published at a rate of 4-5/ yr, from moderately well-known areas in Argentina and Chile.
In most cases these studies have demonstrated that populations assigned to single species based on generalized morphological features and limited field sampling, tend to under-represent biodiversity. Distinct lineages have been revealed by molecular data, many of which are later described as new species (e.g., Breitman et al. 2011a, b). The largest poorly-known areas for the genus are the Andean regions of Bolivia, Peru and northern Chile. Intensive fieldwork and molecular phylogenetic studies have never been systematically carried out in these regions, and species descriptions have traditionally been based on gross comparisons of morphological characters from small sample sizes and limited geographic sampling. So SDL studies are needed in the extreme northern range of Liolaemus (e.g., Peru) based on intensive geographic sampling and large series for collection of new molecular, coloration, and various classes of morphological data.
Currently, 14 species of Liolaemus are known from Peru (L. montanus group, 10 spp; L. alticolor group, 4 spp), but SDL studies based on an integrative approach have not been carried out in either of these groups. Moreover, several areas in the Peruvian Andes remain completely unexplored, and based on recent studies in the southern range of Liolaemus, it is highly probable that the Peruvian Andes harbor many undiscovered species. Here, we use new molecular, morphological, and geographic data from known Peruvian species (L. alticolor Barbour, 1909, L. incaicus Lobo, Quinteros & Gómez, 2007, L. tacnae (Shreve, 1941 and L. walkeri Shreve, 1938), assigned to the L. alticolor group, and three populations morphologically similar to L. walkeri (identified by their regions of occurrence: Ancash, Ayacucho and Cusco), to present the first SDL study based on an IT approach. Our results provide evidence that three new lineages deserve species status, and these are described herein.

Sampling and DNA extraction
Lizards were collected by hand, photographed and sacrificed with an injection of pentobarbital. After liver tissue was collected for DNA samples, whole specimens were fixed in formaldehyde at 10% and transferred to 70% ethanol for permanent storage in museum collections. Tissue samples were collected in duplicate, stored in 96% ethanol and deposited at the Bean Life Science museum at Brigham Young University (BYU) and Museo de Historia Natural de San Marcos (MUSM) (see Data resources below). Total genomic DNA is extracted from liver/muscle tissue following the protocol of Fetzner (1999), and using Qiagen DNeasy kits (Qiagen, Inc., Valencia, CA).

Mitochondrial DNA amplification and sequencing
Forty-eight samples from 40 localities were sequenced for 669bp of the mtDNA cytochrome b (cyt-b) region, with LIO742F 5'-TCGACCTVCCYGCCCCATCA-3' and LIO742R 5'-GAGGGGTTACTAAGGGGTTGGC-3' primers (this study), and all unique cyt-b haplotypes were sequenced for a 12S region (752 bp) using primers 12Stphe 5'AAAGCACRGCACTGAAGATGC-3' and 12SE 5'-GTRCGCTTAC-CWTGTTACGACT-3' (Wiens et al. 1999). Double stranded polymerase chain reactions (PCR) were amplified under the following conditions: 1.0 μL of genomic DNA, 1.0 μL of light strand primer 1.0 of μL of heavy strand primer, 1.0 μL of dinucleotide pairs, 2.0 μL of 5x~ buffer, 1.0 μL of MgCl 10x~ buffer, 0.18 μL of Taq polymerase, and 7.5 μL of diH2O. PCR amplification was executed under the following conditions: initial denaturation at 95°C for 2 min, followed by a second denaturation at 95°C for 35 s, annealing at 52°C for 35 s, followed by a cycle extension at 72°C for 35 s, for 31 cycles. PCR products were visualized on a 10% agarose gels to ensure the targeted products were cleanly amplified, and then purified using a MultiScreen PCR (mu) 96 (Millipore Corp., Billerica, MA) and directly sequenced using the BigDye Terminator v 3.1 Cycle Sequencing Ready Reaction (Applied Biosystems, Foster City, CA). The cycle sequencing reactions were purified using Sephadex G-50 Fine (GE Healthcare) and MultiScreen HV plates (Millipore Corp.). Samples were then analyzed on a ABI3730xl DNA Analyzer in the BYU DNA Sequencing Center.

Phylogenetic reconstruction
All sequences were aligned in MUSCLE (Edgar 2004) plugin, and cyt-b sequences were translated to check for premature stop codons in GENEIOUS®PRO v5.6.6. Cyt b haplotype diversity was estimated using DnaSP (Librado and Rozas 2009), and concatenated cyt-b and 12S regions were edited using GENEIOUS®PRO (Drummond et al. 2011). For ingroups and outgroups we used selected species of the subgenus Liolaemus that are assigned to different species groups and for which cyt-b and 12S sequences are available in GenBank. Our ingroup samples included taxa that have been assigned to the same species group as L. tacnae and L. walkeri (alticolor-bibronii group), including: L. abdalai Quinteros, L. bibronii (Bell), L. gracilis (Bell), L. ramirezae Lobo & Espinoza and L. saxatilis Ávila & Cei (Lobo et al. 2010). To further test for monophyly of the alticolor-bibroni group, we sampled three species assigned to different species groups (robertmertensi, pictus and monticola groups), but nested within the subgenus Liolaemus; these include: L. monticola Müller & Hellmich, L. pictus Duméril & Bribon and L. robertmertensi Hellmich (Lobo et al. 2010). We used L. lineomaculatus Boulenguer, a species belonging to the subgenus Eulaemus (Lobo et al. 2010, Fontanella et al. 2012) as the outgroup. All new sequences were deposited in GenBank (accession numbers KF923633-KF923660 and KF923661-KF923688 for cyt-b and 12s respectively) and a list of all haplotypes, GenBank accession and museum voucher numbers used for the phylogenetic analysis are provided as Supplementary file 1.
Bayesian Information Criteria in JMODELTEST (v 0.01; Posada 2005) identified the best-fit model of evolution for the complete data set of haplotypes as TPM2+ I + Γ. A Maximum-likelihood (ML) search in PHYML (Guindon and Gascuel 2003) was performed with 1000 replicates for bootstrap analyses; we consider strong nodal support for bootstrap values ≥ 70 (Hillis and Bull 1993;with caveats). Because the TPM2+ I + G model is not incorporated in the MRBAYES (Huelsenbeck and Ronquist 2001) plugin of GENEIOUS®PRO v5.6.6, we used a model with the closest likelihood available (GTR + I + Γ). Two parallel runs were performed in MRBAYES using four chains (one cold and three hot) for 1.1 × 10 6 generations and sampling every 200 generations from the Markov Chain Monte Carlo (MCMC). We determined stationarity by plotting the log likelihood scores of sample points against generation time; when the values reached a stable equilibrium and split frequencies fall below 0.01, stationarity was assumed. We discarded 100,000 samples and 10% of the trees as burn-in. A maximum clade credibility (MCC) tree was constructed using TREEANNOTATOR v1.7.5 (Drummond et al. 2012). We consider Bayesian Posterior Probabilities (BPP) >95% as evidence of significant support for a clade (Huelsenbeck andRonquist 2001, Wilcox et al. 2002).

Morphological data and analyses
A total of 199 individuals (see species descriptions and Data resources below) representing three putative different populations and four Peruvian species (L. alticolor, L. incaicus, L. tacnae and L. walkeri) assigned to the L. alticolor group were scored for three classes of morphological characters. We performed a character analysis of 17 discrete binomial characters related to scalation, pattern of coloration and skin folds, including the following: presence/absence of smooth (1) temporal scales and (2) dorsal head scales, contact or not of (3) rostral to nasal scale, presence/absence of (4) mucronate dorsal scales and (5) precloacal pores, (6) preocular scale same or different color as loreal region, presence/absence of (7) spots on dorsal head scales, (8) black line surrounding the interparietal scale, regular spots or marks in (9) paravertebral field and (10) lateral field, presence/absence of dorsolateral stripes (11) and vertebral line (12), marks or spots on throat (13), melanistic belly (14), ringed pattern in ventral tail (15), and presence/absence of antehumeral (16) and neck folds (17). All characters were scored using a stereomicroscope and from photos of live animals taken in the field.
For statistical analyses of these discrete variables we used tolerance intervals as described in the tolerance package of Young (2010), which in a random sample of a univariate population, is an interval expected to contain a specified proportion or more of the sampled population (Krishnamoorthy and Mathew 2009). We used binomial tolerance intervals to estimate the number of individuals that comprise 95% of the population expected to have one state with a 0.05 level of significance (following Wiens andServedio 2000, Zapata andJiménez 2012). One-sided binomial tolerance intervals were estimated using the Wilson method (WS), which is appropriate when the sample sizes are small (n ≤ 40) (Young 2010).
We scored the following 11 morphometric characters: (SVL) snout-vent length, (AGL) axilla-groin length (between the posterior insertion of forelimb and anterior insertion of thigh), (HL) head length (from snout to anterior border of auditory meatus), (HW), head width (at widest point), (FOL) forelimb length (distance from the attachment of the limb to the body to the terminus of the fourth digit), (HIL) hindlimb length (distance from the attachment of the limb to the body, to the terminus of the fourth digit), (SL) snout length (from snout to anterior border of eye), (AMW) auditory meatus width, (AMH) auditory meatus height, (RW) rostral width, and (RL) rostral length. We also scored five meristic characters, including: (MBS) number of midbody scales (counted transversely at the middle of the body), (DTS) dorsal trunk scales (counted from the level of anterior border of ear to anterior border of thighs), (DHS) dorsal head scales (counted from the rostral scale to anterior border of ear), (VS) ventral scales (counted from the mental scales to the cloaca), and (SCI) number of scales in contact with the interparietal. Measurements and counts were taken from the right side of the animal using a stereomicroscope. Morphometric data were only taken for adult males and females.
After testing for normality in all morphometric and meristic characters with the Shapiro-Wilks test (Shapiro and Francia 1972), we summarized means and ranges for all population samples, and performed Principal Component Analyses (PCA) and Correspondence Analyses (CA) separately for each class of characters and by sex, to summarize patterns. Results of PCA and CA were then compared with the analysis of continuous characters by estimating normal tolerance intervals to find gaps or discontinuities in each class of morphological characters. We used normal tolerance intervals to estimate the lowest and highest values of a continuous character that is contained in 95% of the population with a 0.05 level of significance. Two-sided normal tolerance intervals were estimated using the Howe method (HE), which is considered to be extremely accurate, even for small sample sizes (Young 2010).
For comparison with normal tolerance intervals we assessed the morphometric and meristic characters with univariate ANOVA and Mann-Whitney U tests for parametric and non-parametric distributions, respectively. When the assumption of equal-variance was not met for an ANOVA test, the unequal-variance (Welch) version of ANOVA was performed. Each character was tested for intersexual differences, and if present, the sexes were analyzed separately. Results were considered significant when p ≤ 0.05. However, we didn't use the results of the ANOVA and Mann-Whitney U tests in our taxonomic decisions (see Introduction and Discussion). Binomial and tolerance intervals were calculated with the package Tolerance (Young 2010) in R v3.0.1 (R core team 2013). Test of normality, PCA, CA and univariate tests were performed using PASTv. 2.08b, (Hammer et al. 2001).

Distributional models
We used the maximum entropy model implemented in the program MAXENT v3.3.3e ) to predict where the Peruvian lineages of L. walkeri complex are most likely to occur under current climatic conditions. MAXENT generates distributional models (or ecological niche models; ENMs) using presence-only records, contrasting them with background/pseudoabsence data sampled from the remainder of the study area. We chose this approach because of its overall better performance with presence-only data and with small sample sizes (Elith et al. 2006). ENMs were developed from occurrence points used in this study, and records without duplicates are: 22 for Ancash, 31 for Ayacucho, 16 for Cusco, 33 for L. tacnae and 52 for L. walkeri (see Data resources below). For niche predictions, we used the 19-bioclimatic variables from the WorldClim v1.4 dataset with a resolution of 2.5 min (Hijmans et al. 2005). Bioclimatic variables were derived from monthly temperature and precipitation layers, and represents biologically meaningful properties of climate variation (Hijmans et al. 2005). Layers were trimmed to the areas surrounding each species and populations that might represent new species, and then projected over a larger region (-9.828° to -17.839° and -77.486° to -69.811°).
For model calibration we used the default settings with 1000 iterations, and the minimum training value averaged over the 10 replicates as threshold with the default convergence threshold (10 -5 ). Due to our smaller samples sizes, we used for model calibration and performance the cross-validation option with 10 replicates, and average the results to estimate species niche and distributions. For model significance, 25% of localities were randomly set aside as test points and the area under the curve (AUC) was calculated, which summarizes the model's ability to rank presence localities higher than a sample of random pixels (Peterson et al. 2011). AUC values ≤ 0.5 correspond to predictions that are equal or worse than random. AUC values > 0.5 are generally classed into (1) poor predictions (0.5 to 0.7); (2) reasonable predictions (0.7 to 0.9); and (3) very good predictions (>0.90; but see Peterson et al. 2011, for caveats on use of AUC in presence/background data). Model clamping was checked with the "fade by clamping" option available in MAXENT v 3.3.3e. Estimates of bioclimatic variable importance was performed using the Jackknife test. We used the logistic output (probability values) and mapped the distributional models showing areas from the average minimum logistic values (threshold) to 1 as areas suitable for species.
Schoener's D metric was introduced as a measure of niche similarity between pairs of populations (or species) by Warren et al. (2010), and is calculated using the ENMTOOLS package. We calculated these values by comparing the climatic suitability of each grid cell in the projected area obtained with MAXENT. This similarity measure ranges from 0 (niche models have no overlap) to 1 (niche models identical; Warren et al. 2008). We estimated similarity measures and then tested whether the ENMs produced by two populations or species are identical using the niche identity test in ENMTOOLS. One hundred pseudoreplicate data sets were generated to obtain a distribution of D scores, and we reject the hypothesis of niche identity when the empirically observed value for D is significantly lower than the values expected from the pseudoreplicated data set (Warren et al. 2010).

Data resources
The data underpinning the analysis reported in this paper are deposited in the Dryad Data Repository at http://doi.org/10.5061/dryad.0q7pc, and at GBIF, the Global Biodiversity Information Facility, http://ipt.pensoft.net/ipt/resource.do?r=ocurrence_re-cords_liolaemus_walkeri_complex.

Phylogenetic Analysis
A tree with maximum likelihood bootstrap values (logL = -8452.31415, MLB) and Bayesian posterior probabilities (BPP) based on 1421 aligned base pairs is shown in Fig. 1. Differences between both methods are mentioned below. Both ML and Bayesian analyses recovered Ancash, Ayacucho, Cusco, L. tacnae and L. walkeri haplotypes as monophyletic groups with high support. Both also showed a close relationship between Ayacucho and L. walkeri haplotypes, but relationships between Ancash, Cusco and the (L. walkeri + Ayacucho) clade were unresolved and with moderate support in the ML tree (MLB 65%). The Bayesian analysis recovers Ancash as the sister to the (L. walkeri + Ayacucho) clade with low support (BPP 0.5), and Cusco as the sister clade to the ((L. walkeri + Ayacucho) Ancash) clade with moderate support (BPP 0.9). In both analyses, Liolaemus tacnae is recovered as the sister group of the (Ancash + Cusco + (L. walkeri +Ayacucho)) clade with moderate support ( MLB 65%, BPP 0.9). Liolaemus tacnae and L. walkeri are assigned to the alticolor-bibronii group, but the clade (L. tacnae (Ancash + Cusco + (L. walkeri + Ayacucho))) is strongly differentiated from the other species assigned to the alticolor-bibronii group (Fig. 1).
The monophyletic group (L. tacnae (Ancash + Cusco + (L. walkeri + Ayacucho))) is the sister group of a clade comprised of taxa assigned to different species groups in the subgenus Liolaemus, including species of the alticolor-bibronii group. The relationships of these two more inclusive clades showed high MLB, but low BPP values. In this clade, both ML and Bayesian analyses recovered L. alticolor and L. incaicus haplotypes as monophyletic groups with high support. In our ML analysis, the clade (L. alticolor + L. incaicus) has unresolved relationships with L. ramirezae and the clade (L. robertmertensi + (L. gracilis + L. saxatilis)), and this latter clade has high BPP but low MLB support (Fig. 1). Liolaemus abdalai and L. bibronii are recovered as sister taxa with high support, and this clade is sister to the clade (L. ramirezae + (L. incaicus + L. alticolor) + (L. robertmertensi + (L. gracilis + L. saxatilis))) also with high support (Fig. 1). Liolaemus pictus is sister to the clade ((Liolaemus abdalai and L. bibronii) + (L. ramirezae + (L. incaicus + L. alticolor) + (L. robertmertensi + (L. gracilis + L. saxatilis)))), and L. monticola is basal to a clade that includes L. pictus and its sister group.

Binomial discrete characters
Because our phylogenetic analysis did not show a close relationship between (L. alticolor + L. incaicus) and the (L. tacnae (Ancash + Cusco + (L. walkeri + Ayacucho))) clades, we focus our comparisons on these last five taxa. Of the 17 binomial characters, four were useful for species delimitation among these taxa (Table 1). One-sided binomial tolerance intervals (BTI) for 95% of the population with a 0.05 level of significance is indicated below for each of these four characters.
Ancash (n = 12) and L. tacnae (n = 18) males differed from Ayacucho, Cusco and L. walkeri males in lacking precloacal pores ( Fig. 2A and D; vs. presence in panels B, C, and E). Although these differences are fixed in our samples, the BTI tests showed table 1. Binomial characters for females (F) and males (M) of focal populations of Liolaemus lizards sampled for this study. Character states useful for species discrimination are in bold, and states only assessed on adults are indicated with an asterisk.

Morphometric and meristic characters
Our empirical results are summarized in Table 2, and tolerance intervals are given in Tables 3 and 4 for morphometric and meristic variables, respectively. Statistical tests rejected normality for HW, AMW, RW and all meristic characters, but we assumed normality because our sample sizes were too small to implement non-parametric tolerance interval tests. We did not find any diagnostic character or gaps in either data set (Tables 3 and 4).
Principal Component and Correspondence Analyses separated by sex or pooled together did not show any differences, so we present the results of the pooled analyses. Principal Component (PC) Analysis revealed that PC1 and PC2 explained 90% of the variance, and the Correspondence Analysis revealed that Correspondence Axis (CA) 1 and CA2 explained 66% of the similarity for morphometric and meristic data, respectively (see also Supplementary file 4 for corresponding eigenvalues, and percentages of variance and similarity accounted by principal components and correspondence axes). The bivari- ate plot for the morphometric variables revealed extensive overlap of L. walkeri with the remaining four samples, but minimal overlap between the Ancash and Cusco samples, and little overlap between the Ayacucho and Cusco (Fig. 5A). Both of these pairs are differentiated primarily along PC1, for which SVL and AGL contributed the highest loadings (0.85 and 0.47 respectively). The Cusco samples are characterized by shorter SVL and axilla-groin lengths than the Ancash and Ayacucho samples. The PC analyses revealed extensive overlap among all samples along PC2, and the CA for the meristic variables (Fig. 5B) revealed extensive overlap among all five samples along both axes.
Only significant results of ANOVA are mentioned below and the sex of a particular species or population is indicated only if significantly different from the opposite sex. For SVL, there were significant differences between Ancash vs. Cusco, L. tacnae and L. walkeri; Ayacucho vs. Cusco and L. tacnae; Cusco vs. L. tacnae and L. walkeri.
Only significant results of Mann-Whitney U are mentioned below and the sex of a particular species or population is indicated only if significantly different from the opposite sex. For HW, there were significant differences between Ancash males vs. Ayacucho, Cusco, L. tacnae and L. walkeri; Ancash females vs. Cusco and L. tacnae; Ayacucho vs. Cusco and L. tacnae; Cusco vs. L. walkeri.

Distributional models
The predicted distribution in all cases matched the known range of each taxon, although some of these overlap. However, the distributional models of Ayacucho vs L. tacnae ( Fig. 6; C vs. E), as well as those for L. walkeri and L. tacnae ( Fig. 6; E vs. F) are virtually mutually exclusive. All other combinations of distributional models overlapped, but differed in the contribution of bioclimatic variables to each niche envelope, and in predicting the known distribution of particular taxa (Table 5, Fig. 6). For example, the most important bioclimatic variables for the Ancash model were completely different from those for the L. walkeri and Ayacucho models (Table 5). In the same manner, the most important bioclimatic variables contributing to the Ayacucho model were completely different from those for the L. walkeri and Cusco models ( Table 5). The most important bioclimatic variables for the Cusco model were completely different to those for L. tacnae (Table 5). Moreover results from the Niche Identity Test found all pairwise comparison between focal populations and species significantly different, except for Ancash and Cusco ( Table 6).
The Ancash model (Fig. 6B) overlapped the known geographic distributions of Ayacucho, Cusco, L. tacnae, and partially with L. walkeri, but the two most important bioclimatic variables accounting for 94.3% of the contribution to this model were Precipitation of Warmest Quarter (63.3%) and Isothermality (31.0%; Table 5). These were also the most important variables in the permutation and jackknife tests. Thus the Ancash samples are characterized by a niche envelope with relative lower precipitation and more variation in annual temperature. The AUC score for this model = 0.87 (± 0.05), suggesting that the model prediction was reasonable (Fig. 7A).
The Ayacucho model did not overlap known distributions of Ancash, Cusco, L. tacnae, and only partially overlapped L. walkeri (Fig. 6C); the two most important bioclimatic variables accounting for 75.1% of the contribution to this model were Precipitation of Driest Quarter (64.4%) and Maximum Temperature of Warmest Period (10.7%; Table 5). In the permutation and jackknife tests, Precipitation of Driest Quarter was also the most important variable. In other words, the Ayacucho samples are characterized by a relatively wet and warm niche envelope, and the AUC score = 0.76 (± 0.06), suggesting that model prediction was reasonable (Fig. 7B).
The Cusco model did not overlap the known distribution of Ancash, overlapped most of Ayacucho and L. walkeri, and overlapped some of L. tacnae (Fig. 6D). The  (Table 5). In the permutation and jackknife tests, Precipitation of the Wettest Period was also the most important variable indicating a niche envelope with relative more precipitation in the wettest period of the year. The AUC score = 0.91 (± 0.03), suggesting that model prediction was reasonable (Fig. 7C).  The L. tacnae model did not overlap the known distributions of any of the remaining taxa (Fig. 6E); the three most important bioclimatic variables accounting for 64.2% of the contribution to the model are Precipitation of Warmest Quarter, Precipitation of Driest Quarter, and Precipitation of Wettest Quarter (Table 5). In the permutation test, the most important variable was Precipitation of the Coldest Quarter, but in the jackknife tests Precipitation of Warmest Quarter, Precipitation of Wettest Quarter and Annual Precipitation were the most important variables. This indicates that L. tacnae samples are characterized by a drier niche envelope relative to all other populations and the AUC score (0.85 ± 0.06) suggests that this model prediction was reasonable (Fig. 7D).
The L. walkeri model overlaps the known distribution of the Ancash and partially that of the Cusco samples (Fig. 6F); the two most important bioclimatic variables accounting for 84.2% of the contribution to the model are Precipitation of Driest Period and Precipitation of Wettest Period (Table 5). In the permutation and jackknife tests, Precipitation of Wettest Period also was the most important variable. This suggests a relative wetter niche envelope relative to all other populations, and the AUC score (0.77 ± 0.08) suggests that the model prediction was reasonable (Fig. 7E).
The niche identity test results showed that observed values of Schoener's D between all populations and species were significantly lower than null distribution of pseudoreplicates except for Ancash and Cusco (Table 6).

Integrative taxonomy
Results of mitochondrial haplotypes, binary (presence/absence of precloacal pores, spots or regular marks in lateral field, melanistic belly in adult males, ringed ventral tail pattern), morphometric (snout-vent length, axila-groin length and hindlimb length) characters and niche identity tests in various combinations, differentiated Ancash, Ayacucho and Cusco samples from each other, and from L. tacnae and L. walkeri. Despite the fact that binomial tolerance intervals showed the possible presence of polymorphisms even at a frequency cut off of 0.5% in discrete characters, we hypothesize that increasing samples sizes will lower the hypothesized frequencies of the alternative states for each taxon. Normal tolerance intervals and distributional models showed overlap between all paired combinations of samples except for the Ayacucho vs. L. tacnae distributional models and niche identity tests showed statistical differences between all pairwise comparisons but Ancash vs. Cusco. Note that this is an extremely conservative approach; if we simply look at the data and count the number of "fixed" differences between all combinations of samples, we would conclude that the following pairs are unambiguously diagnosed: Ayacucho, Cusco and L. walkeri vs. Ancash (precloacal pores or not), Ancash vs. L. tacnae (melanistic belly or not), Ayacucho vs. Cusco and L. walkeri (ringed pattern in ventral tail or not), Cusco vs. most L. walkeri (lateral markings or not). Based on the integration of molecular, different classes of morphological data, and niche identity test results, we conclude that Liolaemus populations from Ancash, Ayacucho, and Cusco can be delimited as separate species, and we describe these new species below.
Fifteen dorsal head scales (from a line drawn horizontally between anterior edges of external auditory meatus to anterior border of rostral). Dorsal head scales smooth except for the interparietal and surrounding scales, scale organs more abundant in prefrontal, internasal, and supralabial regions. Five scale organs on postrostral. Nasal scale in contact with rostral, separated from first supralabial by one scale, nasal bordered by eight scales; canthus separated from nasal by one scale. Six supralabials. Six lorilabial scales, three in contact with the subocular. Six infralabials. Auditory meatus oval (height 2.3 mm, width 1.2 mm), with three small, projecting scales on anterior margin. Seven convex, smooth temporals. Orbit-auditory meatus distance 4.9 mm. Orbit-anterior margin of rostral distance 6.3 mm. Rostral almost three times wider than high (width 2.9 mm; height 1.2 mm). Mental subpentagonal, about two times as wide as high (width 3.2 mm; height 1.7 mm). Interparietal pentagonal with an elongated posterior apex, bordered by eight scales, the parietal slightly smaller. Frontal quadrangular. Supraorbital semicircles complete on both sides. Semicircles formed by 6 scales. Four enlarged supraoculars. Six distinctly imbricate superciliaries on both sides. Eleven upper and ten lower ciliaries. Subocular elongate, 3.8 mm, longer than eye diameter (2.9 mm), separated from supralabials by a single, but interrupted row of lorilabials. Second supralabial elongate, 1.9 mm. Six lorilabials with single and double rows of scale organs. Sixth, fifth and fourth lorilabials contacting subocular. Preocular small, separated from lorilabial row by one scale. Postocular as large as preocular. Mental in contact with four scales: first infralabials (on each side) and two enlarged chin shields. Chin shields forming a longitudinal row of three enlarged scales separated one from the other by seven smaller scales. Scales of throat round, flat, and imbricate. Twenty-four gulars between auditory meatus. Longitudinal neck fold without keeled scales, that are similar to dorsal in size scales. Antehumeral pocket and antehumeral neck fold well developed. Forty-two scales between auditory meatus and shoulder (counting along postauricular and longitudinal neck fold), thirtytwo scales between auditory meatus and antehumeral neck fold. Gular folds absent.
Color pattern in preservation. Dorsal background color from occiput to base of tail greenish brown. Black continuous vertebral stripe present. Dark paravertebral marks. Paravertebral and vertebral fields of same background color. Dorsolateral stripes distinctly cream-color. Small dark cream-colored markings scattered in lateral field. Cream ventrolateral stripe, beginning on the upper auricular meatus, continuing across the longitudinal neck fold, through the shoulders, ending in the groin. Dark and small cream-colored marks in the ventral field. Black ventral color from about second third of head to femur, tibia and first third of tail. Dark and cream-colored small markings in first third of ventral head and two posterior thirds of tail.
Color pattern in life. Head dorsally brown with black and light brown dots. Subocular cream colored, dorsum bisected by a dark vertebral line. Vertebral field not conspicuous, bordering the vertebral line with a tenuous yellowish line. Paravertebral field with dark marks, bordered dorsally by a yellowish cream dorsolateral stripe. Lateral field with black and yellow reticulated pattern and white dots. Inconspicuous ventrolateral stripe, beginning on upper margin of auricular meatus, continuing from the longitudinal neck fold, through the shoulders, ending in the groin. Ventrolateral similar to lateral field but with more white dots. Fore and hind limbs same color as the paravertebral field, with diffuse dorsal markings. Dark, melanistic ventral color from about second third of head to femur, tibia and first third of tail. Dark and white dots in first third of ventral head and two posterior thirds of tail.
Variation. Variation in characters is summarized in Tables 1-4. There is sexual dichromatism. Adult males exhibit melanistic belly, cloacal region and throat, or mela-nistic belly only; adult females exhibit black and white spots on belly, cloacal region and throat, or yellowish belly and tail.
Etymology. The specific epithet chavin refers to the pre-Inca culture Chavin, which had its center close to the type locality and frequently depicted reptile figures on some of its most remarkable sculptures. The species name is in the nominative singular.
Distribution and natural history. Liolaemus chavin sp. n. is known from four localities in the central Andes, at elevations of 3535-4450 m in Ancash and Huánuco Departments in western central Peru (Fig. 11). It is the northernmost species of the subgenus Liolaemus.
Liolaemus chavin sp. n. was found active and under rocks in grassland and shrubland habitats at higher and lower elevations respectively (Fig. 8). In Pampas de Huamani the new species was usually found basking on grass up to 60 cm above the ground, and when they were disturbed they escaped into the base of grass clumps. Individuals basking on rocks were very rare in all localities. On cloudy days we found this species inactive hidden in the base of grass clumps, although some individuals were also found inactive under rocks. This species is viviparous; one female showed two uterine chambers per side with developed embryos, yolk and no visible shell in either chamber, and three females showed two uterine chambers per side with yolk, without developed embryos and no visible shell in each chamber. At the type locality no sympatric species of reptiles were found, but four amphibians are known: Pleurodema marmoratum (Duméril & Bibron, 1840), Telmatobius mayoloi Salas & Sinsch, 1996, Gastrotheca peruana and Rhinella (Bufo) spinulosa (Wiegmann, 1834) (Lehr, 2002;personal observations). Sympatric species at Catac include the anurans G. peruana, R. (Bufo) spinulosa, Telmatobius rimac Schmidt, 1954, T. mayoloi, andthe lizard Stenocercus chrysopygus Boulenger, 1900;at Carpa, G. peruana (Boulenger, 1900), R. (Bufo) spinulosa and P. marmoratum; at Pampas de Huamani, G. peruana, P. marmoratum and R. (Bufo) spinulosa;and at La Unión, Gastrotheca griswoldi Shreve, 1941, G. peruana, R. (Bufo) spinulosa and S. chrysopygus (Lehr, 2002).
Paratypes. Three males (MUSM 29681, 29687, 29678) and four females (MUSM 29679, 29689, 29680, 29682)  Diagnosis. Small (51.9 mm maximum SVL) Liolaemus closely related to L. chavin sp. n., L. tacnae, L. walkeri, and L. wari sp. n. (described below) (Fig. 1). It differs from L. chavin sp. n. and L. tacnae in having precloacal pores (males). Liolaemus pachacutec differs from L. wari sp. n. in having a partial or complete melanistic belly in adult males and in lacking a ringed pattern in ventral tail. Liolaemus pachacutec differs from most individuals (90%) of L. walkeri in lacking spots in the lateral field. In comparison with other species assigned to the L. alticolor group, L. pachacutec differs from L. chaltin in having precloacal pores in males. It differs from L. paulinae in the presence of a vertebral line and smooth neck scales. It differs from L. puna, L. alticolor and L. incaicus in having a partial or complete melanistic belly in adult males. It differs from L. aparicioi in lacking keeled temporal scales. It differs from L. bitaeniatus and L. pagaburoi in having a smooth dorsal surface of the head. It differs from L. pyriphlogos in the absence of red marks in lateral fields. It differs from L. variegatus in lacking keeled temporal scales, rugose dorsal head scales, and precloacal pores in females.
Dorsal head scales 16, dorsal head scales smooth, scale organs more abundant in loreal and supralabial regions. Two scale organs on postrostral. Nasal scale in contact with rostral, separated from first supralabial by one scale, nasal bordered by six scales; canthus separated from nasal by one scale. Four supralabials. Four lorilabials scales and one in contact with the subocular. Five infralabials. Auditory meatus oval (height 2.0 mm, width 1.0 mm), with two small, projecting scales on anterior margin. Six convex, smooth temporals (counting vertically from buccal commissure to posterior corner of orbit). Orbit-auditory meatus distance 3.9 mm. Orbit-anterior margin of rostral distance 4.3 mm. Rostral about two times wider than high (width 2.3 mm; height 1.0 mm). Mental subpentagonal, about two times as wide as high (width 2.5 mm; height 1.0 mm). Interparietal pentagonal with an elongated posterior apex, bordered by five scales, the parietal of similar size. Frontal trapezoidal.
Supraorbital semicircles complete on both sides. Semicircles formed by six scales. Five enlarged supraoculars. Six distinctly imbricate superciliaries on both sides. Eleven upper and lower ciliaries. Subocular elongate, 2.8 mm, longer than eye diameter (2.1 mm; measured between anterior and posterior commissure of ciliaries), separated from supralabials by a single, but interrupted row of lorilabials. Fourth supralabial elongate, 2.0 mm. Four lorilabials with single row of scale organs. Fourth lorilabial contacting subocular. Preocular small, separated from lorilabial row by one scale. Postocular as large as preocular. Mental in contact with four scales: first infralabials (on each side) and two enlarged chin shields. Chin shields forming a longitudinal row of four enlarged scales separated one from the other by six smaller scales. Scales of throat round, flat, and imbricate. Twenty-two gulars between auditory meatus. Longitudinal neck fold without keeled scales and smaller in size than dorsal scales. Antehumeral pocket and antehumeral neck fold well developed. Thirty-six scales between auditory meatus and shoulder (counting along postauricular and longitudinal neck fold), twenty-six scales between auditory meatus and antehumeral neck fold. Gular folds absent.
Color in preservation. Dorsal background color from occiput to base of tail brownish-green. Black thin continuous vertebral line present. No dark paravertebral marks. Paravertebral and vertebral fields with same background color. Distinct cream dorsalateral stripes. No marks in lateral field. Cream ventrolateral stripes, beginning on the posterior corner of the eye, continuing across the upper auricular meatus, the longitudinal neck fold, through the shoulders, ending in the groin. No marks in the ventral field. Melanistic venter on throat, femur, tibia, and belly. Small and scattered dark marks in chin area and ventrolateraly. Ventral tail melanistic near the cloaca, with a thin longitudinal stripe, first half with small marks lateral to the stripe.
Color pattern in life. Head dorsally brown with scattered black dots. Subocular white. Thin and faint black vertebral line. Paravertebral field without dark marks. Creamy dorsolateral stripes. Lateral field without marks. Faint cream-white ventrolateral stripe, beginning on upper margin of eye, continuing from auricular meatus, the longitudinal neck fold, through the shoulders, ending in the groin. Ventral field yellow. Forelimbs and chin scales white with scattered black dots. Melanistic belly, hind limbs, posterior two thirds of throat. Belly with scattered yellow dots laterally. Tail with a black region close to the cloaca, black longitudinal stripe and dots at each side of the stripe.
Variation. Variation in characters is summarized in Table 1-4. There is sexual dichromatism. Males have a complete or partial melanistic belly and throat, while females have a white or yellow belly and black spots on throat. Some males have orange and yellow dots on lateral belly and yellow dots on chin scales, and ventral field with orange and black dots.
Etymology. The specific epithet pachacutec refers to one of most important Inca rulers, Pachacutec, who built the best known Inca ruins, including Machu Picchu and Pisac, this last site at a higher elevation just above the type locality. The species name is in the nominative singular.
Distribution and natural history. Liolaemus pachacutec sp. n. is known from four localities in the central Andes, at elevations of 4023-4972 m in the departments of Cusco and Apurímac in southeastern Peru (Fig. 11). The species was found under rocks in grassland habitats (Fig. 9). It was found in sympatry at similar elevations with Liolaemus ortizi Laurent, 1982 andTachymenis peruviana Wiegmann, 1835. This species is probably viviparous; two females showed one or two uterine chambers per side, with an embryo and abundant yolk in each chamber, but without a visible shell.
Liolaemus wari sp. n. http://zoobank.org/67A997B8-5854-4D0D-B1E0-77680FF47512 http://species-id.net/wiki/Liolaemus_wari Figure 10 1999 Paratypes. Three males (MUSM 30823, BYU 50184, 50185) and ten females (MUSM 30824,30825,30826,30827,30828,30831,BYU 50186,50187,50191,50243) from the same locality as the holotype. Two males (MUSM 30830,30834) and three females (MUSM 30829, BYU 50188, 50190) from high area above the Historic Sanctuary  Diagnosis. Small (61.4 mm maximum SVL), slender Liolaemus, closely related to L. chavin sp. n., L. pachacutec sp. n., L. tacnae and L. walkeri (Fig. 1). It differs from L. chavin sp. n., L. pachacutec sp. n. and L. walkeri in having a ringed pattern on the ventral tail of adult males. It differs from L. pachacutec sp. n. in having spots in the lateral fields. Liolaemus wari differs from L. tacnae and L. chavin in having precloacal pores in males. In comparison with other species assigned to the L. alticolor group, L. wari sp. n. differs from L. chaltin in having precloacal pores in males. It differs from L. paulinae in lacking keeled neck scales. It differs from L. puna, L. alticolor and L. incaicus in having black spots on belly of adult males. It differs from L. aparicioi in lacking keeled temporal scales. It differs from L. bitaeniatus and L. pagaburoi in having a smooth dorsal surface of the head (rough to slightly dorsal surface of the head). It differs from L. pyriphlogos in the absence of red marks in the lateral field (red marks in the lateral fields present). It differs from L. variegatus in the absence of keeled temporal scales, rugose dorsal head scales and precloacal pores in females.
Dorsal head scales 14, dorsal head scales smooth, scale organs more abundant in loreal and supralabial regions. Five scale organs on postrostral. Nasal scale in contact with rostral, separated from first supralabial by one scale, nasal bordered by seven scales; canthus separated from nasal by one scale. Four supralabials. Five lorilabials scales and two in contact with the subocular. Four infralabials. Auditory meatus oval (height 2.0 mm, width 1.9 mm), with two small, projecting scales on anterior margin. Seven convex, smooth temporals (counting vertically from buccal commissure to posterior corner of orbit). Orbit-auditory meatus distance 4.6 mm. Orbit-anterior margin of rostral distance 7.9 mm. Rostral almost three times wider than high (width 2.7 mm; height 1.0 mm). Mental subpentagonal, about two times as wide as high (width 2.6 mm; height 1.2 mm). Interparietal pentagonal with an elongated posterior apex, bordered by seven scales, the parietal slightly smaller. Frontal trapezoidal. Supraorbital semicircles complete on both sides. Semicircles formed by 6 scales. Four enlarged supraoculars. Five distinctly imbricate superciliaries on both sides. Eleven upper and lower ciliaries. Subocular elongate, 3.2 mm, longer than eye diameter (2.3 mm; measured between anterior and posterior commissure of ciliaries), separated from supralabials by a single, but interrupted row of lorilabials. Second supralabial elongate, 1.6 mm. Five lorilabials with single and double rows of scale organs. Fifth and fourth lorilabials contacting subocular. Preocular small, separated from lorilabial row by one scale. Postocular as large as preocular. Mental in contact with four scales: first infralabials (on each side) and two enlarged chin shields. Chin shields forming a longitudinal row of three enlarged scales separated one from the other by six smaller scales. Scales of throat round, flat, and imbricate. Twenty-one gulars between auditory meatus. Longitudinal neck fold without keeled scales and smaller in size than dorsal scales. Antehumeral pocket and antehumeral neck fold well developed. Twenty-nine scales between auditory meatus and shoulder (counting along postauricular and longitudinal neck fold), 21 scales between auditory meatus and antehumeral neck fold. Gular folds absent.
Color pattern in preservation. Dorsal background color from occiput to base of tail brownish-green. Black continuous vertebral line present. Dark paravertebral marks. Paravertebral and vertebral fields with same background color. Highly distinct creamy-yellow dorsalateral stripes. Large dark and small cream marks in lateral field. Cream ventrolateral stripe, beginning on the posterior corner of the eye, continuing across the upper auricular meatus, the longitudinal neck fold, through the shoulders, ending in the groin. Dark and cream small marks in the ventral field. Black spots on throat, femur, tibia, posterior third of belly and laterally in anterior two thirds of belly. Small and scattered dark marks in chest and anterior two thirds of belly. Tail with dark horizontal rows.
Color pattern in life. Head dorsally brown with black dots. Subocular cream. A black vertebral band with a thin yellow stripe on the middle. The vertebral band has a thin white stripe on each side. Paravertebral field with dark marks with posterior white dots. Creamy-yellow dorsolateral stripes. Lateral field with black marks separated by cream diagonal stripes. Yellowhish-white ventrolateral stripe, beginning on upper margin of eye, continuing from auricular meatus, the longitudinal neck fold, through the shoulders, ending in the groin. Ventrolateral similar to lateral field and same color as the paravertebral field, with diffuse dorsal markings. Forelimbs, chest and belly yellowish-white with scattered and diffuse black dots. Black marks on hind  L. chavin, L. pachacutec, L. tacnae, L. walkeri, and L. wari. limbs, throat, and posterior third of belly. Tail with black horizontal bands separated by white bands.
Variation. The variation in morphological characters is shown in Tables 1-4. There is sexual dichromatism. Males have white or yellow belly and throat covered completely with black spots, yellowish belly and throat with black spots on posterior third of belly, or a melanistic belly on posterior third and cloacal region, with black dots on a white throat; females have white belly and yellowish throat with faint black dots, yellowish belly and throat with faint black spots, or yellowish belly and throat without spots. Adult males have white, yellowish and yellow tails with a conspicuous ringed pattern; adult females have white, yellowish or reddish ventral tails with or without a faint ringed pattern.
Etymology. The specific epithet wari refers to the pre-Inca culture Wari (600-850 AD), which had its center close to the type locality. The species name is in the nominative singular.
Distribution and natural history. Liolaemus wari sp. n. is known from seven localities in the central Andes, at elevations of 3768-4246 m in Ayacucho Department in eastern southern Peru (Fig. 11).
Liolaemus wari sp. n. was active on the ground or found under rocks in grassland (Fig. 10) and shrubland habitats. It was found in sympatry with another Liolaemus species belonging to the L. montanus series and the snake Tachymenis peruviana. This species is probably viviparous; three females each showed three uterine chambers per side; each chamber showed yolk, but with no developed embryos or visible shell.

Phylogenetic relationships
Surprisingly, our phylogenetic analysis showed that the three new species described herein plus L. tacnae and L. walkeri, assigned to alticolor-bibronii group, are strongly separated from the other members of this species group included in this study. Specifically, the species L. alticolor and L. incaicus assigned to the alticolor-bibronii group (Lobo et al. 2010, Quinteros 2013 were not recovered with L. tacnae, L. walkeri, and the three new species. Previous molecular based phylogenies did not include L. alticolor, L. tacnae and/ or L. walkeri , Morando et al. 2007, Schulte and Moreno-Roark 2010 and much of what these different topologies show (including ours) is probably an artifact of incomplete taxon/population sampling. Previous morphology-based phylogenies included better taxon sampling, but all of them recovered clades with low or no statistical support, and relationships of L. tacnae, L. walkeri, L. alticolor, and L. incaicus with each other and other species assigned to the alticolor-bibronii group are ambiguous.

Species delimitation and integrative taxonomy
We take our results based the mtDNA gene tree as a first step in species "discovery" (Carstens et al. 2013), and identify the Ancash, Ayacucho, and Cusco clades as "candidate species" (Morando et al. 2003, Avila et al. 2004. Comparative morphological and niche envelope assessments of these three clades revealed combinations of characters from three different lines of evidence, that unambiguously diagnose these groups as distinct from each other and from L. tacnae and L. walkeri (this is the second step of species delimitation -"validation" -following Carstens et al. 2013). This result highlights the need for using an integrative approach rather than a single line of evidence (e.g. morphology, usually meristic data only) to delimit species.
Our results show that normal tolerance intervals of continuous morphometric and meristic characters could not discriminate between any of these new species nor between L. tacnae and L. walkeri. On the other hand, discrete character analysis revealed some diagnostic characters, including: (1) the presence/absence of pre-cloacal pores in males distinguishing L. chavin and L. tacnae from L. pachacutec, L. walkeri, and L. wari; (2) the presence/absence of a complete or partial melanistic belly in adult males distinguishing L. chavin from L. tacnae; (3) the presence/absence of a ringed ventral tail pattern of adult males distinguishing L. wari from L. pachacutec and L. walkeri; and (4) the presence/absence of regular marks or spots in lateral fields distinguishing L. pachacutec from L. wari and from most (90%) individuals of L. walkeri. However, binomial tolerance intervals showed that all these "fixed" character states in our samples have a high probably of non-fixation when statistical inference is extended to consider large sample sizes. Despite these findings, we encourage the use of these binomial tests to place empirical evidence into a broader context, and to make investigators aware that tolerance intervals will become narrower as sample sizes increase, and that taxonomic decisions should be based on statistical populations not on samples (Zapata and Jiménez 2012). Moreover, samples taken at random are important for strong statistical inferences, but obtaining random samples in observational studies (such as in most taxonomic studies) is often impractical or impossible, and thus potential for bias is a serious concern (Ramsey and Schafer 2002). Besides this limitation, a statistical inference (such as those based on tolerance intervals) is better than no inference at all. However, statistic tests that evaluate differences in central tendencies (e.g., the ANOVA and Mann-Whitney U tests we used here) do not seem relevant as SDL criteria or for practical taxonomic purposes. For instance, most pairwise comparisons of SVL between focal populations (Ancash, Ayacucho, Cusco) and species (L. tacnae and L. walkeri) are significant in an ANOVA test at a confidence level of 0.05, giving the false impression that this character is useful for species delimitation or taxonomic identification, but tolerance intervals indicate that these populations and species completely overlap with respect to this character (Table 2).
Molecular analysis and, in most cases, niche identity tests, support our species units based on these few morphological characters, and in combination provide more robust hypotheses. Our model-based molecular phylogenetic analysis provided the basis for our "candidate species" hypotheses, but molecular phylogenetic analysis relies on the assumption that a chosen evolutionary model is a correct one (Posada 2009), and we recognize that in the absence of corroboration from independent data sets, mtDNA may often over-split species (Miralles and Vences 2013). However, assumptions are also pervasive in morphological and ENM analyses. Discovery of gaps in morphology assumes that discontinuities are not due polymorphisms, ontogenetic variation or phenotypic plasticity (Wiens andServedio 2000, Zapata andJiménez 2012), and ENM (especially those models based on background data and not true absence records) assumes that occupied distribution of a species is not reduced by biotic interactions and dispersal limitations (Peterson et al. 2011). Despite these assumptions, we think that robust hypotheses of species delimitation based on different data sets give stability to scientific names, provide the strongest inference about species boundaries, overcome overlapping character variation in any particular character system, and should be a prioritized research theme in systematics (Balakrishnan 2005, Padial and De la Riva 2006, Padial et al. 2010). In addition, we expect more exciting results when new molecular coalescent-based multi-locus and morphological multivariate methods can be applied to our data (Zapata andJiménez 2012, Camargo andSites 2013).

Northern limits of squamate viviparity in the high Andes
Liolaemus chavin is the northernmost viviparous species of the subgenus Liolaemus. Two recognized Liolaemus species present in the extreme northern range of the genus are L. robustus Laurent, 1992 and L. disjunctus Laurent, 1990 (subgenus Eulaemus). In the case of L. disjunctus, our recent fieldwork in the area of the species' type locality did not locate any specimen. The same result was found when we revisited localities near the type locality of L. disjunctus in 2012, and to our knowledge this species has not been collected at least since its original description and data on its reproductive mode are still lacking (Laurent 1990). On the other hand, the colubrid snake Tachymenis peruviana is another viviparous squamate widely distributed in the high Andes of Argentina, Bolivia, Chile and Peru. Its northern limits are in the department of La Libertad, Peru at about latitude 7°S, and no other viviparous squamate species are present in the high Andes of northernmost Peru, Ecuador, and Colombia.
What selective pressures might have limited the distribution of viviparous squamates in the high Andes? Although there are no field or experimental studies that have addressed this question in particular, one distributional pattern seems to be evident in the northern distributional limit of Liolaemus. For instance, on the Pacific Andean slopes at about latitude 15°S and south in Peru, viviparous Liolaemus species are present in lower, middle and higher elevations (C. Aguilar, personal observations), and oviparous lizards (genera Phyllodactylus, Ctenoblepharys and Microlophus but not Stenocercus) are only present at lower and middle elevations. However, on the Pacific Andean slopes at about latitude 12°S, Liolaemus species are only present at higher elevations and oviparous Stenocercus (Tropiduridae) species become common at lower and middle elevations, together with the above-mentioned oviparous genera. If we consider the actual northern limits of Liolaemus as represented by L. chavin, viviparous lizards in the high Andes do not extend north beyond about latitude 8-9°S. North of latitude 8°S, oviparous Stenocercus, Petracola and Riama (Gymnophthalmidae) species are the only lizard genera present in the high Andes of Peru and Ecuador. One interesting distributional and reproductive pattern that matches this change in reproductive mode in lizards is the distribution pattern of amphibians with direct development (genus Pristimantis). No Pristimantis species have been found in sympatry with northernmost Liolaemus species. At high elevations on the Pacific slopes, the northernmost Liolaemus species (L. chavin and L. robustus) have always been found with anurans having complete (genera Rhinella, Pleurodema and Telmatobius) or partial (Gastrotheca) indirect development.
Direct-development Pristimantis rely on high humidity substrates for egg development (Duellman and Lehr 2009), and what may have limited the distribution of direct-development frogs in the Pacific basin of southern Peru and northern Chile, and the Andean Plateau, is the formation of an Arid Diagonal area due to the interaction of the Humboldt Current and uplift of the Andes. If so, then a working hypothesis for the evolution of viviparity and placentation in some clades of Liolaemus is their relationship to the presence of these arid and hypoxic conditions. Arid environments in hypoxic middle and high elevations might be lethal to the development of oviparous lizard eggs. However, origins of viviparity in Liolaemus seem to be associated with shifts to cold climates (e.g., in the Oligocene; Schulte and Moreno-Roark 2010), thus supporting the cold climate hypothesis (CCH; Tinkle and Gibbons 1977). According to this hypothesis, viviparity has evolved to avoid lethal ambient temperatures in high elevations and latitudes, and through retention of eggs in the uterus coupled with female behavioral thermoregulation, this mode accelerates embryonic development (for a recent review see Sites et al. 2011). The CCH is a special case of a more general maternal manipulation hypothesis (MMH) where females can enhance fitness-related phenotypic attributes in offspring by manipulating thermal conditions during embryogenesis (Shine 1995). However, arid environments may be more important with increasing hypoxic conditions in high altitudes for the evolution of viviparity than cold climates, as has been suggested for Phrynosoma lizards (Hodges 2004, but see Lambert and Wiens 2013). In other words, altitude may be a surrogate of other selective factors important for the evolution of viviparity, not only cold climates (Hodges 2004). High altitude environments tend to be drier and have low oxygen conditions, and viviparous species may be able to provide a better oxygen environment for developing embryos via placental structures (Hodges 2004). Whether shifts in cold climates and/or appearance of arid zones along with Andean uplift are correlated with the origin of viviparity in Liolaemus should be tested with coalescent based multi-locus phylogenetic studies and a time-calibrated hypothesis of species relationships. Spots present in the lateral fields ( Collecting and exportation permits were issued by the DGFFS in Lima, Peru. We thank Ignacio De la Riva and an anonymous reviewer for providing valuable comments to our manuscript.

Key to Peruvian species of the subgenus
appendix 1