Molecular analysis of the Sydney rock oyster (Saccostrea glomerata) CO₂ stress response

Abstract

Background

Human activities have led to a substantial increase in carbon dioxide (CO₂) emission, with further increases predicted. A RNA-Seq study on adult Saccostrea glomerata was carried out to examine the molecular response of this bivalve species to elevated pCO₂.

Results

A total of 1626 S. glomerata transcripts were found to be differentially expressed in oysters exposed to elevated pCO₂ when compared to control oysters. These transcripts cover a range of functions, from immunity (e.g. pattern recognition receptors, antimicrobial peptides), to respiration (e.g. antioxidants, mitochondrial respiratory chain proteins) and biomineralisation (e.g. carbonic anhydrase). Overall, elevated levels of CO₂ appear to have resulted in a priming of the immune system and in producing countermeasures to potential oxidative stress. CO₂ exposure also seems to have resulted in an increase in the expression of proteins involved in protein synthesis, whereas transcripts putatively coding for proteins with a role in cilia and flagella function were down-regulated in response to the stressor. In addition, while some of the transcripts related to biomineralisation were up-regulated (e.g. carbonic anhydrase 2, alkaline phosphatase), a small group was down-regulated (e.g. perlucin).

Conclusions

This study highlighted the complex molecular response of the bivalve S. glomerata to expected near-future ocean acidification levels. While there are indications that the oyster attempted to adapt to the stressor, gauged by immune system priming and the increase in protein synthesis, some processes such cilia function appear to have been negatively affected by the elevated levels of CO₂.

Background

Anthropogenic activities such as deforestation and burning of fossil fuels have led to a 36 % increase in atmospheric carbon dioxide (CO₂) over the last few hundred years. This is expected to further increase to between 540 to 970 ppm by 2100. Of the anthropogenic CO₂ produced, about one third has been taken up by oceans, leading to a 0.1 unit decrease in the ocean’s surface pH, with predictions of further pH reductions of 0.14–0.35 units by 2100. Furthermore, CO₂ uptake affects carbonate chemistry through changes to the saturation state of, for instance, aragonite and calcite [1–8]. During CO₂ uptake, water (H₂O) interacts with CO₂ to produce carbonic acid (H₂CO₃) in the first instance, which then dissociates to form bicarbonate (HCO₃ ^-) and carbonate ions (CO₃ ^2-) that, in the presence of calcium ions (Ca²⁺), eventually leads to the formation of calcium carbonate (CaCO₃) [5, 6]. Calcium carbonate is an important component of the bivalve shell, as well as coral reefs and other marine calcifiers, with its formation and dissolution strongly dependent on the carbonate saturation state of the water column [4, 6]. Ocean acidification, the term that encompasses the effects of CO₂ uptake by the ocean, can affect shell formation, shell growth and thickness with potential flow-on effects of reduced protection from predators and suboptimal environments [3, 6, 9]. In addition, ocean acidification can potentially affect the metabolic rate of benthic invertebrates, protein synthesis and ion exchange [1, 3].

A number of studies examining the effects of ocean acidification on molluscs have been undertaken, with responses varying between species and changing with ontogeny [9]. For instance, within genera comparison of larval shell growth and calcification in Crassostrea virginica and Crassostrea ariakensis in response to elevated pCO₂, show that while C. ariakensis larvae are not affected, C. virginica larvae displayed a slower shell growth and decreased calcification [7]. C. virginica larvae, as well as larvae of the clam Mercenario mercenaria and Argopecten irradians, also showed delayed metamorphosis and decreased larval size at elevated CO₂. In addition, mortality rates of M. mercenaria and A. irradians larvae were higher at elevated CO₂ than ambient CO₂ [10]. A study on the mussel Mytilus galloprovincialis observed no effect on fertilisation in response to elevated levels of CO₂, but fertilised eggs exposed to elevated CO₂ showed delayed development into D-veliger larvae compared with the control fertilised eggs. Moreover, the majority of the fertilised eggs that eventually developed into D-veliger larvae showed a range of morphological abnormalities (e.g. indentation of shell margin) and had also smaller shells [11].

The Sydney rock oyster, Saccostrea glomerata, is native to the east coast of Australia, where it is both ecologically and economically important [12]. Due to its importance, several studies have also been carried out to assess the potential effects of ocean acidification on this species. For instance, elevated pCO₂ has been shown to decrease fertilisation rate and reduce the percentage of S. glomerata gametes that develop to D-veliger stage. As observed in M. galloprovincialis, S. glomerata D-veliger larvae that had been exposed to elevated pCO₂ were smaller and had higher percentages of morphological abnormalities at 24 h [13]. While much of the ocean acidification research has concentrated on the early life-stages of molluscs, studies indicate that not only larvae, but also juvenile and adult molluscs can be affected by ocean acidification. Observed effects were decreased shell strength, reduced growth (lower dry shell and soft-tissue mass) and increased metabolic rate in response to elevated levels of CO₂ [14, 15].

Estuaries, the habitat of many molluscs such as oysters, are believed to be more vulnerable to increases in pCO₂ and its effects than the open ocean as they often already show high levels of CO₂ and are much shallower [7, 10, 16], which may affect the rate in which the local pH changes. Considering the potentially severe and varied impacts ocean acidification can have on all molluscan life-stages, it is essential to gain a better understanding of the mechanisms underlying the observed effects of ocean acidification on estuarine bivalves. To address this question, we carried out a pCO₂ exposure study using adult S. glomerata. In this study, S. glomerata, acclimated to elevated temperature (28 °C) was exposed to elevated temperature and elevated pCO₂ (1000 ppm) and the molecular response of the oysters to these challenges analysed with RNA-Seq.

Methods

Stress exposure and sample collection

CO₂ and temperature

Adult, wild S. glomerata, collected from Cromarty Bay (seawater temperature of 20.5 ± 0.5 °C and 423.36 ± 17.54 μatm CO₂), Port Stephens (NSW, Australia) were slowly brought up to 28 °C (1 °C change in water temperature/day), then acclimated to the elevated temperature (28 °C) for 4 days. After acclimation, the oysters were exposed to 1) 28 °C and ambient (385 ppm or 386.9 μatm) pCO₂ (control) or 2) 28 °C and elevated (1000 ppm or 1108.2 μatm) pCO₂ in three replicate 750 L header tanks per treatment for 4 weeks. This timeframe was chosen as oysters were expected to be affected by the experimental condition within 4 weeks. Six oysters (n = 2 per replicate) were sampled randomly from each treatment at the end of the experiment.

Additional samples for reference transcriptome preparation

In addition, adult, wild Saccostrea glomerata from Cromarty Bay, Port Stephens (NSW, Australia), were exposed to a) salinity and temperature (n = 24 oysters) or b) polycyclic aromatic hydrocarbons (PAH) (n = 12 oysters). Full details of experimental set-up, oyster husbandry and exposure are as described in Additional file 1. Only the results from the CO₂ and temperature exposure are discussed in this paper; however, sequencing reads from the samples of all treatments (n = 48 oysters total) were used to build the reference transcriptome.

Sample preparation and sequencing

Tissues used for RNA-Seq analysis were gill, mantle, adductor muscle, gonad and digestive (PAH and CO₂ experiment) or gill, mantle, adductor muscle and digestive (salinity). Total RNA was extracted from 25 mg (24 mg for salinity samples) of pooled tissue (5 mg per tissue for CO₂ and PAH samples, 6 mg per tissue for salinity samples) from each of the oysters according to the manufacturer’s guidelines, using the Direct-zol RNA MiniPrep kit (Zymo Research Corporation, USA) including the DNase I digestion step. Quality and quantity of the extracted RNA were tested by gel electrophoresis, bioanalyzer (2100 Bioanalyzer, Agilent Technologies, USA) using the RNA 6000 Nano Chip kit (Agilent Technologies) and with the Quantus^TM fluorometer (Promega, Australia). ERCC (External RNA Control Consortium) RNA spike-in control mixes (Ambion, Australia) were diluted 1:100, then 2 μl of the diluted mix1 and mix2 added to 1 μg of total RNA of control and treated samples, respectively. RNA-Seq libraries were then prepared from each of these samples, using the TruSeq RNA sample prep kit-v2 (Illumina, Australia) according to the kit’s preparation guidelines. Quality and quantity of the 48 cDNA libraries were tested with the 2100 Bioanalyzer, using the High Sensitivity DNA chip kit (Agilent Technologies) and with the Quantus^TM fluorometer. Libraries were sent to AGRF (Australia) for 100 bp paired-end sequencing, where the 48 libraries were randomly distributed across four lanes and run on an Illumina HiSeq 2000 sequencer.

Read processing and de novo reference transcriptome assembly

Paired-end read quality was examined pre- and post-cleanup with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads were trimmed for quality with Trimmomatic [17], and adapter sequences and reads with a size below 40 bp after cleanup removed. In addition, reads were assessed for artefacts like E. coli, phiX, human sequences, as well as for the level of ribosomal sequences. Processed reads, along with S. glomerata processed non-strand specific and non-normalised reads from a previous Illumina HiSeq 2000 sequencing run (Additional file 1) were assembled into transcripts, using Trinity [18] with the in silico read normalisation option. Assembled transcripts were aligned to the ERCC and phiX reference sequences, using BLAT [19] and transcripts mapping to either references removed from the S. glomerata reference transcriptome assembly with an in-house script. The removed ERCC transcripts functioned as indicators of how successful the assembly of the reference transcriptome was and were closely examined in the CLC workbench (CLC Bio, USA, version 7.5). In addition, redundancy was removed with cd-hit-est [20], using a 99 % cut-off, and completeness of the assembly was examined with CEGMA (Core Eukaryotic Genes Mapping Approach) [21]. In order to eliminate any potential bias of the analysis, ERCC reference sequences were added to the clean S. glomerata reference transcriptome that was then used as a reference for the differential gene expression analysis.

Differential gene expression analysis of CO₂ samples

Transcripts significantly differentially expressed between control and elevated CO₂ samples were determined with RSEM [22] and EBSeq [23]. Post-processed reads of the 12 CO₂ samples were aligned to the S. glomerata reference transcriptome with Bowtie [24] (Additional file 2: Table S1), using a RSEM internal script. Read counts were estimated with RSEM, after which read count data for transcripts that had received a zero estimated count for all 12 samples was removed from the analysis. Transcripts with mean counts of >0 across all samples were then further analysed with EBSeq in R (version 3.1.1). As a de novo assembled reference transcriptome was used for differential gene expression analysis where the gene-isoform relationship was not known, mapping ambiguity clusters were produced with the RSEM script rsem-generate-ngvector. Median normalisation of the estimated count data was carried out in EBSeq, after which EBTest was run for 12 iterations until convergence had been reached. A false discovery rate (FDR) threshold of 0.05 was applied to determine isoforms/transcripts significantly differentially expressed between control and elevated CO₂ samples. Fold change values used throughout the text were based on the posterior fold change values.

Functional annotation of differentially expressed transcripts

Blastx similarity searches were carried out on transcripts found to be significantly differentially expressed against the NCBI non-redundant (nr) database (downloaded 08.09.14), using an e-value cut-off of 1e^-5 with a hit number threshold of 25. Mapping and functional annotation of the transcripts was carried out with Blast2GO [25], using standard parameters (hit adjusted to 25). In addition, InterProScan searches were run through Blast2GO with the results merged with the already existing annotations. Where domain/family information was available for transcripts with a sequence description of “---NA---“ or “hypothetical protein/uncharacterised protein”, and this information offered an indication as to the potential identity of the transcript, it was added to the respective transcript. Furthermore, transcripts of interest were mapped to the tissue specific S. glomerata transcripts (Additional file 1), using the CLC Genomics Workbench version 7.5 (CLC Bio, USA) with default parameters (minimum length fraction: 0.7, similarity: 0.8), to determine their putative tissue distribution in S. glomerata.

QPCR analysis

QPCR analysis was carried out to validate the differential transcript expression analysis method used in this study. For this, total RNA was extracted from 25 mg of pooled tissue (5 mg of each, gill, mantle, adductor muscle, gonad and digestive) from the six control (28 °C and 385 ppm pCO₂) and six treated (28 °C and 1000 ppm pCO₂) oysters used to prepare the 12 CO₂ RNA-Seq libraries of this study. In addition, total RNA was also isolated from a wild, non-stressed S. glomerata (using 6 mg each of gill, mantle, adductor muscle, gonad and digestive tissue of the oyster) and used as a reference sample in the qPCR analysis. The Direct-zol RNA MiniPrep kit, including the DNase I digestion step, was used for RNA isolation according to the manufacturer’s guidelines. Quality of the extracted RNA was tested with the 2100 Bioanalyzer, using the RNA 6000 Nano Chip kit, and quantity with the Quantus^TM fluorometer. cDNA synthesis was carried out on the reference sample (500 ng of total RNA), as well as on the control and treated samples (1000 ng of total RNA each), using the QuantiTect® reverse transcription kit (Qiagen, Australia) as described in the kit’s guidelines. Furthermore, negative reverse transcription (-RT) reactions were performed on three out of 12 samples, where RNase free water was used instead of the reverse transcriptase to test for genomic contamination of the samples during qPCR analysis.

Primer design and testing

Of the transcripts found to be significantly differentially expressed between the six control and six CO₂ treated S. glomerata, two transcripts, showing an estimated read count of >100 for the lowest sample were randomly chosen for qPCR analysis. Four transcript specific primer pairs per transcript were determined with Primer3Plus (http://primer3plus.com/cgi-bin/dev/primer3plus.cgi), with a primer size of 19–23 bp, melt temperature of 60–61 °C and low to no self-annealing. Potential primer pairs were mapped back to the S. glomerata reference transcriptome with the CLC Genomics Workbench version 7.5 (CLC Bio, USA) to determine primer pairs that either a) only mapped to the transcript of interest, or b) mapped to the target transcript and a very small number of other transcripts with a minimum of two nucleotide mismatch at the 3′ end. Next,.bam mapping files, produced during the differential transcript expression analysis with Bowtie [24] and RSEM [22] were visualised in the Integrative Genomics Viewer (IGV; https://www.broadinstitute.org/igv/) and used to examine potential primer pairs for nucleotide mismatches in the primer sequence across control and CO₂ treated samples. Primer pairs with a) no nucleotide mismatches across either of the primer sequences in all 12 samples, or b) a maximum of two non-crucial (no 3′) nucleotide mismatches across either of the primer sequences in all 12 samples were synthesised by Sigma Aldrich (Australia) for qPCR analysis (Table 1).

Table 1 Primer pairs for qPCR analysis. This table only shows the highest efficiency primer pairs, as well as the primers used to produce the purified PCR products for the standard curves

A temperature gradient PCR was performed to determine the best primer annealing temperature for each of the eight primer pairs. The primers chosen for use in this study are presented in Table 1. Triplicate PCR reactions were prepared for each of the primer pairs, using 0.6 μL of reference sample as template, added to 14.4 μL of mastermix, containing 1.5 μL of 10x PCR reaction buffer, 2 mM of MgCl₂, 200 μM of dNTPs, 1 U of Taq (all reagents from Fisher Scientific, Australia), 200 nM each of the respective forward and reverse primer, and 10.74 μL of RNase free water. Cycling conditions were an initial denaturation for 1 min at 95 °C, followed by 35 cycles of denaturation for 30 s at 94 °C, annealing for 30 s at either a) 58.1 °C, b) 60.6 °C or c) 61.9 °C with the triplicates distributed across the three temperatures, then extension for 45 s at 72 °C. Final extension occurred at 72 °C for 2 min, after which gel electrophoresis was carried out on 3 μL of each PCR product.

Fragment preparation for qPCR standard curves

PCR products of a known length were produced for the two transcripts chosen for qPCR analysis. Tektin-2 and tektin-4-like PCR products were obtained by PCR amplification, using primer pairs T2F_2 and T2R_4, and T4F_1 and T4R_2, respectively. Four replicate PCR reactions were set up for both transcripts, using 1 μL of reference sample as template, 0.2 μM each of the respective forward and reverse primer, 12.5 μL of MyTaq mix (Bioline, Australia) and 10.5 μL RNase free water for a total reaction volume of 25 μL. PCR amplification conditions were: a) initial denaturation for 2 min at 95 °C, followed by 35 cycles of b) denaturation for 30 s at 95 °C, c) annealing for 30 s at 60 °C, d) extension for 2 min at 72 °C, with a final extension for 2 min at 72 °C. Gel electrophoresis was used to determine single banding of the PCR products, after which products for tektin-2 and tektin-4-like were purified with the QIAquick PCR purification kit (Qiagen, Australia) according to the kit’s guidelines. Purified products of the two transcripts were assessed with gel electrophoresis and the four replicates of each primer pair pooled to obtain a single clean product per primer pair. The two clean products were analysed on the NanoDrop2000 spectrophotometer (Thermo Fisher Scientific, USA) three times and the mean concentration value of each product used to determine copy numbers of each product with the copy number calculator (http://cels.uri.edu/gsc/cndna.html).

Primer validation and qPCR analysis of CO₂ samples

Primer efficiency and specificity were determined with qPCR (Table 1). A 1 in 10, 10 point serial dilution was prepared from the 10⁸ copy number stock solution of each of the two purified PCR products. The serial dilutions were then used as qPCR templates to determine primer efficiency and specificity for all primer pairs. For this, 200 nM each per forward and reverse primer were added to 1 μL of template, 5 μL of Platinum®SYBR® Green qPCR SuperMix-UDG (Invitrogen, Australia) and 3.6 μL of RNase free water. Reactions were performed in triplicate, including triplicate no template controls (NTCs), using the Rotor-Gene 6000 thermal cycler (Corbett Research, Australia). Cycling conditions were as follows: initial holding step at 50 °C for 2 min, hold at 95 °C for 2 min, then 40 cycles of 95 °C for 15 s, 60 °C for 15 s and 72 °C for 25 s, with the last step set to acquire to Green. Melt curve analysis was performed by increasing the melting temperature by 1 °C increments from 72 °C to 95 °C. The Rotor-Gene 6000 software, version 1.7.87 (Corbett Research, Australia), was used for quantification and melt curve analysis. Reaction efficiency (E) was calculated for each primer pair by the Rotor-Gene 6000 software, using the following equation: E = [10(-1/M)] -1, where M stands for the slope of the curve.

Absolute transcript expression levels of tektin-2 and tektin-4-like were determined with qPCR in the six control and six CO₂ treated S. glomerata samples, with reaction volumes and cycling conditions as described for primer validation. Reactions were carried out in duplicate, including duplicate NTCs and –RTs, as well as one point of the standard curve and a positive control in triplicate. Primer pairs used for the individual transcripts were: a) T2F_1 and T2R_1 (efficiency: 0.9855) and b) T4F_3 and T4R_3 (efficiency: 0.9698).

Statistical analysis of qPCR data

A two-sample t-test assuming unequal variances (Microsoft Excel Software) was used to determine significant differences between the transcript expression levels of control and CO₂ stressed S. glomerata. In order to allow for a direct comparison between RNA-Seq and qPCR results, the significance level for the t-test was set at p < 0.05.

Gene ontology (GO) enrichment analysis

Functional annotation of the reference transcriptome

Potential coding regions were determined with TransDecoder (v2.0.1, http://transdecoder.github.io/) and annotated by performing sequence similarity searches against the protein databases UniProt, Swiss-Prot (both with a BitScore > 100) and KOBAS (e-value of 1e^-5) using DIAMOND (v0.8.5) [26], and against the HMMER/Pfam protein database (v28.0, DomainScore > 20) [27] using HMMER3.1b [28]. Where multiple annotations were available for a single open reading frame (ORF) or multiple ORFs were predicted for a given transcript, the best matching annotation (highest BitScore, lowest e-value) or ORF with the best annotation were chosen, respectively.

GO enrichment analysis

Using Trinotate (v3.0.1, https://trinotate.github.io/) scripts, the annotation and differential expression files were loaded into a SQLite database (v3.9.2, http://www.sqlite.org/) and a GO enrichment analysis carried out using the ‘goseq’ R package [29]. Taking the length of each transcript into account, GOseq determined GO and COG /eggNOG (Clusters of Orthologous Groups/evolutionary genealogy of genes: Non-supervised Orthologous Groups) terms that are over-represented among the transcripts differentially expressed between control S. glomerata and oysters exposed to elevated CO₂.

Results and discussion

S. glomerata reference transcriptome

For this study, a S. glomerata reference transcriptome was produced and then used to analyse the molecular response of adult, wild S. glomerata exposed to elevated pCO₂ (1000 ppm) and temperature (28 °C), when compared to control S. glomerata that were exposed to elevated temperature and ambient pCO₂ (385 ppm). A pCO₂ concentration of 1000 ppm was chosen based on predictions for the year 2100 [8]. In order to obtain a comprehensive reference transcriptome, tissue (gill, mantle, adductor muscle, gonad and digestive) of S. glomerata exposed to CO₂ and temperature, salinity and temperature, and PAH was extracted and non-strand specific and non-normalised libraries (n = 48) for each individual oyster prepared. Libraries were sequenced with Illumina, resulting in a total of 818,834,356 paired-end reads with a GC content of 42–44 %. A similar GC content has also been found in other molluscs [30, 31]. Raw reads were processed and the resulting processed reads (98.7 %), along with the processed non-strand specific and non-normalised reads of our prior S. glomerata study (Additional file 1), assembled into a reference transcriptome. Assembly statistics before and after redundancy removal are summarised in Table 2. Completeness of the reference transcriptome assembly was assessed with CEGMA [21], as well as by the successful assembly of ERCC and phiX sequences. Even though the CEGMA software was originally developed to assess the completeness of genomes, various studies [32, 33] have also used this software to determine the completeness of transcriptomes. Of the original 92 ERCC reference sequences and one phiX sequence, 89 ERCC’s and one phiX sequence were found in the S. glomerata reference transcriptome, with all sequences close to full length. The N50 value (836 bp) of the non-redundant reference transcriptome is comparable to the N50’s observed for other molluscan de novo transcriptomes [34–36]. Based on these results, the assembly was considered to be of suitable quality for the differential gene expression analysis.

Table 2 Assembly statistics. Statistics of “total transcripts” includes assembled ERCC and phiX transcripts. All other statistics are excluding ERCC and phiX transcripts

Differential transcript expression analysis of CO₂ samples

EBSeq, an R based program was used to determine S. glomerata transcripts differentially expressed between control oysters (exposed to pCO₂ of 385 ppm and 28 °C) and oysters exposed to elevated pCO₂ (1000 ppm and 28 °C). EBSeq uses an empirical Bayesian approach and was chosen in this study as it takes the estimation uncertainty inherent in isoform expression analysis into consideration [23]. Graphical results of the standard diagnostics on the differential transcript expression analysis with EBSeq are presented in Additional file 3: Figure S1, Additional file 4: Figure S2, Additional file 5: Figure S3 and Additional file 6: Figure S4. A total of 1626 S. glomerata transcripts were found to be differentially expressed (DE) between control and elevated pCO₂ S. glomerata, using a false discovery rate (FDR) threshold of 0.05. Functional annotation of the DE transcripts with Blast2GO against NCBI’s non-redundant database with an e-value cut-off of 1e^-5 resulted in the annotation of 75.2 % of the DE transcripts. Annotation information for 73.9 % of DE transcripts could be obtained when transcripts were searched against the InterProScan database through Blast2GO. Of the functionally annotated DE transcripts, GO-terms associated with cellular and metabolic processes, cell and membrane, and catalytic activity and binding contained the most DE transcripts (Fig. 1a, b and c). While the DE transcripts present only a small section of the S. glomerata reference transcriptome, this pattern of GO-terms is comparable to the most common GO-terms found in the transcriptomes of other molluscs [37, 38]. In addition to the main GO-terms, terms associated with response to stimulus, immune system process, biological regulation, signalling and transporter activity were also observed for the DE transcripts (Fig. 1a, b and c), indicating that the DE transcripts are potentially involved in a wide range of processes and functions.

Fig. 1

GO analysis of S. glomerata DE transcripts. GO-terms were determined with Blast2GO, using default parameters. Level 2 GO-terms are associated with (a) biological process, (b) cellular component and (c) molecular function are depicted

GO enrichment analysis

GO enrichment analysis identified a range of GO and COG/eggNOG terms that were over-represented in either control or CO₂ stressed S. glomerata (Table 3). Of the GO terms found to be over-represented in the control group, ‘oxidation-reduction process’ (biological process) and ‘oxidoreductase activity’ (molecular function) were the most enriched GO term in their respective category. Other enriched GO terms were ‘hydrolase activity’ (8) and ‘carbohydrate metabolic process’ (8). In comparison, ‘calcium ion binding’ (molecular function) and ‘regulation of apoptotic process’ (biological process) were the most enriched GO terms found to be over-represented in CO₂ stressed S. glomerata. Apoptosis is an important component of the innate immune system of molluscs [39] and an over-representation of ‘regulation of apoptotic process’, along with an over-representation of ‘alternative oxidase activity’ in CO₂ stressed oysters could suggest that potentially protective mechanisms were induced in response to the stressor. Analysis of the enriched COG/eggNOG terms indicated an over-representation of, for example, ‘collectin sub-family member 12 (1), ‘glutathione S-transferase (2) or ‘RNA binding motif protein’ (2) in control S. glomerata, whereas ‘alternative oxidase’ (2), ‘solute carrier family 17’ (3) or ‘alkaline phosphatase’ (2) were over-represented in CO₂ stressed oysters.

Table 3 GO terms significantly over-represented among the transcripts differentially expressed between control and CO₂ stressed oysters. This table includes GO (GO_BP = Biological process, GO_MF = Molecular function), COG (Clusters of Orthologous Groups) and eggNOG (evolutionary genealogy of genes: Non-supervised Orthologous Groups) terms that were found to be significantly over-represented among the transcripts differentially expressed in S. glomerata exposed to elevated levels of CO₂ when compared to control oysters

Immunity

Closer examination of the 1626 S. glomerata transcripts differentially expressed between control and CO₂ stressed oysters, showed multiple transcripts, potentially involved in innate immunity. Innate immune responses are triggered by the recognition of pathogen-associated molecular patterns (e.g. lipopolysaccharides) or damage-associated molecular patterns (DAMPs) through pattern recognition receptors (PRRs) [40, 41]. DAMPs are molecules such as heat shock proteins, RNA, DNA, galectins, defensins and annexins, that are released from stressed (e.g. hypoxia), injured or necrotic cells [42, 43]. Many of these DAMPs have been found to be differentially expressed in S. glomerata exposed to elevated CO₂ (Additional file 7: Table S2). For instance, S. glomerata DE transcripts of the heat shock Hsp20 family were up-regulated in response to elevated CO₂, while six out of eight Hsp70 transcripts were 4-fold and higher down-regulated. Aside from heat shock proteins, two out of four annexins were also up-regulated (2-fold and higher) in CO₂ challenged oysters. This correlates with the results of the enrichment analysis, where annexin A7 was shown to be over-represented in S. glomerata exposed to elevated CO₂ (Table 3). In addition, the antimicrobial peptide mantle defensin was 4-fold and higher up-regulated and galectins, which also function as a PRR, had one transcript 4-fold and higher up and one transcript 4-fold and higher down-regulated in challenged oysters. Other PRRs found to be differentially expressed in S. glomerata were gram-negative bacteria binding proteins (GNBPs), scavenger receptors (SRs), fibrinogen-related proteins (tenascins, fibrinogen c domain-containing proteins and fibroleukins), c-type lectins, collectins, peptidoglycan recognition proteins (PGRPs), c-type mannose receptors/macrophage mannose receptors and C1q domain containing proteins (Additional file 7: Table S2). Overall, slightly more PRRs (53.8 %) were up-regulated in response to elevated CO₂ than down-regulated, with, for instance, GNBP transcripts 2-fold and higher down-regulated and four out of six collectins 4-fold and higher up-regulated. This suggests that prolonged CO₂ exposure induced the innate immune defence system of S. glomerata, likely due to the action of the mostly up-regulated DAMPs that were also detected in this study. In addition, SRs which were found to be 4-fold and higher up-regulated in response to CO₂, are also linked to phagocytosis [44]. Aside from SRs, other molecules associated with phagocytosis, such as antimicrobial peptides (bactericidal permeability increasing protein and mantle defensin) and lysozyme were also observed to be 4-fold and higher up-regulated in S. glomerata exposed to elevated CO₂ (Additional file 7: Table S2). Phagocytosis is a mechanism of the innate immune system that recognises and removes dead cells, foreign bodies (e.g. bacteria) and environmental debris [44, 45], and its induction in CO₂ stressed S. glomerata might be a preventive mechanism or a response to damage at the cellular level. A similar expression pattern has also been observed in invertebrates exposed to different types of stress. For instance, the impact of injury on the immune responses of Hydra was analysed and showed an up-regulation of antimicrobial peptides (hydramacin and arminin), a lectin (L-rhamnose binding lectin CSL3) that functions as a PRR and small heat shock proteins in response to injury [46]. The authors suggested that the small heat shock proteins might act as cytoprotectors against ROS (reactive oxygen species) and injury stress, while the Hydra antimicrobial peptides could have additional roles in functions such as regeneration [46]. Exposure of C. virginica to elevated levels of CO₂ (800 and 2000 μatm) for 4 weeks in turn showed decreased Hsp70 mRNA levels and an increase in haemocyte lysozyme activity in response to the stressor [47]. In contrast, Hsp70 gene expression of the coral Desmophyllum dianthus to elevated pCO₂ (997 μatm) for 8 months was strongly up-regulated, along with the gene expression of mannose-binding C-type lectin [48]. Allograft inflammatory factor 1-like, another S. glomerata DE transcript that was found to be 2-fold and higher up-regulated in CO₂ stressed oysters (Additional file 7: Table S2) has been shown to be induced by bacterial challenge, tissue injury and shell damage in the pearl oyster Pinctada martensii, suggesting a role in innate immunity [49]. These results indicate that stressors such as pCO₂ and tissue injury appear to result in a similar induction of the affected animal’s immune system, as has been seen in the S. glomerata exposed to elevated pCO₂ of this study. This might be a pre-emptive strategy to protect the oyster from potential infection during stress exposure or a tissue protective measure in response to potential damage at the cellular level due to 4 weeks of continuous pCO₂ challenge. While no apparent tissue damage or infection was observed throughout the experiment and during sample collection, damage at the cellular level could still have occurred, as has been implied by the increase in DAMPs in CO₂ challenged S. glomerata. Furthermore, considering the presence of potentially opportunistic bacteria in the oyster’s natural habitat [50, 51], non-specific priming of the innate immune system might protect stressed S. glomerata from such invading opportunistic bacteria. That stress can affect the composition of the bacterial community in oysters has already been shown in a study in S. glomerata, where individuals infected with the protozoan paramyxean parasite, Marteilia sydneyi, showed a different bacterial community in their digestive gland than non-infected S. glomerata [51]. Interestingly, one transcript putatively coding for a macrophage-expressed gene 1 (MPEG1) protein and multiple laccase transcripts were found to be 2-fold and higher down-regulated in the pCO₂ stressed S. glomerata (Additional file 7: Table S2). MPEG1 of the disk abalone Haliotis discus discus showed antibacterial activity against Gram-positive and –negative bacteria and was up-regulated in response to bacteria and viral hemorrhagic septicaemia virus [52]. Similarly, a laccase of the sponge Suberites domuncula was up-regulated in response to bacterial lipopolysaccharide and also showed antibacterial activity when the laccase mediator ABTS [2,2′-azino-bis(3-ethylbenzothiazoline-6-sulfonic acid)] was present [53]. Considering the function of MPEG1 and laccase in these invertebrates and their down-regulation in the S. glomerata exposed to elevated pCO₂, it is possible that a) both are only induced as a direct response to specific invading pathogens, or b) protection from potential infection might not be the primary reason for the up-regulation of the other immune factors (e.g. antimicrobial peptides and PRRs) in S. glomerata.

Other DE transcripts observed to be differentially expressed in CO₂ challenged S. glomerata are putatively involved in apoptosis. Caspase, which has a role in cell death signalling [54], was 2-fold and higher down-regulated in S. glomerata exposed to elevated pCO₂, while anti-apoptotic transcripts such as Bcl-2, Fas apoptotic inhibitory molecule (FAIM) and some transcripts of the inhibitor of apoptosis protein (IAP) family were up-regulated in response to elevated pCO₂ (Additional file 7: Table S2). FAIM was shown to protect a variety of cell types such as hepatocytes from cell death [55], and Bcl-2 suppressed irradiation-induced apoptosis in transgenic zebrafish [56]. Similar to FAIM, IAPs can protect different cell types (e.g. neurons, macrophages) from stress induced apoptosis [57]. In our study, two out of five IAP transcripts were 4-fold and higher up-regulated in CO₂ stressed S. glomerata, indicating that a) some apoptotic activity was still occurring or b) as different types of IAPs can have a range of functions in different cells (reviewed in Dubrez-Daloz et al. [57]), it is possible that the different expression pattern seen is linked to separate anti-apoptotic functions of the individual IAPs. A discordant expression pattern of IAPs was also observed in the clam Ruditapes philippinarum exposed to ibuprofen for seven days [58]. Aside from anti-apoptotic transcripts, the pro-apoptotic interferon alpha-inducible protein 27 (IFI27) and programmed cell death protein 7, along with TNF (tumor necrosis factor) ligand superfamily member and TNF receptor superfamily member proteins were also differentially expressed in S. glomerata exposed to elevated pCO₂ (Additional file 7: Table S2). IFI27 that was 4-fold and higher up-regulated in CO₂ stressed S. glomerata, has pro-apoptotic functions in vertebrates, which might be blocked by the simultaneous expression of Bcl-2 [59]. While sparse research is available on IFI27, it could act similarly in oysters, suggesting that the concomitantly up-regulated Bcl-2 could function as an IFI27 mediator in S. glomerata. Another pro-apoptotic S. glomerata DE transcript was programmed cell death protein 7 (less than 4-fold down-regulated), which was shown to cause increased cell apoptosis in mice when over-expressed [60]. TNF and TNF receptor superfamily members that were also observed in the S. glomerata DE transcripts (both 4-fold and higher up-regulated), have diverse functions, such as roles in apoptosis and acute inflammation [61]. TNF-mediated cell death signalling is complex, activating different components through varied pathways and resulting either in cell death (dependent on ROS and JNK [Jun NH₂-terminal kinase] signalling) or cell protection (activation of nuclear factor kB) [54]. Similar to our study, a TNF receptor was found to be up-regulated in regenerating tissue of Hydra after injury [46]. While CO₂ exposed S. glomerata had no obvious injuries, it is possible that the increased TNF and TNF receptor expression in the oysters also had a cell protective function. Overall, it appears that apoptosis was decreased in S. glomerata subjected to 1000 ppm pCO₂ for 4 weeks, potentially to protect its cells from the CO₂ apoptotic stimuli. Comparable to our study, brine shrimps (Artemia sinica) exposed to elevated pCO₂ for 14 days up-regulated the anti-apoptotic factor Apoptosis inhibitor 5, with the authors suggesting that this might be a mechanism to deal with the toxic effect of this stressor [62]. Another study, assessing molecular differences between drought intolerant (Pyganodon grandis) and tolerant (Uniomerus tetralasmus) freshwater mussels showed overall an induction of apoptosis in intolerant mussels and an inhibition of apoptosis in drought tolerant mussels, suggesting that apoptosis inhibition might be an important mechanisms in drought tolerance of U. tetralasmus [63]. In summary, priming of the innate immune system and suppressing cell apoptosis appear to be two mechanisms by which S. glomerata exposed to elevated pCO₂ cope with potentially detrimental effects of ocean acidification.

Respiration and antioxidant defence

S. glomerata transcripts potentially involved in antioxidant defence, ROS production and respiration were also observed among the DE transcripts (Additional file 7: Table S2). Dual oxidase, which is linked to phagocytosis and the production of ROS [64], was found to be 4-fold and higher down-regulated in the elevated treatment when compared to control oysters (Additional file 7: Table S2), suggesting that CO₂ stressed oysters aim to limit their ROS production to protect themselves from oxidative stress. This differs from the expression pattern seen for S. glomerata transcripts which are part of the mitochondrial respiratory chain, where NADH dehydrogenase [ubiquinone] iron-sulfur protein 2 (complex I protein), cytochrome b-c1 complex subunit 7 (complex III protein) and alternative oxidase were 4-fold and higher up-regulated and cytochrome c oxidase subunit 5A mitochondrial (complex IV protein) 4-fold and higher down-regulated in the elevated treatment when compared to control oysters (Additional file 7: Table S2). Both, complex I and complex III are considered to be the major sources of ROS production in the mitochondrial respiratory chain [65, 66], suggesting that up-regulation of mitochondrial respiration in CO₂ challenged S. glomerata would lead to an increase in ROS production, elevating the risk of oxidative damage to the tissues of the oysters. However, this appears to be counteracted by two measures in the S. glomerata exposed to the elevated CO₂ treatment. Firstly, alternative oxidase (Additional file 7: Table S2), located in the inner mitochondrial membrane, allows the animal to circumvent complex III by transferring electrons from coenzyme Q to oxygen and therefore limits the amount of ROS produced during cellular respiration [66, 67]. Secondly, ROS originating from complex I are thought to be removed by mitochondrial matrix antioxidants such as glutathione peroxidase and catalase [65]. In S. glomerata exposed to elevated pCO₂, both catalase (4-fold and higher) and glutathione peroxidase (2-fold) were up-regulated (Additional file 7: Table S2) when compared to control oysters, potentially protecting the animal from ROS produced by complex I. While transcripts coding for the antioxidants peroxiredoxin and glutatredoxin were 4-fold and higher down-regulated in CO₂ challenged S. glomerata, nearly half (two out of five) of the transcripts coding for glutathione S-transferase were 2-fold and higher up-regulated in S. glomerata exposed to elevated pCO₂ for 4 weeks (Additional file 7: Table S2). Of the antioxidants found in the DE transcripts of S. glomerata, catalase (CAT), glutathione peroxidase (GPX), peroxiredoxin and glutaredoxin are able to remove the ROS hydrogen peroxide (H₂O₂), with GPX also functioning in the removal of organic hydroperoxides (e.g. fatty acid and phospholipid hydroperoxides) [68, 69]. Glutathione S-transferases, on the other hand, are active against secondary metabolites (e.g. epoxides, hydroperoxides and unsaturated aldehydes), thus offering protection from the effects of oxidative stress [69, 70]. Ferritin, which is involved in hydroxyl radical scavenging [71], was also found to be 2 to 4-fold up-regulated in S. glomerata exposed to elevated pCO₂ (Additional file 7: Table S2). Up-regulation of CAT, GPX, ferritin and glutathione-S-transferase transcripts, as well as transcripts involved in cellular respiration in the CO₂ challenged S. glomerata suggest that the elevated CO₂ treatment led to oxidative stress, which S. glomerata attempted to counteract by up-regulating a range of antioxidants that act on a variety of substrates. Furthermore, S. glomerata increased the expression of alternative oxidase to potentially reduce production of ROS from cellular respiration (mainly complex III), and decreased expression of dual oxidase whose main function is the production of ROS. A similar pattern of expression has been observed in Crassostrea gigas exposed to hypoxia and in C. virginica exposed to elevated pCO₂ [67, 71, 72]. Alternative oxidase and pyruvate kinase, another S. glomerata transcript found to be 4-fold and higher down-regulated in pCO₂ exposed S. glomerata (Additional file 7: Table S2) of our study, were measured in C. gigas removed from water for 3 h, then re-immersed into either normoxic or hypoxic water [71]. They showed an increase in alternative oxidase mRNA levels in C. gigas re-immersed into normoxic water, and a higher level of pyruvate kinase mRNA levels in oysters under normoxic conditions when compared to oysters in hypoxic water [71]. The authors suggested that the increase in alternative oxidase was linked to increased oxygen consumption and would protect the oyster from the resulting ROS production [71]. Similarly, alternative oxidase mRNA levels in the gills and digestive gland of C. gigas exposed to 12 h and 24 h of hypoxia were shown to be up-regulated when compared to normoxic conditions, suggesting that hypoxia could cause changes to the respiratory function of mitochondria in C. gigas [67]. Furthermore, comparable to the S. glomerata DE transcripts, a varied antioxidant defence response was observed in C. virginica (one peroxiredoxin protein down-regulated, three peroxiredoxins up-regulated) in response to 2 weeks of elevated CO₂ exposure [72], and the spider crab Hyas araneus in response to different concentrations (390 μatm, 1120 μatm and 1960 μatm) of CO₂ for 10 weeks, where peroxiredoxin and ascorbate peroxidase were up-regulated and thioredoxin and thioredoxin peroxidase down-regulated in crabs exposed to 1120 μatm of CO₂ [73]. Moreover, GPX was up-regulated in H. araneus exposed to 1960 μatm of CO₂, while superoxide dismutase and thioredoxin were down-regulated under the same CO₂ conditions [73]. These results indicate that stressors like hypoxia and pCO₂ can affect a range of genes that function in ROS protection, with the specific antioxidant response potentially influenced by the concentration of CO₂ that the invertebrates are exposed to. This has also been observed in the S. glomerata of our study, showing that S. glomerata exposed to 1000 ppm of pCO₂ appear to employ a range of mechanisms (e.g. antioxidants, alternative oxidase) to cope with oxidative stress, potentially caused by exposure to elevated CO₂.

Biomineralisation and cytoskeleton

Bivalve shells consist of an outer layer (periostracum) made up of conchiolin, a prismatic layer (ostracum) and a lamellar layer (hypostracum) which is found closest to the bivalve body. The layers below the periostracum consist of calcium carbonate crystals deposited in an organic matrix (proteins, glycoproteins, polysaccharides and lipids), with the lamellar layer containing aragonite and calcite [74–78]. Several studies have been carried out in molluscs such as the mussel Mytilus edulis, the pearl oysters Pinctada margaritifera and Pinctada maxima, the clam Laternula elliptica and the sea snail Patella vulgata to determine genes potentially associated with shell formation [78–82]. Some of the many genes suggested to be involved in biomineralisation are tyrosinase, chitinase, chitin synthase, calponin, carbonic anhydrase, perlucin, nacrein-like, silk-like, perlustrin, lustrin, follistatin, sarcoplasmic calcium binding and calmodulin [79–81]. In our study, three S. glomerata carbonic anhydrase transcripts were found to be differentially expressed, with two carbonic anhydrase 2 transcripts 4-fold and higher up-regulated and carbonic anhydrase 15-like 4-fold and higher down-regulated (Fig. 2). Carbonic anhydrase, which catalyses the reaction of CO₂ to bicarbonate (HCO₃ ^-) [81], has been examined in different invertebrates and showed varying responses to elevated levels of CO₂. In M. edulis (2 months exposure) carbonic anhydrase was shown to be not significantly affected by elevated pCO₂ [83], while a novel carbonic anhydrase of Hyriopsis cumingii (2 weeks exposure) was significantly down-regulated in elevated pCO₂ [84]. In contrast, carbonic anhydrase was found to be up-regulated in juvenile C. gigas (28 days exposure) [85], in the coral D. dianthus (8 months exposure) [48] and in the crab Hyas araneus (10 weeks exposure) [73] exposed to elevated pCO₂. These results suggest that carbonic anhydrase expression in response to CO₂ stress is not consistent across marine species and exposure times, emphasising the necessity to determine the response of individual marine animals to potential future changes in pCO₂. Furthermore, the different expression pattern observed for the two types of carbonic anhydrase in the S. glomerata of this study indicates that there might also be some variation within a species. As carbonic anhydrases also have roles in processes such as respiration, ion transport or pH homeostasis [84, 86], it is possible that the difference in expression seen in our S. glomerata study was due to different roles of the two types of carbonic anhydrase.

Fig. 2

Heatmap of S. glomerata DE transcripts involved in biomineralisation. Heatmap shows the posterior fold change values of the DE transcripts involved in biomineralisation, with transcripts depicted in red up-regulated and transcripts in green down-regulated in CO₂ treated oysters when compared to control oysters. Full list of 1626 DE transcripts, which include the above biomineralisation transcripts, was determined with Bowtie-RSEM-EBSeq, using a FDR threshold of 0.05

Metabolism of chitin, a biopolymer that has been detected in molluscan shells, involves chitin synthase, chitin deacetylase and chitinase [79, 87, 88], which have all been found to be differentially expressed in the S. glomerata of this study (Fig. 2). Chitin synthase (synthesises chitin) and chitinase (degrades chitin) were 2-fold and higher down-regulated in response to elevated CO₂, whereas chitin deacetylase (modifies chitin) was less than 2-fold up-regulated, suggesting that oysters exposed to 1000 ppm pCO₂ protect the chitin shell component from degradation and concentrate on chitin modification instead of chitin synthesis. Additional work is needed to determine the effect (detrimental or beneficial) of chitin modification on S. glomerata shells. A similar pattern of expression has been shown in M. edulis, where elevated levels of CO₂ exposure resulted in a decrease in chitinase; however, chitin synthase showed no response to high CO₂ levels [83]. DE transcripts putatively coding for follistatin, silk like protein, perlucin, alkaline phosphatase, sarcoplasmic calcium-binding protein (one transcript 2-fold and higher up-, one transcript 4-fold and higher down-regulated), calmodulin and calmodulin-like, which have also been associated with biomineralisation in molluscs [79–81, 89], were also found in S. glomerata exposed to elevated pCO₂ (Fig. 2). While perlucin and calmodulin were 4-fold and higher down-regulated, silk like protein, follistatin, calmodulin-like and alkaline phosphatase were up to 4-fold and higher up-regulated in CO₂ stressed S. glomerata. In addition, the enrichment analysis also showed an over-representation of the COG/eggNOG terms ‘follistatin’ and ‘alkaline phosphatase’ in CO₂ challenged S. glomerata (Table 3), further highlighting the effect CO₂ exposure could potentially have on biomineralisation. Contrary to our results, sarcoplasmic calcium-binding protein and silk like protein expression in M. edulis did not change under elevated pCO₂ [83]. Alkaline phosphatase on the other hand, was also found to be up-regulated in A. sinica exposed to two levels of CO₂ for 7 days, but showed a decrease in expression levels at day 14 [62]. A study in Pinctada fucata examined calmodulin (role in calcium metabolism [90]) and calmodulin-like and their potential functions in biomineralisation and found that both had a role in biomineralisation but were differently distributed in the mantle tissue of the oyster [91]. In general, these results suggest that elevated pCO₂ levels might have affected the shell composition of S. glomerata by causing a down-regulation of some of the biomineralisation genes (e.g. perlucin), and up-regulation of others (e.g. chitin deacetylase, calmodulin-like), which could have consequences in regards to shell strength. This could potentially explain the decrease in shell strength that has been observed in P. fucata exposed to low pH [14], as well as the significant changes in shell strength observed in C. gigas exposed to the same elevated CO₂ conditions used within our study [92]. Also, in accord with the molecular results of our study in S. glomerata, research analysing the shell ultrastructure of the mussel M. edulis exposed to projected near-future ocean acidification levels (550, 750 and 1000 μatm pCO₂) showed changes in shell composition and crystal formation [93]. While shell growth did not appear to be affected by exposure to 1000 μatm pCO₂, a significant reduction in aragonite thickness and a significant increase in calcite thickness was observed under this treatment condition [93]. In addition, new calcite crystals formed in the shell of M. edulis exposed to elevated ocean acidification levels were disorientated and could potentially result in reduced shell strength in the exposed mussels [93].

Other S. glomerata DE transcripts found are putative elements of the cytoskeleton and extracellular matrix (ECM). The eukaryotic cytoskeleton consists of microtubules (α- and β-tubulin), microfilaments (actin) and intermediate filaments (proteins of the keratin family), maintains cell shape and has a role in motility (e.g. movement of organelles) [94, 95]. Among the S. glomerata DE transcripts, one beta tubulin and one actin transcript were less than 2-fold up-regulated, while four actin transcripts were 4-fold and higher down-regulated in S. glomerata exposed to elevated pCO₂ (Additional file 7: Table S2). Similar transcript expression has been observed in H. araneus in response to different levels of elevated pCO₂, where alpha and beta tubulin transcripts, as well as actin transcripts were up-regulated in response to the CO₂ treatment [73]. In addition, a β-actin transcript and different actin proteins of D. dianthus and C. virginica, respectively, were also found to be up-regulated in response to elevated pCO₂ [48, 72]. While the majority of actin transcripts were down-regulated in our study, one tubulin and one actin transcript were up-regulated, indicating that the increase in pCO₂ for 4 weeks might have had some impact on the cytoskeleton, forcing the oysters to compensate by slightly increasing the expression of cytoskeletal transcripts. It has been suggested that oxidative stress could affect the cytoskeleton [73], which, considering the increase in complex I and complex III transcripts along with the up-regulation of antioxidant defence mechanisms appears to be supported by the results of our S. glomerata study. Integrin, a transmembrane receptor that connects the ECM with the actin cytoskeleton [96, 97] was 4-fold and higher up-regulated in CO₂ exposed S. glomerata (Additional file 7: Table S2). As integrins are also involved in transmitting mechanical and chemical information to the cytoskeleton [96], it is possible that the up-regulation of the S. glomerata integrin transcript is a coping mechanism of the oyster to the CO₂ stress to protect the cytoskeleton. In addition to cytoskeletal components, transcripts putatively coding for ECM components were also found to be differentially expressed in S. glomerata exposed to elevated pCO₂ for 4 weeks. The ECM in eukaryotes, which is composed of the macromolecules glycoproteins and fibrous proteins, provides a physical scaffold for cells and has other roles such as transmitting mechanical cues [97–99]. Collagens, which are the most abundant proteins of the ECM [97, 99], and an aggrecan core protein that is also a component of the ECM [99] were differentially expressed in the S. glomerata of our study (Additional file 7: Table S2). While the aggrecan core protein transcript was 4-fold and higher down-regulated in CO₂ challenged S. glomerata, five collagen transcripts were 4-fold and higher up-regulated and eight were 4-fold and higher down-regulated in response to elevated pCO₂. In addition, one ADAMTS (a disintegrin and metalloproteinase with thrombospondin motifs) transcript was found to be less than 2-fold down-regulated and two TIMPs (metalloproteinase inhibitors) were 4-fold and higher up-regulated in S. glomerata exposed to elevated pCO₂ (Additional file 7: Table S2). ADAMTS are a family of proteinases that function in ECM degradation [98], whereas TIMPs are regulators of ADAMTS and have a role in the continuous remodelling of the ECM, which is important in maintaining homeostasis during, for instance, injury [97–100]. The differential expression pattern of S. glomerata ECM transcripts suggests that ECM degradation is inhibited, while a small number of collagen transcripts are up-regulated to potentially counteract any detrimental effects caused by the elevated CO₂ stress exposure. A slightly different response was observed in C. virginica exposed to elevated pCO₂ for 2 weeks, where the cytoskeletal component actin was up-regulated but one collagen protein down-regulated [72], indicating that hypercapna negatively affected the ECM but induced actin to potentially protect the cytoskeleton against ROS in this oyster species. Interestingly, studies on TIMP function in oysters surmised that it could have additional roles. For instance, a TIMP in P. martensii was shown to have a putative role in nacre formation as suppression of TIMP expression resulted in disordered growth of the oyster’s nacre [101]. Similarly, TIMP expression in C. gigas was up-regulated in oysters with damaged shells [102]. Based on these studies, it is possible that the up-regulation of S. glomerata TIMPs of this study could also have been in response to putative damage caused by the prolonged CO₂ exposure of the oysters.

Protein synthesis

Nuclear respiratory factor 1 (NRF-1), which is a nuclear transcription factor that is involved in the transcriptional expression of mitochondrial respiratory chain components and other mitochondrion-related genes [103], has been found to be up-regulated (less than 2-fold) in S. glomerata in response to elevated concentrations of CO₂ (Additional file 7: Table S2). Based on the function of NRF-1 in vertebrates, its up-regulation in S. glomerata could potentially be linked to the up-regulation of mitochondrial respiratory chain components, which was also observed in this study, indicating that maintenance of the respiratory apparatus was of importance to the CO₂ challenged S. glomerata. Comparable to NRF-1, transcripts, putatively coding for ribosomal proteins, were found to be 2-fold and higher up-regulated (seven out of eight transcripts) in S. glomerata after CO₂ exposure (Additional file 7: Table S2). Ribosomal proteins have a role in protein synthesis [104] and similar to our study, were up-regulated in D. dianthus and C. virginica following CO₂ exposure [48, 72], showing that protein synthesis was not only affected by CO₂ in our S. glomerata study but also in other marine organisms. Aside from ribosomal proteins and NRF-1, transcripts belonging to the family of mRNA helicases and RNA-binding proteins were also differentially expressed in CO₂ challenged S. glomerata (Additional file 7: Table S2). Of these transcripts, three are putatively coding for eukaryotic translation initiation factors (eIFs), with eIF4A2, eIF4B and eIF3h less than 2-fold up-regulated in the elevated CO₂ treatment. In eukaryotes, eIF4 and eIF3 are involved in translation initiation, where they unwind secondary structures in the mRNA 5′ untranslated region [105]. Similarly, heterogeneous nuclear ribonucleoproteins (hnRNPs) are RNA-binding proteins with functions in the nucleus and cytoplasm of the cells [106] that were mostly less than 2-fold up-regulated in the S. glomerata of this study. The roles of hnRNPs range from transcription to mRNA transport, splicing, 3′-end processing and mRNA stability [106], showing that hnRNPs are an important group of transcriptional regulation proteins. In P. fucata, a hnRNP was cloned and shown to be expressed in the gonad, gill and viscera of the oyster, with its localisation inside the cell restricted to the nucleus [107]. Increasing the expression of transcripts putatively coding for proteins involved in protein synthesis would allow CO₂ stressed S. glomerata to express molecules that might make them more resilient to the stressor. In addition, a slightly higher percentage of the total 1626 DE transcripts were up-regulated (53.3 %) than down-regulated, which might make it necessary to increase the number of proteins involved in translation and post-translational processing.

Ciliary function

In bivalves, cilia are found on a variety of tissues such as gill, mantle or stomach, with important functions in, for instance, filtration, respiration and pseudofeces expulsion [108–111]. Multiple S. glomerata DE transcripts putatively involved in cilia formation or function (e.g. tektin, CCDC65, CCDC176) were found to be less than 4-fold down-regulated in oysters exposed to elevated pCO₂ for 4 weeks (Fig. 3). These RNA-Seq results were confirmed with qPCR, which showed that both, tektin-2 and tektin-4-like were significantly down-regulated in S. glomerata exposed to elevated pCO₂ for 4 weeks (Fig. 4). Tektins are a family of proteins that have been observed in eukaryotic organisms such as mammals, insects and sea urchins and in a wide range of tissues (e.g. testis, brain), and are vital components of cilia and flagella [112]. Similarly, analysis of the putative tissue distribution of the tektin transcripts found to be differentially expressed in our study showed that they were expressed in the haemolymph, gill, mantle, adductor muscle, digestive system and gonad of S. glomerata (mapping data not shown). Comparable to the results of our study, tektin has also been found to be down-regulated in C. gigas larvae exposed to a pH of about 7.5 for 6 days [113]. In addition, another study in C. gigas showed that two forms of tektin were present in the spermatozoa of the oyster [114]. Correspondingly, CCDC135 (coiled-coil domain-containing protein lobo homolog), a flagellar protein, has been shown to localize in the testis and along the sperm flagellum of the fly Drosophila melanogaster, with some sperm motility defects seen in the single and double mutants [115]. Another coiled-coil domain containing protein (CCDC65) was detected in the human sperm tail and the cilia of airway epithelial cells, with a mutation or silencing of the protein negatively affecting cilia motility [116]. Based on these studies and the importance of the flagellum for sperm movement, down-regulation of tektin, CCDC135 and CCDC65 in the S. glomerata of our study could mean a potential impairment of sperm motility in the affected oysters. While some research has examined sperm motility, contradictory results have been obtained in these studies, with C. gigas sperm exposed to elevated CO₂ for 48 h showing a decrease in motility and Strongylocentrotus nudus sperm (20 min exposure) not affected by elevated CO₂ [117, 118]. Although knowledge regarding the effect of increased pCO₂ on gametes is important, exposure of the adult before it reaches gravid stage could already potentially impact on the gametes before they are released into the water column.

Fig. 3

Heatmap of S. glomerata DE transcripts involved in cilia formation and function. Heatmap shows the posterior fold change values of the DE transcripts involved in cilia formation and function, with transcripts depicted in red up-regulated and transcripts in green down-regulated in CO₂ treated oysters when compared to control oysters. Full list of 1626 DE transcripts, which include the above biomineralisation transcripts, was determined with Bowtie-RSEM-EBSeq, using a FDR threshold of 0.05

Fig. 4

QPCR results of differential transcript expression analysis. Transcript expression levels (mean ± SD) of tektin-2 and tektin-4-like were analysed with qPCR in control (n = 6) and CO₂ stressed (n = 6) S. glomerata, and significant differences examined with a two-sample t-test assuming unequal variances (p < 0.05)

Other S. glomerata DE transcripts putatively involved in cilia function and formation are CCDC176 (coiled-coil domain-containing protein 176), CC2D2A (coiled-coil and c2 domain-containing protein 2A), CEP131 (5-azacytidine-induced protein 1) and bardet-biedl syndrome proteins (Fig. 3). A study of CCDC176 in Xenopus laevis showed an up-regulation of the gene during the formation of motile cilia and an involvement in the alignment and maintenance of cilia orientation [119]. The authors also observed that misalignment of the cilia in X. laevis appeared to have negatively affected the flow generated by the beating cilia [119]. Bardet-biedl syndrome (BBS) genes play a role in intraflagellar transport and ciliary function, with BBS2, BBS4 and BBS6 observed to affect flagella formation of sperm in mice and BBS2 and BBS4-BBS8 loss in zebrafish was associated with maintenance or survival of cilia of Kupffer’s vesicle [120]. CEP131 and CC2D2A, on the other hand, have both been implicated in ciliogenesis in vertebrates [121, 122]. While research on the effects of stress on bivalve ciliary motility appears to be sparse, studies in M. edulis showed a decrease in ciliary activity in response to gamma radiation and loss of cilia in mussels exposed to environmental pollution [123, 124]. Overall, the molecular response of CO₂ stressed S. glomerata of our study suggest that elevated pCO₂ could also negatively affect cilia function in oysters, potentially resulting in impairments to the oyster’s ability to generate flow for filtration and sperm movement.

Conclusions

RNA-Seq examination of the molecular response of adult S. glomerata exposed to expected near-future ocean acidification levels for 4 weeks showed that overall slightly more of the DE transcripts were up-regulated than down-regulated. The expression patterns suggest that S. glomerata protects itself from the stressor by priming its immune system, suppressing cell apoptosis and adjusting antioxidant defence mechanisms to compensate for what appears to be an increase in mitochondrial respiration and its respective leakage of ROS. Aside from this protective mechanism, that is similar to what has been seen in response to injury and pathogens in other marine organisms, transcripts associated with maintenance and repair were also up-regulated in S. glomerata. These transcripts have putative roles in translation and post-translational processing and were likely increased to cover the elevated expression of protective molecules such as antioxidants and other immune related proteins. Contrary to this, prolonged exposure to elevated CO₂ seems to have negatively affected structural proteins (e.g. actin) and proteins putatively involved in cilia and flagella function. Cilia are important structures with many functions, such as facilitating filtering and particle transport into the stomach. The decrease in the expression of transcripts implicated in cilia and flagella function suggests that increased pCO₂ might impair processes in the oyster that rely on their optimal action, such as feeding or sperm motility, which in the long-term could be life-threatening. While some ocean acidification studies have been carried out in regards to its effect on sperm motility, these studies only assessed how a very short exposure to elevated pCO₂ impacts the released sperm, ignoring that projected near-future ocean acidification levels could also potentially affect sperm production and maturation, which could have downstream effects on its motility and success in fertilizing an egg. Expression patterns of transcripts putatively involved in biomineralisation, suggest that continuous exposure of S. glomerata to elevated CO₂ resulted in a change of the shell composition with potential downstream effects on shell strength. In contrast to larvae that quite likely have to extend more energy into the initial formation of the shell, adults only need to maintain and slowly grow their shells. As our study assessed the effects of elevated CO₂ on adult S. glomerata, the increase in biomineralisation transcript expression seen might be sufficient to maintain an adequate level of biomineralisation. In summary, this study detailed the complex molecular response of S. glomerata to projected near-future levels of ocean acidification. However, to fully elucidate the molecular and physiological response of bivalves to future ocean acidification levels, long-term studies need to be carried out that include recovery periods to assess the potential of bivalves to reverse any detrimental effects of ocean acidification.