A new species of leaf-toed gecko (Phyllodactylidae, Phyllodactylus) from María Cleofas Island, Nayarit, Mexico

Abstract We describe a new species of leaf-toed gecko of the genus Phyllodactylus from María Cleofas Island, the smallest island of Tres Marías Archipelago, Nayarit, México. Genomic, phylogenomic, and morphological evidence support that the new species presents a unique combination of diagnostic characters. Morphologically, the new species has a high number of tubercles, head to tail (mean 47), longitudinal ventral scales (mean 61), and third labial–snout scales (mean 26). Gene flow tests revealed the genetic isolation of insular populations from mainland counterparts. In addition, we confirmed the non-monophyly of P. homolepidurus and P. nolascoensis, and we show that the taxon P. t. saxatilis is a complex; therefore, we propose taxonomic changes within the saxatilis clade. The discovery of this new insular endemic species highlights the urgency of continued exploration of the biological diversity of island faunas of Mexico.


Introduction
The family Phyllodactylidae currently includes 10 genera and 148 species of geckos (Uetz et al. 2020), with wide geographical distribution mainly in America, Africa, Europe, and the Middle East (Gamble et al. 2008;Vitt and Caldwell 2014). Phyllodactylus is among the 10 gecko genera with the highest number of species (Uetz et al. 2020). Within this diverse genus, particularly in Mexico, recent studies incorporating genomic data indicate a need for taxonomic changes and the description of new species that have remained undetected (or weakly differentiated) on the basis of morphological data alone (Ramírez-Reyes et al. 2020).
In the late nineteenth century, Stejneger (1899) reviewed the specimens of Phyllodactylus collected by E.W. Nelson and E.A. Goldman from the Islas Marías Archipelago (comprising the islands of María Madre, María Magdalena, María Cleofas, and San Juanito), and classified these populations as Phyllodactylus tuberculosus. Sixty years later, Zweifel (1960) collected additional individuals from María Madre Island and María Magdalena Island, and assigned these to the species Phyllodactylus lanei Smith, 1935. Pioneer efforts by James R. Dixon included an extensive taxonomic and systematic review of the genus Phyllodactylus in America (Dixon 1964(Dixon , 1966(Dixon , 1973Dixon and Huey 1970). Based on morphological analysis, Dixon (1964) proposed the name P. tuberculosus saxatilis for gecko populations in Western Mexico, a wide geographic area including the states of Jalisco, Nayarit, Sinaloa, Sonora, Durango, Chihuahua and some islands off the coast of Nayarit. Dixon (1964) also mentioned some morphological differences between insular (from Islas Marías) and mainland populations of P. t. saxatilis, highlighting the morphological variation in the number of dorsal tubercles, tubercles between the axilla and groin, and the number of scales bordering the postmentals. These differences, however, were not statistically significant according to his criteria. Therefore, Dixon (1964) did not believe that this morphological variation justified a taxonomic distinction between insular and mainland populations. This taxonomic classification was subsequently accepted and published in checklists of the herpetofauna of Nayarit and the Isla Marías (McDiarmid et al. 1976;Casas-Andreu 1992;Nolasco-Luna et al. 2016;Woolrich-Piña et al. 2016).
Recently, Ramírez-Reyes et al. (2020) showed that the gecko population from María Cleofas Island was isolated and had differentiated from mainland lineages during the early Pliocene. Later, Ramírez-Reyes et al. (unpublished data) found that Phyllodactylus tuberculosus is a polyphyletic taxon and that none of the three subspecies studied (P. t. tuberculosus, P. t. magnus and P. t. saxatilis) share an exclusive more recent common ancestor with each other (i.e. share MRCA with other Phyllodactylus), and therefore, the evolution of each major clade occurred independently in time and space. In the saxatilis clade, the genetic differentiation of the María Cleofas Island lineage is particularly remarkable (Ramírez-Reyes et al. unpublished data). The gene flow or introgression test (ABBA-BABA) showed that this insular lineage has remained genetically isolated from mainland populations within the saxatilis clade. This last characteristic suggests treating this lineage as a candidate species under the biological species concept. However, there is uncertainty about the phylogenetic relationships of the insular lineage with respect to the continental lineages of the saxatilis clade. According to phylogenetic approximations (ML and Bayesian) based on the supermatrix (or supergene) of genomic data, the insular lineage of the saxatilis clade is related to Phyllodactylus lineages from Sinaloa (Ramírez-Reyes et al. 2020;Ramírez-Reyes et al. unpublished data). A species tree approach with coalescence methods could provide greater certainty about the phylogenetic relationships of lineages and species within the saxatilis clade. Given the set of evidence that allows defining this new species, the main objective of this study is the formal description of the insular lineage from María Cleofas Island, which adds an endemic taxon to the herpetofauna of the Tres Marías Archipelago. In addition, we elaborate a coalescence-based species tree for the saxatilis clade from SNP data, and assess the taxonomic treatment of the lineages and candidate species within the saxatilis clade.

Study area
María Cleofas Island was well described by Emerson (1958), Zweifel (1960), andCasas-Andreu (1992). Of the four islands that comprise the Islas Marías Archipelago Biosphere Reserve, María Cleofas Island is the closest to the mainland, at a distance of 87.5 km to the nearest continental point located in San Blas, Nayarit, and has a surface area of 25 km 2 (Fig. 1). Tropical sub-deciduous and deciduous forests, scrub forests, and a small area of Rhizophora mangle dominate the vegetation (Fig. 2). The climatic conditions are characterized by a dry season from November to April and a rainy season from May to October (Bullock 1986). The island shows an average ambient temperature of 24.9 °C; the minimum and maximum temperatures are 21.1 °C and 28.7 °C in January-February and July, respectively. The mean annual rainfall is 564.2 mm, of which 95% occurs from June to December (CONANP 2007).

Sampling
We performed diurnal and nocturnal surveys in all available habitats to assess the herpetofauna of the María Cleofas Island from May 2017 to August 2018. The survey consisted of 2-3 researchers walking ad libitum through the different types of vegetation on the island, examining every microhabitat that could be potentially utilized by amphibians and reptiles (Fig. 3). Our field observations were conducted during the day from 09:00 to 17:00 h, and at night from 21:00 to 04:00 h. All captured lizards were taken to the María Cleofas Research Boat for sampling and data collection (see below for morphological analyses). We captured 36 individuals of the genus Phyllodactylus during fieldwork under collecting permits SGPA/DGVS/01208/17 and SGPA/ DGVS/010144/18, of which tissue samples of two specimens were used for molecular analyses (Ramírez-Reyes et al. 2020), and 11 specimens were used for the present morphological description. We only sacrificed nine individuals (the type series) by injection of sodium pentobarbital, following the animal ethics guidelines of UNAM.

Genomic data and species tree
We used the generated variant calling format file (VCF) of Ramírez-Reyes et al. (unpublished data), which contains the genetic variants (SNPs) corresponding to the saxatilis clade. This VCF file contains 13 genotypes and 14,271 binary SNPs (8.1% missing data). We analyzed this VCF file in R, mainly using three packages: vcfR (Knaus and Grünwald 2017), adegenet (Jombart and Ahmed 2011), and StAMPP (Pembleton et al. 2013). With these packages, we calculated Nei's distances between populations, as well as pairwise F ST among populations. We visualized Nei's distances as a Neighbor joining (NJ) tree obtained in Rstudio.
We elaborated a species tree in SVDquartets Kubatko 2014, 2015) implemented in PAUP* v. 4 (Swofford 2003) with the 14,271 binary SNPs. The individual samples (n = 13) from the VCF file were assigned to seven species, based upon previous studies about evolution and candidate species in the saxatilis clade (Ramírez-Reyes et al. 2020;Ramírez-Reyes et al. unpublished data). We performed a species tree search using exhaustive quartet sampling (289 quartets) and used the QFM assembly algorithm to assemble the species tree. We also indicated the multispecies coalescent tree model. Branch support was assessed using 100 nonparametric bootstrap replicates.

Morphological analyses
We collected meristic character counts (n = 104) and morphometric measurements (n = 117) of specimens with an electric digital caliper (resolution 0.00005"/0.01 mm). We examined 117 specimens from four scientific collections: Colección Nacional de Anfibios y Reptiles at Instituto de Biología (CNAR-IBUNAM), Museo de Zoología de la Facultad de Ciencias at Universidad Nacional Autónoma de México (MZFC-UNAM), Colección de Anfibios y Reptiles at Escuela Nacional de Ciencias Biológicas All measurements and counts followed the protocol of Ramírez-Reyes and Flores-Villela (2018). Measured characters were: SVL snout-vent length, AG axilla-groin length, HW head width, SE snout-eye length, ED eye diameter, AO auricular opening, IOD interorbital distance, IND internarial distance, LFF length of the fourth finger, and LFT length of the fourth toe. The meristic characters counted were: DLF digital lamellae of the fourth finger and toe (left and right, four total counts), THT tubercles from head to tail, TAG tubercles from axilla to groin, RTD rows of tubercles across the dorsum, IS interorbital scales, SBI scales bordering internasals, TSS third labial-snout scales, SEN scales between nostril and eye, SBP scales bordering postmentals, SAV scales across venter, and LVS longitudinal ventral scales.
We conducted multivariate analyses for morphological characters using PCA and MANOVA to determine significant morphological differences among species within the saxatilis clade. A one-way ANOVA was performed to compare the mean SVL of all species analyzed, as well as the differentiated characters suggested by Dixon (1964): number of dorsal tubercles, tubercles between axilla and groin, and number of scales bordering the postmentals. Finally, the meristic characters were analyzed by multivariate analysis (MANOVA). All statistical analyses (morphometric and meristic characters) were performed in Rstudio.

Diet
We analyzed the stomach content from 36 individuals during fieldwork using stomach flushing (Legler and Sullivan 1979). In contrast with dissection techniques that require sacrificing the animals to obtain stomach contents, the aforementioned technique enables individuals to regurgitate prey items for further analysis (Perkins and Eason 2019). We performed this technique three times consecutively to obtain the largest possible volume of stomach contents from each individual (Barraza-Soltero 2019; Barraza-Soltero and Escobedo-Galván 2020). The captured geckos were taken to the María Cleofas Research Boat for sampling and data collection. After data collection, geckos were under observations between 24 to 36-hours before being returned to the capture sites on the island. During this time, we observed no deaths due to manipulation. Stomach contents were preserved in 70% ethanol for laboratory analysis at the Biodiversity and Ecosystem Services Laboratory of the Universidad de Guadalajara, Puerto Vallarta, Jalisco, Mexico. The stomach contents were separated using a Carl Zeiss Stemi DV4 stereoscopic microscope and were examined under an Olympus optical microscope. Prey items were identified to the lowest taxonomic level possible using specialized literature. In addition, confirmation and identification of some prey species was done by Fabio G. Cupul-Magaña (CUC, Universidad de Guadalajara, México).

Genomic evidence and species tree
The F ST values supported a very high genetic differentiation among the populations studied (F ST = 0.82 on average among all populations; Table 1). In particular, for P. t. saxatilis from the María Cleofas Island population (the focal population in this study), F ST was 0.92, on average, with respect to all other populations. The lowest F ST value (0.36) was calculated between the populations of P. t. saxatilis from Mocuzari and P. homolepidurus (mainland). Nei's distance values were consistent with F ST values; we obtained a distance of 5% on average among all populations, based on SNP data. These differences among individuals and populations were clearly supported and illustrated in the NJ tree (Fig. 4). The NJ topology showed that the most differentiated populations of the saxatilis clade were those from María Cleofas Island and Villa Union. The species tree developed in SVDquartets showed a high percentage of compatible quartets (87.5%) versus incompatible quartets (12.5%). The recovered species tree shows that the mainland populations of P. t. saxatilis (North and South) are rendered non-monophyletic by P. nolascoensis, P. homolepidurus, and P. partidus. This result is broadly consistent with previous studies, in which the non-monophyly of P. t. saxatilis This fundamental difference between populations in southern Sinaloa with respect to those in Sonora (both previously considered as P. t. saxatilis) requires the delimitation of P. saxatilis in order to avoid non-monophyletic groups. We consider that P. saxatilis be restricted to the populations from the type locality (Villa Union, Sinaloa) and surrounding areas (Dixon 1964). Thus, P. cleofasensis sp. nov. (see description) is related to P. saxatilis, while the other insular species P. nolascoensis is related to Phyllodactylus sp. and the two sister species (P. homolepidurus + P. partidus) (Fig. 5).

Morphological analyses
The results of the meristic counts and morphological measurements are summarized in Tables 2 and 3, respectively. All statistical analyses are based on adult individuals. The ANOVA results for the four characters analyzed rejected the null hypothesis of equality of means among the six compared species with high significance (Table 4). Of the four characters analyzed, SVL obtained the highest significance and was supported by Tukey's paired test, of which we highlight the differences between P. magnus-P. cleofasensis sp. nov. (P = 0.009), P. nolascoensis-P. magnus (P < 0.001), P. partidus-P. magnus (P < 0.001), Phyllodactylus sp. (Son)-P. magnus (P < 0.001), and P. saxatilis-P. partidus (P = 0.089). Regarding RTD, Tukey tests showed a greater difference between Phyllodactylus sp. (Son) -P. cleofasensis sp. nov. (P = 0.01). The AGT Tukey tests supported the differences between P. partidus-P. cleofasensis sp. nov. (P = 0.01), P. saxatilis-P. cleofasensis The PCA of morphometric data (excluding SVL) for the six species studied (Fig. 6) showed that about 83% of the variance was explained by the first three components. However, plotting individuals by the main components does not recover any group independently; furthermore, it shows that high overlap exists in the morphometric characters analyzed (Fig. 6). However, as mentioned previously, SVL is the most informative morphometric character; we can even visualize the differences in size of the studied species. For the saxatilis clade, P. cleofasensis sp. nov. and P. saxatilis have the largest body size, while the insular P. partidus and P. nolascoensis are the smallest (Fig. 5). For the meristic characters, we performed the MANOVA excluding the SBI, TSS, and SEN characters since they did not show a normal distribution (or close to it). All MANOVA tests rejected the equality of multivariate means for all species analyzed (P < 0.001), Pillai trace (Trace = 1.51, d.f. = 5, F = 3.6, P = 1.73 × 10 -14 ), Wilk's lambda (λ = 0.14, F = 3.9, P = 5.715 × 10 -16 ), Hotteling test (HL = 2.66, F = 4.2, P = 2.2 × 10 -16 ) and Roy test (R = 1.26, F = 10.5, P = 2.99 × 10 -12 ). Once we determined that the multivariate means were not the same, we proceeded to perform individual ANOVA on paired comparisons of the most significant characters (P < 0.001), but emphasizing the differences of P. cleofasensis sp. nov. with respect to the others. For example, in the case of THT, P. cleofasensis sp. nov. differs significantly with respect to all species (P < 0.001) except P. nolascoensis, a pattern also seen with TAG and RTD characters. On the other hand, the ventral longitudinal scales demonstrated a significant difference between P. cleofasensis sp. nov. and P. magnus (P < 0.001), P. nolascoensis (P < 0.001), P. partidus (P = 0.009) and P. saxatilis (P = 0.002), while the lowest significance value for this character was with Phyllodactylus sp. (P = 0.04).

Taxonomy
The present study constitutes a precedent for the taxonomy of the saxatilis clade (sensu Ramírez-Reyes et al. unpublished data) based on analyses of genomic and phylogenomic data. Here, our taxonomic treatment allows the use of specific categories for the studied populations and avoids non-monophyletic groups and infraspecific categories (P. t. saxatilis, P. homolepidurus). Specifically, in the case of the endemic population from the María Cleofas Island, the genomic evidence suggests that this population is clearly differentiated from its mainland counterparts (Ramírez-Reyes et al. 2020;Ramírez-Reyes et al. unpublished data), and a formal description of this species is warranted.
Description of the holotype. Adult female with SVL 71.17 mm, robust body, head not flattened, neck slightly differentiated from head. Head width 14.4 mm and snout length 12.11 mm. Rostral scale is flat (no grooves or stretch marks) and in contact with the two internasal scales. Nostril in contact mostly with the rostral scale and marginally with the first labial scale on both sides; supranasal scales in contact with 12 scales cross- ing from right to left side; 22 interorbital scales counted from middle of eye, interorbital distance 9.75 mm. 24 scales crossing snout between contralateral second labial scales 27 between third labial scales. Number of loreal scales 14 on right side and 15 on left side; equal number of supralabial scales (14), labial scales (8), and infralabial scales (7) on each side. Auricular opening oval (2.48 mm) smaller than the ocular opening, occipital scales similar in size and shape to interorbitals (not greatly differentiated in size and shape). Mental scale slightly wider than long, forming a "V" but not with pronounced angles; eight postmental scales in contact with first infralabials on both sides (right and left). Body with granular and circular scales interspersed with tubercles of different sizes; the specimen presents a fragmented tail. Internarinal distance 2.67 mm and axilla-groin length 26.19 mm. 14 rows of dorsal tubercles, 20 axilla-groin tubercles on right side and 19 on left side. Presents 29 transversal ventral scales with first ventral scale differentiated from lateral scales (which are small and circular). Ventral scales differentiated in size and shape from lateral and gular scales. Scales imbricate on extremities (anterior and posterior), as well as on dorsal region of the tail. No femoral or precloacal pores. Digital lamella formulae: right posterior (9-11-15-14-13), left posterior (8-10-11-12-13), right anterior (8-11-12-13-10), left anterior (8-10-12-13-11); fourth finger of extremities longer that others (5.55 mm manus, 7.49 mm pes); digital toepads longer than wide on all fingers.
Etymology. Specific epithet is taken from the type locality María Cleofas Island, with the Latin suffix -ensis meaning, "originating from." Specific epithet is masculine, in agreement with the gender of Phyllodactylus.
Variation. All meristic and morphometric characters are presented with mean values and standard deviation in Tables 2 and 3, respectively.
Natural history. Individuals of P. cleofasensis were observed active during night surveys under single rocks, abandoned anthropogenic structures, or in some cases on the trunk of Piranhea mexicana. Some were observed on the ground while they moved between rocks. In some parts of the island, they can be seen in abandoned man-made structures, which are shelters for species such as geckos, anoles, iguanids, and bats. Although predation of P. cleofasensis has not been reported, we suggest that Oxybelis microphthalmus (for taxonomic status see Jadin et al. 2020) could be a potential predator given that we captured two O. microphthalmus at the same site of P. cleofasensis during nocturnal surveys. During our visits (May to August 2017 and April to August 2018), we did not observe reproductive activity or females with a shelled egg in one of their oviducts. Eleven distinct prey items were found in 36 stomach contents. It has been mentioned that the diets of lizards in island environments contain a high percentage of vegetation matter due to the lack of food or scarcity of prey belonging to the Class Insecta (Olesen and Valido 2003). However, in this work, a great variety of arthropods was reported in the stomach contents of geckos (Fig. 8). Prey items of P. cleofasensis were mostly composed of orthopterans (Family Rhaphidophoridae) and coleopterans (Family Passalidae); they also consumed plant matter, arachnids, lepidopterans, scorpions (Centruroides elegans insularis), and cockroaches (Blattodea: Pycnoscelus surinamensis) (Fig. 8). Additionally, during fieldwork in May 2018, we captured two individuals with remains of shed skin in their stomach contents (Barraza-Soltero and Escobedo-Galván 2020).

Discussion
Herein we present genomic, phylogenomic, and morphological evidence that demonstrates the independent species status of the endemic population of Phyllodactylus from María Cleofas Island. To our knowledge, the only work that mentioned morphological differences between the insular and mainland populations of P. t. saxatilis was conducted by Dixon (1964). Fifty-six years later, we confirmed that these differences were supported by genetic isolation of this population from those on the mainland using genomic data (Ramírez-Reyes et al. 2020;Ramírez-Reyes et al. unpublished data). Surprisingly, Phyllodactylus cleofasensis is the third gecko species to be recognized from the islands located off the coast of Nayarit, together with P. isabelae endemic to the Marietas Islands and P. lupitae endemic to El Coral Island (Ramírez-Reyes and Flores-Villela 2018). In the 1960's, there was a consensus among herpetologists that the insular amphibians and reptiles from the Western Mexican Pacific were representatives of taxa found on the adjacent mainland, due to a possible Neotropical influence in their evolution (Casas-Andreu 1992). However, recent molecular techniques have challenged the traditional view of the herpetofauna on both the mainland and the islands of the Central Pacific of Mexico. For example, Card et al. (2016) recently recognized independent lineages of the genus Boa for the Mexican Pacific coast. Similarly, Ramírez-Reyes et al. (2017) analyzed and confirmed divergent lineages within the Phyllodactylus lanei complex. In light of this diversity, increased research efforts in alpha taxonomy should be considered in conjunction with the generation of genetic and genomic resources of insular fauna. The scarcity of current studies reflects minimal interest in the conservation of the insular diversity of Mexico, and at the same time maintains the great uncertainty of the status quo of the herpetofaunistic diversity of insular ecosystems (Blázquez et al. 2018).
The phylogenetic position of P. cleofasensis was unknown until recently. Based on a relaxed molecular clock model, this species diverged from mainland populations during the late Miocene (~7 mya), and shows the highest genetic isolation of the species in the saxatilis clade (Ramírez-Reyes et al. unpublished data). The recovered species tree allows us to propose a more stable taxonomy for this clade, with a clear separation of P. homolepidurus and P. nolascoensis (as was recently proposed by Lemos-Espinal et al. 2019), and some non-monophyletic lineages that were considered part of P. t. saxatilis (Phyllodactylus sp. in this study). Only a broader sampling will establish if there are significant morphological differences between P. homolepidurus and P. nolascoensis, as mentioned by Dixon (1964). Meanwhile, our taxonomic treatment allows the use of specific categories for the studied populations and avoids polyphyletic groups and infraspecific categories (P. saxatilis and P. homolepidurus). Similar to this study, Koch et al. (2016) implemented an integrative taxonomic framework to propose an updated and stable taxonomy for the P. reissi complex in Perú.
The three insular endemic species from Nayarit are differentiated morphologically. Phyllodactylus lupitae (SVL = 64.25) is the largest (mean body size), though very similar in size to P. cleofasensis (SVL = 64.19). The smallest species is P. isabelae (SVL = 45.91). Phyllodactylus lupitae is differentiated from P. cleofasensis by a smaller number of dorsal paravertebral tubercles (mean 28.82 versus 47.89). With respect to the continental species compared, we demonstrate that a set of diagnostic characters (THT, TSS, and SLV) allows a clear morphological differentiation of P. cleofasensis from the other Phyllodactylus species in Mexico (see Table 1 and Diagnosis).
The stomach contents of P. cleofasensis coincide with those reported in other studies on Phyllodactylus species. We found coleopterans, arachnids, and orthopterans, as well as unidentifiable matter (vegetation, rocks, and insect wings) and shed skin remains. Castro-Franco and Gaviño (1990), on El Coral Island, carried out the only study on the diet of the genus Phyllodactylus in insular conditions; they identified coleopterans, hemipterans, hymenopterans, homopterans, orthopterans, and dipterans, as well as larvae of different insects, as part of the diet of P. lupitae. In this study, we found Araneae, Coleoptera, and Orthoptera present in for the dietary habits of the gecko population from María Cleofas Island, as well as Nematoda, Lepidoptera, Isoptera, Scorpiones, and Blattodea. Barraza-Soltero (2019) identified the stomach contents from the lizard community on María Cleofas Island, including Anolis nebulosus, which consumed almost the same number of prey types as P. cleofasensis despite differences in foraging habits. In addition, the identification of shed skin remains in stomach contents from some individuals suggests for the first time keratophagous behavior in insular species from the Central Pacific islands of México.
Finally, a striking result was that we did not find non-native species during explorations on the island, despite María Cleofas Island having suffered ecological disturbances due to anthropogenic activities and the vestiges of man-made constructions visible. To our knowledge, non-native species have been reported on some islands pertaining to Nayarit, specifically the Common gecko Hemidactylus frenatus on María Madre Island (Casas-Andreu 1992), Isabel Island (Valdez-Villavicencio and Peralta-García 2008), and El Coral Island (Ramírez-Reyes et al. 2015). Similarly, some non-native species have been reported along the Pacific Coast of Mexico in recent years, such as the Brown anole (Anolis sagrei, Pazos-Nava et al. 2019), the Mourning Gecko (Lepidodactylus lugubris, Ahumada-Carrillo and Weatherman 2018), and possibly a gecko species of the genus Gehyra (IKB-S pers. obs.). These species are human commensals, frequently transported intentionally or unintentionally to new localities. Therefore, it is urgent to establish a protocol of species monitoring on María Cleofas Island to prevent the entry of non-native species, which may be a direct threat to native species and disturb balance of insular ecosystems.