Automated measurement of upper thermal limits in small aquatic animals

Physiological limits are emerging as a key predictor of organismal sensitivity to environmental change (Huey et al., 2012). Physiological limits can be exemplified by the concept of the thermal performance curve, which typically shows how a physiological variable, but also a life history trait or population vital rate, varies in relation to temperature (Sinclair et al., 2016). These performance curves, and their spatial and temporal variation, are useful for understanding how variation in a given trait can influence the sensitivity of a genotype, population or species to temperature (e.g. Chown et al., 2010). Consequently, in the face of global environmental change, renewed attention has been given to the thermal sensitivity of ‘performance’, and in particular, measurements of the upper boundary on the thermal performance curve (referred to hereafter as the upper thermal limit or UTL), as well as the application of these physiological measures to ecologically relevant questions. For example, UTLs have been used to investigate phenotypic and evolutionary responses to temperature change (Brans et al., 2017; Geerts et al., 2015; Manenti et al., 2014; Phillips et al., 2016; Sandblom et al., 2016; Yampolsky et al., 2014), to understand how species distributions are shaped by the environment (Clusella-Trullas et al., 2011; Kellermann et al., 2012; Overgaard et al., 2014), and to understand how niches and distributions might respond to environmental change (Kearney et al., 2009; Nowakowski et al., 2018). From a methodological perspective, UTLs can be measured in a multitude of ways (summarised by Angilletta, 2009; Terblanche et al., 2011). However, with the exception of direct estimates of survival following thermal stress (e.g. LT₅₀), measurements of UTLs typically involve recording the time to loss of locomotory function following exposure to lethal temperature (e.g. knock-down time) or the temperature at which locomotory function is lost (e.g. knock-down temperature or CT_max). These measurements are typically performed by manual observation, with the observer noting either the time or the temperature at which the experimental animals lose locomotory function. It has previously been recognised that such measures of UTLs are subjective and thus prone to error and bias (Lighton and Turner, 2004), both within and among observers. As such, alternative methods that overcome these issues have been developed (e.g. thermolimit respirometry; Lighton and Turner, 2004), but they are data intensive, requiring simultaneous measurement of individual oxygen consumption and activity data. However, in developing their method, Lighton and Turner (2004) concluded that UTLs can be likely be assessed from individual movement data alone. Here, we describe an automated, high-throughput system for measuring UTLs from movement data. Our system utilises video-based tracking software that simultaneously records the movement of multiple individuals when exposed to high, lethal temperature. We implement an algorithm in R that uses the tracking data to find optimal and objective criteria for detecting the loss of locomotory function at high temperature. To illustrate the utility of our method, we apply these criteria to an independent experimental dataset.

We developed our method using the clonally reproducing zooplankter Daphnia magna Straus 1820 as a study organism. We measured the UTL in this species as the ability to maintain mobility at high temperature, recorded as the ‘time to immobilisation’ (T_imm) at a temperature 37°C, which is lethal after several hours of exposure (Yampolsky et al., 2014). Motility is an established parameter for measuring physiological performance in D. magna (Knie, 1982; Yampolsky et al., 2014; Zeis et al., 2004), and given its small body size (approximate body length at maturity is ∼2 mm), our measurement temperature likely reflected the internal body temperature of the organism. Our method can be summarised in four steps. First, we experimentally induced among-individual variation in T_imm through acclimation to different temperatures (e.g. Cambronero et al., 2017; Moyano et al., 2017). Second, when exposed to 37°C, we filmed the experimental individuals, up to and including the point when they lost swimming function (i.e. became immobile). Third, using animal tracking software, we obtained a time-series of movement data for these individuals. Fourth, we developed a custom R script (https://www.r-project.org/) in conjunction with a quantitative statistical approach, to objectively identify a threshold movement value that could define a state of immobility from the tracking data. Each of these steps is described in further detail below.

Developmentally synchronised clonal animals were obtained by isolating neonates (F₀) produced within a 24 h period from mass cultures of the experimental clone (EF 64, hatched from an ephippium collected at Værøy Island, northern Norway, 67.687°N, 12.672°E) that we maintain in our laboratory at 17°C. F₀ individuals and two subsequent generations (F₁ and F₂), both founded from second-brood offspring, were reared at either 17°C or 22°C in 2.5 l aquaria in temperature-controlled climate cabinets (IPP260, Memmert, Schwabach, Germany), a long-day photoperiod (16 h:8 h light:dark) and initial density of 35 females per aquarium. Culture medium [artificial Daphnia medium (ADaM)] was replaced in the aquaria one to two times per week and commercial shellfish diet (1800, Reed Mariculture, Campbell, CA, USA) was provided at ad libitum levels three times per week. We predicted that individuals from the 22°C acclimation treatment would have a higher T_imm than individuals from the 17°C treatment (Cambronero et al., 2017; Yampolsky et al., 2014).

We exposed individual daphnids from each of the acclimation treatments to lethally high temperature using a custom-built aluminium and glass thermostatic well plate (Fig. 1). Forty-five individual glass wells, open on their upper surface (well diameter 16 mm, depth 21 mm), were inserted in a rectangular 5×9 array on an aluminium plate (length 265 mm, width 125 mm, thickness 3 mm). This plate was fitted (and sealed via a series of screws) on top of a rectangular aluminium frame (depth 25 mm, length 265 mm, width 125 mm, thickness 20 mm). A sheet of glass (thickness 3 mm) of the same dimensions as the aluminium frame was glued to its underside. Water, warmed to 37°C (in a 15.0 l capacity water bath, Grant Instruments, Shepreth, UK), was pumped into the water-jacket via five inlet points (Eheim Compact 600 pump, Deizisau, Germany). The wells were positioned in the plate so that there was a space of approximately 3 mm between the bottom of each well and the floor of the well plate, and a space of approximately 4 mm between adjacent wells and the outer walls of the well plate. This means that water could pass freely from the inlet points, both around and beneath each well, to five outlet points located at the opposite end of the plate. Water then flowed back to the water bath. Pilot measurements confirmed that there was no difference in temperature of medium contained in wells at the periphery or centre of the well plate, nor between wells located immediately adjacent to the inlet or outlet points. The temperature of medium in each well also reflected the set temperature of the water bath.

F₂ animals from each of the acclimation treatments were measured approximately 1, 2, 3, 6, 8, 10 and 13 days after reaching maturity (determined as presence of first clutch eggs in brood pouch). On a given recording day, we filmed n=45 individuals for a unique combination of acclimation treatment and age group. Thus, the 17°C and 22°C acclimation treatments were measured in separate ‘blocks’ and the respective age groups within each block were recorded as sequence of ‘runs’ beginning with the 1 day group on the first day and concluding with the final age group 13 days later. On a given recording day, up to 45 min before filming began, F₂ females from two randomly chosen aquaria of the same acclimation treatment were pooled, mixed and placed in fresh ADaM (without food) of the same temperature as the treatment group from which they originated. Forty-five of these females were pipetted individually into plastic vials along with 3.3 ml of ADaM. These vials were then returned to their climate cabinet until filming began. The 45 individuals from an acclimation group (and medium of the same temperature) were transferred into the wells (i.e. one individual per well), which had been pre-heated to the test temperature by circulating medium from the water bath, noting down the well number and time (in seconds) elapsed from the moment that the first individual was placed in a well (on average, it took 3–4 min to introduce all 45 individuals to the well plate). Video recording commenced after the last individual was introduced to a well. The well plate was backlit with an LED light board (Huion A4 LED light pad, set to maximum intensity; Shenzen, China), to provide contrast between the individual in each well and the background. The well plates were filmed from above with digital video cameras (Basler aCA1300-60gm, fitted with 5–50 mm, F1.4, CS mount lenses; Ahrensburg, Germany) at a working distance of 116 cm (from lens to upper surface of well plate), a frame rate of 30 frames s⁻¹ and a resolution of 1280×960 pixels. Video recording ceased when visual inspection indicated that all individuals were motionless (inspection performed ca. every 2 min). As the animals began to lose swimming ability during video recording, we did not physically stimulate them because of possible interference with subsequent tracking. Thus, the thermal limit identified by our method may differ slightly from one defined by the loss of reaction to mild physical stimulus at the chosen measurement temperature (e.g. gentle prodding). After filming, all animals were photographed under a stereomicroscope, and body size (carapace length, mm) was measured with ImageJ software (version 1.49v, National Institutes of Health, Bethesda, MD, USA). In total, we filmed 315 individuals from each of the acclimation treatments (n=630 total).

The video files produced in step 2 were processed in Ethovision (version XT 11.5, Noldus Information Technology, Wageningen, The Netherlands; settings: greyscale pixel range 0–120, pixel size range 2–130, sample rate 3 observations s⁻¹), to produce a time series of velocity data (in mm s⁻¹, travelled by the center-point of each individual).

We developed an R script based upon a moving median (med.filter function from robfilter package; https://cran.r-project.org/web/packages/robfilter/) that calculates the time taken for an individual's swimming velocity to drop below a specified threshold value (V_thresh). Moving medians are frequently used to detect underlying ‘signals’ or patterns from ‘noisy’ time-series data, such as the tracking data utilised here. The moving median requires the parameter w, which specifies the window width, the number of observations in a time series over which each successive iteration of the moving median is calculated. Given that an individual's T_imm should be positively affected by warm acclimation temperatures (Cambronero et al., 2017; Moyano et al., 2017) and should be lower in larger, older individuals (Brans et al., 2017; Messmer et al., 2017), the optimal choice of values for V_thresh and w should be the combination that maximises the amount of variation in T_imm that can be described by acclimation treatment (17 and 22°C), when also statistically controlling for the variation in body size that was present both among and within each acclimation treatment (see Fig. S1).

To find these optimal values, we used the velocity time-series data for 629 of the individuals produced in step 3 (one animal was misplaced during removal from the well plate) to calculate T_imm for 5400 unique combinations of V_thresh (0.01 to 1 in increments of 0.01 mm s⁻¹) and w (20 to 1080 in increments of 20 observations, corresponding to window widths that ranged from 6.67 to 360 s given that the movement of each individual was tracked at a rate of 3 observations s⁻¹). The upper limits of both V_thresh and w were chosen based upon preliminary tests that revealed that values larger or smaller than these produced estimates of T_imm of which little variation could be attributed to acclimation treatment. For each combination of V_thresh and w, the estimated T_imm of each individual was modelled statistically as a function of acclimation treatment and its body size (general linear model, temperature treatment and body size as additive categorical and continuous predictors, respectively). This procedure was also repeated for the interactive version of the same model.

To test the optimal values for V_thresh and w, we applied them to an independent set of time-series swimming velocity data. The same D. magna genotype used to define optimal values for V_thresh and w was reared at a wider range of acclimation temperatures (12, 17, 20, 24 and 28°C) and nine individuals from each treatment group were tested for T_imm using the same procedure described in step 2, but in a single well-plate ‘run’. Developmentally synchronised clonal animals were obtained by isolating neonates produced within a 24 h period from mass cultures of the same genotype described earlier in the Materials and Methods (EF 64). These individuals were reared in 250 ml jars at the five different temperatures (nine replicate jars per temperature) in temperature-controlled climate cabinets (Memmert IP260) at a long-day photoperiod (16 h:8 h light:dark) and a density of five individuals per jar. A commercial shellfish diet (1800, Reed Mariculture) was provided at ad libitum levels three times per week. Culture medium (ADaM; Klüttgen et al., 1994) in the jars was changed one to two times per week. The establishment of these cultures was staggered, with the lower temperatures established earlier so that all individuals were measured for T_imm at a similar developmental stage (i.e. when all were carrying first clutch offspring in the brood pouch). Up to 45 min before video recording in the well plate commenced, nine mature females from each acclimation treatment (one per jar) were placed individually into plastic vials along with 3.5 ml of fresh ADaM (without food) of the same temperature as the treatment group from which they originated (i.e. 12°C to 28°C). These vials were then returned to their respective climate cabinets until recording began. Again we anticipated a positive relationship between T_imm and acclimation temperature. Moreover, our hypothesis was that the V_thresh and w values that maximised the amount of variation in T_imm that could be explained by acclimation temperature (and body size) in the design employed in Step 2 should be as, if not more, effective in a situation in which estimates of T_imm should be free of variation attributable to ‘run’ and ‘block’ effects.

For each of the 5400 combinations of V_thresh and w that were tested, the amount of variation in T_imm that could be explained by acclimation treatment and body size was maximised for V_thresh=0.03 mm s⁻¹ and w=280 observations (Fig. 2A). When this procedure was performed for the interactive version of the same general linear model (i.e. including body size×acclimation treatment interaction), it yielded the same result. Fig. 2A also reveals how the choice of w had greater influence on the estimation of T_imm than V_thresh. Fig. 2B demonstrates the effect of different w values on the estimation of T_imm: relatively wide w values were less sensitive to the variation inherent in a given time series and vice versa for relatively narrow w value. Consequently, for a given V_thresh, a relatively wide w was more likely to produce estimates of T_imm that were higher than those produced by the optimal w and vice versa for a relatively narrow w (Fig. 2B). The optimal values for V_thresh and w produced T_imm estimates that matched our predictions with respect to acclimation temperature and body size: T_imm was higher in individuals that had been acclimated to 22°C and decreased with body size/age (parameter estimates±s.e. and corresponding t-values for 22°C acclimated individuals compared with 17°C acclimated individuals: 566.58±19.23, t=29.50, P<0.0001; body size: −114.97±32.27, t=−3.56, P<0.001; model R²=0.59; Fig. 3A). Moreover, when applying the optimal combination of parameter values to the independent tracking data obtained from five different acclimation temperatures (12–28°C, n=9 individuals per acclimation temperature), T_imm increased consistently with acclimation temperature, with the latter explaining 88.7% of the total variation in T_imm (Fig. 3B).

We view the method and hardware described here as a general framework that end-users can adapt to objectively measure UTLs in their own study organism(s), but caution that the optimal parameters (V_thresh and w) for identifying UTLs described here are unlikely to be generalizable across study organisms/taxa. Our R script and hardware was developed to measure UTLs in (a) a relatively small, (b) aquatic organism at a (c) static test temperature. However, provided that the focal organism can be filmed and subsequently tracked accurately whilst exposed to warming, our system could be adapted to measure larger or smaller study organisms (both aquatic and terrestrial) by scaling and adjusting the hardware accordingly. Furthermore, via the addition of ‘real-time’ temperature data, our system could also incorporate the ramped changes in test temperature that are employed in other measurements of UTLs, such as CT_max. Nevertheless, our method is unlikely to be applicable to all species. A prerequisite for the successful translation of our method is that the study organism/life stage in question is motile and maintains a certain level of spontaneous activity. The physiological UTL may also be difficult to establish based on our method if the study organism becomes immobile as part of its natural response to an increase in ambient temperature. However, having said that, estimating and explaining variation in the ‘voluntary’ immobility response to thermal heating may be equally interesting owing to the consequences this can have on ecological processes such as feeding rates and predation risk.

In summary, we present an alternative means to eliminate observer subjectivity when estimating UTLs. In comparison with other established techniques such as thermolimit respirometry (Lighton and Turner, 2004), our method can more efficiently facilitate the measurement of larger sample sizes, meaning that weaker biological effects can be identified within shorter periods of time. We acknowledge that a possible criticism of our technique is that it employs proprietary tracking software. We chose the tracking software described here based on its long history in animal behaviour research and its relative ease of use. There are, however, several open-source alternatives that are more coding intensive but perform similar tracking functions (summarised in Robie et al., 2017). Provided these open-source programmes are able to track the focal organisms accurately, we see no reason as to why our algorithm would fail to produce transferrable results. Moreover, driven by rapid developments in the application of computer vision methods to address ecological questions (Weinstein, 2017) and, in particular, automated, quantitative analysis of behaviour, we feel that this number of computer vision methods will undoubtedly increase in the near future, making our method more accessible to a broader range of end-users.

We wish to thank V. Iashchenko for helpful discussion in the design of the well plates, and S. Bete and A. Siversen for their construction.