Plasticity in thermal hardening of the invasive Asian house gecko

The Asian house gecko (Hemidactylus frenatus) is a tropical invasive species that has established and spread throughout several temperate regions around the world. In some invasive species, rapid thermal acclimation (thermal hardening) may contribute to their success in occupying a wide range of climates. In this study, we investigated whether invasive house geckos from southeastern Australia show differing thermal hardening responses than individuals from the native range in Thailand. In the laboratory, we measured the basal heat tolerance (CTmax) of the geckos and their heat hardening response after being subjected to the second thermal stress after 1, 3, 5, 7, 9, or 11 h. When geckos had recovered, we measured their basal cold tolerance (CTmin) and cold hardening response over the same time intervals. We then explored whether hardening responses differed between populations or among time intervals. We found that basal heat tolerances did not differ between populations, but geckos from Australia had lower cold tolerance than geckos from Thailand. The magnitude of the heat hardening was similar between populations, but the introduced geckos had a higher magnitude of cold hardening. The native geckos could maximize their cold tolerance capacity for only 0.6 °C, versus 0.9 °C for the introduced geckos. Also, geckos from Australia exhibited faster responses to thermal stress than did geckos from Thailand. Maximum thermal tolerances as a result of hardening responses peaked within three hours after thermal stress in Australian geckos (adjusted means = 44.0 °C for CTmax and 9.9 °C for CTmin) and at five hours after thermal stress in Thailand geckos (adjusted means = 44.2 °C and 10.2 °C, respectively). The plasticity in the thermal hardening of the invasive gecko should enable it to survive rapid temperature fluctuations, especially winter cold snaps that occur in temperate regions.


Introduction
Invasive species are a global environmental problem due to their ability to disrupt native ecosystems and cause declines or extinctions of native species (Lockwood et al. 2013). In recent decades, research has focused on accurately predicting the spread of invaders in their non-native ranges. Species distribution models (SDMs) have been widely used to predict the distribution limit of invasive species (Jiménez-Valverde et al. 2011;Rödder et al. 2008); nonetheless, many invasive species have surpassed that limit, and have spread further in non-native ranges than initially predicted (Kolbe et al. 2012;Leal and Gunderson 2012;McCann et al. 2014;Vimercati et al. 2018). In some cases, thermal plasticity may be a contributing factor for these unexpected range expansions (Kelley 2014). Several studies on invasive species have shown that individuals from non-native ranges have displayed rapid shifts in thermal biology, allowing them to function over a wide range of temperatures (Braby and Somero 2006;McCann et al. 2014;Zerebecki and Sorte 2011).
Thermal hardening is a rapid thermal acclimation within minutes or hours after a brief exposure to extreme temperature (Angilletta 2009). The exposure to high temperatures can facilitate the upregulation of heat-shock proteins (HSPs), which allow the organisms to increase their heat tolerance and survival at higher temperatures . Likewise, the exposure to cold temperatures can generate the upregulation of coldshock proteins (CSPs) (Ritossa 1962;Seebacher 2005;Thieringer et al. 1998) or induce metabolic adjustments without additional protein syntheses (Lee et al. 1987;Teets and Denlinger 2013;Teets et al. 2020), which allow organisms to survive exposure to colder temperatures. These rapid responses to thermal stresses should facilitate the survival of introduced species during the transport and early introduction phases of the invasion when they encounter unfamiliar climates (Chown et al. 2007;Nyamukondiwa et al. 2010). Dipteran insects have served as models to study thermal hardening studies for decades (Berrigan and Hoffmann 1998;Hoffmann et al. 2003;Hu et al. 2014;Loeschcke and Hoffmann 2007;Overgaard and Sørensen 2008;Sejerkilde et al. 2003;van Heerwaarden et al. 2016). By contrast, until recently, there have been fewer studies on hardening responses in reptiles (Deery et al. 2021;Gilbert and Miles 2019;Phillips et al. 2016;Refsnider et al. 2021).
The Asian house gecko (Hemidactylus frenatus) is one of the most successful species of tropical invasive reptiles that has spread throughout tropical and temperate regions worldwide (Carranza and Arnold 2006). As a tropical species, the gecko should be vulnerable to thermal extremes since it has evolved in a relatively thermal invariant climate (Ghalambor et al. 2006;Janzen 1967). Despite this, the gecko has established populations in temperate regions of Mexico, Australia, and East Asia (Farr 2011;Hoskin 2011;Kurita 2013), where ambient temperatures fluctuate widely daily, and thermal extremes (both heat and cold) are greater than those experienced in its native range. In southeastern Australia, the New South Wales (NSW) population of H. frenatus is the most southern population of the species (Hoskin 2011). This part of Australia experiences greater thermal fluctuations than does southern Thailand (Fig. 1). Southeastern Australia also experiences frequently prolonged heatwaves (Cowan et al. 2014); therefore, we predicted that this introduced NSW population of H. frenatus should exhibit higher heat hardening ability than native conspecifics. The geckos in southeastern Australia should also experience colder temperatures in winter than those experienced in the tropics; therefore, we expected to see more significant cold hardening in the introduced population. Furthermore, based on previous studies on insects, we predicted that the time course of hardening responses would differ between geckos from the introduced population and the native range. To test these predictions, we collected geckos from a native population in Thailand and an introduced population in southeastern Australia and measured their thermal hardening capacity in the laboratory.

Study populations
In this study, we chose a population from southern Thailand to represent the native population because it is located at the center of the native distribution (Ota and Whitaker 2010). For the introduced population, we chose NSW since it is the southernmost population of this species (Hoskin 2011). Although the Asian house gecko has established populations in northern Australia since the 1960s, they only reached NSW in the 2000s, making the NSW population approximately 20 years old (Hoskin 2011). We collected 60 H. frenatus from southern Thailand (Hat Yai; 7.006278,100.498871, and Satun; 6.831708, 99.5363708) from November 2018 -January 2019 to represent the native population. We collected 60 geckos from New South Wales (NSW), southeastern Australia (Yamba;-29.436890, 153.357986, and Coffs Harbour;-30.292685, 153.119707) in February 2019 to serve as the introduced population. Between 1900 and 2100 h, we walked around buildings and caught the geckos by hand, with the help of a laser pointer to attract them to reachable heights (Cole 2004). Geckos from Thailand were transported to the Prince of Songkla University (PSU), while geckos from NSW were transported to the University of Technology Sydney (UTS). We measured snout-to-vent length (SVLs) of all geckos to 0.01 mm using a vernier caliper. We started the experiments the day after arrival to the lab to minimize the effect of acclimatization. In the lab, we kept each gecko in a separate 2L ventilated plastic cage (200 × 150 × 60 mm, Sistema®, New Zealand) containing tissue paper, a water dish to reduce desiccation, and a cardboard tube as a shelter. In Thailand, we put the cages in a room with ambient temperatures and glass windows to provide natural light. In Australia, since we transferred the geckos to location further south, we set up the environments to resemble those in northern NSW. We put the cages on a shelf containing a heating cable (set to 32 °C) in a temperature-controlled room (23 °C) to allow the geckos to thermoregulate. This created a thermal gradient (23-32 °C) inside the cages which encompasses the thermal preferences of the Asian house gecko in northern New South Wales, which is approximately 29 °C (T set = 26-32 °C) (Lapwong et al. 2020). We set the light cycle at 12:12 h. We fed each gecko with five crickets every third day between 1700 and 1900 h. After the experiments, we euthanased the geckos from NSW using MS222 (Conroy et al. 2009) due to their invasive status, and released the geckos from Thailand to their sites of capture.

Thermal hardening measurement
We applied the heat hardening measurement method developed by Phillips et al. (2016). Firstly, we put a single gecko into a cylindrical plastic tube with a plastic cap and acclimated it at 23 °C for 10 min. Then we changed the cap with another one with a thermocouple already inserted before partially submerging the tube into a water bath to moderate the temperature in the tube. The thermocouple was connected to an electric thermometer (OMEGA® Thermistor thermometer-450 ATH, accuracy ± 0.1 °C) for real-time temperature measurements. We used the same tubing in both heat and cold tolerance measurements. We used a water heater (Anova Precision Cooker 2.0-Bluetooth, China) to increase, or ice to decrease, water temperatures. We controlled the change of the temperature in the tube at the rate of 1 °C per minute. We regularly rolled the tube to check the gecko's righting reflex, i.e., the ability to rotate itself over after being flipped over. When the gecko lost the righting reflex, we stopped the experiment and recorded the last temperature at which the gecko could right itself as CT max or CT min . To determine thermal hardening, we measured the critical thermal limitations of the same animal twice, which we assigned as the basal and the final CTs. While the basal CTs indicate the initial thermal tolerances of the animals, the final CTs show the maximized thermal tolerances as a result of hardening responses. Additionally, the change of CTs after the second thermal shock (i.e., differences between basal CTs and final CTs; ∆CTs) represent the animals' thermal hardening capacities. We varied the time intervals between each measurement as 1, 3, 5, 7, 9, and 11 h. There were equal numbers of male and female geckos in each treatment group (5:5), except for the 1 h interval treatment for native geckos, which had 6 males and 4 females. Each animal was measured only twice to avoid any carry-over effects. All of the measurements were conducted solely by Lapwong.

Data analyses
In the analyses, we aimed to determine the effects of locations (Thailand vs. Australia) and time intervals between two thermal shocks as treatment groups on basal CTs, final CTs, and ∆CTs. We used general linear models (GLM) to determine the effects of locations and time intervals on basal CTs, final CTs, and ∆CTs (each in a separated model). For basal CTs, we included two covariates in the model; sexes and SVLs. On the other hand, we added basal CTs as another covariate to analyze final CTs and ∆CTs because it could affect final CTs and ∆CTs as baselines. If any covariates affect the CTs, GLM will eliminate those covariates' effects and provide adjusted values for better comparisons (e.g., adjusted final CTs) (McNeil et al. 1996). Then, we used a GLM to determine the effects of locations and sexes on SVLs. Additionally, we used the Pearson correlation coefficient to evaluate the correlation between SVLs and ∆CTs, and between basal CTs and ∆CTs. Before the analyses, we used the Kolmogorov-Smirnov test to check the normality of basal CTs and the residuals of the CTs. We found that 95% of our data sets were normally distributed (P > 0.05). The Levene's test of equality of variances also confirmed the homogeneity of variances (P > 0.05). We performed all statistical analyses using IBM® SPSS® Statistics Version 24.

Heat tolerance
The general linear model revealed no effect of locations (F 1,106 = 1.520, P = 0.22), treatment groups (F 5,106 = 0.954, P = 0.45), and interaction between locations and groups (F 5,106 = 0.179, P = 0.97), on the basal CT max . Sexes and SVLs did not have significant effects as covariates (F 1,106 = 0.593, P = 0.44, and F 1,106 = 1.403, P = 0.24, respectively). The mean ± SD basal CT max were 43.6 ± 0.7 °C and 43.6 ± 0.5 °C for native and introduced geckos, respectively. There was a significant negative correlation between basal CT max and ∆CT max (r =− 0.474, P < 0.05). The model revealed no effect of locations (F 1,105 = 0.453, P = 0.50), but a significant effect of time interval (F 5,105 = 2.416, P < 0.05), and a significant interaction between time and location (F 5,105 = 3.762, P < 0.05) on final CT max . That is, the time course of heat resistance differed between locations (Fig. 2a).). Among three covariates, only basal CT max had a significant effect on final CT max (F 1,105 = 51.794, P < 0.05), but not sexes (F 1,105 = 0.270, P = 0.61) and SVLs (F 1,105 = 0.470, P = 0.50). There was also no effect of locations (F 1,105 = 0.192, P = 0.67), but a significant effect of time interval (F 5,105 = 2.914, P < 0.05), and a significant interaction between time and location (F 5,105 = 5.083, P < 0.05) on ∆CT max . That is, the time course of heat hardening differed between locations (Fig. 3a). Among three covariates, only basal CT max had a significant effect on ∆CT max (F 1,105 = 31.920, P < 0.05), but not sexes (F 1,105 = 0.039, P = 0.84) and SVLs (F 1,105 = 0.547, P = 0.46). For instance, the Thailand geckos had the highest final CT max at the 5-h interval (adjusted final CT max = 44.2 °C, adjusted ∆CT max = 0.6 °C), whereas the NSW geckos had the highest final CT max at the 1 h interval (adjusted final CT max = 44.0 °C, adjusted ∆CT max = 0.4 °C).

Cold tolerance
The general linear model revealed a significant effect of location on basal CT min (F 1,106 = 18.414, P < 0.05), and there was no effect of treatment group (F 5,106 = 1.151, P = 0.34) or the interaction between location and group (F 5,107 = 0.284, P = 0.92). While sexes had a significant effect as a covariate (F 1,106 = 10.483, P < 0.05), SVLs did not 1 3 (F 1,106 = 0.140, P = 0.71). The mean ± SD basal CT min were 11.2 ± 0.8 °C for Thai geckos and 10.4 ± 0.8 °C for Australian geckos. Male geckos had higher basal CT min and female geckos in both locations (11.4 °C vs. 11.0 °C in Thailand, and 10.8 °C vs. 10.1 °C in Australia). There was a significant positive correlation between basal CT min and ∆CT min (∆CT min ; r = 0.554, P < 0.05). The model also showed significant effects of locations (F 1,106 = 5.945, P < 0.05), time intervals between cold shocks (F 5,106 = 2.494, P < 0.05), and interaction between time intervals and locations (F 5,106 = 5.537, P < 0.05) on final CT min . That is, the time course for cold hardening differed between native and introduced geckos   (Fig. 2b). Among three covariates, only basal CT min had a significant effect on final CT min (F 1,105 = 13.511, P < 0.05), but not sexes (F 1,105 = 0.092, P = 0.76) and SVLs (F 1,105 = 1.918, P = 0.17). There were also effects of locations (F 1,105 = 5.945, P < 0.05), time intervals between cold shocks (F 5,105 = 2.494, P < 0.05), and interaction between time intervals and locations (F 5,105 = 5.537, P < 0.05) on ∆CT min . That is, the time course of cold hardening differed between locations (Fig. 3a). Among three covariates, only basal CT min had a significant effect on ∆CT min (F 1,105 = 60.414, P < 0.05), but not sexes (F 1,105 = 0.092, P = 0.76) and SVLs (F 1,105 = 1.918, P = 0.17). For instance, while the native geckos had the lowest final CT min at the 5-h interval (adjusted final CT min = 10.2 °C, adjusted ∆CT min = − 0.6 °C), the introduced geckos had the lowest final CT min at the 3 h interval (adjusted final CT min = 9.9 °C, adjusted ∆CT min = − 0.9 °C).

Protein kinetics
In Phillips et al. (2016), the author suggested that their heat hardening graphs fitted quadratic equations. However, in other works, researchers showed that heat hardening consists of two stages. First, after being stimulated, the organism quickly produced heat shock proteins. Then, the protein level slowly declined after reaching the peak (Richter et al. 2010;Simões-Araújo et al. 2003). Therefore, the first part of the graph (before the peak) should be explained by a quadratic equation, while the second part of the graph (after the peak) should be explained by a logarithm equation. Accordingly, we used quadratic equations to explain the thermal hardening responses of the native population because the geckos have increased their thermal tolerances during the first 5 h before decreased after that. On the other hand, we used logarithms equations to explain the introduced populations' responses to thermal shocks because their thermal hardenings have peaked within the first hour and gradually decrease after that. The equations fitted well with the graphs with R 2 ranging from 0.6058 to 0.9186 (Fig. 3).

Discussion
The Asian house geckos have established in the temperate region of southeastern Australia for more than 20 years (Hoskin 2011), so we expected them to have more substantial thermal hardening capacity, especially cold hardening, than their native conspecifics. As predicted, we found a greater degree of cold hardening, but not heat hardening, of the introduced populations. For cold tolerance, geckos from NSW had lower basal CT min and final CT min than geckos from Thailand. The analysis of covariance confirmed that the lower final CT min was beyond the effect of the lower basal CT min . The geckos from NSW could shift their cold tolerance to a greater degree than the geckos from Thailand. Our findings for both basal cold tolerance and cold hardening agree with the results reported for other species of invasive lizards that have shifted cold tolerance downwards following successful spread to colder regions (Angetter et al. 2011;Kolbe et al. 2012;Leal and Gunderson 2012). The ability to withstand low temperatures of this invasive gecko also suggests that they do not always rely on the warm climate inside the buildings and can spread into natural habitats. Furthermore, we found that female geckos had lower basal CT min . Although the study by Cameron et al. (2018) found that male H. frenatus tended to have a wider performance breadth, they also observed a higher level of activities in males, presumably induced by mating and territory defense behaviors. Therefore, we suggest that male geckos might keep thermoregulating, so the ability to tolerate low temperatures has become less vital. Another study in the Italian Wall Lizard (Podarcis siculus) found a similar phenomenon to ours, whereby females could tolerate colder temperatures than males, suggesting that this sexual divergence in cold tolerance could correlate with different microhabitat use or physiology between sexes (Liwanag et al. 2018).
For heat tolerance, basal CT max and ∆CT max were similar in geckos from Thailand and NSW. These results contrast with those from a study on another tropical reptile, Lampropholis coggeri, which found an interpopulation divergence in heat hardening (Phillips et al. 2016). In that study, the magnitude of heat hardening diverged among populations and was higher for skinks from localities with higher seasonal variation in daily maximum undercanopy temperatures. By contrast, we found little interpopulation variation in either basal heat tolerance or the magnitude of heat hardening. This finding could be due to the already high initial heat tolerance, so there is no evolutionary reason for the geckos to uplift their heat tolerance. Interestingly, while the magnitude of the hardening response (around 0.7 °C) is similar to that reported for other geckos (e.g., Amalosia lesueurii, Abayarathna et al. (2019)), the basal CT max of H. frenatus (43.6 °C) is much higher than that reported for most other gekkonids (i.e., mean of 40.8 °C, (Clusella-Trullas and Chown 2014)). Alternatively, perhaps the CT max of geckos may have already reached an upward limit, such that there is little opportunity for further upward shifts. In support of this idea, there was a negative correlation between the basal CT max and ∆CT max , a finding that was also reported for tropical skinks (Phillips et al. 2016) and anoles (Deery et al. 2021). That is, individuals with low basal CT max showed higher hardening responses and vice versa. This pattern mirrors finding in other taxa, such as Drosophila, and suggests that there is a hard upper limit to shift thermal tolerance upwards (Hoffmann et al. 2003;van Heerwaarden et al. 2016).
Besides the magnitude of the hardening response, we documented apparent differences in the time courses of plastic responses to thermal tolerance. Interestingly, geckos from NSW responded faster to both heat shocks and cold shocks than did geckos from Thailand. This result mirrors findings from comparative studies on invasive insects and their congeneric species with limited distributions. For example, the cosmopolitan Mediterranean fruit fly, Ceratitis capitata, responded to the thermal stresses faster than the less successful invasive congener, C. rosa (Nyamukondiwa et al. 2010). In another study on fruit flies, 3 rd instar larvae of the widespread invasive Bactrocera dorsalis and more geographically restricted B. correcta were allocated to groups subjected to different temperatures (25, 30, 35, 37, 39, 41 °C) followed by exposure to 45 °C. Interestingly, larvae of the invasive fly had higher survival after exposure to milder temperatures (35 °C and upwards). In contrast, the non-invasive B. correcta only showed a heat shock response after exposure to temperatures of 39 °C and above (Hu et al. 2014). While these results document clear differences in the plasticity of thermal hardening responses between invasive and non-invasive species, as far as we are aware, few studies have compared hardening responses in the introduced populations to the native ones. Thus, it is possible that invasive flies may develop those thermal sensitivities as the post-invasion response to novel climates.
In the current study, we suggest that plasticity in thermal hardening in house geckos has arisen in response to the colder climate and the strong predictability in thermal extremes in southeastern Australia. Even in temperate NSW, thermal fluctuations can range from 12.9 to 43.3 °C during summer (Fig. 1b). In the native range, geckos are rarely exposed to critical temperatures (annual thermal fluctuation ranges from 18.2 to 39.2 °C, Fig. 1a), so there may be little benefit to shifting down their cold tolerance or reacting quickly to exposure to thermal stress. On the other hand, there is a higher chance for the introduced geckos to experience critical or near-critical temperatures, so faster physiological responses should enhance their survival. It is also possible that plasticity for thermal hardening occurred during the transport phase of the invasion. For example, a study on marine bivalves found that exposure to high thermal stress during simulated transport promoted strong selection for enhanced survival upon exposure to the second thermal stress (Lenz et al. 2018). Irrespective of when the shift in plasticity in hardening occurred in H. frenatus, it is likely that invasive geckos would benefit from rapid hardening responses during regular heatwaves or cold snaps, both of which frequently occur in southeastern Australia. However, without further molecular studies, we do not know whether the shifts in thermal hardening responses of the Asian house geckos in NSW were the result of natural selection or just non-inherited plasticity.
In summary, we found that H. frenatus from NSW had greater cold tolerance, and responded faster to both heat and cold stresses, than did geckos from Thailand. Such plasticity to thermal stress is likely to influence the survival of individuals and should facilitate further invasions, especially in the temperate zone. However, since this study included only two populations, additional studies of H. frenatus from different geographic ranges are required to test this prediction.