Thermal tolerance breadths among groundwater crustaceans living in a thermally constant environment

In the face of global warming, major progress has been made in the development of species distribution models to predict the influence of temperature change on biodiversity distribution and extinction risks (reviewed in Guisan and Thuiller, 2005). Model predictions are commonly generated by means of correlative approaches linking species' occurrence data with climatic variables (e.g. Raxworthy et al., 2003; Bourg et al., 2005). These correlative models are largely founded on the ecological niche concept (Grinnell, 1917), and more precisely on the realized niche of species in their present-day environment (Hutchinson, 1959). However, as the realized niche depends on both abiotic and biotic interactions as well as historical factors, simple correlative distribution modeling is unlikely to provide realistic predictions if climatic variables do not predominantly shape species distribution patterns (Kearney, 2006). Improving predictions of species distribution requires that models incorporate mechanistic links between an organism's fitness and its environment (e.g. Kearney and Porter, 2009; Chown et al., 2010). The integration of thermal physiology in biogeographical studies is a pertinent tool to evaluate the likelihood of a species to maintain itself in a landscape and its vulnerability to future warming (Calosi et al., 2008; Deutsch et al., 2008).

The thermal tolerance of species is generally considered to play a pivotal role in their geographical distributions (Addo-Bediako et al., 2000; Compton et al., 2007; Calosi et al., 2010). Indeed, the climate variability hypothesis is probably the most widely accepted theory to explain the increase in species range extent, with increasing latitude in the Northern Hemisphere [i.e. Rapoport's rule (Stevens, 1989)]. This hypothesis suggests that the thermal tolerance breadth of species becomes larger with increasing latitude because the annual amplitude of temperature increases towards the poles. Larger seasonal variation in temperature selects for eurythermal species, which may colonize a wider range of habitats in the absence of strong physical or biological barriers to dispersal. Recently, a global analysis of thermal physiology in ectotherms gave support to this hypothesis by showing a relationship between the thermal tolerance breadth of species and the temperature variability they experienced on land and in the ocean (Sunday et al., 2011). Accordingly, ectotherms living in environments with little annual variation of temperature would be expected to exhibit a narrow thermal tolerance range, and thus should be very vulnerable to future warming. For instance, Antarctic marine species living at stable temperatures below 5°C cannot survive a temperature increase of 2°C (Peck et al., 2004). These species are considered to be highly stenothermal (Somero and De Vries, 1967; Peck and Conway, 2000; Pörtner et al., 2007), with thermal tolerance ranges that are two to four times smaller than those for lower latitude species (Peck et al., 2004).

Similarly to Antarctic marine environments, groundwater deeper than 5 m below the soil surface shows almost no seasonal variation of temperature [i.e. annual amplitude of temperature typically less than 2°C (Ginet and Decou, 1977; Freeze and Cherry, 1979)]. Yet groundwater ecosystems are not less sensitive to global warming, because the annual mean temperature of groundwater closely tracks that of air temperature (Freeze and Cherry, 1979). According to the climatic variability hypothesis and the low seasonal variation of temperature in their environment, groundwater-obligate species should exhibit narrow breadths of thermal tolerance because they should maximize their physiological performance over a very narrow thermal range [stenothermal profile (Huey and Kingsolver, 1989)]. Physiological experiments on the thermal physiology of groundwater invertebrates are scarce, and none of them were specifically designed to test the predictions of the climate variability hypothesis (Pattée, 1965; Ginet and Mathieu, 1968; Mathieu, 1983; Issartel et al., 2005a; Issartel et al., 2005b; Issartel et al., 2006; Colson-Proch et al., 2009; Colson-Proch et al., 2010). Several studies examined the physiological responses of groundwater amphipods to cold temperatures [<3°C (Issartel et al., 2005a; Issartel et al., 2006; Colson-Proch et al., 2009)] or heat shocks (Colson-Proch et al., 2010) but did not attempt to determine their thermal tolerance breadths. Issartel et al. (Issartel et al., 2005b) examined the physiological response of two groundwater amphipods, Niphargus virei and Niphargus rhenorhodanensis, over a wide range of temperatures (−2 to 28°C). Yet the stenothermy and/or eurythermy of these groundwater amphipods was defined by comparison of their physiological responses to a third surface-water amphipod (i.e. Gammarus fossarum), without mention of the thermal variability experienced by the three species in their natural habitats.

In this study, we determined the thermal tolerance breadth of three groundwater obligate isopods colonizing habitats characterized by annual temperature amplitude of less than 1°C. The behavioral and metabolic responses of these isopods were examined over a range of temperatures (from 2 to 26°C) by measuring five quantitative variables. Four of them – survival, locomotor activity, oxygen consumption and the concentrations of total free amino acids and sugars – have been commonly used to determine the influence of temperature on the fitness of ectotherms (Huey and Kingsolver, 1989; Angilletta et al., 2002; Pörtner, 2002; Issartel et al., 2005b; Colson-Proch et al., 2009; Dell et al., 2011). The fifth variable, the immune defense, has rarely been used to study the thermal physiology of ectotherms, although it is also tightly linked to organism fitness (Moret and Schmid-Hempel, 2000), and was recently shown to respond quickly to thermal stress in crustaceans (Matozzo et al., 2011).

The three selected species of isopods belonged to the genus Proasellus Dudich 1925 (Asellidae). The first two species, namely Proasellus n. sp. 1 and Proasellus n. sp. 2, belonged to the morphospecies Proasellus cavaticus (Henry, 1976). These two cryptic species were recently identified by Morvan et al. (Morvan et al., 2013) using the molecular species delineation method (Lefébure et al., 2006), and the general mixed Yule–coalescence method (Pons et al., 2006). These two molecular methods also enabled to corroborate the taxonomic validity of the third species, Proasellus valdensis (Chappuis 1948) (Morvan et al., 2013). The three species were collected in caves and springs of the Jura Mountains, France (Table 1). The annual mean temperature of the sites ranged from 6.3 to 11.7°C and the seasonal variation of temperature in each site was less than 1°C. Only one population of Proasellus n. sp. 1 and Proasellus n. sp. 2 (N=140 individuals per population) was sampled because these two species were restricted to a single limestone aquifer (geographic ranges <50 km²; F.M. et al., unpublished data). Two populations of P. valdensis (Pv1 and Pv2) were sampled, because this species occupied a much larger geographic range in the Jura, the Prealps and the Alps [ca. 10,000 km² (Henry, 1976)]. These two populations were used to assess differences in response to temperature variation among populations (Gaston and Spicer, 1998). All individuals were collected delicately from stones using a paintbrush. They were transported to the laboratory within 6 h of capture in aerated and refrigerated plastic containers filled with freshwater from the collection sites. At the laboratory, individuals were kept in aquaria in constant darkness at 10°C (±0.3°C) in a thermo-regulated room for 15 days. Aquaria were previously filled with chemically controlled water [96 mg l⁻¹ NaHCO₃, 39.4 mg l⁻¹ CaSO₄·2H₂O, 60 mg l⁻¹ MgSO₄·7H₂O and 4 mg l⁻¹ KCl, pH 7.5, specific conductance=250–260 μS cm⁻¹ (US Environmental Protection Agency, 1991)] continuously aerated with an air pump to maintain concentrations of dissolved oxygen between 9.5 and 10.5 mg l⁻¹. During these 15 days of acclimation to laboratory conditions, animals were fed ad libitum with soaked alder leaves. We used an acclimation temperature of 10°C, which was very close to the habitat temperature of three out of the four studied populations. The Pv2 population lives at lower temperatures (5.9–6.7°C; Table 1), but previous experiments showed that populations of P. valdensis could be reared at a temperature of 10°C for several months in the laboratory without significant mortality (Colson-Proch, 2009). Therefore, we assumed that 10°C was a comfortable temperature for the four populations tested.

The thermal tolerance breadth of each population was determined by measuring the physiological responses of the populations at eight temperatures (2, 4, 7, 10, 13, 16, 19 and 22°C). For each population and temperature tested, one set of 14 individuals was placed in a 1 liter glass beaker containing soaked alder leaves and filled with 700 ml of chemically controlled aerated water (see composition above). Beakers were placed in thermoregulated water baths (TECO-TC20, precision ±0.5°C, Ravena, Italy), and their water was cooled or heated by conduction through glass walls. Starting from a rearing temperature of 10°C, water temperature in the beakers was increased or decreased by 1°C every 4.5 days. Although recent investigations (Rezende et al., 2011; Santos et al., 2011) proposed the use of fast ramping rates of temperature (1°C min⁻¹ or 1°C h⁻¹) to measure the thermal tolerance of organisms, we applied a slower ramping rate (0.22°C day⁻¹), which was more appropriate and ecologically relevant for groundwater species experiencing low thermal variability in their native habitats (Terblanche et al., 2011). The ramping rate used in the present study was a compromise between the fast (1°C day⁻¹) and the slow (1°C month⁻¹) rate used in previous experiments with Antarctic ectotherms (Peck et al., 2009; Barnes et al., 2010), which also live in thermally buffered environments. Once the desired temperature was reached, the temperature ramp was stopped and animals were maintained at this temperature for the next 5 days for acclimation. During this 5 day period, alder leaves were removed from beakers to ensure that the digestive tract would be empty during subsequent O₂ consumption measurements. At the end of the 5 day acclimation period, live asellids in the beakers were counted (measure of survival) and collected for other physiological measurements (locomotor activity, respiratory activity, immune defense and concentrations of total free amino acids and sugars, see below). As our protocol led to longer experimental times for isopods tested at high and low temperatures (2°C, 22°C) compared with isopods tested at temperatures close to 10°C (7°C, 13°C), we verified that time per se had no influence on survival rates: 14 individuals per population were maintained at 10°C for 59 days (the longest exposure duration) to evaluate their survival rates in the absence of temperature variation.

Locomotor activity was measured in darkness using near infrared (NIR) light. Prior to the measurement, we checked that NIR light did not significantly affect the behavior of the three Proasellus species. Two rectangular aquaria (30 cm long, 10 cm wide) were filled to a height of 3 cm with chemically controlled aerated water and covered at 50% with a screen opaque to NIR light. For each population, 10 individuals per aquarium were introduced at the boundary between the covered and uncovered parts of the aquaria. After 2 h under NIR light or under control conditions (in darkness without NIR light), the distribution of the individuals was recorded. This test was repeated four times for each population and each condition (NIR, control). Results showed that the proportion of individuals recovered below the opaque screen was not significantly different between NIR and control conditions (Student's t-tests, P>0.05).

For each population, live animals collected at the end of the 5 day acclimation period were introduced into glass beakers with a bottom area adapted to video recording (i.d. 6.5 cm). Glass beakers were placed in thermoregulated baths to perform video recordings at the acclimation temperature. A maximum of seven individuals were recorded at the same time under NIR light using an infrared light video camera (SONY, HDR-XR550VE with Nightshot mode, Tokyo, Japan). Videos of 30 min were recorded after 7 h of acclimation of animals to measurement conditions. Videos were performed between 01:00 and 01:30 h to minimize the potential influence of building vibrations on animal activity. Analysis of video recordings was performed with Ctrax 1.5.2 [http://ctrax.sourceforge.net/ (Branson et al., 2009)] and MATLAB R-2008.a (The MathWorks, Natick, MA, USA). We measured the track length of each individual monitored during the 30 min video and calculated its average velocity. Depending on mortality rate, up to 14 average velocity values were obtained per population for each tested temperature.

Oxygen consumption of Proasellus was measured using a micro-respiration system (UNISENSE, Aarhus, Denmark) according to Navel et al. (Navel et al., 2010). After video recording of locomotor activity, each individual was transferred into a respiration chamber (volume: 4.2–4.5 ml) previously filled with chemically controlled aerated water. Chambers contained a stirring system (magnetic stirrer) that prevented any vertical gradient of dissolved oxygen into water. Control chambers without animals were also prepared to determine background microbial respiration in the experimental system. All chambers were immersed into a thermoregulated bath to control temperature. Animals were inserted in respiration chambers 1 h before starting the measurements to reduce the potential influence of animal stress on respiration measurements. Dissolved O₂ was measured using a microsensor that was inserted into each chamber without affecting the transfer of atmospheric O₂ into chambers. Dissolved O₂ was measured every 1.5 h during 9 h. The rate of linear decrease of oxygen concentration with time (from 0 to 9 h) was used to calculate the respiration rate of each individual in μg O₂ consumed h⁻¹. The respiration rates were corrected for microbial respiration measured in control chambers and are reported relative to animal dry mass.

The immune defense of Proasellus was assessed by measuring phenoloxidase activity in haemolymph according to Rigaud and Moret (Rigaud and Moret, 2003). After oxygen respiration, each individual was dried on a paper towel, killed by inserting a fine inox tip into their brain, and pricked in the abdomen to extract haemolymph. Animals were then transferred into pierced 0.5 ml centrifuge tubes that were individually inserted into 1.5 ml pre-weighted centrifuge tubes. The tubes were centrifuged at 800 g during 5 min at 4°C to collect haemolymph. The mass of the 1.5 ml tubes containing haemolymph was determined to evaluate the volume of collected haemolymph for each individual (assuming a haemolymph density of 1 mg μl⁻¹). After centrifugation, haemolymph was used for immune defense analyses. Remaining individuals after haemolymph extraction were frozen, freeze-dried and weighted for oxygen consumption calculations (see above) and analyses of total free amino acids and sugars (see below). For each individual haemolymph extract, activity of naturally activated phenoloxidase enzymes was measured using a spectrophotometric assay (Rigaud and Moret, 2003). Haemolymph was mixed with 9.5 μl of cold phosphate buffered saline (PBS; 8.74 g l⁻¹ NaCl, 1.78 g l⁻¹ Na₂HPO₄·2H₂O, pH 6.5) and centrifuged at 1300 g for 5 min at 4°C. Five microliters of the supernatant (or 5 μl of PBS for the blanks) were deposited in a well of ELISA microplate and mixed with 140 μl of ultrapure water and 20 μl of PBS. Then, 20 μl of l-DOPA (3,4-dihydroxyphenyl-l-alanine, 4 mg ml⁻¹) were added to each well to initiate the reaction between phenoloxidase and l-DOPA, which produces a quinone. The reaction was allowed to proceed at 30°C in a microplate reader (VersaMax ELISA microplate reader, Molecular Devices, Sunnyvale, CA, USA) for 1.5 h. Readings were obtained every 30 s at 490 nm. Enzyme activity was measured as the slope of the reaction curve during the linear phase of the reaction (linear regression with R²>0.8) and reported as the increase in optical density (ΔOD) per minute and per microliter of haemolymph (Vogelweith et al., 2011).

The concentration of total free amino acids and sugars was measured in the four populations using individuals acclimated at 4, 10 and 16°C, except for P. n. sp. 1, for which analyses were only performed at two temperatures (10 and 16°C) due to the low survival of this species at 4°C. For each population, two to three individuals were pooled in each replicate (four to seven replicates per experimental condition, depending on survival and individual masses) to obtain sufficient material for the quantification of free amino acids (glycine, serine, ornithine, leucine, valine, isoleucine, proline, glutamic acid, phenylalanine, putrescine and threonine), free sugars (ribose, glucose, maltose and trehalose) and three other metabolites (ethanolamine, inositol and glycerol). These compounds were extracted from lyophilized organisms according to Laparie et al. (Laparie et al., 2012). Pools of lyophilized animals were weighed and then homogenized in 300 μl of a mix of methanol/chloroform (2:1) using a bead-beating device (Retsch MM301, Retsch, Haan, Germany) at 25 Hz for 1.5 min with two tungsten beads. The samples were transferred at −20°C for 10 min, 200 μl of ice-cold ultrapure water was added, and each sample was vortexed. Following centrifugation at 4000 g for 10 min at 4°C, 300 μl aliquots of the upper aqueous phase, which contained polar metabolites, were transferred to microtubes and vacuum-dried (Speed Vac Concentrator, MiVac, Genevac, Ipswich, UK) at 30°C.

The polar phase aliquots were resuspended in 30 μl of 20 mg ml⁻¹ methoxyaminehydrochloride (Sigma-Aldrich, St Louis, MO, USA) in pyridine prior to incubation under orbital shaking at 40°C for 90 min. Following incubation, 30 μl of N-methyl-N-(trimethylsilyl) trifluoroacetamide (MSTFA; Sigma-Aldrich, #394866) was added, and derivatization was conducted at 40°C for 45 min. We measured the concentrations of free amino acids, free sugars and three other metabolites (indicated above) from the whole crustacean by gas chromatography-mass spectrometry (GC–MS). Our GC–MS system comprised a CTC CombiPal autosampler (GERSTEL, Mülheim an der Ruhr, Germany), a Trace GC Ultra chromatograph and a Trace DSQII quadrupole mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA). One microliter of each sample was injected using the split mode (25:1). The injector temperature was held at 250°C. The oven temperature ranged from 70 to 147°C at 9°C min⁻¹, from 147 to 158°C at 0.5°C min⁻¹, from 158 to 310°C at 5°C min⁻¹, and remained for 4 min at 310°C. A 30 m fused silica column (TR5MS, i.d. 25 mm, 95% dimethyl siloxane, 5% phenyl polysilphenylene–siloxane, Thermo Fisher Scientific, Courtaboeuf, France) was used with helium as the carrier gas at a rate of 1 ml min⁻¹. MS detection was achieved using electron impact. The ion source temperature was set to 250°C, and the MS transfer line temperature was set to 300°C. All samples were run under the SIM mode (electron energy: −70 eV). GC–MS peaks were accurately annotated using both mass spectra (two specific ions), and a retention index specific to each compound. Randomized sample sequences were established for sample injection, and each sequence was initiated with a quality control. Calibration curve samples for pure reference compounds at 5, 10, 20, 50, 100, 200, 500, 750 and 1000 μmol l⁻¹ concentrations were run. Arabinose was used as an internal standard. Chromatograms were deconvoluted using XCalibur 2.0.7 (Thermo Fisher Scientific). Then, metabolite concentrations were quantified using the quadratic calibration curves for each reference compound and are reported relative to animal dry mass.

For each population, survival values, expressed as a percentage, followed a normal distribution. A Gaussian curve was adjusted to survival data (all R² were >80%) in order to determine the experimental temperatures delimiting 95% of survival for each population. The adjustment was performed with Minitab software (Minitab, version 13.32, State College, PA, USA). For each population, the influence of temperature on locomotor activity (mean individual velocity) and immune defense (phenoloxidase activity) was tested by one-way ANOVA. When a significant influence of temperature was detected, we determined the temperature range around 10°C (temperature of rearing at the start of the experiment) for which no significant change in locomotor activity (or immune defense) was measured. For that, data obtained at a temperature of 10°C were compared with data obtained at each temperature using Tukey's post hoc tests. The normality and the homoscedasticity assumptions were verified using the Shapiro–Wilk test and Levene's test, respectively. When necessary, data were ln(x+1)-transformed to fit these assumptions. Statistical analyses were performed using Statistica 7 (Statsoft, Tulsa, OK, USA).

The influence of temperature on concentrations of total free amino acids and sugars was tested for each population using one-way ANOVA except for P. n. sp. 1, for which we used a Student's t-test due to the lack of data obtained at a temperature of 4°C. A between-group principal component analysis (Dolédec and Chessel, 1991) was run to determine the free amino acids and sugars that were specifically affected by temperature change in the four populations. This multivariate analysis aimed to maximize differences in mean concentrations of free amino acids and sugars among groups defined here as the combinations between all populations and temperatures. ADE 4 software (Thioulouse et al., 1997) was used to run the between-group principal component analysis.

No mortality was observed in control animals maintained at 10°C during the course of the experiment. We obtained contrasted survival curves among the three species of Proasellus (Fig. 1). More than 90% of the specimens of P. valdensis survived exposures at temperatures ranging from 2 to 13°C (Pv1), or from 7 to 16°C (Pv2). In contrast, the survival of P. n. sp. 1 and P. n. sp. 2 rapidly declined below 90% when temperature departed slightly from 10°C. Only 67 and 47% of the individuals survived at 7°C for P. n. sp. 1 and P. n. sp. 2, respectively. At 13°C, the survival rates of these two species were 84% (P. n. sp. 1) and 52% (P. n. sp. 2). The thermal range over which 95% of individuals survived exceeded 11°C in P. valdensis whereas it was lower than 4°C for P. n. sp. 1 and P. n. sp. 2 (Table 2).

The locomotor activity of the four Proasellus populations was significantly influenced by temperature (one-way ANOVAs, F_4,41=4.9 and P<0.01 for P. n. sp. 1, F_6,42=6.1 and P<0.001 for P. n. sp. 2, F_7,84=4.4 and P<0.001 for Pv1, F_6,70=4.3 and P<0.01 for Pv2). Yet locomotor activity was differentially altered among populations (Fig. 2). One population of P. valdensis (Pv1) maintained a similar activity at 7°C and 10°C (Tukey's post hoc test, P=0.87), and decreased its locomotor activity at 4 and 13°C compared with that at 10°C (Tukey's post hoc tests, P<0.05 for the two temperatures). For this population, velocity rose up again at 19°C (Fig. 2) when thermal conditions became stressful, as indicated by much lower survival rates (<70%; Fig. 1). The second population of P. valdensis (Pv2) reduced its locomotor activity by more than 40% when temperature diverged from 10°C (Tukey's post hoc test, P<0.05 for comparison between 7 and 10°C, P<0.01 for comparison between 10 and 13°C) and it did not show any increase in activity at highest temperatures. Locomotor activity was the highest at 10°C and 13°C for P. n. sp. 1 and was significantly reduced at 7 (−53%) and 16°C (−55%) in comparison with activity measured at 10°C (Tukey's post hoc tests, P<0.02 for comparison between 7 and 10°C, P<0.04 for comparison between 10 and 16°C). Proasellus n. sp. 2 showed the largest temperature range, with no significant change in locomotor activity from 7 to 16°C when compared with that measured at 10°C (Tukey's post hoc tests, P>0.8).

Increased acclimation temperatures led to a significant rise in oxygen consumption in the four populations (Fig. 3). An exponential relationship between oxygen consumption and temperature was observed for all surviving individuals of the Pv2 population of P. valdensis, P. n. sp. 1 and P. n. sp. 2. In contrast, individuals of the other population of P. valdensis (Pv1), which survived at temperatures of 19 and 22°C, no longer increased their metabolic rates above 16°C. Arrhenius activation energies were determined for all populations using data fitting the exponential curves. Pv1 had the highest activation energy (86.4±18.2 kJ mol⁻¹, mean ± s.d.) whereas Pv2 showed the lowest value (42.4±5.1 kJ mol⁻¹). The two species P. n. sp. 1 and P. n. sp. 2 presented intermediate activation energies (53.0±9.9 kJ mol⁻¹ for P. n. sp. 1 and 63.2±9.5 kJ mol⁻¹ for P. n. sp. 2). Pv1 had a significantly higher activation energy than the other populations (ANCOVAs, P<0.01 for all pairwise comparisons between Pv1 and other populations), which exhibited comparable activation energies among them (ANCOVAs, P>0.4 for pairwise comparisons between Pv2, P. n. sp. 1 and P. n. sp. 2).

Temperature had contrasting influences on the immune response of the four Proasellus populations (Fig. 4). The phenoloxidase activity measured in the haemolymph of P. valdensis was not significantly affected by temperature (one-way ANOVAs, F_7,80=1.6 and P>0.15 for Pv1, F_6,71=1.0 and P>0.44 for Pv2). In contrast, temperature induced an immune response in P. n. sp. 1 and P. n. sp. 2 (one-way ANOVAs, F_4,42=3.9 and P<0.01 for P. n. sp. 1, F_6,42=7.2 and P<0.001 for P. n. sp. 2). For these two species, phenoloxidase activity was the lowest at 10°C and increased by more than 60% when temperature decreased to 7°C (P. n. sp. 1 and P. n. sp. 2) or increased to 13°C (P. n. sp. 2) or 16°C (P. n. sp. 1).

There was no significant effect of temperature on the concentration of total free amino acids in the two populations of P. valdensis (one-way ANOVAs, F_2,15=2.0 and P>0.16 for Pv1, F_2,10=0.8 and P>0.45 for Pv2; Fig. 5). In contrast, their concentration was significantly influenced by temperature in P. n. sp. 1 and P. n. sp. 2 (Student's t-test, t₁₀=3.1 and P<0.02 for P. n. sp. 1, one-way ANOVAs, F_2,10=28.3 and P<0.001 for P. n. sp. 2), and decreased when specimens were acclimated at 4 or 16°C compared with 10°C (Fig. 5). Temperature change also significantly influenced the concentration of total free sugars in the two populations of P. valdensis (one-way ANOVAs, F_2,15=13.9 and P<0.001 for Pv1, F_2,10=5.1 and P<0.05 for Pv2) and in P. n. sp. 2 (one-way ANOVAs, F_2,10=6.5 and P<0.05). However, we found no significant effect of temperature on the concentration of free sugars in P. n. sp. 1 (Student's t-test, t₁₀=1.5, P>0.15), for which no data could be obtained at 4°C due to high mortality (Fig. 5).

The C1 and C2 axes of the between-group principal component analysis explained 80.3% of the inter-group variability contained in the data matrix, with 50.3% on the first axis (Fig. 6). Temperature changes from 10°C to 4 or 16°C increased the scores of the two populations of P. valdensis along axis C2 (Fig. 6A). This was associated with a decrease in the concentrations of glucose, maltose, trehalose and glycerol and an increase in the contents of phenylalanine, glutamic acid and ribose (Fig. 6B). In contrast to P. valdensis, temperature changes from 10°C to 4 or 16°C increased the scores of P. n. sp. 1 and P. n. sp. 2 along axis C1. For these two species, increasing scores along axis C1 reflected a decrease in the concentrations of several free amino acids including glycine, serine, ornithine, leucine valine, isoleucine, proline and threonine.

The present study demonstrates that physiological responses to temperature variation differ greatly among species of subterranean isopods belonging to the Proasellus genus. The thermal breadth delineating a 95% survival was threefold higher for the two populations of P. valdensis (>13.0°C for Pv1 and 11.4°C for Pv2) compared with P. n. sp. 1 and P. n. sp. 2. The climate variability hypothesis (Stevens, 1989) was supported by the low ranges of thermal tolerance observed in P. n. sp. 1 and P. n. sp. 2. As expected for organisms living in thermally buffered environments (Pörtner et al., 2007), P. n. sp. 1 and P. n. sp. 2, which experience stable thermal conditions in their natural habitats, exhibited stenothermy and maximized their performance over a narrow range of temperature (Huey and Kingsolver, 1989). In contrast, P. valdensis exhibited a wide breadth of thermal tolerance that was comparable to that of the subterranean amphipod Niphargus rhenorhodanensis (Issartel et al., 2005b). The ‘conservatism’ of eurythermy in two populations of P. valdensis indicated that the process of local adaptation to thermally buffered habitats (e.g. Alleaume-Benharira et al., 2006; Sanford and Kelly, 2011) did not occur in this subterranean isopod. This absence of local adaptation could result from dispersal and gene flow among populations of P. valdensis living in thermally contrasting habitats. Indeed, Colson-Proch (Colson-Proch, 2009) reported that habitat temperature among 11 populations of P. valdensis living at altitudes ranging from 288 to 1639 m varied from 4.1 to 11.8°C. Although each population lives in a thermally buffered habitat, the recent dispersal and gene flow among populations might have swamped the establishment of locally adaptive alleles (Bridle and Vines, 2007). This hypothesis was supported by a recent phylogeographic study that showed no clear genetic structuration in COI and 16S genes among populations of P. valdensis (F.M., D. Eme, C. Colson-Proch, P. Jean, S. Calvignac, L. Konecny-Dupré, F. Hervant and C.J.D., submitted). However, another plausible hypothesis would be that the three species investigated in this study differed significantly in their effective population size, with P. valdensis having much smaller effective population sizes. Under this scenario, genetic drift would hamper adaptive selection.

This study also showed that all physiological variables were not equally efficient in terms of evaluating the thermal sensitivity of our subterranean organisms. The response of immune defense to temperature was congruent with the survival rates of the two species P. n. sp. 1 and P. n. sp. 2: the phenoloxidase activity increased at temperatures inducing mortality higher than 25% (Fig. 4). Phenoloxidase activity is known to play a key role in recognition of and defense against microbial infections in crustaceans (Söderhäll and Cerenius, 1998). Matozzo et al. (Matozzo et al., 2011) measured significant increases in phenoloxidase activity in the haemolymph of the crab Carcinus aestuarii during exposures to stressful temperatures (4 and 30°C). Similarly, the measurement of phenoloxidase activity in the haemolymph of Proasellus species appeared to be a pertinent method with which to detect their level of thermal stress. The response of total free amino acids to temperature changes in the four populations was congruent with survival rates and immune responses. Data on survival rates and phenoloxidase activity showed that P. n. sp. 1 and P. n. sp. 2 did not perform well at temperatures of 4 and 16°C, for which the amounts of free amino acids and sugars were significantly reduced in comparison with the amounts measured at 10°C in these two species. The large cold- and warm-induced decreases of amino acids and sugars might have resulted from an alteration of the metabolic network and/or from a reduced energetic investment to reproduction (Harshman and Zera, 2006; O'Brien et al., 2002) at these stressful temperatures for both P. n. sp. 1 and P. n. sp. 2. In contrast, the two populations of P. valdensis were less vulnerable to temperature changes, and the amounts of free amino acids remained similar among specimens exposed to 4, 10 and 16°C. A decrease in total free sugars was detected at 4 and 16°C compared with specimens maintained at 10°C for three of the four studied populations. Free sugars represent essential compounds that fuel the energetic metabolism. From 10 to 16°C, the aerobic metabolism was enhanced, increasing the energetic demand at the whole-organism level. The reduction of sugar amounts, and more particularly glucose, suggests that these metabolites were more rapidly catabolized and that the intermediary metabolism was not able to efficiently fulfill the changing physiological needs of the crustaceans. Reduced sugar amounts were also detected at 4°C in comparison with 10°C, which is consistent with the results obtained for glucose and glycerol in the cold-hardy groundwater amphipod N. rhenorhodanensis exposed to 10°C and 3°C (Colson-Proch et al., 2009). These compounds, and more particularly trehalose, which protects biological structures from cold injuries, have been found to be accumulated in cold-acclimated insects (Lee et al., 2002; Koštál et al., 2011), but their present diminutions likely prevent any role for the enhancement of the basal cold tolerance of Proasellus populations. Glycerol, the most commonly known polyol that is accumulated in cold-acclimated insects (Koštál et al., 2007; Lalouette et al., 2007), also decreased when acclimation temperature was cooled down to 4°C. These metabolic changes in the 4°C-exposed crustaceans were indicative of an alteration of the energetic metabolism, as already evidenced in some cold-acclimated insects (Colinet et al., 2012; Foray et al., 2013). The variations of intermediary compounds of the tricarboxylic acid cycle, which represent useful biomarkers of physiological adjustments in cold-exposed arthropods (Michaud et al., 2008; Foray et al., 2013), together with the activity of the ATP-dependent sodium–potassium ion pump, should be examined in further studies to precisely document this cold-induced metabolic alteration.

Locomotor activity measured for the three species did not respond to temperature in a manner similar to survival rate, immune defense and concentrations of total free amino acids and sugars. For instance, the endemic species P. n. sp. 2, which showed little capacity for survival at temperatures different than 10°C, maintained a similar locomotor activity at temperatures ranging from 7 to 16°C. In addition to the technical challenge of accurately measuring the movement of these 5-mm-long arthropods in the dark, and despite numerous studies using locomotor activity as a proxy of isopods' performance (Lagerspetz and Vainio, 2006; Dailey et al., 2009; Schuler et al., 2011), mean velocity did not appear to be a pertinent indicator of the performance in Proasellus species for at least two reasons. First, mean velocities were extremely low (between 1 and 8 cm min⁻¹), probably reflecting the low mobility of subterranean isopods. In comparison, video recordings of the surface aquatic isopod Asellus aquaticus reported locomotor velocities varying between 7 and 48 cm min⁻¹ in relation with temperature (Lagerspetz, 2003). Second, the three species of Proasellus apparently exhibited different behavioral strategies to cope with temperature changes. Indeed, the two populations of P. valdensis and P. n. sp. 1 presented a narrow peak of locomotor activity around 10°C (7–10°C for Pv1, 10°C for Pv2 and 10–13°C for P. n. sp. 1), whereas the locomotor activity of P. n. sp. 2 was constant between 7 and 16°C despite high mortality rates (i.e. ca. 50%). Hervant and Renault (Hervant and Renault, 2002) distinguished between two distinct behavioral strategies among aquatic isopods exposed to constraining conditions. Accordingly, P. valdensis and P. n. sp. 1 reduced their locomotor activities to maximize their trade-off between foraging and survival (‘sit-and-wait’ strategy), whereas P. n. sp. 2 maintained a high locomotor activity to avoid stressful temperatures (‘escape’ strategy). Therefore, mean locomotor velocity was not a pertinent variable to compare the thermal sensitivities of the four populations, and other behavioral measures might have been more indicative of individual performance under contrasted temperatures (e.g. Harris et al., 2011).

Like locomotor activity, the change in oxygen consumption during temperature variation was not linked to thermal tolerance breadth as determined by survival rates, immune defense and concentrations of total free amino acids and sugars. The prediction of a negative relationship between activation energy and thermal tolerance breadth was not validated as the most tolerant population (Pv1) according to survival rates exhibited the highest activation energy. Moreover, there was a twofold difference in activation energy between the two populations of P. valdensis. Nevertheless, the exponential increase in oxygen consumption with temperature observed in Proasellus species was coherent with the Boltzmann–Arrhenius relationship from chemical reaction kinetics (Wittmann et al., 2008; Dell et al., 2011; Huey and Kingsolver, 2011). This exponential relationship suggested that the survival of Proasellus species at low and high temperatures depended on the ability of these isopods to satisfy their aerobic metabolic demand (Pörtner, 2010). The unique exception to this relationship was represented by population Pv1, the aerobic metabolic rate of which fell below the exponential curve of activation energy at high temperature (19 and 22°C). However, one should note that a temperature of 22°C was lethal to all individuals except those of this population. The physiological model developed by Pörtner (Pörtner, 2010) would attribute this departure from a Boltzmann–Arrhenius relationship to a shift from aerobic to anaerobic metabolism, when energy demand can no longer be satisfied by environmental oxygen availability. This temperature-mediated shift in metabolism could be checked by analyzing major end-products of anaerobic metabolism (e.g. lactate, succinate) in individuals of Pv1 exposed to high temperatures.

In summary, differences in physiological traits among groundwater species was best evaluated using survival as a surrogate for measuring their thermal tolerance breadth, and immune defense and concentrations of free amino acids and sugars to detect sub-lethal effects of temperature variation.

Based on these three variables, the main result of this study was to demonstrate that populations of three groundwater species could exhibit distinct responses to temperature variation, although they all colonized habitats showing little seasonal temperature fluctuations. This indicates that the climate variability hypothesis is unlikely to explain differences in thermal niche breadth among groundwater species, unless it is expanded to account for by the effect of the two interacting processes controlling local adaptation, namely, the strength of divergent selection and gene flow among populations (Räsänen and Hendry, 2008). This expanded view of the climate variability hypothesis necessarily leads to the need to consider for each population of a given species the seasonal amplitude of temperature in its natural habitats, the range of temperatures experienced by populations in the species geographic range, and the significance of dispersal among populations. First, the strength of divergent selection among populations of a groundwater species may not necessarily be strong if they colonize shallow subterranean habitats characterized by strong seasonal variation of temperature (Culver and Pipan, 2011). In this case, populations would be expected to exhibit a large thermal niche breadth. Second, most populations of a widely distributed species colonizing deep groundwater habitats would typically encounter little seasonal variation of temperature, although they may live under different mean annual temperatures. Although this situation should favor divergent selection, thermal niche narrowing among populations may be counteracted by maladaptive gene flow among populations (Räsänen and Hendry, 2008; Sexton et al., 2009; Geber, 2011; Hardie and Hutchings, 2010). Finally, populations of narrowly distributed species colonizing deep groundwater habitats would be expected to have a narrow thermal niche breadth because they experience little seasonal or spatial variation of temperature. We explored the likelihood of these predictions using available data for five groundwater species (Table 3). Thermal niche breadth of populations did not significantly relate to the seasonal amplitude of temperature in the habitats (Pearson correlation, R²=0.24, P=0.31). However, thermal niche breadth increased significantly with increasing temperature heterogeneity when the latter was measured as the maximum value between seasonal amplitude and spatial heterogeneity of temperature over the species' geographic range (Pearson correlation, R²=0.95, P<0.001, thermal niche breadth=1.18×temperature heterogeneity + 0.22). Clearly, additional data from multiple species are critically needed to validate this relationship, but it provides both an expanded perspective of the climate variability hypothesis and a working methodology for explaining differences in thermal niche breadth among groundwater organisms.

The present study also has implications for assessing the risk of species extinction due to global warming. Although our results derive from an experimental rise of temperature that is considerably higher than the present-day increase in air temperature, they emphasize the double sensitivity of narrowly distributed species colonizing deep groundwater habitats. Indeed, populations of P. n. sp. 1 and P. n. sp. 2 not only have little physiological capacity to respond to fast temperature changes, but they also benefit from a reduced number of thermally suitable habitats that could act as thermal refuges during climate change. Thus, the extinction risk of these groundwater endemics is critically higher than that of widely distributed species, the populations of which exhibit wider thermal niche breadth and dispose of many potential thermal refuges. Yet a major difficulty in determining the fate of groundwater endemics in response to climate warming stands in our uncertainty about the temporal scale at which key physiological traits evolve. Nevertheless, as pointed out by several authors (e.g. Gaston et al., 2009; Wiens et al., 2009; Chown et al., 2010; Kelly et al., 2012), we greatly encourage the use of climatic models that integrate differences in thermal physiology among species for evaluating extinction risks under various scenarios of global warming.

We are very grateful to Yannick Moret (UMR CNRS 6282 BioGéoSciences, Dijon) for his frequent advice concerning the analysis of phenoloxidase activity in our tiny organisms and Vanessa Larvor (UMR CNRS 6553, ECOBIO, Rennes) for her helpful assistance in metabolite analyses. We also thank the staff of ‘Les Grottes de la Balme’ for the authorization to collect specimens of the population Pv1. The manuscript greatly benefited from insightful comments raised by two anonymous reviewers.