Geckos were collected at night by opportunistic searching in potential habitats. Voucher specimens were euthanized by injection with ketamine solution or MS-222, 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) and Aurélien Miralles (MirZC). Vouchers were deposited in the Museo Regionale di Scienze Naturali (MRSN), the Parc Botanique et Zoologique de Tsimbazaza (PBZT), the collections of the Mention Zoologie et Biodiversité Animale of the Université d’Antananarivo (UADBA-R) and the Zoologische Staatssammlung München (ZSM). Additional specimens were examined from the collections of the Field Museum of Natural History, Chicago (FMNH), the Muséum National d’Histoire Naturelle, Paris (MNHN), the Zoological Museum Amsterdam (ZMA), and the Zoologisches Forschungsmuseum Alexander 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.

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 = maximum 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).

Genomic DNA was extracted from tissue samples using proteinase K (10 mg/ml) digestion followed by a standard salt extraction protocol (Bruford et al. 1992). To delineate species, for all available tissue samples, we amplified and sequenced a DNA fragment of a mitochondrial (mtDNA) gene (cytochrome oxidase subunit 1, CO1) and fragments of two protein-coding nuclear genes (nDNA: oocyte maturation factor, CMOS; and leucine-rich repeat and WD repeat-containing protein 1239, KIAA1239). Furthermore, in order to resolve the deep phylogenetic (interspecific) relationships within the P. bastardi species complex, we selected seven samples representative of all species-level lineages, plus P. picta as out-group, and amplified and sequenced for each of them four additional fragments of mitochondrial genes: 12S rRNA (12S), 16S rRNA (16S), and NADH-dehydrogenase subunits 2 and 4 (ND2, ND4); as well as seven additional protein-coding nuclear genes: acetylcholinergic receptor M4 (ACM4), matrix remodelling-associated protein 5 (MXRA5), phosducin (PDC), prolactin receptor (PRLR), recombination activating gene 1 (RAG1), sacsin (SACS), and titin (TTN). Standard polymerase chain reactions were performed in a final volume of 12.5 μl containing 0.3 μl each of 10 pmol primer, 0.25 μl of total 10 mM dNTP (Promega), 0.1 μl of 5 U/ml GoTaq, and 2.5 μl of GoTaq Reaction Buffer (Promega). Primer sequences and PCR conditions are given in Appendix 2. PCR products were purified with ExoSAP-IT (Thermo Fisher Scientific, Waltham, MA, USA) and directly sequenced on an ABI 3130 capillary sequencer. New sequences were checked, corrected and trimmed with the software CODONCODE ALIGNER (CodonCode Corporation), and submitted to GenBank (www.ncbi.nlm.nih.gov/genbank/; accession numbers: MW311329–MW311370, MW315522–MW315525, MW318983–MW318987, MW319200–MW319361; cf. Appendix 3). Sequences were aligned using the Clustal algorithm in MEGA7 (Kumar et al. 2016).

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.

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 (ACM4, 380 bp; CMOS, 426 bp; KIAA1239, 869 bp; MXRA5, 795 bp; PDC, 409 bp; PRLR, 534 bp; RAG1, 1041 bp; SACS, two fragments of 975 and 1032 bp, respectively; and TTN, 849 bp); (2) an analysis based on five concatenated mitochondrial markers (12S, 1072 bp; 16S, 603 bp; ND2, 697 bp; ND4, 852 bp; and CO1, 597 bp); and (3) an analysis combining both nuclear and mitochondrial datasets (cf. Appendix 3B). These three datasets represent a total of 7310 bp, 3821 bp and 11133 bp, respectively. For most genes, the following samples were used: P. ibityensis, ZMA 19642 (FGMV 2002-0990); P. neglecta, ZSM 163/2006 (FGZC 991), P. tanjaka, ZSM 40/2006 (FGZC 750); and to represent the three lineages in P. bastardi sensu lato: ZSM 180/2004 (FGZC 332), ZSM 189/2004 (FGZC 354), and ZSM 849/2010 (ZCMV 12740). For some taxa and gene fragments, alternative samples were used, and in some cases, complemented by previously published sequences downloaded from GenBank, 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.

To compare trees obtained from mtDNA, nDNA, and the combined data sets, we calculated the Icong congruence index (de Vienne et al. 2007) with the Icong website tool (http://max2.ese.u-psud.fr/icong/index.help.html).

Mitochondrial DNA phylogenetic tree (CO1). The deepest nodes of the tree are unresolved, and the lineages considered 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. ibityensis, 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, 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) examination of the chromatograms of MirZC 40, MirZC 91, MirZC 101, MirZC 104 sequences identified the presence of an unambiguous stop codon; (2) sequences of three specimens (MirZC 40, MirZC 91, MirZC 104) for the 12S rRNA fragment clearly place them in the P. bastardi C lineage, the only lineage otherwise known from Kirindy; (3) sequences of three specimens from Isalo (MH063295: MRSN R3736; MH063296: FAZC 14631; and MH063297: MRSN R3745) for the 16S rRNA fragment unambiguously place them in the P. bastardi C lineage; and (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 lineages. 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. ibityensis sister to the (P. neglecta, 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.

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 morphology-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.

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).

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 GenBank 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. bastardi A (an adult, ZSM 189/2004; a subadult, ZSM 187/2004; and a juvenile, ZSM 188/2004) and two genotyped vouchers of P. bastardi B (a juvenile, ZSM 178/2004; and an adult, ZSM 180/2004).

Adults and subadult specimens in Tranoroa. Despite the limited sample size and lack of molecular assignment of one specimen (ZSM 58/2004), morphological qualitative examinations and multivariate statistical analyses appear to confirm the morphological differentiation between the two lineages occurring in this locality. Specimens 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.

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 and 1900.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 1899.0337, 1899.0339, 1900.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.

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), Blaesodactylus (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 divergent 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.