Gene expression profiling in gills of the great spider crab Hyas araneus in response to ocean acidification and warming

Abstract

Background

Hypercapnia and elevated temperatures resulting from climate change may have adverse consequences for many marine organisms. While diverse physiological and ecological effects have been identified, changes in those molecular mechanisms, which shape the physiological phenotype of a species and limit its capacity to compensate, remain poorly understood. Here, we use global gene expression profiling through RNA-Sequencing to study the transcriptional responses to ocean acidification and warming in gills of the boreal spider crab Hyas araneus exposed medium-term (10 weeks) to intermediate (1,120 μatm) and high (1,960 μatm) PCO₂ at different temperatures (5°C and 10°C).

Results

The analyses reveal shifts in steady state gene expression from control to intermediate and from intermediate to high CO₂ exposures. At 5°C acid–base, energy metabolism and stress response related genes were upregulated at intermediate PCO₂, whereas high PCO₂ induced a relative reduction in expression to levels closer to controls. A similar pattern was found at elevated temperature (10°C). There was a strong coordination between acid–base, metabolic and stress-related processes. Hemolymph parameters at intermediate PCO₂ indicate enhanced capacity in acid–base compensation potentially supported by upregulation of a V-ATPase. The likely enhanced energy demand might be met by the upregulation of the electron transport system (ETS), but may lead to increased oxidative stress reflected in upregulated antioxidant defense transcripts. These mechanisms were attenuated by high PCO₂, possibly as a result of limited acid–base compensation and metabolic down-regulation.

Conclusion

Our findings indicate a PCO₂ dependent threshold beyond which compensation by acclimation fails progressively. They also indicate a limited ability of this stenoecious crustacean to compensate for the effects of ocean acidification with and without concomitant warming.

Background

Increasing anthropogenic emissions of CO₂ induce ocean warming and acidification. These changes in environmental conditions may have adverse effects on marine organisms [1–5]. However, the responses to ocean acidification (OA) are highly variable between organisms [4, 5] based on the fact that various animals differ in their capacities to compensate for acid–base disturbances caused by elevated seawater CO₂ and resulting blood hypercapnia for review see [5]. Organisms with low compensation abilities show depressed metabolism, altered energy budgets, and as a result, lower rates of growth or development [6–9]. In contrast, organisms compensating for acid–base disturbances through active ion transport, such as fish, cephalopods and some crustaceans are projected to be more tolerant towards OA [3, 10]. In parallel to these differential capacities, sensitivities within a phylum seem to be related to differences in lifestyle and associated energy turnover [5, 11]. Furthermore, species or populations from highly variable environments with natural variations in PCO₂ may have evolved to be more tolerant than species from relatively stable environments. As an extreme example, the shallow living crab Cancer magister can compensate within 24 h for hypercapnia-induced acidosis, while the extracellular acidosis in the deep-sea crab Chionoecetes tanneri remains mostly uncompensated during this time [12]. However, such short-term studies have limited value if it comes to the projection of long-term ocean acidification effects.

The great spider crab Hyas araneus is an osmoconforming, slow-moving and inactive species living in relatively stable physical conditions and is thus an excellent candidate to study the medium to long-term effects of abiotic stressors. A number of physiological studies have already investigated the effects of elevated seawater PCO₂ on this species: CO₂ induced decreases in growth rates and fitness of larvae were demonstrated in a North Sea population, whereas an Arctic population seemed more sensitive towards thermal stress [13]. In the Spitsbergen population elevated seawater PCO₂ (1,100 μatm) caused an increase in metabolic rate during larval development pointing to higher metabolic costs [14]. Adult H. araneus became more heat intolerant under elevated CO₂ with potential consequences for biogeographical distribution [15]. In the Arctic population synergistic effects of increased temperature and PCO₂ adversely influenced the capacities for activity associated with disturbances in acid–base status [16].

To understand organismal sensitivities and tolerance-limits to OA with and without concomitant warming it is important to identify and differentiate between the mechanisms that shape an organism’s capacity to cope with the projected changes. At the whole organism level, crustaceans are impacted by OA with and without concomitant warming with effects ranging from changes in acid–base homeostasis [12, 16, 17], metabolism [9, 18, 19], growth [7, 13, 20, 21], to development [13, 14, 22, 23] and even survival [21, 24]. These processes are highly interdependent. While active acid–base regulation is an energy-consuming process [25], eventually leading to enhanced metabolic requirements [26], uncompensated extracellular pH can elicit metabolic depression [27] via effects on transmembrane ion exchange [25]. Furthermore, low pH can trigger a decrease in protein synthesis [28]; this may result in reduction of growth under hypercapnic conditions [6]. These previous studies provide us with important insights into the mechanistic background of responses to ongoing OA and warming, but also highlight the complexity of the processes involved. To elaborate the sensitivities and potential tolerance limits further, it is important to investigate the key regulatory mechanisms shaping affected processes and the tradeoffs between them.

A transcriptomic approach makes it possible to simultaneously investigate the genetic response of a wide range of cellular processes, and thus to identify the early responses to environmental changes [29]. Gene expression analyses can be used to characterize the molecular phenotype and the cellular changes that underpin physiological responses. They can also be used to uncover molecular mechanisms that might define physiological plasticity. Transcriptomic analysis can further reveal the connections between response mechanisms to environmental changes such as OA or temperature that may otherwise be overlooked [30]. Due to technological advances in recent years, analyses of the whole transcriptome have become increasingly attractive to study non-model (marine) organisms, and their molecular responses to a variety of environmental changes such as warming [31, 32], salinity fluctuations [33], hypoxia [34, 35] or OA [29, 36, 37] in marine organisms.

In the present study, we used gene expression profiling to explore the molecular response in gills of H. araneus exposed to hypercapnia at different temperatures. In marine crustaceans, gills are the first line of defence against acid–base disturbances of body fluids and thus the most important regulatory tissue for CO₂ induced acidification of the hemolymph [38, 39]. We used a quantitative transcriptomic approach based on direct cDNA sequencing using high-throughput Illumina sequencing [40]. Since the present study focuses on the mechanisms involved and the potential sensitivity of H. araneus to climate changes, we selected CO₂ concentrations projected for the year 2100 and 2300 by the Intergovernmental Panel on Climate Change (IPCC) as well as two different temperatures (5°C as the summer control temperature for the Arctic population and 10°C as the median habitat temperature of the species considering the whole distribution range. This study provides comprehensive insights into the transcriptional changes involved in the responses to warming and OA.

Methods

Animals, experimental treatments and tissue sampling

Adult specimens of the Arctic spider crab H. araneus (Linnaeus 1758) were collected by scientific divers in May 2009 in Kongsfjord at the west coast of Spitsbergen at 7–12 m depth (N 78°58.635'; E 11°29.454') and transferred to the Alfred Wegener Institute, Bremerhaven, Germany. Animals were maintained at 5°C in flow through aquaria with natural seawater prior to experimentation. During this period, seawater was aerated with ambient air and animals were fed ad libitum twice per week with frozen mussels and cockles (Mytilus edulis and Cerastoderma edule).

For the medium-term experiment, male spider crabs with a carapace width of 26 to 42 mm were divided into six groups and each group was randomly assigned to the different treatments. Animals were exposed to three different CO₂ concentrations (390 μatm as control, 1,120 μatm as intermediate concentration and 1,960 μatm as high CO₂ treatment) and two different temperatures (5°C as control and 10°C as elevated temperature) for 10 weeks (Figure 1). For each treatment, 5–7 animals were individually placed in 2 l wide-mouth containers (Kautex, Bonn, Germany).

Figure 1

Overview of the experimental design used in the differential expression . Animals were collected in 2009 and acclimated to 5°C until the start of experimentation. Subsequently, exposure experiments were conducted for a time period of 10 weeks for all treatments. After exposure, tissue samples were taken and total RNA was extracted for analyses by quantitative real-time polymerase chain reaction (qRT-PCR) and Sequencing. Sequencing data were used for differential expression analysis by DESeq.

Experiments were carried out in recirculating seawater CO₂ manipulation systems of 1 cubic meter volume each. Seawater of the storage tank was pumped to a header tank at a rate of 20 l min⁻¹, which supplied the wide-mouth container by gravity feed at a flow rate of 200 ml min⁻¹. Water of experimental containers was retained in a collection tank and pumped back to the storage tank at a flow rate of 20 l min⁻¹. Ambient air temperature in the experimental rooms was thermostated to keep water temperature constant. Seawater CO₂ manipulation was accomplished by constantly aerating the storage and header tanks with a defined air/CO₂ mixture using an automatic mass flow controller (HTK 6 channel, HTK Hamburg GmbH, Germany). A light–dark cycle of 12:12 h was established. Water was partly changed every week by refilling the storage tank with pre-equilibrated seawater (PCO₂ and temperature). Experimental animals were fed once a week ad libitum with frozen mussels (C. edule).

To monitor water physicochemistry, seawater samples were collected in airtight glass vials to prevent exchange with the atmosphere, and total dissolved carbon (DIC) concentration was immediately measured with a Seal QuAAtro SFA Analyzer (Seal Analytical, Mequon, United States of America). Temperature, salinity and pH were measured at the time of collection and, together with DIC, used to calculate the PCO₂ in seawater using CO2SYS [41]. Seawater pH was measured using a pH electrode (ProfiLine pH 3310, WTW Wissenschaftlich-Technische Werkstätten GmbH, Weilheim, Germany) calibrated at the respective temperature with National Institute of Standards and Technology (NIST) standard pH buffer and salinity with a conductivity meter (ProfiLine Cond 1970i, WTW Wissenschaftlich-Technische Werkstätten GmbH, Weilheim, Germany). For CO2SYS, NBS (NIST) scale of seawater pH and constants of Mehrbach et al. [41] refitted by Dickson and Millero [41] were used. A summary of water physicochemistry data is given in Table 1.

Table 1 Summary of the seawater physiochemical conditions during experiments with Hyas araneus

After experimental exposure, all 6 gill pairs were collected from 5–7 animals in each treatment. Tissue samples were immediately frozen in liquid nitrogen and stored at −80°C until usage.

Hemolymph sampling and measurements

Directly before tissue sampling, around 1 ml of hemolymph was extracted at the coxa of the third walking leg using a 1 ml sterile syringe (Henke-Sass, Wolf GmbH, Tuttlingen, Germany). Hemolymph was immediately transferred to a 1.5 ml tube (AG Eppendorf, Hamburg, Germany), placed in a thermostatted water bath and pH was measured at acclimation temperature using a pH microelectrode (PHM 93 Reference pH meter, Radiometer, Copenhagen, Denmark; InLab Micro, Mettler Toledo GmbH, Germany). The pH meter was calibrated at the respective temperature with NIST standard pH buffer. A hemolymph subsample was withdrawn using a gas-tight 200 μl syringe (Hamilton Company, Reno, United States of America) and total dissolved inorganic carbon (CCO₂) of extracellular fluid was analysed according to the modified gas chromatographic method [42, 43]. Extracellular fluid was injected in gas tight glass vials containing 3 ml of air equilibrated 0.1 M hydrogen chloride (HCl) and analysed by gas chromatography in an Agilent 6890 N GC System (Agilent Technologies, Santa Clara, United States of America). The bicarbonate (HCO₃⁻) concentration of the extracellular fluid was calculated from CCO₂ and pH using equations derived from the Henderson-Hasselbalch equation. PCO₂ was calculated as PCO₂ = CCO₂ * (10^ph-pkIII * αCO₂ + αCO₂)⁻1 and HCO₃⁻ as HCO₃⁻ = CCO₂- αCO₂ * PCO₂, with CCO₂ being the total CO₂ concentration in mM, αCO₂ the physical solubility of CO₂, PCO₂ the partial pressure of CO₂ in kPa and pK the apparent dissociation constant of the CO₂/apparent HCO₃⁻ system. αCO₂ and pK were calculated according to Pörtner et al. [44]. Raw data of hemolymph sampling and measurements are available at http://doi.org/10.1594/PANGAEA.833705.

RNA extraction and sequencing

Total tissue RNA of gills was extracted using the RNeasy Mini Kit according to the Purification of Total RNA from Animal Tissue protocol (QIAGEN, Hilden, Germany). RNA quantities were determined by a NanoDrop 2000c spectrometer (PeqLab, Erlangen, Germany), and RNA was analysed for quality by microfluidic electrophoresis in an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, United States of America). Total RNA from all gill pairs of 4 animals was pooled for each treatment, and used for library constructions and sequencing by GATC Biotech (Konstanz, Germany). The cDNA libraries for each treatment were constructed according to the SMART protocol for Illumina sequencing (Clontech, Mountain View, USA) and after adapter ligation pooled into two samples. To obtain appropriate deep sequencing results, samples were sequenced at least twice. Illumina sequencing was performed on a HiSeq 2000 Sequencer by GATC Biotech (Konstanz, Germany). Raw reads were quality controlled by FastQC (Babraham Institute, Cambridge, UK) and cleaned using the FastX-Toolkit (Hannon Lab - Cold Spring Harbor Laboratory, New York, USA). Quality control and trimming was performed using the following parameters: Minimum quality score of 20, minimum percentage of bases within the quality score of 90 and a minimum length of 25 bases. The cleaned raw data of the Illumina sequencing were deposited in the European Nucleotide Archive (ENA) at the European Molecular Biological Laboratory – European Bioinformatics Institute (EMBL-EBI) (http://www.ebi.ac.uk/ena/data/view/ERP002128). A summary of the cleaned sequencing results for all samples is given in Additional file 1: Table S1.

Mapping and identification of differentially expressed genes

Short reads of each sample were separately aligned against the annotated H. araneus transcriptome [40], using the Burrows-Wheeler Aligner (BWA) (version 0.5.9) with default parameters [45]. Obtained files were processed into bam files for further analysis, using SAMTools (version 0.1.18) [46]. An overview of the mapping and efficiency is described in Additional file 1: Table S1. Differential expression analysis was conducted with the R statistic software [47]. Read counts were summed up for all sequencing runs of each sample and used for the differential expression analysis without biological replicates. Differential expression of genes was evaluated using a test based on the negative binomial distribution as integrated in the Bioconductor R package DESeq [48], with a standard level of p ≤ 0.05 indicating significance. Control (control PCO₂/control temperature) was compared to five treatments: (I) elevated temperature; (II) intermediate PCO₂ at control temperature; (III) high PCO₂ at control temperature (IV) intermediate PCO₂ at elevated temperature; (V) high PCO₂ at elevated temperature. The previously annotated transcriptome made a Gene Ontology enrichment analysis possible to test for particular affected terms, using Fisher’s exact test (FDR ≤ 0.05) as implemented in the Blast2GO software (version 2.6.0) [49, 50]. All subsets of significantly regulated genes identified by the binominal distribution test were tested against the full set of annotated sequences of the H. araneus transcriptome. To cut down on redundancy, GO terms were summarized into a more representative subset of terms using the web-based clustering tool REVIGO [51].

Validation by quantitative real-time polymerase chain reaction (qRT-PCR)

A set of transcript sequences known to be involved in acid–base regulation and/or transcripts that showed differential expression in one or more treatments was selected for validation of RNA-Seq results. Primers were designed using the PrimerExpress software (version 3.0) (Applied Biosystems, Darmstadt, Germany) with the Taq-Man MGB Quantification method and default parameters (Additional file 2: Table S2). Primer specificity was given by using sequences of the annotated H. araneus transcriptome [40]. All primer pairs were tested for performance and efficiency across a series of cDNA dilutions (1:20; 1:40; 1: 100; 1:200; 1:1000; 1:2000). Primers used displayed a suitable per cycle amplification rate, with an efficiency (E) of 2.0 ± 0.1 and R² > 0.98. Efficiency was calculated as E = 10^(−1/S), with s being the slope of linear regression.

Total RNA was extracted from gills as described above. Ten micrograms of total RNA per sample was treated with DNAse for DNA digestion using the Turbo DNA-free kit (Ambion, Darmstadt, Germany) and 0.4 μg DNA free RNA was transcribed into cDNA with the High-Capacity cDNA Reverse Transcription kit (Applied Biosystems, Darmstadt, Germany). Real-time PCR was performed on a 7500 Real-time PCR System (Applied Biosystems, Darmstadt, Germany) and SYBR® Green PCR master mix (Applied Biosystems, Darmstadt, Germany). All genes were finally analysed in a 40-fold dilution and amplified with 300 nM of primer. To verify the amplification specificity of fragments a melting curve analysis was performed for each reaction.

Gene expression calculation was based on the C_T-threshold. Absolute mRNA quantities were calculated as Q_X = E^(CT) and normalized with the formula Q_N(X) = Q_X/Q_X(HK), with Q_X(HK) being the absolute mRNA quantity of the housekeeping gene sodium bicarbonate cotransporter (NBC). The housekeeping gene was determined using geNorm implemented in the software qbasePlus (version 2.1) (Biogazelle, Zwijnaarde, Belgium) with a relative expression stability of M ≤ 0.42 (high reference stability is given at an average geNorm of M ≤ 0.5). To ensure consistency with the differential expression results of the DESeq analysis, gene expression results of the qRT-PCR were calculated as log₂ fold change (log₂FC) of mean normalized quantities of treatment and control.

Statistics

To identify significant differences in the sum of all significantly up- and down-regulated transcripts between treatments, significantly changed transcripts in one or more treatment, identified by the DESeq analysis, were transformed into a matrix with 1 = significantly up-regulated, −1 = significantly down-regulated and 0 = not significantly regulated transcript. Treatments were analysed for statistical differences applying the Wilcoxon matched pairs test as implemented in SigmaPlot (Version 12.0, Systat Software Inc., San Jose, USA) with p < 0.05. Data from each treatment were tested against each other.

The correlation between the differential expression results of the DESeq analysis and the corresponding gene expression results of the qRT-PCR was determined by Pearson Correlation as implemented in SigmaPlot 12.0 (Systat Software Inc., San Jose, USA).

A one-way ANOVA was used to identify the effect of seawater PCO₂ on hemolymph pH and bicarbonate (HCO₃⁻). Data obtained under various PCO₂ levels were tested against each other for each temperature separately. A Holm Sidak test for multiple comparisons was used for a posteriori analyses. Tests were performed in SigmaPlot (Version 12.0, Systat Software Inc., San Jose, United States of America) with p < 0.05 indicating significant differences.

Results and discussion

A total of 55 million reads (56%) from initial Illumina sequencing passed the quality filter and was used for the differential expression analysis. After processing, an average of 9.2 million high quality reads were produced for each sample from 2–3 sequencing runs per sample (Additional file 1: Table S1). To obtain the differential expression of each gene, high quality reads were aligned on the H. araneus transcriptome [40]. An average of 5.2 million reads for each sample produced distinct alignments. The alignment process yielded an average efficiency of 56% for the high quality reads (Additional file 1: Table S1). The achieved mapping efficiency is actually higher than in a comparable study of a non-model organism, which used an analogous approach for differential expression analysis (41% efficiency) [52]. Furthermore, 96.5% of all transcript sequences in the H. araneus transcriptome were detected in the RNA-Seq data. However, the occurrence of a large amount of unmapped reads might result from sequencing errors, repetitive sequences or inadequate quality filtering of the Illumina reads. Furthermore, it has to be considered that transcripts supported by only a small number of aligned reads (≤10) may reflect incompletely assembled transcripts in the reference. Those poorly supported transcript sequences were excluded from the subsequent analysis, and a final test-set of 16,201 transcripts sequences was used for the differential expression analysis. As there were no biological replicates, it has to be considered that the variance for genes can only be estimated by comparing the mean-variance relationship between samples/treatments, as if they were replicates, resulting in an overestimation of the variance and thus make this approach more conservative. Furthermore, expression levels of stable and highly expressed genes, based on the RNA-Seq data, were analysed by quantitative real-time PCR (qRT-PCR) confirming the RNA-Seq methodology used in this study (Additional file 3: Figure S1 and Additional file 4: Figure S2).

We could identify 864 (5.3%) out of the 16,202 tested transcripts to be differentially expressed after medium-term (10 weeks) exposure to the abiotic effectors (Additional file 5: Table S3). Out of these differentially expressed genes (DEG) 40.0% and 31.3% were differentially expressed under intermediate CO₂ (treatment I) and high CO₂ (treatment II; Additional file 6: Figure S3A, B), respectively. For the high temperature (treatment III), 41.0% showed significantly different expression levels (Additional file 6: Figure S3C). The combination of factors intermediate CO₂ (treatment IV) and high CO₂ (treatment V) with elevated temperature revealed 38.7% and 29.4% DEG, respectively (Additional file 6: Figure S3D, E). While the total amount of significant DEGs was similar in all treatments, individual genes displayed large differences in up- or down-regulation (Figure 2). Intermediate CO₂ (I) led to strongly up-regulated transcript levels, while high temperature (III) caused strong down-regulation.

Figure 2

Differences in the sum of all significantly regulated transcripts of Hyas araneus after exposure experiments. For each treatment, up- and down-regulated transcripts were counted and represented as bars (black bars = up-regulated; grey bars = down-regulated). Treatments were tested for statistical differences using the Wilcoxon signed rank test (p < 0.05). Letters denote the significant differences. Differing letters indicate significant differences; identical letters indicate no significant difference.

These general results, especially the strong up-regulation of gene expression under intermediate CO₂ (I), suggest that the regulatory capacity is high at moderately elevated PCO₂ in adult H. araneus but reduced at higher CO₂ levels. Up-regulation was also reduced at intermediate CO₂ when combined with elevated temperature (IV) rather than at control temperature. Down-regulation in DEG predominated under warming alone (treatment III, high temperature) when compared to control temperature (I). 137 of the 177 genes that displayed significantly changed transcripts in the intermediate CO₂/elevated temperature treatment (IV) were also differentially expressed in the high temperature treatment (III) supporting a strong temperature-dependent response, in line with the hypothesis that temperature affects most biochemical processes. Faced with a long-term temperature change, ectothermal organisms retain physiological homeostasis by several acclimation strategies, which can be of quantitative, qualitative or modifying nature [53]. The enhanced down-regulation seen in the high temperature treatments suggest that H. araneus may adopt a quantitative strategy to maintain physiological rates by down-regulating the concentrations of enzymes between 5 and 10°C [54].

To identify processes actually responding to elevated PCO₂ and temperature, a first analysis was carried out using gene ontology (GO) terms. The set of GO-annotated differentially expressed genes was statistically tested for the over- and underrepresentation of GO terms to identify molecular functions, cellular components and biological processes affected most by the experimental treatments (Additional file 7: Table S4). The GO enrichment analysis revealed a variety of significantly over-represented GO terms that can primarily be summarized underneath the more generic categories ‘metabolism’ and ‘cell structure’. 23 and 25 over-represented GO terms, respectively, could be associated with these categories. It is important to mention that both intermediate CO₂ treatments (I, IV) constitute 57% and 27% of the over-represented GO terms, respectively, representing compensation mechanisms mentioned above. The only over-represented GO term under high CO₂, trehalose metabolic process, was evaluated for the down-regulated genes of the combined stress treatment (V). The fact that it was over-represented in all treatments indicates that enhanced expression of trehalose metabolism can be rated as a unifying response to both elevated CO₂ and elevated temperature.

Within the ‘cell structure’ related GO-terms, intermediate CO₂ (I, IV) led to significant up-regulation of genes concerning cell surface, extracellular matrix, structural molecule activity or brush border membrane, suggesting a structural modification of the gills. Some GO terms were found underrepresented at high temperature and in combined, intermediate CO₂ and high temperature treatments (III, IV), and are assigned mainly to intracellular structures such as organelles. The over-representation of ‘cell structure’ related GO terms suggests that gill epithelial structure is adjusted in response to PCO₂ disturbances. Gills are the principal organs for gas exchange and, together with the excretory organs, responsible for osmotic and ionic regulation in crustaceans [55]. As passive ion transport is influenced by the conductivity of gill epithelia [56], their structural modification might lead to a change in conductivity and would change the diffusion rate of ions. Structural changes were in fact identified in gills of Carcinus maenas during salinity exposure, with a modification of the apical plasma membrane system and an enlargement of the subcuticular compartment [57]. As environmental hypercapnia and salinity changes can cause similar mechanistic responses [39], similar transcriptomic modifications may occur. This is supported by an over-representation of the GO term response to salt stress in the intermediate CO₂ treatment (I).

Although a GO enrichment analysis offers initial insights into processes affected by hypercapnia exposure and elevated temperature, a strong bias exists towards conserved and well-characterized processes, functions and cellular components in model organisms. This bias particularly applies to H. araneus, with a lack of GO annotation for about 76% of the transcripts. Additionally, many genes are grouped into more than one GO term depending on their resolution and are thus difficult to interpret. In light of these contraints, GO analysis can only provide a general overview of possibly affected processes and make a more detailed look indispensible.

For a more comprehensive understanding of the mechanisms responding to PCO₂ and temperature changes, we performed a second analysis. Here, all genes included in the most affected categories ‘metabolism’ and ‘structural modification’, identified by the GO enrichment analysis, were considered. Additionally, we integrated all genes related to ‘acid–base and ion regulation’ and ‘response to stress’ into our analysis, as adjustments in these mechanisms are likely relevant in shaping resistance to hypercapnia exposure or heat stress.

For the interpretation of transcriptomic results, it has to be considered that expression profiles represent one regulatory level in the response to environmental changes and thus do not necessarily reflect the changes of other regulatory levels, e.g. protein abundances and activities. At least 25% of the proteome cannot be covered by gene expression profiling [30]. However, there was a positive correlation between transcription and translation in 87% of genes that changed ≥ twofold in living cells of yeast strains [58]. Certainly, a lower correlation has to be considered for more complex organisms, however an additional proteomics study on H. araneus revealed a comparable correlation between gene expression and protein abundance of ~70% for the intermediate treatment at 5°C (unpublished). A transcriptomic approach can provide first insights into the regulatory processes responding to environmental changes such as OA or temperature.

Response of specific groups

Acid–base regulation

The extracellular pH (pH_e) measured in the hemolymph of adult H. araneus showed partial compensation under intermediate CO₂ levels (I, IV) involving an increase in bicarbonate (HCO₃⁻) concentration (Figure 3A,B). Under high CO₂ (II, V), reduction in pH_e was greater and the increase in ${\mathrm{HCO}}_3^{‒}$ was reduced (Figure 3A,B). These findings suggest a limited capability to compensate for pH_e disturbances caused by high seawater PCO₂. According to a crustacean model by Freire et al. [55], proton (H⁺) excretion is generated by apical vacuolar-type (H⁺)-ATPase (V(H⁺)-ATPase) and/or sodium/proton exchanger (NHE), the latter dependent on sodium/potassium-ATPase (Na⁺/K⁺-ATPase). HCO₃⁻ is enriched in the hemolymph by basolateral anion exchangers. Enzymes that support active ion transport are the intracellular carbonic anhydrase (CA) and, in terms of a general support of energy consuming mechanisms, arginine kinase (AK). CA is assumed to accelerate the dissociation of carbonic acid (H₂CO₃) and provide the substrate for H⁺ and HCO₃⁻ transporters [55]. AK catalyses the reversible dephosphorylation of phosphoarginine, contributing to the restoration of adenosine triphosphate (ATP) used in energy consuming processes [59]. The expression of corresponding genes, V(H⁺)-ATPase, AK and partial sequences of two alpha CAs was significantly up-regulated at intermediate CO₂, whereas such mRNA concentrations were only moderately increased at high CO₂ (II) (Table 2). These expression levels follow the course of the hemolymph HCO₃⁻ parameters and the more effective acid–base regulation of adult H. araneus under moderately elevated CO₂. However, only a few sequences encoding for Na⁺/K⁺-ATPase were up-regulated under CO₂, in contrast to their response to elevated temperature treatment (III). Under the combined effect of temperature and CO₂ Na⁺/K⁺-ATPase was also up-regulated. Another enzyme, DOPA decarboxylase, which catalyses the biosynthesis of dopamine by decarboxylation of L-DOPA, was found down-regulated at moderate (intermediate) CO₂ elevations (I). Elevated dopamine leads to increased sodium influx and concomitantly, increased Na⁺/K⁺-ATPase activity in gills of C. maenas[60]. In light of this finding constant mRNA levels of Na⁺/K⁺-ATPase at intermediate CO₂ (I) indicate that the down-regulation of DOPA decarboxylase may even prevent an activation of Na⁺/K⁺-ATPase (Table 2). The gene expression of other transporters, supposed to be involved in acid–base regulation, such as NHE and/or bicarbonate/chloride co-transporter were not influenced at all by elevated PCO₂ values or temperature. Acid–base regulation predominantly via the V(H⁺)-ATPase (as seen in the present transcriptome) might involve minimal disturbance to ionic composition, e.g. cellular sodium homeostasis which would be affected by strong involvement of Na⁺/K⁺-ATPase. In the sipunculid S. nudus, extracellular acidosis induced a shift in ion transporters during hypercapnia from high to low energy cost mechanisms of acid–base regulation, resulting in decreased Na^+/K⁺-ATPase activity due to lower requirement for sodium regulation [25, 27]. In line with our present observations these results suggest a pH regulation system independent of Na⁺/K⁺-ATPase being used under hypercapnic exposure (Figure 4). No significant up-regulation of apical bicarbonate anion exchangers was observed. In C. sapidus, CO₂ induced acid–base disturbances (PCO₂ 10,000 μatm) were mainly compensated by an uptake of HCO₃⁻ from the surrounding seawater via the gill epithelia [39]. However, our results with no up-regulation of respective transporters suggests hemolymph buffering may be achieved through the dissociation of respiratory CO₂ resulting in an accumulation of HCO₃⁻ in the hemolymph. This might be exceptionally for crustaceans with low compensatory capacities. In H. araneus, also no significant up-regulation of basolateral bicarbonate anion exchangers was observed. However, expression of the anion channel bestrophin was significantly enhanced (Table 2, Figure 4). Bestrophin is activated by calcium (Ca²⁺) and enhances membrane permeability for anions, such as chloride and HCO₃⁻[61]. Even though bestrophin is commonly found in the retinal pigment epithelium, ion channels of this family have been identified in different tissues and are involved in a variety of cellular processes [62]. Its role in responses to elevated CO₂ remains to be investigated.

Figure 3

Hemolymph acid–base status Hyas araneus (A, pH-values; B, bicarbonate levels; C, P CO ₂ levels) in response to different exposure experiments. Squares represent the means with error bars depicting the standard error for each treatment. Grey squares refer to treatments at 5°C and black squares to treatments at 10°C. One-way ANOVAs were used to identify the effect of seawater PCO₂ concentration on hemolymph pH, bicarbonate (HCO₃⁻) and PCO₂. A Holm-Sidak test for multiple comparisons was used for posteriori analysis (p < 0.05). Differing letters indicate significant differences; identical letters indicate no significant difference. Capital letters denote differences for 5°C treatments and lower cases for 10°C treatments.

Table 2 Regulation of transcripts of specific interest in the hypercapnia and elevated temperature experiments on Hyas araneus

Figure 4

Schematic description of proposed processes in the gill epithelium of Hyas araneus in response intermediate P CO ₂ exposure. Medium-term hypercapnia acclimation leads to a shift to a new acid–base equilibrium by accumulation of hemolymph bicarbonate (HCO₃⁻). CO₂ is hydrated into H⁺ and HCO₃⁻ by cytoplasmic carbonic anhydrase (CA). Protons are actively pumped out of the epithelial cell by an apical vacuolar proton ATPase (V(H⁺)-ATPase), followed by a transport of HCO₃⁻ via a basolateral anion exchanger (AE) and/or ion channel, such as bestrophin (Best1). Increased energy demand is in part met by an enhanced expression of complex I (CI) and complex IV (CIV) of the electron transport system and possibly triggered by a soluble adenylyl cyclase (sAC) induced signalling pathway. sAC is stimulated by HCO₃⁻ and increases the formation of cyclic adenosine monophosphate (cAMP), which activates protein kinase A (PKA) that subsequently leads to an induced expression of CI and CIV. While enhanced aerobic metabolic processes increase the generation of oxidants, the cellular defence against oxidative stress is adversely affected by a lower production of NADPH due to a decrease of the pentose phosphate pathway enzyme glucose-6-phosphate dehydrogenase (G6PD). NADPH acts as reducing agent for the regeneration of reduced glutathione (GSH) to oxidised glutathione (GSSH). The symbols (+) and (−) mark significantly up- and down-regulated genes, respectively modified after [55, 64].

Energy metabolism

Acid–base regulation compensating for CO₂ induced disturbances likely entails a metabolic cost due to the ATP demanding movement of H⁺ by ion transporters [25]. Elevated metabolic rates have already been observed under hypercapnia exposure [26, 63], however, the mechanisms causing such cost increments remain unidentified. The differential expression analysis revealed an up-regulation of complex I and complex IV of the electron transport system (ETS) during exposure to intermediate CO₂ levels (i,iv) (Table 2). Transcripts encoding for NADH dehydrogenase and cytochrome c oxidase subunits were significantly up-regulated above control levels possibly leading to an increase in mitochondrial density to meet the increased energy demand. Interestingly, complex II and III were not affected by hypercapnia exposure. Increased activities of complex I and IV activated by protein kinase A (PKA) resulted in an increased oxidative phosphorylation and ATP synthesis in human kidney cells [64]. The activation was triggered by a soluble adenylyl cyclase (sAC) induced signalling pathway that implies phosphorylation of PKA being stimulated by cyclic adenosine monophosphate (cAMP), which, in turn, is formed from sAC. In mammals and in elasmobranchs it is known that sAC is stimulated by HCO₃⁻[65, 66], which would thereby support an increased activity of mitochondrial electron transport system under hypercapnia exposure. cAMP and PKA regulate key enzymes, such as complex IV, by alterating gene expression [66, 67]. Considering the higher extracellular bicarbonate levels in animals under hypercapnia (Figure 4) this signalling pathway possibly led to the increased up-regulation in ETS related genes and/or increased activity and thus might explain how the organisms meet the suggested increase in ATP demand (Figure 4).

An increased ATP production by the ETS would also lead to an elevated demand for metabolic substrates and turnover of the resulting reduction equivalents (NADH, FADH₂). In the intermediate CO₂ treatments (I, IV), the glycolytic pathway only experienced a significant up-regulation of glyceraldehyde-3-phosphate dehydrogenase (GAPDH) (Table 2), whereas the rate-limiting enzymes phosphofructokinase or pyruvate kinase remained unaffected, suggesting no general up-regulation of glycolysis. However, the enhanced expression of sodium glucose transporter indicates increased capacity for glucose transport from the hemolymph into the gill cells (Table 2). Trehalose is the major hemolymph sugar in insects and decapod crustaceans, with higher levels than glucose [68, 69]. Enzymes for trehalose synthesis were found in crustacean tissues, including gills [70]. Trehalose and its fast transport into cells and consecutive transformation into glucose reflect its immediate availability to meet sudden bouts of energy demand. Accordingly, a 6.7 fold increase of trehalose concentration was measured in the hemolymph of C. maenas over 10 days under osmotic stress [71]. In H. araneus, however, the expression of trehalose-6-phosphate synthase was significantly down-regulated in all treatments (I-V) suggesting suppressed synthesis of the already depleted trehalose stores (Table 2). Significant up-regulation of transcripts encoding for glycogen phosphorylase (I, IV) and alpha glucosidase (IV) was found instead. Both enzymes catalyse the glucose releasing steps of glycogenolysis indicating the use of glycogen as a glucose source during long-term increased demand.

Enhanced demand for glucose is paralleled by a down-regulation of glucose-6-phosphate dehydrogenase (G6PD) in all treatments (I-V) indicating the potentially reduced production of NADPH (Table 2). G6PD is the key enzyme of the oxidative phase of the pentose phosphate pathway, the main source of NADPH for biosynthetic pathways in the cells (Figure 4). Furthermore, the cytosolic (NADP dependent) isocitrate dehydrogenase (IDH) was significantly down-regulated under intermediate CO₂ levels at high temperature (IV), also suggesting lowered biosynthetic rates such as lipid biosynthesis under combined exposure (Table 2).

Oxidative stress

Besides being involved in lipid biosynthesis, NADPH is an important reducing agent in cellular antioxidative defence, e.g. by regenerating reduced glutathione, a major cellular antioxidant (Figure 4). Thus, besides a general down-regulation of anabolic reactions it seems conceivable that H. araneus encounters a reduced capacity to counteract oxidative stress under hypercapnic and thermal stress.

Significant changes in the expression of several genes involved in cellular antioxidant defence, including several peroxidases, indicate potential oxidative stress in H. araneus under intermediate CO₂ exposure (I). These changes in expression level were less pronounced under high CO₂ (II), suggesting a decreasing acclimation capacity of the H. araneus with increasing external PCO₂(Table 2, see above). Among up-regulated genes, especially genes associated with the detoxification of hydrogen peroxide (H₂O₂) were affected under hypercapnia exposure. An ascorbate peroxidase (APX) was significantly up-regulated under intermediate CO₂ concentrations (I). APX is a peroxidase that utilizes ascorbate as electron donor to detoxify H₂O₂ into water (Table 2). Additionally, a glutathione peroxidase (GPX) that reduces H₂O₂, using glutathione as substrate, was up-regulated under high CO₂ exposure (II). This contrasts the down-regulation of two thiol-specific peroxiredoxin-1, at intermediate CO₂ concentration (I), indicating a balanced response possibly by differential transcription of different splice variants. Peroxiredoxins are ubiquitous enzymes detoxifying peroxides, such as H₂O₂, by oxidising their active cysteine site using peroxide as substrate and are regenerated by oxidation of a thiol-containing electron donor, commonly thioredoxin [72]. However, only two sequences encoding for a thiol-containing protein (thioredoxin-1) were significantly down-regulated (II, IV).

The up-regulation of genes for anti-oxidants, such as glutathione peroxidase and peroxiredoxin, may indicate compensation for enhanced ROS (reactive oxygen species) production and concomitantly oxidative stress in the gill tissue of H. araneus under CO₂ exposure. This is further supported by a significant up-regulation of a ribosomal cytochrome p450 like gene in the intermediate CO₂ treatment (I). Cytochrome p450 is involved in the oxidative metabolism of a variety of organic substrates and incomplete catalytic processes can result in a continuous release of ROS [73, 74]. In contrast, a transcript encoding for urate oxidase (uricase) was significantly down-regulated in all treatments (I-V) (Table 2). Uricase is commonly located in the peroxisomes of the hepatopancreas tissue, however, uricase activity has also been detected in gill tissue of the kuruma shrimp Marsupenaeus japonicas[75]. Uricase catalyses the reaction from urate to allantoin and contributes to the generation of H₂O₂ by the oxidation of uric acid [76]. Consequently, a down-regulation of uricase may contribute to alleviate the generation of ROS. Interestingly, CO₂ exposure also led to an increase in several vitellogenin like transcripts (Table 2). Although expression of vitellogenins is generally sex- and tissue-specific, the expression in both sexes of the mud shrimp Upogebia major[77], revealed a positive effect on oxidative stress resistance regardless of the developmental stage. Vitellogenin is also beneficial for oxidative stress resistance in honeybees, Apis mellifera[78]. Even if the function of vitellogenin in oxidative stress resistance is far from being completely understood and further investigations are needed to validate this hypothesis, the strong up-regulation of vitellogenin under CO₂ (I, II, IV) may indicate that vitellogenin is an important protein in the resistance to CO₂-induced oxidative stress (see below).

Hypercapnia induced enhancement of oxidative stress defence was recently demonstrated in the Eastern oyster, Crassostrea virginica[79]. In a proteomic approach, an up-regulation of several proteins, e.g. superoxide dismutase and several peroxiredoxins, was detected after exposure to high PCO₂ (~3,520 μatm) for 2 weeks. The authors suggested several ways of how increased CO₂ levels could directly or indirectly cause oxidative stress. On the one hand, a reaction of CO₂ with peroxynitrite, a ROS formed through the reaction between superoxide anions and nitric oxide, resulting in the formation of reactive carbonate and nitrogen species, can lead to oxidative stress by oxidizing molecular compounds [80]. On the other hand, an indirect influence of elevated CO₂ and/or pH could adversely affect mitochondrial functions and/or the non-enzymatic production of ROS [79]. However, our findings cannot support a direct influence of CO₂ on oxidative stress generation, as one would assume an increase in the response to oxidative stress with increasing seawater CO₂. Our data rather indicate an indirect influence of elevated CO₂ on ROS production, likely by enhanced ROS production due to metabolic stimulation. The suggested increase in oxidative metabolic processes might cause enhanced ROS production and would also explain why the oxidative stress response was higher under intermediate CO₂ than in the high CO₂ treatment.

Cell structure

It is well known that the formation of ROS can damage lipids, proteins and DNA e.g. [81]. A large group of genes that belong to the functional category ‘cytoskeleton’ include several actins and tubulins which are up-regulated together with the antioxidant genes under CO₂ exposure. The cytoskeleton is one major target for oxidative stress (Table 2) when the exposed cysteine component of actin forms oxidized derivates, such as intermolecular disulfide bridges [82]. This presumably has adverse effects on the interaction between actin and actin binding proteins and leads to changes in the structure of the actin cytoskeleton. The up-regulation of two transcripts encoding for actin in both intermediate CO₂ treatments (I, IV) and the high CO₂ treatment (II) may counter the damages caused by oxidative stress. An additional up-regulation of two actin binding proteins, nesprin and adducin, further supports the need for structural adaptation under oxidative stress triggered by CO₂ exposure (Table 2). There is strong evidence that ROS induce the expression of tubulin [83], which could explain the up-regulation of β-tubulin and three transcript sequences encoding forα-tubulin in both intermediate CO₂ treatments (I, IV) and the high CO₂ (II) treatment (Table 2). Although the effect of oxidative stress on the cytoskeleton is still poorly understood and needs further investigation, an interaction of oxidative stress and adaptive changes in the cytoskeleton is well recognized [79, 82] and supported by our findings.

Besides enhanced antioxidative defence, a reorganisation of the cytoskeleton may occur in gill epithelia during the initial stage of hypercapnia acclimation. In C. maenas, a reorganisation of gill epithelia was observed after short-term exposure to very high PCO₂ (~4,340 μatm; 7 days), but not after medium-term exposure to high PCO₂ (~2,270 μatm; 11 weeks) [36]. Similarly, a reorganization of the cytoskeleton may not occur during medium-term CO₂ acclimation of H. araneus (10 weeks). The response of genes that belong to the functional category ‘cytoskeleton’ was higher at intermediate CO₂ (I, IV) than at high CO₂ (II), in line with the findings in antioxidative defence. This also suggests that the changes observed result from oxidative stress via enhanced metabolic rate at intermediate CO₂ with the cytoskeleton as a direct target of ROS formation.

Cytoskeletal genes also responded to elevated temperature (III). Several genes encoding for actin-binding proteins as well as membrane and cuticle proteins were significantly up-regulated (Table 2). Three transcripts of gelsolin, an actin-binding protein and regulator of actin filament assembly, were significantly up-regulated. The expression of an actin-binding nesprin, a membrane-binding ankyrin and an actin-related cytoskeletal structure protein spectrin was also increased. Warming affects various cellular processes, including stability and/or dynamics of the cytoskeleton for review see [84]. An increase in cytoskeletal gene expression despite the general down-regulation of expression levels indicates a requirement for the stabilization of the cytoskeleton, possibly elicited by the warm-induced stimulation of metabolic rate and associated ROS formation. Interestingly, such up-regulation of cytoskeletal gene expression could not be seen in under combined intermediate PCO₂ and high temperature (IV, V), suggesting that the general down-regulation induced by high CO₂ also affects cytoskeleton assembly.

Commonly, stresses such as temperature extremes, cellular energy depletion, and extreme concentrations of ions, osmolytes and gases induce the synthesis of heat-shock proteins (HSP) for review see [85]. However, HSP expression levels being unchanged in the high temperature treatment (III) suggest that 10°C is not a temperature extreme for the Arctic H. araneus population. In contrast, the expression of several HSPs was up-regulated under CO₂ exposure (Table 2). One HSP70 was up-regulated in both the intermediate and the high CO₂ treatments (I, II, IV). Four transcripts encoding for HSP90 were significantly up-regulated in the intermediate CO₂ treatment (I) and two in the intermediate CO₂ plus high temperature treatment (IV). In treatments I, III, and V one HSP90 showed decreased expression. Acting as molecular chaperones, HSPs can play an important role in maintaining proteins in a folded or unfolded state, controlling the unintended aggregation of proteins or target proteins for degradation [86]. Their induced expression, especially under intermediate CO₂ concentrations, paired with increases in mRNA transcripts for proteins involved in antioxidant defence suggest an increased capacity to defend the cell against cellular damage. The up-regulation specifically of HSP90 suggests HSP90 to be of particular importance in maintaining protein homeostasis.

Conclusions

Based on a comprehensive expression analysis of genes related to acid–base regulation, metabolism, cell structure and their coupling to the stress response, this study has identified moderate, but distinct responses to ocean acidification in gills of adult H. araneus. We could demonstrate that the molecular response strongly depends on the CO₂ concentration. At PCO₂ values proposed for the end of the century, changes in expression suggest elevated metabolism caused by stimulated acid–base regulation and associated with increased oxidative stress. This up-regulation was attenuated at even higher PCO₂ (business as usual scenario for Year 2300) with expression levels closer to control values (Figure 5), indicating down-regulation of these processes. Responses in transcripts related to acid–base regulation and energy metabolism were in line with those observed at elevated temperature. Observations under hypercapnia in the blue mussel M. edulis are in line with these conclusions. An increase in aerobic metabolism under intermediate PCO₂ exposure (<2,400 μatm) was followed by a decrease at higher PCO₂[87, 88]. In earlier studies on marine invertebrates, high PCO₂ exposure (~3,200-5,000 μatm) induced a decrease in metabolic rates. In the mussel Mytilus galloprovincialis and the decapod crab Necora puber, a decrease in respiration rates was seen in response to long-term exposure at a PCO₂ of approx. 5,000 and 3,200 μatm, respectively [6, 9]. Our findings suggest a PCO₂ dependent threshold where metabolic stimulation might turn into metabolic depression. In H. araneus this threshold may be reached at a PCO₂ of approximately 2,000 μatm. Although further examinations besides transcriptomics covering a broader range of PCO₂ concentrations are certainly needed to confirm a possible threshold, a recent meta-analysis revealed that 50% of all marine crustacean species were negatively affected at this PCO₂[5].

Figure 5

Schematic description of the response of Hyas araneus to different exposure experiments for specific categories. For each treatment, log₂-fold changes of all up-regulated transcripts within one specific category were summed and plotted as a line (gray line = 5°C; black line = 10°C). Scales are harmonized for all categories. The figure shows how acclimation to warming reduces the transcriptomic response of ‘cell organization’ and ‘response to stress’ to elevated CO₂ tensions.

Changes in the expression of genes related to energy metabolism paralleled the expression of those related to acid–base regulation and both paralleled the changes in extracellular pH (pH_e), indicating a feedback regulation between pH and metabolic rate (Figure 5). In the marine worm S. nudus and the mussel M. galloprovincialis, metabolic depression under CO₂ was caused by a drop in pH_e[6, 27, 89]. H. araneus was able to partially compensate for the acidosis via bicarbonate (HCO₃⁻) accumulation, but capacities were not sufficient to fully compensate. At both experimental temperatures (5°C and 10°C) hemolymph acidosis was compensated to the same extent. At high CO₂ concentrations, the capacity to fully compensate may be even more limited related to decreased ion transporter expression and activities. Since acid–base regulation is an ATP-consuming process, the decreased expression of metabolism related genes might at least in part be explained by a pH_e-triggered inhibition.

These shifts cause changes in energy balance which can therefore reflect the limits in stress tolerance [90]. A bioenergetic framework [90] builds on the physiological concepts of oxygen- and capacity-limited thermal tolerance (OCLTT) [91–93] and the dynamic energy budget (DEB) [94]. Two cases distinguished between physiological responses to moderate and extreme environmental stress. Moderate stress induces additional ATP turnover along with a higher metabolic rate, in H. araneus exemplified in enhanced acid–base regulation and cellular damage repair from oxidative stress at intermediate PCO₂. These mechanisms support long-term persistence at the expense of a shift in energy budget (pejus range) [95]. Even if the enhanced energy demand of stress resistance can be met by increased feeding, it may still result in a reallocation of energy from fitness-related functions such as reproduction and growth to maintenance and damage repair [90]. More extreme stress exacerbates the disturbances in homeostasis, exemplified in stronger disturbances in acid–base status as seen in H. araneus favouring a suppression of metabolic rate [27] which preserve energy resources and lessen the generation of detrimental metabolic by-products [90]. In H. araneus the respective response to high PCO₂ was paralleled by a lower expression of transcripts associated with energy metabolism, stress response and cellular organization. Depressed metabolism would result in insufficient energy balance and time-limited tolerance to environmental stress making long-term survival impossible (pessimum range) [90, 93].

To our knowledge, we report the first gene expression profiling study analysing the responses to a combination of two drivers in an osmoconforming crustacean. The study confirms the interdependence of physiological processes affected by elevated seawater PCO₂ and temperature. Furthermore, the study demonstrates the importance of considering projected climate scenarios in experimental work, as responses to increasing seawater CO₂ concentrations are not necessarily linear, as presented here.