Multiple stressors, allostasis and metabolic scaling in developing zebrafish

A recent report from the Intergovernmental Panel on Climate Change (IPCC) on the ocean and the cryosphere highlighted the alarming increase of deoxygenation, warming and acidification in aquatic ecosystems (IPCC, 2019). Although some aquatic organisms may be able to respond behaviorally to these changes, most will have to physiologically adjust to reach energetic stability through change, a process known as allostasis (Sterling and Eyer, 1988; Sterling, 2012). We use allostasis instead of homeostasis because, as Sterling (2012) stated, ‘the goal of regulation is not to preserve constancy of the internal milieu…but to continually adjust the milieu to promote survival’. Adult and juvenile stages can make allostatic adjustments when exposed to stressors by incurring minimal energetic costs or allostatic loads (McEwen and Wingfield, 2003, 2007; Sterling and Eyer, 1988). Embryos and larvae, however, face greater challenges because of inherently higher growth rates, smaller aerobic scopes (Rombough, 1988, 2006), elevated metabolic costs (including transport costs; Kauffman, 1990), tight energy budgets and limited food supplies (yolk) (Kamler, 1992; Rombough, 2006). This makes early stages less energetically plastic and more vulnerable to stressors than older stages and more likely to suffer unsustainable cumulative energetic costs or ‘wear and tear’ (Ellis and Del Guidice, 2014).

From past studies, we know that fish larval growth rates typically range from 25–30% day⁻¹ to 100% day⁻¹ (Kamler, 1992; Wieser, 1994; Conceiçaõ et al., 1997; Pedersen, 1997), while juveniles and adults grow at much slower rates ranging from 1 to 3% day⁻¹. In a suite of studies in the 1980s and 1990s, Wolfgang Wieser and colleagues determined that high larval growth rates were possible because of compensatory energy budgeting in which the limited energy available from yolk is reduced (or traded off) from one physiological function (e.g. metabolic activities) and redirected to another (e.g. growth) (Wieser, 1994, 1989; Rombough, 2006; see Pan et al., 2015, 2018 and Griffiths et al., 2021 for metabolic allocation models developed for marine invertebrates). This energetic strategy may provide early stages with the potential to attain novel stable physiological or allostatic states when exposed to stressors (Mauss et al., 2015; McEwen and Wingfield, 2003; Glazier, 2005) and thus improve survival. This strategy differs from that of juveniles and adults in which energy budgets are additive and additional energy needs (e.g. such as those brought on by stressors) can be obtained through processes such as behavioral modification and/or feeding. Although balancing the energy budget during any life stage is difficult, acclimation in the face of external stressors becomes even more challenging in early life stages when development (i.e. rates of cellular differentiation) and growth (increases in somatic mass) are rapid, and must occur in concert with maintaining energetic stability or allostasis.

Multi-stressor experiments are important, because interactions may be additive, synergistic, multiplicative or antagonistic, but given the technical difficulties of conducting experiments on the effects of multiple stressors on fish early life stages, most past studies have selected to examine only the effects of a single stressor (e.g. high temperature or oxygen concentration). Only a few studies have investigated the effects of multiple stressors on metabolic rates (Garside, 1966; Hamor and Garside, 1976; Rombough, 1994; Conceiçaõ et al., 1997) and results were variable, ranging from stressors acting in an additive, synergistic, multiplicative or antagonistic manner (see reviews in Crain et al., 2008; Tekin et al., 2020). Therefore, to further understand how multiple stressors act on developing fish larvae, we conducted a sister study (Pan and Hunt von Herbing, 2017) where we examined the chronic and acute effects of two stressors (temperature and P_O2) on larval morphological development and metabolic thermal sensitivity (Q₁₀ values). As in the present study, short-term or acute stress was defined as a transient change in metabolism in response to stressors (Xu et al., 2006; McEwen and Wingfield, 2003), and long-term or chronic stress was defined as exposure to either multiple or constant stressors (Santos et al., 2010). Results from Pan and Hunt von Herbing (2017) showed that although development and morphology were minimally impacted, metabolic compensation (a synergistic response) was clearly evident at the yolk-sac larval stage. However, it was unclear what effect these two stressors would have on the energetic costs of acclimation (allostasis) or on metabolic-mass scaling relationships. Therefore, the present study investigated both chronic and acute effects of multiple stressors (temperature and P_O2) on % net growth efficiency (K_N), net costs of growth (C_r) and metabolic scaling during the earliest period of larval development (e.g. the yolk-sac period) to better determine how physiologically plastic this early stage might be to rapidly changing aquatic environments.

Metabolic–mass scaling relationships are most commonly expressed as those in which body mass (M) affects metabolic rate (R) according to the power function R=aM^b, where metabolism and mass are expressed in log–log space, R is standard (R_s) or routine (R_r) metabolic rate, a is the scaling coefficient or mass-independent constant (antilog of the intercept in the linear log–log plot) and b is the scaling exponent (slope of the log–log plot). Although the value of the scaling exponent b has been actively debated, it is assumed to be fixed at a value of around 0.75 (explained by variation of fractal nutrient transport networks; West et al., 1997) across all organisms, representing a universal scaling ‘law’ (Hemmingsen, 1960; Savage et al., 2004), and underlies the metabolic theory of ecology (Gillooly et al., 2001; Gillooly and Allen, 2007). However, other studies have shown that the relationship between body mass and metabolic rate can change through ontogeny and does not follow a power function (Moran and Wells, 2007; Kolokotrones et al., 2010; Seymour et al., 2013; White et al., 2022). Further, an alternative to the metabolic theory of ecology called the ‘metabolic level boundaries hypothesis’ (MLBH) was suggested by Glazier (2005, 2008, 2009), in which b ranged between 0.67 (scaling limited by surface-area limits on the fluxes of resources, wastes and heat) and 1.0 (mass or volume limits on energy use or power production), such that scaling can often vary with changes in ambient temperature (Killen et al., 2010; Tan et al., 2019) and lifestyle (e.g. mobility) or thermoregulatory mode (Glazier, 2010). Although b has received much attention because of its potential predictive power, the scaling coefficient a, or metabolic level (L), has received less as it is considered to be highly variable across taxa (White et al., 2006; Seibel, 2007; Glazier, 2009) and substantial covariation has been found to exist in b with elevation (a, L) (Glazier, 2008, 2009, 2010; West and West, 2013). Moreover, most scaling studies have been conducted on juvenile or adult stages because obtaining reliable scaling relationships under optimal conditions in rapidly developing and growing fish larvae is more difficult. As a result, metabolic scaling exponents (b) of early stages, although more variable than adults, generally ranged from 0.75 to 1.0 (Giguère et al., 1988; Hunt von Herbing and Boutilier, 1996; Hunt von Herbing, 2006; Kauffman, 1990; Wieser, 1989; Rombough, 1988, 1994; Post and Lee, 1996; Bochdansky and Leggett, 2001). Further, unlike in adults, a and L, and the relationships between scaling variables (e.g. between b and a or L), remain unexamined. For life stages in which energy budgets are tight, such as in fish larvae, examining scaling relationships under stress is therefore valuable because added energetic costs from stressors may contribute to better informed energetic or allostatic models (Glazier, 2005; Finn et al., 1995; Post and Lee, 1996; Glencross, 2008; Killen et al., 2010; Glencross and Bermudas, 2011, 2012).

Although allostatic models have helped explain physiological stress responses in mammals (Korte et al., 2007; McEwen and Stellar, 1993; McEwen and Wingfield, 2003; Sterling, 2012; Sterling and Eyer, 1988), including humans (Ellis and Del Guidice, 2014; Edes and Crews, 2016), they only recently gained acceptance for use in adult and juvenile fish stress physiology (Glencross and Bermudes, 2012; Ogoshi et al., 2012; Schreck and Tort, 2016; Samaras et al., 2018; Oldham et al., 2019; Leeuwis et al., 2019). The concept of allostasis provides a way to understand how developing embryos and larvae may differentially allocate limited energy from yolk to growth and metabolic activities in response to stressors, but it has not yet been applied to determine how multiple stressors could alter metabolic–mass scaling relationships. We build on our earlier zebrafish (Danio rerio) study (Pan and Hunt von Herbing, 2017) in which the unexpected finding that metabolic thermal sensitivity (Q₁₀) remained at around 2.2, when larvae were exposed to a 3°C increase in rearing temperature (28 to 31°C). In the present study, we expected that zebrafish larvae might be energetically tolerant (exhibit allostasis) to combinations of high temperature (31°C) (which drive metabolic demand) and hypoxia (11 kPa) (which drive supply; see Rombough, 2006) because zebrafish in the wild lay demersal eggs in shallow tropical ponds and in slow-moving streams, where temperatures are likely to be high and O₂ levels low (Spence et al., 2008). Therefore, the major goal was to determine how larvae adjusted their energetic partitioning including the metabolic cost of growth (C_r) and % net growth efficiency (K_N), when they experienced chronic and acute combinations of high temperature and/or hypoxia. Given the much higher mass-specific daily growth rates and tight energy budgets generally found in fish yolk-sac larvae (Kamler, 1992), we hypothesized that changes in energetic partitioning (C_r and % K_N) would depend upon whether these two stressors acted synergistically and/or in an additive manner (see Pan and Hunt von Herbing, 2017). If synergistic or additive, enhanced thermal responses by hypoxia on energetic costs or allostatic loads would result in increased C_r and decreased % K_N, possibly altering scaling relationships between body mass and oxygen consumption rate (Ṁ_O2) because of cumulative energetic costs.

Danio rerio (Hamilton 1822) is especially suited for studies that examine multiple stressors because this species appears to operate close to the upper limit in terms of energy expenditure (Rombough, 2006), produces small larvae (<200 µg total mass at hatch), with one of the highest routine Ṁ_O2 reported for fish larvae (∼100 mol O₂ g⁻¹ h⁻¹; Pelster and Burggren, 1996), and optimally develops at high temperatures between 27 and 29°C. Adult zebrafish (D. rerio; 0.6–0.9 g) were an outbred stock purchased from local commercial suppliers in Denton, TX, USA. Before breeding, male and female fish were kept separately in a semi-circulating system under a 14 h:10 h light:dark cycle at 28±1°C and were fed twice daily with commercial flake food (Tetra, VA, USA) as well as fresh brine shrimp. A total of 13 cohorts over a period of 6 months obtained from eggs produced by several breeding pairs were used. Embryos from each adult pairing were collected from the bottom of the breeding tanks within 2 h post-fertilization (hpf) (1–8 cell stage), divided into four groups (∼200 eggs each) and transferred to 20-l experimental tanks, each containing water pre-equilibrated to one of the four chronic treatments (see experiment 1 below). Larvae were exposed to the same photoperiod as the adults. After hatching at ∼48 hpf, larvae relied only on yolk until ∼144 hpf, after which the anus opened, and while exogenous feeding was possible, it was not observed and is usually poor at this early transition stage (Hunt von Herbing et al., 1996). All measurements were conducted in the first week after hatching (48 to 168 hpf).

Two different experiments were conducted: experiment 1 – chronic exposure, and experiment 2 – acute exposure. Both experiments used a 2×2 factorial combination of two temperatures, control (28°C) and high temperature (31°C), and two oxygen concentrations, normoxia (21 kPa) and hypoxia (11 kPa) (Table 1). Experiment 1 was designed to determine effects of chronic exposure to the four treatment conditions on the use of yolk energy for larval production of somatic tissue mass (G), routine metabolic rate (R_r), % net growth efficiency (K_N), mass-specific cost of growth (C_r), and the scaling relationships between R_r and dry mass (M_D) over the first week after hatching (48–168 hpf). In experiment 2, embryos and larvae raised under the four different chronic conditions (see experiment 1) were further exposed to acute stressors of high temperature (31°C) and hypoxia (11 kPa) either singly or in combination, to estimate larval metabolic resilience to acute stress following exposure to chronic stress.

Embryos were raised at the combinations of temperatures and oxygen concentrations as follows: (1) control temperature–normoxia (T28O21); (2) control temperature–hypoxia (T28O11); (3) high temperature–normoxia (T31O21); and (4) high temperature–hypoxia (T31O11) (Table 1). Chronic treatments were repeated using different cohorts. The selected oxygen concentration of ∼11 kPa (∼50% reduction in oxygen concentration relative to control) is slightly below the critical partial pressure of oxygen (P_O2 at which organismal O₂ consumption becomes limited) for zebrafish embryos (Barrionuevo and Burggren, 1999), and is physiologically stressful during embryonic development. Hypoxia was generated by directly bubbling N₂ into the water and water temperature and P_O2 levels were monitored daily using a hand-held optical dissolved oxygen and temperature meter (Model ProODO, YSI Inc., Yellow Springs, OH, USA). Metabolic rates under chronic conditions were measured at the temperature and P_O2 under the four (see above) larval rearing conditions. A total of 66 measurements were made to investigate effects of chronic exposure on larval metabolic rates during the first week of development. Additionally, a subset of larvae was used to measure larval growth and yolk utilization.

Embryos were raised (0–48 hpf) until hatching under the specific combinations of temperature and hypoxia listed in experiment 1 (i.e. T28O21, T31O21, T28O11, T31O11). Rates of oxygen consumption (routine metabolic rate, R_r) were first measured at the embryonic rearing conditions and then larvae were exposed to acute changes in temperature and P_O2 (either singly or in combination) in a respirometer to evaluate the effects of added acute stressors on larval metabolic rates. For each acute measurement, larvae from one of the four chronic conditions (T28O21, T31O21, T28O11 and T31O11) were transferred to respirometry chambers in which water was first equilibrated to the combinations of temperature and P_O2 different from the chronic conditions, and larvae were exposed to these acute conditions for 1 h after initiation, during which a series of measurements were taken (see respirometry section for details). As there were a high number of combinations of chronic and acute exposure treatments (four chronic treatment, each measured at two temperatures and two P_O2 levels), Ṁ_O2 under different environmental conditions was measured at different times throughout the first week of development (48–168 hpf). In total, 174 measurements were made to investigate effects of chronic versus acute exposure on larval metabolic rates. For the calculation of energetic budgets, measurements of oxygen consumption rates at different developmental time points were pooled for each chronic treatment group.

To determine the effects of chronic exposure on the use of yolk energy between growth (G) and metabolism (R_r), larvae (n=7–10) were sampled for changes in size every 24 h from 48 to 168 hpf and preserved in 10% buffered formalin. Determination of the changes in growth and yolk-sac utilization included measurement of the following variables; dry mass (M_D), standard length (SL, measured from the tip of the snout to the posterior end of the last vertebra), surface area of the yolk-sac (Y_A, as measured by the largest optical section across the yolk). Images were recorded for each preserved larva from the left side of the body using a Zeiss Axis Observeor inverted microscope (model: Axio Observer.Z1, Carl Zeiss, Oberkochen, Germany) equipped with an AxioCam HRm digital camera mounted on the microscope. All measurements were made on calibrated images using AxioVision Software (Version 4.8.2, Carl Zeiss). M_D was determined to the nearest 0.01 mg using a Cahn microbalance (precision ±1.0 μg) (Thermo Electron Corp., MA, USA), after drying larvae at 65°C overnight in a drying oven. In addition to the samples (n=7–10 at each time point) collected for the morphometric measurements, M_D measurements were also made on larvae used in respirometry, and those measurements were included in the analyses.

AGR (dM_D/dt) was calculated as the first derivative of Eqn 4, dM_D/dt=(ab)e^−bt. % SGR was calculated as dM_Df/M_Didt×100%, where dM_Df is final dry mass and M_Di is initial dry mass.

Measurements of changes in Ṁ_O2 during chronic and acute exposures were made using a Loligo microrespirometer system (Loligo Systems, Viborg, Denmark). For each measurement, three larvae from one of the four rearing conditions were placed in a 2-ml glass respirometer chamber. Each respirometer chamber was placed on a mini magnetic stirrer to allow proper mixing of the water in the chamber. After a 30-min acclimation period, P_O2 in the chamber was continuously monitored using a fiber-optic O₂ sensor connected to an OXY-4 O₂ meter (PreSens, Germany) for up to 1 h after initiation. Control measurements of Ṁ_O2 were performed without an animal in the chamber and these blanks were subtracted from every animal's measurement of Ṁ_O2. P_O2 was recorded by an automated data acquisition system (Model DAQ-M, Loligo Systems), and analyzed using AutoResp Software. Decline of oxygen over time was linear in all assays with larvae, with a regression of r²>0.90. Ṁ_O2 was calculated from P_O2 between the start and end of each measurement period.

Two-factor ANOVAs and Bonferroni multiple comparisons tests were used to determine the effects of chronic treatments on growth variables: larval dry mass at hatching (M_Dh) and final dry mass at the end of the experiment (M_Df), standard length at hatching (SL_h) and final standard length (SL_f) as well as yolk utilization rates measured as changes in yolk-sac area (Y_A) over time. A significance level of P<0.05 was used. For metabolic–mass scaling relationships, least-squares regression (LSR) was used because LSR is more appropriate when the x variable is thought to be influencing the y variable, rather than the reverse (Smith, 1997), and when y is associated with significantly more measurement error than x, which is the case for all variables (M_D, L_S and Y_A) as well as oxygen consumption rates (Ṁ_O2) (after Glazier, 2010; White, 2011). log₁₀ values were used to normalize the data variation and to produce linear, proportional relationships that could be compared (see Kerkhoff and Enquist, 2009; White, 2011; Glazier et al., 2011). ANCOVAs (using body mass as a covariate) determined the significance of differences among b, scaling exponents (slopes) and a, coefficients or y-intercepts (Systat Software, San Jose, CA, USA, Version 13) by comparing these parameters using 95% confidence intervals (CI). A mean value outside the 95% CI of another mean value was considered significantly different (P<0.05) (see Smith, 1997). Complete non-overlap of the 95% CI of two means indicated a highly significant (P<0.01) difference (Belia et al., 2005; Cumming, 2008).

The log–log scaling slope (b) was plotted against two complementary values of metabolic level (L): the intercept (log a) of the scaling relationship (L_a) and the mass-specific oxygen consumption rate at the midpoint of the scaling relationship (L_m) (after Glazier, 2009, 2010). Two different values of L were used because the standard estimate of the elevation of the regression line L_a may be partially autocorrelated with the slope (b). This autocorrelation effect is greatest when the scaling relationships are close and their intercepts are far from the midpoint of the sample body mass range (see Glazier, 2009, 2010). Intercepts (log a) at the low end of the body mass range will scale negatively with b, while intercepts at the upper end will scale positively with b (Peters, 1983; Glazier, 2009). Using L_m as well as L_a helped address the autocorrelation problem because L_m is the pivotal midpoint of the regression line and thus independent of the slope (b), mathematically (see Glazier, 2009). If the relationships of L_a and L_m with b are similar, then the effect of L on b is likely to be biologically significant and not just a mathematical or statistical artifact (see also results and discussion of Glazier, 2009 and the present study).

All scaling analyses were based on sample sizes of metabolic rates (n>15) and growth (n>20), reducing the possibility that treatment differences were due to measurement error, given the short total time of the experimental period (48–68 hpf). All values are shown as means±95% CI.

Temperature had a significant effect on M_Dh, which in the high temperature (T31O21) treatment resulted an 11% increase in dry mass and a 9% increase in standard length (M_Dh, F_3,52=9.16, P<0.05; SL_h, F_3,52=9.16, P<0.05) compared with control treatments (T28O21), while P_O2 had no effect (P>0.05) and there were no interactions between temperature and P_O2 (Table 2, Table S1). From hatching to the end of yolk-sac stage (168 hpf) M_D and standard length (SL) increased in a hyperbolic pattern in all treatments (P<0.0001) (Fig. 1A,B). Peak growth rates occurred in all treatments at 72 hpf (Table 2, Fig. 1A). High temperature reduced growth rates (AGR and SGR) by 55–57% at 72 hpf compared with controls, whereas hypoxia had little effect on growth rates (AGR and SGR), yet when combined with high temperature exhibited the lowest growth rates among all the treatments (Table 2). The domination of temperature effects on M_Dh, SL_h and SGR were also evident in the combined high temperature (T31O11 and T31O21) versus combined low temperature (T28O11 and T28O21) values (Table 2). At the end of the experiment (168 hpf), P_O2, not temperature (P>0.05), had a significant effect on final dry mass (M_Df) (F_3,79=10.144, P<0.05) (Table 2), and there was an interaction between temperature and P_O2 (F_3,79=5.433, P<0.05) (Table S1), but there was no effect of either temperature or P_O2 on final larval standard length (SL_f, P>0.05) (Table S1).

At hatching, yolk surface area (Y_Ah) was significantly different among treatments (F_3,30=4.56, P<0.05), with the smallest Y_Ah occurring at high temperature (31°C) (P<0.05) (Fig. 2, Table S1). For all treatments, yolk surface area (Y_A) declined exponentially (P<0.05) from hatching to 168 hpf (Fig. 2), but rates of yolk utilization differed among temperatures (F_3,34=8.35, P<0.05). By 96 hpf, only 11% of yolk remained at high temperature (31°C), compared with 40.8% at the control temperature (28°C). For hypoxia and high temperature/hypoxia treatments, 38.4% and 51.9% yolk remained at 96 hpf, respectively. Overall, the lowest yolk utilization rates occurred in the combined high temperature/hypoxia treatment (T31O11) in which 2.6% yolk remained at 120 hpf (Fig. 2), while there was no visible yolk in the other treatments.

Based on Eqns 1 and 2 (see Materials and Methods), % net efficiency of growth (K_N) was highest at the control temperature of 28°C, when 46% of the total assimilated yolk energy (C_y, see Eqn 1) was directed to growth (Table 3). At high temperature, K_N was reduced compared with controls and only 20.3% was directed to growth, while hypoxia did not appear to be synergistic with high temperature and had little effect on the K_N of high-temperature larvae (Table 3). Mean net mass-specific cost of growth (C_r; nmol O₂ μg⁻¹) was calculated from Eqn 3 (see Materials and Methods) for each chronic treatment over the endogenous feeding period of 48–168 hpf. Chronic exposure to high temperature increased C_r by 3.8 times compared with control temperature, and as a consequence, the final tissue dry mass (M_Df) that formed more than doubled (2.7×) at the lower temperature. Chronic exposure to hypoxia decreased C_r by 50% compared with control temperatures, and by 63% combined with high temperature, correspondingly reducing total tissue (M_Df) by 56% and 75% at 28°C and 31°C, respectively (Table 3). Thus, hypoxia had less of an effect on K_N and C_r at higher temperatures.

Yolk-sac larvae raised under control conditions (T28O21) showed a significant scaling relationship between log M_D and Ṁ_O2. The relationship between metabolism and mass was negative and had an allometric (b<1) scaling exponent (b) of 0.83±0.66 μmol O₂ h⁻¹ and an intercept (a) of 2.11±0.94 μmol O₂ h⁻¹ (±95% CI; Table 4) and was not significantly different from 0.75 or 1.0 (P>0.05). For the rest of the chronic treatments (T31O21, T28O11 and T31O11), metabolic scaling slopes (b) were not significantly different from zero (P>0.05; Table 3) and Ṁ_O2 was independent of dry mass (M_D) (Fig. 3, Table 4).

With the exception of two cases, i.e. high-temperature treated larvae exposed to acute hypoxia and hypoxia-treated larvae exposed to acute high temperature, scaling exponents (b) in metabolic–mass scaling relationships for acute exposures of larvae raised under a combination of chronic high temperature and/or hypoxia were not significantly different from zero (Table 4). For all larvae raised under control conditions and acutely exposed to any one of the following, high temperature, hypoxia or a combination of high temperature and hypoxia, metabolic–mass scaling relationships were significant (Table 4, Fig. 4A–D). When scaling exponents or slopes (b) were compared with the control slope of b=0.83±0.68, there were no differences (ANCOVA, F_3,77=1.39, P=0.25) (Table 4, Fig. 4A) and the combined b (1.19±0.25) was different from b=0.69 and 0.75, but not isometry (b=1.0) (testing for adherence to the 3/4 rule or unity). In addition, scaling exponents of control larvae acutely exposed to high temperature (b=0.81±0.77) or for control larvae acutely exposed to combined high temperature and hypoxia (b=1.25±1.08) were not different from either b=0.75 or b=1.0. But for control larvae acutely exposed to hypoxia, the scaling exponent b=1.88±1.08 was significantly higher than b=0.75, but not b=1.0.

In contrast to scaling slopes (b), the scaling coefficients (a) or y-intercepts (log a) from all the acute experiments differed when compared with scaling coefficients of larvae raised and tested under control conditions (T28O21) (log a=2.11±0.97; F_3,77=20.84, P<0.0001). There were significant differences between control larvae and larvae acutely exposed to hypoxia (log a=3.35±1.52; t=2.798, P<0.01) and between control larvae and larvae acutely exposed to combined high temperature/hypoxia (log a=2.60±1.55; t=2.881, P<0.005). However, there was no difference between larvae acutely exposed to high temperature (log a=2.17±1.10) and controls (P>0.05) (Table 4).

Finally, following Glazier (2008, 2009) and West and West (2013), we plotted scaling exponents (b) versus log metabolic level (L_a) or allometric coefficients (log a) and included two additional significant scaling relationships – from high-temperature larvae acutely exposed to hypoxia and hypoxic larvae acutely exposed to high temperature (Table 3) – which resulted in a significant positive linear relationship (r²=0.996, P<0.05) with a mean slope of 0.73 and a mean log a of −0.70 (Fig. 5A). To test for possible autocorrelation, we also plotted the mid-point (0.04 mg; based on the log dry mass that ranged from −1.2 to −1.6 mg) against b, which resulted in a negative allometric relationship (r²=0.92, P<0.05) with a mean b of −3.53 and a mean L_m of 4.34 (Fig. 5B). The two additional data points resulted in their own negative relationship (Fig. 5B). Thus, the relationships between b and L were not the same when a second estimate of metabolic level was used (i.e. the mid-point or mass-specific metabolic rate of the regression line), and autocorrelation may account for relationship between b and L, at least in these experiments.

After hatching and during the first week of life, yolk fuels growth and metabolic activities and larvae act as closed systems with respect to energy input. Results from the present study showed that zebrafish embryos raised from fertilization to the end of the yolk-sac stage (168 hpf) under two stressors (high temperature and hypoxia) experienced imbalances in their energy budgets, because of redirection of yolk energy from growth to metabolic activities in the face of higher metabolic demands of high temperature. For chronic treatments, although high temperature resulted in larger (M_Dh, SL_h) larvae at hatch, this size advantage disappeared by the end of the yolk-sac larval period because of reduced SGR, elevated net costs of growth (C_r; by almost 3× compared with controls) and decreased % net growth efficiency (% K_N), resulting in lower total mass (M_TD) at the end of the yolk-sac period. But, a synergistic response was observed as hypoxia suppressed the effects of high temperature (3°C above controls) on larval metabolism, as previously reported for Q₁₀ from our sister study (Pan and Hunt von Herbing, 2017). Therefore, yolk-sac larvae can serve as instructive models for understanding partitioning of limited energy reserves when exposed to multiple stressors (see also Crain et al., 2008; Darling and Côté, 2008; Côté et al., 2016; Tekin et al., 2020).

Imbalances in energy budgets for larvae raised under multiple stressors resulted in disruptions in metabolic scaling, because there were no significant metabolic–mass relationships (i.e. the scaling exponent b was not significant) in any chronic treatments of high temperature (i.e. 3°C higher than optimal) or the hypoxia treatments. Of the 16 different treatments (with the exception of two treatments, i.e. larvae raised at T31O21 and tested at T28O11, and larvae raised at T28O11 and tested at T31O21; see Table 4), only control larvae (T28O21) showed significant metabolic–mass scaling relationships. Thus, it is likely that yolk-sac larvae may be unable to physiologically adjust or reach energetic stability (allostasis) when raised under conditions other than those considered to be optimal (i.e. T28O21).

We considered that the reason for the large number of insignificant scaling relationships might have been because of the narrow range of body sizes tested (M_D range=25.0–50.0 μg, SL range=2.5–4.0 mm). Thus, given that in some of the treatments, scaling exponents had low precision [i.e. the confidence intervals often included zero and a range of other values, e.g. for hypoxic larvae (T28O11) exposed to high temperature (T31O11), b=−0.11±0.83; for hypoxic, high-temperature larvae (T31O11) exposed to normoxia (T31O21), b=−0.23±1.56], future work should expand the mass range as well as sample size to improve precision and reduce uncertainty. But although size ranges and sample sizes were small, scaling relationships were still significant and allometric for all control treatments (T28O21), including those in which controls were acutely exposed for 1 h in the respirometer to the two stressors. For these acute treatments, it is interesting to note that although scaling exponents (b) did not differ among treatments, scaling coefficients [log a, or metabolic level (L_a)] of acute high temperature/hypoxia treatments (2.60±1.55) were significantly lower than those of hypoxia alone (3.35±1.52), but higher than either high temperature (2.17±1.10) and controls (2.11±0.97), providing further support for synergistic or additive interactions. These physiological responses to stressors at the yolk-sac larval stage and the ability for resilience confirm results from previous studies (Schnurr et al., 2014; Beaman et al., 2016; Best et al., 2018) that the ability to regulate energetic costs through change (allostasis) exists much earlier in development than previously thought. Moreover, some of the plasticity to stressors observed in older fish stages may be inherited from the integration of genomic and epigenetic responses occurring early in development (Schnurr et al., 2014; Beaman et al., 2016; Best et al., 2018).

Detailed studies of the effects of multiple stressors on developing fish in which morphological and physiological markers change rapidly are difficult to interpret because effects of a single variable, such as temperature, on developmental and physiological timing are not linear, as was reported for Atlantic cod (Gadus morhua) and haddock (Melanogrammus aeglefinus) (Hunt von Herbing et al., 1996; Martell et al., 2005). In the present study for D. rerio, in embryos and yolk-sac larvae raised from embryos under two stressors (high temperature and hypoxia) there were no detectable differences in timing of developmental stages at hatch (48 hpf) [see our previous study (Pan and Hunt von Herbing, 2017) in which staging was based on Nüsslein-Volhard and Dahm (2002)], or on the morphometric measurements of myotomal height, body surface area and fin surface area (Pan and Hunt von Herbing, 2017). In an earlier study, Schmidt and Starck (2010) concluded that zebrafish embryos responded with an unexpected flexibility in developmental staging, as well as body mass and standard length, when exposed to multiple abiotic factors (e.g. temperature, oxygen concentration, pH and salinity). The results of Schmidt and Starck (2010) support our results that the observed developmental and physiological plasticity of zebrafish may be indicative of internal developmental mechanisms, which can foster evolutionary robustness or canalization and can result in resilience in developmental stages (Schmidt and Starck, 2010; Waddington, 1942).

The importance of multiple stressors on the timing of developmental and physiological markers as well as on the rates of energy (yolk-sac) utilization is mentioned here because of their implications on the evolution of change in developmental timing, referred to as heterochronic change (Bonner, 1988; Kuratani, 2009). By exposing embryos and yolk-sac larvae to stressors of high temperature and hypoxia, developmental modularization may act as a phenotypic buffering mechanism (i.e. conferring resilience) against extreme developmental deviation that has previously been observed as a function of abiotic change (Lukina, 1973; Kaushik et al., 1982; Schmidt and Starck, 2010; Ho and Burggren, 2012; Pan and Hunt von Herbing, 2017). In the present study, yolk-sac larvae raised under control conditions (T28O21) transformed approximately 45% (% K_N) of the energy from yolk into somatic tissue to fuel very high SGRs of 31% day⁻¹ (at 72 hpf). These values are close to the % K_N of 44% (Kamler, 1992) and SGR of 25–30% day⁻¹ (Rombough, 2006) recorded for most larvae of oviparous fishes, and typical of the endogenous feeding period in which no energy is lost as feces and active metabolism is low. In D. rerio, the mouth and anus do not open until 72 and 144 hpf, respectively (Nüsslein-Volhard and Dahm, 2002), energy is derived entirely from yolk in the first days after hatching, and the effects of stressors on the yolk-utilization rates after hatching can be estimated without considering exogenous feeding. Even after hatching, exogenous feeding is at best inefficient because, in addition to an incomplete digestive tract, more energy is spent in chasing prey than can be gained (Hunt von Herbing et al., 2001). For yolk-sac larvae from embryos raised under chronic conditions of combined high temperature and hypoxia (T31O11), the effects of high temperature dominated after hatching, and by 72 hpf, an increase of 3°C above the control temperature (28°C) resulted in decreased efficiencies (% K_N), faster yolk utilization rates and a 60% decrease in SGR because of higher growth costs (C_r). By the end of the yolk-sac stage (168 hpf), SGR in these yolk-sac larvae was 35% lower than that of controls. Based on the results from the present study and previous studies, D. rerio yolk-sac larvae may be more physiologically plastic in adjusting to increased metabolic demand (high-temperature) than to decreased metabolic capacity (hypoxia) (see also Pan and Hunt von Herbing, 2017).

Although fish yolk-sac larvae grow at much higher rates than juveniles and adults (1–3% day⁻¹) (Houlihan, 1991), larvae have intrinsically smaller aerobic scopes and higher C_r, resulting in tighter energy budgets compared with older stages. Thus, an increase in metabolic rates such as in response to an increase in temperature would result in a negative energy balance and a reduction in growth rates by the redirection of energy away from growth to metabolic activities to maintain homeostasis or to establish new homeostatic (Glazier, 2005) or allostatic states (Sterling, 2012). In the present study, both AGR and SGR were 50% lower at high temperature compared with control temperature and hypoxia. Decreases in growth rates are likely a result of almost equal reductions in net growth efficiency or % K_N, especially notable under high temperature owing to elevated metabolic rates. This rebalancing within the energy budget is possible for fish larvae under stress because there is good evidence that fish larvae exhibit compensatory energy allocation (Wieser, 1989, 1994; Rombough, 2006), unlike older stages (juveniles and adults), which operate using additive energy budgets. In compensatory partitioning, total metabolic rate remains constant even though activity levels change, and energetic costs may be met by an increase or decrease of the amount of energy devoted to another activity (Wieser, 1994; Rombough, 2006). This energetic rebalancing may only work to a point, after which energy budgets will reach an energetic imbalance. In the present study, this imbalance was observed when embryos and larvae, raised under chronic stressor conditions (e.g. high temperature or hypoxia), exhibited a reduced ability to adjust to further acute exposures, resulting in the breakdown of physiological relationships (see section on metabolic scaling). Energetic imbalances to this extent, or duration, produce unsustainable relationships between growth and metabolism, making allostasis difficult, if not impossible. Yet, the fact that rapidly developing larvae under multiple chronic stressors can still grow at rates significantly higher than adults, as demonstrated in the present study, is truly remarkable. This suggests that compensatory energy partitioning could underlie some of the reasons for heightened developmental and physiological variation, surprising phenotypic plasticity, and energetic resilience, sometimes observed in early stages (see Pan and Hunt von Herbing, 2017; Pan et al., 2015, 2018; Griffiths et al., 2021).

Energetic imbalance can also be represented by the non-linear patterns of metabolic scaling exponents (b) during ontogenetic shifts (larva to juvenile), when b shifts from 1.0 to <1.0, and in some cases from >1.0 to <1.0, when factors other than growth are responsible (Post and Lee, 1996; Bochdansky and Leggett, 2001; Glazier, 2005). In the present study, control yolk-sac larvae had a mean b of 0.83±0.34 (CI=0.49–1.17) and did not differ from the value of b=0.83 (for larvae up to 40 dpf) reported for this species (Pelster and Burggren, 1996; Bagatto et al., 2001; Barrionuevo and Burggren, 1999; Rombough and Drader, 2009; Lucas et al., 2014). But it was not significantly different from b=0.75 (Gillooly et al., 2001; Gillooly and Allen, 2007) or from isometry (b=1.0), which is often reported for developing fishes and amphibians (Giguére et al., 1988; Hunt von Herbing, 2006; Moran and Wells, 2007; Mueller et al., 2011) and for pelagic larval invertebrates (Glazier, 2006). In contrast, all log–log metabolic-scaling relationships of zebrafish yolk-sac larvae raised at chronic high temperature, hypoxia or high temperature/hypoxia were not significant. Although scaling exponents (b) might not be helpful for understanding the effects of chronic stressors on metabolic scaling in developing fish, our results agree with predictions of the MLBH (Glazier, 2005, 2008, 2014, 2020), which holds that temperature may not affect the scaling exponent (b). But Glazier (2018) asserts that for early staged ectotherms, ‘relationships of temperature on ontogenetic metabolic scaling appears to be quite complex’, and therefore some of the variation in temperature and hypoxia (negative, absent and sometimes even positive) may be a result of small sample sizes and body size ranges (see Glazier, 2005, 2020; Pan and Hunt von Herbing, 2017).

Results from the present study agree with predictions of MLBH that while scaling exponents (b) did not differ, scaling coefficients (log a) did, and larvae acutely exposed to hypoxia (log a=3.35±1.52) and high temperature/hypoxia (log a=2.61±1.55) had higher values of a (or log L_a) than controls (log a=2.11±0.97). Although we cannot explain the reason for the differences among the values of a (or log L_a), one possible source of variation may be that stressors (i.e. hypoxia and high temperature) have different effects on larval activity and on metabolic costs of activity (Kauffman, 1990; Hunt von Herbing, 2006). Activity levels have been shown to affect scaling relationships (Glazier, 2008, 2010), but were not measured in the present study as they are difficult to quantify in larval fish and we assumed routine activities. Thus for D. rerio, the effects of temperature and hypoxia are mixed and the strong positive isometric metabolic scaling between b and log L_a across two temperatures (28 and 30°C) and two oxygen concentrations (11 and 21 kPa) may be as Tan et al. (2019) suggested, owing to high sustained growth and metabolic rates characteristic of larval fish. However, in careful consideration of making general conclusions for scaling, Glazier (2009) states this result may be because of a statistical autocorrelation between b and a (or L_a), and examination of the relationship between b and L_m (i.e. the mid-point or mass-specific metabolic rate of the regression line) is needed.

In our study, autocorrelation is possible as scaling slopes were positive for b and L_a, but negative for b and L_m. Yet, aside from observations of Gould (1966) and Heusner (1991), it is only relatively recently that the elevation (a) or metabolic level (L_a or L_m) of the fitted data in metabolic–mass scaling relationships is considered important owing to systemic variation with slope (Glazier, 2005, 2008, 2009; Tan et al., 2019). We suggest that for early life history stages with complex and plastic developmental timing and patterning, it is insufficient to consider only the variation in the exponent b with size while neglecting the allometry coefficient (a) (see West and West, 2013). It is clear from current and past work on fish early life stages, as well as for the 179 animal and plant species evaluated in Glazier (2020), that metabolic scaling is plastic and not a result of intrinsic (physical) constraints related to body design (Glazier, 2020; Moyano et al., 2018). Acute or chronic exposure to stressors will be expressed in a stress response whose duration and magnitude is determined by the duration and magnitude of the consequences on the physiology and incurred allostatic load, which in turn will affect energy budgets and influence scaling relationships. Nowhere is this more evident than in rapidly developing and growing, energy-limited early life history stages.

Using a multifactorial experimental design matrix of temperature (28 and 31°C) and oxygen concentration (P_O2 of 11 and 21 kPa), we determined that the dynamics of the stress response differed between the demand side (i.e. temperature) and supply side (i.e. oxygen concentration) of the energy cascade, with the effects of high temperature dominating that of hypoxia in D. rerio yolk-sac larvae. Differential responses to environmental stressors may result in acceleration of epigenetic strategies enhancing developmental and metabolic plasticity, leading to unexpected resilience and adaptive capacity in fish early life stages. This plasticity may become fixed by genetic assimilation and act as a ‘pacemakers’ (see Glazier, 2015) or as an epigenetic factor for the evolution of development and function (see Waddington, 1942, 1953; West-Eberhart, 2003), and may be critical for species longevity under long-term environmental change. The occurrence of multiple non-significant scaling relationships when larvae were exposed to chronic high temperature and hypoxia suggest that fish larvae, with their inherently high development and growth rates, small aerobic scopes and limiting energy supplies may be unable to cope with long-term stress because of cumulative energetic costs (allostatic load), but can, to some extent, adjust to acute or short-term stress. A need for a different, perhaps novel non-linear or fractional probability allometric approach (see West and West, 2013) may be more suited to quantify how complex adaptive systems (exemplified by rapidly developing fish life stages) may respond to what seem to be stage-specific selective pressures, while at the same time coping with rapid self-organizational, developmental and functional change.

We are grateful for help from Drs Warren Burggren and Jose Mendez-Sanchez for technical assistance with the respirometry and metabolic measurements. We also thank Dr Doug Glazier for discussion and comments on preparation of the manuscript.