What makes Hemidactylus invasions successful? A case study on the island of Curaçao

Hemidactylus spp. (House geckos) rank among the most successful invasive reptile species worldwide. Hemidactylus mabouia in particular has become ubiquitous across tropical urban settings in the Western Hemisphere. H. mabouia’s ability to thrive in close proximity to humans has led to the rapid displacement of native geckos in urban areas, however the mechanisms driving this displacement remain understudied. Here we combine data from nitrogen and carbon 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. Consistent with expectations of direct food resource competition, we demonstrate substantial overlap of invertebrate prey resources between the species. Additionally, we found strong evidence from both diet content and stable isotope analyses that H. mabouia acts as a vertebrate predator, preying upon P. martini as well as other native and non-native reptiles. Finally, we 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 suggest the successful establishment of H. mabouia likely involves a combination of both exploitative interspecific competition and predation. Given the ubiquity of H. mabouia, illuminating the role of this species as both a competitor and a predator casts new concerns on the ecological and demographic impacts of this widespread urban invader.

Despite the prevalence of invasive reptiles around the world (Kraus 2015), most attention has 85 been devoted to the loss of biodiversity following the spread of a few larger bodied species such as 86 Burmese pythons (Smith et  urban and suburban environments to its ability to capitalize on the aggregation of insects around human 101 light sources (Hughes et al. 2015). Restricting the access of native geckos to these clustered food 102 resources is thought to represent a competitive advantage for H. mabouia that promotes high densities of are noted for being particularly aggressive (Short and Petren 2011), the ability of this species to 105 aggressively restrict access to spatially clustered food resources suggests interference competition, 106 whereby high densities of aggressive competitors fuel the displacement of native gecko species. However, 107 alternate hypotheses remain untested. 108 In addition to their impact as competitors, two aspects of H. mabouia invasions that have received 109 particularly little attention are locomotor morphology and role of H. mabouia as potential predators. The 110 feeding mode of H. mabouia combines ambush tactics (Vitt 1983)  individual geckos collected in areas where both species co-occur. This allowed us to test whether pooling 209 across habitat types potentially masked differences or overlap in prey items. Additionally, differences in 210 parasite load between species were compared using a Welch's t-test.

212
Comparisons of morphology 213 We compared absolute differences in log snout-vent length (SVL) between species using an 214 ANOVA and created raincloud plots (Allen et al. 2019) to visualize differences. These plots combine 215 classic boxplots with violin raw data plots to simultaneously visualize data, the difference in size 216 quartiles, and a kernel density estimate of the probability density of the SVL data. We conducted a 217 principal components analysis (PCA) to visualize the overall morphospace occupied by both species. In 218 geckos, size has been shown to covary with our target morphological measurements (Dornburg et al. 219 2016). As such, we first regressed all the measurements per species against SVL (supplemental materials) 220 and used the residual values of individual traits regressed against SVL as data for the PCA. To assess if 221 differences in morphospace occupancy were mostly driven by uneven sample sizes, we randomly sampled 222 equal numbers of both gecko species from our data 200 times in intervals of 5 between 10 and 55. For 223 each of these 2000 datasets, we conducted a PCA and computed the mean and quantiles (25% & 75%) of 224 the ratio of H. mabouia to P. martini morphospace. 225 While morphospace visualization is advantageous for assessing the overall overlap of phenotypic 226 variation, it is possible that allometric slopes are identical between species and simply have different 227 intercepts (i.e., at a given body size a focal trait in one species is larger in one species than the other). To 228 further scrutinize our data, we used an analysis of covariance (ANCOVA) to test for differences in each 229 morphological trait between species. For each analysis, we kept log transformed SVL as the covariate and 230 treated each log transformed morphometric measurement (e.g., jaw length, limb length, etc) as the 231 response. This approach allowed us to test the potential correlation for each measured trait and SVL as 232 well as the possibility of significant differences between species that take trait covariation with SVL into 233 account. We repeated analyses with non-significant interactions removed, as inclusion or omission of 234 non-significant interactions can potentially impact ANCOVA analyses. 235 In many lizard species, including geckos, head size is a sexually dimorphic trait with males often 236 having larger heads relative to females (Kratochvíl et al. 2002; Scharf and Mieri 2013; Iturriaga and 237 Marrero 2013). As such, we used an ANCOVA to assess whether morphological differences for each trait 238 were potentially masked when pooling sexes by species. For all analyses, we again kept log transformed 239 SVL as the covariate and treated each log transformed morphometric measurement as the response. 240 Finally, we assessed potential differences in total limb lengths (humerus length + radius length; femur 241 length + tibia length) between species and sexes using log transformed limb length as the response and 242 log transformed SVL as the covariate in an ANCOVA. This additional analysis facilitated additional 243 comparisons of expectations of gecko locomotion as studies often discuss differences in total limb 244 lengths.

245
Prior work has suggested large hind limbs compensate for large heads in the locomotion of 246 Hemidactylus spp. geckos (Cameron et al. 2013). As such, we examined scaling relationships between 247 head size and hindlimb length for both species by constructing a set of generalized linear models (GLMs). 248 We built models using sex, species, SVL, and head size as explanatory variables, with one set of models 249 using head length to quantify head size and another set using post-orbital width. All models except the 250 intercept-only null models contained an interaction term between SVL and the head size term, so that the 251 effects of head size on limb length would be controlled for overall body size. Additional candidate models 252 included (1) sex, (2) species identity, and (3) sex, species identity, and an interaction term between 253 species identity and head size. Model fit was evaluated with the Akaike information criterion with a 254 correction for small sample size (AICc). This method of model selection identifies models that predict the 255 data well while penalizing overparameterization (Burnham and Anderson 2004). 256 257 Results

258
Differences in feeding ecology 259 Analysis of ∂15N revealed a significant (Welch's t-test: p < 0.004; t= 3.123 df = 34.272) shift 260 towards a higher mean trophic level in H. mabouia versus P. martini (Figure 2A) individual prey items and frequency within these major groupings (Figure 3 & Table1) and ANOSIM 280 results supported significant differences between groups (R=0.213, p=0.001). Visualizations of diet data 281 based on NMDS analyses of invertebrate prey items for the two species support a large degree of overlap 282 in diet with H. mabouia utilizing the same resources as P. martini, but with H. mabouia also utilizing 283 more resources not exploited by P. martini ( Figure 4A). Repeating analyses for just individuals residing 284 in areas of co-occurrence again supported significant differences between species in an ANOSIM analysis 285 when vertebrates were included as a prey category (R=0.020, p=0.040). All instances of vertebrate 286 predation by H. mabouia were found in areas where the two species overlap (supplemental materials). 287 Restricting an ANOSIM analyses to just invertebrate prey items supported no significant differences in 288 diet between the two species (R=0.018, p=0.070). Visualizations of both the raw (supplemental materials) 289 and NMDS analyses of the invertebrate prey item diet data for the two species areas further depicted a 290 large degree of overlap ( Figure 4B). In addition to prey contents, parasitism infestations by nematodes 291 were significantly different between the two species (Welch's t-test: p < 0.001; t= -3.768 df = 71), 292 suggesting higher parasite pressure within P. martini (Table 1).

294
Differences in morphology 295 We found a significant overall size difference between H. mabouia versus P. martini (F= 10.61; p 296 = 0.00143), with H. mabouia generally being larger ( Figure 5A). Three axes of a principal components 297 analysis (PCA) of morphological traits collectively capture 64.1% of the measured variation (PC1: 298 34.84%; PC2: 16.90%; PC3: 12.37%). PC1 largely captures differences in limb lengths (~39% total 299 hindlimb, 18% total front limb) and variation in the postorbital width (~24%). In contrast, PC2 mostly 300 captures variation in cranial measurements with over 70% of the loadings belonging to a combination of 301 head length (~29%), jaw length (~17%), temporalis width (~13%), and postorbital width (~13%). PC3 302 largely captured further variation in cranial morphology (Supplemental materials). Visualization of these 303 PC axes revealed a high degree of overlap between species, with H. mabouia occupying more 304 morphospace overall. Between PC1 and PC2 ( Figure 5B) Figure 5C) and 307 PC2 & PC3 (H. mabouia CHA = 16.532; P. martini = 5.474; Figure 5D) the CHAs of H. mabouia were 308 larger. Results of our dataset resampling analyses support that these differences were not due to sample 309 size differences alone (Supplemental materials). SVL was significantly correlated with all measured 310 morphological traits (Table 2; Supplemental materials) and ANCOVA results further support significant 311 differences between residual trait variation after accounting for SVL scaling between species for all traits 312 ( Table 2; Supplemental materials). The only exception to this general trend of a significant relationship 313 between species identity and trait was head height (F= 3.232; p= 0.075). These results were consistent 314 whether non-significant interactions were included in the analysis or not (Supplemental materials). Tests 315 for sexual dimorphism for no evidence for trait differences between male and female P. martini. In 316 contrast, head width was significantly different between male and female H. mabouia, suggesting H. 317 mabouia males have wider heads than females (Supplemental materials).

318
GLM analyses of the relationship between head size and hind limb length reveal largely 319 concordant patterns regardless of which metric (head length or post-orbital width) is used to quantify head 320 size (Table 3; supplemental materials). For both measurements, the top model (lowest AICc score) was 321 the one containing a different intercept of the relationship between head size and limb length for the two 322 species, but without a difference in slope (i.e., no interaction between species identity and the head 323 size/SVL relationship). These top models also include no effect of sex on the relationship between head 324 size and limb length, but in both cases the model that did include sex was also within or nearly within the 325 set of credible models (deltaAIC of 1.52 for head length, and deltaAIC of 2.2 for post-orbital width presumably larger prey items such as roaches, beetles, and spiders were often found in the stomachs of H. 355 mabouia in comparison with P. martini. Additionally, high numbers of isopods were found in some 356 individuals. This suggests that H. mabouia could be opportunistically feeding on larger prey as well prey 357 encountered in daytime refugia. The latter could also explain the finding of a blind snake within an 358 individual H. mabouia. While partially digested fragments of invertebrate body parts prohibit further 359 testing of whether H. mabouia is more effectively harnessing larger prey, this hypothesis raises several 360 possibilities of how the natural history of these species influences differential patterns of foraging and 361 prey capture.

362
There are different responses to bright lighting between these species with Hemidactylus mabouia 363 readily Bursey 2000), further testing of differences in parasite frequencies between Hemidactylus and its native 374 competitors represents an exciting direction additional research of high relevance to animal health. 375 In addition to having an advantage in light tolerance, our analyses of trait morphological variation 376 suggest that H. mabouia has a size advantage over P. martini, possessing overall larger size, as well as 377 larger heads, hind limbs, and other traits (Fig. 4 & Table 2). Increases in head height and head length are Functionally, this advantage is thought to arise by the combination of increasing space to accommodate 381 increases in mandible adductor muscle sizes as well as changes in attachment angle that provide force 382 advantages (Herrel et al. 2001 Hemidactylus mabouia has been subjected to an unintentional experiment of introduction to 393 human mediated landscapes across the new world for centuries (Goeldi 1902 surfaces such as walls that is coupled with large hindlimbs for sprinting. But, the significant negative 411 scaling relationship between forelimb length and body size for P. martini also highlights the potential that 412 additional major differences in locomotor mode and performance between these species exist. Currently, 413 the foraging mode and activity patterns of P. martini remain little studied, as do those of H. mabouia in 414 their native range. As comparative studies of gecko functional locomotor morphology and performance 415 continue to illuminate the role of forelimbs in gecko locomotor morphology, future comparisons of 416 locomotor morphology and performance between and within these species offers a promising and exciting 417 research frontier.

419
The role of predation in Hemidactylus invasions? 420 Superiority in food resource competition has repeatedly been hypothesized as a major factor 421 facilitating the establishment of H. mabouia at the expense of native geckos (