Unexpected Morphological Diversity in New Zealand’s Large Diplodactylidae Geckos

Prehistoric anthropogenically-mediated extinctions have impacted global biodiversity; however effects on herpetofauna are poorly-documented. New Zealand’s Diplodactylidae geckos exhibit high species-level diversity, largely independent of discernible osteological changes (cryptic). Consequently, taxonomic anities of isolated skeletal elements (fossils) are primarily determined by relative size, particularly in the identication of Hoplodactylus duvaucelii; New Zealand’s largest extant gecko species. Here, three-dimensional geometric morphometrics of maxillae (a common fossilized element) was used to determine whether consistent shape and size differences exist between genera, and if cryptic extinctions have occurred in ‘Hoplodactylus cf. duvaucelii’. Sampling included 13 Diplodactylidae species from ve genera, and 11 Holocene ‘H. cf. duvaucelii’ subfossil individuals. We found phylogenetic history was the most important predictor of maxilla morphology among extant Diplodactylidae genera. Relative size comparisons could only differentiate Hoplodactylus from other genera, with the remaining genera exhibiting variable degrees of overlap. Six subfossils were positively identied as H. duvaucelii, conrming their proposed Holocene distribution throughout New Zealand. Conversely, ve subfossils showed no anities towards any modern Diplodactylidae genera, implying either increased morphological diversity in mainland ‘H. cf. duvaucelii’ or the presence of at least one extinct, large, broad-toed Diplodactylidae species. These results highlight the impact of anthropogenic disturbances on insular reptile diversity.

Herein, three-dimensional geometric morphometrics was used to characterise and quantify both shape and size variation in the maxilla of modern Diplodactylidae genera (Dactylocnemis, Hoplodactylus, Mokopirirakau, Naultinus and Woodworthia), for comparison with Holocene 'H. cf. duvaucelii' subfossils. Three main research questions were tested: (a) can revised Diplodactylidae genera be distinguished based on maxilla shape; (b) is relative-size a reliable  Table 4) across all PC axes suggest shape similarities with Dactylocnemis (E, K), Hoplodactylus (A, C, D, G, I, J) and Woodworthia (B, F, H), with no a nities towards Mokopirirakau or Naultinus.

(b) Predictors of shape and size
Procrustes ANOVA (Supplementary Table 5) revealed that phylogenetic a liation (i.e. genus) is a highly signi cant predictor (F (4,38) = 9.01, p < 0.001) of maxillae shape, accounting for 45.2% of the shape variation. Multivariate pairwise post-hoc tests found differences to be signi cant between most genera (p < 0.05), excluding Dactylocnemis-Hoplodactylus (p = 0.229), and Hoplodactylus-Mokopirirakau (p = 0.056) comparisons (Supplementary Table 6). A weak but signi cant relationship also exists between maxillae shape and centroid size (F (1,41) = 5.39, p = 0.020), and their interaction (F (4,38) = 1.35, p = 0.023; Supplementary Table 5), suggesting a small proportion of the shape diversity (6.8%) is due to allometry. Table 7  Phylogenetic position is a highly signi cant predictor of maxilla shape diversity in New Zealand Diplodactylidae, with all genera (Dactylocnemis, Hoplodactylus, Mokopirirakau, Naultinus and Woodworthia) being morphologically distinct. These results contrast previous long-held notions of skeletal conservatism in New Zealand's geckos (e.g. (27,37)) through identi cation of taxonomically informative morphological variation within a single skeletal element. This retention of genus-level phylogenetic signal is remarkable given pronounced ecological species radiations since the early Miocene (14). However, similar trends of reduced disparity in maxillae (relative to rate of evolution) are observed across both extant and extinct squamates (excluding snakes; (38)), suggesting constrained evolution in this cranial region.

One-way ANOVA (Supplementary
Diplodactylidae maxilla shape is predominantly characterized by two character-states, described by the rst two axes of both PCA and CVA: (1) posterior extension/reduction of the nasal margin; and (2) increase/decrease in dorsoventral extent of the facial process. The separation of genera along PC1 appears to re ect broad habitat use of the New Zealand Diplodactylidae, with terrestrial-arboreal (Dactylocnemis, Hoplodactylus and Woodworthia) and exclusively arboreal (Naultinus) genera occupying positive and negative regions respectively (39,40). This morphological signature of habitat use extends to species-level comparison, most notably in the discrimination of the terrestrial-arboreal M. 'southern North Island' from the arboreal M.granulatus (41), characterized by a shift to more positive values.
In gekkotans, arboreal forms tend towards broad, pointed and dorsoventrally shallow skulls, enabling faster climbing speeds on non-horizontal surfaces (42,43). While cranial modi cations associated with habitat use are undocumented in the New Zealand Diplodactylidae, extension of the nasal margin in arboreal species appears to be linked to two super cial morphological changes in the adjacent prefrontal margin: (1) a reduction in anterior extent (observed in other Gekkota; (44)); and (2) formation of a thickened ridge along the prefrontal orbital margin (Supplementary Figure 7). While the function of these structures remains unclear, association with arboreality provides strong evidence for ecomorphological convergence between phylogenetically independent lineages.
In addition to habitat use, skull-shape evolution in lizards is strongly in uenced by diet, with shape variation concentrated in the premaxilla, nasal and jaw joint, re ecting their roles in rostral prey capture and feeding biomechanics (38,42). Herbivorous lizard skulls tend towards reduced snouts and high temporal regions relative to carnivorous lizards, contributing to an increased bite strength required for processing brous and tough foliage (45)(46)(47). Conversely, omnivorous gekkotans represent intermediate forms not specialized to particular feeding behaviors, and consequently lack unique morphological adaptations (48). New Zealand geckos are predominantly omnivorous, consuming a wide variety of food items including plant matter (fruit, honeydew and nectar) and arthropods (40). Such extensive dietary overlap effects the performance of diet as an explanatory variable of maxilla shape diversity, given categories (omnivorous and insectivorous) are not discrete.

(b) E cacy of size-based discrimination
Maxilla size was signi cantly correlated with phylogenetic a nity, however, only Hoplodactylus could be fully differentiated (under post-hoc comparisons), with the remaining Diplodactylidae genera exhibiting variable degrees of overlap. This highlights the ine ciency of previous size-based taxonomic identi cation of non-Hoplodactylus Holocene subfossil geckos, especially intermediate-sized genera (Dactylocnemis, Mokopirirakau and Naultinus), which exhibit complete size overlap. Similarly, while large relative size proves reliable in discriminating extant H. duvaucelii, applications in Holocene subfossil identi cation are limited given assumptions of temporal taxonomic homogeneity (or "covert biases"; (49)).

(c) Increased Holocene diversity of large geckos
Our results provide evidence for increased morphological diversity of large geckos during the Holocene in New Zealand, with declines in both shape and size variation following Polynesian and European colonization.
Combined Procrustes and Mahalanobis distance comparisons provide support for previous size-based classi cation of ve Holocene subfossils (A, D, E, J, K) as H. duvaucelii, con rming assumed prehuman distribution across both the North and South Islands. The remaining six Holocene subfossil specimens (B, C, F, G, H, I) exhibited classi cation discrepancies and/or reduced assignment probabilities (below relevant thresholds), re ected in their unique position across CV1/CV2. These distinct Holocene subfossil maxillae ("unknown taxa") are not re ective of differential adaptation to mainland and island habitats (see above), therefore re ecting either increased morphological diversity of mainland large species (not encompassed by extant populations) or the presence of at least one extinct, large, broad-toed Diplodactylidae species.
Based on digit morphology, the extinct giant H. delcourti was positioned within the broad-toed clade, sister to H. duvaucelii (24), suggesting these "unknown taxa" could potentially represent small or even juvenile H. delcourti (with respect to the latter hypothesis). However, this seems unlikely given the paucity of reported subfossil remains of H. delcourti (56), despite extensive collections of other Diplodactylidae taxa (37). Accurate phylogenetic a nities of both H. delcourti and "unknown taxa" could be determined through future ancient DNA analysis.
During the Holocene, mainland H. duvaucelii (and "unknown taxa") reached larger sizes than extant populations, re ected in a reduction in maximum maxilla size (a proxy for body size; e.g. (50)). Such sized-biased extinction is well-documented in Quaternary lizards globally (51,(57)(58)(59), including the extinction of two large-bodied Eugongylinae skink species (Oligosoma northlandi and Oligosoma sp.) in northern New Zealand (12,31,60). This re ects the inherent vulnerability of New Zealand's large-bodied, nocturnal herpetofauna towards high-predation rates and ecological displacement by exotic mammals (including the Paci c rat (kiore); (61,62)), particularly in forest-cleared environments (63). Smaller lizards can escape predation during periods of inactivity through utilizing narrow retreats, given limited overlap in body diameter with small mammalian predators (39). Conversely, refugia utilized by large-bodied lizards can be accessed by mammalian predators, evidenced by reductions in body weight, tail width and recruitment of H. duvaucelii on kiore-inhabited islands (35,64).
Similar to extant H. duvaucelii populations (65), Holocene subfossil H. duvaucelii also exhibit a latitudinal cline in maxilla size opposing Bergmann's rule (i.e. increased size at high latitudes), with individuals from northern localities being noticeably larger than those from southern localities. For diurnal lizards, reduced body size appears to be an advantageous thermoregulatory strategy in cooler climates, with high surface-area to volume ratio permitting rapid heat gain whilst sun-basking (66,67). Despite being nocturnal, H. duvaucelii occasionally emerge from retreats to thermoregulate through cryptic sun-basking (68,69), suggesting small body size provided an adaptive advantage at high latitudes.

Conclusions
New Zealand Diplodactylidae genera can be fully differentiated based on maxilla shape, which exhibits strong correlations with phylogenetic history. Additional species-level discrimination based ecomorphological adaptations highlights the potential application of geometric morphometrics to more functionally variable elements (or whole skulls) in taxonomic descriptions of extant Diplodactylidae species. Previous sized-based identi cation of Holocene subfossils is ineffective and grossly underestimates extinct diversity, suggesting global assemblages of insular reptiles are depauperate in comparison to prehuman diversity.

Methods (a) Specimen Selection
To capture extant morphological variation, we examined both left and right maxillae (sensu (70) Shape variation in maxillae of the extant species was assessed using principal component analysis (PCA); with intergeneric differences (shape ~ genus * size) tested using a Procrustes analysis of variance (ANOVA; (73)), and visualized using canonical variate analysis (CVA; (74)) with cross-validations, based on a reduced set of PC scores (75,76). Three-dimensional surface warps (77) representing minimum and maximum shapes along both principal component (PC) and canonical variate (CV) axes were generated using the thin-plate spline (TPS) method (76,78).
Holocene subfossil maxillae were then projected into these two-dimensional morphospaces (i.e. PCA and CVA) through matrix multiplication with respective eigenvectors (e.g. (79)). Phylogenetic classi cation of Holocene subfossil specimens was performed through Procrustes and Mahalanobis distance comparisons (to the mean maxilla shape of each genus), with the latter used to calculate typicality (80,81) and posterior (82) probabilities Variation in size of the maxilla (represented as centroid-size of the landmark con guration) between genera was examined using a one-way ANOVA and Tukey's honestly signi cant difference (HSD) post-hoc tests (83), and visualised using a barplot. All statistical analyses were performed in the R statistical environment v. 3