Comparing the dietary niche overlap and ecomorphological differences between invasive Hemidactylus mabouia geckos and a native gecko competitor

Abstract Hemidactylus mabouia is one of the most successful, widespread invasive reptile species and has become ubiquitous across tropical urban settings in the Western Hemisphere. Its ability to thrive in close proximity to humans has been linked to the rapid disappearance of native geckos. However, aspects of Hemidactylus mabouia natural history and ecomorphology, often assumed to be linked with this effect on native populations, remain understudied or untested. Here, we combine data from ∂15N and ∂13C stable isotopes, stomach contents, and morphometric analyses of traits associated with feeding and locomotion to test alternate hypotheses of displacement between H. mabouia and a native gecko, Phyllodactylus martini, on the island of Curaçao. We demonstrate substantial overlap of invertebrate prey resources between the species, with H. mabouia stomachs containing larger arthropod prey as well as vertebrate prey. We additionally show that H. mabouia possesses several morphological advantages, including larger sizes in feeding‐associated traits and limb proportions that could offer a propulsive locomotor advantage on vertical surfaces. Together, these findings provide the first support for the hypotheses that invasive H. mabouia and native P. martini overlap in prey resources and that H. mabouia possess ecomorphological advantages over P. martini. This work provides critical context for follow‐up studies of H. mabouia and P. martini natural history and direct behavioral experiments that may ultimately illuminate the mechanisms underlying displacement on this island and act as a potential model for other systems with Hemidactylus mabouia invasions.


| INTRODUC TI ON
Hemidactylus mabouia (Tropical House Gecko) is perhaps the most pervasive and formidable gecko to invade the Western Hemisphere (Agarwal et al., 2021;Weterings & Vetter, 2018). This species will readily capitalize on the aggregation of insects around human light sources (Hughes et al., 2015), a foraging strategy which promotes high densities of individual H. mabouia that use aggressive tactics to restrict access to these spatially clustered food resources (Short & Petren, 2012;van Buurt, 2004;Williams et al., 2016). While this suggests competition to be an important aspect of H. mabouia invasions, dietary niche overlap with native species is often not known, and other hypotheses concerning possible ecological advantages over native species remain little explored. For example, the feeding mode of H. mabouia combines ambush tactics (Vitt, 1983) with active pursuit of nearby prey (Dornburg et al., 2016). Such a foraging mode could have been selective for larger hind limb or shorter fore-limb proportions that, respectively, offer a locomotor advantage over native species when accelerating or decelerating on sheer vertical surfaces common in urbanizing landscapes (Zaaf & Van Damme, 2001).
Given that both dietary overlap and morphological advantages have been invoked as major drivers of displacement in the wake of invasions by the distantly related Hemidactylus frenatus in the Pacific (Bolger & Case, 1992;Case et al., 1994;Petren & Case, 1996, 1998Short & Petren, 2012), case studies from across the range of H. mabouia are warranted if we are to gain a broader understanding of the factors that promote the success of H. mabouia invasions.
The gradual decline of the native Phyllodactylus martini (Dutch Leaf Tailed Gecko) from urbanizing areas on the Lesser Antillean island of Curaçao provides an exceptional opportunity for developing further hypotheses of what allows H. mabouia to replace ecologically similar geckos. Similar to H. mabouia, P. martini readily colonizes walls and takes advantage of insect and arachnid prey drawn to artificial lights (Dornburg et al., 2016;van Buurt, 2004). Although P. martini is predicted to fare well in suburban conditions (Dornburg et al., 2016;van Buurt, 2004), established populations have been rapidly declining following the introduction of H. mabouia in the 1980s (Dornburg et al., 2016;Hughes et al., 2015;van Buurt, 2004van Buurt, , 2006. This decline has been attributed to H. mabouia possessing a superior competitive ability based primarily on behavioral observations (Hughes et al., 2015), with no quantification of dietary overlap between these species. Therefore, it is unclear if there are differences in morphological traits relevant to foraging between H. mabouia and P. martini.
Variation in locomotor morphology and size plays vital roles in the feeding ecology of lizards. However, a quantitative comparison of these traits between H. mabouia and P. martini have yet to be conducted. Analyses of limb lengths relevant to acceleration or deceleration on vertical surfaces would provide evidence for whether H. mabouia possess locomotor traits that enable more efficient prey pursuit and capture. Likewise, head size represents another key ecological trait for investigating the invasion biology of these species.
Increased head sizes are correlated with more efficient prey capture in lizards (Dufour et al., 2018;Verwaijen et al., 2002). On the one hand, if aspects of H. mabouia cranial morphology are larger than those in P. martini, this could confer an advantage to energy-efficient prey capture for H. mabouia. On the other hand, changes in head size that are correlated with shifts in prey would suggest a segregation of dietary niches, illuminating a potential mechanism of coexistence.
Habitat type and species occupancy vary across sampling locations.
For example, both species co-occur in Lagun and Westpunt. At these sites, we restricted our sampling to suburban areas near natural habitats to maximize the potential of both species co-occurring as P. martini has been found to be absent far from edge habitats in the presence of H. mabouia (Hughes et al., 2015). In contrast, Shete Boca is a natural area in which H. mabouia are absent, while Saint Anna Bay and Willemstadt are urban areas in which P. martini are absent.
Prior to preservation, muscle biopsies were taken from each individual and dehydrated for analysis of stable isotopes. Additionally, leaf samples from each locality were collected and dehydrated for use as baselines in isotopic analyses. Specimens were then fixed in 10% formalin and later transferred to 70% ethanol and deposited in the Yale Peabody Museum of Natural History (Table S1).
Adult H. mabouia (n = 30) and P. martini (n = 48) specimens collected in 2011 had their stomach contents preserved in 10% formalin and dissected, with contents identified and enumerated under a dissecting MVX10 microscope (Olympus Corp.; http://www. olymp us-lifes cience.com/). During dissection, the sex of each individual was determined by visual inspection of the gonads. Prey items were identified to the taxonomic groupings similar to those in other studies of Caribbean lizards (e.g., Perry, 1996): Arachnida (scorpiones), Arachnida (Araneae), Blattaria (Blattodea), Chilopoda, Coleoptera, Diptera, Ephemeroptera, Hemiptera, Hymenoptera, Isopoda, Lepidoptera, Orthoptera, and "other" for unidentifiable digested fragments. Any vertebrate remains encountered were additionally identified to the highest taxonomic resolution possible, and we additionally identified any parasites encountered in the stomach.
As formalin and alcohol preservation can have heterogeneous effects on the volume of invertebrate organisms (Donald & Paterson, 1977), enumeration of diet contents was restricted to the relative frequency. For each species, the relative frequency of each prey item was calculated based on the total number of prey items encountered across all individuals of that species within a given set of sites (i.e., sympatric sites, Table S1; Figure 2). We further collected measurement data on 79 H. mabouia for 10 morphological traits associated with feeding and locomotion: snout-vent length (SVL), postorbital width, temporalis width, head length (distance between anterior margin of tympanum to the tip of the rostrum), jaw length (distance between posterior margin of last supralabial scale to the tip of the rostrum), head height (distance between the ventral and dorsal surfaces of the head at the eye), brachium length, antebrachium length, thigh length, and shin length. All measurements were taken to the nearest 0.01 mm using digital calipers (Fowler Promax). Both stomach content and morphological data were integrated with the dataset of Dornburg et al. (2016)

Leg muscle biopsies from 21 individual Hemidactylus mabouia and 17
Phyllodactylus martini as well as 8 plant stems and leaf baseline samples were used in nitrogen (∂15N) and carbon (∂13C) stable isotope analysis. Although tail tips are commonly used in lizard isotope studies (Delibes et al., 2015), variation in lipid storage within gecko tails has the potential to confound analysis (Rode et al., 2016). Therefore, leg tissue was chosen as an alternative that also provided consistent yields of tissue required for analysis. Skin was removed from each muscle biopsy, and individual muscle and plant baseline samples were dehydrated at 40°C degrees for 48 h. Following dehydration, samples were powdered using a bead beater (MP FastPrep24 Hyland Scientific). From each sample, 1.5 mg of powder was loaded into 3 × 5 mm tins. Samples were analyzed at the University of California Davis Stable Isotope Facility using an isotope ratio mass spectrometer. As nitrogen enrichment can vary over spatial or temporal periods, quantification of trophic position for each individual was standardized using primary producer baseline samples from plant leaves and stems collected at each locality (Des Roches et al., 2016;Vidal & Sabat, 2010). To account for ∂15N values not reflecting primary producer-level values (Marshall et al., 2007), baseline samples were compared across sites with aberrant samples (i.e., primary producer ∂15N > consumer ∂15N) removed. Nitrogen values were standardized following Post (2002), in subtracting the mean ∂15N of the primary producers from ∂15N of each individual lizard and assuming fractionation of 3.4% per trophic level (Post, 2002). ∂15N values for each species were visualized using raincloud plots (Allen et al., 2019). We tested for differences between the mean ∂15N values of H. mabouia and P. martini using a Welch's t-test and additionally used Levene's test to assess whether there was a significant difference in ∂15N variance between species. To test for potential differences in ∂13C, we used the same statistical approaches as those used in the analysis of ∂15N, assuming carbon fractionation to be 0% (Post, 2002). In this case, non-significant differences in ∂13C would support the expectation that these species forage in similar habitats. To reduce the likelihood of a Type 1 error, p-values for all analyses involving multiple comparisons were adjusted using a Benjamini-Hochberg (BH) procedure (p adjusted = q; Thissen et al., 2002). All analyses were conducted in R, v. 3.4.3 (R Development Core Team, 2018).

| Stomach content analysis
Stomach contents were analyzed to test for dietary niche overlap between species that would support a hypothesis of resource competition. Differences in stomach contents between species were visualized using a principal components analysis. Relative frequencies were calculated as stated above and compared between species in sympatry and also within species between sympatric and allopatric sites. This allowed us to assess if P. martini is altering its prey base in urbanizing settings where it co-occurs with H. mabouia, and also to visualize dietary niche overlap. To test for significant differences between species, we first used an analysis of similarity (ANOSIM; Chapman & Underwood, 1999;Clarke, 1993) based on Bray-Curtis distances and 999 permutations in the vegan software package (Oksanen, 2011;Oksanen et al., 2007). Differences in mean ranks were quantified using the R statistic, for which values close to 0 indicating high similarity and values close to 1 indicating high F I G U R E 2 Map of sampling locations and representative images of habitats. Location and map of the island of Curaçao. Circles indicate sampling locations and species sympatry. (a) Sympatric area (yellow circles). (b) Allopatric "mondi" habitat with no Hemidactylus mabouia presence (blue circle). (c) Allopatric urban habitat with no Phyllodactylus martini presence (pink circles) dissimilarity (Chapman & Underwood, 1999). Additionally, we quantified the Schoener's D index (α) to assess dietary overlap between species when these co-occur with values close to 0 indicating little dietary overlap and values close to 1 indicating high dietary overlap (Schoener, 1974). A threshold of 0.60 was used to determine significant overlap (Wallace, 1981). Stomach content data were pooled from sites where species co-occur, as visualizations of dietary overlap indicated similar patterns of dietary overlap and prey consumption between sites. We also tested for differences within species by comparing their dietary overlap in sympatry versus allopatry. Males and females were pooled; Dornburg et al. (2016) previously found no significant dietary differences between P. martini males and females using these data, so only H. mabouia males and females were compared. As an ANOSIM analysis can be misled when small sample sizes are coupled with high dispersion between samples (Anderson & Walsh, 2013), these were not conducted in instances where sample sizes were lower than 20. Given that internal parasites are common in geckos (see review in Dornburg et al., 2019), we tested for differences in nematode prevalence between species using a Welch's t-test.

| Comparisons of morphology
We compared absolute differences in log snout-vent length (SVL) between species using an ANOVA and created raincloud plots (Allen et al., 2019) to visualize differences. These plots combine classic boxplots with violin raw data plots to simultaneously visualize data, the difference in size quartiles, and a kernel density smoothed estimate of the frequency distribution of the SVL data. We conducted a principal components analysis (PCA) to visualize the overall morphospace occupied by both species. In geckos, size has been shown to covary with our target morphological measurements (Dornburg et al., 2016). To account for this, we first regressed all the measurements per species against log-transformed SVL and used the residual values of individual traits regressed against log-transformed SVL as data for the PCA. To assess if differences in morphospace occupancy were mostly driven by uneven sample sizes, we randomly sampled equal numbers of both gecko species from our data 200 times across a series of datasets that reflect a range of sample sizes.
We resampled our data to assemble a series of datasets that span intervals of 5 additional samples per species starting with 10 and ending with 55 individuals per species. For each of these 2000 datasets, we conducted a PCA and computed the mean and quantiles (25% and 75%) of the ratio of H. mabouia-to-P. martini morphospace.
While morphospace visualization is advantageous for assessing the overall overlap of phenotypic variation, it is possible that allometric slopes are identical between species and simply have different intercepts (i.e., at a given body size, a focal trait in one species is larger in one species than the other). To further scrutinize our data, we used an analysis of covariance (ANCOVA) to test for differences in each morphological trait between species. For each analysis, we kept the log-transformed SVL as the covariate and treated each log-transformed morphometric measurement (e.g., jaw length, limb length, etc.) as the response. This approach allowed us to test the potential correlation for each measured trait and SVL as well as the possibility of significant differences between species that take trait covariation with SVL into account. We repeated analyses with nonsignificant interactions removed, as inclusion or omission of nonsignificant interactions can potentially impact ANCOVA analyses.
In many lizard species, including geckos, head size is a sexually dimorphic trait with males often having larger heads relative to females (Iturriaga & Marrero, 2013;Kratochvíl et al., 2002;Scharf & Mieri, 2013). Therefore, we used an ANCOVA to assess whether morphological differences for each trait were potentially masked when pooling sexes by species. For all analyses, we again kept the log-transformed SVL as the covariate and treated each logtransformed morphometric measurement as the response. Finally, we assessed potential differences in total limb lengths (brachium length + antebrachium length; thigh length + shin length) between species and sexes using log-transformed limb length as the response and log-transformed SVL as the covariate in an ANCOVA. This additional analysis facilitated additional comparisons of expectations of gecko locomotion as studies often discuss differences in total limb lengths.
Prior work has suggested large hind limbs compensate for large heads in the locomotion of Hemidactylus spp. geckos (Cameron et al., 2013). To examine scaling relationships between head size and hind limb length for both species, we constructed a set of generalized linear models (GLMs) using sex, species, SVL, and head size as explanatory variables, with one set of models using head length to quantify head size and another set using postorbital width. All models except the intercept-only null models contained an interaction term between SVL and the head size term, so that the effects of head size on limb length would be controlled for overall body size. Additional candidate models included (1) sex, (2) species identity, and (3) sex, species identity, and an interaction term between species identity and head size. Model fit was evaluated with the Akaike information criterion with a correction for small sample size (AICc). This method of model selection identifies models that predict the data well while penalizing overparameterization (Burnham & Anderson, 2004).

| Differences in feeding ecology
Analysis of ∂15N revealed a significantly (Welch's t-test: p < .004, q = 0.007, t = 3.123, df =34.272) higher mean trophic level in  (Dornburg et al., 2011). Comparing the invertebrate prey found in H. mabouia to P. martini revealed the two species to consume similar prey items, but differ in the overall percentages of prey items consumed ( Figure 4; Table S2). In comparing areas where the two spe-  Figure S2). In addition to prey contents, parasitism infestations by nematodes were significantly different between the two species (Welch's t-test: p < .001, t = −3.768, df = 71), suggesting higher parasite pressure within P. martini.

| Differences in morphology
We found a significant overall size difference between H. mabouia PC1 largely captures differences in limb lengths (~39% total hind limb and 18% total front limb) and variation in the postorbital width (~24%). In contrast, PC2 mostly captures variation in cranial measurements with over 70% of the loadings belonging to a combination of head length (~29%), jaw length (~17%), temporalis width (~13%), and postorbital width (~13%). PC3 largely captured further variation in cranial morphology (Table 4). Visualization of these PC axes revealed a high degree of overlap between species, with H. mabouia occupying more morphospace overall. Between PC1 and PC2 (Figure 5b larger. Results of our dataset resampling analyses support that these differences were not due to sample size differences alone ( Figure S3) and that mean trait values (Table 5) vary between species. SVL was significantly correlated with all measured morphological traits and ANCOVA results further support significant differences between residual trait variations after accounting for SVL scaling between species for all traits ( Table 6). The only exception to this general trend of a significant relationship between species identity and trait was head height (F = 3.232, p = .075, q = 0.075). These results were consistent whether non-significant interactions were included in  (Table S3). Tests for sexual dimorphism show no evidence for trait differences between male and female P. martini ( Figures S4 and S5). In contrast, head width was significantly different between male and female H. mabouia, suggesting H. mabouia males have wider heads than females ( Figure S4).
GLM analyses of the relationship between head size and hind limb length reveal largely concordant patterns regardless of which metric (head length or postorbital width) is used to quantify head size ( intercept of the relationship between head size and limb length for the two species, but without a difference in slope (i.e., no interaction between species identity and the head size/SVL relationship).
These top models also include no effect of sex on the relationship between head size and limb length, but in both cases the model that did include sex was also within or nearly within the set of credible models (deltaAIC of 1.52 for head length and deltaAIC of 2.2 for postorbital width).

| DISCUSS ION
The dominant hypothesis that explains the success of invasive H. mabouia populations is that they restrict access to food resources (Rocha et al., 2011;Williams et al., 2016), thereby promoting the extirpation of both native and even other non-native gecko species (Meshaka, 2000;Short & Petren, 2012

| Dietary overlap between Hemidactylus mabouia and Phyllodactylus martini
Our analyses demonstrate overlap of major invertebrate prey categories between H. mabouia and P. martini when the two species co-occur ( Figure 4; Figure S1). However, our analyses also reveal significant differences in diet between P. martini populations that are sympatric or allopatric with H. mabouia. In sympatry, the most common prey categories consumed by P. martini largely reflect common groups of invertebrates associated with human dwellings and artificial lighting in Curaçao (Dornburg et al., 2016). These categories are consistent with the diet of H. mabouia in other urbanizing areas (Bonfiglio et al., 2006;Drüke & Rödder, 2017;Iturriaga & Marrero, 2013), but contrast with the diet of P. martini in the dry mondi habitats of Curaçao. As such, dietary differences between populations of P. martini could reflect the movement of P. martini into a new foraging niche on human structures that places this species into direct contact and possible conflict with H. mabouia.
Urbanizing habitats have a pronounced effect on the foraging strategy of geckos, as the energetic cost of finding prey is reduced through utilization of artificial lights as a lure for attracting large prey resources (Gaston et al., 2013) that simultaneously confer a thermal advantage (Perry et al., 2008). This spatial clustering of food resources increases the probability of interaction and food resource competition between H. mabouia and P. martini and raises the question of whether H. mabouia has a morphological advantage for either capturing larger prey items or defending these aggregated food resource centers. We observed larger fragments of prey items such as roaches in the stomachs of H. mabouia in comparison with P. martini.
Although partially digested fragments of invertebrate body parts prohibit further testing of whether H. mabouia is exploiting larger prey, our morphological analyses lend some insight that either of these hypotheses are possible. H. mabouia possesses overall larger body sizes, as well as larger heads, hind limbs, and other traits relative to P. martini ( Figure 5; Tables 5 and 6). Increases in head height and head length are associated with increases in bite force and more efficient prey capture in geckos (Cameron et al., 2013;Massetti et al., 2017), as well as other lizard species (Dufour et al., 2018;Verwaijen et al., 2002). Functionally, this advantage is thought to arise by the combination of increasing space to accommodate increases in mandible adductor muscle sizes as well as changes in attachment angle that provide force advantages (Herrel et al., 2001). These changes in size and potential bite force poise H. mabouia to both acquire prey that may be more difficult for P. martini to capture, and better defend ambush sites from competitors. Indeed, the finding of P. martini as a prey item for H. mabouia suggests that such interactions could have very asymmetric fitness consequences. More behavioral observations in the field are warranted.
It is additionally possible that habitat-specific changes in prey base are a major determinant of the spatial distribution and the interactions of H. mabouia and P. martini that could explain our observed differences in total prey count. For example, both high numbers of terrestrial isopods in some individual H. mabouia and a higher level of nematodes likely transmitted from arthropod vectors (Dornburg et al., 2019) in P. martini suggest local variation in prey. However, our results also reveal a broad degree of dietary niche overlap within each site and a consistent pattern of overlap between sympatric habitats. This broad dietary overlap in sympatry may reflect increased homogenization of the prey base, a phenomenon commonly associated with urbanizing areas (Bang & Faeth, 2011;Fenoglio et al., 2021). Investigating the changes in species composition of Curacao's ecological communities in response to urbanization represents an urgent research frontier. Such efforts will be crucial not only for gaining additional insights into the biology of these geckos but also

| The locomotor morphology of Hemidactylus mabouia and Phyllodactylus martini
It is well known that larger heads confer prey capture and bite force advantages to lizards (Cameron et al., 2013;Massetti et al., 2017).
However, larger heads also come at a cost to locomotion as they negatively impact sprinting speed in lizards (Cameron et al., 2013;Lailvaux et al., 2019;Lopez & Martin, 2002). Our analyses suggest that increased hind limb lengths in both sexes of H. mabouia relative to P. martini may reflect the species partially mitigating a fundamental locomotor trade-off (Tables 5 and 6; Figure S5). In lizards, longer hind limbs are often correlated with increased sprint speeds and forward propulsion (Bonine & Garland, 1999;Cameron et al., 2013;Winchell et al., 2018), thereby providing an advantage for an ambush predator such as H. mabouia relying on a combination of ambush and pursuit to capture prey. A similar increased scaling in hind limb lengths has been reported in Hemidactylus frenatus (Cameron et al., 2013) and our results suggest this is potentially a general feature of Hemidactylus locomotor morphology that allows these geckos to maintain sprinting speeds and prey capture advantages.
It is intuitively appealing to view the larger front limb size of H. thereby aiding in maintaining speed and stance in downward movements (Birn-Jeffery & Higham, 2016). Quantifications of limb morphology across major lineages of geckos suggest shorter front limbs relative to hind limbs to be a hallmark of gecko locomotor morphology, with all species having between a 10 and 35% reduction in front limb proportions (Hagey et al., 2017), a finding consistent with our analysis of both H. mabouia and P. martini limb proportions (Table 5; Figure S5). Our finding of shorter front limbs relative to hind limbs in H. mabouia is consistent with expectations of selection for locomotion on steeply inclined surfaces such as walls that are coupled with large hind limbs for sprinting; but the significant negative scaling relationship between forelimb length and body size for P. martini also highlights the potential that there are additional major differences in locomotor mode and performance between these species.
Currently, the functional morphology and associated performance of P. martini remain little studied, as does that of H. mabouia in their native range in sub-Saharan Africa. These studies are of particular importance as H. mabouia is increasingly found in non-urban areas throughout its invaded range (Rocha et al., 2011), raising the question of which habitats offer the highest locomotor advantages to this competitor. While our study suggests that the vertical surfaces of dry forests in Curaçao provide habitat potentially suitable, if not advantageous, for H. mabouia, this species has yet to be found outside of urbanizing landscapes on that island. We speculate one potential reason for the absence of H. mabouia from the native bush habitat on Curaçao, and other similar desert habitats, is that all Hemidactylus geckos possess a basal toe pad system that may not be capable of successfully gaining traction on the loose, and dusty, rocky soil of the island (Russell & Delaugerre, 2017). Future analyses testing performance between and within H. mabouia and P. martini on different substrates may offer a particularly promising and exciting research frontier with high conservation applicability. For example, if H. mabouia is indeed a poor locomotor on sandy substrates, native gecko populations could be preserved integrating continual corridors of native semi-arid and arid habitat into urban planning efforts. Such efforts would yield "enemy-free" space and thereby increase the probability of the long-term persistence of native gecko species such as P. martini (Cole et al., 2005).

ACK N OWLED G EM ENTS
We Herpetology at the Yale Peabody Museum of Natural History). We would like to thank Connor Neagle for his assistance in conceptualizing the early stages of this work, and E. Ferraro, Enie Hensel, and two anonymous referees for excellent feedback on previous versions of this manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.