Elasmobranch Responses to Experimental Warming, Acidification, and Oxygen Loss—A Meta-Analysis

Santos, Catarina Pereira; Sampaio, Eduardo; Pereira, Beatriz P.; Pegado, Maria Rita; Borges, Francisco O.; Wheeler, Carolyn R.; Bouyoucos, Ian A.; Rummer, Jodie L.; Frazão Santos, Catarina; Rosa, Rui

doi:10.3389/fmars.2021.735377

SYSTEMATIC REVIEW article

Front. Mar. Sci., 01 October 2021

Sec. Global Change and the Future Ocean

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.735377

Elasmobranch Responses to Experimental Warming, Acidification, and Oxygen Loss—A Meta-Analysis

$\nCatarina Pereira Santos$ Catarina Pereira Santos¹^*

Eduardo Sampaio^1,2,3

Beatriz P. Pereira¹

Maria Rita Pegado¹

Francisco O. Borges¹

Carolyn R. Wheeler^4,5

Ian A. Bouyoucos⁶

Jodie L. Rummer^5,7

Catarina Frazão Santos^1,8

Rui Rosa¹

¹Laboratório Marítimo da Guia, Marine and Environmental Sciences Centre (MARE), Faculdade de Ciências, Universidade de Lisboa, Cascais, Portugal
²Department of Collective Behaviour, Max Planck Institute of Animal Behavior, University of Konstanz, Konstanz, Germany
³Centre for the Advanced Study of Collective Behaviour, University of Konstanz, Konstanz, Germany
⁴School for the Environment, The University of Massachusetts Boston, Boston, MA, United States
⁵ARC Centre of Excellence for Coral Reef Studies, James Cook University, Townsville, QLD, Australia
⁶Department of Biological Sciences, University of Manitoba, Winnipeg, MB, Canada
⁷Marine Biology, College of Science and Engineering, James Cook University, Townsville, QLD, Australia
⁸Environmental Economics Knowledge Center, Nova School of Business and Economics, New University of Lisbon, Carcavelos, Portugal

Despite the long evolutionary history of this group, the challenges brought by the Anthropocene have been inflicting an extensive pressure over sharks and their relatives. Overexploitation has been driving a worldwide decline in elasmobranch populations, and rapid environmental change, triggered by anthropogenic activities, may further test this group's resilience. In this context, we searched the literature for peer-reviewed studies featuring a sustained (>24 h) and controlled exposure of elasmobranch species to warming, acidification, and/or deoxygenation: three of the most pressing symptoms of change in the ocean. In a standardized comparative framework, we conducted an array of mixed-model meta-analyses (based on 368 control-treatment contrasts from 53 studies) to evaluate the effects of these factors and their combination as experimental treatments. We further compared these effects across different attributes (lineages, climates, lifestyles, reproductive modes, and life stages) and assessed the direction of impact over a comprehensive set of biological responses (survival, development, growth, aerobic metabolism, anaerobic metabolism, oxygen transport, feeding, behavior, acid-base status, thermal tolerance, hypoxia tolerance, and cell stress). Based on the present findings, warming appears as the most influential factor, with clear directional effects, namely decreasing development time and increasing aerobic metabolism, feeding, and thermal tolerance. While warming influence was pervasive across attributes, acidification effects appear to be more context-specific, with no perceivable directional trends across biological responses apart from the necessary to achieve acid-base balance. Meanwhile, despite its potential for steep impacts, deoxygenation has been the most neglected factor, with data paucity ultimately precluding sound conclusions. Likewise, the implementation of multi-factor treatments has been mostly restricted to the combination of warming and acidification, with effects approximately matching those of warming. Despite considerable progress over recent years, research regarding the impact of these drivers on elasmobranchs lags behind other taxa, with more research required to disentangle many of the observed effects. Given the current levels of extinction risk and the quick pace of global change, it is further crucial that we integrate the knowledge accumulated through different scientific approaches into a holistic perspective to better understand how this group may fare in a changing ocean.

Introduction

The rise in anthropogenic greenhouse gas emissions since the Industrial Revolution, particularly carbon dioxide (CO₂), is likely to shape the future of life on Earth (Pecl et al., 2017; Bindoff et al., 2019). The ocean has been assuaging the surface-level repercussions of increased emissions, namely by absorbing over 90% of the energy trapped by the excess greenhouse gases and up to 30% of the CO₂ released to the atmosphere. However, the consequences of this key service are brewing beneath the surface, with ocean warming (OW) and acidification (OA) emerging as major threats to marine life and dependent communities (Pecl et al., 2017; Bindoff et al., 2019; Collins et al., 2019). According to the most pessimistic forecasts by the Intergovernmental Panel on Climate Change [IPCC; representative concentration pathway (RCP) 8.5], mean sea surface temperatures are projected to increase up to nearly 4°C, while mean pH may potentially drop by over 0.3 units by the end of the century (Abram et al., 2019). Meanwhile, rising temperatures have been instigating an additional facet of climate change with potentially widespread consequences: ocean deoxygenation (OD). Indeed, increasing temperatures represent a major driver of oxygen loss by reducing oxygen solubility, promoting stratification, slowing ocean circulation, and increasing heterotrophic respiration (Oschlies et al., 2018; Laffoley and Baxter, 2019). Under the same RCP scenario and timeframe, mean dissolved O₂ in the upper ocean layers is expected to decrease by up to 5% (Abram et al., 2019).

Alongside gradual changes in mean conditions, extreme events associated with the same three factors—warming, acidification, and deoxygenation—are also expected to become more frequent and widespread (Bijma et al., 2013; Bindoff et al., 2019; Collins et al., 2019; Burger et al., 2020). Over the past decades, marine heatwaves (MHWs) frequency and duration increased by over 30% and nearly 20%, respectively, in a trend expected to continue in the near future (Oliver et al., 2018; Collins et al., 2019). Likewise, hypoxic (HEs) and extreme acidification (EAEs) events are expected to become more prevalent over the present century due to climate change (Bijma et al., 2013; Bakun et al., 2015; Burger et al., 2020), compounding the impacts associated with nutrient runoff (Bijma et al., 2013; Duarte et al., 2013; Gissi et al., 2021). These phenomena can feature transient shifts in temperature, pH, and oxygen levels that far surpass IPCC long-term projections, with potentially devastating consequences (Leggat et al., 2019; Schulz et al., 2019; Sampaio et al., 2021). While acting at distinct temporal and spatial scales, both extreme and long-term environmental shifts intertwine and have the potential to reshape marine ecosystems, with considerable socio-economic consequences (Smale et al., 2019; Burger et al., 2020; Cheung and Frölicher, 2020). Thus, it is increasingly urgent that we understand how marine life will respond to this “deadly trio” (Bijma et al., 2013; Sampaio et al., 2021).

With early records dating to the Lower Silurian, Chondrichthyes represent one of the most ancient and diverse vertebrate lineages, having radiated to occupy a wide range of ecological niches and play key functional roles in aquatic ecosystems (Grogan et al., 2012; Hammerschlag et al., 2019). However, while the long evolutive history of Chondrichthyes—and their most speciose subclass, elasmobranchs—vouches for their resilience toward changing environmental conditions, the array of pressures brought by the Anthropocene have been pushing this group toward a worldwide decline (Dulvy et al., 2014; Queiroz et al., 2019; Pacoureau et al., 2021). Currently, sharks and their relatives represent one of the most threatened vertebrate groups (Dulvy et al., 2014) and are at disproportional risk of functional diversity and evolutionary history loss (Stein et al., 2018; Pimiento et al., 2020). While the most consequential driver of chondrichthyan extinction risk is definitely overfishing (Dulvy et al., 2014), global change may further jeopardize this group. In fact, climate change is likely to affect this group both through direct impacts over their biological responses (Chin et al., 2010; Rosa et al., 2017; Bouyoucos et al., 2019; Wheeler et al., 2020) and indirect effects due to the disruption of ecological dependences, including changes in habitat quality, ocean productivity, and prey availability (Chin et al., 2010; Schlaff et al., 2014; Bindoff et al., 2019). Moreover, the expected changes may increase their exposure to other anthropogenic threats, namely fishing pressure (Vedor et al., 2021), or compromise the efficacy of conservation measures (Davies et al., 2017). Additionally, their k-selected life strategy, with long generation times and low fecundity, reduces their margin for adaptation to rapid environmental change (Wheeler et al., 2020). On the other hand, sharks' high position in marine trophic webs renders them as potential climate change multipliers (Zarnetske et al., 2012), with previous studies showcasing their potential to moderate ecological responses to extreme events (Nowicki et al., 2019, 2021). Hence, understanding how sharks and their relatives will respond to warming, acidification, and deoxygenation, represents a priority in conservation physiology.

Experimental exposure to warming (EW), acidification (EA), and deoxygenation (ED) in a controlled setting represents a precious tool to understand the effects of climate change over marine life. While entailing a substantial trade-off with environmental realism and sidestepping ecological dependences, EW, EA, and ED allow us to disentangle the influence of the corresponding factors over marine biota from a myriad of confounding variables, and better understand the mechanisms behind the observed effects (Baumann, 2019; Sampaio et al., 2021). While experimental research addressing the effects of these three factors on elasmobranchs has lagged behind more amenable taxa, a considerable body of research has been building for EW and EA over the past decade. In fact, an extensive array of effects has been reported in the context of exposure to both EW and EA, from biochemical and metabolic adjustments to altered behavior and survival rates (e.g., Green and Jutfelt, 2014; Heinrich et al., 2014; Rosa et al., 2014; Di Santo, 2016; Pegado et al., 2018), although trends may not always be consistent (Rosa et al., 2017). Sustained ED, meanwhile, has been mostly neglected despite the potential for steep consequences, as ubiquitously observed for aerobic organisms (Sampaio et al., 2021). In fact, across all factors, several vital biological responses have seldom been addressed, with experimental research being hampered by a wide range of constraints, from logistic limitations to conservation and ethical concerns. The consequences of data paucity are amplified by unique physiology and life traits of elasmobranchs, which complicate extrapolation from other groups (Rosa et al., 2017). As such, the identification of knowledge gaps, along with the targeted focus toward the most pressing concerns, is particularly critical for this group. On the other hand, while experimental research provides a central stepping-stone to understand the impacts of climate change over this taxon, it is crucial that we further integrate these insights with the information collected through different approaches (e.g., Chin et al., 2010; O'Brien et al., 2013; Lighten et al., 2016; Birkmanis et al., 2020; Vedor et al., 2021).

Meta-analyses facilitate the synthesis and integration of available research, allowing us to quantitatively realize the consistent trends and overarching effects derived from a pool of highly context-specific empirical datasets (Rosa et al., 2017; Sampaio et al., 2021). In this context, we systematically surveyed the currently available literature for studies featuring the sustained (>24 h) and controlled exposure of sharks and their relatives to EW, EA, and/or ED. Using a comprehensive comparative framework, we employed sequential meta-analyses to assess the effects of these factors over elasmobranchs. However, the adaptations necessary to thrive in a specific climate, the traits accumulated by distinct lineages over geological time, the specialization associated with different lifestyles and habitats (O'Brien et al., 2013; Nadeau et al., 2017; Baumann, 2019), and the contingencies related to distinct reproductive modes (Wheeler et al., 2020), may all conceivably modulate species susceptibility toward rapidly changing environmental conditions. Likewise, different ontogenetic stages often explore distinct habitats and express different levels of tolerance to environmental change, thus, the identification of ontogenetic bottlenecks of performance represents an important undertaking (Baumann, 2019; Dahlke et al., 2020). As such, in addition to comparing the overall magnitude of change imposed by the different treatments over elasmobranchs in general, we looked for potential differences between groups with distinct attributes, namely climate of origin, order, lifestyle, reproduction mode, and life stage. On the other hand, only by scrutinizing the directional effects of these factors over specific biological responses can we understand the nature of the changes taking place. Hence, we assessed the effects of a comprehensive array of biological responses, specifically survival, development, growth, aerobic metabolism, anaerobic metabolism, oxygen transport, feeding and digestion, behavior, acid-base status, thermal tolerance, hypoxia tolerance, and cell stress. With this approach, we aim to provide a quantitative overview of the effects observed in a context-specific manner and communicate priorities for future research.

Materials and Methods

Literature Search

We surveyed the literature for peer-reviewed studies featuring the controlled and sustained (>24 h) exposure of elasmobranchs to increased temperature (EW), reduced pH through increased pCO₂ (EA), and/or reduced oxygen levels (ED) at the whole-body level. We searched across all databases in ISI Web of Knowledge (last search January 21st, 2021) using an array of pertinent keywords (“shark” OR “elasmobranch” AND “climate change” OR “ocean warming” OR “ocean acidification” OR “ocean deoxygenation” OR “temperature” OR “hypercapnia” OR “hypoxia” OR “°C” OR “CO2” OR “pH” OR “O2”). This search was complemented using Google scholar and reviewing the references from relevant reviews (Rosa et al., 2017; Bouyoucos et al., 2019; Wheeler et al., 2020).

Distinct exposure events and different species or populations tested within the same paper were classified as different studies, each including one or more experimental treatments contrasted with the same set of control animals. Excluded were studies (or individual responses) that (i) lacked a suitable control treatment, namely no treatment within the natural range of the species or the presence of confounding factors (e.g., distinct batches); (ii) lacked any treatment deviating from control conditions in the direction expected in a climate change context (e.g., temperature decrease only); (iii) exposure lasted less than 24 hours; (iv) featured overly invasive procedures with no ecological relevance (e.g., surgeries); (v) featured cyclical or intermittent treatments; (vi) did not report quantitative data, namely dispersion measures, or data could not be confidently retrieved. For details regarding the number of studies screened and removed at each stage see the flowchart provided in Supplementary Figure 1, prepared according to PRISMA (Preferred Reporting Items for Systematic Reviews and Meta-Analyses) guidelines (Moher et al., 2009; Page et al., 2021). The dataset retrieved builds on an earlier compilation by Rosa et al. (2017), broadening its scope to all elasmobranchs and implementing a more comprehensive analysis to address the effects of both warming, acidification, and deoxygenation - in a framework akin to Sampaio et al. (2021).

Data Compilation

Suitable studies (Table 1; Figure 1A) were organized according to the experimental treatments tested: EW, EA, ED, and combined treatments. Specifically, combined treatments consisted of EWA (experimental warming and acidification) and EWD (experimental warming and deoxygenation), with no studies experimentally addressing the combined effects of all three factors (EWAD), nor the combination of acidification and deoxygenation (EAD).

TABLE 1

Table 1. Summary of the references including in the analysis across the different treatments, namely experimental warming (EW), acidification (EA), deoxygenation (ED), and their combination (EWA and EWD).

FIGURE 1

Figure 1. (A) Number of suitable studies (and control-treatment contrasts) per treatment, featuring experimental warming (EW), acidification (EA), and deoxygenation (ED), along with their combination (EWA, EWD, and EWAD). (B) Elasmobranch overall response to each experimental treatment standardized according to the stressor ratio (SR) implemented in each study. Effect sizes represent the absolute value of the response ratio (-|LnR|) and 95% confidence intervals (CI). Significant effects, marked with asterisks, correspond to the absence of overlap between the CIs and the threshold used to test the null hypothesis (dashed line; effect size = 0). Different letters showcase differences between treatments; numbers indicate sample sizes (control-treatment contrasts). Statistical results are detailed in Supplementary Table 1. With only two studies available for ED, the obtained results should be interpreted with special caution. With only one study, WD was entirely excluded from the analysis. (C) Number of studies per treatment featuring a lower (light) or higher (dark) SR than the worst-case scenario (Representative concentration pathway 8.5) predicted by the Intergovernmental Panel on Climate Change for the end of the century.

For each suitable study, controls were set as the lowest temperature/ highest pH/ highest O₂ treatment within the species natural range. Measures of central tendency and dispersion, sample sizes, and treatment conditions, were collected or computed from raw data, text, tables, or plots, using the software Im2graph (v1.21) to obtain data points when necessary. Alternative central and dispersal measures were converted to means and standard deviations prior to analyses using the applicable mathematical formulas; specifically, the QQ model was used to retrieve estimates from boxplots, when necessary (McGrath et al., 2020). In the case of unclear sample sizes per treatment, the most conservative interpretation of the information provided was used to compute the final effect sizes. Whenever dispersal was zero in the control treatment (e.g., 100% survival across replicates), the average dispersal observed for the same response across the other treatments in the same study was applied to the control, thus allowing the conservative computation of effect sizes. The endpoints analyzed in each suitable study were screened and classified into the response categories, as specified in Figure 2 and Supplementary Data 1. Also retrieved, both through author-provided information and alternative sources [namely: Froese and Pauly (2021); Ebert et al. (2013); Nakaya et al. (2020); and Sternes and Shimada (2020)], was information regarding the model organism used in study, specifically pertaining to their (i) taxonomy, (ii) climate of origin, (iii) ontogeny stage, (iv) lifestyle, and (v) reproduction mode (Figure 2 and Supplementary Data 1).

FIGURE 2

Figure 2. Number of suitable studies assessing each factor [experimental warming (EW), acidification (EA), deoxygenation (ED), along with their combinations (EWA and EWD)] across different attributes of the model organisms used, namely their climate (A), order (B), lifestyle (C), reproduction mode (D), and life stage (E), as well as (F) the response categories assessed.

Statistical Analysis

Analysis Framework

A series of meta-analyses were performed using mixed-effects models to test the effect of distinct moderator variables over different subsets of data depending on the research question. First, the overall response of elasmobranchs was compared across the different treatments, combining the responses from all categories (Stage I). Then, the overall response to each specific stressor was compared across different attributes (order, climate, lifestyle, reproduction mode, and life stage) of the model organisms (Stage II). Lastly, each biological response was individually compared across treatments (Stage III).

All analyses were carried out in R (R Core Team, 2021; version 4.0.5) and R studio (RStudio Team, 2019; version 1.2.5019). Meta-analyses were implemented using the metafor package (Viechtbauer, 2010), in a framework adapted from Sampaio et al. (2021). Exploration, integration, and visualization of the data largely relying on the ecosystem of packages within tidyverse for the (Wickham et al., 2019). Additional packages are detailed below and fully listed along with the script used in the analysis (see Supplementary Code 1).

Effect sizes and variance estimates were calculated for each control-treatment comparison using the escalc function and the ROM measure, which calculates the natural log-transformed ratio of the means (LnR) between control and treatment (Hedges et al., 1999; Lajeunesse, 2011; Viechtbauer, 2020). Depending on the analysis stage, treatment or the different attributes were used as categorical moderators (detailed in corresponding sections). To fairly compare between treatments and account for the different stressor deltas used across studies in a standardized framework, a stressor ratio (SR) was included in the analysis. Specifically, the delta between control and treatment conditions for the factors tested in each study (EW, EA, ED) was calculated and standardized as the percentage of change in relation to the most pessimistic IPCC projections (EW: + 3.5°C; EA: −0.317 pH units; ED: −5.0% change [O₂]; RCP8.5; Abram et al., 2019). In the case of treatments with factor combinations, the average of both stressor ratios was used. The natural logarithm of the SR (LnSR) was computed and included as an interacting moderator across all models (LnSR calculation and values used are detailed in Supplementary Data 1). Additionally, “-1” was included along with the categorical moderator to generate a dummy variable to which each level of the moderator was compared, directly testing the null hypothesis. To deal with the uncertainty associated with small sample sizes, the significance of effect sizes and confidence intervals, calculated by restricted maximum likelihood, was verified using t-statistic instead of the default z-statistic, retrieving more conservative results (Knapp and Hartung, 2003; Viechtbauer, 2020).

To minimize pseudo-replication, in studies where several condition deltas on the same treatment were assessed, control-treatment contrasts were run solely using data obtained under the treatment conditions corresponding to the SR closest to IPCC predictions. Furthermore, while several responses within each response category were retrieved and separately analyzed in Stage III, to provide a more comprehensive overview of the literature, only one response per response category was included in Stage I and Stage II, preferentially the most frequently reported metric within each category (e.g., routine oxygen consumption, within aerobic metabolism). Likewise, for experiments reporting several time points, only the results obtained following the longest exposure period were considered. To further mitigate the potential effects of pseudo-replication, treatment, response category, study, and paper were all introduced as random effects whenever appropriate (detailed below, in the stage-specific sections). These effects were hierarchically structured, with high complexity structures prioritized to achieve high independence between stressors, responses, and studies. Most models were run using an unstructured variance/covariance matrix (UN), and the lowest structure level implemented was compound symmetry (CS; Viechtbauer, 2020; Sampaio et al., 2021). A thorough sensitivity analysis was also conducted (see section Publication Bias and Sensitivity Analysis).

For all stages, pairwise post hoc comparisons were drawn between the different levels of the categorical moderator (e.g., between different treatments; Gurevitch and Hedges, 1999; Ferreira et al., 2015; Sampaio et al., 2021). Tukey's honest significance tests were employed using general linear hypotheses across a contrast matrix to run pairwise comparisons between levels; this was performed using the multcomp package (Hothorn et al., 2008), as described by Sampaio et al. (2021).

Stage I

One response per response category for each study and treatment (EW, EA, ED, and EWA) were pooled and compared across treatments by including the different treatments as distinct levels of the treatment moderator. With only one study to date, EWD was excluded from the analysis. In addition to the inclusion of “1 | Paper” as a random effect, two additional random effects were hierarchically included in the analysis, specifically “Treatment | Study” and “Treatment | Response category,” to account for the use of several treatments from the same study (which make use of the same set of control animals) and response category (which may be differently affected by the distinct treatments). Moreover, given that several responses were included in the analysis, the absolute value of the response ratio was calculated (-|LnR|) and used in the model, thus providing an absolute measure of deviation from control, regardless of direction of response (e.g., temperature increase typically results in a decrease in development time and an increase in routine metabolism). In this context, effect sizes significantly differing from 0 indicate a significant impact of treatment over the biological response.

Stage II

The dataset used in Stage I was subset according to treatment and a series of meta-analyses were run using different attributes of the model organisms as categorical moderators, specifically climate, taxonomic order (as a proxy for phylogenetic lineage), lifestyle, reproduction mode, and life stage. In this case, given that only one treatment was tested in each model and each study featured only one level of each moderator, the random effects included in the analysis were “1 | Study,” “1 | Response category,” and “1 | Paper.” With several responses once again included in the analysis, the absolute value of the response ratio was used in these models and effect sizes should be interpreted as described in Stage I.

Stage III

Each biological response was separately compared across treatments. Treatment was once again included as the moderator and, given that only one response per study and treatment was inherently included, only “Treatment | Study” and “1 | Paper” were included as random effects. Individual analyses were run for multiple responses within the same category (e.g., length and weight, within growth) whenever possible, to provide a more comprehensive overview of the literature. Given that only one response type is, for the most part, included in each analysis, the response ratio was directly inputted in the models, allowing the evaluation of the direction of effect. Hence, in this stage, effect sizes above 0 indicate a response stimulation while values below 0 correspond to response inhibition. An exception was, nonetheless, made in the case of the behavior response, which features distinct types of behavioral responses and, thus, the absolute response ratio was used in the analysis. In the latter case, effect sizes should be interpreted as described in Stages I and II.

Publication Bias and Sensitivity Analysis

Potential publication bias was assessed through the estimation of the Rosenthal fail-safe number, which estimates how many effect sizes with no effect would be necessary to alter the significance value provided by the model (Rosenberg, 2005; Viechtbauer, 2010; Supplementary Table 1). Furthermore, for each level of the moderator, we computed Duval and Tweedie's Trim and Fill operations, which represents a non-parametric method to estimate the number of missing effect sizes due to the suppression of the most extreme points on each side of the funnel plot. This provides conservative estimates, given they do not account for the random effects included in the main analysis (Duval and Tweedie, 2000; Viechtbauer, 2010; Supplementary Code 1). Additionally, to assess the robustness of the significant results obtained, forest plots were visually screened for potentially influential data-points, and leave-one-out diagnostics were performed using the Cook's distance and DFBETAS statistic [see Viechtbauer and Cheung (2010)]. Potentially influential datapoints were tentatively removed, one-by-one, to assess their effect on the overall mean effect size and direction of the trends observed; this was performed for all models with N > 5. Additionally, whenever sufficient studies were available (corresponding to a N_paper > 3), papers contributing with more than one study were removed one-by-one. Finally, to account for the potential effects of non-independence stemming from the use of multiple responses from the same study, models were run using only one weighed effect size per study and stressor, obtained from mixed effect models computed using “study per moderator level” as the categorical moderator and “Treatment | Paper” as a random effect (Ferreira et al., 2015). This analysis was performed for all models that featured multiple responses per study (Stage I and II) and included at least three studies. No changes in the significance of each level in relation to the main models were observed in most cases. The exceptions were Climate EWA, Order EA, and Life stage EA; although, even here, a very close and similar trend was observed. Hence, the results of the main models are provided and discussed here, and the independence analysis is discriminated in the Supplementary Code 1.

Results and Discussion

Stage I: General Response

Temperature has long been acknowledged as a key modulator of biological processes and its effects have been experimentally studied for decades (Britton, 1924; Neale et al., 1977). Not surprisingly, of the climate change-related factors considered here, the most often investigated in an experimental setting is EW (Figure 1A). However, the first studies explicitly addressing the effects of EW over elasmobranchs in a climate change-context were rather recent (Dabruzzi et al., 2013; Rosa et al., 2014). When evaluating the overall response to each treatment in a comparative framework, warming arises as the factor with the most extreme effects over elasmobranchs (Figure 1B; Supplementary Table 1).

Shark competence at acid-base regulation has been previously established in the literature, leading to an earlier assumption of minimal direct effects of OA over this group (Evans et al., 2004; Chin et al., 2010). However, earlier studies assessing elasmobranchs' responses to high CO₂ levels (e.g., Claiborne and Evans, 1992) mostly addressed the physiological effects of acute and often extreme hypercapnia, using experimental designs with little ecological significance in a climate change-context (Baumann, 2019; Bindoff et al., 2019). Research aiming to understand EA effects over this group in an OA-context only became available in recent years (Green and Jutfelt, 2014; Heinrich et al., 2014; Rosa et al., 2014; Di Santo, 2015), revealing potential secondary effects of EA (Rosa et al., 2017). Indeed, albeit modest in comparison to the effects observed for EW and EWA, considerable effects of EA over elasmobranchs emerge from the pool of experimental research conducted to date (Figure 1B; Supplementary Table 1).

Studies implementing factorial experimental designs are also restricted to the last decade and far more sparce than studies featuring isolated factors (Figure 1A). With both factors expected to co-occur over the foreseeable future, these studies provide crucial insights into what lies ahead. Based on the current analysis, the combined effects of EW and EA (EWA) closely match those observed under the sole exposure to EW (Figure 1B; Supplementary Table 1). While the results of these studies are known to be particularly context-dependent (e.g., Pistevos et al., 2015, 2017), such a close match further pins temperature rise as the leading driver of change.

Oxygen loss is by far the most overlooked factor (Figure 1A). Indeed, currently available studies mostly feature acute (a few hours or less) or intermittent exposures (e.g., Rytkönen et al., 2012), which, similarly to early studies on high CO₂ effects, are not suitable for comparison with the remaining factors within the present framework. In fact, only two studies featured a sustained exposure to reduced levels of oxygen (Morash et al., 2020; Musa et al., 2020), one of them still fairly short in exposure duration (48H; Morash et al., 2020). While severely neglected, the effect of ED over elasmobranchs seems well within the range of the more well-studied EA. However, as emphasized in Figure 1B (Supplementary Table 1), the model outcome should be considered with special caution, given the current data limitations. While no significant effects were detected, a detrimental physiological response to lower levels of O₂ is to be expected for aerobic biota in general, namely at more extreme levels (Sampaio et al., 2021). This is particularly pertinent in a climate change context, given the projected increase in HEs pervasiveness and the expansion of oxygen minimum zones alongside long-term OD (Bijma et al., 2013; Bakun et al., 2015; Sampaio et al., 2021). Nonetheless, as ED inclusion in the model resulted in no substantial changes in the output regarding the remaining treatments (detailed in Supplementary Code 1), these results are provided for an overview of the current body of literature and to underscore the magnitude of this knowledge gap.

Data paucity extends to most multi-stressor treatments. While the interaction between ED and other factors has been on the radar of researchers for several decades (Butler and Taylor, 1975; Perry and Gilmour, 1996), only one study addressing EWD (Musa et al., 2020) and none regarding EDA or EWDA would be suitable for comparison within the present framework. This contrasts with EWA, which has been addressed by 17 studies between 2010 and 2020 (Figure 1A; Supplementary Table 1). Moreover, while roughly half of EW and EA studies implement SRs below the worst-case IPCC scenario for the end of the century (Figure 1C), the only two suitable ED studies to date feature SRs far beyond IPCC predictions, aligning with HEs rather than gradual OD. While extreme phenomena relating to warming (i.e., MHWs), acidification (i.e., EAEs), and deoxygenation (i.e., HEs) are all expected to become more pervasive over the next century, this approach disparity complicates the comparison between experimental treatments. This issue is, nonetheless, prevalent across all taxa in climate change-related experimentation (Sampaio et al., 2021) and compounds the general paucity of studies addressing the effects of sustained ED over elasmobranchs.

It should be noted that, by using the absolute value of change, the framework implemented at both this (Stage I), and the next (Stage II) stages of the analysis disregards the direction of the observed effects, compounding effects that move in opposite directions for the same response. On the other hand, this provides a more reasonable ground for comparisons across different responses (e.g., comparing the effects of EW over metabolism and development). In Stage III, directional analyses were performed for each specific response, to reveal the directional effects of these factors.

Stage II: Specific Attributes

Climate Region

Regarding climate affinity, subtropical elasmobranchs are by far the most studied to date (Figure 2A). Meanwhile, only one study featuring a single boreal species is currently available (Schwieterman et al., 2019), which does not allow for a meaningful grasp as to how species from high latitudes may be affected. This represents a considerable knowledge gap as, while relatively few elasmobranch species inhabit such high latitudes (Ebert et al., 2013), they may be particularly sensitive to climate change. Trapped by rising temperatures, these species are likely to see their suitable habitat restricted, having to cope with both climate change-induced physiological impacts and increasing competition from species moving in from lower latitudes (Berteaux et al., 2004; Somero, 2010; Poloczanska et al., 2016). On the other end, tropical climates are typically more stable, animals often living closer to their thermal limits. Hence, tropical species may also be particularly impacted by rising temperatures (Somero, 2010; Poloczanska et al., 2016). However, EW produced a widespread effect, with no detectable differences across the different climate affinities (Figure 3A; Supplementary Table 2). This may possibly reflect a low resolution provided by the available body of literature and the approach implemented here rather than the true absence of differences in the biological responsiveness across climate groups. On the other hand, while the effects of EWA were also ubiquitous across groups, the effects of EA significantly impacted both tropical and subtropical animals but not temperate animals (Figures 3B,C; Supplementary Table 2), in a similar pattern to that observed by Sampaio et al. (2021) across all life stages. This disparity may be related to the natural patterns of pH fluctuations naturally experienced by species in their habitat, with conditions at temperate regions being typically more variable and, particularly in coastal habitats, unpredictable (Duarte et al., 2013; Baumann, 2019). Nonetheless, with only two species of temperate elasmobranchs tested under EA to date, more research is required, particularly given that these regions represent a major carbon sink and experience a quicker pace of acidification (Bindoff et al., 2019).

FIGURE 3

Figure 3. Elasmobranch response to each experimental treatment [warming (EW); acidification (EA); and warming + acidification (EWA)] depending on the attribute [climate (A–C), order (D–F), lifestyle (G–I), reproduction mode (J–L), life stage (M–O)] of the model species used and standardized according to the stressor ratio (SR) implemented in each study. Effect sizes represent the absolute value of the response ratio (-|LnR|) and 95% confidence intervals (CI). The dashed line (R = 0) represents the threshold used to test the null hypothesis, with overlapping CI indicating the lack of treatment effect; asterisks indicate a significant effect. Different lowercase letters showcase differences between treatments; numbers indicate sample sizes (number of control-treatment contrasts). Statistical results are detailed in Supplementary Table 2.

Lineage

Within the class Chondrichthyes, the bulk of research to date has focused on the sub-class Elasmobranchii. Indeed, only one study featuring the exposure of chimeras (sub-class Holocephali) to EW would be considered suitable under the current framework (Lyon et al., 2011). Holocephali represents a rather evolutionarily distinct clade, with a basal position within the vertebrate tree of life, rendering it a particularly valuable model for future research (Venkatesh et al., 2014; Stein et al., 2018). Across Elasmobranchii, all experimental research conducted in a climate change context has focused on either the super-order Galeomorphii or Batoidea. Hence, although the current research body has been fairly split between sharks and batoids, an entire super-order of sharks has seldom been contemplated (Squalomorphii; Bockus et al., 2020). Additionally, no suitable study has evaluated the effects of these stressors in any lamniform species—although, whole-organism experimentation is, for the most part, out of reach for these large, active, and exceptionally threatened pelagic animals (Ebert et al., 2013; Dulvy et al., 2014). Hence, alternative approaches are vital to assess the climate change vulnerability of this clade. Within batoids, currently, no study has assessed the effects of either stressors over Torpediniformes, and only one investigated the effects of EW over Rhinopristiformes (Figure 2B; Du Preez et al., 1988).

The effects of EW appear to be pervasive across all orders assessed, except Rajiformes—although a close trend was still observed for this group (Figure 3D; Supplementary Table 2). Regarding the effects of EA, all orders except Carcharhiniforms showed a significant response (Figure 3E; Supplementary Table 2). This latter clade represents the most speciose and diverse shark order with a wide diversity of traits, from lifestyle to reproductive mode (Ebert et al., 2013). In fact, most benthopelagic and pelagic species tested to date belong to this order, with both viviparous and oviparous species represented (Supplementary Data 1). Similar to EW, all orders show a response to EWA, although Heterodontiformes showcase higher responsiveness to these conditions (Figure 3F; Supplementary Table 2).

Lifestyle

Few benthopelagic, and even fewer pelagic species, have been studied in comparison to benthic species (Figure 2C; Supplementary Data 1), the latter being typically more amenable to captive rearing. Hence, data from pelagic and benthopelagic species was merged and compared against benthic species. Once again, EW produced a substantial impact over both groups (Figure 3G; Supplementary Table 2), although the effects observed over pelagic/benthopelagic species were less robust, highlighting a higher heterogeneity within this group (see Sensitivity analysis; Supplementary Code 1). On the other hand, pelagic/benthopelagic species showed no response to EA, while a considerable effect was observed for benthic animals (Figure 3H; Supplementary Table 2). Given species' inherent dispersal capacity and the spatial and temporal heterogeneity they experience in relation to pH and pCO₂, one would expect an increase in sensitivity from benthic to benthopelagic and pelagic species, as previously observed for teleost fishes (Cattano et al., 2018). However, all pelagic and benthopelagic species studied to date are known to explore estuaries, some often performing freshwater incursions, thus likely representing particularly plastic species (Ebert et al., 2013; Nadeau et al., 2017; Baumann, 2019). Moreover, while pelagic species' association with more stable environments may underscore a higher degree of sensitivity toward acidification (Cattano et al., 2018), there are currently not enough data to empirically assess how they may be affected. In fact, contrary to coastal regions, which undergo frequent and often unpredictable pH/pCO₂, the open ocean is far more stable (Duarte et al., 2013; Hannan et al., 2020). Oceanic species may, thus, potentially represent the most EA sensitive group of elasmobranchs, with temperature-induced range shifts likely pushing these species toward higher latitudes, where they may experience OA at an increased rate (Bindoff et al., 2019). In this context, understanding how these species may respond to EA and EWA represents a rather challenging but important undertaking. Likewise, more research is necessary to understand how these, more active, species will be affected by EWA (Figure 3I; Supplementary Table 2).

Reproduction Mode

The diverse reproductive modes displayed by elasmobranchs can influence species vulnerability to anthropogenic pressures (García et al., 2008; Wheeler et al., 2020). In the context of climate change, it is worth to consider that the embryos of oviparous species are naturally exposed to high environmental fluctuations, while viviparous species benefit from an extended maternal protection. Indeed, while the progeny of viviparous and oviparous species was equally affected by EW (Figure 3J; Supplementary Table 2), viviparous species seem impervious to EA (Figure 3K; Supplementary Table 2). From an experimental design perspective, viviparous species are typically exposed for shorter periods. However, when running the analysis with the dataset restricted to individuals not exposed during embryogenesis, the same pattern was obtained (Supplementary Code 1). As such, this outcome may reflect the effects of increased maternal investment. Nonetheless, it should be considered that all but one study investigating the effects of EA over viviparous species features benthopelagic or pelagic species (Supplementary Data 1; Bouyoucos et al., 2020c). As such, it is currently not feasible to fully disentangle the effects of EA over different lifestyles and reproduction modes. More research is also necessary regarding the effects of EWA over viviparous species, with clear effects observed for oviparous species (Figure 3L; Supplementary Table 2).

Life Stage

Juveniles are typically more resilient to EW in comparison to both adults, burden by the energetic costs associated with reproduction, and embryos, still amidst cardiovascular development (Pörtner and Farrell, 2008; Dahlke et al., 2020). However, EW produced no detectable differences across elasmobranch developmental stages, affecting all groups regardless of exposure onset (Figure 3M; Supplementary Table 2). Significant effects of EA, on the other hand, were only detected for adults and animals exposed over their embryogenesis (Figure 3N; Supplementary Table 2). As previously mentioned, elasmobranchs strive to produce a particularly robust progeny, which, following extended embryogenesis and considerable maternal provisioning, are born as smaller versions of their adult counterparts (Carrier et al., 2004). On the other hand, elasmobranchs' rebound potential is closely tied with adult numbers and overall fitness (Smith et al., 1998). Meanwhile, despite the small number of studies under EWA precludes further analysis, animals exposed as embryos appear to be particularly affected by the combination of these factors—although a close trend was also observed for adults (Figure 3O; Supplementary Table 2). It should, nonetheless, be noted that elasmobranch studies featuring embryonic exposure to climate change-related conditions are restricted to oviparous species. These species have their energetic reserves limited by female provisioning during oviposition (i.e., the volume of the yolk-sac) and are essentially sessile and, thus, particularly susceptible to extreme events (Carrier et al., 2004; Wheeler et al., 2020).

Stage III: Specific Responses

Survival

Survival rates are ultimately an indicator of the strength of the selection pressure caused by environmental change. As no significant effects were observed over survival rates (Figure 4A; Supplementary Table 3) and hatching success (Figure 4B; Supplementary Table 3) across all analyzed treatments. Hence, projected levels of EW and EA for the end of the century may not directly impose a strong selection pressure over elasmobranchs, particularly considering their low generation turn-over (Wheeler et al., 2020). However, it should be noted that EW consistently resulted in lower survival rates across studies, many of which were not included in this analysis due to the lack of dispersal data for this specific response (e.g., Pouca et al., 2018; Hume, 2019; Musa et al., 2020). Additionally, in the only sustained study featuring ED, a considerable decline in survival was observed (Musa et al., 2020). Nonetheless, as k-strategists, sharks and their relatives are particularly susceptible to the levels of overexploitation currently taking place (Wheeler et al., 2020), which may further camouflage the minor selection pressure likely imposed by climate change. As such, potential sublethal effects are more likely to persist across generations.

FIGURE 4

Figure 4. Directional effect size (LnR) of experimental warming (EW), acidification (EA), and their combination (EWA) over elasmobranchs, considering specific biological responses [(A,B) survival; (C) development; (D–F) growth]. The dashed line (R = 0) represents the threshold used to test the null hypothesis, with overlapping confidence intervals (CI) indicating the lack of significant treatment effect. R values above 0 indicate a stimulation and below indicates an inhibition of the response; asterisks indicate a significant effect. Different letters showcase differences between treatments; numbers indicate sample sizes (number of control-treatment contrasts). Statistical results are detailed in Supplementary Table 3.

Development

Temperature is known to strongly influence ectotherm development and, consistently, both EW and EWA resulted in a significant decline in the duration of embryogenesis (Figure 4C; Supplementary Table 3). Indeed, rising temperatures are likely to result in a reduced “sessile stage” for oviparous species (Wheeler et al., 2020), during which they are particularly susceptible to predation (Cox and Koob, 1993). While the sustained effects of EW over viviparous species' development have yet to be investigated in a controlled setting, elasmobranch females often display habitat use patterns suggestive of maternal thermophily. Indeed, there are several reports of species with different reproductive modes exploiting environmental thermal heterogeneity during pregnancy or oviposition, presumably to decrease pregnancy energy costs and accelerate development (Wallman and Bennett, 2006; Jirik and Lowe, 2012; Salinas-de-Leon et al., 2018). However, the effects of EW over female oviposition rates and pregnancy costs for females have yet to be addressed and represents a challenging but important frontier. Another issue worth considering is the potential community-level impacts of these changes (e.g., early hatching) in resource availability for the neonates and pregnant females, particularly in regions and species with consistent seasonal and phenological cycles [match-mismatch hypothesis, first proposed by Cushing (1975); Durant et al., 2007; Anderson et al., 2013]. Meanwhile, although not significant, EA resulted in a slight trend in the opposite direction that should be elucidated with further research. Similar effects have been observed in some teleost fishes (Heuer and Grosell, 2014). However, the observed decline in embryogenesis duration in response to EA conditions appears to be population-specific (Di Santo, 2015). This highlights the variability associated with the response to EA and underscores the importance of intra-species variation in response to climate change-related drivers (Di Santo, 2015; Gervais et al., 2021).

Growth and Feeding

Shark body size is tightly linked to its diet (and consequent predation pressure on both sides of the equation; Heupel et al., 2014), and temperature has long been established to influence adult body size (temperature-size rule; Atkinson, 1994; Rosa et al., 2012). Hence, understanding how climate change may affect growth patterns has fittingly represented a research priority (Figure 2F). Based on the currently available literature on elasmobranchs, no significant effect is observed across all metrics analyzed (Figures 4D–F; Supplementary Table 3). A downward trend for body condition in response to EA and EWA should, nonetheless, be remarked. Indeed, while EW and EA have been observed to influence growth in fishes (Audzijonyte et al., 2020; Sampaio et al., 2021), results appear to be widely context-dependent for both elasmobranchs (e.g., Pistevos et al., 2015) and teleosts (Cattano et al., 2018). Moreover, most studies feature ad libitum feeding, with food availability known to modulate growth under EW and EA (Baumann, 2019; Cominassi et al., 2020; Otjacques et al., 2020). Indeed, feeding and digestion responded positively to EW, with a close trend observed for EAW (Figure 5A; Supplementary Table 3), which aligns with hypothesized increased energetic needs. Meanwhile, by having animals reared in a mesocosm, Pistevos et al. (2015) observed a decline in growth in animals exposed to EA and the disappearance of the EW effects observed in animals maintained under standard laboratory settings, likely related to increased energetic needs coupled with the costs of foraging. Moreover, biochronological reconstructions have revealed a negative correlation between growth and temperature during the breeding season in the same species (Heterodontus portusjacksoni; Izzo and Gillanders, 2020). All things considered, growth patterns in the wild are likely to respond differently from what is observed under normal laboratory conditions.

FIGURE 5

Figure 5. Directional effects size (LnR) of experimental warming (EW), acidification (EA), and their combination (EWA) over elasmobranchs, considering specific biological responses [(A) Feeding and digestion; (B–D) Aerobic metabolism; (E) Oxygen transport; (F) Anaerobic metabolism]. The dashed line (R = 0) represents the threshold used to test the null hypothesis, with overlapping confidence intervals (CI) indicating the lack of significant treatment effect. R values above 0 indicate a stimulation and below indicates an inhibition of the response; asterisks indicate a significant effect. Different letters showcase differences between treatments; numbers indicate sample sizes (number of control-treatment contrasts). Statistical results are detailed in Supplementary Table 3.

Aerobic Metabolism

Performance declines with temperature have been widely attributed to an aerobic ceiling stemming from the species' systemic capacity for oxygen delivery [oxygen- and capacity-limited thermal tolerance (OCLTT) hypothesis; Pörtner and Farrell, 2008]. Under this hypothesis, aerobic scope and the associated biological performance are presumed to increase until the optimum temperature is reached, above which both performance and oxygen uptake decline relatively rapidly, triggering the onset of unsustainable anaerobic pathways. Other environmental and physiological stressors, such as hypercapnia and hypoxia, are thought to act as limiting factors, narrowing the thermal tolerance window (Pörtner and Farrell, 2008). In this context, a relatively large number have assessed oxygen uptake rates, used as a proxy for aerobic metabolism, in elasmobranchs under climate change-related scenarios (Figure 2F; Bouyoucos et al., 2019). Still, available data is sparce. As such, given the current data limitations, here we merge some metrics that do not strictly adhere to the conventional definitions of maximum/standard metabolic rates and aerobic scope (Lefevre, 2016), yet represent relevant proxies in the context of the present work and the discussion of future avenues of research (see Supplementary Data 1; e.g., Di Santo, 2015). To safeguard this distinction, we conservatively implement a more accommodating nomenclature, namely routine metabolism (RM), active metabolism (AM), and metabolic scope (MS).

Analyzing the data obtained across all elasmobranchs, EW and EWA robustly produced an expectable increase of RM (Figure 5B; Supplementary Table 3), likely matching an increase in basal energetic demands, which is further supported by the increase in food consumption. AM, on the other hand, was also positively affected by EW, with a close trend observed for EWA—although the significance of effects observed is still quite sensitive to case deletion (Figure 5C; Supplementary Table 3; Sensitivity analysis, Supplementary Code 1). Meanwhile, MS appears to be positively affected by EWA, with a close trend also observed for EW (Figure 5D; Supplementary Table 3). The latter effect is, however, restricted to batoids and largely driven by Di Santo (2015, 2016) studies (see Sensitivity analysis, Supplementary Code 1). Regardless, while some of the differences observed in response patterns may reflect species' physiology, differences may likely also reflect methodological disparities [e.g., extensive air exposure (Schwieterman et al., 2019) or the use of anesthesia to achieve a resting state (Di Santo, 2015)] necessary due to specificities in batoid body plan and behavior. Practical limitations in protocol standardization across these two groups may be mitigated by the combined use of molecular and biochemical indicators (e.g., Tullis and Baillie, 2005; Rosa et al., 2016a). No significant effects were observed for EA across all three metrics of aerobic metabolism (Figures 5B–D; Supplementary Table 3).

Oxygen Transport

In addition to increasing oxygen demands, oxygen affinity typically decreases with temperature (Nikinmaa et al., 2019), and several elasmobranch species have been shown to adjust their hematological profile in response to environmental challenges, namely temperature (Butler and Taylor, 1975; Neale et al., 1977; Bouyoucos et al., 2020a; Pegado et al., 2020b). However, no significant overall effect was detected for EW or EA for the endpoints analyzed here, with insufficient data to conduct the analysis for EWA (Figure 5E; Supplementary Table 3). While this may suggest that no consistent directional changes in oxygen transport are made to support increasing aerobic requisites, it should be noted that different species may implement distinct strategies, with further research required. On the other hand, the erythrocyte properties of elasmobranchs may inherently support much of the induced shifts (Nikinmaa et al., 2019). Still, how the cumulative action of these factors affects this response, however, remains seldomly addressed (Bouyoucos et al., 2020c).

Anaerobic Metabolism

Although the long-term reliance on anaerobic metabolism is not sustainable, such pathways are key for ecological persistence—namely by supporting short-term survival under physiological extremes and for sustaining burst swimming, which is essential for hunting and predator evasion (Hochachka and Somero, 2002). Meanwhile, without enough studies assessing either biochemical or whole-body proxies for anaerobic metabolism, we combined both approaches to gauge the effects of EW, EA, and EWA over these pathways, with no significant effects being observed regardless of treatment (Figure 5F; Supplementary Table 3). More standardized research is recommended, however, and would benefit from the parallel use of anaerobic biochemical markers, allowing us to better understand how animals may change their reliance on different metabolic pathways as temperatures move away from their thermal optima (Pörtner, 2002).

Tolerance and Stress Indicators

Thermal tolerance showed an overall positive response to EW (Figure 6A; Supplementary Table 3). It is nonetheless worth noting that, while thermal tolerance does seem to increase, thermal safety margins are likely to decline (Bouyoucos et al., 2020a), which is particularly important in the context of gradual OW and simultaneously more pervasive MHWs. In this context, while animals may cope with OW and mild MHWs, their capacity to withstand severe or extreme MHWs may be undercut. Moreover, the effects of EA and EWA have not been addressed and should be considered in future studies. On the other hand, no significant effects of either treatment were observed over hypoxia tolerance (Figure 6B; Supplementary Table 3), despite previous indication that these two responses may be intertwined (Butler and Taylor, 1975; Ely et al., 2014; Bouyoucos et al., 2020a). Likewise, cell stress levels, particularly oxidative damage, appear not to be affected by EA (Figure 6C; Supplementary Table 3), with elasmobranchs relying on non-enzymatic ROS-scavengers in addition to enzymatic antioxidants, which may be a potential explanation for their resilience (Lopes et al., 2018; Pegado et al., 2020a). However, this same response has rarely been addressed under EW and EWA, despite more evidence suggesting potentially considerable effects of these treatments over elasmobranchs and other taxa (Lesser, 2006; Rosa et al., 2016a).

FIGURE 6

Figure 6. Directional effect sizes (LnR) of experimental warming (EW), acidification (EA), and their combination (EWA) over elasmobranchs, considering specific biological responses [(A) Thermal tolerance; (B) Hypoxia tolerance; (C) Cell stress; (D,E) Acid-base status; (F) Neuroendocrinology; (G–I) Behavior]. The dashed line (R = 0) represents the threshold used to test the null hypothesis, with overlapping confidence intervals (CI) indicating the lack of significant treatment effect. R values above 0 indicate a stimulation and below indicates an inhibition of the response; asterisks indicate a significant effect; close trends are also noted. Different letters showcase differences between treatments; numbers indicate sample sizes (number of control-treatment contrasts). Statistical results are detailed in Supplementary Table 3.

Acid-Base Status

As previously established (e.g., Evans et al., 2004; Green and Jutfelt, 2014; Heinrich et al., 2014), elasmobranchs promptly resume their acid-base balance when exposed to EA (Figure 6D; Supplementary Table 3), specifically through a competent buffering capacity (Figure 6E; Supplementary Table 3). Still, the underlying biochemical mechanism for elasmobranch is still poorly understood and the energetic consequences of this process are not well-established (Lefevre, 2016). Moreover, this process has also been hypothesized to be the catalyst for some of the behavioral changes observed in teleost fishes (GABA_A hypothesis; Nilsson et al., 2012), which remains to be tested in the context of elasmobranch behavior.

Behavior

Behavior is at the forefront of individual response to environmental stress, thus representing an important facet of animal response to changing ocean conditions (Nagelkerken and Munday, 2016). However, the diversity of behaviors analyzed, along with the lack of standardization in procedures, complicates this analysis. To provide a general overview of the effects of these factors on behavior, the different endpoints were analyzed in terms of absolute response (Figure 6G; Supplementary Table 3). While all treatments significantly affected behavior, the effect of EWA was significantly lower than its isolated counterparts. This may reflect an antagonist effect of both factors (e.g., Pistevos et al., 2017), although the specific implications of this finding are more difficult to evaluate because of the paucity of studies assessing behavioral effects with comparable approaches. In fact, despite multiple individual reports of behavioral shifts (e.g., Green and Jutfelt, 2014; Pistevos et al., 2015; Pegado et al., 2018), no net changes in swimming patterns and foraging behavior are detectable for both EA and EW (Figures 6H,I; Supplementary Table 3). Meanwhile, studies into the neurologic and endocrine underpinnings of behavioral effects under climate change are scant, although effects at multiple levels are to be expected (Servili et al., 2020). Indeed, this specific response was only possible to analyze in the context of EW, for which no significant net effects were detected (Figure 6F; Supplementary Table 3). However, it is worth noting that two of the studies included in this analysis measured the effects of EW over testosterone concentrations, both reporting significant differences although in opposite directions (Neale et al., 1977; Mull et al., 2008). This disparity reflects the context specificity of this response, highlighting the need for more research, particularly considering potential repercussions over reproduction. Indeed, no studies have specifically addressed the effects of these factors over shark reproduction, which potentially constitutes one of the most consequential knowledge gaps in experimental research and a challenging frontier.

Conclusions, Limitations, and Future Directions

One of the most evident insights from experimental research, as examined here, is that rising temperatures will likely represent the most prevalent direct driver of elasmobranch response, with a consistent effect observed across attributes (see summarized results in Figure 7). Moreover, although considerable knowledge gaps still linger, directional effects are relatively well-understood, with a clear reduction of development time and an increase in aerobic metabolism, feeding and thermal tolerance. Meanwhile, although EA did induce a considerable overall effect, the responses observed appear to be more context-specific and difficult to pin-down. In fact, there were no directional effects over specific biological responses, apart from the buffering adjustment required to re-establish acid-base balance. The effects observed for EWA mostly match those observed for EW, although different levels of responsiveness are apparent across distinct attributes and considerable knowledge gaps preclude further assessment. On the other hand, while no significative effects were observed here, sustained ED has mostly been disregarded in experimental research, although consequential effects are to be expected (Sampaio et al., 2021).

FIGURE 7

Figure 7. Summary the results obtained in the different stages of the analysis the effects of experimental warming (EW), acidification (EA), deoxygenation (ED), and the combination of warming and acidification (EWA) over elasmobranchs. Asterisks (*) indicate significant impacts over combined biological response in terms of absolute change (Stage I and II). Arrows indicate the direction of effects [stimulation (↑) or inhibition (↓)] over specific biological responses (Stage III). Equals sign (=) indicates the lack of significant effects. Empty cells identify responses not assessed due to data paucity.

While the present analysis shows how climate change-related pressures have the potential to strongly influence elasmobranch biology, whether these effects will translate into detrimental effects to populations, species, and communities, will likely depend on a wide range of ecological dependencies. For example, shifts in ocean productivity and associated changes in prey distribution may further constraint this group (Schlaff et al., 2014; Bindoff et al., 2019), particularly considering increased energetic demands. Indeed, the current knowledge gaps and the interweaved nature of biological and ecological processes in the wild, precludes a direct assessment. This uncertainty is compounded by indirect (and potentially interacting) effects of other drivers of global change, as well as by the sheer diversity within elasmobranchs. For example, the logistic constraints associated with controlled exposure in a laboratory context enfold an inherent bias toward the use of species and life stages amenable to captivity rearing. On the flip side, observations obtained from particularly robust species are likely to represent an underestimation of the general effects, representing conservative indicators (Heinrich et al., 2014; Wheeler et al., 2021). Moreover, the effects of anthropogenic CO₂ emissions extend far beyond decreasing pH, rising temperatures, and oxygen loss. Global warming is further driving changes in precipitation patterns, ocean circulation and stratification, and sea-level rise, entailing the potential disruption of key habitats for these animals (Chin et al., 2010; Bindoff et al., 2019). At the same time, OD, which is mostly driven by OW, is expected encompasses the expansion of oxygen minimum zones and the concurrent compression of the available habitat for pelagic sharks, altering ecosystem dynamics while increasing the exposure of these animals to fishing pressure (Vedor et al., 2021). Moreover, these effects are further likely to cascade into considerable ecological consequences, such as changes in productivity, shifts in food availability, deterioration of foundation species, and the potential increased prevalence of harmful algae blooms (Bindoff et al., 2019; Leggat et al., 2019; Smale et al., 2019)—all of them likely to affect elasmobranchs to an extent.

Finally, while we have emphasized important research gaps throughout this discussion (see Figure 2), some final considerations ought to be remarked. First, one crucial research target, which should provide important insights for integration with different approaches is the identification of the condition thresholds that induce relocation. Although some efforts have been made in this direction (e.g., Gervais et al., 2018; Nay et al., 2021), most experiments rely on a set-up that precludes animal behavioral regulation. This is important across multiple spatial and temporal scales, from transient shifts in the use of specific coastal habitats in response to extreme events to the general poleward movement of species. Second, while understanding elasmobranchs' current capacity to tolerate the levels of change expected in the future is key to assess their resilience, transient extreme events are likely to shape the overall system response and may represent bottlenecks for species' coping capacity (Grant, 2017). Hence, experimental designs aligning with transient extreme events represent a priority from both ecological and evolutionary perspectives. Third, the biochemical and molecular mechanisms underpinning elasmobranch response to these factors are still very poorly understood. This represents a major concern, given that this group's physiology differs widely from other vertebrate lineages (Speers-Roesch and Treberg, 2009). Likewise, knowledge gaps relating to the effects of these factors over several additional biological traits essential for fitness, namely reproduction and immune competence, should also represent a research priority. Furthermore, given the long generation times typical of elasmobranchs, the quick pace of global change, and the apparently weak selection pressure placed by these factors over this group, the potential for genetic adaption is quite slim. In this context, mechanisms such as epigenetic changes (Lighten et al., 2016) and transgenerational acclimation (Donelson et al., 2018) may represent key components to this group's coping capacity. These mechanisms may broaden elasmobranch phenotypic plasticity and provided a buffer for slower adaptation processes, meriting further investigation. Moreover, in addition to multi-stressor treatments featuring the factors explored here, studies exploring potential interactions with other global change drivers (e.g., pollutants) would be of great value. Finally, while some lineages with somewhat amenable species have seldom been investigated, research featuring certain lineages and stages is limited (or entirely out of reach) due to logistical, ethical, and conservation constraints. In these cases, the use of alternative approaches (e.g., in vitro studies, telemetry, and environmental modeling) becomes increasingly important. Integrating different approaches is therefore key to robustly assess elasmobranchs' vulnerability to climate change—and to the many other challenges of the Anthropocene.

As one of the oldest lineages of vertebrates, sharks and their relatives have experienced climate changes many times over, however, not without important consequences (O'Brien et al., 2013; Pimiento et al., 2017). This time around, not only are conditions shifting at a remarkably quick pace (Bindoff et al., 2019), but these challenges are compounded by a much more dire type of pressure for k-strategists—overexploitation (Dulvy et al., 2014; Queiroz et al., 2019; Pacoureau et al., 2021). With the clock ticking for the preemptive design and implementation of strategies directed toward the mitigation of climate change impacts over this threatened group, it is crucial that knowledge gaps are bridged and that multiple lines of research are developed, integrated, and consolidated.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

CS, RR, and ES conceptualized the article. CS, BP, and MP retrieved the data. CS, FB, and CW reviewed the dataset for accuracy. CS and ES analyzed the data. CS produced the first draft of the manuscript. RR, CF, JR, and IB provided supervision and editorial comments. All authors reviewed and contributed to the submitted version.

Funding

This work was supported by the Fundação para a Ciência e Tecnologia (FCT), through the strategic project granted to MARE strategic project (UID/MAR/04292/2019), the PhD grants attributed to CS (SFRH/BD/117890/2016), ES (SFRH/BD/131771/2017), and FB (SFRH/BD/147294/2019), financed by national and community funds from FCT and the European Social Fund (ESF), through the Human Capital Operating Programme and Regional Operation Programme (Lisboa 2020). BP acknowledges funding in the context of the MAR2020 programme through the project MAR-01.04.02-FEAMP-0006. IB was supported by NSERC/Manitoba Hydro Industrial Research Chair grant (479350/1) and an NSERC Discovery grant (05328). JR was supported by funding from the Australian Research Council Center of Excellence for Coral Reef Studies. CW was supported by an American Australian Association Graduate Scholarship and a James Cook University Postgraduate Research Scholarship.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.735377/full#supplementary-material

References

Abram, N., Gattuso, J.-P., Prakash, A., Cheng, L., Chidichimo, M. P., Crate, S., et al. (2019). “Framing and context of the report - intergovernmental panel on climate change special report on oceans and cryosphere in a changing climate,” in IPCC Special Report on the Ocean and Cryosphere in a Changing Climate, eds H.-O. Pörtner, D.C. Roberts, V. Masson-Delmotte, P. Zhai, M. Tignor, E. Poloczanska, K. Mintenbeck, A. Alegría, M. Nicolai, A. Okem, J. Petzold, B. Rama, and N. M. Weyer.

Anderson, J. J., Gurarie, E., Bracis, C., Burke, B. J., and Laidre, K. L. (2013). Modeling climate change impacts on phenology and population dynamics of migratory marine species. Ecol. Model. 264, 83–97. doi: 10.1016/j.ecolmodel.2013.03.009