One hundred specimens of the tropical sacoglossan sea slug E. crispata (41.1 ± 3.8 mm of total length) were collected off the Florida Keys coastline and shipped to Laboratório Marítimo da Guia (LMG, Cascais, Portugal) by Tropical Marine Centre (TMC, Iberia, Portugal), a marine aquarium wholesaler recognized for its efforts on the sustainable collection and trade of reef organisms and promotion of animal welfare. One hundred and forty-four specimens of the temperate sacoglossan sea slug E. viridis (13.3 ± 0.9 mm of total length) were hand collected during low tides, in Cabo Raso (38° 42′ 34.67″ N, 9° 29′ 12.38″ W; Cascais, Portugal).

Upon arrival to the LMG aquatic facilities, organisms were randomly distributed in recirculating life support systems (RAS) according to Dionísio et al. (2013). Each RAS was composed by a 250-L holding aquaria, filled with 0.2 μm altered natural seawater (NSW), and equipped with mechanical (100 μm, TMC Iberia, Portugal), physicochemical (REEF-Skim Pro 400, TMC Iberia, Portugal) and biological (Fernando Ribeiro Lda, Portugal) filtration. All RAS were additionally equipped with UV irradiation (Vecton 600, TMC Iberia, Portugal). Ammonia (<0.5 mg/L) and nitrite (<0.05 mg/L) levels were daily checked using colorimetric test kits (Aquamerk, Merck Millipore, Germany). Overhead tank illumination was provided through dimmable LED illumination apparatus (Aquabeam 1500 Ultima NP Ocean Blue, TMC Iberia, Portugal), consisting of five white XP-G LEDs (9000 K) and five XP-E blue LEDs (50000 K). Photosynthetically active radiation (PAR) was measured (FluorPen FP100 light meter, Photo System Instruments, Czechia) and maintained at 150 ± 15 μmol photons m⁻² s⁻¹ at the water surface, and photoperiod was set to 14 h light: 10 h dark. The siphonaceous macroalgae Codium tomentosum and Bryopsis plumosa (previously acclimated for 2 days to the same conditions of stocked sea slugs) were provided ad libitum as feed source. During the first 2 weeks of laboratory acclimation, sea slugs were kept at control conditions, corresponding to the ambient temperature and pH conditions at collection sites, i.e., 26°C and pH 8.0 for E. crispata, and 18°C and pH 8.0 for E. viridis). After laboratory acclimation, E. crispata individuals were randomly divided into five 5^−L tanks per treatment (n = 5 individuals per tank, n = 25 individuals per treatment), and E. viridis into three 5^−L tanks per treatment (n = 12 individuals per tank, n = 36 individuals per treatment). Subsequently, organisms were exposed for 5 days to a gradual increase of pCO₂ and temperature levels. After this period, organisms were exposed for 8 weeks to 4 different experimental conditions, namely: (i) Control scenario – normocapnia (pH 8.0) and control temperature (26 and 18°C for E. crispata and E. viridis, respectively); (ii) hypercapnia/high CO₂ scenario – hypercapnia (pH 7.6) and control temperature; (iii) thermal stress scenario (+4°C, i.e., 30 and 22°C for E. crispata and E. viridis, respectively) and normocapnia; and (iv) thermal stress + high CO₂ combined scenario – the warming and hypercapnia scenarios.

Seawater temperature and pH were adjusted automatically by using a Profilux control system (GHL, Germany) connected to individual temperature and pH probes (GHL, Germany). The temperature was automatically upregulated by submergible heaters and downregulated using cooling systems (HC-1000A, Hailea, China). Monitoring of pH values was automatically performed (every 2 s) and adjusted via a solenoid valves system, being downregulated through the injection of a certified CO₂ gas mixture (Air Liquide, Portugal) or upregulated by aerating the tanks with atmospheric filtered air (soda lime, Sigma-Aldrich). Salinity was measured with a refractometer (V2 Refractometer, Tropical Marine Centre, Portugal) and kept at 35 ± 1 μS cm⁻¹. Seawater carbonate system speciation (Table 1) was calculated weekly based on total alkalinity (Sarazin et al., 1999), pH, temperature, and salinity measurements using the CO₂SYS software (Lewis and Wallace, 1998), with dissociation constants accordingly (Mehrbach et al., 1973).

Survival at each treatment was daily checked throughout the entire experimental period (i.e., 60 days). The integrity of the symbiosis was evaluated at the initial (T = 0), mid (T = 30), and final (T = 60) days of exposure, based on the presence of green kleptoplasts inside the digestive tubules of the sacoglossan sea slugs. Kleptoplasts were qualitatively evaluated using morphological features, namely color, symmetry and their distribution in the tubules. Images were taken using a binocular microscope (DM1000, Leica, Germany) equipped with a digital camera (DFC 450, Leica, Germany).

Variable chlorophyll a fluorescence was measured at day 60 using a PAM (Pulse Amplitude Modulated) fluorometer, comprising a computer-operated PAM-control unit (JUNIOR-PAM, Walz Heinz GmbH, Germany) and a WATER-EDF emitter-detector unit (Gademann Instruments GmbH, Germany). The actinic and saturating light was provided by a blue LED-lamp (450 nm peak and 20 nm half-band width) and supplied through a plastic fiber optic bundle (1.5 mm diameter) perpendicularly positioned to the surface of the sea slug parapodia. A saturation pulse of 2500 μmol photons m⁻² s⁻¹ with a duration of 0.8 s was applied to at least 8 slugs per treatment (previously anaesthetized as described in Cruz et al. (2012), in order to determine the fluorescence at both dark and light conditions. Sea slugs were dark-adapted for 30 min and the minimum (F_o) and maximum fluorescence (F_m) in the dark-adapted state were used to determine the variable fluorescence (F_v = F_m – F_o) and the maximum quantum yield of PSII (F_v/F_m). Sea slugs were then light-adapted at 150 ± 15 μmol photons m⁻² s⁻¹ for 30 min. The minimum (F) and maximum fluorescence (F_m’) in the light-adapted state were used to determine the variable fluorescence (ΔF = F_m’ – F) and the PSII maximum quantum yield (ΔF/F_m’) in the light-adapted state. The relETR was then calculated as:

where PAR is the photosynthetic active radiation and 0.5 compensates for irradiance being split between two photosystems.

Oxygen consumption was determined 60 days after exposure to experimental scenarios according to previously established methods (Rosa et al., 2009, 2012, 2013). Sea slugs (n = 6 per treatment) were individually incubated in sealed water-jacketed respirometry chambers (Strathkelvin, United Kingdom) containing 1 μm filtered and UV-irradiated NSW derived from the respective experimental treatments. Water volumes were adjusted in relation to animal mass (up to 3 mL) in order to minimize locomotion and stress but still allow for spontaneous and routine activity rates. Respiration chambers were immersed in Lauda water baths (Lauda-Brinkmann, Germany) to control temperature. Oxygen concentrations were recorded with Clark-type O₂ electrodes connected to a multi-channel oxygen interface (Model 928, Strathkelvin, United Kingdom). Controls (blanks) were used to correct for possible bacterial respiratory activity. Two runs of 3 h were made per individual, one exposed to light and the other in complete darkness to inhibit photosynthesis. Light or dark incubations were performed within the respective photoperiod of the animals. Oxygen concentration measurements (μmol O₂ L⁻¹) were transformed into μmol O₂ g⁻¹ L⁻¹ h⁻¹ by taking into consideration the volume of the chamber and the wet weight of the slug. Respiration was determined as the oxygen consumption rate in complete darkness, while net primary photosynthesis was determined as the oxygen production rate in the light exposed conditions, according to Baker et al. (2015).

The oxidative stress response was analyzed based on both the HSP production (HSP70/HSC70) and the activity of the antioxidant enzyme (GST). A total of 3 samples (each one containing 3 slugs) were analyzed per treatment. Samples were homogenized using an Ultra-Turrax (Staufen, Germany) in phosphate-buffered saline (PBS), pH 7.4: 0.14 M NaCl (≥99%), 2.7 mM KCl (≥99%), 8.1 mM Na₂HPO₄ (≥99%), and 1.47 mM KH₂PO₄ (≥99%), Sigma-Aldrich, United States) and centrifuged at 10,000 ×g for 15 min at 4°C. Afterward, homogenized samples were frozen at -80°C until further analyses.

Total protein measurements were determined according to Bradford (1976) adapted to 96-well microplates. Briefly, 20 μL of each sample and 200 μL of 5 % Bradford reagent solution (Sigma-Aldrich, United States) were added to a 96-well microplate and the absorbance read at 595 nm (Asys UVM 340, Biochrom, United States). Albumin bovine serum (BSA, Sigma-Aldrich, United States) dilutions (0–1 mg) were used as standards. Bradford results were then used to normalize HSP and GST results to total protein content.

The HSP70/HSC70 content was assessed by Enzyme-Linked Immunosorbent Assay (ELISA), by adapting the protocol from Njemini et al. (2005) (see more details in Supplementary Methods). Briefly, 10 μL of the homogenate supernatant was diluted in 250 μL of PBS. Afterward, 50 μL of the diluted sample was added to 96-well microplates (Nunc- Roskilde, Denmark) and allowed to incubate overnight at 4°C. After 24 h, the microplates were washed in PBS containing 0.05% Tween-20 (≥40%, Sigma-Aldrich, United States). A total of 100 μL of blocking solution (1% bovine serum albumin, Sigma-Aldrich, United States) was added to each well and left to incubate at room temperature for 2 h. After washing the microplates, 50 μL of a solution of 5 μg mL⁻¹ of primary antibody (anti-HSP70/HSC70, Acris, United States) was added to each well and then incubated at 37°C for 90 min. According to the manufacturer details, the primary antibody (anti-HSP70/HSC70) has a broad range of reactivity. The primary antibody reactivity for the species E. crispata and E. viridis was validated by Western blot. The non-linked antibody was removed by an additional washing step of the microplates. The alkaline phosphatase-conjugated anti-mouse IgG (Fab specific, Sigma-Aldrich, United States) was then used as a secondary antibody, by adding 50 μL of a solution at 1 μg mL-1 to each well and incubating the microplates for 90 min at 37°C. After three additional washing steps, 100 μL of substrate (SIGMAFAST^TM p-nitrophenyl phosphate tablets, Sigma-Aldrich, United States) was added to each well and incubated for 10–30 min at room temperature. Subsequently, 50 μL of stop solution (3 M NaOH (≥98%), Sigma-Aldrich, United States) was added to each well, and the absorbance was read at 405 nm in a 96-well microplate reader (Asys UVM 340, Biochrom, United States). The concentration of HSP70/HSC70 in the samples was calculated from a curve of absorbance based on serial dilutions (between 0 and 2 μg mL⁻¹) of purified HSP70 active protein (Acris, United States). Results were expressed in relation to the protein content of the samples, which was determined according to Bradford (1976).

The activity of the antioxidant enzyme GST was determined according to Rosa et al. (2012) and Lopes et al. (2013) and optimized for a 96-well microplate. This assay uses 1-chloro-2,4-dinitrobenzene (CDNB) as substrate, which conjugates with the thiol group of the glutathione (GSH) causing an increase in absorbance. A total of 180 μL of substrate solution (composed by 200 mM L-glutathione reduced in Dulbecco’s PBS and 100 mM CDNB (≥99%)) was added to each well of a 96-well Nunclon microplate (Thermo Scientific Nunc, ıUnited States), along with 20 μL of GST standard (≥25 units/mg protein, Sigma-Aldrich, United States) or sample. Equine liver GST was used as a positive control to validate the assay. The enzyme activity was determined spectrophotometrically at 340 nm by measuring the formation of the conjugate of GSH (≥99%, Sigma-Aldrich) and CDNB (≥99%, Sigma-Aldrich). The absorbance was recorded every minute for 6 min, using a plate reader (BioRad, United States). The increase in absorbance per minute was estimated and the reaction rate at 340 nm was determined using the CDNB extinction coefficient of 0.0053ϵμM, as follows:

where TV is the total volume, ST is the sample volume and DF is the dilution factor. Results were expressed in relation to the protein content of the samples, which was determined according to the Bradford method (Bradford, 1976).

All data were analyzed using generalized linear mixed models (Zuur et al., 2009). The distributional family used was Binomial (logit link function) for proportions (i.e., survival), Gaussian (identity link function) for quantities (i.e., F_v/F_m and relETR), and Gamma (log link function) for positive quantities with a severe positively skewed distribution (i.e., R and NPP). The sample size of oxidative stress variables was not enough to model HSP and GST as response variables, which were therefore analzsed only through descriptive statistics. The initial mixed models included the species, temperature and pH as fixed effects, the corresponding second and third order interactions, and the tank as a random effect to account for possible dependency within tanks. Following the recommendation from Barr et al. (2013), the random effects were kept in the models irrespectively of the amount of variation they explained.

The most parsimonious models were selected based on the Akaike Information Criterion. Model residuals were checked for departures from the assumed distributions and no significant deviations were found. For Binomial models, odds ratios and confidence limits were determined to allow a more informative discussion of the results. Considering that odds define the ratio of the probability of success and the probability of failure, odds ratios were built by the ratio of odds between the two species (E. crispata vs. E. viridis), temperatures (control temperature vs. warming) or pH (normocapnia vs. hypercapnia).

All statistical analyses were implemented in R, using the lme4 (Bates et al., 2015) and nlme (Pinheiro et al., 2018) packages. Results were considered statistically significant at a significance level of 0.05.

Sea slug survival was significantly affected by temperature (p = 0.005) but not by high CO₂ (p = 0.624) (Figure 1). The odds of survival under control temperature were more than “7 times higher” than the odds of survival under warming conditions. No mortality was observed under control conditions for both species. E. crispata survival decreased under warming conditions to 40 ± 34.6 and 53.3 ± 30.6% under normocapnia and acidification, respectively. E. viridis survival decreased under warming conditions down to 72.2 ± 9.6 and 41.7 ± 8.3%, under normocapnia and acidification, respectively. Moreover, no significant differences were found between species (p = 0.152), although the interaction between species and pH was found to be significant (p = 0.004). While we cannot detect a pH effect over E. crispata survival, E. viridis survival decreased under hypercapnia by 69.4 and 30.6 percentage points, under control temperature and heat conditions, respectively (see more statistical details in Supplementary Table S1).

Under control conditions, E. viridis’ kleptoplasts were packed tightly in the tubule cells surrounding the terminus of the tubule (Figure 2a). Neither high temperature nor high CO₂ affected the color or the morphology of kleptoplasts (Figure 2b). In contrast, E. crispata kleptoplasts were mainly distributed in the tip of the tubule cells under control conditions (Figure 2c). Bleaching was observed in all the slugs exposed to warming condition, with the majority of the host tubule cells being unfilled or displaying degraded kleptoplasts (Figure 2d).

The photosynthetic efficiency of kleptoplasts (Figure 3) was significantly affected by temperature (p < 0.001 for F_v/F_m and relETR) and pH (p = 0.026 for F_v/F_m and p = 0.001 for relETR), but these effects varied between species (p < 0.001 for F_v/F_m and p = 0.003 for relETR). Moreover, the interaction between species and temperature was also significant (p = 0.022 for F_v/F_m and p = 0.030 for relETR). While the F_v/F_m and relETR of E. viridis varied little among treatments, the photosynthetic efficiency of E. crispata decreased under both hypercapnia and heat conditions. More specifically, F_v/F_m decreased 35.9 and 55.3%, while relETR decreased 48.9 and 53.0% under hypercapnia and heat, respectively. However, under the combined effect of hypercapnia and heat, the negative impact of these variables was not cumulative, resulting in a significant interaction between temperature and pH (p = 0.012 for F_v/F_m and p = 0.005 for relETR; see more statistical details in Supplementary Table S1).

Sea slug metabolism was significantly affected by pH (p < 0.001), but not by temperature (p = 0.673) (Figure 4). Moreover, the interaction between species and pH was also significant (p < 0.001 for R and NPP). NPP was significantly affected by both pH (p < 0.001) and temperature (p = 0.043). Moreover, the interaction between species and pH (p < 0.001) and between species and temperature (p = 0.045) was also significant (see more statistical details in Supplementary Table S1). While E. viridis metabolism varied little or even increased, the metabolism of E. crispata decreased to values near zero under acidification and/or warming conditions.

Elysia viridis revealed significantly higher HSP content than E. crispata (Figure 5A). However, when exposed to the combined scenario, E. crispata showed the highest value observed (418.6 μg HSP70/mg of total protein. In E. viridis, such increase was from 232.4 to 385.3 μg HSP70/mg of total protein. Besides the recurrent interspecific differences, GST levels increased 58.8% with warming in E. crispata (Figure 5B); yet, such response was not observed under the combination of both stressors (with a down-regulation of 20% compared to control treatment. Regarding E. viridis, GST levels were significantly higher under both warming and warming+acidification treatments.