Beyond survival experiments: using biomarkers of oxidative stress and neurotoxicity to assess vulnerability of subterranean fauna to climate change

Accurate assessments of species vulnerability to climate change need to consider the physiological capacity of organisms to deal with temperature changes and identify early signs of thermally induced stress. Oxidative stress biomarkers and acetylcholinesterase activity are useful proxies of stress at the cellular and nervous system level. Such responses are especially relevant for poor dispersal organisms with limited capacity for behavioural thermoregulation, like deep subterranean species. We combined experimental measurements of upper lethal thermal limits, acclimation capacity and biomarkers of oxidative stress and neurotoxicity to assess the impact of heat stress (20 ◦ C) at different exposure times (2 and 7 days) on the Iberian endemic subterranean beetle Parvospeonomus canyellesi . Survival response (7 days of exposure) was similar to that reported for other subterranean specialist beetles (high survival up to 20 ◦ C but no above 23 ◦ C). However, a low physiological plasticity (i.e.incapacitytoincreaseheattoleranceviaacclimation)andsignsofimpairmentatthecellularandnervoussystemlevelwere observed after 7 days of exposure at 20 ◦ C. Such sublethal effects were identified by significant differences in total antioxidant capacity, glutathione S-transferase activity, the ratio of reduced to oxidized forms of glutathione and acetylcholinesterase activity between the control (cave temperature) and 20 ◦ C treatment. At 2 days of exposure, most biomarker values indicated some degree of oxidative stress in both the control and high-temperature treatment, likely reflecting an initial altered physiological status associated to factors other than temperature. Considering these integrated responses and the predicted increase in temperature in its unique locality, P. canyellesi would have a narrower thermal safety margin to face climate changethanthatobtainedconsideringonlysurvivalexperiments.Ourresultshighlighttheimportanceofexploringthermally sensitive processes at different levels of biological organization to obtain more accurate estimates of the species capacity to face climate change.


Introduction
Species vulnerability to climate change depends on the capacity of individuals to maintain current populations or to shift geographical ranges to future suitable environments (Williams et al., 2008). These processes are ultimately related to the niche breadth (Kearney, 2006) and dispersal capacity. However, most studies forecasting species responses to climate change do not account for their actual environmental tolerances, as relevant physiological traits are not explicitly considered (Kearney et al., 2010). Experimental measurements of such traits, as proxies of the species' fundamental niche breadth (sensu Hutchinson, 1957), have the potential to improve species vulnerability assessments (Somero, 2010;Arribas et al., 2012).
Amongst the multiple dimensions of species niches, the 'thermal niche' (i.e. the range of body temperatures maintaining positive population growth; Gvoždík, 2018) is key for small ectotherms, as ambient temperature affects most aspects of individual performance and fitness (Chown and Nicolson, 2004;Chown, 2012;Gonzalez-Tokman et al., 2020). Some studies have combined climatic and experimental data to calculate thermal safety margins (TSM), i.e. the difference between habitat temperature and the species' critical thermal maxima (CT max ) (e.g. Diamond et al., 2012;Kellermann et al., 2012;Kellermann and van Heerwaarden, 2019). Such CT max are typically estimated from survival experiments and assumed to be reliable proxies for the species' capacity to deal with warming. However, thermal sensitivity often occurs in a hierarchical manner, such that processes most sensitive to environmental change can limit the overall fitness of an organism, and survival is often possible over a wider range of temperatures than locomotion, development or reproduction (Buckley and Kingsolver, 2012;Evans et al., 2015). Otherwise, capacity for thermal acclimation may extend critical limits for performance or survival (Seebacher et al., 2015).
Many molecular and biochemical mechanisms are triggered under thermal stress in ectotherms (Tattersall et al., 2012). Amongst them, high temperature stress increases the production of cellular reactive oxygen species (ROS), leading to oxidative stress (Tomanek, 2015). Oxidative stress is defined as an imbalance of oxidant molecules and antioxidant mechanisms in favour of the former, which leads to disruption of redox signalling and control as well as biomolecule damage (Sies and Jones, 2007). The organism has evolved complex molecular and enzymatic defence systems against excessive ROS production to maintain the cellular redox homeostasis. It is generally assumed that a limited antioxidant capacity leads to poor acclimation potential (Madeira et al., 2018), and ultimately, limited health status and low capacity to cope with environmental stressors (e.g. Yang et al., 2010;Madeira et al., 2013;Enzor and Place, 2014). Thermal stress has also a significant impact on locomotion, an energetically costly trait and major component in the individual fitness and organisms' ability to cope with different environments (Domenici et al., 2007). In locomotion, cholinergic synapses are essential in the neuromuscular transmission, and the enzyme acetylcholinesterase (AChE, EC 3.1.1.7) plays a pivotal role in the regulation of the nervous impulse via hydrolysis of the neurotransmitter acetylcholine. This highly temperature-sensitive enzyme (Baldwin, 1971;Scaps and Borot, 2000;Pfeifer et al., 2005;Dimitriadis et al., 2012;Attig et al., 2014) can be used as an indicator of physiological stress at the nervous system level. Together with oxidative stress parameters, it provides a potentially useful biomarker to identify early signs of thermal stress and altered metabolic processes. Such estimates may help to interpret outcomes from survival experiments in thermal stress research and to refine assessments of species vulnerability to climate change (e.g. Deschaseaux et al., 2010;Goh and Lai, 2014;Madeira et al., 2016;Campos et al., 2019).
One of the main limitations for predicting species responses to climate change is the difficulty to account for species dispersal capacity and the ability to buffer from stressful temperatures by behavioural adjustments (Jiménez-Valverde et al., 2008;Sánchez-Fernández et al., 2011, 2012. Exceptionally, such limitations inherent to most environments are minimized in deep subterranean habitats . These habitats are characterized by extremely reduced environmental variability across both time and space (Lauritzen, 2018); so, deep subterranean organisms have a limited behavioural capacity to exploit different microhabitats compared to most surfacedwelling species. Furthermore, these species are generally poor dispersers, usually with restricted ranges and isolated populations (Culver and Pipan, 2019;Ribera et al., 2018). This means that their fate under climate change will depend almost entirely on their physiological thermal tolerance and the capacity to persist in their current locations, which somehow simplifies vulnerability assessments, but places subterranean fauna as a potentially highly endangered  (Sánchez-Fernández et al., 2016, Mammola et al., 2019a. Accordingly, experimental estimates on the temperature effects at the physiological level and on individual performance are essential if we aim to get accurate predictions on the impact of global warming in subterranean fauna. Previous studies have shown that ectotherm subterranean species have narrow thermal tolerance ranges and limited acclimation capacities compared to surface-dwelling ones, some of them being extremely sensitive to thermal changes (e.g. Mermillod-Blondin et al., 2013;Di Lorenzo and Galassi, 2017;Mammola et al., 2019c). However, the upper lethal thermal limits of a number of subterranean species are well above their current habitat temperatures and those experienced through their evolutionary history (e.g. Issartel et al., 2005;Rizzo et al., 2015;Sánchez-Fernández et al., 2016;Pallarés et al., 2019). Yet, little is known about the underlying biochemical mechanisms that modulate thermal tolerance in these species as well as the physiological effects of exposure to sublethal temperatures (but see Colson-Proch et al., 2010;Bernabò et al., 2011). In this study, we combine experimental measurements of lethal thermal limits, acclimation capacity and biomarkers of oxidative stress (total antioxidant capacity, glutathione S-transferase activity, glutathione concentration and lipid peroxidation) and neurotoxicity (AChE activity) to assess the impact of heat stress on the subterranean beetle species Parvospeonomus canyellesi (Lagar, 1974) (fam. Leiodidae). We also accounted for the effect of exposure time in biomarker responses to identify transient and long-term biological adjustments to thermal stress. Some studies have shown that certain oxidative stress biomarkers (e.g. antioxidant enzyme activities) respond faster to environmental stressors than others (e.g. lipid peroxidation or DNA damage) in both vertebrates (Bouchard, 2003) and invertebrates (Colacevich et al., 2011), and the dynamic of their responses depends on the level of exposure (time × dose) to the stressors. After first short-term defence activation, an organism under thermal stress may either exhaust the cellular protection mechanisms or acclimate (Madeira et al., 2016b). Furthermore, transport and the time of maintenance of subterranean organisms in laboratory conditions may be an additional source of variation that could affect outcomes when measuring highly stress-sensitive parameters. We hypothesize that exposure to high non-lethal temperatures in subterranean insects triggers oxidative stress and locomotor impairment, consequently reducing their effective TSM under climate change, and that such effects might be time-dependent.
The study species shows an extremely limited geographic distribution, as it is found in a single cave in NE Spain. Beyond the interest in this particular, potentially endangered species, we seek to take advantage of the potential of subterranean environments for ecological, biogeographical and conservation research (Sánchez-Fernández et al., 2018;Mammola et al., 2019b) to both improve the understanding of the mechanisms underlying tolerance of insects to heat stress and to obtain more accurate estimates of the response of biodiversity to climate change.

Study site and species
Parvospeonomus canyellesi is a troglobiont (i.e. obligate subterranean) species endemic to Spain and cited in a single cave ('Forat de las Pedreres') (Fresneda and Salgado, 2016). This cave is located in a mining area in the Montseny mountainous massif, a Catalan pre-coastal mountain range designated as Natural Park and UNESCO Biosphere Reserve. The species is neither considered of conservation concern nor has any national or international legal protection despite of its extremely reduced geographical distribution and high habitat specialization. It shows several traits associated with the subterranean specialization in Leptodirini beetles, e.g. complete loss of eyes, elongated appendages and a modified life cycle with two larval instars, observed only in strictly subterranean species and different from the three-instar cycle typical in most Coleoptera (Cieslak et al., 2014a, b;Ribera et al., 2018). However, it does not exhibit other extreme modifications observed in related species specialized in deep subterranean habitats (e.g. increased body size and the most modified life cycle with a single larval instar), which may suggest an intermediate degree of subterranean specialization for this species.
Temperature in the deep parts of caves is relatively constant through the day and year (Badino, 2010;Mammola & Isaia, 2016), and it can be indirectly estimated from the mean annual temperature of the surface, as both values are highly correlated (Poulson and White, 1969;Badino, 2004;Mammola et al., 2017;Sánchez-Fernández et al., 2018). Accordingly, we obtained the current surface mean annual temperature at the study site from the WorldClim database (WorldClim v.1.4; http://www.worldclim.org). We also obtained the predicted mean annual temperature for future climatic scenarios (2070) under the most optimistic (2.6) and pessimistic (8.5) representative concentration pathway (RCP). We averaged the estimates from all the general circulation models available in WorldClim for each RCP, to account for the uncertainty associated to different models. Both current and future temperatures were obtained at 30-arc-second spatial resolution. The current surface annual temperature at the cave location is 13.7 • C, and it is predicted to increase up to 15.5 (RCP 2.6) or 17.4 • C (RCP 8.5) by 2070.

Animal collection and maintenance
Adult specimens of P. canyellesi (90 specimens approx.) were collected in Forat de las Pedreres and transported to the laboratory at cool and wet conditions in a portable fridge, using the substrate from the cave and moss to keep a high humidity (>90% RH). Sampling permits were previously approved by the corresponding local authority ( permit numbers SF/994-SF/1001). In the laboratory, groups of 10 specimens were placed in different 10 × 15-cm plastic boxes with a white plaster substratum, 2-3 small volcanic stones and tissue paper (e.g. Rizzo et al., 2015;Pallarés et al., 2019). These were covered with holed plastic film to maintain RH values close to saturation but allow aeration. Containers were placed in a laboratory incubator (Radiber ERF-360, Radiber S.A., Barcelona, Spain) at the approximate temperature of the study site (13 • C) and permanent dark for 2 days previous to the experiments. For the entire duration of the experiments, specimens were fed ad libitum with freshly frozen Drosophila melanogaster, and relative humidity was kept >90% by wetting the substrate, stones and papers daily and placing trays with distilled in the incubators. Temperature and relative humidities inside the containers were recorded every 5 min using HOBO MX2301 data loggers (Onset Computer Corporation, Bourne, MA, USA).

Estimation of upper lethal temperature and acclimation capacity
To determine the upper lethal temperature (ULT) of P. canyellesi, we measured survival at different temperatures, which represent its current natural conditions (13 • C, used as a control treatment), and potentially sublethal (20 and 23 • C) or extreme temperatures (25 • C) according to the ULT reported for species of the same clade (between 20 and 23 • C) (Rizzo et al., 2015). Ten specimens, placed in a plastic box were subjected to each treatment in individual climatic incubators, for 1 week. The rest of the experimental conditions were maintained as described above. Survival was monitored every 24 h. After 1-week exposure, specimens from the 13 and 20 • C treatments (non-lethal temperatures) were exposed to 23 • C to look at acclimation capacity (i.e. to assess whether survival under heat stress is influenced by previous thermal conditions), and survival was monitored daily until all specimens were dead (i.e. 8 days).

Biomarker assays
Other set of specimens was exposed to the control (13 • C) or high temperature (20 • C) for either 2 (short-term exposure) or 7 days (long-term exposure) (n = 10 specimens per each temperature × time treatment), in order to identify potential sublethal effects of heat stress and the time dependence of the responses. The selected temperatures closely represented current and future extreme climate scenarios in the study site, respectively (see above). The short-term, 2-day exposure is a typical time point used in oxidative stress bioassays to assess the first antioxidant defence activation under stress (Saint-Denis et al., 1999;Bouchar, 2003;Colacevich et al., 2011). The 7 days of exposure was used as a longer-term exposure to relate biomarker responses with acclimation capacity, because thermal acclimation responses have been found in some subterranean invertebrates within such time (Pallarés et al., 2019 and unpublished data).
Specimens were immediately frozen at −80 • C at the end of each treatment and individually homogenized in 300 μl of ice-cold 10-mM sucrose buffer (pH = 7.4) containing 1 mM EDTA, using a plastic tissue grinder (kept on ice) and an ultrasonic cell disruptor. The homogenates were centrifuged at 10 000 × g for 5 min at 4 • C, and the post-mitochondrial fraction was aliquoted and stored at −80 • C until biomarker analyses. Total protein content was determined using Bradford method (Bradford, 1976) adapted to a 96-well microplate, and using bovine serum albumin as the standard.
We measured the following biomarkers of oxidative stress: the total antioxidant capacity (TAC), glutathione Stransferase activity (GST, EC 2.5.1.18), concentration of reduced (GSH) and oxidized (GSSG) forms of glutathione and lipid peroxidation (LPO). TAC is a global measure of the antioxidant compounds present in the sample. This biomarker was measured following the protocol by Erel (2004). GST activity participates not only in the metabolism of xenobiotic compounds (van der Oost et al., 2003), but also in the ROS inactivation (Hayes et al., 2005). It was measured by the method described by Habig et al. (1974). Glutathione is a tripeptide highly reactive against ROS (van der Oost et al., 2003). Its oxidation by ROS contributes to reduce the cellular oxidative stress, so a decrease in the GSH:GSSG ratio is commonly considered as an indicator of oxidative stress (Storey, 1996;Asensi et al., 1999;Amaral et al., 2012). GSH and GSSG concentrations were measured following the method by Rahman et al. (2006). LPO was measured as an indicator of oxidative damage because it is the result of the interaction between ROS and polyunsaturated fatty acids, leading to severe injury of plasma membranes (Halliwell and Chirico, 1993;Halliwell and Gutteridge, 1999). LPO was estimated by the formation of thiobarbituric acid reactive substances (TBARs) according to the method by Agarwal and Chase (2002). We also measured AChE activity as a biomarker of stress at the nervous system level, according to the method by Ellman et al. (1961). Detailed information of the protocols used to measure each parameter is provided as Supplementary Material (Appendix S1).

Data analysis
We used Kaplan-Meier survivorship curves (Altman, 1992) for a visual comparison of survival at the different tested temperatures. Right censored data were specified for those individuals that were alive at the end of the 7 days of exposure (see Therneau, 2015). The overall effect of treatment on survival time was assessed by a log-rank test (Harrington & Fleming, 1982). Differences amongst treatments were determined by pairwise comparisons using also log-rank tests with Bonferroni-adjusted P values. The same procedure was used to determine acclimation capacity, i.e. to determine the effect of acclimation temperature (13 vs. 20 • C) on survival time under the subsequent exposure to 23 • C. To estimate ULTs, survival data at the end of the 7 days of exposure at 13-25 • C were fitted to a logistic regression model from which LT 50 values (median lethal temperature, i.e. the temperature at which 50% of the exposed individuals have died) were obtained.
To assess the impact of temperature on the biomarkers measured at short-(2 days) and long-term exposure (7 days), we used the integrated biomarker response (IBR) index, the IBRv2, which is a modified version by Sanchez et al. (2013) of the original IBR index (Beliaeff and Burgeot, 2002). The index integrates the response of multiple biomarkers to provide a measurement of the organism' physiological status at specific environmental conditions. Star-plots of the A i -score outcome (individual biomarker deviation index) obtained from IBR calculation allowed the visual assessment of the response (inhibition or induction) of each biomarker (see Supplementary Material, Appendix S2). The specific responses of individual biomarkers to temperature and exposure time were also explored by performing a two-way ANOVA test for each biomarker, including temperature, time and their interaction as predictors. When significant effects were detected, pairwise post hoc tests with Bonferroni correction were used to compare biomarker values between treatments.

Upper lethal limits and acclimation capacity
Temperature had a significant effect on survival time (logrank test: χ 2 = 47.6, df = 3, P < 0.001). Most of the specimens survived 7 days after exposure to both 13 and 20 • C (Bonferroni-adjusted P = 0.31, df = 1). Survival times in the 23 and 25 • C treatments were significantly different from the other treatments (all Bonferroni-adjusted Ps < 0.001, df = 1). Exposure to 23 • C caused a progressive decrease of animal survival, which led to 100% mortality at the end of the survival experiment. No specimens survived longer than 1 day at 25 • C (Fig. 1). The estimated LT 50 at 7 days was 21.13 ± 0.49 • C.

Biomarker responses
IBRv2 values were 1.04 and 4.40 for 2 and 7 days of exposure, respectively, showing a higher impact of temperature on biomarker responses at longer than short-term exposure. The star plots of the A i -scores evidenced an increase of TAC and GST activity at both exposure times, as well as an increase of AChE activity and a marked and slight decrease in the GSH:GSSG ratio and LPO, respectively, at longer-term exposure (Fig. 2).
Individually, most biomarkers (GST, GSH:GSSG ratio and AChE) showed a time-dependent response to temperature (significant interaction temperature × time, see Supplementary Material, Table S1). TAC increased with both temperature and exposure time (Fig. 2a, Table S1). GST and AChE activity increased with temperature, but such effect was only significant after 7 days of incubation ( Fig. 2b and e, Table S1). The GSH:GSSG ratio was significantly lower at short-than long-term exposure, and at 7 days, it was lower at 20 • C than 13 • C (Fig. 3c, Table S1). TBAR levels were lower at short than long-term exposure, but not significantly affected by temperature at any exposure time (Fig. 2d, Table S1).

Discussion
Adult specimens of P. canyellesi showed similar heat tolerance in the survival experiment than other species of the same clade, whose ULT were measured under the same approach (i.e. no survival above 23 • C after 7 days of exposure) (Rizzo et al., 2015). This species can withstand temperatures much higher (around >6 • C) than that of its habitat, which is also in agreement with other studies of heat tolerance in several subterranean taxa (e.g. Issartel et al., 2005;Mermillod-Blondin et al., 2013;Rizzo et al., 2015;Pallarés et al., 2019). However, acclimation and biomarker experiments indicated that P. canyellesi might have difficulties in coping with rising temperature due to a combination of low physiological plasticity and signs of heat-induced impairment at the cellular and nervous system level.

Acclimation capacity and sublethal effects of heat stress
P. canyellesi did not show an enhancement of heat resistance after previous acclimation at 20 • C during 7 days. In agreement with such acclimation failure, biomarker responses revealed sublethal effects of high temperature exposure, of much greater magnitude at long than short-term exposure (4-fold higher IBR v2 values). High IBRv2 values, reflecting altered physiological status, have been associated with low acclimation abilities (Madeira et al., 2016(Madeira et al., , 2018Campos et al., 2019). Specifically, after 7 days of incubation, the significant increase of TAC and GST activity in parallel with the decrease of GSH:GSSG ratio in the 20 • C treatment, compared to controls, suggest activation of both molecular and enzymatic antioxidant mechanisms. Despite the lack of acclimation capacity of the species, such antioxidant response appeared to be effective in controlling irreversible damage to membranes, since LPO levels did not differ between the control and stressful treatment at this exposure time. It has been long reported that exposure to chemical stressors, such as environmental contaminants, cause an imbalance in the cellular oxidative status that lead to a decreased GSH:GSSG ratio and activation of GSH, amongst other enzymatic antioxidant mechanisms (Storey, 1996;Asensi et al., 1999;Amaral et al., 2012). Therefore, our findings suggest a quick response of glutathione-dependent antioxidant mechanisms against heat stress; the sensitivity and magnitude of such response could be critical for the adaptation potential of the study species to global warming. Whether this response is similar in other subterranean species remains to be investigated, as glutathione metabolism appears to be highly species-specific (Enayati et al., 2005;Domingues et al., 2010;Rane et al., 2019).
Response of AChE activity indicated heat-induced stress also at the nervous system at the long-term, as the activity of this enzyme was higher at 20 than 13 • C at 7 days. Such increase of AChE activity may reflect higher locomotor activity under high temperature (e.g. Kjaersgaard et al., 2010). Indeed, previous studies have found a positive correlation between AChE activity and locomotor behaviour parameters in invertebrates (e.g. Engenheiro et al., 2005;Xuereb et al., 2009), showing the significant physiological role of this esterase in animal locomotion activity. Although locomotor parameters were not systematically registered in our experiment, we observed that specimens were generally much more active in the 20 • C treatment, supporting such hypothesis. However, such up-regulation of AChE under thermal stress contrasts with other studies, which have shown the opposite pattern (i.e. an inhibition of AChE activity) (e.g. Scaps and Borot, 2000;Dimitriadis et al., 2012;Attig et al., 2014). It could not be discarded that a more prolonged exposure or higher temperatures would down-regulate the AChE activity in P. canyellesi.
Recent studies on various subterranean taxa suggest that the physiological responses underlying thermal tolerance could vary between species with different degree of specialization to the subterranean environment (e.g. Raschmanová et al., 2018;Mammola et al., 2019c). Amongst cave beetles, some deep subterranean specialist Leiodidae appear to lack thermal plasticity, like P. canyellesi, (i.e. they show no capacity for acclimation to high temperature), whilst related species less specialized to subterranean life (facultative subterranean) have some degree of acclimation capacity (Pallarés et al., unpublished data). Bernabò et al. (2011) demonstrated the occurrence of a heat shock response in two obligate subterranean leiodids, but its intensity was lower in the species confined to the internal-more thermally stable parts of caves than in its congener that lives closer to cave entrances, exposed to higher thermal fluctuations. Therefore, in the evolutionary process of specialization to the highly stable deep subterranean environments, the physiological capacity to cope with heat stress might have  been reduced. The capacity to onset antioxidant responses under heat stress observed here may allow cave beetles to survive at temperatures well-above that of the habitat for a limited exposure time, as seen here and in related species (Rizzo et al., 2015). However, given the high energetic cost of such responses (Krebs and Loeschke, 1994;Tomanek, 2010), it is likely that prolonged exposure to heat stress in resource-limited environments such as caves decreases fitness in the longer time scale (Monaghan et al., 2009). Further work exploring these issues (e.g. transcriptome and proteome-level responses, metabolic rates and the fitness consequences of long-term exposure to heat stress), in species with different degrees of specialization to the subterranean environment, would contribute not only to understand the complex biochemical processes underlying heat tolerance in ectotherms, but also their evolution associated with the colonization of these extreme and stable environments.

Exposure time effect
It is well-established that biomarker responses are timedependent (Madeira et al., 2016). After a first defence activation, antioxidant biomarkers may decrease reflecting either exhausted responses or acclimation. We observed here that the cellular defence against thermal stress was not yet onset at 2 days of exposure, as no significant effect of temperature was found for any biomarker. At 7 days of exposure, antioxidant mechanisms were active and not exhausted, but were not effective in allowing thermal acclimation (as mentioned above).
On the other hand, high and low LPO and GSH:GSSG ratio values, respectively, indicated some signs of oxidative stress in both the control and heat treatment at 2 days, if such values are compared with those at 7 days. Oxidative stress in specimens exposed to a temperature that is a priori not stressful (i.e. the control temperature in this study), and its subsequent recovery after longer exposure time, might reflect an initial altered physiological status associated to factors other than temperature, likely as a consequence of the recent collection, transport and handling in the laboratory. This side effect may have important implications from a methodological point of view. When using field-collected specimens for experiments, exposure time should be set so that specimens can recover from any stress associated to collection and transport but also avoiding stress caused by long-term holding at laboratory, non-natural conditions. In the light of our results, controlling for these potentially confounding effects is especially critical for species inhabiting highly stable habitats such as subterranean environments.

Implications in a climate change context
According to our survival experiments, the studied caveadapted species would have a relatively wide TSM to cope with the predicted temperature increase in its current (and unique) locality (a difference of ca. 3.5-5.6 • C between the future predicted temperature at its locality under a pessimistic and optimistic scenario, respectively, and the estimated LT 50 values). However, by combining survival, acclimation and molecular stress biomarkers, we showed that exposure at temperatures below lethal limits induced oxidative stress and alteration of the activity of a key nervous system enzyme (AChE), and prevented acclimation under a subsequent thermal challenge. Long-term persistence of an organism in a given location is more likely to be defined by thermal constraints on physiological performance than thresholds for heat-induced mortality (Sará et al., 2011;Evans et al., 2015). Therefore, in the basis of our results, persistence of P. canyellesi in a scenario of maintained temperatures close to its lethal limit could not be guaranteed. For most taxa, ULTs or other CT max are frequently the only physiological end point available (or the easiest to measure) to estimate TSM in assessments of vulnerability to climate change. For example, Sánchez-Fernández et al. (2016) predicted that species of a lineage of Leptodirni beetles will have suitable habitat conditions (i.e. wide TSM) under climate change scenarios considering experimental data on thermal limits for survival, which was a more optimistic prediction than that obtained with other approaches (bioclimatic models). However, results in our study suggest that such forecasts should be also taken with caution, as they might overestimate the actual sensitivity to thermal stress. Further work exploring the responses of these biomarkers within a range of temperatures between 13 and 20 • C would allow to estimate the exact threshold temperature that onsets cellular and biochemical alterations and get more biologically relevant TSM for the study species. It is also noteworthy that longer exposure times would be needed to detect effects of temperature in longer-term processes, and that measurements on other traits such as fecundity, or on other life-cycle stages, might reduce even more the thermal window for species persistence.
These results stress the need of exploring thermally sensitive processes at different levels of biological organization, beyond typical survival experiments. This approach would allow to obtain more accurate estimates of the capability of poor dispersal species to cope with temperatures outside those they currently experience, and consequently, better estimates of species capacity to face climate change.
Despite the fact that the fragility of the subterranean world is widely recognized (Reboleira et al., 2011;Mammola et al., 2019a), subterranean habitats and species remain largely neglected in conservation programmes. Thus, it is time to go a step forward into the effective protection of this unique and valuable component of biodiversity. The identification and protection of areas with high values of subterranean biodiversity and the application of conservation measures to minimise the threats affecting it must be a priority.

Author contributions
S.P, D.S-F and J.C.S-H conceived the idea and designed the experiments, P.B-G and J.C. collected the specimens, S.P., J.C.S-H and R.C. performed the experiments and S.P. analysed the data and led the manuscript writing. All authors discussed the idea and results, contributed to the manuscript drafts and gave final approval for publication.

Funding
This work was supported by the Agencia Estatal de Investigación (Spain), the Spanish Ministry of Economy and Competitiveness and the European Regional Development