A statistical reanalysis of morphological differentiation among island night lizards (Xantusia riversiana) from the California Channel Islands

This study re-analyzes morphometric and meristic data among island night lizards, Xantusia riversiana, from the California Channel Islands of San Clemente, Santa Barbara, and San Nicolas in order to ascertain whether the implementation of different statistical methods can recover different results that could potentially alter biological interpretations. Our results concur with a recent previous study demonstrating that the three island populations differ morphologically from one another and that the San Nicolas Island population is the most divergent. Several important aspects, however, of the previous study depart significantly from those recovered here. Our analyses found sexual dimorphism within each population for both morphometric and meristic characters to be relatively uncommon whereas the previous study found nearly all characters to be sexually dimorphic for all island populations. The previous study also recovered significant differences among the three island populations for all morphometric characters whereas far fewer differences were recovered in the present study. Both studies found few significant inter-island differences among the meristic characters. The discordances between these two studies stem from differences in the a priori treatment of the raw character data and the different downstream statistical analyses and visualization techniques used on those data. This was particularly relevant with the use here of an allometric growth algorithm for size-correcting the morphometric data not used in the previous study and by treating all three populations as independently evolving groups. We did not conduct analyses where data from the San Clemente and Santa Barbara island populations were conflated based on their subspecific designation (X. r. reticulata) and then compared to data from the independently evolving San Nicolas Island population. This imprudent use of taxonomy violates the assumptions of statistical independence. We emphasize that explicit justification for the use of particular statistical analyses should occur in all studies—especially if the results bear on the implementation of effective and efficient resource management programs.


Introduction
The California Channel Islands comprise the western limits of the California Borderlands-one of the most tectonically active areas in North America. Distributed across these eight islands are a suite of endemic vertebrates that continue to serve as exemplary models for studying the evolutionary patterns and processes of island radiations (Johnson 1972;Eggert et al. 2004;Floyd et al. 2011;Funk et al. 2016;Pergams et al. 2000). One of these vertebrates, the island night lizard, Xantusia riversiana, occurs on three of the four southernmost islands-San Clemente, San Nicolas, and Santa Barbara (Fig. 1). These islands, unlike the four northern islands-Anacapa, San Miguel, Santa Cruz, and Santa Rosa-are separated by deep ocean basins and have never been connected to each other, making overwater dispersal the only likely mode of their colonization (Junger and Johnson 1980;Wenner and Johnson, 1980;Vedder and Howell 1980;Ward and Valensise 1996). Whereas the continental species of Xantusia are relatively small (adult SVLs < 70 mm) microhabitat specialists (e.g. Zweifel and Lowe 1966;Bezy 1967;1989;Webb 1970;Lee 1975;Grismer and Galvan 1986;Levitt et al. 2007;Noon et al. 2013;Grismer 2021), the insular X. riversiana is much larger (adult SVLs 70-117 mm) and more of a habitat generalist but with possible inter-island differences in microhabitat preference, diet, and other aspects of their ecology (Fellers et al. 1998(Fellers et al. , 2008Fellers and Drost 1991;Mautz et al. 1993;Drost et al. 2018). Many studies across multiple lineages of lizards have demonstrated that morphology reflects habitat use as well as a number of other ecological traits (e.g. Williams 1972;Melville and Swain 2000;Bergmann and Irschick 2010;Losos 2010;Lee et al. 2013;Meiri 2013, 2015;Grismer and Grismer 2017;Grismer et al. 2015Grismer et al. , 2018Grismer et al. , 2020aKaatz et al. 2021;Pelegrin et al. 2021;Pianka et al. 2017;Toyama 2017) and that "ecomorphology" may represent an average of long-term responses to multiple ecological parameters (Pianka et al. 2017). In fact, the application of ecomorphology to uncover, explain, and predict ecological patterns has become formalized into a hypothetico-deductive framework that reveals causal roles of morphological traits in organismal ecology (Ricklefs and Miles 1994;Wainwright and Reilly 1994;Wainwright 1996;Kingsolver and Huey 2003;Feilich and López-Fernández 2019;Kaatz et al. 2020).
Considerable high-quality field research on the ecology of the three populations of Xantusia riversiana has been ongoing since the 1990's (e.g., Fellers et al. 1998Fellers et al. , 2008Fellers and Drost 1991;Mautz et al. 1993;Drost et al. 2018) but it was only recently that the first focused, statistically based, comparative morphometric and meristic analyses among the three populations were conducted. Adams et al. (2018) employed a battery of poorly explained and poorly justified statistical analyses on morphometric and meristic characters using large sample sizes from all three islands to investigate the morphological differences among the island populations and between the subspecies X. r. riversiana from San Nicolas Island and X. r. reticulata from San Clemente and Santa Barbara islands (Smith 1946). Using their raw data, we employed a standardized analytical framework for morphological comparisons proposed by Chan and Grismer (2021a) that is more streamlined and considerably different from the framework employed by Adams et al. (2018) to determine whether the implementation of different statistical methods can recover different results that could potentially alter biological interpretations. Although both studies concluded that all island populations were morphologically distinct and that the San Nicolas Island population was the most divergent, there were notable differences in several important details. For consistency, we compare and contrast the results of different sets of analyses and visually interpret the efficacy of these analyses using many of the same multivariate visualization techniques employed by Adams et al. (2018).

Samples and data
All morphological analyses were based on the raw data set of Adams et al. (2018) which included 172 adults: 37 females and 32 males from San Clemente Island, 28 females and 12 males from Santa Barbara Island, and 32 females and 31 males from San Nicolas Island. The data set included nine mensural (i.e. morphometric) characters: snoutvent length (SVL), head length (HL), head width (HW), head depth (HD), snout length (SNT), interorbital distance (IO), forelimb length (FLL), hind limb length (HLL) and pelvic width (PW). Tail length was omitted from our analyses owing to a high number of samples with regenerated or broken tails. Five meristic characters examined included ventral scales (VS), gular scales (GS), preanal (i.e. precloacal used here; reptiles do not have an anus) scales (PA), femoral pores (FP), and fourth toe subdigital lamellae (TL4). Methods for measuring body dimensions and counting scales are detailed in Adams et al. (2018).

Statistical analyses
We followed a standardized statistical protocol for morphological analyses (Chan and Grismer 2021a) and employed multivariate analyses and data visualization methods recently used for other species of Xantusia (Grismer 2021) and commonly used in lizard taxonomy. Given that the three island populations are allopatric and there is no evidence of gene flow among them (Bezy et al. 1980;Noonan et al. 2013;Rice 2017), they were treated as independent evolutionary units regardless of their taxonomy (Smith 1946) so as to avoid violating the assumptions of statistical independence.
Because preliminary Student t-tests on males and females from each island detected sexual dimorphism among some of the adjusted morphometric (see below) and raw meristic characters within each population, males and females were analyzed separately in all subsequent analyses unless stated otherwise. Adams et al. (2018) did not log-transform their meristic data but without explanation size-corrected them by dividing each character by its geometric mean. We conducted no a priori treatment of the raw meristic data because scale patterns form early in embryogenesis (Alibardi 1996) and do not change during ontogeny (Chang et al. 2009). However, to ensure that allometric biases (variation due to differences in body size) in the raw morphometric data were appropriately removed prior to analysis, these data were adjusted using the allometric formula: X adj =log(X)-β[log (SVL)log(SVL mean )], where X adj =adjusted value; X=measured value; β=unstandardized regression coefficient for each population; and SVL mean =overall average SVL of all populations (Thorpe 1975(Thorpe , 1983Turan 1999;Lleonart et al. 2000)-implemented through the R package GroupStruct ( available at https://github.com/chankinonn/ GroupStruct). Mor phometric adjustments were conducted separately on each sex from each island. The data were then concatenated by sex in order to construct separate data frames composed of only males and only females from each island. This ensured there was no intersexual or interisland conflation of variation (Reist 1985;McCoy et al. 2006). To remove the effects of allometry, Adams et al. (2018) log-transformed all morphometric data and divided each measurement (except for SVL) by its geometric mean.
One-way analyses of variance (ANOVA) were performed on the morphometric and meristic data sets to test for the presence of statistically significant interisland differences. If detected, Tukey HSD post hoc tests were performed to ascertain which pairs of island populations had significantly different mean values for which character(s). All statistical tests were performed at a significance threshold of α = 0.05. Boxplots for the discrete meristic characters and violin plots embedded with boxplots for continuous morphometric characters were generated to visualize the distribution of variation among the characters across islands. Summary statistics (mean, range, and ±1 standard deviation) were generated for all meristic and adjusted morphometric characters.
Separate and concatenated principal component analyses (PCA) of the morphometric and meristic data were employed to visualize and assess the degree of morphospatial separation among the three populations. PCA is an unsupervised analysis that does not group individuals a priori according to species/population/sex. A non-parametric permutation multivariate analysis of variance (PERMANOVA) from the vegan package in R (Oksanan et al. 2020) was used to determine if the centroid locations of each population in PCA were statistically different (Skalaski et al. 2018). The analysis is based on the prior calculation of the distance between any two data points in a Euclidean (dis)similarity matrix using 5000 permutations, not on the output of the PCA. A pairwise post hoc test calculates the differences between all combinations of population pairs, generating a Bonferroni-adjusted p-value and a pseudo-F ratio (F statistic). p < 0.05 is considered significant and larger F statistic values indicate more pronounced group separation. A rejection of the null hypothesis (i.e. centroid positions and/or the spread of the data points are no different from random) signifies a difference between populations.
A subsequent supervised analysis, discriminant analysis of principal components (DAPC; Jombart and Collins 2015) that does group individuals a priori, was also performed. DAPC relies on data calculated from a PCA as a prior step to ensure that the variables analyzed are not correlated and number fewer than the sample size. Dimension reduction of the DAPC prior to plotting is accomplished by retaining the first set of principal components that account for 90-95% of the variation (Jombart and Collins 2015) as determined from a scree plot generated as part of the analysis. The DAPC analyses were run using the adegenet package in R. All analyses were performed in R (R Core Team 2014).

Results and comparisons with Adams et al. (2018)
Comparing the results recovered here with those of Adams et al. (2018) was, in many cases, hampered by the fact that Adams et al. (2018) often reported significant differences resulting from their ANOVAs on the size-corrected SVLs and analyses of covariance (ANCOVAs) on all other size-corrected characters using size-corrected SVL as a covariate. However, they presented no evidence of any post hoc tests so it remains unclear how their comparisons were obtained.

Sexual dimorphism
Overview. Overall, the following analyses found far less sexual dimorphism within all three island populations compared to that of Adams et al. (2018) as summarized in Table 1. The differences between these studies stem from their ambiguous justification for using MANOVAs, ANOVAs, and ANCOVAs on incorrectly size-adjusted morphometric characters and inappropriately size-corrected meristic characters (see below). Whereas Adams et al. (2018) found all morphometric characters to be sexually dimorphic in all populations, only snout-vent length (SVL), head depth (HD), and pelvic width (PW) in the San Clemente Island population and PW in the San-ta Barbara Island population were recovered as sexually dimorphic. No sexual dimorphism was recovered in the San Nicolas Island population. Both Adams et al. (2018) and the analyses here recovered sexual dimorphism in the meristic characters of precloacal scales (PA) and femoral pores (FP) in the San Clemente Island population and gular scales (GS) and PA in the San Nicolas Island population. The analysis here recovered no sexual dimorphism in the Santa Barbara Island population whereas Adam's et al. (2018) found this population to be sexually dimorphic for GS, PA, and FP.
Morphometric characters. All characters conformed to parametric test assumptions of normality (Shapiro-Wilk test, p < 0.05) and homogeneity of variances (F-test, p > 0.05). Student t-tests between the sexes for each character for each island population recovered minimal instances of sexual dimorphism (Table 1) while two-way ANOVAs and Tukey post hoc tests using island and sex as independent variables recovered no sexual dimorphism. Student t-tests indicated that males from San Clemente Island were significantly smaller (SVL), had flatter heads (HD), and narrower pelves (PW) than females (Tables 1 and 2). Males from Santa Barbara Island also had significantly narrower pelves than females. No sexual dimorphism was recovered in the San Nicolas Island population. The PCA and DAPC analyses concur, showing wide overlap and slight separation between the sexes from each island (Fig. 2). For the Santa Barbara Island population, principal component (PC1) accounted for 31.8% of the variation and loaded most heavily for head depth (HD), head length (HL), forelimb length (FLL), and pelvic width (PW) ( Table 3). The second principal component (PC2) accounted for 15.9% of the variation and loaded most heavily for head width (HW). For the San Clemente Island population PC1 accounted for 28.4% of the variation and loaded most heavily for head length (HL) and head width (HW) ( Table 3). PC2 accounted for 17.9% of the variation and loaded most heavily for hind limb length (HLL). For the San Nicolas Island population PC1 accounted for 26.8% of the variation and loaded most heavily for snout length (SNT) and PC2 accounted for 19.4% of the variation and loaded most heavily for forelimb length (FLL) ( Table 3). The PERMANOVA analysis recovered no statically significant differences (p adjusted values for Santa Barbara, San Clemente, and San Nicolas islands = 0.42, 0.10, and 0.21, respectively) in centroid placement in the PCA (Table 4).
Using an ANOVA on the size-corrected SVLs and analyses of covariance (ANCOVAs) on all other size-corrected characters with size-corrected SVL as a covariate, Adams et al. (2018) stated that all morphometric characters (see below), except tail length, to be sexually dimorphic within all populations (i.e. their "Xantusia riversiana") although PW was not significantly different within the combined data of the San Clemente and Santa Barbara island populations (i.e. X. r. recticulata) nor in the San Nicolas Island population. However, no post hoc tests were discussed so it is not known how these comparisons were obtained. Table 1. p-values for statistically significant mean differences in the characters among males and females from each island population based on Student t-tests. Character abbreviations are in the Materials and methods. Orange shaded cells denote sexually dimorphic characters reported by Adams et al. (2018), green shaded cells denote sexually dimorphic characters recovered herein and by Adams et al. (2018), and gray shaded cells denote characters for which neither study found sexual dimorphism.  Meristic characters. All characters conformed to parametric test assumptions of normality (Shapiro-Wilk test, p < 0.05) and homogeneity of variances (F-test, p > 0.05). Student t-tests indicated that males from the San Clemente Island population have significantly more precloacal scales and femoral pores (PA and FP, respectively) than females (Tables 1 and 5). No sexually dimorphic meristic characters were recovered in the Santa Barbara Island population and the males of the San Nicolas Island population had significantly more gular (GS) and PA scales (Tables 1 and 5). The PCA and DAPC analyses mirrored those of the morphometric characters in showing wide overlap between the sexes from each island (Fig. 3). In the Santa Barbara Island population, PC1 accounted for 36.5% of the variation and loaded most heavily for ventral scales (VS) ( Table 6). PC2 accounted for 22.7% of the variation and loaded most heavily gular scales (GS) and fourth toe lamellae (TL4). In the San Clemente Island population, PC1 accounted for 30.3% of the variation and loaded most heavily for ventral scales (VS) and gular scales (GS) ( Table 6). PC2 accounted for 24.7% of the variation and loaded most heavily preanal scales (PA). In the San Nicolas Island population, PC1 accounted for 36.9% of the variation and loaded most heavily for ventral scales (VS) and gular scales (GS) and PC2 accounted for 19.1% of the variation and loaded most heavily for precloacal scales (PA) and femoral pores (FP) ( Table 6). The PERMANOVA analysis recovered statically significant differences in centroid placement in the PCAs for the San Clemente and San Nicolas island populations (p adjusted = 0.0002 and 0.04, respectively) but not for the Santa Barbara population (p adjusted = 0.55) ( Table 4).
Using ANCOVA analyses with size-corrected SVL as a covariate, Adams et al. (2018) found all size-corrected meristic characters except femoral pores, to be sexually dimorphic between the two subspecies but in an apparent contradiction stated there were no "dramatic sexual dimorphisms" in the size-corrected meristic characters except for males having more "preanal" scales (PA) than females in the San Nicolas and San Clemente Island populations. However, their Table 1 of the raw meristic data shows all means for each character to be nearly equal with widely overlapping ranges.

Inter-island differences
Overview. The results below recover statistical morphometric and meristic differences among all three island populations using the data types separately or combined as summarized in Table 7. The San Nicolas Island population is considerably different from the San Clemente and Santa Barbara Island populations which are morphologically similar but also statistically different for many characters. Adams et al. (2018) found similar results but across a different suite of characters based on different and sometimes inappropriately used analyses (see above).
Female morphometrics. ANOVA and Tukey HSD post hoc analyses recovered a number of statistical differences between San Nicolas Island and the other two islands (Tables 2 and 7). San Clemente and Santa Barbara islands differed only in head width (HW). The PCA recovered modest separation of the San Nicolas Island population and broad overlap of the San Clemente and Santa Barbara island populations (Fig. 4A). PC1 accounted for 53.6% of the variation in the data set and loaded most heavily for snout length (SNT), forelimb length (FLL), interorbital distance (IO), and less so for head length, width, and depth (HL, HW, and HD, respectively) (Fig. 5). PC2 accounted for 9.8% of the variation in the data set and loaded most heavily for snout-vent length (SVL). The DAPC recovered complete separation of the San Nicolas Island population and modest overlap of the 66% confidence ellipses between the San Clemente and Santa Barbara island populations (Fig. 4B). Quantitative differences among the characters among the three populations are visualized in Figure 4C illustrating the higher mean value of all characters in the San Nicolas Island population. The PERMANOVA analysis recovered statistically significant differences in centroid placement among all island pairs with those comparing San Nicolas Island being magnitudes greater in difference compared to that between San Clemente and Santa Barbara islands (F statistic 58.34 and 60.42 versus 4.25; Table 4).
Using an ANOVA on the size-corrected SVLs and AN-COVAs on all other size-corrected morphometric characters with size-corrected SVL as a covariate, Adams et al. (2018: Table 2) recovered significant differences among the three island populations for all characters for both sexes (Table 8). However, no post hoc tests were indicated to differentiate which island pairs differed for which characters. Adams et al. (2018) stated that "X. r. reticu-lata have relatively longer heads". Although their data were not shown, our data indicate that lizards from the San Nicolas Island population have the longest, widest, and thickest heads. Additionally, their PCA analyses recovered no separation among any island populations for either sex.
Female meristics. ANOVA and Tukey HSD post hoc analyses recovered a number significant differences among the females from all islands for ventral and gular scales (VS and GS, respectively) and subdigital lamellae (TL4) (Tables 5 and 7). The ANOVAs and Tukey HSDs post hoc tests also recovered significant differences among the Santa Barbara and San Nicolas populations for precloacal scales (PA). The PCA recovered much less separation among the three populations compared to that recovered by the morphometric data (Fig. 6A). PC1 accounted for 41.4% of the variation and loaded most heavily for gular scales (GS), fourth toe lamellae (TL4) and ventral scales (VS) (Fig. 7). PC2 accounted for 22.0% of the variation and loaded most heavily for femoral pores (FP). The DAPC mirrored that of the female morphometrics in recovering complete separation of the San Nicolas Island population and notable overlap of the 66% confidence ellipses between the San Clemente and Santa Barbara island populations (Fig. 6B). Quantitative differences in the characters among the three populations are illustrated in Figure 6C. The PERMANOVA analysis recovered statistically significant differences in centroid placement among all island pairs for females with those compared to San Nicolas Island being magnitudes greater in difference compared to that between San Clemente and Santa Barbara islands (F statistic 17.89-28.63 versus 7.98; Table 4).
Conducting ANCOVAs on the size-corrected meristic data set using the size-corrected SVL as a covariate, Adams et al. (2018: Table 2) recovered significant differences among the three island populations for ventral scales (VS), gular scales (GS), precloacal scales (PA), and fourth toe lamellae (TL4) but only among GS and PA between the sexes. However, no post hoc tests were indicated to differentiate which island pairs differed for which character. Their combined data set from the San Clemente and Santa Barbara islands (i.e. Xantusia river- siana reticulata) recovered significant differences in all meristic characters between it and the San Nicolas Island population (Adams et al. 2018).
Female concatenated data. The concatenated morphometric and meristic data set recovered the same pattern as the individual data sets of each data type but with even clearer separation of the San Nicolas Island population from the other island populations in the PCA (Fig. 8A). PC1 accounted for 41.4% of the variation and loaded most heavily for the morphometric characters of snout length (SNT), head depth (HD), head width (HW), head length (HL), interorbital distance (IO), forelimb length (FLL), and pelvic width (Fig. 9). PC2 accounted for 10.8% of the variation and loaded most heavily for fourth toe lamellae (TL4) and ventral scales (VS). These analyses and the concatenated analysis in particular, indicate that overall morphological variation in females is disproportionately driven by body shape (compare Figs  4A, 6A, and 8A). The PERMANOVA analysis recovered statistically significant differences in centroid placement among all island pairs with those among San Nicolas Island being magnitudes greater in difference compared to that between San Clemente and Santa Barbara islands (F statistic 41.91 and 47.10 versus 5.50; Table 4).
Male morphometrics. The pattern of morphometric variation in the males mirrored that of the females. ANOVA and Tukey HSD post hoc analyses recovered a number of significant differences between San Nicolas Island and the other two islands for all characters (Tables 2 and 7).
San Clemente and Santa Barbara islands differed only in head width (HW). The PCA recovered nearly complete separation of the San Nicolas Island population and broad overlap of the others (Fig. 10A). PC1 accounted for 59.9% of the variation in the data set and loaded most heavily for snout length (SNT), interorbital distance (IO), head length (HL), and head depth (HD) (Fig. 5). PC2 accounted for 9.6% of the variation and loaded most heavily for snout-vent length (SVL). The DAPC recovered nearly complete separation of all three populations with minimal overlap of the 66% confidence ellipses between the  San Clemente and Santa Barbara island populations (Fig.  10B). Quantitative differences in the characters among the three populations are illustrated in Figure 10C, again showing the higher mean values across all characters in the San Nicolas Island population. The PERMANOVA analysis recovered statistically significant differences in centroid placement among all island pairs with those compared to San Nicolas Island being magnitudes greater in difference compared to that between San Clemente and Santa Barbara islands (F statistic 46.06 and 66.56 versus 2.96; Table 4).
Using an ANOVA on the size-corrected SVLs and AN-COVAs on all other size-corrected morphometric characters with size-corrected SVL as a covariate, Adams et al. (2018: Table 2) recovered significant differences among the three island populations for all characters. However, no post hoc tests were indicated to differentiate which island pairs differed for which character. Additionally, their PCA analyses recovered no separation among any island populations for either sex.
Male meristics. The pattern of meristic variation in the males also mirrored that of the females. ANOVA and Tukey HSD post hoc analyses recovered a number of statistical differences among all islands for gular scales (GA), between the San Clemente and San Nicolas island populations and between the Santa Barbara and San Nicolas island populations for subdigital lamellae (TL4), between the San Clemente and San Nicolas island populations for precloacal scales (PA), and between the San Clemente and San Nicolas island populations for ventral scales (VS) (Tables 5 and 7). The PCA recovered much less separation among the three populations than did the PCA of morphometric data (Fig. 11A). PC1 accounted for 39.5% of the variation and loaded most heavily for gular scales (GS), fourth toe lamellae (TL4), and ventral scales (VS) (Fig. 7). PC2 accounted for 20.0% of the variation and loaded most heavily for precloacal scales (PA). However, the DAPC recovered complete separation of the San Nicolas Island population with notable overlap between the San Clemente and Santa Barbara island populations (Fig. 11B). Quantitative differences in the characters among the three populations are illustrated in Figure 11C. The PERMANOVA analysis recovered statistically significant differences between San Nicolas Island and the other islands (p adjusted = 6.00E -5 ) but not between San Clemente and Santa Barbara islands for males (p adjusted = 0.062) ( Table 4).
Conducting ANCOVAs on the size-corrected meristic data set using the size-corrected SVL as a covariate, Adams et al. (2018: Table 2) recovered significant differences among the three island populations for ventral scales (VS), gular scales (GS), precloacal scales (PA), and fourth toe lamellae (TL4) but only among gular scales and precloacal scales between the sexes. However, no post hoc tests were indicated to differentiate which island pairs differed for which character.
Male concatenated data. The PCA of the concatenated morphometric and meristic data from all islands was consistent with the individual PCAs and that of the female concatenated PCA in that it recovered complete separation of the San Nicolas Island population from the other island populations which overlapped greatly (Fig. 8C). PC1 accounted for 43.8% of the variation and loaded most heavily for snout length (SNT), head depth (HD), Table 7. p-values for statistically significant mean values of the adjusted morphometric and raw mensural data between pairs of island populations determined from Tukey HSD post hoc tests following ANOVAs. Island abbreviations are SB = Santa Barbara Island, SC = San Clemente Island, and SN = San Nicolas Island. Character abbreviations are in the Materials and methods. Orange cells denote significant differences reported by Adams et al. (2018), green cells denote significant differences recovered herein and by Adams et al. (2018), and gray cells denote characters for which neither study found significant differences. head width (HW), head length (HL), interorbital distance (IO), forelimb length (FLL), hind limb length (HLL), and pelvic width (PW) (Fig. 12). PC2 accounted for 10.0% of the variation and loaded most heavily for fourth toe lamellae (TL4). The separation was less clear among all three populations in the DAPC with wide overlap between the San Clemente and Santa Barbara island populations (Fig.  8D). As with females, this analysis indicates that overall variation in the males is disproportionately driven by body shape (Figs 4A, 6A, and 8A). The PERMANOVA analysis recovered statistically significant differences in centroid placement among all island pairs with those compared to San Nicolas Island being magnitudes greater in difference compared to that between San Clemente and Santa Barbara islands (F statistic 31.25 and 40.68 versus 3.20; Table 4).
Concatenation of all data and both sexes. The PCA of the concatenated data set that included all characters, all islands, and males and females, recovered the same clear separation of the San Nicolas Island population from the other two island populations which again, showed complete overlap (Fig. 13A). These data also continue to underscore the weak separation between males and females from their respective islands. PC1 accounted for 41.6% of the variation and loaded most heavily for the morphometric characters snout length (SNT), head depth (HD), head width (HW), head length (HL), interorbital distance (IO), forelimb length (FLL), hind limb length (HLL), and pelvic width (PW) (Fig. 14). PC2 accounted for 10.4% of the variation and loaded most heavily for fourth toe lamellae (TL4), ventral scales (VS), and preanal scales (PA). The DAPC recovered complete separation of the 66% confidence ellipses of San Nicolas Island males from both sexes of the San Clemente and Santa Barbara island populations. The San Nicolas Island females slightly overlapped with both sexes from San Clemente Island. Both sexes from San Clemente and Santa Barbara islands broadly overlapped (Fig. 13B). As with the other concatenated analyses, this analysis also indicates that the overall variation among the island populations is disproportionately driven by body shape. These results are consistent with those of Grismer et al. (2018Grismer et al. ( , 2020aGrismer et al. ( , 2021 indicating that concatenated data sets generally outperform other data sets in PCA and DAPC analyses by showing greater separation among the operational taxonomic units (OTUs). The PERMANOVA analysis recovered statistically significant differences in centroid placement among all possible combinations of island pairs and sexes except for those between Santa Barbara Island males and females and San Nicolas Island males and females. As before, F statistic values were magnitudes greater in difference when comparing any population to the San Nicolas Island population (Table 4). Table 8. Comparisons of a priori treatment of raw morphological data and subsequent analyses and visualization techniques between Adams et al. (2018) and those herein. Only inter-island analyses are listed. Subspecies analyses are not included.

Adams et al. (2018) This study
a priori treatment of raw data a priori treatment of raw data

log-transformation of morphometric data
Size correction of log-transformed morphometric data using geometric means size-correction of raw morphometric data using the GroupStruct equation Size correction of meristic data presumably using geometric means raw meristic data untreated

Sexual dimorphism Sexual dimorphism
ANOVA on log-transformed SVLs Student t-test for each island population ANCOVA on log-transformed morphometric data and meristic data with island and sex as factors and size-corrected log-transformed SVL as a covariate.
two-way ANOVA using species and sex as independent variables MANOVA including all islands Discriminant analyses. We did not perform a linear discriminant analysis (LDA) on any data set as there were no questions as to the insular provenance of any sample. Furthermore, LDA does not always return the correct provenance probabilities (see Grismer et al. [2019] and the data for Santa Barbara Island of Adams et al. [2018]). As with the PCA plots of Adams et al. (2018), their LDA plots showed no separation among any of the populations. We did perform discriminate analyses (DAPC) on the PC scores of the most influential PCs which generally showed modest to complete separation of the San Nicolas Island population and overlap between the other two populations across the varying data sets (see above).

Comparison of the statistical treatment of data
Differences between the results of Adams et al. (2018) and those recovered here stem from differences in the a priori treatment of the raw data (i.e. body-size correction and data transformation), the different statistical analyses and visualization techniques used on those data, their lack of post hoc tests to explain putative significant differences, and our decision not to violate the assumption of statistical independence by conflating data from the San Clemente and Santa Barbara island populations based on their taxonomy given there is no evidence of gene flow between them (Bezy et al. 1980;Noonan et al. 2013;Rice 2017).
A priori treatment of morphometric data. Adams et al. (2018) created male and female data sets for the morphometric data. One data set contained the log-transformed values of all body measurements and the other contained log-transformed vales of body measurements that were then size-corrected "by dividing each measurement by its (i.e. the character, [our parentheses]) geometric mean across specimens." However, it is not clear if the character's geometric means were calculated separately for each island population and then concatenated into a single data frame, or if the geometric means for each character were calculated from the combined measurements from each population. If the latter, inter-populational differences in body shape-which we know exists-would be conflated and potentially affect the outcome of the analyses. Additionally, a geometric mean is calculated from measurements of multiple characters (not just one) and then used to scale parameters of other characters within localized anatomical systems (e.g. different components of the cranium or pelvis) (Jungers et al. 1995;Tsegai et al. 2013). We size-corrected the morphometric data using a method that is theoretically derived from the equations of allometric growth (Lleonart et al. 2000) that has been Reist 1985)-neither study evaluated the geometric mean method. As a crosscheck, we ran the size-corrected data of Adams et al. (2018) using geometric means calculated separately from each population and then concatenated into a single data frame as well as geometric means calculated from pooled samples from all populations. In both cases, the PCAs (not shown) for both males and females showed no separation among any of the populations and were consistent with those PCAs generated by Adams et al. (2018: Fig. 2). This indicates that the allometric equation may also outperform the size-correction method using geometric means.
A priori treatment of meristic data. In their Table 2, Adams et al. (2018) state that their meristic data were size-corrected. We assume this meant the characters were divided by their geometric means as were the morphometric characters but this was not explained. As above, we do not know if the characters were size-corrected separately for each population and concatenated or if the populations were combined prior to calculating the geometric means. They did not generate PCA plots for the meristic data. We did not size-correct the meristic data nor do we understand the justification for doing so as scale counts generally do not scale isometrically with body size (i.e. the number of scales a lizard is born with does not change as it grows). We did generate PCA plots and boxplots for the meristic data which generally mirrored those plots of the morphometric data in showing the overall separation of the San Nicolas Island population and broad overlap of the San Clemente and Santa Barbara Island populations (Figs 4,6,8,10,11,13).
Differences in morphometric analyses. Adams et al. (2018) stated that the "MANOVA on the complete data set recovered significant differences among islands (P<2.2 X 10 -16 ) and between subspecies (P=3.19 X 10 -

16
; Table 3)." We are uncertain what the "complete data set" was as they went on to state they conducted separate MANOVAs on males and females for morphometric and meristic characters, comparing all islands and subspecies. By our calculations this would potentially total to 30 data sets-two data types for three islands and two sexes = 12, plus two sexes for two data types for two subspecies = 6, and if they used both log-transformed morphometric data and size-corrected log-transformed morphometric data for each sex (and additional 12 data sets), that would total 30. We conducted ANOVAs and Tukey HSD post hoc tests separately for each sex across all island populations. We find these analyses to be justified as univariate ANOVA is considered a post hoc test for MANOVA to determine which dependent variables differ significantly.  06512 -0.12049 -0.19088 -0.46185 -0.33753 -0.35656 -0.25408  Differences in meristic analyses. In some of their analyses, Adams et al. (2018) also conflated data from the San Clemente and Santa Barbara island populations based on subspecies taxonomy and employed ANCOVAs for size-corrected meristic characters with log-transformed SVL as a covariate. We conducted ANOVAs and Tukeys HSD post hoc tests on the raw meristic data separated by sex for all three populations only. We did not conflate data from separate populations.
Sexual dimorphism. Adams et al. (2018) used ANCO-VAs with island and sex as factors and log-transformed size-corrected SVL as a covariate. We used Student t-tests and two-way ANOVAs and Tukey post hoc tests using species and sex as independent variables.
Why statistics matter. Morphology is inextricably linked with many essential organismal functions, such as feeding and locomotion, which in turn influence ecological  interactions, resource use, and ultimately, survival and reproduction (Bock and von Wahlert 1965;Arnold 1983;Garland and Losos 1994). If effective and efficient resource management programs are to be implemented, the morphological data and subsequent ecological interpretations on which these programs are based, must be generated appropriately. Different interpretations from different statistical analyses often generate conflicting results. Therefore, explicit explanations of what analyses are used, what they mean, and why they are used is paramount to any thorough investigation. In some cases, a judicious use of taxonomy must also be taken into consideration. Treating separate populations as a single evolving entity (fide Adams et al. 2018) when in fact they are evolving independently, not only violates the assumptions of statistical independence, but can lead to spurious biological interpretations.

Ecomorphological interpretations
Our data indicate that the majority of overall morphological variation among the three island populations is reflected in body shape (Figs 4,6,8,10,11,and 13). Among all the data sets, lizards from San Nicolas Island are consistently differentiated from those of the other two islands and among the morphometric data sets, they consistently have the largest body metrics in being longer (SVL) and having relatively longer heads and snouts (HL and SNT, respectively); wider heads (HW) with larger interorbital distances (IO); thicker heads (HD); longer forelimbs and hind limbs (FLL and HLL, respectively), and wider pelves (PW) than do lizards from the other two populations (Tables 2 and 6; Figs 4 and 10) which overlap greatly in most characters.   All data and sexes concatenated A B We presume that based on their log-transformed data corrected for size, Adams et al. (2018) hypothesized that the longer heads in lizards from Santa Barbara and San Clemente islands may be an adaptation to having a higher proportion of vegetation in their diet. Their raw data, however, show that the mean head lengths (and all other head proportions) of lizards from those populations are shorter than those of lizards from San Nicolas Island not longer, which is consistent with our results. They did not list the size-corrected data for each island. Additionally, robust head dimensions have been reported to be an adaptation for eating a larger proportion of hard-bodied prey not softer prey (e.g. Stayton 2005;Kohlsdorf et al. 2008;Sagonas et al 2014). Bite force, in particular, has proven to be an informative performance metric in numerous taxa, which is often correlated with head morphology (principally skull width and depth), as well as diet (e.g., Binder and Van Valkenburg 2000;Herrel et al. 2001Herrel et al. , 2005Erickson et al. 2004;Dumont et al. 2005, Anderson et al. 2008). This is an example where the results from different analyses can lead to different biological interpretations. Furthermore, Charles Drost (pers com. in lit. 2021) stated "these lizards are decidedly generalized in their feeding, and the studies that have been reported from the different islands differ substantially in terms of methods. And further, the limited information from San Clemente and San Nicolas are almost certainly strongly biased in terms of both time of year and location where the specimens were collected-both of which would strongly affect diet contents in a limited sample. In short, I don't think there is any basis for serious considerations of diet differences among the islands." This statement indicates we are far from understanding the dietary nature of these populations and that diet analyses throughout the year from various parts of the islands may shed light on reasons for the significantly larger head dimensions of lizards from the San Nicolas Island. Adams et al. (2018) noted that "San Nicolas Island largely consists of sandstone, whereas Santa Barbara and San Clemente islands are composed mainly of volcanic rocks. These geological differences likely result in numerous habitat differences that could relate to the observed morphological divergence of X. r. riversiana and X. r. reticulata" although the particular differences to which they were referring were not stated. Charles Drost (pers com. in lit. 2021) however, stated that this "has no practical effect, as these are not rock-dwelling lizards; the great majority of Island night lizards on all three islands never see rocks larger than pebbles. As the studies of Island night lizard ecology have documented, the primary  habitats on all of the islands are cactus and boxthorn, on silt/clay substrates. Comparing among the islands, San Clemente stands out as most different, because of the significantly more arid conditions on that island; Santa Barbara and San Nicolas are more similar in all respects that I can see (but that does not fit the analyses presented by Adams et al., of course)." Our data demonstrated that lizards from the San Nicolas island population had relatively longer limbs than those from the other two islands. Many studies have shown that shorter limbs are more advantageous for locomoting through restrictive types of vegetation (e.g. Garland and Losos 1994;Van Damme et al. 1998;Bonine and Garland 1999;Melville and Swain 2000;Kohlsdorf et al. 2001;Irschick et al. 2005;Goodman et al. 2008;Siler and Brown, 2011;Grismer et al. 2018) and that lizards with longer hind limbs spend more time in the open away from cover (e.g. Pianka 1969Pianka , 1986. However, this has not been quantified for the San Nicolas Island population. These and other functional ecomorphological scenarios will be tested using a statistically robust quantifiable framework (i.e. Feilich and López-Fernández 2020).