Fish embryo vulnerability to combined acidification and warming coincides with a low capacity for homeostatic regulation

Embryonic development is a critical period in the lifecycle of most organisms (Hamdoun and Epel, 2007). This could be particularly true for ectothermic species that release their eggs into the ocean (Przeslawski et al., 2015), which is expected to warm and acidify at an unprecedented rate over the coming decades (IPCC, 2019). The vulnerability of embryos and early larvae to climatic changes (and environmental factors in general) is often attributed to a lack of adult-like organ systems involved in energy and ion homeostasis (Melzner et al., 2009b; Hamdoun and Epel, 2007). However, knowledge about the development of homeostatic capacity and its contribution to environmental tolerance is still incomplete (Burggren and Bautista, 2019), making it difficult to identify lifecycle bottlenecks and climatic risks (Burggren and Mueller, 2015; Esbaugh, 2018).

The eggs of aquatic ectotherms such as fish are permeable to dissolved gases (e.g. O₂ and CO₂) and thermally equilibrated with their environment (Finn and Kapoor, 2008). Embryos are therefore directly exposed to changes in water temperature and CO₂-driven acidification without having regulatory (defensive) organ systems like gills (Finn and Kapoor, 2008). Instead, early (cleavage) stages are thought to be protected by maternally provisioned (passive) defenses such as non-bicarbonate pH buffering and constitutive heat-shock proteins (Melzner et al., 2009b; Hamdoun and Epel, 2007). These mechanisms may support embryonic resilience to natural environmental variability (Hamdoun and Epel, 2007), but the level of innate robustness is probably species specific and, in some cases, insufficient to cope with the challenges posed by anthropogenic climate change (Przeslawski et al., 2015; Dahlke et al., 2017). After the cleavage stage, developmental control and defense are transferred from maternal factors to those synthesized from the embryonic genome (Tadros and Lipshitz, 2009), and it is expected that the progressive differentiation of specialized cells (ionocytes) and tissues promotes active homeostatic regulation and thus improved environmental tolerance (Alderdice, 1988; Rombough, 1997; Melzner et al., 2009b).

As inferred from studies on adult fish, maintenance of homeostasis in a thermally dynamic environment includes adjustments to the structure and functioning of cell membranes and enzymes involved in energy (ATP) production (Somero et al., 2017). Within temperature limits that are typically narrow in embryos and larvae (Rombough, 1997), such responses (i.e. thermal acclimatization) may support normal development and optimal use of energetic resources (Schnurr et al., 2014; Scott et al., 2012). Developmental defects can result from a mismatch between ATP demand and supply capacity at critically high or low temperatures (Sokolova et al., 2012; Dahlke et al., 2017; Leo et al., 2018), as well as from thermal damage to proteins at extreme temperatures (Somero, 2010). CO₂-driven acidification hampers the diffusive release of metabolic CO₂ across epithelial surfaces, which causes an increase in internal P_CO2 and, consequently, a potentially harmful decline in pH of extracellular/intracellular body fluids (Brauner et al., 2019). Restoring acid–base balance requires ATP-intensive ion transport mechanisms (Pörtner, 2008), including proton (H⁺) excretion and bicarbonate (HCO₃⁻) accumulation (Brauner et al., 2019). The additional ATP demand for CO₂ compensation is expected to reduce embryonic/larval growth efficiency and thermal tolerance by tightening energy supply constraints at critical temperatures (Pörtner, 2008; Dahlke et al., 2017; Dahlke et al., 2018). Increased climate vulnerability of embryos and larvae compared with adults probably represents a common feature among vertebrate and invertebrate taxa living in different climate zones and habitats (Rombough, 1997; Przeslawski et al., 2015). Knowledge regarding the ontogeny of regulatory functions in relation to environmental tolerance, acclimation potential and energy budgeting may help advancing the concept of early-life vulnerability and its potential for directing future eco-physiological research.

Here, we used a marine cold-water fish (Atlantic cod, Gadus morhua) to investigate (i) whether vulnerability to end-of-century acidification (increase in P_CO2 from 400 to 1100 µatm, −0.4 pH) and warming (+3.5°C; Fig. 1A,B) during early embryogenesis coincides with low homeostatic capacity, and (ii) whether this environmental challenge modifies the development of energy-intensive homeostatic functions, which can inform about acclimatization potential (Burggren, 2018). Embryo vulnerability was quantified based on daily mortality rates until hatching. Homeostatic capacity was assessed at five stages until yolk sac absorption (Fig. 1C) through measurements of ion transport and ATP synthesis capacity [in vitro enzyme activity and/or protein expression of mitochondrial F₁F_O-ATP-synthase (hereafter ATP-synthase), Na⁺/K⁺-ATPase (NKA), V-Type H⁺-ATPase (VHA), Na⁺/HCO₃⁻ cotransporter 1 (NBC1) and Na⁺/K⁺/2Cl⁻ cotransporter 1 (NKCC1)] in combination with immunohistological analyses of ionocyte morphology. Cod represents a suitable model because rearing methodology (Puvanendran et al., 2015; Dahlke et al., 2017) and protocols for biochemical analyses and immunolocalization of relevant ion transporters are well established (e.g. Kreiss et al., 2015; Michael et al., 2016b). Furthermore, our previous experiments confirmed the vulnerability of cod embryos to acidification and warming (Dahlke et al., 2017), and data available on ion-regulatory mechanisms in gill tissue of adult cod (e.g. Kreiss et al., 2015; Michael et al., 2016b) allow for quantitative and qualitative comparisons between embryos, larvae and adults.

This experiment was conducted in Norway in 2014 according to local regulations of the Norwegian Animal Research Authority (Forsøksdyrutvalget, permit: FOTS ID 6382).

Mature Atlantic cod, Gadus morhua Linnaeus 1758, were caught by longlining in the southern Barents Sea (Tromsøflaket, approximately 70.5°N, 18°E) in March 2014. After being transported to the Centre for Marine Aquaculture of Nofima outside Tromsø (Norway), the fish were held in a flow-through tank (25 m³, supplied with fjord water) under ambient photoperiod, salinity (34 ppt), pH (8.1) and temperature (5–6°C) conditions.

Gametes used for in vitro fertilization were obtained by means of strip spawning from randomly selected females (n=6, 66–94 cm length) and males (n=12, 59–91 cm length; Table S1). All fertilizations were conducted within 30 min of stripping according to a standard protocol (Brown et al., 2003). In brief, each of six egg batches was divided into two equal portions to be fertilized under two different P_CO2 conditions (400 and 1100 µatm, pH_FreeScale 8.1 and 7.7) using filtered (0.2 µm) seawater adjusted to 6°C. To maximize genetic diversity, the eggs of each female were fertilized with a sperm mix from two different males. Fertilization success (Table S1) was determined in subsamples of 100 eggs (n=3 subsamples per batch and P_CO2 treatment), which had been incubated in sealed Petri dishes (20 ml water volume) for 12 h at 6°C before being photographed with a digital camera mounted to a binocular. Eggs (8 or 16 cell stage) with a clear and regular cleavage pattern were considered fertilized.

A full-factorial design with two levels of P_CO2/pH (400 µatm, pH 8.1 and 1100 µatm, pH 7.7) and two temperatures (6 and 9.5°C) was used for incubation of eggs and larvae (Fig. S1). Eggs previously fertilized at different P_CO2 conditions were subdivided into different temperature groups with the same P_CO2, resulting in four experimental treatments: (i) 6°C, 400 µatm; (ii) 6°C, 1100 µatm; (iii) 9.5°C, 400 µatm; and (iv) 9.5°C, 1100 µatm. Eggs for treatments iii and iv were warmed from 6°C to 9.5°C at a rate of ∼0.5°C h⁻¹. Increased temperature and P_CO2 conditions reflect end-of-century climate projections according to the IPCC high-emission scenario RCP8.5 (Van Vuuren et al., 2011). Treatment combinations were established in a flow-through incubation system, which consisted of 24 upwelling incubators (6 egg batches×2 P_CO2×2 temperatures=24) with a volume of 25 l each. The flow rate was set to 1.5 l min⁻¹ to ensure even distribution and mixing of eggs within the incubator. Two flow-through header tanks (150 l volume, one for each temperature) were connected to the water supply pipes and equipped with a multi-channel feedback system (IKS-Aquastar, IKS, Karlsbad, Germany) to adjust (and control) elevated pH/P_CO2 values online via infusion of pure CO₂. Automatic recordings (every 30 min) of temperature and pH values within the incubation systems were referenced against daily measurements of seawater pH/temperature with a laboratory-grade pH-electrode (InLab Routine Pt 1000, Mettler Toledo, Columbus, OH, USA) connected to a WTW 3310 pH-meter (Weilheim, Germany). Prior to each measurement, the electrode was recalibrated against tempered Tris–HCl seawater buffers (Dickson et al., 2007) to convert pH_NBS readings to the free proton concentration scale of seawater pH (pH_F) (Waters and Millero, 2013). The CO2SYS program (Lewis et al., 1998) was used to calculate P_CO2 values based on total alkalinity and dissolved inorganic carbon determined in water samples (n=3) from each treatment combination (taken during the running experiment). Seawater parameters are summarized in Table 1.

Egg batches (n=6, 200–520 ml; Table S1) were equally distributed across treatments (50–130 ml per incubator). Egg mortality was determined volumetrically every 24 h until the onset of hatching by draining dead (sunken) eggs into a graduated cylinder (±0.5 ml). Fertilization success data (Table S1) were used to estimate the volume of fertilized and unfertilized eggs within each incubator. Daily mortality rates were calculated as percentages relative to the estimated volume of eggs that was present in the incubator on the previous day. To better resolve changes in embryo vulnerability in relation to developmental age, egg mortality was also displayed as a function of degree-days (days post-fertilization multiplied by incubation temperature) (Trudgill et al., 2005). Total embryo mortality was calculated as the percentage of fertilized eggs that died until hatching stage.

Subsamples of eggs and larvae used for analyses of enzyme activity and protein expression were obtained from each treatment at developmental stages I–V (Fig. 1C). These stages were chosen because they represent developmental landmarks that can be reliably identified under a stereomicroscope: stage I: end of cleavage period, stage II: end of gastrulation period, stage III: 50% eye pigmentation, stage IV: peak-hatching, stage V: complete yolk sac absorption (Hall et al., 2004). Developmental progress was monitored in each incubator every 12 h. Prior to sampling, the aeration and water supply of the incubators was turned off so that eggs and larvae accumulated at the water surface. Eggs and larvae were then concentrated within a small kitchen sieve, pipetted into 1.5 ml cryovials (∼500 individuals per vial) and immediately frozen in liquid nitrogen after excess water was removed by pipetting. Larvae were previously centrifuged to the bottom of the vial (∼3 s at 500 rpm). Additional samples of larvae (stage IV) reared under control conditions (6°C, 400 µatm) were fixed in 4% buffered formaldehyde (pH 7.4) and stored in 70% phosphate-buffered saline (PBS)-buffered isopropanol (pH 7.4) for immunolocalization of NKA, VHA, ATP-synthase, NBC1 and NKCC1.

Crude extract preparation followed the same protocol as described elsewhere (Kreiss et al., 2015; Melzner et al., 2009a; Michael et al., 2016a,b), with the exception that different sample-to-extraction buffer dilutions were applied for egg (stage I–III) and larval (stage IV and V) stages to ensure similar concentrations of biologically active tissue in all extracts. This adjustment was necessary because a given fresh mass of egg-shelled embryos contains approximately half the number of individuals (tissue) than the same mass of frozen larvae (many individuals compressed into one). Given that soluble protein content of egg and early larval stages was reported to vary by less than 3% (Finn et al., 1995), we normalized all measurements against protein content of the crude extract to account for potential differences in the number of individuals per sample. In brief, crude extracts were produced by homogenizing ∼100 mg of frozen sample suspended in 5 (eggs) or 10 (larvae) volumes of ice-cold extraction buffer, which contained 50 mmol l⁻¹ imidazole (pH 7.8), 250 mmol l⁻¹ sucrose, 5 mmol l⁻¹ Na₂-EDTA, 0.1% sodium deoxycholate, 5 mmol l⁻¹ β-mercaptoethanol and 0.2 ml proteinase inhibitor cocktail (Sigma-Aldrich, Taufkirchen, Germany; catalog no. P8340). Samples were homogenized in a temperature-controlled tissue homogenizer (Precellys 24, Bertin Technologies, Aix-en-Provence, France) set to 0°C, 6000 rpm (3×15 s). Cell fragments were removed by centrifugation (1000 g for 10 min at 2°C) and the supernatant (crude extract) was taken for analysis. Final protein concentration of the crude extracts were determined according to Bradford (1976), using bovine serum albumin (BSA) as standard. Differences in protein concentration between treatments were not significant (Fig. S2A, Table S3). Half of each sample was used for the enzyme assay and the other half for immunoblotting procedures.

Maximum total ATPase activity and fractional (inhibitor-sensitive) activity of NKA, VHA and ATP-synthase (reversed catalysis) were measured in crude extracts by means of a coupled enzyme assay based on pyruvate kinase (PK) and lactate dehydrogenase (LDH) as described elsewhere (Kreiss et al., 2015; Melzner et al., 2009a; Michael et al., 2016a,b). Transporter specific enzyme activity was determined as the difference between total ATPase (TA) activity and inhibitor-insensitive (residual) ATPase activity, using inhibitor concentrations previously applied to gill tissue samples of adult cod: 5 mmol l⁻¹ ouabain for NKA (Michael et al., 2016a), 60 μmol l⁻¹ oligomycin for ATP-synthase (Michael et al., 2016b) and 0.1 μmol l⁻¹ bafilomycin A1 for VHA (Kreiss et al., 2015). The assay was conducted in a micro-plate reader format under temperature-controlled conditions. Samples from 6.0 or 9.5°C incubations were assayed at both temperatures to assess acclimation effects on enzyme activity. The reaction process (oxidation of NADH, hydrolysis of ATP) was initiated by the addition of crude extract to 20 volumes of reaction buffer containing 100 mmol l⁻¹ imidazole, pH 7.8, 80 mmol l⁻¹ NaCl, 20 mmol l⁻¹ KCl, 5 mmol l⁻¹ MgCl₂, 5 mmol l⁻¹ ATP, 0.24 mmol l⁻¹ Na-NADH₂, 2 mmol l⁻¹ phosphoenolpyruvate and 12 U ml⁻¹ PK with 17 U ml⁻¹ LDH in a PK/LDH enzyme mix (Sigma-Aldrich). All samples were arranged in a randomized order and measured in quadruplicate with 10 readings over a time period of 10 min at λ=339 nm. Inhibitor-insensitive activity was calculated based on the extinction coefficient of NADH (ε=6.31 l mmol⁻¹ cm⁻¹) and expressed as micromoles of ATP consumed per mg protein in the crude extract per hour (µmol ATP mg⁻¹ protein h⁻¹).

The specificity of all primary antibodies used for protein quantification and/or localization was confirmed by western blotting (Fig. S3). Antibodies against NKA, NKCC1 and NBC1 were previously used to study ion regulation processes in gill tissue of cod (Kreiss et al., 2015; Melzner et al., 2009a; Michael et al., 2016b). The polyclonal antibody against NBC1 was developed specifically for cod (Michael et al., 2016b). Details on primary antibodies are given in Table S2. A mix of subsamples from randomly selected crude extracts (developmental stage IV) was split into membrane and cytosolic fractions by ultra-centrifugation (350,000 g for 30 min at 4°C). A 15 μl sample from each fraction was separated by SDS-PAGE on 8–10% polyacrylamide gels according to Laemmli (1970), and transferred onto PVDF membranes (Immuno-BlotTM, Bio-Rad, Munich, Germany) using a tank blotting system (Bio-Rad). To prevent non-specific protein binding, blots were quenched for 1 h at room temperature (RT) in TBS-Tween buffer [TBS-T, 50 mmol l⁻¹ Tris-HCl, pH 7.4, 0.9% (w/v) NaCl, 0.1% (v/v) Tween20] containing 5% (w/v) non-fat skimmed milk powder. Incubation of blots with primary antibodies was done overnight at 2°C (dilutions in TBS-T are given in Table S2). After washing with TBS-T, blots were incubated for 1 h (RT) with horseradish-conjugated goat anti-rabbit/anti-mouse IgG antibody (Pierce, Rockford, IL, USA, diluted 1:50,000 in TBS-T). Protein bands were visualized with ECL Advanced Western Blotting Detection Reagent (GE Healthcare, Munich, Germany), and imaged with a cooled charge-coupled device camera (LAS-1000: Fuji, Tokyo, Japan). Protein bands (Fig. S3) were referenced to pre-labelled SDS-PAGE standards (Bio-Rad).

A 48-fold slot-blot filtration system (Hoefer PR 648, Amersham Biosciences, Freiburg, Germany) was used to quantify the expression of NKA, ATP-synthase, NBC1 and NKCC1. VHA was excluded because of insufficient antibody reactivity. After activation in 100% methanol, PVDF membranes were equilibrated for 30 min in transfer buffer [10 mmol l⁻¹ NaHCO₃, 3 mmol l⁻¹ NA2CO₃, 20% (v/v) methanol, 0.025% (w/v) SDS, pH 9.5–9.9]. Crude extracts were diluted 1:10 in electrophoresis running buffer [25 mmol l⁻¹ Tris, 192 mmol l⁻¹ glycine, 0.1% (w/v) SDS], and 80 µl of each sample were applied to the system followed by repeated rinsing with transfer buffer (3×500 µl). A dilution series of pooled samples was used as a reference standard on each membrane (Fig. S3F), which were always loaded with randomly ordered samples from two egg batches (2×20 samples) yielding 12 runs in total (3×2 batches×4 proteins). After the loading process, membranes were immediately blocked for 1 h (RT) with TBS-T buffer containing 5% (w/v) non-fat skimmed milk powder. Methods for protein detection and imaging were the same as described above (see Antibody specificity). Signal intensity was analyzed using AIDA Image Analyzer software (version 3.52, Raytest, Straubenhardt, Germany). Values were recalculated from the reference curve and expressed as arbitrary units per mg protein content of the original sample (a.u. mg⁻¹ protein).

Fixed larvae (stage IV) were rehydrated in 0.1 mol l⁻¹ PBS (pH 7.4) and incubated for 30 min in 3% BSA to block non-specific binding during immunolocalization of NKA, ATP-synthase, VHA, NBC-1 and NKCC. Incubation with primary antibodies diluted 1:300 in PBS was done overnight at 2°C (see Table S2 for details on primary antibodies). After being carefully rinsed with PBS, larvae were incubated for 2 h (RT) with two secondary antibodies (DyLight^© 488 anti-mouse and DyLight^© 594 anti-rabbit, Jackson ImmunoResearch, West Grove, PA, USA, diluted 1:300 in PBS) for co-localization of NBC1/NAK, NBC1/NKCC1, NBC1/ATP-synthase and NKA/VHA, respectively. Negative controls were performed for each pair without application of the primary antibody (Fig. S3G). Finally, larvae were rinsed once more with PBS and placed on a fluorescence slide prior to image acquisition with an inverse confocal laser-scanning microscope (Leica TCS SP5 II, Leica, Wetzlar, Germany). The total number of mitochondria (ATP-synthase)-rich ionocytes on the body surface of five larvae was evaluated using a cell-counter plugin for ImageJ^© (ATCN 1.6). Ionocyte density was estimated from the number of cells divided by body area (cells mm⁻²).

Statistical analyses were conducted with the open source software R (www.r-project.org). If normality and homoscedasticity were not violated (assessed via Q–Q plots), linear mixed models (LME, package ‘lme4’; Bates et al., 2014) were applied (total mortality, enzyme activity-to-expression ratios). Otherwise, generalized linear mixed-effect models (GLMM, package ‘lme4’) were used (daily mortality, enzyme activity, protein expression) and appropriate probability distributions were assessed using the ‘MASS’ package (https://CRAN.R-project.org/package=MASS; Venables & Ripley, 2002). In all cases, different levels of temperature, P_CO2 and developmental stage were treated as fixed factors while egg batch was included as a random factor. The package ‘lsmeans’ (Lenth, 2016) with Tukey's P-value corrections was used to conduct post hoc pairwise comparisons of single model factors. All data are presented as means (±s.e.m.) and statistical tests with P<0.05 were considered significant. Statistical results are summarized in Table S3.

The difference in mean (±s.e.m.) fertilization success under control conditions (88.4±4.3% at 6.0°C, 400 µatm CO₂) and elevated P_CO2 (87.0±4.5% at 6.0°C, 1100 µatm CO₂) was not significant (paired t-test, P=0.08, n=6). Temperature-dependent development times until stages I–V (Fig. 2A,B) were highly synchronous among egg batches (Fig. S2); differences between P_CO2 treatments were not detectable by visual inspection (data not shown).

The influence of elevated P_CO2 on daily mortality rates of cod embryos was a function of incubation temperature and developmental stage (Fig. 2A,B), as indicated by an interaction between age, P_CO2 and temperature (GLMM: P=0.0233). Increased mortality due to elevated P_CO2 occurred primarily during early development (cleavage and gastrulation, stages I−II) and in combination with warming. Few losses were observed during organogenesis, segmentation and hatching (stages II−IV). Total mortality until hatching (Fig. 2C) of embryos exposed to elevated P_CO2 and warming (28.2±0.5%) was 2-fold higher than in the control group (13.1±1.4%, post hoc test: P<0.0001). Differences in mortality under either elevated P_CO2 (14.1±1.3%) or warming (16.1±1.2%) were not significant. These results clearly demonstrate that early embryogenesis (and particularly gastrulation) is a critical bottleneck with respect to the combined effects of elevated P_CO2 and temperature on Atlantic cod (Fig. 2D).

TA activity and specific activity of VHA, NKA and ATP-synthase increased with developmental stage in a similar, sigmoidal way in all treatment combinations (GLMM, P<0.0001; Fig. 3A,B). The activity of NKA, VHA and ATP-synthase was extremely low at stage I and II (end of gastrulation). Thereafter, activity increased exponentially until hatching (stage II–IV) while a less rapid increase was observed between hatching and yolk sac absorption (stage IV−V; Fig. 3B). Between the end of gastrulation and yolk sac absorption (stage II–V), enzyme specific activity increased 40- to 60-fold.

Incubation at elevated P_CO2 had no effect on TA activity (GLMM, P=0.145), but caused a significant increase (20–30% at stage III–V) in the activity of NKA (P=0.021), VHA (P<0.001) and ATP-synthase (P<0.001; Fig. 3B). Accordingly, there was a significant decrease in residual ATPase activity (the difference between TA activity and the sum of NKA, VHA and ATP-synthase activity) at elevated compared with control P_CO2 (GLMM, P<0.001; Fig. 3D−G). The effect of elevated P_CO2 on enzyme activity (TA, NKA, VHA and ATP-synthase) did not differ between temperature treatments (Fig. 3A,B). In both CO₂ treatments, the summed contribution of NKA, VHA and ATP-synthase activity to TA activity decreased with warming, resulting in larger fractions of residual ATPase activity at 9.5°C compared with 6.5°C (GLMM, P=0.023; Fig. 3D−G).

The protein expression of NKA, ATP-synthase and secondary ion transporters NKCC-1 and NBC1 (Fig. 4A−D) increased from low levels during cleavage and gastrulation (stage I and II) to a maximum at hatching (stage IV), followed by a 30–40% decrease during yolk sac absorption (GLMM, stage effect: P<0.001). A trend towards increased protein expression at elevated P_CO2 was statistically not significant (Table S3). Incubation at increased temperature (9.5°C) led to an overall reduction in protein expression of ∼40% (stage III–V) compared with that at 6.0°C (GLMM, P<0.05).

Lower protein expression of NKA and ATP-synthase at 9.5°C compared with 6.0°C did not result in a difference in enzyme activity when assayed at a common temperature (Fig. 5A,B; Table S3). Moreover, enzyme activity of NKA and ATP-synthase increased throughout yolk sac absorption (stage V) despite a reduction in expression between stage IV and stage V (Fig. 5A,B). As a result, enzyme activity-to-expression ratios (indicating enzyme catalytic power) increased significantly due to warming and developmental progress (GLMM, P<0.0001; Fig. 5C,D).

All targeted ion transporters (NKA, VHA, NBC1 and NKCC1) and mitochondrial ATP-synthase were localized in epidermal ionocytes of newly hatched larvae (stage IV; Fig. 6A−F). Ionocytes located on the yolk sac and trunk were larger (∼50 µm diameter, Type I) than densely clustered ionocytes (∼30 µm, Type II) within the primordial gill cavity (Fig. 6A, top). Ionocyte density on the yolk sac ranged between 200 and 300 cells mm⁻², while more than 1000 cells mm⁻² were counted within the gill cavity (Fig. 6G). NBC1 was also localized around the apical pit of mucous cells (Fig. 6A,B), which occur across the entire body surface (see Ottesen and Olafsen, 1997). Immunoreactivity of adhesive mucous droplets (Fig. 6A, bottom) was probably caused by non-specific binding of antibodies.

Embryonic life stages are considered particularly vulnerable to environmental factors because of a lack of adult-like organ systems and associated capacity for homeostatic regulation (Hamdoun and Epel, 2007; Melzner et al., 2009b). In line with this hypothesis, we found that Atlantic cod embryos are most vulnerable to the combined effects of ocean acidification and warming during early development, especially during gastrulation. The extremely low capacity for homeostatic regulation (activity and expression of ion transporters) throughout cleavage and gastrulation suggests that early embryogenesis is mainly protected by inherited (maternal) defense mechanisms. By the time of hatching, however, cod larvae were less sensitive and had ion transport and ATP production capacities similar to those of specialized adult gill tissue. Moreover, enhanced enzyme activity and modulation of protein expression levels in response to acidification and warming imply that post-gastrulation stages are increasingly capable of responding to changing environments via plasticity in the development of homeostatic mechanisms. Based on these results, a conceptual model of fish early-life vulnerability and homeostatic plasticity is proposed in Fig. 7.

The combined effects of future ocean acidification and warming on offspring survival demonstrated in this study are consistent with previous work on Atlantic cod (Dahlke et al., 2017, 2018), Antarctic dragon fish Gymnodraco acuticeps (Flynn et al., 2015), sand lance Ammodytes dubius (Murray et al., 2019) and many other marine organisms (Przeslawski et al., 2015). Collectively, these data point to considerable risks for marine species and ecosystems if CO₂ emissions continue to rise over the coming decades (Hoegh-Guldberg et al., 2018). It should be noted, however, that our approach did not account for the possibility that species could adapt to the expected long-term changes in acidity and temperature through selection and/or transgenerational plasticity, i.e. parental effects (Chevin et al., 2010; Burggren, 2018). Several studies have shown that parent–environment interactions, including epigenetic mechanisms, can alter offspring responses (reaction norms) to environmental factors (e.g. Dahlke et al., 2016, Shama et al., 2014, Schunter et al., 2018), potentially enabling physiological adaptation to moderate climate change (but see Byrne et al., 2019). However, phylogenetic analyses of past adaptation rates (Quintero and Wiens, 2013; Comte and Olden, 2017) and observed distribution shifts of marine species to higher latitudes and/or greater depths in response to global warming (Burrows et al., 2019; Pinsky et al., 2013) indicate that, for many species, adaptation of physiological thresholds may be too slow to keep pace with ongoing anthropogenic climate change (Hoegh-Guldberg et al., 2018). Our results primarily emphasize the importance of understanding the development of homeostatic capacity and associated lifecycle bottlenecks as a basis for reliable climate risk assessment.

Critical periods during development often coincide with functional transitions (Burggren and Mueller, 2015), such as those occurring during gastrulation. This early period is characterized by complex morphogenetic processes, including the maternal-to-zygotic transition (MZT), where developmental control is handed over from maternally provisioned factors (e.g. mRNAs, enzymes and chaperones) to those synthesized from the embryonic genome (Schier, 2007; Kimmel et al., 1995). Moreover, gastrulation leads to the formation of germ layers which give rise to different tissues and organ systems (Kimmel et al., 1995). Any defect at this stage may cause disproportionately serious damage, leading to either instantaneous mortality or deformities. In fact, high vulnerability of gastrula stages to diverse abiotic stressors (e.g. temperature, CO₂, UV radiation, hypoxia or toxicants) was demonstrated in different gadoid species (Dahlke et al., 2017; Kouwenberg et al., 1999; Nahrgang et al., 2017; Skjærven et al., 2013; Wieland et al., 1994), as well as in model organisms like zebrafish, Danio rerio (Ali et al., 2011; Jesuthasan and Strähle, 1997; Sawant et al., 2014), fruit fly, Drosophila melanogaster (Bergh and Arking, 1984), and clawed frog, Xenopus laevis (Heikkila et al., 1985; Metikala et al., 2018; Degitz et al., 2000). Overall, these observations identify the gastrulation period as a critical lifecycle bottleneck for many, if not all, metazoans.

Proposed mechanisms underlying embryo mortality from external stressors include impairment of functions related to cell division and cell motility, such as pH-sensitive and energy-dependent microtubule dynamics (Cheng et al., 2004; Zalik et al., 1999). It is thus likely that the increasing capacity for pH and energy homeostasis caused the observed decrease in cod embryo vulnerability after gastrulation (Fig. 7), although additional mechanisms may have been involved. For instance, studies on vertebrate and invertebrate model organisms suggest that inducible heat shock responses and apoptotic pathways occur after MZT (Rupik et al., 2011; Hamdoun and Epel, 2007), possibly contributing to reduced stress sensitivity of post-gastrulation stages.

The capacity of fish embryos and larvae for homeostatic regulation is thought to be linked to the differentiation of extrabranchial (epidermal) ionocytes (Burggren and Bautista, 2019; Varsamos et al., 2005). Proportional changes in the expression and activity of primary (NKA, VHA) and secondary ion transporters (NBC1, NKCC1), as well as their localization in mitochondria (ATP-synthase)-rich ionocytes suggest that early larvae and adults utilize similar ion regulation mechanisms to defend pH homeostasis against respiratory and environmental acidosis (Fig. 6H). Notably, the high enzyme activity of cod larvae (whole-animal extracts, stage V) was similar to that in specialized adult gill tissue, indicating that these fragile life stages are already equipped with powerful homeostasis systems. From an evolutionary perspective, this appears plausible given that fish larvae are among the fastest growing vertebrates (30% day⁻¹ in cod; Finn et al., 2002) with correspondingly high requirements for efficient removal of metabolic CO₂ and acid–base regulation (Brauner and Rombough, 2012). In fact, an early study (Ishimatsu et al., 2004) showed that (short-term) survival of marine fish larvae is possible at P_CO2 levels (>10,000 µatm) 10-fold higher than those projected for the end of this century.

While all investigated ion transporters are highly expressed by ionocytes of newly hatched cod larvae (Fig. 6), it was not possible to determine the contribution of these cells to whole-organism homeostatic capacity. However, observed changes in whole-organism enzyme activity and protein expression due to elevated P_CO2 and temperature are consistent with the idea that developmental (or homeostatic) plasticity allows species to maintain fitness within a limited range of environmental conditions (Burggren, 2018). Early mortality of embryos with insufficient homeostatic capacity (i.e. selection of beneficial genotypes) provides an alternative explanation for the observed treatment effects on NKA, VHA and ATP-synthase activity and expression. However, selection through treatment-related mortality would have narrowed the data frequency distribution (e.g. at elevated P_CO2 compared with control conditions), which was clearly not the case (Fig. S2,D).

Higher catalytic power of enzymes (Fig. 5) in combination with higher kinetic energy at warmer temperatures may support maintenance of ionic balance despite lower expression (and energy expenditure) of ion transporters (Fields et al., 2015; Somero, 1995). This view is supported by an inverse relationship between incubation temperature and ionocyte density previously demonstrated in cod larvae (Dahlke et al., 2017). At extreme temperatures, however, energetic benefits associated with warm temperature acclimation are increasingly outweighed by rising maintenance costs and constraints on mitochondrial energy production (Dahlke et al., 2017; Leo et al., 2018).

Additional homeostatic requirements of cod embryos at elevated P_CO2 (1100 µatm) may have triggered the establishment of increased enzyme activity, as was previously demonstrated for gill NKA of adult cod (Melzner et al., 2009a) and eelpout Zoaces viviparus (Deigweiher et al., 2008), although at much higher CO₂ concentrations (6000 and 10,000 µatm, respectively). Increased activity of VHA, NKA and ATP-synthase in cod embryos and larvae incubated at elevated P_CO2 implies that pH homeostasis and normal development under these conditions involves additional ATP-dependent ion transport in specialized ionocytes (Fig. 6H), but probably also in other cell types and tissues. Increased ATP demand by ionocytes under OA conditions may primarily relate to the activity of NKA, which facilitates electroneutral proton export via apical Na⁺/H⁺ exchanger (NHE) proteins, as well as bicarbonate buffering via basolateral NBC1 and NKCC1 by establishing and maintaining the required Na⁺ and Cl⁻ gradient (Brauner et al., 2019; Fig. 6H). The observed increase in VHA activity may not be directly related to proton export by ionocytes, as our results and previous studies on marine fish do not provide evidence for an apical orientation of this transporter (Allmon and Esbaugh, 2017; Brauner et al., 2019). Instead, it is possible that the contribution of VHA to the maintenance of membrane potentials and various other cellular functions (e.g. acidification of lysosomes) (Tresguerres, 2016) led to an increased VHA activity and thus ATP demand at elevated P_CO2.

In line with previous work on echinoderms (Pan et al., 2015), our results (e.g. smaller residual enzyme activity; Fig. 3D–G) indicate that additional costs for acid–base regulation and cellular maintenance at elevated P_CO2 are met through energy reallocation. In cod and many other fish species, altered energy budgets in response to elevated P_CO2 are reflected by reductions in developmental growth (Cattano et al., 2018; Esbaugh, 2018; Dahlke et al., 2018), sometimes still detectable at the juvenile stage (Murray et al., 2016). Accordingly, our results link CO₂-related enzyme adjustments with developmental trade-offs at the animal level (growth deficits), which in turn may contribute to increased susceptibility to natural sources of mortality, i.e. starvation and predation (Garrido et al., 2015).

Although it is clear that in vitro analyses of enzyme activity and expression cannot resolve the actual (in vivo) biochemical responses of an intact organism (Somero et al., 2017), we consider the presented data as reliable proxies. Firstly, the relative increase in total ATPase activity (in vitro metabolic capacity) between stage I and stage V is directly proportional to the increase in oxygen consumption (in vivo metabolic intensity) determined in cod embryos (Finn et al., 1995) over the same developmental period (Fig. S2,E). Secondly, the fractional activity of whole-larvae NKA (28.5% of TA activity at stage V) corresponds with the relative amount of available ATP (30%) that is typically allocated to regulate sodium–potassium fluxes in metabolically active tissues (Somero et al., 2017), including gill tissue (29–36%) of adult cod (Kreiss et al., 2015; Michael et al., 2016b).

A low capacity to maintain pH and energy homeostasis of early cod embryos corresponds with the concept (Fig. 7) that maternally provisioned defense mechanisms protect initial development against natural environmental variability (Hamdoun and Epel, 2007). Innate defense levels differ between locally adapted species or populations, sometimes due to parental pre-exposure (Byrne et al., 2019). Future climate changes may exceed the range of natural variability and thus innate defense levels. The establishment of improved tolerance after the critical gastrulation period most likely involves increasing capacity for homeostatic regulation associated with the differentiation of specialized cells (ionocytes) and organ systems (Varsamos et al., 2005). Regulatory mechanisms of early cod larvae are more sophisticated than previously expected, possibly reflecting a physiological prerequisite for highly active and rapidly growing fish larvae (Rombough, 2011). Capacity adjustments and modifications of regulatory mechanisms potentially support short-term (within-generation) acclimatization to environmental change, but the extent of homeostatic plasticity may cause additional energetic costs and trade-offs.

We acknowledge the technical support of Silvia Hardenberg, Elettra Leo, Felix Mark, Martina Stiasny, Catriona Clemmensen and Gwendolin Göttler. Special thanks are dedicated to the staff of the Tromsø Aquaculture Research Station and the Centre for Marine Aquaculture.