Intraspecific variation in physiological performance of a benthic elasmobranch challenged by ocean acidification and warming

Temperature can vary considerably across the geographic range of a temperate species (Angilletta, 2009). This thermal gradient often favors the evolution of locally adapted populations that exhibit different ecological, morphological and physiological responses to changes in temperature. The accelerating changes in ocean temperature and chemistry (i.e. acidification) experienced by marine ecosystems and the pressing need to better predict the consequent shifts in distribution and resilience of species underscore the urgency to identify the response of organisms challenged by climate change-related stressors at a finer geographic scale. In fact, as carbon dioxide (CO₂) and temperature increase in the oceans, the likelihood of individual species of marine ectotherms surviving local extinction may depend on the level of intraspecific variation in physiological responses and the capacity for acclimatization and rapid adaptation to these major environmental stressors (Baumann and Conover, 2011; Di Santo, 2015; Donelson et al., 2011, 2012; Pörtner and Farrell, 2008; Rummer et al., 2014; Somero, 2010). Temperature is an important ecological factor, known to have a profound effect on many metabolic processes in aquatic ectotherms, such as fishes (Di Santo and Bennett, 2011a; Fry, 1971; Magnuson et al., 1979). At the same time, although some studies report little to no effect (or even a positive effect) of CO₂ on fishes (see, for example, Green and Jutfelt, 2014; Heinrich et al., 2014; Jutfelt and Hedgärde, 2015; Rummer et al., 2013), several studies have presented data suggesting that increasing ocean acidification has the potential to exert several adverse effects on fish life history traits, such as malformation during skeletogenesis (Chambers et al., 2013), reduced survival (Baumann and Conover, 2011), reduced body condition and increased developmental time (Di Santo, 2015), as well as on several other important behavioral and physiological traits, such as alertness, predator avoidance and hunting (Ferrari et al., 2012a,b; Hamilton et al., 2014; Jutfelt et al., 2013; Munday et al., 2009; Näslund et al., 2015).

Predator evasion is considered a key factor in fish survival, especially during early life stages; it is a fundamental process that drives community dynamics (Dahlgren et al., 2006), and individuals with effective anti-predatory responses are able to survive, grow and reproduce. In elasmobranch fishes, the effect of warming and acidification on predator escape performance, especially for juveniles, may pose additional challenges because they grow slowly, thus increasing and prolonging the chance of being predated upon before reaching sexual maturity (DiBattista et al., 2007). Although the consequences of rapid ocean acidification and warming on escape performance could be high, no studies to date have evaluated their combined effect on an elasmobranch fish. To address this issue and improve predictions of species responses to environmental challenges, it is necessary to implement multistressor studies that quantify responses of individuals to simultaneous warming and acidification (Todgham and Stillman, 2013). Temperature is a crucial modulator of fish activities and ocean acidification could exacerbate the potential negative effect of warming on performance. Moreover, as individuals are likely to respond to abiotic stressors depending on their environmental history, it is increasingly being recognized that it is important to test individuals from different localities raised in common garden conditions to account for local adaptation of metabolic functions (Angilletta et al., 2004; Angilletta, 2001; Baumann and Conover, 2011; Di Santo, 2015; Lombardi-Carlson et al., 2003; Norin and Malte, 2012).

The little skate Leucoraja erinacea (Mitchill 1825) is a small benthic oviparous elasmobranch, reproducing year round, and inhabiting near-shore waters along the northwestern Atlantic from Cape Hatteras to the Gulf of Maine (GoM). Frisk and Miller (2006, 2009) found that growth of L. erinacea differed along its geographic range, with size-at-age increasing with latitude. In particular, L. erinacea from the GoM show significantly larger body size than conspecifics in Georges Bank (GB) and the Mid-Atlantic, perhaps suggesting strong philopatry and low population exchange in the fish, although genetic data are still needed to designate population-level differences. However, this anecdotal evidence has been supported by empirical data showing that skates from the northern locality, the GoM, exhibited larger size at hatching when compared with individuals from the neighboring and more southern GB when reared in the same conditions, even though the trade-off was a reduced body condition in the larger skates (Di Santo, 2015). Interestingly, although inhabiting geographically contiguous areas, the GoM and GB L. erinacea show low mixing, which might have favored local thermal adaptation of metabolism and growth (Di Santo, 2015). In fact, although high-resolution data on pH and temperature at a fine geographic scale are scarce, especially for deeper waters typically inhabited by skates, the GB and the GoM present some potentially crucial abiotic differences. The GoM is much deeper, exceeding 200 m in depth, while the GB is shallower, with frequent upwelling (Signorini et al., 2013). The GB may therefore experience more frequent short-term variation in P_CO2, salinity and temperature, but higher resolution data from the two locations are needed to confirm this (Pershing et al., 2001; Rossby, 1996; Signorini et al., 2013). As perturbation of metabolic rates may have consequences for fitness and therefore could affect resilience and survival of organisms, differences in whole-organism physiological responses along their geographic range may establish individual vulnerability to ocean warming and acidification (Barnes et al., 2010; Di Santo, 2015; Somero, 2010; Yamahira and Conover, 2002), thereby providing crucial information for implementation of effective conservation plans.

The maximum and routine metabolic rates (MMR and RMR, respectively) are key measures of performance that relate to Darwinian fitness (Bennett and Huey, 1990). In fact, maximal and routine performance may influence reproductive and competitive success, and the capacity to procure food and escape predators (Bennett and Huey, 1990; Pough, 1989). Although various methods are used to induce MMR, experiments that use chasing protocols have been found to elicit exhaustion and higher metabolic rates than classic critical swimming protocols in benthic fishes (Cutts et al., 2002; Roche et al., 2013; Svendsen et al., 2012). Moreover, manually chasing fish to exhaustion is regarded as a commonly used and appropriate method (Cutts et al., 2002; Roche et al., 2013; Svendsen et al., 2012) and approximates the energy needed to escape an intense chase by a predator. Here, I quantified the effects of ocean warming, acidification and location of origin on the performance of juvenile L. erinacea developmentally acclimatized to different pH and temperature conditions. To distinguish between the single and combined effect of each factor, I employed a fully crossed experimental design that compared (i) MMR and RMR, (ii) absolute aerobic scope (AAS, the metabolic capacity to sustain the escape response), (iii) intensity of exercise (turns min⁻¹) and (iv) recovery time after exhaustion in juvenile skates from the GoM and GB. I hypothesized that, similar to responses previously observed in embryos, juvenile L. erinacea from two locations would exhibit local adaptation in metabolic performance when developmentally acclimatized to different temperatures and acidification conditions.

The research was conducted under the approved Institutional Animal Care and Use protocol no. 11-041 at Boston University. Leucoraja erinacea (total N=84, N=42 per location, N=7 replicates per treatment per location, 6 treatments) were obtained from different females (i.e. one embryo per female, N=84 females) caught at two localities, the GoM (43°N, 68°W; ∼12°C) and GB (41.21°N, 67.38°W; ∼15°C) and maintained at 15–16°C. Newly laid embryos (<1 week old) were transported in a constant-temperature tank, and were randomly assigned to a treatment group (1 embryo per replicate tank, 7 replicate tanks per treatment), and reared in common garden conditions in a temperature-controlled cold room (Harris Environmental Systems; 12°C) at Boston University throughout their embryonic development (∼5 months) and after hatching (∼3 months).

Temperature treatments (15, 18, 20°C) were chosen to test the performance at present and projected temperatures by the end of the century (Di Santo, 2015; IPCC, 2013; Palm et al., 2011). Each tank (150 l) was supplied with clean artificial seawater (Instant Ocean and deionized water) and had completely independent filtration, temperature and CO₂ control (i.e. closed system, 84 independent tanks). Constant temperature (15, 18 or 20°C) was maintained in each tank by a submersible titanium heater unit (Finnex 300 W) controlled by a digital thermostat (Aqua Logic Inc.), and CO₂ levels (with corresponding pH of 8.1 or 7.7) to simulate current and projected (year 2100) conditions according to the climatic model RCP 8.5 (IPCC, 2013; Riebesell et al., 2010). To maintain the appropriate pH, each tank was independently supplied with a mix of air:CO₂ (water pH 7.7 which resulted in P_CO2 ∼1100 μatm) or ambient air (water pH 8.1 which resulted in P_CO2 ∼400 μatm) and controlled by an independent Aqua Medic pH computer. Each CO₂ system had an electronic solenoid that triggered the delivery of CO₂ into a diffuser within the tank, if the pH increased above the set value. Each tank was covered with a clear Plexiglas^® sheet to reduce fluctuations in CO₂. Total alkalinity was estimated using titration and reference material (Andrew Dickson, Scripps Institute of Oceanography, personal communication). Water parameters were calculated in CO2SYS (Pierrot et al., 2006) using suggested constants (Dickson and Millero, 1987) (Table 1). Temperature and pH levels were monitored every 12 h using a Traceable NIST calibrated thermometer and a portable Hach pH meter. Each pH probe was calibrated weekly. Water quality analyses (ammonia, nitrites, nitrates, salinity, dissolved oxygen) were performed and 30% water changes were made weekly to ensure excellent water quality (no ammonia and nitrites), and constant salinity (33 ppt) and dissolved oxygen (>90%). Skates were reared on a 14 h:10 h light:dark photoperiod and, once hatched, were fed daily a diet of frozen mysis shrimp ad libitum but were fasted for 24 h prior to experiments.

All experimental values are reported as means±s.e.m. Two sets of experiments (N=4 and N=3) were conducted about 6 months apart. To ensure there was not an effect of time on the results, the two data sets were compared using an ANOVA but no significant difference was observed (all P>0.05) so data were pooled together. The effects of temperature, pH and location of origin on metabolic rates were explored using a 3-way ANOVA after metabolic rates were mass-adjusted with the scaling coefficient used for elasmobranchs, followed by a Tukey–Kramer multiple comparison test (MCT) to identify statistical differences between treatment group means. The effects of temperature, pH, location and mass on exercise metrics were explored using a 3-way ANCOVA, with mass as the co-variate, followed by a Tukey–Kramer MCT to identify statistical differences between treatment group means. If interactions between factors were detected, they were reported following the analysis. Metabolic rates following fatigue were compared with RMR (controls) using a repeated measures ANOVA followed by a Dunnett's test. All statistical comparisons were based on α=0.05. All analyses were performed in JMP Pro version 11 (SAS Institute Inc.).

The wet mass of juvenile (3 month old) L. erinacea from the GoM was significantly larger than that of same-age GB individuals, regardless of temperature and pH (24.61±0.38 and 12.34±0.16 g, respectively; 3-way ANOVA: F_7,76=136.60, P<0.0001; Fig. 1). Overall, temperature and pH had no significant effect on the mass of juvenile skates (P>0.05), with the exception that low pH exacerbated the effect of high temperature (20°C) on the mass of juvenile L. erinacea from the GoM, resulting in lower body mass (low pH: 21.54±1.04 g, control pH: 26.11±0.39 g; Dunnett's test, P=0.04). The validity of the mass exponent 0.67 was tested using an ANCOVA on metabolic rates and resulted in similar statistical outcomes to the 3-way ANOVA on mass-adjusted rates, so results using the mass exponent 0.67 are reported instead.

MMR following exhaustion differed between the two locations (3-way ANOVA, F_7,76=31.24, P<0.0001; Fig. 2A), with a significant effect of pH (P=0.0004) but no significant effect of temperature (P=0.07). There was a significant interaction of the three treatments (P=0.0005). Routine oxygen consumption rates were significantly affected by temperature, acidification and location (3-way ANOVA, F_7,76=19.04, P<0.0001; Fig. 2A). In particular, RMR significantly increased with temperature (P<0.0001) and pH (P=0.002), and was higher in the GB skates (P<0.0001). Furthermore, there was a significant interaction between location and temperature (P=0.008). AAS differed between individuals from different locations under different conditions (3-way ANOVA; F_7,76=18.01, P<0.0001; Fig. 2B) and decreased with temperature (P=0.0004). However, only the AAS of the GoM skates was significantly reduced by acidification and warming (P=0.002), and pH levels alone did not have a statistically significant effect on AAS (P>0.05). AAS had a significant effect on time to fatigue (P=0.0001) and intensity of exercise (P<0.0001), but had no effect on overall number of turns during exercise (P=0.3). Thermal sensitivity quotients (Q₁₀) in skates from both localities were high for RMR (Q₁₀>5), while they were relatively insensitive for MMR (Q₁₀≈1, with the exception of GB skates at control pH conditions, which had a typical Q₁₀=2.83; Table 2). Consequently, Q₁₀ values for AAS showed a negative effect of temperature on rates (all Q₁₀<0.5, with the exception of GB skates at control pH conditions, which had a typical Q₁₀=1.98; Table 2).

Overall, all exercise metrics (endurance, number of turns to reach fatigue, and intensity of exercise) were affected by treatments across locations (3-way ANCOVA, all P<0.0001; Fig. 3). There was a significant interaction between location and temperature in time to fatigue (P=0.0005) and number of turns to fatigue (P=0.005; Fig. 3). Additionally, significant interactions were found in the intensity of exercise, in particular between location (or mass) and pH (P=0.02), and location (or mass), temperature and pH level (P=0.02). In this study, body mass was closely linked to location of origin, regardless of treatment (P<0.0001); therefore, it is not surprising that the two factors show a similar effect on the physiological responses tested. Juvenile skates from GB recovered faster from exercise than skates from the GoM (Fig. 4). Even though acidification increased recovery time across locations, skates from the GoM exhibited elevated metabolic rates for as much as double the time of individuals from GB under higher temperatures and acidification conditions (Dunnett's test, P<0.05; Fig. 4). Metabolic rates returned to RMR levels within 50 min of exercise (repeated measures ANOVA, Dunnett's test, P<0.05; Fig. 4).

This study is the first to investigate the effects of ocean acidification and warming on the performance of an elasmobranch sampled at different localities and developmentally acclimatized in a common garden experimental set-up. The present multistressor study demonstrated that even within a small geographic distance, it is possible to find locally adapted individuals of mobile marine organisms that respond differently to two major climatic stressors. Growth performance in L. erinacea from two locations raised in common garden conditions exhibited the same thermal optimum but GoM skates grew to a larger body size at the same age compared with those from GB across the experimental temperatures used in this study, suggesting countergradient variation in performance. These findings confirm previous experimental work on L. erinacea embryos that found strong evidence for local adaptation in body size (Di Santo, 2015) as well as field observations that described the same pattern but could not discern the effect of plasticity from adaptation (Frisk and Miller, 2006, 2009). Here, fish from two geographic locations were raised under the same conditions as newly laid embryos, therefore eliminating the effect of different acclimatization levels on individuals. Although embryos were obtained from wild-caught females held in different laboratories, these were kept at about 15–16°C, reducing the potential effect of different transgenerational acclimation in the fish (Donelson et al., 2012; Ho and Burggren, 2009). Furthermore, each egg was obtained from a different female and yet had similar yolk size, suggesting that parental body condition was not a significant factor in development (Di Santo, 2015; Vleck and Vleck, 1986).

Interestingly, embryonic little skates from both the GoM and the GB exhibit a metabolic performance thermal optimum around 18°C (Di Santo, 2015), while the results from the present study showed that the thermal optimum for performance is shifted ontogenetically towards lower temperatures in the GoM skates but not in the GB skates, as they move from the embryonic to the juvenile stage. These results seem to suggest that GoM skates may be even more vulnerable to warming and acidification as juveniles than embryos, albeit the effect on fitness may be only sub-lethal. Similar patterns of variation in life history traits have developed in several fishes along latitudinal gradients as the result of thermal adaptation in growth and metabolic rates. For instance, a similar pattern of variation was observed in the growth of the Atlantic silverside Menidia menidia along the coast in the northwestern Atlantic (Baumann and Conover, 2011) and in the respiratory performance of several coral reef fishes in the Great Barrier Reef (Gardiner et al., 2010). However, a larger body size could potentially have important trade-offs (Fangue et al., 2009). In fact, AAS was greater in skates from GB while it was significantly depressed in the GoM skates across the same temperature range. Conversely, GoM skates were able to endure the escape response for a longer period of time and righted themselves more times before fatiguing. The greater AAS measured in the GB skates may provide the capacity to recover from a more intense, but shorter, escape response. Performance in juveniles from GB peaked at 18°C, while GoM skates showed a decline in all performance parameters with temperature increase. Rates of physiological processes typically show a 2- to 3-fold increase for each 10°C increase in ambient temperature (Schmidt-Nielsen, 1990). However, in this study, RMR between the low and high temperatures increased dramatically (Q₁₀>5). This increase in RMR Q₁₀ values and the relative thermal insensitivity of MMR explain the decline of AAS with temperature (with the exception of GB skates at control pH; Fig. 2, Table 2). High Q₁₀ values for RMR have been attributed to fish species that live in stable environments (Rummer et al., 2014) as well as species that exploit thermal variability in their environment to enhance physiological processes (Hopkins and Cech, 1994). Often, Q₁₀ values for RMR do not directly correlate with the ability to enhance specific physiological processes (e.g. swimming or digestion efficiency; Di Santo and Bennett, 2011b), but as the results from this study show, they might explain the reduction in AAS with increasing temperature. Consequently, even though a 3–5°C increase in average temperature was not lethal for juvenile skates from the GoM, it reduced endurance during the escape response to chasing, prolonged recovery time after exhaustion and lowered aerobic scope with the potential to decrease fitness and local persistence in near-future projected climate change. Moreover, in both localities, low pH increased the recovery time after exhaustive exercise. Possibly, acidification may have caused an elevation in P_CO2 and HCO₃⁻ in the plasma (Green and Jutfelt, 2014), thereby increasing the costs associated with maintaining homeostasis. It is therefore possible that these additional costs might have significantly prolonged the time to recovery as young skates had to ‘pay off’ the oxygen debt after exhaustive anaerobic exercise (Svendsen et al., 2012; Webb, 1994).

High-energy output is required during anaerobic bursts (Webb, 1994) and these are typically associated with fast escapes; thus, reduced capacity for intense activity is likely to impact an individual's ability to survive predation (Allan et al., 2014; Johnson and Bennett, 1995). Although righting and escape bursts involve anaerobic pathways, these will be negatively impacted by a reduced aerobic scope (Svendsen et al., 2012). In fact, intensive anaerobic exercise results in an oxygen debt that must be accounted for at the expense of other important activities such as foraging. It is possible, however, that some fishes at higher latitudes would benefit from allocating most of their energy towards growth at the expense of other metabolic activities. For example, this might allow skates from the GoM to develop faster to compensate for a shorter growth season. The drop in escape performance of GoM skates would manifest at temperatures that are just a few degrees (∼3°C) higher than those currently experienced by individuals at this particular location, thereby pushing performance beyond its thermal optimum. It has been suggested that the negative effects of climate change stressors on species may be mitigated through geographic shifts toward cooler areas (Burrows et al., 2011; Chen et al., 2011; Genner et al., 2010; Somero, 2010). However, generally, skates show strong site fidelity, and some species do not seem to recolonize nearby areas where other populations were extirpated (Dulvy and Reynolds, 2002; Dulvy et al., 2003), but more empirical data on this particular species are needed. In addition, another coping mechanism could involve a shift in phenology, particularly for reproduction, to adjust to warmer temperatures (Gardner et al., 2011). However, L. erinacea has a long embryonic stage (Di Santo et al., 2016) and lays eggs year-round, so temperature is not a strong factor in determining reproductive timing and behavior (Palm et al., 2011).

Although increasing ocean acidification and warming are likely to have substantial effects across species at different latitudes, we need a better understanding of which particular traits can confer an advantage to more tolerant physiotypes (Lauder and Di Santo, 2015). The results from this study underscore the importance and benefits of employing multistressor experimental designs, in which crucial physiological traits are evaluated in organisms from different latitudes and populations (Todgham and Stillman, 2013). Furthermore, by reducing aerobic scope and escape endurance, increases in temperature and acidification beyond optimal conditions are likely to compromise vital activities such as predator evasion in fishes. This study shows that skates from the GoM tend to be more sensitive to acidification than those from GB, and their performance declines at temperatures above 15°C. The responses of the GB skates to increased acidification and temperature are more promising; skates from this locality recover much faster after exhaustive exercise, and acidification exerts a less detrimental effect on metabolic rates and therefore escape responses. It is possible that L. erinacea from GB might be already ‘pre-adapted’ to acidification because of the elevated and frequent fluctuations in pH in their environment as a consequence of strong upwelling in the area (Pershing et al., 2001). Alternatively, it is possible that given that GoM skates have a reduced aerobic capacity, they might be more strongly affected by elevated plasma P_CO2 and HCO₃⁻ (as seen in other elasmobranch species; Green and Jutfelt, 2014) as a consequence of increased acidification.

Lastly, variation in size between individuals complicates the interpretation of the results, as body mass is closely linked to location of origin, regardless of treatment, and larger fish (from the northern locality) generally tend to exhibit higher sensitivity to sub-optimal environmental conditions than smaller conspecifics (from the more southern locality) (Pörtner, 2010). Sensitivity to warming and acidification could perhaps be a function of local adaptation of body size, rather than metabolic adaptation of the escape response to local environmental conditions. A shift towards smaller body size has been invoked as one of the universal responses to global warming (Gardner et al., 2011). In fact, smaller fishes are known to be more tolerant and to perform better at higher temperatures than larger conspecifics (Di Santo and Lobel, 2016; Pörtner, 2010). In this study, skates from two locations exhibited differences in body size as well as physiological responses to environmental challenges that are not plastic but rather seem to be genetically fixed as developmental acclimatization to common garden conditions did not eliminate them. However, as body mass may affect physiological responses of fish species to ocean warming and acidification, future studies should also include analyses of different size classes rather than a single one, if the objective is to forecast more realistic responses. Ultimately, the differences in escape performance observed in L. erinacea from two locations underscore the importance of investigating intraspecific variation in physiological responses of species challenged by climate change stressors.

I thank George Lauder for critical logistic and conceptual support, James Sulikowski, NOAA and the Marine Biological Laboratory for providing some of the skates used in this study, and Phil Lobel for sharing laboratory space. Anna Tran helped during trials. Chris Pomory provided helpful suggestions during data analysis. Eric Widmaier, John Mandelman, Phil Lobel, Pam Templer, Jud Crawford and two anonymous reviewers provided useful suggestions and comments on a previous version of this manuscript.