Skip to main content


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



Human activities have led to a substantial increase in carbon dioxide (CO2) 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 pCO2.


A total of 1626 S. glomerata transcripts were found to be differentially expressed in oysters exposed to elevated pCO2 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 CO2 appear to have resulted in a priming of the immune system and in producing countermeasures to potential oxidative stress. CO2 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).


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 CO2.


Anthropogenic activities such as deforestation and burning of fossil fuels have led to a 36 % increase in atmospheric carbon dioxide (CO2) over the last few hundred years. This is expected to further increase to between 540 to 970 ppm by 2100. Of the anthropogenic CO2 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, CO2 uptake affects carbonate chemistry through changes to the saturation state of, for instance, aragonite and calcite [18]. During CO2 uptake, water (H2O) interacts with CO2 to produce carbonic acid (H2CO3) in the first instance, which then dissociates to form bicarbonate (HCO3 -) and carbonate ions (CO3 2-) that, in the presence of calcium ions (Ca2+), eventually leads to the formation of calcium carbonate (CaCO3) [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 CO2 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 pCO2, 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 CO2. In addition, mortality rates of M. mercenaria and A. irradians larvae were higher at elevated CO2 than ambient CO2 [10]. A study on the mussel Mytilus galloprovincialis observed no effect on fertilisation in response to elevated levels of CO2, but fertilised eggs exposed to elevated CO2 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 pCO2 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 pCO2 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 CO2 [14, 15].

Estuaries, the habitat of many molluscs such as oysters, are believed to be more vulnerable to increases in pCO2 and its effects than the open ocean as they often already show high levels of CO2 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 pCO2 exposure study using adult S. glomerata. In this study, S. glomerata, acclimated to elevated temperature (28 °C) was exposed to elevated temperature and elevated pCO2 (1000 ppm) and the molecular response of the oysters to these challenges analysed with RNA-Seq.


Stress exposure and sample collection

CO2 and temperature

Adult, wild S. glomerata, collected from Cromarty Bay (seawater temperature of 20.5 ± 0.5 °C and 423.36 ± 17.54 μatm CO2), 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) pCO2 (control) or 2) 28 °C and elevated (1000 ppm or 1108.2 μatm) pCO2 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 CO2 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 CO2 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 CO2 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 QuantusTM 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 QuantusTM 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 ( 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 CO2 samples

Transcripts significantly differentially expressed between control and elevated CO2 samples were determined with RSEM [22] and EBSeq [23]. Post-processed reads of the 12 CO2 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 CO2 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 pCO2) and six treated (28 °C and 1000 ppm pCO2) oysters used to prepare the 12 CO2 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 QuantusTM 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 CO2 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 (, 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; and used to examine potential primer pairs for nucleotide mismatches in the primer sequence across control and CO2 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 MgCl2, 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 (

Primer validation and qPCR analysis of CO2 samples

Primer efficiency and specificity were determined with qPCR (Table 1). A 1 in 10, 10 point serial dilution was prepared from the 108 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 CO2 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 CO2 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, 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, scripts, the annotation and differential expression files were loaded into a SQLite database (v3.9.2, 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 CO2.

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 pCO2 (1000 ppm) and temperature (28 °C), when compared to control S. glomerata that were exposed to elevated temperature and ambient pCO2 (385 ppm). A pCO2 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 CO2 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 [3436]. 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 CO2 samples

EBSeq, an R based program was used to determine S. glomerata transcripts differentially expressed between control oysters (exposed to pCO2 of 385 ppm and 28 °C) and oysters exposed to elevated pCO2 (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 pCO2 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 CO2 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 CO2 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 CO2 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 CO2 stressed oysters.

Table 3 GO terms significantly over-represented among the transcripts differentially expressed between control and CO2 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 CO2 when compared to control oysters


Closer examination of the 1626 S. glomerata transcripts differentially expressed between control and CO2 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 CO2 (Additional file 7: Table S2). For instance, S. glomerata DE transcripts of the heat shock Hsp20 family were up-regulated in response to elevated CO2, 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 CO2 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 CO2 (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 CO2 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 CO2 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 CO2, 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 CO2 (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 CO2 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 CO2 (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 pCO2 (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 CO2 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 pCO2 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 pCO2 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 pCO2 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 CO2 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 pCO2 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 pCO2, 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 CO2 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 pCO2, 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 pCO2 (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 CO2 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 pCO2 (Additional file 7: Table S2). IFI27 that was 4-fold and higher up-regulated in CO2 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 NH2-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 CO2 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 pCO2 for 4 weeks, potentially to protect its cells from the CO2 apoptotic stimuli. Comparable to our study, brine shrimps (Artemia sinica) exposed to elevated pCO2 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 pCO2 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 CO2 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 CO2 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 CO2 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 pCO2, 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 CO2 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 pCO2 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 (H2O2), 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 pCO2 (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 CO2 challenged S. glomerata suggest that the elevated CO2 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 pCO2 [67, 71, 72]. Alternative oxidase and pyruvate kinase, another S. glomerata transcript found to be 4-fold and higher down-regulated in pCO2 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 CO2 exposure [72], and the spider crab Hyas araneus in response to different concentrations (390 μatm, 1120 μatm and 1960 μatm) of CO2 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 CO2 [73]. Moreover, GPX was up-regulated in H. araneus exposed to 1960 μatm of CO2, while superoxide dismutase and thioredoxin were down-regulated under the same CO2 conditions [73]. These results indicate that stressors like hypoxia and pCO2 can affect a range of genes that function in ROS protection, with the specific antioxidant response potentially influenced by the concentration of CO2 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 pCO2 appear to employ a range of mechanisms (e.g. antioxidants, alternative oxidase) to cope with oxidative stress, potentially caused by exposure to elevated CO2.

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 [7478]. 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 [7882]. 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 [7981]. 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 CO2 to bicarbonate (HCO3 -) [81], has been examined in different invertebrates and showed varying responses to elevated levels of CO2. In M. edulis (2 months exposure) carbonic anhydrase was shown to be not significantly affected by elevated pCO2 [83], while a novel carbonic anhydrase of Hyriopsis cumingii (2 weeks exposure) was significantly down-regulated in elevated pCO2 [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 pCO2. These results suggest that carbonic anhydrase expression in response to CO2 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 pCO2. 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 CO2 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 CO2, whereas chitin deacetylase (modifies chitin) was less than 2-fold up-regulated, suggesting that oysters exposed to 1000 ppm pCO2 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 CO2 exposure resulted in a decrease in chitinase; however, chitin synthase showed no response to high CO2 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 [7981, 89], were also found in S. glomerata exposed to elevated pCO2 (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 CO2 stressed S. glomerata. In addition, the enrichment analysis also showed an over-representation of the COG/eggNOG terms ‘follistatin’ and ‘alkaline phosphatase’ in CO2 challenged S. glomerata (Table 3), further highlighting the effect CO2 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 pCO2 [83]. Alkaline phosphatase on the other hand, was also found to be up-regulated in A. sinica exposed to two levels of CO2 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 pCO2 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 CO2 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 pCO2) showed changes in shell composition and crystal formation [93]. While shell growth did not appear to be affected by exposure to 1000 μatm pCO2, 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 pCO2 (Additional file 7: Table S2). Similar transcript expression has been observed in H. araneus in response to different levels of elevated pCO2, where alpha and beta tubulin transcripts, as well as actin transcripts were up-regulated in response to the CO2 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 pCO2 [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 pCO2 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 CO2 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 CO2 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 pCO2 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 [9799]. 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 CO2 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 pCO2. 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 pCO2 (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 [97100]. 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 CO2 stress exposure. A slightly different response was observed in C. virginica exposed to elevated pCO2 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 CO2 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 CO2 (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 CO2 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 CO2 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 CO2 exposure [48, 72], showing that protein synthesis was not only affected by CO2 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 CO2 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 CO2 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 CO2 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 [108111]. 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 pCO2 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 pCO2 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 CO2 for 48 h showing a decrease in motility and Strongylocentrotus nudus sperm (20 min exposure) not affected by elevated CO2 [117, 118]. Although knowledge regarding the effect of increased pCO2 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 CO2 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 CO2 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 CO2 stressed S. glomerata of our study suggest that elevated pCO2 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.


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 CO2 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 pCO2 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 pCO2 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 CO2 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 CO2 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.


A. irradians :

Argopecten irradians

A. sinica :

Artemia sinica


2’-azino-bis(3-ethylbenzothiazoline-6-sulfonic acid)


A disintegrin and metalloproteinase with thrombospondin motifs


Australian Genome Research Facility


Bardet-biedl syndrome

C. ariakensis :

Crassostrea ariakensis

C. gigas :

Crassostrea gigas

C. virginica :

Crassostrea virginica




Coiled-coil and c2 domain-containing protein 2A


Coiled-coil domain containing protein


Coiled-coil domain-containing protein lobo homolog


Core Eukaryotic Genes Mapping Approach


5-azacytidine-induced protein 1

D. dianthus :

Desmophyllum dianthus


Damage-associated molecular patterns


Differentially expressed


Extracellular matrix


Eukaryotic translation initiation factors


External RNA Control Consortium


Fas apoptotic inhibitory molecule


False discovery rate


Gram-negative bacteria binding proteins


Glutathione peroxidase

H. araneus :

Hyas araneus


Heterogeneous nuclear ribonucleoproteins


Heat shock protein


Inhibitor of apoptosis protein


Interferon alpha-inducible protein 27


Integrative genomics viewer


Jun NH2-terminal kinase

M. edulis :

Mytilus edulis

M. galloprovincialis :

Mytilus galloprovincialis

M. mercenaria :

Mercenario mercenaria


Macrophage-expressed gene 1


National Center for Biotechnology Information




Nuclear respiratory factor 1


New South Wales


No template controls

P. fucata :

Pinctada fucata


Polycyclic aromatic hydrocarbon


Peptidoglycan recognition proteins


Pattern recognition receptors


Quantitative polymerase chain reaction


Reactive oxygen species


Negative reverse transcriptions

S. glomerata :

Saccostrea glomerata


Sequence read archive


Scavenger receptors


Metalloproteinase inhibitors


Tumor necrosis factor

U. tetralasmus :

Uniomerus tetralasmus

X. laevis :

Xenopus laevis


  1. 1.

    Harley CDG, Hughes AR, Hultgren KM, Miner BG, Sorte CJB, Thornber CS, et al. The impacts of climate change in coastal marine systems. Ecol Lett. 2006;9:228–41.

  2. 2.

    Poloczanska ES, Brown CJ, Sydeman WJ, Kiessling W, Schoeman DS, Moore PJ, et al. Global imprint of climate change on marine life. Nat Clim Chang. 2013;3:919–25.

  3. 3.

    Przeslawski R, Ahyong S, Byrne M, Wörheides G, Hutchings P. Beyond corals and fish: the effects of climate change on noncoral benthic invertebrates of tropical reefs. Glob Chang Biol. 2008;14:2773–95.

  4. 4.

    Zeebe RE, Zachos JC, Caldeira K, Tyrrell T. Carbon emissions and acidification. Science. 2008;321:51–2.

  5. 5.

    Lee S-W, Park S-B, Jeong S-K, Lim K-S, Lee S-H, Trachtenberg MC. On carbon dioxide storage based on biomineralization strategies. Micron. 2010;41:273–82.

  6. 6.

    Doney SC, Fabry VJ, Feely RA, Keypas JA. Ocean acidification: the other CO2 problem. Ann Rev Mar Sci. 2009;1:169–92.

  7. 7.

    Miller AW, Reynolds AC, Sobrino C, Riedel GF. Shellfish face uncertain future in high CO2 world: Influence of acidification on oyster larvae calcification and growth in estuaries. PLoS ONE. 2009;4(5):e5661.

  8. 8.

    IPCC. Climate change 2007: the physical science basis. Cambridge: Cambridge University Press; 2007.

  9. 9.

    Gazeau F, Parker LM, Comeau S, Gattuso J-P, O’Connor WA, Martin S, et al. Impacts of ocean acidification on marine shelled molluscs. Mar Biol. 2013;160:2207–45.

  10. 10.

    Talmage SC, Gobler CJ. The effects of elevated carbon dioxide concentrations on the metamorphosis, size, and survival of larval hard clams (Mercenaria mercenaria), bay scallops (Argopecten irradians), and Eastern oysters (Crassostrea virginica). Limnol Oceanogr. 2009;54(6):2072–80.

  11. 11.

    Kurihara H, Kato S, Ishimatsu A. Effects of increased seawater pCO2 on early development of the oyster Crassostrea gigas. Aquat Biol. 2007;1:91–8.

  12. 12.

    O’Connor W, Dove MC. The changing face of oyster culture in New South Wales, Australia. J Shellfish Res. 2009;28(4):803–11.

  13. 13.

    Parker LM, Ross PM, O’Connor WA. The effect of ocean acidification and temperature on the fertilization and embryonic development of the Sydney rock oyster Saccostrea glomerata (Gould 1850). Glob Chang Biol. 2009;15:2123–36.

  14. 14.

    Welladsen HM, Southgate PC, Heimann K. The effects of exposure to near-future levels of ocean acidification on shell characteristics of Pinctada fucata (Bivalvia: Pteriidae). Molluscan Res. 2010;30(3):125–30.

  15. 15.

    Beniash E, Ivanina A, Lieb NS, Kurochkin I, Sokolova IM. Elevated level of carbon dioxide affects metabolism and shell formation in oysters Crassostrea virginica. Mar Ecol Prog Ser. 2010;419:95–108.

  16. 16.

    Lannig G, Eilers S, Pörtner HO, Sokolova IM, Bock C. Impact of ocean acidification on energy metabolism of oyster, Crassostrea gigas - changes in metabolic pathways and thermal response. Mar Drugs. 2010;8:2318–39.

  17. 17.

    Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

  18. 18.

    Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nature Biotechnolgy. 2011;29:644–52.

  19. 19.

    Kent WJ. BLAT - the BLAST-like alignment tool. Genome Res. 2002;12:656–64.

  20. 20.

    Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22(13):1658–9.

  21. 21.

    Parra G, Bradnam K, Korf I. CEGMA: a pipeline to accurately annotate core genes in eukaryotic genomes. Bioinformatics. 2007;23(9):1061–7.

  22. 22.

    Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.

  23. 23.

    Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BMG, et al. EBSeq: an empirical Bayes hierarchical model for interference in RNA-Seq experiments. Bioinformatics. 2013;29(8):1035–43.

  24. 24.

    Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.

  25. 25.

    Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.

  26. 26.

    Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12(1):59–60.

  27. 27.

    Finn RD, Bateman A, Clements J, Coggill P, Eberhardt RY, Eddy SR, et al. Pfam: the protein families database. Nucleic Acids Res. 2014;42:D222–D30.

  28. 28.

    Eddy SR. Accelerated profile HMM searches. PLoS Comput Biol. 2011;7(10):e1002195.

  29. 29.

    Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.

  30. 30.

    Cantacessi C, Prasopdee S, Sotillo J, Mulvenna J, Tesana S, Loukas A. Coming out of the shell: building the molecular infrastructure for research on parasite-harbouring snails. PLoS Negl Trop Dis. 2013;7(9):e2284.

  31. 31.

    Deng Y, Lei Q, Tian Q, Xie S, Du X, Li J et al. De novo assembly, gene annotation, and simple sequence repeat marker development using Illumina paired-end transcriptome sequences in the pearl oyster Pinctada maxima. Biosci Biotechnol Biochem. 2014;78(10):1685–92.

  32. 32.

    Nakasugi K, Crowhurst RN, Bally J, Wood CC, Hellens RP, Waterhouse PM. De novo transcriptome sequence assembly and analysis of RNA silencing genes of Nicotiana benthamiana. PLoS ONE. 2013;8(3):e59534.

  33. 33.

    Hu H, Bandyopadhyay PK, Olivera BM, Yandell M. Characterization of the Conus bullatus genome and its venom-duct transcriptome. BMC Genomics. 2011;12:60.

  34. 34.

    Wit PD, Palumbi SR. Transcriptome-wide polymorphisms of red abalone (Haliotis rufescens) reveal patterns of gene flow and local adaptation. Mol Ecol. 2013;22:2884–97.

  35. 35.

    Wang S, Hou R, Bao Z, Du H, He Y, Su H, et al. Transcriptome sequencing of Zhikong scallop (Chlamys farreri) and comparative transcriptomic analysis with Yesso scallop (Patinopecten yessoensis). PLoS ONE. 2013;8(5):e63927.

  36. 36.

    Gerdol M, Moro GD, Manfrin C, Milandri A, Riccardi E, Beran A, et al. RNA sequencing and de novo assembly of the digestive gland transcriptome in Mytilus galloprovincialis fed with toxinogenic and non-toxic strains of Alexandrium minutum. BMC Res Notes. 2014;7:722.

  37. 37.

    Zhang L, Li L, Zhu Y, Zhang G, Guo X. Transcriptome analysis reveals a rich gene set related to innate immunity in the Eastern oyster (Crassostrea virginica). Mar Biotechnol. 2014;16(1):17–33.

  38. 38.

    Huan P, Wang H, Liu B. Transcriptomic analysis of the clam Meretrix meretrix on different larval stages. Mar Biotechnol. 2012;14(1):69–78.

  39. 39.

    Sokolova IM. Apoptosis in molluscan immune defense. Invertebrate Surviv J. 2009;6:49–58.

  40. 40.

    Lee MS, Kim Y-J. Signaling pathways downstream of pattern-recognition receptors and their cross talk. Annu Rev Biochem. 2007;76:447–80.

  41. 41.

    O’Neill LAJ, Golenbock D, Bowie AG. The history of toll-like receptors - redefining innate immunity. Nat Rev Immunol. 2013;13:453–60.

  42. 42.

    Boone BA, Lotze MT. Targeting damage-associated molecular pattern molecules (DAMPs) and DAMP receptors in melanoma. In: Thurin M, Marincola FM, editors. Molecular diagnostics for melanoma: methods in molecular biology, vol. 1102. New York: Humana Press; 2014. p. 537–52.

  43. 43.

    Vos P, Lazarjani HA, Poncelet D, Faas MM. Polymers in cell encapsulation from an enveloped cell perspective. Adv Drug Deliv Rev. 2014;67-68:15–34.

  44. 44.

    Underhill DM, Goodridge HS. Information processing during phagocytosis. Nat Rev Immunol. 2012;12(7):492–502.

  45. 45.

    Flannagan RS, Jaumouillé V, Grinstein S. The cell biology of phagocytosis. Annu Rev Pathol. 2012;7:61–98.

  46. 46.

    Wenger Y, Buzgariu W, Reiter S, Galliot B. Injury-induced immune responses in Hydra. Semin Immunol. 2014;26:277–94.

  47. 47.

    Ivanina AV, Hawkins C, Sokolova IM. Immunomodulation by the interactive effects of cadmium and hypercapnia in marine bivalves Crassostrea virginica and Mercenaria mercenaria. Fish Shellfish Immunol. 2014;37:299–312.

  48. 48.

    Carreiro-Silva M, Cerqueira T, Godinho A, Caetano M, Santos RS, Bettencourt R. Molecular mechanisms underlying the physiological responses of the cold-water coral Desmophyllum dianthus to ocean acidification. Coral Reefs. 2014;33:465–76.

  49. 49.

    Li J, Chen J, Zhang Y, Yu Z. Expression of allograft inflammatory factor-1 (AIF-1) in response to bacterial challenge and tissue injury in the pearl oyster, Pinctada martensii. Fish Shellfish Immunol. 2013;34:365–71.

  50. 50.

    Labreuche Y, Lambert C, Soudant P, Boulo V, Huvet A, Nicolas J-L. Cellular and molecular hemocyte response of the Pacific oyster, Crassostrea gigas, following bacterial infection with Vibrio aestuarianus strain 01/32. Microbes Infect. 2006;8(12-13):2715–24.

  51. 51.

    Green TJ, Barnes AC. Bacterial diversity of the digestive gland of Sydney rock oysters, Saccostrea glomerata infected with the paramyxean parasite, Marteilia sydneyi. J Appl Microbiol. 2010;109:613–22.

  52. 52.

    Bathige SDNK, Umasuthan N, Whang I, Lim B-S, Won SH, Lee J. Antibacterial activity and immune responses of a molluscan macrophage expressed gene-1 from disk abalone, Haliotis discus discus. Fish Shellfish Immunol. 2014;39:263–72.

  53. 53.

    Li Q, Wang X, Korzhev M, Schröder HC, Link T, Tahir MN, et al. Potential biological role of laccase from the sponge Suberites domuncula as an antibacterial defense component. Biochim Biophys Acta. 1850;2015:118–28.

  54. 54.

    Ashe PC, Berry MD. Apoptotic signaling cascades. Prog Neuropsychopharmacol Biol Psychiatry. 2003;27:199–214.

  55. 55.

    Huo J, Xu S, Lin B, Chng W-J, Lam K-P. Fas apoptosis inhibitory molecule is upregulated by IGF-1 signaling and modulates Akt activation and IRF4 expression in multiple myeloma. Leukemia. 2013;27:1165–71.

  56. 56.

    Langenau DM, Jette C, Berghmans S, Palomero T, Kanki JP, Kutok JL, et al. Suppression of apoptosis by bcl-2 overexpression in lymphoid cells of transgenic zebrafish. Blood. 2005;105(8):3278–85.

  57. 57.

    Dubrez-Daloz L, Dupoux A, Cartier J. IAPs: more than just inhibitors of apoptosis proteins. Cell Cycle. 2008;7(8):1036–46.

  58. 58.

    Milan M, Pauletto M, Patarnello T, Bargelloni L, Marin MG, Matozzo V. Gene transcription and biomarker responses in the clam Ruditapes philippinarum after exposure to iboprofen. Aquat Toxicol. 2013;126:17–29.

  59. 59.

    Cheriyath V, Leaman DW, Borden EC. Emerging roles of FAM14 family members (G1P3/ISG 6-16 and ISG12/IFI27) in innate immunity and cancer. J Interferon Cytokine Res. 2011;31(1):173–81.

  60. 60.

    Nagaoka-Yasuda R, Matsuo N, Perkins B, Limbaeck-Stokin K, Mayford M. An RNAi-based genetic screen for oxidative stress resistance reveals retinol saturase as a mediator of stress resistance. Free Radic Biol Med. 2007;43:781–8.

  61. 61.

    Locksley RM, Killeen N, Lenardo MJ. The TNF and TNF receptor superfamilies: integrating mammalian biology. Cell. 2001;104:487–501.

  62. 62.

    Zheng C-q, Jeswin J, Shen K-l, Lablche M, Wang K-j. Detrimental effect of CO2-driven seawater acidification on a crustacean brine shrimp, Artemia sinica. Fish Shellfish Immunol. 2015;43:181–90.

  63. 63.

    Luo Y, Li C, Landis AG, Wang G, Stoeckel J, Peatman E. Transcriptomic profiling of differential responses to drought in two freshwater mussel species, the giant floater Pyganodon grandis and the pondhorn Uniomerus tetralasmus. PLoS ONE. 2014;9(2):e89481.

  64. 64.

    Bedard K, Krause K-H. The NOX family of ROS-generating NADPH oxidases: physiology and pathophysiology. Physiol Rev. 2007;87(1):245–313.

  65. 65.

    Chen Q, Vazquez EJ, Moghaddas S, Hoppel CL, Lesnefsky EJ. Production of reactive oxygen species by mitochondria. J Biol Chem. 2003;278(38):36027–31.

  66. 66.

    Munro D, Pichaud N, Paquin F, Kemeid V, Blier PU. Low hydrogen peroxide production in mitochondria of the long-lived Arctica islandica: underlying mechanisms for slow aging. Aging Cell. 2013;12:584–92.

  67. 67.

    Sussarellu R, Dudognon T, Fabioux C, Soudant P, Moraga D, Kraffe E. Rapid mitochondrial adjustments in response to short-term hypoxia and re-oxigenation in the Pacific oyster, Crassostrea gigas. J Exp Biol. 2013;216:1561–9.

  68. 68.

    Matés JM, Pérez-Gómez C. Castro INd. Antioxidant enzymes and human diseases. Clin Biochem. 1999;32(8):595–603.

  69. 69.

    Birben E, Sahiner UM, Sackesen C, Erzurum S, Kalayci O. Oxidative stress and antioxidant defense. World Allergy Organ J. 2012;5(1):9–19.

  70. 70.

    Hayes JD, McLellan LI. Glutathione and glutathione-dependent enzymes represent a co-ordinated regulated defence against oxidative stress. Free Radic Res. 1999;31(4):273–300.

  71. 71.

    Sussarellu R, Fabioux C, Sanchez MC, Goïc NL, Lambert C, Soudant P, et al. Molecular and cellular response to short-term oxygen variations in the Pacific oyster Crassostrea gigas. J Exp Mar Biol Ecol. 2012;412:87–95.

  72. 72.

    Tomanek L, Zuzow MJ, Ivanina AV, Beniash E, Sokolova IM. Proteomic response to elevated PCO2 levels in eastern oysters, Crassostrea virginica: evidence for oxidative stress. J Exp Biol. 2011;214:1836–44.

  73. 73.

    Harms L, Frickenhaus S, Schiffer M, Mark FC, Storch D, Held C, et al. Gene expression profiling in gills of the great spider crab Hyas araneus in response to ocean acidification and warming. BMC Genomics. 2014;15:789.

  74. 74.

    Healy JM. The mollusca. In: Anderson DT, editor. Invertebrate zoology. 2nd ed. Victoria, Australia: Oxford University Press; 2001. p. 120–71.

  75. 75.

    Ruppert EE, Fox RS, Barnes RD. Invertebrate zoology: a functional evolutionary approach. 7th ed. United States: Thomson Brooks/Cole; 2004.

  76. 76.

    Lucas J. Bivalves. In: Lucas JS, Southgate PC, editors. Aquaculture: farming aquatic animals and plants. Oxford: Fishing News Books; 2003. p. 443–66.

  77. 77.

    Gosling E. Bivalve molluscs: biology, ecology and culture. Oxford: Fishing News Books; 2003.

  78. 78.

    Werner GDA, Gemmel P, Grosser S, Hamer R, Shimeld SM. Analysis of a deep transcriptome from the mantle tissue of Patella vulgata Linnaeus (mollusca: gastropoda: patellidae) reveals candidate biomineralising genes. Mar Biotechnol. 2013;15:230–43.

  79. 79.

    Joubert C, Piquemal D, Marie B, Manchon L, Pierrat F, Zanella-Cléon I, et al. Transcriptome and proteome analysis of Pinctada margaritifera calcifying mantle and shell: focus on biomineralization. BMC Genomics. 2010;11:613.

  80. 80.

    Freer A, Bridgett S, Jiang J, Cusack M. Biomineral proteins from Mytilus edulis mantle tissue transcriptome. Mar Biotechnol. 2014;16:34–45.

  81. 81.

    Clark MS, Thorne MAS, Vieira FA, Cardoso JCR, Power DM, Peck LS. Insights into shell deposition in the Antarctic bivalve Laternula elliptica: gene discovery in the mantle transcriptome using 454 pyrosequencing. BMC Genomics. 2010;11:362.

  82. 82.

    Gardner LD, Mills D, Wiegand A, Leavesley D, Elizur A. Spatial analysis of biomineralization associated gene expression from the mantle organ of the pearl oyster Pinctada maxima. BMC Genomics. 2011;12:455.

  83. 83.

    Hüning AK, Melzner F, Thomsen J, Gutowska MA, Krämer L, Frickenhaus S, et al. Impacts of seawater acidification on mantle gene expression patterns of the Baltic Sea blue mussel: implications for shell formation and energy metabolism. Mar Biol. 2013;160(8):1845–61.

  84. 84.

    Ren G, Wang Y, Qin J, Tang J, Zheng X, Li Y. Characterization of a novel carbonic anhydrase from freshwater pearl mussel Hyriopsis cumingii and the expression profile of its transcript in response to environmental conditions. Gene. 2014;546:56–62.

  85. 85.

    Wei L, Wang Q, Ning X, Mu C, Wang C, Cao R, et al. Combined metabolome and proteome analysis of the mantle tissue from Pacific oyster Crassostrea gigas exposed to elevated pCO2. Comp Biochem Physiol Part D Genomics Proteomics. 2015;13:16–23.

  86. 86.

    Henry RP. Multiple roles of carbonic anhydrase in cellular transport and metabolism. Annu Rev Physiol. 1996;58(1):523–38.

  87. 87.

    Bai Z, Zheng H, Lin J, Wang G, Li J. Comparative analysis of the transcriptome in tissues secreting purple and white nacre in the pearl mussel Hyriopsis cumingii. PLoS ONE. 2013;8(1):e53617.

  88. 88.

    Tetreau G, Cao X, Chen Y-R, Muthukrishnan S, Haobo J, Blissard GW et al. Overview of chitin metabolism enzymes in Manduca sexta: identification, domain organization, phylogenetic analysis and gene expression. Insect Biochem Mol Biol. 2015;62:114–26.

  89. 89.

    Xiao R, Xie L-P, Lin J-Y, Li C-H, Chen Q-X, Zhou H-M, et al. Purification and enzymatic characterization of alkaline phosphatase from Pinctada fucata. J Mol Catal B Enzym. 2002;17:65–74.

  90. 90.

    Li S, Xie L, Zhang C, Zhang Y, Gu M, Zhang R. Cloning and expression of a pivotal calcium metabolism regulator: calmodulin involved in shell formation from pearl oyster (Pinctada fucata). Comp Biochem Physiol B Biochem Mol Biol. 2004;138:235–43.

  91. 91.

    Fang Z, Yan Z, Li S, Wang Q, Cao W, Xu G, et al. Localization of calmodulin and calmodulin-like protein and their functions in biomineralization in P. fucata. Prog Nat Sci. 2008;18:405–12.

  92. 92.

    Wright JM, Parker LM, O’Connor WA, Williams M, Kube P, Ross PM. Populations of Pacific oysters Crassostrea gigas respond variably to elevated CO2 and predation by Morula marginalba. Biol Bull. 2014;226:269–81.

  93. 93.

    Fitzer SC, Phoenix VR, Cusack M, Kamenos NA. Ocean acidification impacts mussel control on biomineralisation. Sci Rep. 2014;4:6218.

  94. 94.

    Liang P, MacRae TH. Molecular chaperones and the cytoskeleton. J Cell Sci. 1997;110:1431–40.

  95. 95.

    Janmey PA. The cyotskeleton and cell signaling: component localization and mechanical coupling. Physiol Rev. 1998;78(3):763–81.

  96. 96.

    Morse EM, Brahme NN, Calderwood DA. Integrin cytoplasmic tail interactions. Biochemistry (Mosc). 2014;53(5):810–20.

  97. 97.

    Frantz C, Stewart KM, Weaver VM. The extracellular matrix at a glance. J Cell Sci. 2010;123:4195–200.

  98. 98.

    Brown BN, Badylak SF. Extracellular matrix and an inductive scaffold for functional tissue reconstruction. Transl Res. 2014;163(4):268–85.

  99. 99.

    Mouw JK, Ou G, Weaver VM. Extracellular matrix assembly: a mutiscale deconstruction. Nat Rev Mol Cell Biol. 2014;15:771–85.

  100. 100.

    Murphy G. Tissue inhibitors of metalloproteinases. BMC Genome Biology. 2011;12:233.

  101. 101.

    Yan F, Jiao Y, Deng Y, Du X, Huang R, Wang Q, et al. Tissue inhibitor of metalloproteinase gene from pearl oyster Pinctada martensii participates in nacre formation. Biochem Biophys Res Commun. 2014;450:300–5.

  102. 102.

    Montagnani C, Roux FL, Berthe F, Escoubas J-M. Cg-TIMP, an inducible tissue inhibitor of metalloproteinase from the Pacific oyster Crassostrea gigas with a potential role in wound healing and defense mechanisms. FEBS Lett. 2001;500(1-2):64–70.

  103. 103.

    Huo L, Scarpulla RC. Mitochondrial DNA instability and peri-implantation lethality associated with targeted disruption of nuclear respiratory factor 1 in mice. Mol Cell Biol. 2001;21(2):644–54.

  104. 104.

    Korobeinikova AV, Garber MB, Gongadze GM. Ribosomal proteins: structure, function, and evolution. Biochemistry (Mosc). 2012;77(6):562–74.

  105. 105.

    Parsyan A, Svitkin Y, Shahbazian D, Gkogkas C, Lasko P, Merrick WC, et al. mRNA helicases: the tacticians of translational control. Nat Rev Mol Cell Biol. 2011;12:235–45.

  106. 106.

    Krecic AM, Swanson MS. hnRNP complexes: composition, structure, and function. Curr Opin Cell Biol. 1999;11:363–71.

  107. 107.

    Xiong X, Feng Q, Xie L, Zhang R. Cloning and characterization of a heterogeneous nuclear ribonucleoprotein homolog from pearl oyster, Pinctada fucata. Acta Biochim Biophys Sin. 2007;39(12):955–63.

  108. 108.

    Beninger PG, Veniot A. The oyster proves the rule: mechanisms of pseudofeces transport and rejection on the mantle of Crassostrea virginica and C. gigas. Mar Ecol Prog Ser. 1999;190:179–88.

  109. 109.

    Millar RH. Notes on the mechanism of food movement in the gut of the larval oyster, Ostrea edulis. Q J Microsc Sci. 1955;96(4):539–44.

  110. 110.

    Ward TJ. Effect of cadmium on particle clearance by the Sydney rock oyster, Saccostrea commercialis (I. & R.). Mar Freshw Res. 1982;33(4):711–5.

  111. 111.

    Micallef S, Tyler PA. Effect of mercury and selenium on the gill function of Mytilus edulis. Mar Pollut Bull. 1990;21(6):288–92.

  112. 112.

    Amos LA. The tektin family of microtubule-stabilizing proteins. Genome Biol. 2008;9(7):229.

  113. 113.

    Dineshram R, Wong KKW, Xiao S, Yu Z, Qian PY, Thiyagarajan V. Analysis of Pacific oyster larval proteome and its response to high-CO2. Mar Pollut Bull. 2012;64:2160–7.

  114. 114.

    Kingtong S, Kellner K, Bernay B, Goux D, Sourdaine P, Berthelin CH. Proteomic identification of protein associated to mature spermatozoa in the Pacific oyster Crassostrea gigas. J Proteomics. 2013;82:81–91.

  115. 115.

    Yang Y, Cochran DA, Gargano MD, King I, Samhat NK, Burger BP, et al. Regulation of flagellar motility by the conserved flagellar protein CG34110/Ccdc135/FAP50. Mol Biol Cell. 2011;22:976–87.

  116. 116.

    Horani A, Brody SL, Ferkol TW, Shoseyov D, Wasserman MG, Ta-shma A, et al. CCDC65 mutation causes primary ciliary dyskinesia with normal ultrastructure and hyperkinetic cilia. PLoS ONE. 2013;8(8):e72299.

  117. 117.

    Barros P, Sobral P, Range P, Chícharo L, Matias D. Effects of sea-water acidification on fertilization and larval development of the oyster Crassostrea gigas. J Exp Mar Biol Ecol. 2013;440:200–6.

  118. 118.

    Sung C-G, Kim TW, Park Y-G, Kang S-G, Inaba K, Shiba K, et al. Species and gamete-specific fertilization success of two sea urchins under near future levels of pCO2. J Mar Syst. 2014;137:67–73.

  119. 119.

    Chien Y-H, Werner ME, Stubbs J, Joens MS, Li J, Chien S, et al. Bbof1 is required to maintain cilia orientation. Development. 2013;140:3468–77.

  120. 120.

    Yen H-J, Tayeh MK, Mullins RF, Stone EM, Sheffield VC, Slusarski DC. Bardet-Biedl syndrome genes are important in retrograde intracellular trafficking and Kuppfer’s vesicle cilia function. Hum Mol Genet. 2006;15(5):667–77.

  121. 121.

    Wilkinson CJ, Carl M, Harris WA. Cep70 and Cep131 contribute to ciliogenesis in zebrafish embryos. BMC Cell Biol. 2009;10:17.

  122. 122.

    Veleri S, Manjunath SH, Fariss RN, May-Simera H, Brooks M, Foskett TA, et al. Ciliopathy-associated gene Cc2d2a promotes assembly of subdistal appendages on the mother centriole during cilia biogenesis. Nat Commun. 2014;5:4207.

  123. 123.

    Cappello T, Mauceri A, Corsaro C, Maisano M, Parrino V, Paro GL, et al. Impact of environmental pollution on caged mussels Mytilus galloprovincialis using NMR-based metabolomics. Mar Pollut Bull. 2013;77:132–9.

  124. 124.

    Karpenko AA, Ivanovsky YA. Effect of very low doses of γ radiation on motility of gill ciliated epithelia of Mytilus edulis. Radiat Res. 1993;133:108–10.

Download references


We would like to thank the following people and organisations: Mr John Wright for giving us access to his CO2 and temperature experiments for the samples; Mr Stephan O’Connor, Mr Kyle Johnston and the rest of the DPI hatchery team for their help and advice with oyster husbandry; Dr David Schoeman for statistical advice and writing the R based script used to produce the heatmaps; Dr Ido Bar for carrying out the GO enrichment analysis; CSIRO for a postgraduate studentship for NGE; the National Computational Infrastructure Specialised facility in Bioinformatics (Barrine@UQ) for Barrine support.


The project was funded by the Fisheries Research & Development Corporation (FRDC –, the Australian Seafood Cooperative Research Centre (CRC – (2011/718 to NGE, AE and WAO), and the University of the Sunshine Coast, Australia. NGE was supported by an Australian Seafood CRC and University of the Sunshine Coast scholarship, as well as a scholarship from Queensland Education and Training International (QETI). The funding parties had no role in study design, data collection, analysis and interpretation of data, decision to publish, or preparation of the manuscript.

Availability of data and materials

The raw CO2 reads supporting the conclusions of this article are available from the sequence read archive and can be accessed under the SRA study accession number SRP055052.

Authors’ contributions

NGE: designed the study, collected samples from the stress experiments, carried out all relevant laboratory and bioinformatics work, as well as the differential gene expression analysis, drafted manuscript. WAO: designed study, advised on all experiments, provided oysters and necessary equipment for experiments, reviewed manuscript. ANW: wrote the in-house script used to remove ERCC and phiX transcripts from the transcriptome and advised on linux based errors. AE: designed study, reviewed manuscript. All authors have read and accepted the manuscript.

Competing interests

The authors declare that they have no competing interests.

Author information

Correspondence to Abigail Elizur.

Additional files

Additional file 1:

Detailed methods of experimental exposure trials and previous sequencing. (DOCX 23 kb)

Additional file 2: Table S1.

Bowtie alignment statistics. Table shows the number of raw reads and reads surviving the processing, with the total alignment percentage based on the post-processed reads aligning to the S. glomerata reference transcriptome, using Bowtie. (DOCX 12 kb)

Additional file 3: Figure S1.

Variance versus mean plot for each Ng group (C1). This plot shows the mean-variance relationship (using polynomial regression) for each isoform (Ng) group of condition 1 (control samples). Mapping ambiguity clusters were produced with RSEM (rsem-generate-ngvector), while the plot was visualised in R using EBSeq’s PolyFitPlot function. (PDF 260 kb)

Additional file 4: Figure S2.

Variance versus mean plot for each Ng group (C2). This plot shows the mean-variance relationship (using polynomial regression) for each isoform (Ng) group of condition 2 (elevated CO2 samples). Mapping ambiguity clusters were produced with RSEM (rsem-generate-ngvector), while the plot was visualised in R using EBSeq’s PolyFitPlot function. (PDF 268 kb)

Additional file 5: Figure S3.

Quantile-quantile plot. QQ-plots show the fitted Beta prior distributions within each condition and each Ig group (uncertainty group) and were visualised in R using EBSeq’s QQP function. (PDF 337 kb)

Additional file 6: Figure S4.

Density plot. Plot shows the prior distribution fit within each condition and each Ig group, visualised in R using EBSeq’s DenNHist function. (PDF 38 kb)

Additional file 7: Table S2.

S. glomerata DE transcripts. List of 1626 DE transcripts determined with Bowtie-RSEM-EBSeq, using a FDR threshold of 0.05. Sequence descriptions are based on blast homology searches against the NCBI nr database (e-value cut-off: 10-5, hit number threshold: 25), and on InterProScan domain/family information. Posterior fold change (FC) was based on the normalised data, whereas real FC was based on the raw data. C1 stands for control, C2 for treatment condition. (XLSX 124 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Ertl, N.G., O’Connor, W.A., Wiegand, A.N. et al. Molecular analysis of the Sydney rock oyster (Saccostrea glomerata) CO2 stress response. Clim Chang Responses 3, 6 (2016).

Download citation


  • Saccostrea glomerata
  • Sydney rock oyster
  • Molluscs
  • RNA-seq
  • Stress
  • Carbon dioxide
  • Immunity
  • Biomineralisation