Correction: 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 1.4% of the anemone transcripts, but only ~0.5% of the Symbiodinium sp. transcripts were differentially expressed. Processes enriched at high seawater CO2 were mainly linked to cellular stress, including significant up-regulation of protective cellular functions and deregulation 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 CO 2 [1], yet non-calcified cnidarians appear to be more resilient [2][3][4][5][6]. Since  expression profiling analysis presented below. We infer that the host-symbiont relationship is not disrupted at low seawater pH conditions.

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).
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 fourand 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.
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 ). 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. Ocean acidification effects on a symbiotic sea anemone

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.
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 upregulation 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 upregulation 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 hostsymbiotic 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 CO 2 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 CO 2 in the sea anemone, but not in the dinoflagellate. Even though high CO 2 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 CO 2 , 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 CO 2 . 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 CO 2 (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 CO 2 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 CO 2 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 CO 2 produced a stress reaction in A. viridis, whereas Symbiodinium sp. was less affected by increased levels of CO 2 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 CO 2 , but those that do can be abundant. The present study serves as a reminder that so-called 'winners' in a high CO 2 world may incur significant costs.

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 CO 2 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 pCO 2 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 cm 3 ) 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 inhouse 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 (evalue < 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). 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.  Table. 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) S3 Table. 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) S4 Table. 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 DEtranscripts 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)  Table. 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) S8 Table. 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) S9 Table. 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) S10 Table. 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) S11 Table. Differentially expressed SYMBIOTIC transcripts at low seawater pH 7.6 compared to normal seawater pH 8.2 in Anemonia viridis. (PDF) S12 Table. Differentially expressed APOSYMBIOTIC transcripts at low seawater pH 7.6 compared to normal seawater pH 8.2 in Anemonia viridis. (PDF) S13 Table. Independent measurements recorded at the sampling site. Average (S.D.) values of carbonate chemistry parameters at sampling locations off Vulcano Island CO 2 seeps. On 13 th and 14 th May 2013, we performed daily measurements (am and pm) of pH NBS , 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) S14 Table. 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 14 th 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)