Completing a taxonomic puzzle: integrative review of geckos of the Paroedura bastardi species complex (Squamata, Gekkonidae)

The Paroedura bastardi clade, a subgroup of the Madagascan gecko genus Paroedura, currently comprises four nominal species: P. bastardi, supposedly widely distributed in southern and western Madagascar, P. ibityensis, a montane endemic, and P. tanjaka and P. neglecta, both restricted to the central west region of the island. Previous work has shown that Paroedura bastardi is a species complex with several strongly divergent mitochondrial lineages. Based on one mitochondrial and two nuclear markers, plus detailed morphological data, we undertake an integrative revision of this species complex. Using a representative sampling for seven nuclear and five mitochondrial genes we furthermore propose a phylogenetic hypothesis of relationships among the species in this clade. Our analyses reveal at least three distinct and independent evolutionary lineages currently referred to P. bastardi. Conclusive evidence for the species status of these lineages comes from multiple cases of syntopic occurrence without genetic admixture or morphological intermediates, suggesting reproductive isolation. We discuss the relevance of this line of evidence and the conditions under which concordant differentiation in unlinked loci under sympatry provides a powerful approach to species delimitation, and taxonomically implement our findings by (1) designating a lectotype for Paroedura bastardi, now restricted to the extreme South-East of Madagascar, (2) resurrecting of the binomen Paroedura guibeae Dixon & Kroll, 1974, which is applied to the species predominantly distributed in the South-West, and (3) describing a third species, Paroedura rennerae sp. nov., which has the northernmost distribution within the species complex.


Introduction
The gekkonid genus Paroedura is endemic to Madagascar and the Comoro islands Vences 2007, Hawlitschek andGlaw 2013). Its monophyly has been recovered by several molecular phylogenetic studies (e.g. Jackman et al. 2008, Glaw et al. 2018. Currently, 22 species are recognized, with more than half of the species scientifically named in the course of the 21st century (Uetz et al. 2020). While some species of Paroedura, like for example P. masobe or P. lohatsara, can be immediately recognized by their strikingly distinct external morphology Raxworthy 1994, Glaw et al. 2001), other species are morphologically more similar to close relatives and their identification based on morphology alone remains challenging. The application of molecular genetics furthermore has revealed the existence of deep phylogenetic lineages within the genus, many of them apparently qualifying as undescribed candidate species (Jackman et al. 2008, Nagy et al. 2012, 2018, indicating that the species diversity within the genus is still insufficiently documented. Integrative taxonomic approaches already helped to revise species diversity within the P. oviceps clade, a group of largely limestone-adapted species occurring on karstic outcrops, resulting in the description of four new species (Glaw et al. , 2018. Another major clade, the P. stumpffi clade, has been revised by Hawlitschek and Glaw (2013) and has resulted in the discovery and description of a new species from the Comoran Island of Mayotte.
Another subgroup within Paroedura whose monophyly is supported by molecular data (Jackman et al. 2008, 2018, Köhler et al. 2019) as well as by chromosomal data (a 2n=34 karyotype, compared to 2n=38 or 36 in other species of the genus, Aprea et al. 2013, Koubová et al. 2014) contains four currently recognized nominal species: P. bastardi (Mocquard, 1900), supposed to be widely distributed in southern and western Madagascar; P. ibityensis Rösler & Krüger, 1998, a montane species endemic to the region of Itremo and Ibity; and P. tanjaka Nussbaum & Raxworthy, 2000 and P. neglecta Köhler, Vences, Scherz & Glaw, 2019, both of the latter known only from the karstic Tsingy de Bemaraha formation in western Madagascar. Three recently published studies provided molecular phylogenetic trees in which the taxon P. bastardi was paraphyletic, with P. ibityensis nested among at least three strongly divergent mitochondrial lineages of P. bastardi and thus rendering P. bastardi sensu lato paraphyletic , 2018, Köhler et al. 2019. The data accumulated so far strongly suggest the existence of several cryptic or morphologically similar species contained in the name P. bastardi, particularly as some of these divergent lineages appear to occur in sympatry at certain localities (Köhler et al. 2019). We here use the term Paroedura bastardi clade to refer to the entire monophyletic group of species (P. bastardi, P. ibityensis, P. neglecta, P. tanjaka) and associated putative species-level lineages, and Paroedura bastardi complex to refer to P. bastardi sensu lato, i.e., including the cryp-tic lineages currently associated with this species. The present work aims at elucidating species limits and thus revising the species diversity within the P. bastardi complex, by combining a multilocus molecular dataset (one mitochondrial and two nuclear markers) with a detailed morphological study in an integrative taxonomic framework. Furthermore, we provide a phylogenetic analysis involving a representative set of samples for an extended set of molecular markers, to infer the relationships within this clade of geckos.

Sampling
Geckos were collected at night by opportunistic searching in potential habitats. Voucher specimens were euthanized by injection with ketamine solution or MS-222, Figure 1. Localities of specimens of the Paroedura bastardi clade analyzed in the present study. Colored shapes represent distinct lineages (species-level taxonomic units, anticipating the results of the present study) and open white rectangles highlight co-occurrence of two lineages. Black crosses superimposed on color shapes indicate specimens for which molecular data were not available, and whose taxonomic assignment is based exclusively on morphology. The localities of Marofandilia and Miandrivazo are based on additional genetic data (two 16S sequences) from Aprea et al. (2013). fixed in 90% ethanol or 5% formalin, then transferred to 70% ethanol for long-term storage. Tissue samples were stored in 99% ethanol. Field numbers refer to the collections of Angelica Crottini (ACZC), Franco Andreone (FAZC), Frank Glaw and Miguel Vences (FGMV), Frank Glaw (FGZC)  Koenig, Bonn (ZFMK). A map showing the geographic provenance of samples examined for DNA sequences and/or morphology is represented in Figure 1. The list of specimens examined for morphology is presented in the Appendix 1.

Morphological analyses
Measurements were taken using a digital Vernier calliper to the nearest 0.1 mm (taken by TB, Fig. 2): SVL = snout-vent length; TaL = tail length; HL = maximum head length from the anterior margin of the ear opening to the tip of the snout; HW = maximum head width; HH = maximum head height behind the eyes; distE = minimum distance between the bony edges of the orbits in dorsal view; AGL = axilla-groin distance; ED = maxi- Figure 2. Illustration of the morphological characters used in this study. Quantitative characters: SVL = snout-vent length; TaL = tail length; HL = maximum head length from the anterior margin of the ear opening to the tip of the snout; HW = maximum head width; HH = maximum head height behind the eyes; distE = minimum distance between and the bony edges of the eyeballs in dorsal view; AGL = axilla-groin distance; ED = maximum eye diameter; EO = maximal ear opening; TIBL = distance between the ankle and the knee (tibia length, knee flexed at 90°, left side). Qualitative characters: IO = minimum number of interocular scales separating the eyes (above the center of the eyes); SO = number of granular scales across the upper eyelid (transversally) ; SnoutS = arrangement of the mediodorsal scale rows of the snout tip: mostly forming two transverse rows of granules in contact (c), separated by a third median rows (s) or intermediate pattern (i). Pictures of P. bastardi taken by AM. mum eye diameter; EO = maximal ear opening; FoL = foot length (from ankle to tip of the 3rd toe, left side), HAL = hand length, distance between the wrist and the tip of the longest finger (until the insertion of the claw, which is not included, left side); TIBL = distance between the ankle and the knee (tibia length, knee flexed at 90°, left side). Qualitative characters (collected by AM) are the following: IO = minimum number of interocular scales separating the eyes (above the center of the eyes); SO = number of granular scales across the upper eyelid (transversally); SnoutS = arrangement of the mediodorsal scale rows of the snout tip: mostly forming two transverse rows of granules in contact (c), separated by a third median rows (s) or intermediate pattern (i); DigC = coloration of the toes, uniform (u) or bicolor (b). Multivariate statistical analyses were performed in R 3.5.2 (R Core Team 2018) with the R studio graphic interface R studio (RStudio Team 2016). To accommodate mixed data (here, qualitative, meristic, and mensural data), we used the function PCAmix from the package PCAmixdata (Chavent et al. 2017), which combines classical principal component analysis (PCA) with multiple correspondence analysis, and the package ade4 (Dray and Dufour 2007) for visualization of principal components. Sample sizes were too small to make application of allometric growth curves for full allometry corrections, so to remove the main impact of body size on metric measurements, all measurements were size corrected by dividing them by SVL. This procedure also has the advantage of providing simple ratios that can be easily calculated from raw measurements, even in the field, and used for diagnosis if consistent interspecific differences are detected. The PCAmix function was subsequently used to normalize the dataset (centered and normed) prior to the component analyses. Specimens displaying missing data and juveniles were removed from the dataset. Analyses were performed on 22 adult or subadult specimens from both sexes with 12 quantitative and one qualitative characters (DigC was excluded from the analysis as this trait presented too many missing data).

Molecular analyses for species delimitation
To delimit species among the P. bastardi group, analyses based on three independent datasets and involving all the tissue samples available for specimens in the P. bastardi clade (n = 57), plus samples representing 16 different species of Paroedura as out-group (see details in Appendix 3A) were carried out: one phylogenetic tree was inferred based on the mitochondrial DNA dataset (mtDNA: CO1 fragment) and two haplotype networks were reconstructed based on phased nuclear loci (nDNA: CMOS and KIAA1239).
To reconstruct a phylogenetic tree from mtDNA (CO1), the best fitting substitution model was determined in MEGA 7 (Kumar et al. 2016) based on the Bayesian Information Criterion (GTR+I+G). Phylogenetic inference under the Maximum Likelihood optimality criterion was carried out, with SPR branch swapping, and the robustness of nodes was assessed with 500 bootstrap replicates, in MEGA 7 (Kumar et al. 2016).
To assess the amount of allele sharing between populations and evaluate the amount of gene flow in contact zones, we built haplotype networks using statistical parsimony, as implemented in the program TCS v.1.21 (Clement et al. 2000) implemented in PopART (http://popart. otago.ac.nz). We used PHASE (Stephens et al. 2001) as implemented in DnaSP v.5 (Librado and Rozas 2009) to infer haplotypes from both sets of nuclear DNA sequences (CMOS and KIAA1239), based on alignments that were trimmed to ensure that all sequences have the same length (detailed list of the haplotypes is available in the Appendix 4). Networks were imported in Adobe Illustrator to add colors (according to the corresponding mitochondrial assignment) and additional connections representing co-occurrences of different haplotypes at a given locality.

Interspecific phylogenetic relationships within the P. bastardi clade
In addition to the previous analyses carried out for species delimitation purposes, three complementary analyses aimed at resolving deep phylogenetic relationships among species were undertaken on a reduced sub-sample (only one specimen per delineated species of the P. bastardi clade) but with a significantly greater number of markers: (1) a phylogenetic analysis based on nine concatenated nuclear markers ( . For some taxa and gene fragments, alternative samples were used, and in some cases, complemented by previously published sequences downloaded from Gen-Bank, or in the case of the outgroup taxon P. picta, extracted from the full genome of this species (Hara et al. 2018). See Appendix 3B for a complete list of samples used.
Phylogenetic analyses were carried out by partitioned Bayesian Inference, using MrBayes 3.2 (Ronquist et al. 2012) with best-fit substitution models and partitions calculated using Partition Finder 2.1 (Lanfear et al. 2016) (cf. Appendix 5). For each dataset, we performed one run of 50 million generations (started on random trees) and four incrementally heated Markov chains (using default heating values) each, sampling the Markov chains at intervals of 5000 generations. The first 12.5 million generations (burn-in = 25%) were conservatively discarded and the remaining trees were retained post burn-in and summed to generate a 50% majority-rule consensus tree. For the same individual sampling we reconstructed haplotype networks for each nuclear marker, using the same methodology described in the previous section.

Molecular analyses
Mitochondrial DNA phylogenetic tree (CO1). The deepest nodes of the tree are unresolved, and the lineages con-sidered to be part of the Paroedura bastardi clade are recovered as a monophyletic group with very weak support (bootstrap value = 75%). Nevertheless, three of the four nominal species currently recognized -i.e. P. ibi tyensis, P. neglecta, and P. tanjaka -are recovered as distinct and strongly supported monophyletic groups (≥ 99%). Together they form an unsupported clade (36%). In contrast, Paroedura bastardi sensu lato is not recovered as a monophyletic unit. Instead, this nominal species is divided into four distinct deep mitochondrial lineages, hereafter referred to as P. bastardi A, P. bastardi B, P. bastardi C, and P. bastardi D (which is represented by a single individual). In the CO1 tree, these form a basal polytomy within which the ibityensis-neglecta-tanjaka clade is nested. Of these lineages, P. bastardi A corresponds to P. bastardi Ca02 and Ca03 of Cocca et al. (2018) and Köhler et al. (2019), while the other lineages were all referred under the name P. bastardi in Köhler et al. (2019). Although the relationships among these four lineages are unresolved given the lack of support, it is worth noting that P. bastardi C and D are topologically recovered (with no support) more closely related to the ibityensis-neglecta-tanjaka clade than to P. bastardi A and B (Fig. 3A). The respective monophyly of P. bastardi B and C is fully supported (100%) whereas support for the monophyly of P. bastardi A is very low (60%). Within P. bastardi A, the subclade excluding samples from Tranoroa is supported by 96%. The monophyly of P. bastardi D cannot be discussed further given that it is here represented by a single specimen only.
To ensure clarity and the consistency of the comparisons between the different data sets, we arbitrarily choose to use the clustering suggested by the CO1 tree as an interpretative framework, and in the following report individuals based on their mitochondrial assignment as P. bastardi A, B, C or D.

Nuclear DNA networks (KIAA1239 and CMOS).
In terms of overall grouping of individuals, the KIAA1239 haplotype network (Fig. 3C) shows important similarities with the mtDNA tree. The specimens of each of the four main mitochondrial lineages also show closely related haplotypes in the KIAA1239 network, respectively. In each of these groups, KIAA1239 haplotypes differ most often by only one to two mutational steps (up to a maximum of five), whereas the groups are differentiated from each other by a minimum of four mutational steps (between two haplotypes of P. bastardi A and P. bastardi D) but most often by at least 10 steps. In terms of overall structuring, the only major difference between CO1 and KIAA1239 concerns the specimen ZCMV 12790 from Anja (the only representative of P. bastardi D in the CO1 tree), which clusters with individuals of P. bastardi A (Isalo, Toliara, and Tranoroa) in the KIAA1239 network.
The CMOS haplotype network ( Fig. 3B) is more conserved than the KIAA1239 network (22 distinct haplotypes, versus 37 for KIAA1239). From a structural perspective, the only notable difference between the clusters suggested by the CMOS network and lineages recovered in the CO1 tree is the sharing of a haplotype between P. ibityensis and three specimens of P. bastardi A (from Tranoroa). Moreover, the number of mutational steps between the specimens belonging to different mitochondrial lineages is lower than in KIAA1239.

Identity of P. bastardi D and of further mitochondrial variants.
Paroedura bastardi D is represented by a single sample in the molecular data set only. This sample exhibits discordant phylogenetic signal among markers, Figure 3. Summary of molecular results from one mitochondrial and two nuclear gene fragments in individuals of the Paroedura bastardi clade. A Maximum Likelihood phylogenetic tree based on the CO1 dataset, B, C haplotype networks inferred from the phased sequences of CMOS and KIAA1239, respectively. Circles represent haplotypes inferred by phasing (size proportional to their frequency in the individuals sequenced) and crossbars indicate the number of mutational steps. The colors assigned to the different clades in the CO1 tree are reported on the haplotypes of the corresponding specimens to facilitate comparisons. Gray lines and boxes highlight haplotypes co-occurring in Anja, Isalo, Tranoroa, and Bemaraha. and no voucher specimen was available for morphological comparison. Although we are confident that this clade represents a species distinct from P. bastardi C (no shared haplotype between these two clades co-occurring in Anja, suggesting an absence of gene flow), it is impossible, with the limited data, to determine whether this lineage is conspecific with P. bastardi A (affinities suggested by KIAA1239), conspecific with P. bastardi B (affinities suggested by CMOS), or if it represents a fourth distinct evolutionary lineage (as suggested by the CO1 tree). Further investigations, involving a larger number of samples and an examination of their morphological characteristics, are needed to clarify its status.
In previous analyses of the P. bastardi clade (e.g., Köhler et al. 2019), a further mitochondrial variant was included in phylogenetic trees based on CO1 sequences, i.e., the samples MRSN R3736, MRSN R3745, MRSN R2553, and FAZC 14631, all from Isalo (published by Cocca et al. 2018; GenBank accession numbers MH063359-MH063362) considered as candidate species Paroedura bastardi Ca01 by Cocca et al. (2018) and Köhler et al. (2019). In our study, four additional samples yielded CO1 sequences that clustered with this variant: MirZC 40, MirZC 91, MirZC 101, and MirZC 104, all from Kirindy. We have excluded all these CO1 sequences from our analysis and do not further discuss them because we are convinced they do not represent a biological entity but rather a segment of nuclear DNA of mitochondrial origin (NUMT), i.e., a nuclear pseudogene rather than a genuine CO1 sequence. This is based on four lines of evidence: (1) (4) the re-sequencing of these three specimens for the CO1 fragment (this study) confirm that the COI GenBank sequences MH063359-MH063362 are different and correspond to a NUMT, despite lacking stop codons.

Phylogenetic relationships of species and main line ages.
Our multi-gene phylogenetic analysis included single representative samples of all nominal species of the P. bastardi clade, plus P. bastardi A, B, and C. We did not include P. bastardi D because we were only able to sequence two nuclear loci (out of nine) for the single sample at our disposal.
The separate analyses of concatenated nuclear versus concatenated mitochondrial datasets (nine and five markers respectively) recovered incongruent topologies (Icong = 1.14 with P = 0.31, meaning that both trees are less congruent than or as congruent as expected by chance: de Vienne et al. 2007;Fig. 4). Both trees agree in placing the sympatric species P. neglecta and P. tanjaka as sister species. The mitochondrial tree places P. ibi ty ensis sister to the (P. neglec-ta, P. tanjaka) clade, and P. bastardi A sister to these three species, followed by P. bastardi C, and P. bastardi B. In contrast, the nuclear tree places a clade with P. bastardi C and P. ibityensis sister to the (P. neglecta, P. tanjaka) clade, and a clade with P. bastardi A and B sister to the remainder of species. Whereas the nuclear tree is very well supported (PP = 1.0 at each nodes), the mitochondrial tree is less robustly supported, with two nodes having PP < 0.9 (Fig. 4).
The tree involving the complete concatenated dataset (nuclear and mitochondrial) is relatively congruent with the mtDNA tree (Icong = 1.42 with P = 0.01, meaning that both trees are more congruent than expected by chance), but incongruent with the nDNA tree (Icong = 1.14 with P = 0.31, not more congruent than expected by chance). This suggests that the mtDNA phylogenetic signal is predominantly contributing to the topology of the combined mtDNA+nDNA tree (Fig. 4). As a main difference, in comparison with the mtDNA tree, in the combined tree the position of P. bastardi A and C is reversed, whereas the relationships among P. ibityensis, P. neglecta, and P. tanjaka, as well as the position of P. bastardi B as sister to all other lineages, are identical (Fig. 4).
Although this phylogeny and the underlying mito-nuclear discordance remains to be confirmed by more comprehensive phylogenomic studies, it is worth noting that the suggested relationships of P. ibityensis with the (P. neglecta, P. tanjaka) clade coincide with the occurrence of these two groups at localities relatively distant from one another in central Madagascar: P. ibityensis at high elevations on rocky mountain tops in the central high plateau, and P. neglecta and P. tanjaka in the central-west (Fig. 1). The intervening area between the high plateau and the western massifs such as the Tsingy de Bemaraha is poorly studied, with few reptile records (Glaw and Vences 2007), but it is possible that these apparently related species occur in geographical proximity to one another. Also, the relationships of P. bastardi C with the (P. ibityensis (P. neglecta, P. tanjaka)) clade is biogeographically plausible as this is apparently the lineage of P. bastardi sensu lato that occurs furthest north (Fig. 1). The discordance might be compatible with an ancient hybridization event between the P. bastardi and P. ibityensis lineages, but such a scenario will require more extensive future investigation to ascertain.

Morphological comparisons
Within the P. bastardi clade, the two species P. neglecta and P. tanjaka can easily be diagnosed because both of them have the nostril in contact with the rostral scale (Köhler et al. 2019). We therefore focused our morphological comparisons on the other species and lineages, mainly relying on PCA. Raw morphological data are presented in Appendix 6 and the contributions of each variable to the first Principal Components (PCs) in Appendix 7. The three groups within the Paroedura bastardi complex revealed by molecular analyses (P. bastardi A, B, and C) were relatively well separated from each other in the first two Principal Components (PCs) of our morphol-ogy-based PCA, but poorly separated from P. ibityensis (Fig. 5). Relative head size (HL/SVL, HH/SVL, and HW/ SVL; PC1) and body size (SVL; PC2) represent the traits that contribute the most (> 0.50) to the variation expressed along the first two PCs. Nevertheless, although the different values related to head shape (HL/SVL, HH/SVL and HW/SVL) contribute, along with body size, to unambiguously differentiate few pairs of species in the PCA (i.e., no overlap in the two-dimensional space of PC1 and PC2 comparing P. bastardi B and C with P. bastardi A), neither of these single variables is unambiguously diagnostic for any species. For instance, P. bastardi B and C reach larger body sizes than P. bastardi A and P. ibityensis (Fig.  5B), but with considerable overlap. Relative head length (HL/SVL) in P. bastardi B is > 0.3 for all nine adult and subadult individuals measured (0.302-0.336) and < 0.3 in P. bastardi A for five out of eight individuals measured (0.289-0.299); however, three P. bastardi A have values in the range of P. bastardi B (0.307-0.322) (see Appendix 7 for measurements and ratios of all individuals). Therefore, the value of these ratios for reliable, practical field identification of specimens is limited, although some statistical differences appear to exist and we expect relative head size having higher discriminatory power when comparing lineages of the P. bastardi complex with several other Paroedura species.

Integrative species delimitation: Molecular evidence of reproductive isolation in contact zones
While P. ibityensis is a highland species apparently not occurring in sympatry with any other lineage of the P. bastardi clade, among the other lineages there are several instances of co-occurrence without genetic admixture that are informative to infer reproductive isolation, as examined in the following.
(1) Sympatry in Isalo. Several specimens from Isalo are placed in P. bastardi A (ACZC 1828, 6464, 7932, and 7934, plus MH063363-69 available on GenBank) in the mitochondrial (CO1) tree, and those sequenced for the nuclear genes consistently show affinities to other P. bastardi A specimens (from Toliara and Tranoroa) in the KIAA1239 and CMOS haplotype networks. In contrast, three other specimens from Isalo are placed in P. bastardi C (ACZC 6438, 6534, and 7938) in the CO1 tree and consistently show affinities to other P. bastardi C specimens (from Kirindy and Anja) in the KIAA1239 and CMOS networks.
(2) Sympatry in Tranoroa. Three specimens of P. bastardi A in the CO1 tree (FGZC 352, 353, and 354) present affinities to other samples of this lineage in the KIAA1239 network, whereas in the CMOS network, they share a haplotype with specimens of P. ibityensis. In contrast, two other individuals of Tranoroa placed in the P. bastardi B clade in the CO1 tree (FGZC 332 and 327) have their nuclear haplotypes consistently recovered as closely related (or identical) to those of other P. bastardi B specimens (from Tolagnaro).
(3) Sympatry in Anja. Haplotypes of two specimens placed in P. bastardi C in the CO1 tree (ZCMV 12789 and 12791) consistently show affinities to other specimens of P. bastardi C (from Kirindy and Isalo) in both the KIAA1239 and CMOS networks. In contrast, the only representative of P. bastardi D (ZCMV 12790), shows either affinities to specimens of P. bastardi A (KIAA1239 network) or to P. bastardi B (CMOS network).
(4) Sympatry in Bemaraha. Paroedura neglecta and P. tanjaka are sympatric in this locality. Both species represent distinct mtDNA lineages and do not share any nuclear haplotypes (CMOS, KIAA1239).
Such patterns, consistently involving the same clusters of specimens for each of the three sequenced markers, strongly support the hypothesis of at least four occurrences of sympatry with reproductive isolation between different pairs of lineages. To summarize, our dataset supports reproductive isolation between: P. bastardi A and C in Isalo, P. bastardi A and B in Tranoroa, P. bastardi C and D in Anja, and between Paroedura neglecta and P. tanjaka in Bemaraha (Fig. 3).

Integrative species delimitation: Morphological divergence in contact zones
The molecular evidence for reproductive isolation between different pairs of sympatric lineages of the P. bastardi complex in Anja, Isalo, and Tranoroa also offers an opportunity to understand more precisely the morphological variability encountered, by differentiating the features whose variability is due to intraspecific polymorphism within the same species (and which are most often variable across the distribution range), from those which characterize each of the considered species. The distinction between these two cases of polymorphism (intra-versus interspecific) could even be facilitated if the interspecific morphological differentiation has been exacerbated in sympatry by reinforcement mechanisms (i.e. character displacement). Unfortunately, the heterogeneity of the material at our disposal, and the lack of multiple adult individuals, did not always allow us to objectively link molecular and phenotypic data. For instance, for Anja, we had only two voucher specimens available for morphological examination (ZSM 850/2010, member of the clade C, and ZSM 779/2009, also a member of the clade C, but a subadult; 16S sequence available in Gen-Bank confirm species identity: MT981881); and for Isalo, we had only one voucher specimen available (ZFMK 59808), too old to be genotyped using classical PCR approaches.
Luckily, our sampling from the population of Tranoroa was richer and provided us with several specimens unambiguously belonging to each of the two co-occurring lineages: three genotyped vouchers of P.

Adults and subadult specimens in Tranoroa.
Despite the limited sample size and lack of molecular assignment of one specimen (ZSM 58/2004), morphological quali-tative examinations and multivariate statistical analyses appear to confirm the morphological differentiation between the two lineages occurring in this locality. Spec- Figure 5. Morphological differentiation among four lineages of the Paroedura bastardi clade, here considered as representing distinct species. A Scatterplot of first two principal components (PC1, PC2) from a Principal Component Analysis using the morphological variables (adults and subadults only, males and females not analysed separately). Genotyped specimens are highlighted in bold italics (other specimens tentatively assigned to one of the molecular clusters based on morphological examination). B Violin plot of SVL (in mm), illustrating lower maximum body sizes in P. ibityensis and especially, in P. guibeae despite overlap of size ranges among all lineages. See appendices 6 and 7 for details. imens of P. bastardi A are indeed smaller in size than those of P. bastardi B, have fewer interocular scales (3 to 4 versus 5 to 6); and visually differ by less prominent tubercles, by a paler coloration, and by the presence of a banded pattern with alternating dark and light stripes on fingers and toes (versus a darker coloration, more spiny tubercles and a uniform light coloration of digits in P. bastardi B; Fig. 6A).
Juvenile specimens in Tranoroa. Juvenile coloration can be informative in lizard taxonomy, as their patterns tends to be more pronounced (more contrasted and differentiated) than in adults, and especially in Paroedura, juveniles often have a distinct and species-specific color pattern. For instance, Glaw et al. (2001) emphasized the distinct juvenile coloration of P. lohatsara, and Köhler et al. (2019) used the unique number of dorsal bands in juvenile P. neglecta in the diagnosis of this species from other Paroedura not belonging to the P. bastardi clade. It is nevertheless important to take into consideration the fact that in Paroedura, the coloration intensity changes dramatically (and thus likely rapidly) during ontogeny. Differences between juveniles in color and pattern should therefore be interpreted with great caution, as it can be difficult to distinguish ontogenetic from inter-individual or inter-specific variation. We had the opportunity to examine two juvenile specimens of approximately identical size, and therefore supposedly at relatively similar developmental stages, from Tranoroa ( Fig. 6B; ZSM 178/2004, genetically assigned to P. bastardi B, SVL = 30.0 mm, and ZSM 188/2004, genetically assigned to P. bastardi A, SVL = 26.4 mm). The latter specimen (cluster A) has a dull dorsal coloration whereas the former specimen (cluster B) shows a highly contrasted color pattern consisting of a dark brown dorsal background with two very light transverse bands. This pattern is also observed in five other juveniles of similar size: ZSM 42/2004 from Tolagnaro; one syntype of P. bastardi, MNHN 1899.0338, from Tolagnaro; and ZFMK 48434-48436, all from Berenty), all from a part of Madagascar where only the presence of P. bastardi B has been confirmed by molecular data. These six juvenile specimens (the one from Tranoroa and the five others from the region of Tolagnaro) also share a similar pattern in the parietal region, which roughly evokes a butterfly or a bat-like silhouette perforated in its center (Fig.6B, 6C). The shape of this pattern is remarkably constant across these six specimens and seems to characterize the juveniles of P. bastardi B.

Taxonomic conclusions
The three molecular datasets (Fig. 3) consistently suggest that Paroedura bastardi A, B, and C represent independent and reproductively isolated evolutionary lineages, i.e., distinct biological species, whereas the status of Paroedura bastardi D remains to be assessed in the future. This conclusion is supported by morphological data, despite limited genotyped material available for examination. From a nomenclatural point of view, two taxon names unambiguously referring to populations of the Paroedura bastardi complex are available: (1) Phyllodactylus bastardi Mocquard, 1900, whose type series is composed of five syntypes: three females (including two subadults) collected by Grandidier (MNHN 1899.0337, "environs de Tuléar" (= vicinity of Toliara) and MNHN 1899.0338 and 1899.0339, both from "Fort-Dauphin" (= Tolagnaro), and two specimens collected by M. Bastard (MNHN 1900.0006 and1900.0007, both from the "pays Mahafaly" (= Mahafaly country), which were larger in size and darker in coloration at the time the species was named (according to Mocquard 1900). In all likelihood, this type series is mixing different clades: MNHN 1899.0337 is likely a member of P. bastardi A, the only lineage of the P. bastardi complex known to be present in Toliara, and MNHN 1899.0338 and 1899.0339 are very likely members of P. bastardi B, the only lineage of the P. bastardi complex known from Tolagnaro. Although the PCA shows morphological resemblances between the specimen MNHN 1900.0006 and those assigned to P. bastardi B, we consider that reliable assignment of MNHN 1900.0006 and 1900.0007 to a given genetic lineage based on morphology alone is not possible because the "Pays Mahafaly" designates a vast region of the South of Madagascar where we suspect several lineages occur.
We therefore elect to designate the syntype MNHN 1899.0338 (see Fig. 6C and Fig. 7 for photos of the head and the entire body, respectively) from Tolagnaro as lectotype of the name Phyllodactylus bastardi Mocquard, 1900. Although this specimen is a juvenile, we consider it the best choice because (1) all the genotyped specimens from this locality are unambiguously assigned to P. bastardi B, (2) this specimen presents the typical "butterfly-shaped" pattern (albeit discolored and hardly visible) observed in another juvenile unambiguously genetically assigned to P. bastardi B (ZSM 178/2004, Tranoroa) and in all other juveniles (not genotyped) from the extreme South East of Madagascar (ZSM 42/2004 from Tolagnaro, ZFMK 48434-48436 from Berenty) where only P. bastardi B is known to occur (Fig. 6B, 6C). The choice of this lectotype is also motivated by our aim to promote nomenclatural stability, by limiting the establishment of new names (it enables us to erect one of the former synonyms, guibeae, for P. bastardi A, rather then erecting a new name; see below). As a consequence of our lectotype designation, the name Paroedura bastardi (Mocquard, 1900) becomes restricted to the specimens corresponding to P. bastardi B (= Paroedura bastardi sensu novo), and four of the original syntypes (i.e. MNHN 1899MNHN .0337, 1899MNHN .0339, 1900MNHN .0006, and 1900.0007) lose their status of name-bearers and become paralectotypes.
(2) Paroedura guibeae Dixon & Kroll, 1974, whose holotype is an adult male (FMNH 73049) collected at "10 km S Betroka 23° 18' S 46° 06' E, Madagascar". This nomen has been synonymized with P. bastardi by Nussbaum and Raxworthy (2000). We consider that FMNH 73049 very likely belongs to P. bastardi A because it shows banded patterns on digits, a trait that characterizes most of the specimens assigned to this lineage and that is always absent in the specimens assigned to P. bastardi B and C (cf. Fig. 6A). As a consequence, we here resurrect the name Paroedura guibeae Dixon & Kroll, 1974 and unambiguously apply it to P. bastardi A.
Note: The spelling of the epithet of Paroedura guibeae Dixon & Kroll, 1974 has been subsequently changed into guibei by Michels and Bauer (2004), who supported this decision by the fact that this species was originally dedicated to Mr. (not Ms.) Jean Guibé. Nevertheless, from a nomenclatural perspective, this emendation of epithet is incorrect (unjustified emendation according to the articles 33.2 and 33.4 (ICZN 1999); see also Dubois 2007 for this particular case) and for this reason we retain the original spelling.
No earlier name is available for the species corresponding to P. bastardi C, and consequently, this lineage is, in the following, described as a new species. Figure 6. Differences in coloration between juveniles of lineages provisionally named P. bastardi A and P. bastardi B (corresponding to the species P. guibeae and P. bastardi, according to the taxonomic hypothesis proposed herein), with a special emphasis on those from Tranoroa. A Digit coloration. Juveniles of P. bastardi A present a very characteristic banded pattern on fingers and toes (finger schematically represented in green). The same pattern is also present in the holotype of P. guibeae, and this name is therefore assigned to P. bastardi A. B Juvenile of P. bastardi A showing a dull dorsal coloration, whereas the juveniles of P. bastardi B show a highly contrasted pattern color consisting in a dark brown dorsal background with two very light transverse bands. C Detail of the dorsal side of the head of the newly designated lectotype of P. bastardi (above, a schematic drawing represent the "butterfly" pattern characterising juveniles of P. bastardi B, which is also present, although hardly distinguishable (probably faded) in the lectotype of P. bastardi). Scale bars = 5 mm.
Paroedura rennerae sp. nov. can be distinguished from most other currently recognized Paroedura species by the presence of only three broad light crossbands on the dorsum in juveniles and subadults (the first one between forelimbs, the second one at midbody, and the third one between hindlimbs) versus four light crossbands in all other species except those of the P. bastardi clade (P. bastardi, P. guibeae, P. ibityensis, P. neglecta, and P. tan jaka, which all have three crossbands) and P. ovi ceps and P. vahiny (in which the juvenile coloration is still unknown). It can be distinguished from P. gracilis by larger dorsal scales, absence of a white tip to the original tail, absence of a raised vertebral ridge on the dorsum and shorter forelimbs, which do not extend forward beyond tip of snout; from P. masobe by much smaller eyes and absence of a dorsal row of paired spines on the tail; from P. fasciata, P. ho ma lor hina, P. hordiesi, P. vahiny, and P. spelaea by presence of spines on the original tail (versus absence); from P. gracilis, P. homalorhina, P. kloki, P. maingoka, P. masobe, P. oviceps (from its type locality Nosy Be), P. picta, P. spelaea, most P. tanjaka, and P. vahiny by the presence of prominent dorsal tubercles arranged in regular longitudinal rows (versus rather irregular rows of dorsal tubercles).
Within the P. bastardi clade, the species can easily be distinguished from P. tanjaka and P. neglecta by the absence of contact between the nostril and the rostral scale (versus presence). It can be distinguished from P. ibityensis by larger maximum SVL (> 70 mm versus 61 mm). In comparison with P. bastardi sensu novo and P. guibeae, the new species can be distinguished by the presence of a very sharp and contrasting dark transverse pattern, evoking the shape of a thin curly-bracket ( { ), in the occipital region and delimiting the skull from the neck. Moreover, Paroedura rennerae sp. nov. is unambiguously larger in size than P. guibeae (adult SVL > 70 mm versus < 60 mm in P. guibeae), and its dorsal tubercles are more prominent. It also lacks striped fingers (versus striped in P. guibeae), and the light patch on its head lacks concave anterior edge and central vacuity in juveniles (versus both present in P. bastardi).
Description of the holotype. Adult female in very good condition, with the exception of the regenerated tail tip, which is amputated (ca. 10 mm missing). Head distinctly wider than neck, as wide as the body. Canthal ridges relatively well developed with a marked median depression. Ear opening is a vertical slit. Tail regenerated, nearly round (slightly flattened dorso-ventrally) in cross section in its proximal part; ventral pygal section of tail with a pair of poorly developed postcloacal sacs. Digits distinctly expanded at tips. Rostral scale rectangular, more than two times wider than tall and barely wider than mental. Nostrils separated from the rostral by prenasals. The two enlarged prenasals in contact with rostral and first supralabials, both separated by a single small granular scale. 12/12 (left/right) smooth supralabials, followed by two carinated tubercules above the mouth commissure. Eyes desiccated. Scales covering canthal ridges, loreal, temporal and periphery of the parietal region distinctly enlarged, spiny and tuberculate. Scales covering the dorsolateral side of neck and body heterogeneous, with enlarged, spiny, carinate and tuberculate scales regularly separated from each other by one (most often transversally) to three (most often longitudinally) rows of small, flat and juxtaposed scales or, along the vertebral line, by a single distinct row of smaller spiny tubercles. Seventeen longitudinal rows of tuberculate scales at midbody. Dorsal scales of forelimbs and hindlimbs mostly tuberculate and keeled, with a tetrahedral outline. Ventral scales of forelimbs distinctly smaller than surrounding ventral scales of the body. Three transverse rows at the base of the tail with six very spiny pygal scales per row. Ventrally, six rows of pygal scales squared and flat. Tail segments with irregular transverse row of spiny tubercles. Mental triangular, bordered posteriorly by a pair of elongated, irregular hexagonal postmentals. Each postmental in contact with six scales: other postmental, mental, first infralabial, one enlarged lateral gular, one smaller posterolateral gular, and one larger central gular. First three infralabials slightly larger (taller) than others. Gulars small, slightly granular. Ventrals of chest and abdomen flat and roundish. Proximal subdigitals in rows of mostly two. One pair of squarish terminal lamellae. Claws curving downwards between terminal pads of digits.
After nine years in alcohol (Fig. 9), head dorsally ochre colored with a pair of dark temporal bands, run-ning from the eyes to contact each other in the nuchal region. The contrast between the ochre dorsal side of the head and the darker temporal/nuchal bands is amplified by a dark blackish curly-bracket shaped transverse stripe in the occipital region (Fig. 8). Area along the upper lip alternating taupe-gray and cream. Body dorsally brown with three distinct lighter ochre (strongly contrasting thanks to very dark anterior and posterior borders) crossbands fading at the flanks: one transverse light crossband below forelimb insertion (width along the vertebral axis 6.1 mm), one distinctly broader light bow-tie-shaped crossband at midbody (width along the vertebral axis 9.6 mm), and one slightly less distinct band between the hindlimbs (6.3 mm). Dorsal surfac-es of forelimbs and hindlimbs slightly marbled with brown and ochre (hindlimbs not darker than forelimbs). Flank coloration lighter than dorsum, fading gradually towards the ventral surface. Ventral coloration (throat, chest, abdomen, ventral parts of forelimbs and hindlimbs) cream (very slightly pigmented on the throat and chest).
Coloration in life (Fig. 9). The coloration of the preserved specimen is very similar to that of the living in- dividual, although it is slightly duller (the contrasts are a little less strong and the colors a little less warm).
Variation. Both paratypes, from Anja, present a lighter and more contrasted color pattern, with sharper dark lines (anterior and posterior margin of the light dorsal cross bands and dark curly-brackets delimiting the occipital region) (Figs 7-8). The tail of ZSM 779/2009, which is not regenerated, has seven pairs of regularly alternating light brown and cream crossbands delimited by dark (brown) transverse stripes, whereas the tail of ZSM 850/2010 (regenerated) is cream with five thin transverse zig-zagging dark brown stripes (Fig. 7). The specimen ZFMK 59808, juvenile (but not neonate) with its original tail, is relatively similar in coloration to the paratype ZSM 779/2009, although slightly paler. In contrast with adults with a regenerated tail, younger specimens (ZSM 779/2009, subadult and ZFMK 59808, juvenile) present very regular rows of spiny tubercules all along the tail (around 20 rows). See also Table 1 for the variation in measurements.
Etymology. This new species, elegant and prickly, is dedicated to Susanne Renner, eminent botanist and evolutionary biologist, and Professor Emeritus of the University of Munich, in recognition of her substantial contributions to taxonomy and her invaluable collaboration in the framework of the "Taxon-omics" priority program of the German Research Foundation, DFG.
Habitat, habits, and distribution. Paroedura rennerae is reliably known from five localities, some of them relatively distant from each other, suggesting this species is widely distributed in the central/southern region of Madagascar. In the dry forest of Kirindy CNFEREF, specimens have been observed on vertical surfaces (tree trunks, wooden walls of the CNFEREF camp huts), around 1 to 2 m above the ground. Like other members of the P. bastardi species complex, it is quick to bite when handled. In Anja, several specimens have been collected on granitic boulders. ZSM 779/2009 was found in a large cavity below two large granitic boulders, in a quite humid environment. In this cavity, ZSM 779/2009 and other individuals were found on the walls. In Isalo, specimens belonging to this species were found at two sites (Zahavola and Namazaha Valley). These individuals were found within the canyons of the sandstone Massif in shaded areas and in close proximity to a small cave or a small waterfall, again in quite humid microhabitats. Two additional 16S sequences confirm the presence of this species also in Marofandilia and Miandrivazo (GU129005 and GU128989, Aprea et al. 2013). All specimens have been observed at night or near dusk.

Diversity, biogeography and species delimitation in the Paroedura bastardi clade
The genus Paroedura has seen a remarkable increase in the number of recognized species. Only nine species were recognized by Dixon and Kroll (1974), whereas Köhler et al. (2019) distinguished 22 species, and already pointed to the probable existence of additional species in the P. bastardi clade. By resurrecting P. guibeae and naming P. rennerae, the genus now contains 24 species, and we suspect that additional unnamed species still exist -for instance the enigmatic P. bastardi D identified herein.
As seems to be typical for many other reptiles in Madagascar, Paroedura contains several regional endemics with moderately large distributions, as well as a handful of extremely range-restricted species. For instance, several Paroedura species specialized on karstic limestone often inhabit caves, and are restricted to particular limestone massifs (Glaw et al. 2018). The Ankarana Massif in northern Madagascar harbors microendemic gecko species of the genera Paroedura (i.e., P. homalorhina; Glaw et al. 2018), Blaesodactlyus (B. microtuberculatus; Jono et al. 2015), Geckolepis (G. megalepis; Scherz et al. 2017), Phelsuma (P. roesleri; Glaw et al. 2010), and Uroplatus (U. fetsy; Ratsoavina et al. 2019), and at least in the latter case, the closest relative of the microendemic species is more widespread. The P. bastardi clade also contains microendemic species (P. neglecta and P. tanjaka, only known from the Tsingy de Bemaraha limestone massif) and species spread over wider ranges in southern Madagascar (P. bastardi, P. guibeae, and P. rennerae), but a further example for microendemism in the genus might be P. bastardi D from Anja, whose status we could not reliably determine due to the lack of material. Anja Reserve is characterized by a specific habitat of large granitic boulders where range-restricted species of ground chameleons (Brookesia brunoi: Crottini et al. 2012) and geckos (Phelsuma gouldi: Crottini et al. 2011; Paragehyra felicitae: Crottini et al. 2014) occur. We recommend extended sampling in this area to clarify the status of P. bastardi D, and to identify possible additional microendemic species occurring at this site. The improved knowledge on the taxonomy of the P. bastardi complex will, in the future, allow specifically targeting questions on possible ecological or behavioral specialization of the taxa involved. Especially in cases of sympatric occurrence, we assume that possibly, the taxa involved may prefer different substrates. We have found P. bastardi mostly on tree trunks and other vertical wooden surfaces, and the same is true for P. rennerae in Kirindy, but not in Anja, where at least ZSM 779/2009 was found on the walls of a cave-like large cavity. In Isalo, P. rennerae and P. guibeae occur syntopically at least at one site (Zahavola). However, while P. rennerae was found in quite humid microhabitats (and always inside the canyons), P. guibeae was mostly found on rock surfaces in open grasslands near canyon entrances.
Instances of sympatric occurrence of lineages may not only serve to understand their ecological specialization; they can also provide one of the most reliable lines of evidence to delineate species, and this has been applied both by Köhler et al. (2019) and in this study for the Paroedura bastardi complex. While numerous approaches to species delimitation use geographic sorting of mitochondrial haplotypes as argument for species distinctness (see Sites and Marshall 2004), we emphasize that sympatry can be one of the most powerful arguments for species distinctness, if used properly. The essential aspect of sympatry is that two lineages co-occurring at the same site at the same time while maintaining their genetic identity, in principle must be reproductively isolated. However, there are several caveats that need to be taken into account, and the Paroedura example serves to recapitulate these. (1) First, in highly mobile species or long-distance migrants, conspecific individuals representing different geographic variants may regularly occur in sympatry outside of the breeding season, thus not representing evidence for reproductive isolation. Such situations however are extremely unlikely in less mobile, small-sized taxa, and thus can be excluded for geckos. (2) Secondly, exceptional events such as human translocation can bring specimens of different geographic variants into situations of immediate co-occurrence where genetic admixture will only become apparent after multiple generations. While human translocation is common in commensal geckos, it is unlikely to be a major factor in Paroedura, especially given that many of our collections were made in natural areas and a distinct phylogeographic structure was obvious in several species such as P. guibeae. (3) Sympatric variation in one marker alone is insufficient for sound conclusion. Different alleles of both nuclear and mitochondrial genes can co-exist in the same population, and especially in mitochondrial DNA, examples of fast divergence in geographic isolation, subsequent admixture and thus co-occurrence of substantially diverged haplotypes at the same site are common. In such cases, it is paramount to assess the variation in other, independent sets of characters (Padial et al. 2010), which can be morphological, ecological, or behavioral, or other unlinked molecular markers such as nuclear gene sequences.
In the P. bastardi clade, there are at least two examples that will require future scrutiny: in P. tanjaka, three mitochondrial haplogroups of substantial divergence co-occur in the Tsingy de Bemaraha (Fig. 3; see also Köhler et al. 2019), and in P. guibeae, two haplogroups co-occur in the Isalo Massif ( Fig. 3; treated as different candidate species Ca02 and Ca03 by Cocca et al. 2018). In these cases, we found no evidence for divergence in the nuclear genes studied, and had too limited material available for a thorough morphological comparison, and therefore treated the respective individuals as conspecifics (representing deep conspecific mitochondrial lineages sensu Vieites et al. 2009). (4) Lastly, it has to be considered that particularly closely related lineages in a state of incipient speciation may be connected by hybrid zones, and the width of these is informative about the species status of the taxa involved (e.g., Dufresnes et al. 2020). Across such hybrid zones, limited sampling from one site (i.e., few samples sequenced for few molecular markers) involves the risk -even if unlikely -to choose individuals where these markers show a concordant signal and thus suggest reproductive isolation. In the case of the P. bastardi clade, however, this situation is very unlikely, especially because in Tranoroa we could assess a strong concordance not only of unlinked molecular markers, but also of molecular with morphological differentiation.
To conclude, we advocate that sympatry of lineages without genetic admixture is one of the most immediate means to delimit species, even with limited sample size, if several other biological phenomena are appropriately considered and can be excluded. It is important to emphasize the need for concordance of various characters or unlinked markers; by no means should new species be based on co-occurrence of different, even strongly diver-gent mitochondrial haplotypes alone. The probability of recovering by chance concordant differentiation among different unlinked markers, or between molecular markers and morphology, decreases drastically with increasing numbers of markers and sampled individuals, and we suggest that this could be taken into account by probabilistic approaches to species delimitation.

Funding
The Portuguese National Funds through FCT -Fundação para a Ciência e a Tecnologia -supported the Investigador FCT (IF) grant to AC (IF/00209/2014). This study would not have been possible without the support to AM and TB in the framework of the Taxon-Omics priority program of the Deutsche Forschungsgemeinschaft (SPP 1991 -RE 603/29-1).