Ocean acidification at a coastal CO2 vent induces expression of stress-related transcripts and transposable elements in the sea anemone Anemonia viridis.

Ocean acidification threatens to disrupt interactions between organisms throughout marine ecosystems. The diversity of reef-building organisms decreases as seawater CO2 increases along natural gradients, yet soft-bodied animals, such as sea anemones, are often resilient. We sequenced the polyA-enriched transcriptome of adult sea anemone Anemonia viridis and its dinoflagellate symbiont sampled along a natural CO2 gradient in Italy to assess stress levels in these organisms. We found that about 3.1% of the anemone transcripts, but <1% of the Symbiodinium sp. transcripts were differentially expressed. Processes enriched at high seawater CO2 were linked to cellular stress and inflammation, including significant up-regulation of protective cellular functions and down-regulation of metabolic pathways. Transposable elements were differentially expressed at high seawater CO2, with an extreme up-regulation (> 100-fold) of the BEL-family of long terminal repeat retrotransposons. Seawater acidified by CO2 generated a significant stress reaction in A. viridis, but no bleaching was observed and Symbiodinium sp. appeared to be less affected. These observed changes indicate the mechanisms by which A. viridis acclimate to survive chronic exposure to ocean acidification conditions. We conclude that many organisms that are common in acidified conditions may nevertheless incur costs due to hypercapnia and/or lowered carbonate saturation states.

Introduction

Reef-forming cnidarians are in global decline due to rapidly increasing levels of atmospheric CO2 [1], yet non-calcified cnidarians appear to be more resilient [2-6]. Since the Industrial Revolution, average surface ocean pH has decreased by 0.1 units and is projected to drop by further 0.2–0.4 units by the end of this century, depending on CO2 emission scenarios [7].

Ocean acidification is expected to lead to widespread marine biodiversity loss and cause major shifts in coastal ecosystems as some species are vulnerable whereas others are more hardy [2, 8]. Many calcifying corals are affected as seawater carbonate saturation levels fall and reefs can be corroded if carbonate levels become undersaturated [9-11]. Some cnidarians grow well in high CO2 environments, but corrosive waters can weaken the skeletons of those that are calcified and hypercapnia can increase metabolic cost as organisms allocate resources when coping with the changes [12-15]. Surveys in areas with naturally high levels of CO2 have shown that non-calcified cnidarians (e.g. soft corals, sea anemones and jellyfish) are more resilient to the effects of ocean acidification [3, 5]. However, the molecular mechanisms underlying this resilience are largely unknown.

Physiological measurements of sea anemones exposed to high CO2 conditions both in situ and during laboratory studies show that sea anemones with symbiotic algae increase their primary productivity at high CO2 due to enhanced carbon availability [3, 16–20]. To better understand the physiological plasticity and adaptive potential of cnidarians to high CO2, several transcriptome-level studies have been performed, mostly using corals [21-24]. The results of these studies imply that the strongest response is to an acute increase in CO2 levels [20, 23, 24], and that cnidarians exposed long-term to high CO2 conditions appear to be more resilient [3, 20, 22, 24]. Although the impact of high CO2 can be assessed under laboratory conditions, investigations at marine CO2 seeps are useful for studies of the long-term effects of ocean acidification in situ in the whole ecosystem [25], and therefore complement laboratory experiments.

Beyond changes in expression of protein-coding genes, environmental stress can induce the activity of transposable elements that are able to influence both coding and non-coding potential of a genome through genomics rearrangements [26-29]. Transposable elements contribute to the ability to tolerate both biotic and abiotic stresses and could potentially assist the organism to acclimatize or adapt to new environmental conditions [27, 29]. A relatively high proportion of the sea anemone genomes is composed of repeats and transposable elements [30, 31]. However, the involvement of transposable elements in response to stress has not been extensively studied in cnidarians. An increase in transposable element activity upon bleaching in the coral Montastraea faveolata is the only observed case so far [32].

As with many organisms that are able to live intertidally, Anemonia viridis is physiologically robust and can cope with stressful physicochemical conditions. It has a wider variety of antioxidant enzymes, such as catalases, peroxidases and superoxide dismutases (SOD) than non-symbiotic anemones, and these seem to make this species more resilient to stressors [33]. In addition, the symbiotic cnidarians show an ancient expansion of transposable elements compared to non-symbiotic cnidarians [30], which could contribute in adaptation to a greater range of environmental changes.

The sea anemone A. viridis exposed to natural ocean acidification conditions appears to be acclimatized to high CO2 [3, 18]. The goal of this research was to elucidate the changes in global gene expression patterns in the host and its symbiont that might be associated with the acclimatization of A. viridis to high CO2 conditions in situ. We recently reported whole genome and mitochondrial genome sequencing of A. viridis, and assessed the expression of small RNAs (miRNAs and piRNAs) and mitochondrial RNAs at in situ ocean acidification conditions [34, 35]. Here, we report transcriptomic studies performed on A. viridis and its symbiont Symbiodinium sp., living at a natural CO2 gradient off Vulcano Island, Italy, where these organisms are very common in seawater between pH 7.6 and pH 8.2 [3, 18]. We investigated the transcriptomic response of A. viridis living in high seawater CO2 conditions (pH 7.6) to find out whether these conditions were stressful for A. viridis.

Results

Anemonia viridis reference transcriptome

Adult polyps of the sea anemone A. viridis were sampled from three different locations along a natural CO2 gradient near a set of volcanic seeps off Vulcano Island, Sicily (Italy), representing seawater pH conditions 8.2, 7.9, and 7.6 (Fig 1A and 1B). Twelve individual polyps, four from each location, were subjected to total RNA extraction, mRNA isolation and subsequently whole transcriptome sequencing using the Ion Torrent PGM platform (Fig 1C). We generated 3.4 billion nucleotide (nt) sequence data (S1 Table). Sequences obtained from all the individuals were pooled, quality filtered and then used for transcriptome assembly by Trinity [36], yielding 244,294 contigs with an average length of 538 nt (including 41,197 isoforms). Contigs were clustered based on 90% sequence similarity, resulting in 154,015 contigs (including 9,510 isoforms) after filtering for contaminants. About 40% of the contigs had blast hits (38% in nr database and 40% in Swiss-Prot/UniProtKB, e-value < 10-5) and 57% of these were annotated with the Blast2GO pipeline B2G4Pipe [37]. This assembly was used as a reference transcriptome in further analyses (see additional assembly metrics in S1 Table).

Fig 1

Sampling site overview and data analysis workflow.

(A) The sea anemone Anemonia viridis growing at high CO2 (mean pH 7.6) at 1 m depth. (B) Sampling site location at Levante Bay (Vulcano Island, Italy); 1 – low seawater pH site with average pH 7.6; 2 – low seawater pH site with average pH 7.9 and 3 – reference site with average seawater pH 8.2. (C) General overview of laboratory methods and data analyses performed. V – vent site. (Photo credit: A - Demetris Kletou, B - Ilona Urbarova).

The reference transcriptome was divided into two fractions using PSyTranS tool (https://github.com/sylvainforet/psytrans) [38], one representing A. viridis and the other its symbiotic dinoflagellate Symbiodinium sp. (hereafter called the symbiont). This resulted in 90,535 A. viridis contigs with average GC-content 41% and 63,480 symbiont contigs with average GC-content 57% (Fig 2, S1 Table), which were subsequently used in separate transcriptome analyses for each species. The presence of symbiont at all sample locations was confirmed by PCR amplification of the symbiont-specific nuclear apx gene [39] (S2 Table), and from gene expression profiling analysis presented below. We infer that the host-symbiont relationship is not disrupted at low seawater pH conditions.

Fig 2

Separation of host and symbiont contigs using PSyTranS software.

Host sequences (Anemonia viridis) were separated from the symbiont sequences (Symbiodinium sp.) using BLAST and Support Vector Machine (SVM) classification implemented in PSyTranS software tool.

Anemonia viridis and its symbiont respond differently to low seawater pH

Differential gene expression and its variability among specimens from the same conditions was assessed separately for the host and the symbiont. Quality filtered sequence reads were first aligned to the reference transcriptome, followed by extraction of mapped reads belonging to host and symbiont contigs. We filtered the data set, keeping only transcripts with minimum 10 aligned reads each in at least four samples, and performed differential gene expression analysis accounting for two different variables; pH condition and day of sampling (glm edgeR; S3 Table) [40]. Because the host/ symbiont contig mapping ratio was found relatively consistent among samples (1.567 ± 0.159), we deemed it unnecessary to normalize our data to this parameter. We found 1233 (~1.4%) of the A. viridis transcripts and 291 (~ 0.5%) of the symbiont transcripts to be differentially expressed (DE) among the assessed pH conditions (FDR < 0.05). Nearly all the DE-transcripts appeared to be affected by the pH condition variable only (S4 Table). Hierarchical clustering clearly showed for most individuals that A. viridis DE-transcripts from specimens at the lowest pH condition studied (pH 7.6) appeared on a separate branch compared to that of pH 7.9 and pH 8.2 (Fig 3B). In the symbiont, however, a closer association could be observed between symbiont specimens at pH 7.6 and pH 8.2 (Fig 3B). We infer that A. viridis and the symbiont each responded to pH 7.6 by a different set of DE-transcripts compared to that of pH 7.9 (Fig 3, S1 Fig).

Fig 3

Differential gene expression profiles of Anemonia viridis and Symbiodinium sp. at the sampling sites.

Differential expression (DE) pipeline using a glm edgeR approach was applied to account for both the day of sampling and the different pH where samples were taken. The DE analysis was performed separately for A. viridis and Symbiodinium sp. (A) Principal component analysis (PCA) plots show clustering similarity of individual samples. Numbers in the PCA plots represent different individuals sequenced. (B) Shown are heatmaps with hierarchically clustered, significantly differentially expressed (DE) transcripts between the sampling sites separately for A. viridis and Symbiodinium sp. (C) Venn diagrams visualize the number of private and shared DE-transcripts at pH 7.6 and pH 7.9 compared to normal conditions (pH 8.2). A. viridis contained 526 private DE-transcripts at pH 7.6 and 318 private DE-transcripts at pH 7.9. The symbiont contained 61 and 123 private DE-transcripts at pH 7.6 and pH 7.9, respectively. Venn diagrams were created using venneuler in R software.

Down-regulation of transcripts from seawater pH 8.2 to either pH 7.6 or pH 7.9 dominated in A. viridis and symbiont data sets. A. viridis showed more transcripts as DE at pH 7.6 compared to pH 7.9, which was contrary to that of the symbiont (Fig 3C). In A. viridis, we noted 294 up-regulated and 315 down-regulated transcripts at pH 7.6, and 165 up-regulated and 236 down-regulated transcripts at pH 7.9. In the symbiont however, 28 transcripts were up-regulated and 56 were down-regulated at pH 7.6, versus 36 up-regulated and 110 down-regulated transcripts at pH 7.9. Significant amounts of DE-transcripts in A. viridis and the symbiont were found to be specific (private) to only one pH condition when compared to pH 8.2 (Fig 3C). From these analyses it therefore appears that low seawater pH has a more significant effect on gene expression in the sea anemone compared to that of the symbiont.

LTR-retrotransposons are heavily up-regulated at low seawater pH

Homology search for transposable elements in our reference transcriptome resulted in 1397 significant contig hits in the RepBase database [41] (1209 in A. viridis and 188 in the symbiont). The majority (~70%) of these hits showed similarities to known retrotransposons, while DNA transposons were found to constitute only about 30% of all transposons. We found 10 transposable elements (> 200 nt) from eight A. viridis contigs as differentially expressed at low seawater pH 7.6 compared to that of normal pH 8.2 conditions (glm edgeR, FDR < 0.05), and most of them were up-regulated (Table 1). Several different transposable element sub-classes were represented, but the most dramatic changes were noted for the long terminal repeat (LTR) retrotransposons. Here, one element of the BEL-family showed extreme up-regulation at pH 7.6 (> 100-fold; Table 1). Interestingly, at pH 7.9 only three transposable elements were found to be differentially expressed. Two Gypsy elements (LTR retrotransposons) were four- and six-fold up-regulated and the chicken repeat 1 (CR1) element (annotated in N. vectensis) was 11-fold down-regulated. No transposable element could be observed as differentially expressed in the symbiont. We conclude that several classes of transposable elements, and especially the LTR retrotransposons, are strongly activated in A. viridis at low seawater pH conditions.

Table 1

Differentially expressed transposable elements (TEs) at pH 7.6 compared to normal pH 8.2 conditions in Anemonia viridis.

Transcript	TE region	Fold change	p-value	False Discovery Rate (FDR)	Transposon name	Class of transposon	Specie
TR33113|c0_g2_i1	219–2079	112.55	1.53E-07	2.14E-04	BEL11-I_AG	LTR Retrotransposon	Anopheles gambiae
TR82105|c3_g1_i1	2757–4658	5.75	1.29E-04	1.87E-02	Gypsy-43-I_NV	LTR Retrotransposon	Nematostella vectensis
TR15677|c3_g2_i1	1636–2329	6.78	7.32E-06	3.07E-03	EnSpm-5_CCri	DNA transposon	Chondrus crispus
TR15677|c3_g2_i1	1372–1599	6.98	1.14E-07	1.69E-04	Copia-33_GM-I	LTR Retrotransposon	Glycine max
TR15677|c3_g2_i1	1177–1326	6.08	1.01E-05	3.76E-03	Copia-33_GM-I	LTR Retrotransposon	Glycine max
TR82105|c3_g1_i1	3–2337	6.04	3.47E-07	3.85E-04	Gypsy-21_Adi-I	LTR Retrotransposon	Acropora digitifera
TR53336|c0_g1_i1	339–737	5.95	6.07E-05	1.19E-02	Polinton-2_PH	DNA transposon	Parhyale hawaiensis
TR41995|c0_g2_i1	60–509	3.33	5.01E-04	4.35E-02	EnSpm-5_CCri	DNA transposon	Chondrus crispus
TR82351|c2_g27_i1	2–301	3.39	5.31E-04	4.50E-02	BEL1-I_SM	LTR Retrotransposon	Schmidtea mediterranea
TR41961|c21_g6_i3	313–784	-10.88	7.01E-05	1.27E-02	DIRS-35_NV	LTR Retrotransposon	Nematostella vectensis

Taxonomically restricted genes are mostly down-regulated at low seawater pH

A significant part of our DE-transcript data set included contigs lacking annotation. We found about 55% (672 DE contigs) for A. viridis and 35% (102 DE contigs) for the symbiont with no hits in protein databases (nr and Swiss-Prot/UniProtKB, e-value < 10-5). The great majority contained ORFs, but we identified 14 DE contigs for A. viridis lacking ORFs, and without any hit to protein databases. These may represent putative long non-coding RNAs (lncRNAs) or may originate from untranslated regions (UTRs). We designated the contigs without annotation as taxonomically restricted genes (TRGs), and most of these were found to be down-regulated at low pH conditions in both A. viridis and the symbiont. For A. viridis, 44% TRGs were differentially expressed at pH 7.6 compared to pH 8.2, with 58% down-regulated. At pH 7.9, 58% of TRGs were down-regulated. In the symbiont, most TRGs (56%) were differentially expressed at pH 7.9, with 77% down-regulated. Only 50 putative protein domains in 34 TRGs of A. viridis, and 20 putative protein domains in 12 TRGs for the symbiont were identified by Hidden-Markov model search (hmmscan, e-value < 10-3).

Up-regulation of stress-related transcripts and deregulation of metabolism-related transcripts at low seawater pH

A. viridis showed significant up-regulation of stress-related transcripts at low pH, and an increase in oxidative stress response (Table 2). Furthermore, we also noted deregulation of metabolism-related transcripts at the same conditions. The symbiont showed a similar trend, but only few transcripts could be annotated and therefore results did not appear significant in a gene set enrichment analysis.

Table 2

Selected genes with significant differential expression at low seawater pH 7.6 compared to normal seawater pH 8.2 in Anemonia viridis.

Transcripts 1	Transcript	Fold Change	p-value	False Discovery Rate (FDR)	Transcript length	e-value	Blast similarity [%]
Stress-response genes
heat shock 70 kda protein (Hsp70)	TR56459|c7_g9_i1	146.23	3.72E-09	1.24E-05	1900	0	89.05
heat shock protein 90 (Hsp90)	TR15677|c3_g2_i1	6.80	4.22E-06	2.13E-03	2783	0	90.40
glucose-regulated protein 78 (Grp78)	TR44825|c0_g1_i1	4.70	1.75E-06	1.14E-03	364	4,96E-15	72.30
glucose-regulated protein 94 (Grp94)	TR41995|c0_g2_i1	3.04	5.61E-04	4.96E-02	2534	0	85.80
hypoxia up-regulated protein 1-like (HYOU1)	TR82131|c1_g7_i3	3.87	7.81E-05	1.43E-02	2474	0	76.35
cyclic AMP-responsive element-binding protein 3-like protein 3-A (CREB3L3A)	TR67696|c1_g12_i1	4.70	1.71E-05	5.27E-03	304	2.00E-10	48.04%
Antioxidant response
nuclear factor erythroid 2-related factor (Nrf2)	TR506|c0_g1_i1	5.11	2.24E-06	1.36E-03	2220	1,42E-31	59.90
Heme oxygenase 1 (HO-1)	TR9534|c0_g1_i1	8.41	3.42E-05	8.62E-03	401	3.00E-08	65.91%
Metabolism
heme oxygenase 2-like (HO-2)	TR2095|c0_g3_i1	3.92	9.77E-05	1.67E-02	299	2.90E-43	75.20
niemann-pick c 2 (Npc2d) – like	TR4599|c1_g6_i3	-3.44	8.88E-05	1.56E-02	818	4,97E-29	55.20
Cell growth and survival genes
NF-kappa-b (NF-κB) p100 subunit	TR20187|c0_g1_i1	4.65	2.49E-04	2.96E-02	2584	1,02E-71	56.55
Bcl-2-like protein 2 (BCL2L2)	TR18835|c6_g1_i1	5.88	5.28E-05	1.12E-02	1132	7,01E-58	78.00
eukaryotic translation initiation factor 2-alpha kinase 3 (EIF2AK3)	TR5221|c2_g3_i2	5.76	1.56E-05	4.94E-03	2885	0	52.65

1 Selection of the most important differentially expressed genes at pH 7.6 compared to pH 8.2 in A. viridis. If isoforms exist, the most abundant isoform was selected.

A. viridis and symbiont contigs with hits in either of the protein databases were annotated with Gene Ontology (GO) terms in order to study biological processes affected by the acidified conditions. Using the GOseq method (S5 Table) [42], a number of GO categories were found over-represented among the DE-transcripts (FDR < 0.05; S6 Table). GOseq terms enriched at pH 7.6 were manually classified into presented ancestor categories (with the most important DE-transcripts shown in Table 2 and extended listing in S7 Table). As indicated by hierarchical clustering (Fig 3), the largest expression changes of the A. viridis transcriptome were found at pH 7.6 compared to pH 7.9 and pH 8.2. GO categories enriched at pH 7.6 were mainly associated with cellular stress (S6 Table). The most significantly up-regulated transcripts were heat shock proteins from Hsp70 and Hsp90 family, molecular chaperones, and other stress response genes (Table 2). We also observed many deregulated transcripts linked to metabolic enzymes, and expression of genes enhancing cell survival was increased (Table 2, S7 Table). Signal transduction pathway members were mostly up-regulated, though transport channels appeared both up- and down-regulated in our analysis. Repression of global synthesis was supported by up-regulation of EIF2AK3 (also known as PERK), leading to reduction in translation, and by highly up-regulated ubiquitin (S7 Table). In addition, an anti-apoptotic response was observed by up-regulation of Bcl-2-like protein, as also reported in corals by Moya et al. [24] (Table 2). Taken together, these gene features indicated stress-responses at low seawater pH conditions.

For the symbiont, no GO categories appeared significantly over-represented (FDR < 0.05), though the symbiont seemed to respond similarly as the host (S8 Table). Interestingly, we found up-regulation of autophagy 8i transcript at pH 7.6 (S8 Table) and down-regulation of putative chlorophyll A-B binding protein domain by InterPro at both low pH conditions (S9 Table), potentially resulting from oxidative stress response [43]. Despite their successful detection in the reference assembly (S10 Table), the important symbiont stress-response genes (i.e. Hsp70 and 90, SODs, GR, TRX, APX and CYP450) were surprisingly not observed differentially expressed in our study (S8 Table).

The expression of genes related to symbiosis is affected at low seawater pH

Many genes involved in innate immunity, lipid metabolism, cell signalling, oxidative stress, apoptosis, autophagy and phagocytosis implicated previously in regulation of host-symbiont relationship in sea anemones [44] appeared differentially expressed at pH 7.6 in our study (Table 2, S7 Table). However, we targeted our analysis on important candidate genes for host-symbiotic relationship maintenance identified specifically in A. viridis [39]. To elucidate the involvement of these genes in response to low seawater pH, we performed BLASTN homology search of the DE-transcripts identified in our study to the A. viridis EST database [45]. The output was filtered for transcripts found exclusively up-regulated in the symbiotic and the aposymbiotic A. viridis [39]. Here, we identified 4 and 5 transcripts significantly differentially expressed at low pH from the symbiotic and aposymbiotic state (S11 and S12 Tables). The majority of the symbiotic genes showed down-regulation, and all the aposymbiotic genes showed up-regulation in our data set. Three of these transcripts were from the ‘Kern set’ genes (S11 Table), including carbonic anhydrase (CA). Significant down-regulation of CA transcript expression was observed at pH 7.6 compared to pH 8.2 (S11 Table), a result consistent with the findings by Ventura et al. from the same sampling site [20]. Interestingly, we also found a significant down-regulation of Niemann-Pick C2-like protein transcript (specifically its Npc2d variant) at pH 7.6 (Table 2).

Verification of significantly differentially expressed genes by quantitative PCR

In order to verify key findings from the RNA-seq experiment, we re-examined the expression of nine selected A. viridis DE-transcripts involved in stress responses and metabolism by the qPCR approach (S2 Fig). From initial testing of six candidates, we successfully established the ribosomal protein L12 (RPL12), beta-actin and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) transcripts as suitable reference transcripts with relatively stable expression levels between seawater pH 7.6, pH 7.9 and pH 8.2 (S3 Fig). We then examined the expression levels of the selected transcripts from the DE analysis. The main aim was to verify the observed RNA-seq expression profiles of transcripts at pH 7.6 compared to control pH 8.2, as was validated for six transcripts (Hsp70, NF-kappaB, HO-2, HYOU1, IRF2 and NPC2d) by qPCR (S2 Fig, S7 Table). Majority of assessed transcripts did not show significantly differential expression at pH 7.9 compared to pH 8.2 in RNA-seq experiment, as was also confirmed by the qPCR analysis (S2 Fig).

Discussion

Initial observations by Suggett et al. [3] showed that both the abundance of A. viridis and the photosynthetic rates of its symbiont (Symbiodinium sp. type A19) increased as seawater CO2 levels rose. There is growing evidence that ocean acidification is a threat to many hard corals, but that their less calcified relatives, e.g. soft corals, jellyfishes and sea anemones, are more resilient [4, 5]. Here, we conducted a transcriptomic response study of the sea anemone and its symbiont. We found high levels of stress-response genes at high CO2 in the sea anemone, but not in the dinoflagellate. Even though high CO2 conditions might not be benign or beneficial to some organisms, those organisms may proliferate in hypercapnic conditions because their competitors are less tolerant than they are [46]. At high CO2, A. viridis showed up-regulation of cellular stress response gene transcripts and deregulation of transcripts involved in metabolism. Furthermore, we observed high up-regulation of several LTR retrotransposons. The symbiont also responded, but fewer genes were differentially expressed. This observation is consistent with previous findings by Barshis et al. [47] and Leggat et al. [48] of symbiotic cnidarians in response to heat stress, where the symbiont was found less affected than the host. Low pH conditions neither appear to represent stressful conditions for the symbiont physiologically. Under normocapnia, the symbiont is apparently inorganic carbon (iC) limited and the increased iC availability at low pH promotes autotrophy, as discussed by Suggett et al. [3], and in more detail presented by Horwitz et al. [18].

One of our important findings was a dramatic increase in the amount of stress-related transcripts at low pH. A hallmark of cellular stress is up-regulation of heat shock proteins that are involved in many important cellular processes [49]. Specifically, Hsp70 expression has been reported to be induced in response to various environmental and physiological stressors [50]. Furthermore, we observed increased expression of glucose-regulated proteins (GRPs; GRP78 and GRP94) at pH 7.6, indicative of endoplasmic reticulum (ER) stress and unfolded protein response (UPR). GRPs assist in folding of damaged proteins, suppress caspase activation, and have positive effects on proliferation [51]. At 0.6 pH units decrease (pH 7.6) we observed a significantly higher expression of Nrf2 and NF-κB in A. viridis, two key transcription factors that regulate cellular responses to stressful conditions and promote cell survival [52, 53]. In stressful conditions, incorrectly folded proteins accumulate and aggregate, stimulating inflammatory responses and cellular responses to oxidative injury [54]. Oxidative stress plays apparently a critical role in the cnidarian bleaching cascade by modulating cell death and survival pathways [55, 56]. Respiration, photosynthesis, and autotrophy have previously been reported to increase at low seawater pH in A. viridis (and its symbiont) at the same sampling site [3, 18]. Rather interestingly, respiration and photosynthetic rates were found unaltered during laboratory-induced short-term ocean acidification conditions in A. viridis [16]. Surprisingly, we did not observe any significant changes in gene expression of antioxidant enzymes catalase or glutathione S-transferase. However, these observations are in agreement with other studies of symbiotic sea anemones under stressful conditions [33, 44, 57].

Besides observed cellular stress, several metabolism-related enzymes were deregulated in A. viridis. We observed an up-regulation of heme oxygenase 2 (HO-2), implicated in protection against lipid peroxidation [58], in A. viridis at pH 7.6. Importantly, one of the Niemann-Pick type c proteins regulating intracellular sterol-trafficking (NPC2d) was observed down-regulated at pH 7.6 in A. viridis. NPC2d has been proposed as a key regulator of host-symbiont relationship maintenance [59], localized exclusively in vacuole structures containing the symbionts [59]. The NPC2d gene from A. viridis, Exaiptasia pulchella and Exaiptasia pallida have all been observed up-regulated in the symbiotic state, and its expression was negatively affected in sea anemones undergoing various stresses including bleaching [39, 57, 59]. In addition to decrease in NPC2d expression, our data also indicate that the host-symbiont relationship appears affected at high seawater CO2. However, we observed no physical sign of sea anemone bleaching at the same sampling sites.

One of the most significant findings in our study was the dramatic up-regulation of LTR retrotransposons at high seawater CO2 (pH 7.6). Activation of retrotransposons, which constitute one of two main classes of transposable genetic elements, may impact the genome structure and genome function [60]. These include genome expansion and chromosomal rearrangements, as well as gene expression abnormalities. Transposable elements appear to play an important biological role in response to biotic and abiotic stresses [27], and an increased activity has been reported during bleaching of the reef building coral Montastraea faveolata [32]. Furthermore, substantial activation of LTR retrotransposons has been observed in the marine diatom Phaeodactylum tricornutum upon nitrate starvation [61]. But what could be the biological role of LTR retrotransposon activation in A. viridis at high seawater CO2 conditions? One possibility is that the activation can be a pure selfish response of the transposable element in order to escape from the host environment during stress. This resembles stress-induced proliferations of certain viruses previously reported in corals [62]. Alternatively, the activation could induce genome rearrangements that are of advantage to the host as a response to environmental changes and could lead to local adaptation to the low pH environment. We cannot exclude that some of the observed acclimatization processes involved in response to high CO2 environment could be influenced by local population genetic structure.

Conclusion

We confirm that the sea anemone Anemonia viridis and its algal symbiont Symbiodinum sp. are highly resilient to ocean acidification and that these organisms tolerate the changes in carbonate chemistry that occur as pH falls from 8.2 to 7.6. Despite no bleaching of the host, or any other obvious outward signs of stress, at the transcriptome level we observed the following responses in the sea anemone to chronic hypercapnia: 1) up-regulation of protective cellular functions; 2) deregulation of metabolic pathways; and 3) activation of transposable elements. Our data show that high CO2 produced a stress reaction in A. viridis, whereas Symbiodinium sp. was less affected by increased levels of CO2 at the transcriptome level. Studies at natural analogues for the effects of ocean acidification have shown that many organisms are unable to cope with long-term elevations in CO2, but those that do can be abundant. The present study serves as a reminder that so-called ‘winners’ in a high CO2 world may incur significant costs.

Materials and methods

Sampling

The temperate sea anemone Anemonia viridis (the Snakelocks Anemone) was collected in Levante Bay, Vulcano Island, Sicily - Italy (Fig 1). Here, low seawater pH conditions are created by CO2 release from a natural vent site at -1 m depth [63, 64]. Sampling was performed on May 13 and 14, 2013, at a depth of 1–2 meters from two different locations >350 m from the vent site along a gradient of decreasing pH (~ pH 7.6 and pH 7.9), and at a reference site ~800 m from the vent site with pH corresponding to ambient seawater levels (~ pH 8.2). For simplicity, we are referring to average pH values throughout this work as reported in Johnson et al. [63], but we know the sea anemones were exposed to variable pH conditions due to shallow currents influenced mostly by dominant western winds, as previously documented [63, 64]. Therefore, we also independently collected measurements of temperature, salinity, pH (NBS scale), and water samples for total alkalinity, at each sampling site to note variability between different days of sampling (S13 Table). In addition, sampling days and times were noted for all the individual samples taken (S14 Table). The sampling site is not only well described in regards to various geochemical measurements (e.g. [63, 64]), several physiological measurements of A. viridis have been previously performed at this site [3, 18, 20, 65]. These show that long-term exposure to high pCO2 has no impact on Symbiodinium density (cells mg protein-1), as well as on protein concentrations between sea anemones from the different sampling locations [18, 20, 65]. In addition, the clade of Symbiodinium has been previously determined for this sampling site [3]. Symbiodinum populations across all sites/ anemones were identified to be the same ITS2 ‘type’ of clade A (A19). A total of twelve A. viridis specimens were sampled in two consecutive days (four from each location). Small pieces of tissue (approximately 0.5 cm3) were collected from body wall, tentacles, and oral disc of each specimen. Each tissue type from each individual was subsequently stored separately at 4°C in RNAlater (ThermoFisher Scientific, Waltham, MA, USA), and transported from the sampling site to the laboratory. RNAlater solution was removed and all samples were stored at -80°C prior to further analysis.

RNA extraction

One ml of ice cold TRIzol reagent (ThermoFisher Scientific, Waltham, MA, USA) was added to the frozen tissue sample after removing excess RNAlater, and then quickly crushed using Precellys lysis homogenizer at 6000 rpm for 30 seconds (Stretton Scientific, Stretton, UK) in order to minimize degradation of RNA. RNA was twice extracted by chloroform, subsequently precipitated in isopropanol at -20°C overnight, centrifuged at 4°C, washed with 70% ethanol, and rehydrated in Nuclease-Free Water (ThermoFisher Scientific, Waltham, MA, USA). The RNA quality was examined and RNA amount quantified using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and Qubit 2.0 fluorometer (ThermoFisher Scientific, Waltham, MA, USA), respectively. Total RNA was isolated from each tissue type of an individual separately. Only high quality samples with an RNA integrity number (RIN) value equal to 7 or higher were used in further sequencing experiments.

Whole transcriptome sequencing

Whole transcriptome sequencing of A. viridis fragment libraries was performed using our in-house Ion Torrent Personal Genome Machine (PGM) sequencing platform (UiT, Norway). Equal amounts of total RNA from three different tissue samples of each specimen were pooled to prepare the sequencing libraries. PolyA enrichment was performed using Dynabeads mRNA DIRECT Purification Kit (ThermoFisher Scientific, Waltham, MA, USA). A total of twelve transcriptome libraries were prepared corresponding to four biological replicates at each location using Ion Total RNA-Seq Kit v2 (ThermoFisher Scientific, Waltham, MA, USA). Here, two replicates were prepared from each sampling day. Libraries were barcoded using Ion Xpress RNA-Seq Barcode 01–16 Kit (ThermoFisher Scientific, Waltham, MA, USA) and were sequenced once or twice on six different PGM 318 chips using the Ion PGM 200 Sequencing Kits (ThermoFisher Scientific, Waltham, MA, USA).

Transcriptome assembly and analysis pipeline

Raw Ion Torrent PGM sequencing reads from the same specimens were pooled and trimmed using cutadapt [66]. Barcodes and adapters were removed and reads were quality filtered using Phred quality score > 20. Quality of reads was monitored before, during and after trimming and filtering using the FastQC program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Based on quality graphs from FastQC program, the first 10 nucleotides of each sequence read were removed, and all sequences less than 30 nucleotides were discarded. De novo assembly was performed by the Trinity assembler software v2.0.6 using sequences from all individual samples [36]. To reduce contig redundancy, the assembled contigs were further processed by CD-HIT-EST clustering tool using 90% sequence similarity [67]. The assembly was explored for putative contamination by BLASTing against univec and refseq for bacteria and fungi (BLASTN, e-value < 10-5), where 352 contigs with higher similarity to these databases than to nucleotide databases of Cnidaria and Dinoflagellata, were filtered away from the assembly. This assembly was used as a reference transcriptome for further analyses. Quality filtered sequence reads from each specimen were aligned to the reference transcriptome using BWA mapping software v0.7.12 (BWA-MEM) [68], and an estimation of relative expression levels of transcripts was performed using eXpress v1.5.1 [69]. The assembled reference transcriptome was divided into host and symbiont fraction using BLASTX algorithm for training data sets and Support Vector Machine (SVM) classification for separating the two fractions; both are implemented in the PSyTranS tool (https://github.com/sylvainforet/psytrans). In this way, all contigs can be assigned to either of the fractions, not only those that have a BLAST match. A total of 27,273 published and predicted protein gene sequences from the sea anemone N. vectensis [31] and 47,014 predicted protein gene sequences from the symbiont Symbiodinium minutum [70], the closest relatives with published predicted protein sequences at the time of the analysis, were used as training data sets. Successful separation of the host and symbiont contigs was assessed by a GC-content plot, which displayed a clear bimodal distribution. We performed additional analysis to further validate the performance of the PSyTranS tool (S1 Table). After separation of mapped sequence reads into host and symbiont fraction, host and symbiont were assessed individually for statistically significant DE-transcripts by the edgeR package [71]. Raw read counts from eXpress were filtered for lowly expressed transcripts (with minimum 10 reads in at least four samples). Data were normalized using TMM (edgeR) [71] and differential gene expression was inferred using a two-factor, negative binomial generalized linear model (glm edgeR) [40], accounting both for the different sampling sites (pH 7.6, 7.9 and 8.2) and various days of sampling (May 13 and 14, 2013). Furthermore, DE-transcripts were filtered with two-fold change difference and a FDR cutoff of 0.05 (S3 Table). DE-transcripts were then hierarchically clustered into heatmap based on normalized expression values (rpkm) and scaled by row.

Annotation pipeline

The transcriptome assembly was assessed for homologous protein-coding genes by BLAST (Basic Local Alignment Search Tool) [72]. Contigs were blasted against locally installed protein databases (nr and Swiss-Prot/UniProtKB) on a high-performance cluster using BLASTX (e-value < 10-5). Functional GO annotation of our reference transcriptome was carried out using a local Blast2GO pipeline B2G4Pipe [37]. GO analysis was performed using GOseq R package [42] to further explore over- and underrepresented GO categories among different pH conditions. Transcripts with a likely protein-coding capacity were extracted from reference transcriptome by Transdecoder, included in the Trinity package [73], and translated queries were further searched locally for homologous sequences using the hidden Markov model (HMM) [74] against Pfam database [75] (hmmscan, e-value < 10-3). InterPro domain search was performed on the Transdecoder output using RunIprScan (http://michaelrthon.com/runiprscan/). Sequences that lacked hits in protein databases were searched for ORF features using OrfPredictor (http://bioinformatics.ysu.edu/tools/OrfPredictor.html) [76].

Identification of transposable elements in the reference transcriptome

Searches for transposable elements in the complete assembly (154,015 contigs) were performed using RepBase database v20.10 [41], which was queried against our reference transcriptome using TBLASTX (e-value < 10-20). A tabular output from the search was further run through two perl scripts, “blast92gff3.pl” with additional options -lowscore 0.0001 -alignmax 9999 -exonType exon (http://arthropods.eugenes.org/EvidentialGene/evigene/scripts/blast92gff3.pl.) and the “overbestgene2.pl” (http://iubio.bio.indiana.edu/gmod/tandy/perls/) to create a gff file from blast results and to remove overlapping blast hits, respectively. Only sequence regions corresponding to transposable elements in our transcriptome reference assembly were then extracted from our contigs using BLAST fastacmd tool and substituted the corresponding contigs in the reference transcriptome. Subsequently, an analogous glm edgeR DE analysis was performed, where expression changes of transposable elements were investigated.

Quantitative PCR

Selected transcripts reported by glm edgeR analysis as DE between conditions were further verified by quantitative PCR (qPCR). The ribosomal protein L12 (RPL12), beta-actin and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) were used as endogenous controls for normalization. 20 bp primers with low self-complementarity were designed by Primer 3 [77] to produce 90 – 110 bp products (S2 Table), and amplicons were tested by regular PCR and gel electrophoresis prior to qPCR. PCR products were sequenced by Sanger sequencing to confirm the amplified sequences. RNA for the qPCR analysis was isolated from the same samples that were used for preparation of transcriptome libraries. An equal amount of RNA from each specimen (250 ng) was DNase I treated, and cDNA was subsequently prepared using SuperScript III Reverse Transcriptase (ThermoFisher Scientific, Waltham, MA, USA) with combination of random hexamer primers. cDNA was 10 x diluted and 2 μl was used together with 2x SYBR mix and 2.5 μl of 1 μM primer mix, in 10 μl reactions as input into Roche Light Cycler 96 (Roche, Basel, Switzerland). Each sample was analysed in duplicate by qPCR. Mean Ct (cycle threshold) values for samples were normalized to Ct values of the reference genes (RPL12, beta-actin and GAPDH) and plotted as ΔΔCt values compared to the control condition (pH 8.2).

Gene expression profiles of Anemonia viridis and Symbiodinium sp. at low pH.

Shown are differential gene expression profiles of (A) the host and (B) the symbiont from two sampling locations of decreasing pH (pH 7.6 and pH 7.9) compared to normal seawater pH 8.2. The differential expression values are plotted as logarithmic fold change (logFC) values, and were calculated from the average expression values of four individuals per sampling location. The heatmap serves only as a visualization of transcripts that are significantly up- and down-regulated at the two low pH conditions compared to normal seawater pH 8.2.

(PDF)

Click here for additional data file.

Verification of significantly differentially expressed genes by quantitative PCR.

We assessed the expression of nine transcripts reported as differentially expressed in our RNA-seq analysis by quantitative (q)PCR. Three individuals per each condition studied (pH 7.6, pH 7.9 and pH 8.2) were examined, and values were normalized to the expression of our reference transcripts: ribosomal protein L12 (RPL12), beta-actin and glyceraldehyde 3-phosphate (GAPDH). Gene expression at two low seawater pH sites (pH 7.6 and pH 7.9) is shown as a relative measure with standard deviations compared to the gene expression at normal seawater pH 8.2.

(PDF)

Click here for additional data file.

Expression levels of six candidate reference genes at the sampling sites.

The expression levels of six potential reference transcripts: glyceraldehyde 3-phosphate dehydrogenase (GAPDH), ribosomal protein L12 (RPL12), beta-actin, adenosylhomocysteinase (AHCY), senescence-associated protein and NADH dehydrogenase (NDH), were assessed among the different sampling sites. Presented are the normalized values of these transcripts from the RNA-seq transcript expression matrix. Three transcripts were successfully established as reference transcripts after qPCR testing: GAPDH, RPL12 and beta-actin.

(PDF)

Click here for additional data file.

Sequencing information and various metrics of the assembled reference transcriptome of Anemonia viridis.

(PDF)

Click here for additional data file.

Primer pairs used for regular and quantitative PCR analysis.

Shown are primer sequences used for PCR amplification of symbiont-specific nuclear apx gene, together with 12 primer pairs used for verification of results from differential expression analysis by qPCR.

(PDF)

Click here for additional data file.

R script for differential expression analyses.

Shown is the R script of glm edgeR analyses that was used to find differentially expressed genes between the various pH conditions separately for the host and for the symbiont, using four different individuals per condition. Additionally, PCA plots and heatmaps were created.

(PDF)

Click here for additional data file.

R script to show the influence of each variable in glm edgeR analyses.

Shown is the R script used to investigate in detail the influence of each variable on the differential gene expression analyses separately for the host and for the symbiont.

Differentially expressed genes identified in our glm edgeR analysis were separated based on the influence of the two variables in our study, pH condition and day of sampling. By forming contrasts from the design matrix in the glm edgeR pipeline, we were able to identify 1124 DE-transcripts for the host Anemonia viridis (91.2%) and 249 DE-transcripts for the symbiont Symbiodinium sp. (85.6%) affected by the pH condition variable only. We therefore concluded that day of sampling did not have significant effect on the identification of the differentially expressed genes between the individual pH conditions.

(PDF)

Click here for additional data file.

R script for gene set enrichment analyses.

Shown is the R script used to find enriched gene ontology (GO) categories at pH 7.6 compared to normal seawater pH 8.2.

(PDF)

Click here for additional data file.

Enriched Gene Ontology (GO) categories at pH 7.6.

A list of up- and down-regulated gene ontology (GO) categories in Anemonia viridis detected at low pH 7.6 compared to normal seawater pH 8.2.

(PDF)

Click here for additional data file.

Selected differentially expressed transcripts at low pH compared to normal seawater pH in A. viridis.

Shown is an extended list of significantly differentially expressed transcripts at low seawater pH 7.6 compared to normal seawater pH 8.2 in A. viridis from glm edgeR analysis (FDR < 0.05).

(PDF)

Click here for additional data file.

Selected differentially expressed transcripts at low pH compared to normal seawater pH in Symbiodinium sp.

Shown is a list of selected differentially expressed transcripts at (A) pH 7.6 and (B) pH 7.9 compared to normal seawater pH 8.2 in Symbiodinium sp. from glm edgeR analysis (FDR < 0.05).

(PDF)

Click here for additional data file.

Analysis of protein signatures in the symbiont DE-data set.

Protein-coding regions within differentially expressed transcripts were searched for domains or other functional signatures using InterPro database. Depicted is a list of recognized signatures that appeared up- or down-regulation in the symbiont at low seawater pH 7.6 compared to normal seawater pH 8.2.

(PDF)

Click here for additional data file.

Presence of symbiont stress-response genes in the reference transcriptome assembly.

The presence of important Symbiodinium sp. stress-response genes in our reference transcriptome assembly have been assessed. Eight different stress-response proteins, corresponding to 62 different proteins from the genus Symbiodinium present in the NCBI database, have been queried to the symbiont reference transcriptome assembly (e-value < 10-3). 29 symbiont contig hits were then BLASTed to the NCBI’s nr database and the best alignments are presented below.

The following stress-response genes have been assessed: heat shock proteins 70 and 90 (Hsp70 and Hsp90), superoxide dismutases (SODs), glutathione reductase (GR), thioredoxin (TRX), catalase peroxidase (katG), ascorbate peroxidase (APX) and cytochrome P450 (CYP450). Neither of these proteins was observed as DE in our study. We only observed up-regulation of certain InterPro domains at pH 7.6 compared to normal seawater pH 8.2 (S9 Table).

(PDF)

Click here for additional data file.

Differentially expressed SYMBIOTIC transcripts at low seawater pH 7.6 compared to normal seawater pH 8.2 in Anemonia viridis.

(PDF)

Click here for additional data file.

Differentially expressed APOSYMBIOTIC transcripts at low seawater pH 7.6 compared to normal seawater pH 8.2 in Anemonia viridis.

(PDF)

Click here for additional data file.

Independent measurements recorded at the sampling site.

Average (S.D.) values of carbonate chemistry parameters at sampling locations off Vulcano Island CO2 seeps. On 13th and 14th May 2013, we performed daily measurements (am and pm) of pHNBS, salinity and temperature (n = 3–4) at the study site during sea anemone sampling, using a 556 MPS YSI (Yellow Springs, USA) probe. The pH sensor was calibrated using NBS scale standards buffers. Three replicate sub-samples of seawater were analyzed at 25° C for total alkalinity (TA) using a titration system (Mettler Toledo, Inc.).

(PDF)

Click here for additional data file.

Dates and times of sampling of Anemonia viridis from the individual sampling locations.

Adult polyps of the sea anemone were sampled from the individual locations accessed from shore (sampling site pH 7.6 and pH 8.2) or from a boat (pH 7.9). The sampling times at the individual dates could not be synchronized due to bad weather forecast reported for the following days of the stay. Therefore, samples had to be taken as soon as possible on the 14th May. However, we adjusted for the day of sampling in our statistical model (glm edgeR) and found that only small amount of the transcripts were affected by the day of the sampling (see more information in the S4 Table).

(PDF)

Click here for additional data file.