Skeletal variation in extant species enables systematic identification of New Zealand’s large, subfossil diplodactylids

New Zealand’s diplodactylid geckos exhibit high species-level diversity, largely independent of discernible osteological changes. Consequently, systematic affinities of isolated skeletal elements (fossils) are primarily determined by comparisons of size, particularly in the identification 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 subfossil ‘Hoplodactylus cf. duvaucelii’. Sampling included 13 diplodactylid species from five genera, and 11 Holocene subfossil ‘H. cf. duvaucelii’ individuals. We found phylogenetic history was the most important predictor of maxilla morphology among extant diplodactylid genera. Size comparisons could only differentiate Hoplodactylus from other genera, with the remaining genera exhibiting variable degrees of overlap. Six subfossils were positively identified as H. duvaucelii, confirming their proposed Holocene distribution throughout New Zealand. Conversely, five subfossils showed no clear affinities with any modern diplodactylid genera, implying either increased morphological diversity in mainland ‘H. cf. duvaucelii’ or the presence of at least one extinct, large, broad-toed diplodactylid species.


Background
New Zealand's lizard fauna is characteristic of isolated archipelagos, exhibiting high species endemism, extensive in situ radiations [1,2] and insular gigantism [3]. Despite this diversification, the New Zealand Diplodactylids (compared with the New Caledonian radiation) are relatively conserved in form; exhibiting reduced morphological divergence at a similar evolutionary depth [2,4].
Formal taxonomic descriptions of New Zealand diplodactylid species [5,6], similar to those of Australian [7,8] and New Caledonian representatives [9][10][11], have been based exclusively on external morphological characters (e.g. coloration and scalation), with interspecific skeletal variation rarely analysed. Early anatomical studies [12], however, noted osteological differences between the two then recognized genera: Hoplodactylus ('brown geckos') and Naultinus ('green geckos'); focussing on the neotenic condition of Naultinus and the 'primitive' osteology of the New Zealand Diplodactylidae. Despite these indepth comparisons, re-examination of generic level skeletal variation in the context of current nomenclature is required, given considerable taxonomic fluidity over the last 65 years [13][14][15].
Descriptions of these revised genera were based exclusively on external morphology [2]; with osteological differences only recognized for the frontal [21,22], which has led to assertions of skeletal uniformity at the specieslevel [23,24]. In comparison, the identification of interspecific skeletal variation in cranial elements of New Zealand's large eugongyline skinks has enabled both differentiation from smaller congeners [25], and description of the extinct 'giant' Oligosoma northlandi [26,27] from Holocene subfossil remains. Classification of isolated diplodactylid remains however, has been restricted to size comparisons with extant representatives (in reference to now outdated taxonomy; [23]), particularly in the identification of 'H. cf. duvaucelii' , New Zealand's largest extant diplodactylid.
Hoplodactylus duvaucelii ('Duvaucel's gecko') is a large, nocturnal species, with a pseudoendemic (realized) distribution on predator-free islands in the Cook Strait and off the north-eastern coast of the North Island ( Fig. 1A; [28]). Prior to Polynesian (~ 1280 AD; [29]) and European (effectively the late 1700s) arrival, 'H. cf. duvaucelii' was widely distributed throughout the North Island [25], and the northwest and eastern South Island ( Fig. 1a; [23,[30][31][32][33]), evidenced by the presence of large, isolated diplodactylid remains in Holocene predator (laughing owl and falcon) middens. Subsequent range contractions occurred as a result of the synergistic effects of competitive exclusion, direct predation by introduced mammals, and degradation of forest habitat [25,34]. This extensive distribution across multiple biogeographic regions [23,35], given pronounced regional endemism recognised (or proposed) in other New Zealand diplodactylid genera [1,2], combined with the extinction of other large lizards Surface models of a diplodactylid maxilla in lateral (top), dorsal (middle) and medial (bottom) views demonstrating placement of fixed landmarks (black circles) and equally spaced semilandmarks (white circles). Numbers and C-prefixed numbers correspond to anatomical landmark descriptions (Additional file 1: Table S2. 3) in New Zealand [26,27] implies unrecognized diversity may exist within 'H. cf. duvaucelii'; perhaps detectable through fine-scale osteological analysis. Geometric morphometrics [36], a method of statistical shape analysis that enables improved detection and visualisation of subtle morphological differences (compared with traditional linear-based morphometrics [37,38]), has been widely applied in herpetofaunal studies; including the successful discrimination of closely-related species [39][40][41] and classification [42][43][44] of isolated cranial elements. We therefore predict osteological differences, if examined appropriately (using geometric morphometrics), will be sufficient for discriminating between extant diplodactylid genera (and potentially species), enabling the identification of isolated subfossil material, or reveal unidentified taxa that require description and diagnosis.
Herein, three-dimensional geometric morphometrics is used to characterise and quantify both shape and size variation in the maxilla of modern New Zealand diplodactylid genera (Dactylocnemis, Hoplodactylus, Mokopirirakau, Naultinus and Woodworthia), for comparison with Holocene 'H. cf. duvaucelii' subfossils. Three main research questions are tested: (a) can recognised diplodactylid genera be distinguished based on maxilla shape; (b) is size a reliable method for generic-level identification of isolated cranial elements; and (c) have cryptic extinctions occurred in the Diplodactylidae (with a focus on 'H. cf. duvaucelii')?
PC3 (8.0% of variance) describes shape differences in both the anterolateral lappet and prefrontal margin (Additional file 1: Figure S4); with dorsoventral thickening (and lateral thinning) of the anterolateral lappet and plateauing of the prefrontal margin at negative values. Conversely, at positive values, the prefrontal margin forms a near continuous curve with the adjacent orbital margin (Additional file 1: Figure S4). Shape differences along PC4 (7.2%) are primarily associated with increased curvature of the tooth row along the negative sector of the axis (Additional file 1: Figure S4).
Visually, Holocene subfossil specimens cluster in the intermediate region of the positive zones of PC1, PC2 and PC4; overlapping the morphospaces of extant genera ( Fig. 2a; Additional file 1: Figure S3). Conversely, Holocene subfossil specimens (excluding H) occupy increasingly positive regions of PC3, with some individuals (B, E, I, J) exhibiting no overlap with extant genera (Additional file 1: Figure S3). Procrustes distances of the Holocene subfossil specimens (Additional file 1: Table S4) 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 shape similarities with Mokopirirakau or Naultinus.

Phylogenetic shape differences
Phylogenetic signal in maxilla shape is statistically significant (K mult = 0.828, p < 0.001); however, species resemble each other less than expected under a model of Brownian motion (given K mult is less than 1). This is reflected in the distribution of species throughout the phylomorphospace (Additional file 1: Figure S6), with several overlapping branches (e.g. Dactylocnemis and Hoplodactylus) and non-adjacent closely related species (e.g. M. granulatus and M. 'southern North Island').
Canonical variate (CV) analysis (Fig. 2c) and Mahalanobis distance probabilities (Additional file 1: Table S9) show all genera form significantly different groups, with a cross-validation accuracy of 100%. Canonical function 1 (CV1; 53.9% among-group variance) clearly distinguishes Naultinus and Woodworthia, which occupy opposite extremes of the morphospace (Fig. 2c). The positive sector of CV1 describes shortening of the nasal margin and adjacent medial flange, with corresponding shortening in the prefrontal margin (similar to PC1; Additional file 1: Figure S7). Hoplodactylus occupies the extreme positive region of canonical function 2 (CV2; 30% among-group variance), characterized by a relative slope decrease of the nasal margin and consequent shortening of the orbital margin ( Fig. 2c; Additional file 1: Figure S7).
The Holocene subfossil specimens are broadly distributed throughout the positive zones of CV1 and CV2 (Fig. 2c), with some individuals visually falling within the 95% confidence-interval of the morphospaces of extant genera (Hoplodactylus: D, E, J, K; Woodworthia: B). Typicality probabilities of Mahalanobis distances across all CVs (Table 1) show that while many Holocene subfossil specimens strongly associate with Hoplodactylus (A, D, E, F, J, K), other specimens (B, C, G, H, I) show no clear phylogenetic affinities, indicating that Holocene subfossil specimens display greater variation in maxillary form than that encompassed by the extant genera. Conversely, despite posterior probabilities (Table 1) showing similar significant support for Holocene subfossil Hoplodactylus classification (A, C, D, E, F, H, J, K), specimens with lower typicality probabilities were assigned to Woodworthia (B, G, I).

Variation and morphological convergence in diplodactylid maxillae
Phylogenetic position is a significant predictor of maxilla shape diversity in New Zealand diplodactylids, with all genera (Dactylocnemis, Hoplodactylus, Mokopirirakau, Naultinus and Woodworthia) being morphologically distinct. Identification of taxonomically informative morphological variation within a single skeletal element contrasts with previous assertions of skeletal uniformity in New Zealand's geckos (e.g. [23,26]).
Variation in diplodactylid maxilla shape is predominantly explained by two characters, described by the first two axes of both PCA and CVA: (a) posterior extension/reduction of the nasal margin; and (b) increase/ decrease in dorsoventral extent of the facial process. Separation of genera along PC1 appears correlated with broad habitat use of New Zealand diplodactylids, with terrestrial-arboreal (Dactylocnemis, Hoplodactylus and Woodworthia) and exclusively arboreal (Naultinus) genera occupying positive and negative regions respectively [45,46]. This morphological signature of habitat use extends to species-level comparison, most notably for discrimination of the terrestrial-arboreal M. 'southern North Island' from the arboreal M. granulatus [47], characterized by a shift to more positive values.
In lizards, arboreal forms tend towards broad, pointed skulls, and, similar to saxicoline species, tend to be dorsoventrally flattened, presumably enabling faster climbing speeds on non-horizontal surfaces [48][49][50]. While cranial modifications associated with habitat use are undocumented for New Zealand diplodactylids, extension of the nasal margin in arboreal species appears to be linked to two superficial morphological changes in the adjacent prefrontal margin: (a) a reduction in anterior extent (observed in other Gekkota; [51]); and (b) formation of a thickened ridge along the prefrontal orbital margin (Additional file 1: Figure S8). While the function of these features remains unclear, association with arboreality may indicate ecomorphological convergence between phylogenetically independent lineages. Despite describing similar shape change, separation of genera along CV1 reflects broad phylogenetic relationships, distinguishing broad (Hoplodactylus, Woodworthia) and narrow (Dactylocnemis, Mokopirirakau, Naultinus) toed clades by positive and negative values respectively; supporting previous morphological classification [16].
In addition to habitat use, skull-shape evolution in lizards is strongly influenced by diet, with shape variation concentrated in the premaxilla, nasal and jaw joint, reflecting their roles in jaw-based prehension and feeding biomechanics [48,52]. Herbivorous lizard skulls tend towards reduced snout lengths and high temporal regions relative to carnivorous lizards, contributing to increased bite strength required for processing fibrous and tough foliage [53][54][55]. Conversely, omnivorous gekkotans represent intermediate forms not specialized for particular feeding behaviors, and consequently lack unique morphological adaptations [56]. New Zealand geckos are predominantly omnivorous, consuming a wide variety of food items, including plant matter (fruit, honeydew and nectar) and arthropods [46]. Such extensive dietary overlap affects the use of diet as a variable of maxilla shape, given that the categories (omnivorous and insectivorous) are not discrete. Finally, as lizard maxillae are evolutionarily conserved, exhibiting reduced disparity relative to rate of evolution (compared with other cranial elements; [52]); similar analyses of elements critical to cranial biomechanics (e.g. quadrate, which also supports the auditory system) may enhance detection of stronger species-level morphological signals in the New Zealand Diplodactylidae.

Efficacy of size-based discrimination
Maxilla size is significantly correlated with phylogenetic affinity; however, only Hoplodactylus can be fully differentiated (under post-hoc comparisons), with the remaining diplodactylid genera exhibiting variable degrees of overlap. This highlights the inefficiency of previous sizebased taxonomic identification of non-Hoplodactylus Holocene subfossil geckos, especially intermediate-sized genera (Dactylocnemis, Mokopirirakau and Naultinus), which exhibit complete size overlap. Similarly, while large size proves reliable in discriminating extant H. duvaucelii, application to Holocene subfossil identification requires assumptions of temporal taxonomic homogeneity (or "covert biases"; [57]).

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. This suggests that well-characterized biodiversity reductions (and extinctions) observed across insular avifauna [64,65] extend to lineages comprised of taxa of smaller body size, including herpetofauna.
Combined Procrustes and Mahalanobis distance comparisons provide support for previous size-based assignment of five Holocene subfossils (A, D, E, J, K) to H. duvaucelii, confirming assumed prehuman distribution of this species across both the North and South Islands. The remaining six Holocene subfossil specimens (B, C, F, G, H, I) exhibit classificatory discrepancies and/or reduced assignment probabilities (below relevant thresholds), reflected by their unique position across CV1/CV2. These distinct Holocene subfossil maxillae ("unknown taxa") do not reflect differential adaptation of H. duvaucelii to mainland and island habitats (given shape overlap of mainland and island populations) but reflect either increased morphological diversity of mainland large species (not encompassed by extant populations) or the presence of at least one extinct, large, broad-toed diplodactylid species.
Based on digit morphology, the extinct giant H. delcourti was positioned within the broad-toed clade, sister to H. duvaucelii [20], suggesting these "unknown taxa" could potentially represent small or even juvenile H. delcourti. However, this seems unlikely given the paucity of reported subfossil remains of H. delcourti [66], despite extensive collections of other diplodactylid taxa [26]. More precise phylogenetic affinities 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, reflected in a reduction in maximum maxilla size (a proxy for body size; e.g. [58]). Such sized-biased extinction is well-documented for Quaternary insular lizards globally [59,[67][68][69], including the extinction of two largebodied eugongyline skink species (Oligosoma northlandi and Oligosoma sp.) in northern New Zealand [25,27,70]. This reflects the inherent vulnerability of New Zealand's large-bodied, nocturnal herpetofauna to high predation rates and ecological displacement by exotic mammals (including the Pacific rat Rattus exulans (kiore); [71,72]), particularly in forest-cleared environments [73]. Smaller lizards can escape predation during periods of inactivity through utilizing narrow retreats, given limited overlap in body diameter with small mammalian predators [45]. Conversely, refugia utilized by large-bodied lizards can be accessed by mammalian predators, as evidenced by reductions in body weight, tail width and recruitment of H. duvaucelii on kiore-inhabited islands [34,74].
Similarly to extant H. duvaucelii populations [75], Holocene subfossil H. duvaucelii also exhibit a latitudinal cline in maxilla size, in opposition to Bergmann's rule (i.e. increased size at higher 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 sunbasking [76,77]. Despite being nocturnal, Duvaucel's gecko occasionally emerges from retreats to thermoregulate through cryptic sun-basking [78,79], suggesting that small body size provides an adaptive advantage at high latitudes.

Conclusions
The majority of New Zealand diplodactylid genera can be differentiated from each other based exclusively upon the shape of the maxilla, which exhibits strong correlations with phylogenetic relationships. Additional species-level discrimination based on ecomorphological adaptations highlights the potential application of geometric morphometrics to the morphological characterization of highly functionally variable elements (or whole skulls) in taxonomic descriptions of extant diplodactylid species. Previous sized-based identification of Holocene subfossils is ineffective and underestimates extinct diversity, suggesting global assemblages of insular reptiles are depauperate in comparison to their prehuman diversity.

Geometric morphometrics
Geometric morphometric analyses were performed on a total of 94 maxillae (see Additional file 1: Methods for additional analytical details). Three-dimensional rendered surface models were generated from micro-CT reconstructions of maxillae, with shape characterized by 15 landmarks and 40 sliding semi-landmarks ( Fig. 1b; Additional file 1: Figures S2, Table S2) digitized in Checkpoint (Stratovan Corporation, Davis, CA). Landmark coordinates were aligned using a generalized least-squares Procrustes superimposition [38], with semi-landmark position optimized using the Procrustes distance criterion [83] and paired elements symmetrized (following mirroring of left maxillae coordinates; Additional file 1: Table S3).
Holocene subfossil maxillae were then projected into these two-dimensional morphospaces (i.e. PCA and CVA) through matrix multiplication with respective eigenvectors (e.g. [93]). Phylogenetic classification 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 [94,95] and posterior [96] probabilities Variation in size of the maxilla (represented as centroid-size of the landmark configuration) between genera was examined using a one-way ANOVA and Tukey's honestly significant difference (HSD) post-hoc tests [97], and visualised using a barplot. All statistical analyses were performed in the R statistical environment v. 3.6.1 [98] using the packages geomorph v. 3.1.2 [92] and Morpho v. 2.7 [99].
Additional file 2. Three-dimensional landmark coordinates for modern and subfossil maxilla.
Additional file 3. R code for geometric morphometric and statistical analyses.
Otago) for coding assistance; and both Aaron Bauer (Villanova University) and Daniel Paluh (University of Florida) for access to diplodactylid cranial scans. In addition, we would like to thank both Alexander Verry and Kerry Walton (both University of Otago) for manuscript review. We would also like to thank reviewers 1 (Anonymous) and 2 (Jendrian Riedel) for their constructive comments on our manuscript.

Authors' contributions
LS and NR conceived the study; LS carried out data collection and analyses with assistance from ES; ES, RH and NR assisted with data interpretation; LS drafted the manuscript, and all authors edited the manuscript; NR provided funding. All authors read and approved the final manuscript.

Funding
This work was supported by the Departments of Geology and Zoology, University of Otago; and a Royal Society of New Zealand Marsden FastStart Grant (16-UOO-096).

Availability of data and materials
The dataset (i.e. raw landmark coordinates and R-code) supporting the conclusions of this article is included within the article and its Additional file 2 and Additional file 3.