Genes Related to Ion-Transport and Energy Production Are Upregulated in Response to CO2-Driven pH Decrease in Corals: New Insights from Transcriptome Analysis

Since the preindustrial era, the average surface ocean pH has declined by 0.1 pH units and is predicted to decline by an additional 0.3 units by the year 2100. Although subtle, this decreasing pH has profound effects on the seawater saturation state of carbonate minerals and is thus predicted to impact on calcifying organisms. Among these are the scleractinian corals, which are the main builders of tropical coral reefs. Several recent studies have evaluated the physiological impact of low pH, particularly in relation to coral growth and calcification. However, very few studies have focused on the impact of low pH at the global molecular level. In this context we investigated global transcriptomic modifications in a scleractinian coral (Pocillopora damicornis) exposed to pH 7.4 compared to pH 8.1during a 3-week period. The RNAseq approach shows that 16% of our transcriptome was affected by the treatment with 6% of upregulations and 10% of downregulations. A more detailed analysis suggests that the downregulations are less coordinated than the upregulations and allowed the identification of several biological functions of interest. In order to better understand the links between these functions and the pH, transcript abundance of 48 candidate genes was quantified by q-RT-PCR (corals exposed at pH 7.2 and 7.8 for 3 weeks). The combined results of these two approaches suggest that pH≥7.4 induces an upregulation of genes coding for proteins involved in calcium and carbonate transport, conversion of CO2 into HCO3 − and organic matrix that may sustain calcification. Concomitantly, genes coding for heterotrophic and autotrophic related proteins are upregulated. This can reflect that low pH may increase the coral energy requirements, leading to an increase of energetic metabolism with the mobilization of energy reserves. In addition, the uncoordinated downregulations measured can reflect a general trade-off mechanism that may enable energy reallocation.


Introduction
Coral reefs are key ecosystems characterized by a high level of biodiversity, ecological complexity and primary productivity. Reef-building corals are key species in coral reefs, providing physical and biological foundations for the ecosystem. The ability of corals to build reefs in oligotrophic tropical oceans is mainly explained by the mutualistic association they form with zooxanthellae (photosynthetic dinoflagellates of the genus Symbiodinium), which produce autotrophically most of the energy required for the calcification process [1].
Global climate change is impacting on coral reefs, with approximately 19% of reefs worldwide being permanently degraded, 15% showing symptoms of imminent collapse, and another 20% at risk of becoming critically affected in the next few decades [2]. This alarming level of reef degradation is mainly due to an increase in frequency and intensity of natural and anthropogenic disturbances affecting reefs [3]. Among, these disturbances is ocean pH decrease, also termed ''ocean acidification'' [4,5].
Several experimental studies have shown that low pH decreases the net calcification rate and skeletal growth of various coral species [6][7][8][9][10][11][12][13][14][15]. However, in other studies no effects on gross calcification (or an increase) have been reported [16][17][18]. In ecosystems naturally subject to low pH conditions (e.g. CO 2 vents), coral net calcification has been reported to not decrease, or to decrease only slightly [19,20], but at the ecosystem level low pH induces marked community changes [19]. In addition to its impact on calcification rates, low pH has been shown to: (i) modify the composition, size, shape and orientation of aragonite crystals in primary polyps [7]; and (ii) increase the coral tissue biomass [11,21]. Contradictory findings have been reported in relation to the effects of low pH on the symbiont component of the coral holobiont. These include findings of a decrease in the density of zooxanthellae [17,22], and effects on their chlorophyll concentration [11] and photosynthetic efficiency [9]. On the other hand, these results have been challenged in recent studies [9,14,16]. To date only 2 recent works have explored the transcriptomic response of corals exposed to low pH: a whole transcriptomic analysis on early life history stages [23] and a cDNA microarrays analysis on adults [17]. Both of these studies were performed on the coral Acropora millepora which phylogenetically belongs to the coral ''complex'' clade which has a less heavily calcified skeletons, and a more porous construction of corallite walls than corals belonging to the ''robust'' clade [24,25].
The aim of the present study was to initiate a wholetranscriptome analysis in the adult scleractinian tropical coral Pocillopora damicornis exposed to low pH. Such whole-transcriptome analysis has already successfully been used in studies on the impact of global climate change on coastal organisms [23,[26][27][28]. It is also an interesting prerequisite for studies that will target specific biological functions or phenotypic plasticity across populations facing different pH conditions. Finally, RNAseq experiments offer the possibility to highlight unexpected physiological pathways. In the present study we choose to work on the scleractinian tropical coral P. damicornis for several reasons: i) it is a major reef-building coral widely distributed in the Indo-Pacific ocean [29]; ii) it is naturally confronted to diverse values of ocean pH, ranging from 7.8 in the particular case of volcanic CO 2 vents, to 7.9-8.0 in the Eastern tropical Pacific and to 8.1 in other locations [19,30]; iii) it is known to be highly sensitive to a wide range of natural and anthropogenic disturbances [31][32][33]; iv) it belongs to the ''robust'' clade which has not yet been studied by transcriptomic approaches [24,25]. This ''robust'' clade is characterized by a highly calcified and low porous skeleton compared to corals belonging to the ''complex'' clade. Taking into account the impact of low pH on calcifying organisms, these last characteristics can lead to different responses to low pH conditions between clades. To address our aims we first assembled the transcriptomes of P. damicornis using the result of several RNAseq experiments. In order to maximize the amplitude of the transcriptomic response, an exposure of three weeks to an extreme value of today (volcanic CO 2 vents) or near future natural variability (lower value on the daily variation scale) of pH (7.4) was used and the transcriptome obtained was compared to the control transcriptome (nubbins maintained for three weeks at a current seawater pH of 8.1). After the identification of the biological functions of interest responding to this low pH, a candidate gene approach (by q-RT-PCR) was performed on mRNA from nubbins exposed to a predicted pH of 7.8 and to a very low pH of 7.2 (3 weeks of exposure).

Coral collection, maintenance, and experimental stress
The Pocillopora damicornis (Linnaeus, 1758) isolate used in this study was collected in Lombok, Indonesia (CITES Management Authority, CITES number 06832/VI/SATS/LN/2001-E; Direction de l'Environnement, CITES number 06832/VI/SATS/LN/ 2001-I). The site of collection corresponds to an usual fringing-reef subjected to the usual day/night pH variation (8.3/7.85) observed on such reefs [34,35]. P. damicornis colonies were then maintained since the year 2001 in an open water aquarium system where pH fluctuates between day and night due to photosynthesis and respiration between maxima and minima values of 8.3 and 7.9 respectively for day and night. For experiments, coral explants (3 cm height, 3 cm diameter) from the same parent colony were left to allow tissues to grow over skeleton (recover) for a period of 1 month prior to use in experiments.
Following recovery, coral nubbins were incubated during 3 weeks into experimental tanks. Temperature (2560.5uC) was kept constant in each tank using heaters connected to electronic controllers (IKS, Karlsbad, Germany), and seawater was mixed using 2 submersible pumps. Corals received a constant saturating irradiance of 175610 mmol photons m 22 s 21 (photoperiod was 12 h:12 h light/dark) using HQI-10000K metal halide lamps (BLV-Nepturion). They were fed twice a week with Artemia salina nauplii in order to reach a natural heterotrophic rate. Tanks, pumps, sensors and electrodes were rigorously cleaned every week to prevent the growth of epiphytic algae and fouling communities or the accumulation of detritus. Water renewal rate in each tank was 60% per hour (continuous water flow).
Carbonate chemistry was manipulated in 3 of 4 tanks by bubbling with CO 2 to reduce pH to target values and the fourth tank (control tank) was not bubbled with CO 2 (see values below). pH electrodes and temperature sensors installed in each tank were connected to a pH-stat system (IKS, Karlsbad) that continuously monitored pH (calibrated to NBS scale), temperature and controlled CO 2 bubbling. On a daily basis, additional pH checks were carried with a Plastogel-Ponsel pH probe (Ponsel), calibrated to pH total scale (pH T ). Moreover biweekly pH T measurements were made using the indicator dye m-cresol purple (mCP Acros 199250050) adapted from Dickson et al. [36]; the absorbance was measured using spectrophotometer (UVmc 2 , Safas, Monaco). According to these results, the pH-stat system was adjusted. The control tank was maintained at ambient pH T : 8.0760.07 (396 matm pCO 2 ), the other tanks at pH T : 7.7560.12, (856 matm pCO 2 ) ; pH T : 7.4260.21 (2181 matm pCO 2 ) ; pH T : 7.1960.18 (3880 matm pCO 2 ). For easier reading, throughout the manuscript, these four values of pH T will be respectively indicated as pH 8.1 (control tank), pH 7.8, pH 7.4 and pH 7.2. After 3 weeks of exposure, 3 nubbins in each treatment tank were sampled at the same time and stored in liquid nitrogen until analyzed.
Water samples were collected in scintillation vials (Wheaton) and stored refrigerated no more than 1 week prior to measurement. Biweekly total alkalinity (TA) was measured via titration with 0.01 N HCl containing 40.7 g NaCl l 21 using a Metrohm Titrando 888 dosimat controlled by Tiamo software to perform automated normalized Gran titrations of 1 ml samples. For each sample run, certified seawater reference material supplied by the lab of A.G. Dickson (Scripps Institution of Oceanography) was used to check quality of measurements.
Parameters of carbonate seawater chemistry were calculated from pH T , mean AT, temperature, and salinity using the free access CO 2 Systat package [37] using constants from Mehrbach et al. [38] as refit by Dickson and Millero [39]. Parameters of carbonate seawater chemistry are given in Table 1 and 2.

cDNA library construction and high-throughput sequencing
We assembled a transcriptome de novo from 80-nucleotide paired-end short sequence reads (SSRs) based on 6 lanes of Illumina sequencing. To maximize the coverage of potential transcriptomes, the cDNA from corals exposed to various environmental conditions (pH 7.4 and its dedicated control (see above), thermal stress, bacterial challenge, bacterial infection and a dedicated control for these last conditions [40]) were sequenced (one illumina lane per condition). Despite the recent publication of a reference transcriptome for P. damicornis [41] we preferred to assemble de novo a reference transcriptome with the SSRs generated by our RNAseq approach rather than map our reads on the one previously published. Indeed, this reference transcriptome (Roche 454 data) was obtained from Hawaiian P. damicornis, a population that is geographically very distant from the one we studied (Lombok, Indonesia). Such a distance may lead to a high number of single nucleotide polymorphisms that can skew the mapping step. In addition, the Hawaiian study did not use specimens exposed to low pH, thatpH that could lead to the absence of transcripts or spliced transcripts specifically expressed under such environmental condition.
Total RNA from 3 nubbins per condition was extracted, purified using TRIzol reagent (Invitrogen), and pooled together (one pool per condition) to reach a sufficient amount of biological material for the RNAseq procedure, as described previously [42]. A cDNA library was subsequently constructed following a published protocol [43]. A minimum of 1 mg of cDNA was used to generate each of the paired-end Illumina sequencing libraries. The libraries were prepared using Illumina adapter and PCR primers according to previously published protocols [44,45]. Libraries with an average insert size of 400-500 bp were isolated and the concentration was adjusted to 10 nM. The samples (7 pM per sample) were loaded into separate channels of an Illumina GAIIx sequencer, and sequenced at the Oregon State University Center for Genome Research and Biocomputing.

Reference transcriptome assembly and annotation
The SSRs obtained were processed using an Illumina pipeline (RTA1.6_CASAVA1.6), and SSRs from 'sequence.txt' files were filtered to select only those SSRs that matched high quality standards. Velvet software version 1.1.04 was subsequently used, evoking velveth with parameters [k-mer] -fasta -shortPaired, and velvetg with parameter -read_trkg yes [46]. The k-mer value was tested between 31 and 61 and optimized using a subset of data for maximum size of longest contig, average contig length, contig number, n50 value and total transcriptome length. Optima were obtained at k-mer 55 and 57, and we decided for k-mer 55 based on the higher number of included reads. Using the above parameters, 45,596,723 of 131,859,606 SSRs were incorporated into the assembly, producing 102,571 contigs. These served as the input for the de novo transcriptome assembler software Oases 0.1.21 [47], using the input parameter -ins_length 400 and default parameters -cov_cutoff 3 -min_pair_cov 4 -paired_cutoff 0.1 and with or without -scaffolding. This assembly generated 84,511 contigs. The scaffolding had no effect. Since oases 0.1.x was known to generate a number of exactly identical sequences (duplicates), a reassembly was performed using TIGR assembler [48] based on the default parameters in the original run_TA script (ftp://ftp.jcvi.org/pub/ software/assembler/), which reduced the number of contigs to 72,924. The remaining exact duplicates and reverse complement exact duplicates were removed using prinseq-lite.pl (http://prinseq. sourceforge.net; parameter -derep 14), yielding a total of 72,890 transcript contigs.
Functional annotation of this reference transcriptome was performed using Blast2GO version 2.4.2 [49], which enabled a semi-automated functional annotation of all contigs using a set of similarity search tools. The default parameters was used and included: i) an initial annotation with BLASTX (against the nonredundant NCBI database; e-value at 1610 -3 ); ii) a protein domain search using InterProscan; iii) an enzyme annotation using the Kyoto Encyclopedia of Genes and Genomes (KEGG); and iv) assignment of a Gene Ontology term (GO; http://www. geneontology.org/).
The taxonomic assignment of contigs (sequences from the coral host or zooxanthellae symbiont) were predicted using 3 approach-es: i) top hit species from the BLASTX results (the species with the best/first sequence alignment for a given BLASTX result); ii) the percentage of GC bases, GC content of corals is approximately 42%, and 54% for zooxanthellae [50]; and iii) species specific sequence signatures based on CLaMS analysis [51].

Differential gene expression analysis
To assess changes in gene expression regulated by low pH, SSRs from both low pH and control libraries were mapped against the reference transcriptome and counted using BWA [52]. Because RNAseq data are functions of both the molar concentration and the transcript length, the results of the mapping step were corrected and expressed as reads per kilobase per million mapped reads (RPKM; [53]). This enabled comparison of SSRs counts between transcripts and from samples derived from differing conditions [53]. To identify significantly different gene expression among conditions, the MARS method (MA plot-based methods using a random sampling model) of the R package in DEGseq was used [54]. Differences in transcripts levels between low pH and control conditions were considered statistically significant at p,0.0001.

q-RT-PCR
Quantitative real-time PCR (q-RT-PCR) was used to validate expression profiles obtained from the MARS DEGseq analysis of the RNAseq data [27]. It was also used to measure expression levels of selected 48 candidate genes in P. damicornis subjected to complementary pH conditions (7.8 and 7.2). Total RNA was extracted, treated with DNase, and the poly(A) RNA was purified as described above. Approximately 50 ng of purified poly(A) RNA were reverse transcribed with hexamer random primers using ReverTAid H Minus Reverse Transcriptase (Fermentas). The q-RT-PCR experiments were performed on cDNAs obtained from 3 nubbins per treatments as described previously [40]. For each candidate gene the level of transcription was normalized using the mean geometric transcription rate of 3 reference sequences encoding ribosomal protein genes from P. damicornis (60S ribosomal protein L22, GenBank accession numbers HO112261, 60S ribosomal protein L40A, HO112283 and 60S acidic ribosomal phosphoprotein P0, HO112666). These housekeeping genes had ,10% variation in expression between control and low pH conditions and were previously shown to be stable in P. damicornis subjected to different environmental conditions [40]. The primers used for amplification are provided in File S1.

Assembly of the reference transcriptome
To quantify changes in gene expression between the pH 7.4 treatment and the control, we assembled a reference transcriptome de novo. To maximize the diversity of transcripts, the cDNA samples derived from corals exposed to various environmental conditions (pH 7.4 and its dedicated control, thermal stress, bacterial challenge, bacterial infection and a dedicated control for these last conditions) were sequenced and assembled together. After the assemble work flow, a total of 72,890 transcript contigs were generated. The functional annotation performed using Blast2GO produced the following result: i) 22,419 contigs had significant similarities to proteins of known function; ii) GO term was assigned to 23,170 contigs; iii) KEGG enzymatic codes were found for 2,837 contigs; and vi) 46,295 conserved protein domains were detected by InterProScan. Using a combination of BlastX top hit species results, GC% calculation of each contig (perl script) and sequence signature analysis (CLaMS program; [51], 27.7% and 69.8% of the contigs were predicted to belong to the symbiont and the host transcriptome, respectively. The remaining contigs could not be taxonomically attributed. The assembled and annotated reference transcriptome is publicly available at http:// 2ei.univ-perp.fr/telechargement/transcriptomes/blast2go_fasta_Pda mv2.zip. The raw data (untreated SSRs) for the pH treatment and its corresponding control are publicly available at http://www.ncbi.nlm. nih.gov/sra (study accession number SRP011059.1).

Expression in 16% of the genes changes after exposure to low pH
The RNAseq approach resulted in the sequencing of 25.2 and 22.3 million Short Sequence Reads (SSRs) for the pH 8.1 and the pH 7.4 conditions, respectively. A total of approximately 10 and 9 million SSRs (,40%) passed the quality filter and were successfully mapped on the reference transcriptome for the control and the low pH conditions, respectively. Nearly 80% (57,259 contigs) of the reference transcripts were mapped by at least one SSR from either 1 of the 2 experimental conditions.
To highlight the genes significantly up or downregulated by the pH 7.4 treatment, the mapped results (RPKM) were analyzed using the MARS method in the DEGseq R package. In response to pH 7.4, 6.0% (3,204) of the transcripts were significantly upregulated and 10.0% (5,758) were downregulated (p,0.0001; Figure 1). Among the upregulated genes, 60.4% showed no significant similarity to known proteins, 2.7% showed significant similarity to proteins of unknown function, and 36.9% showed significant similarity to proteins of known function (Significant similarity: evalue ,1610 -3 ). Among the downregulated genes, 69.5% of their products showed no significant similarity to known proteins, 3.5% showed significant similarity to proteins of unknown function, and 27.0% showed significant similarity to proteins of known function.
To assess the accuracy of the quantification analysis a validation step involving an alternative method of quantification (q-RT-PCR) was conducted on 25 transcripts arbitrarily selected along the gradient of expression (from highly upregulated to highly downregulated gene) obtained from the RNAseq data. There was a significant correlation (r 2 = 0.86; p,0.0001; Figure 2) between the log 2 fold change in expression (relative to the control) of the RNAseq and the q-RT-PCR data which confirmed the accuracy of the results obtained with the RNAseq data.
Specific biological processes are altered after exposure to pH 7.4 for three weeks To study the transcriptomic changes resulting from exposure to pH 7.4 for 3 weeks, a GO term enrichment analysis (Figure 3) was performed [55]. It showed that 32 biological process categories were enriched in the subset of sequences associated with significant upregulation (Fig. 3A). Several of the functions involved were of particular interest in the context of this study, including transport, photosynthesis, gluconeogenesis, generation of precursor metabolites and energy and pyruvate metabolic process. Analysis of the subset of transcripts of significantly downregulated genes showed that only those for the GO term 'cellular macromolecule metabolic process' were represented at a higher rate than expected by chance (p,0.01; Fig. 3B). The GO terms 'post-translational protein modification' and 'cellular protein metabolism' were more significantly represented (p,0.05) among the downregulated set than in the reference sequence set. The numbers of genes belonging to each biological process are indicated in File S2.
Response of some biological functions to an exposure to pH 7.8, 7.4 or 7.2 for 3 weeks Based on GO term enrichment analysis and DEGseq results (Figure 3), 48 candidate genes significantly regulated (p,0.0001) were identified as corresponding to biological functions of interest in the context of our study and were selected (see File S1) and investigated further. Their transcript abundance was measured  Figure 4 and the numerical values are detailed in File S3. Briefly, the GO term enrichment analysis ( Figure 3) showed a significant increase in gene expression involved in the biological processes 'transport', which includes transport of ions and proteins. As calcification involves Ca 2+ and inorganic carbon [56] we investigated genes encoding proteins involved in the transcellular transport of Ca 2+ and HCO 3 2 ions, and carbonic anhydrases involved in converting CO 2 into HCO 3 2 . Among these candidates: i) The six HCO 3 2 transporters were upregulated at pH 7.8 and 7.4 (between 17.75 and 1.52 fold increase) and were downregulated at pH 7.2 (between 3.81 and 1.57 fold decrease); ii) a Ca 2+ plasma membrane ATPase was upregulated at pH 7.8 and 7.4 (1.47 and 1.43 fold increase, respectively), but was downregulated at pH 7.2 (3.01 fold decrease); iii) an extracellular (upregulated at all pH between 2.51 and 3.34 fold increase) and a cytosolic carbonic anhydrase (upregulated at pH 7.8 and 7.4, 2.51 and 1.85 fold increase respectively). The Ca 2+ plasma membrane ATPase and the two carbonic anhydrase present significant similarities with proteins that were previously shown to be involved in Stylophora pistillata calcification [57][58][59]. In addition to CaCO 3 , the skeleton of scleractinian corals contains organic matrix proteins [56], and their production was affected by low pH treatment. Indeed, seven genes encoding skeleton organic matrix proteins were identified in the present study. One of these genes regulated at pH 7.4 (3.5 fold increase) presents significant similarities to galaxin from Galaxea fascicularis [60], and the four other transcripts to the galaxin-like (only 2 of them were subjected to q-RT-PCR experiment, and were upregulated at all pH between 70.03 and 4.11 fold increase ) proteins family from Acropora millepora [61]. The remaining two genes had significant similarities to the bone morphogenetic protein superfamily [62] and were upregulated at pH 7.8 and 7.4 (between 2.91 and 1.93 fold increase), but downregulated at pH 7.2 (between 2.68 and 1.56 fold decrease). For a complete description of candidate genes expression see Figure 4 and File S3.
Biological functions involved in energy production (photosynthesis, glycolysis, Krebbs cycle, oxidative phosphorylation, lipolysis and beta oxidation; Figure 4) were also significantly enriched in the set of up-regulated genes (GO term enrichment analysis, Figure 3). As calcification represents approximately 13-30% of the daily coral energy demand [1], we also investigated the regulation of these metabolic related genes by q-RT-PCR for the 3 pH treatments. Among the eight candidates genes involved in photosynthesis, seven were upregulated (between 1.15 and 5.54 fold increase) and one, the PS I p700 chlorophyll a apoprotein a2like was downregulated at pH 7.4 (1.96 fold decrease). Three of the four candidate genes involved in glycolysis were upregulated at pH 7.8 and 7.4 (between 6.68 and 1.22 fold increase) and two were downregulated at pH 7.2 (between 2.11 and 1.22 fold decrease). Among the eight the candidates genes involved in the Krebs cycle six shows upregulations at pH 7.8 and/or 7.4 (between 3.36 and 1.33 fold increase) and five were downregulated at pH 7.2 (between 2.36 to 1.1 fold decrease). Four of the genes involved in oxidative phosphorylation were upregulated at pH 7.8 and 7.4 (between 10.85 and 1.28 fold increase), but one become downregulated at pH 7.2 (1.89 fold decrease). One of the genes involved in lipolysis was upregulated at pH 7.8 and 7.4 (1.15 and 1.95 fold increase respectively) but was downregulated by a factor of 2 in the pH 7.2 treatment. Finally, the three candidate genes involved in beta oxidation were upregulated at all pH (between 5.98 and 1.37 fold increase), except the Acyl-CoA dehydrogenaselike that was downregulated by a factor of 2.17 at pH 7.2. For a complete description of candidate genes expression see the Figure 4 and File S3.
Finally, very few biological processes were significantly enriched in the genes downregulated at pH 7.4 (GO term enrichment analysis, Figure 3). However, these included the category 'posttranslational protein modification' (GO term enrichment analysis; p,0.05), which contains genes encoding kinases. Because of their important role in signal transduction, four of those were selected, Figure 1. Differential gene expression between control samples and samples exposed to pH 7.4 for three weeks. The detection of genes differentially expressed in the two conditions was performed using the DEGseq R package and plotted as a MA plot. The M axis is the log2 fold change for the pH-treated sample compared with the control sample (log2 fold change = log2 [RPKM acid stress/RPKM control]), and the A axis is the average log2 normalized counts in both samples. Each point represents a single contig of the reference transcriptome. The red points correspond to contigs significantly differentially represented between the two conditions (p,0.0001). Among the 57,259 contigs analyzed, 5758 were significantly downregulated (10%) and 3204 (6%) were upregulated. doi:10.1371/journal.pone.0058652.g001 Figure 2. Validation of the RNAseq approach using q-RT-PCR. Twenty-five genes that whose expression was significantly different in the control and treatment conditions were arbitrary selected from highly upregulated to highly downregulated contig. Their levels of expression were quantified by q-RT-PCR, and the results were compared with those obtained using the RNAseq approach. The log2 change in expression of q-RT-PCR and RNAseq was closely correlated (r 2 = 0.86; p,0.0001), indicating the accuracy of the RNAseq approach for quantification. doi:10.1371/journal.pone.0058652.g002 two encoded tyrosine protein kinase while the other two genes encoded serine/threonine protein kinases. These four genes were downregulated at pH 7.4 and 7.2 (between 22.01 and 1.82 fold decrease). For a complete description of candidate genes expression see the Figure 4 and File S3.

Discussion
Validation of the RNAseq approach to the study of coral responses to low pH An RNAseq approach was used to compare the transcriptomes of corals maintained under seawater pH conditions with those of corals exposed to pH 7.4 (an extreme value of today or near future natural variability) for a 3-weeks period. Based on our reference transcriptome, the constant expression of 84% of the genes evidenced that Pocillopora damicornis transcriptome is not completely disturbed by the treatment and that its response to pH 7.4 is mediated by the modification of 16% of its transcriptome. Complementary analysis of these 16% of genes allowed us to identity several molecular pathways and candidate genes putatively involved in the response to low pH. In order to better understand the response of these affected biological functions at a pH 7.8 and 7.2, we measured the transcript abundance of genes belonging to these affected functions by q-RT-PCR method.
The quantification of gene expression by RNAseq requires a reference genome/transcriptome for the mapping step [53]. However, such molecular data are scarce, especially for nonmodel species. This explains why few global transcriptomic analyses have been undertaken using non-model species. Among scleractinian corals such analyses were only possible for a few species such as Acropora millepora and A. digitifera because their transcriptome or genome data were available [27,63,64]. Recent progress in the development of next generation sequencing made it possible to deal with this problem. In a two-step process, the sequences obtained under different physiological conditions were first used to assemble the transcriptome. Secondly, the sequences for each condition were used for quantification of gene expression. The use of transcripts from different environmental conditions allowed us also to cover the largest possible range of mRNAs resulting from differential transcription and alternative splicing [65,66]. Comparison of the quantitative data obtained using the RNAseq and the q-RT-PCR approaches (the latter performed on 25 arbitrarily selected genes) showed a significant correlation between these two methods of quantification. The correlation obtained (r 2 = 0.86) is comparable to those reported in other studies by Oliver and collaborators [67], r 2 = 0.83; Wang and collaborators [68], r 2 = 0.83; Castruita and collaborators [69], r 2 = 0.95; and Meyer and collaborators [27] r 2 = 0.74). This confirms that the two-step RNAseq approach we used is suitable for measuring gene expression in corals exposed to various environmental conditions.

Specific upregulation of biological functions under low pH
The overall pattern of gene expression showed that 6% of the genes were significantly upregulated during exposure to pH 7.4 compared to pH 8.1. One important finding is that genes coding for organic matrix proteins and proteins involved in the transport of Ca 2+ and HCO 3 are modulated under low pH conditions, suggesting that low pH affects the calcification process at the transcriptomic level. We showed that genes coding for HCO 3 2 and Ca 2+ transporters and carbonic anhydrases were upregulated at low pH values (7.8 and 7.4) but downregulated at the extreme level of pH 7.2 ( Figure 4). A similar effect was observed for genes encoding skeletal organic matrix proteins (similar to BMP1 and BMP7 of the bone morphogenic protein superfamily, and the galaxin superfamily).
These results contrasts with some of those obtained in two recent works: one studying by RNAseq the short term (3 days exposure) response of Acropora millepora primary polyps to pH 7.96 and 7.86 compared to pH 8.16 [23], and the second one studying by microarray the short and midterm (1 day and 28 days) responses of adults A. millepora to pH 7.85 and 7.65 compared to pH 8.1 [17]. Among these three studies, if the two using A. millepora as a model species are the most realistic in terms of pH exposure (pH,7.9), ours offer the opportunity to study the response to pH below 7.6 (7.4 and 7.2), which can exacerbate the transcriptomic response. Concerning the comparison with the study of Moya and collaborators [23], in condition where pH between studies is the closest (7.8), their results are in agreement with ours for the genes coding for Ca 2+ transporters and Galaxinlike protein. However, carbonic anhydrase and some organic matrix encoding genes are downregulated in their experiment while it is the contrary in ours. These contradictions are not surprising since we worked at a different time-scale, on a different species and on adult corals whereas they worked on early aposymbiotic life stages. However, even when considering the same species A. millepora, Moya and collaborators [23] have obtained results different from the ones obtained by Kaniewska and collaborators [17]. For example, the genes coding for organic matrix proteins (such as galaxin) clearly shows that effect of pCO 2 is dependent on the life stage. Then differences obtained between our study and the study of Kaniewska and collaborators [17] mainly concern the opposite pattern of expression for the metabolic related genes. Indeed, we measured upregulation while they highlight downregulation. In this last case the contradiction reported may be the results of i) phylogenetic differences (Acropora belongs to the ''complex'' coral clade whereas Pocillopora belongs to the ''robust'' coral clade; [24,25]), ii) and/or experimental protocol differences (light and nutrition were different in the two studies) iii) and/or differences in previous life history (corals are from different geographical origin). Such differences in the response of marine invertebrates belonging to the same phyla were also previously described. In the sea urchin group, Paracentrotus lividus and Strongylocentrotus purpuratus larvae showed mixed transcriptomic responses for ortholog genes when larvae are submitted to a comparable exposition to low pH. In S.purpuratus the genes involved in calcification processes msp130 and SM30 were down-regulated at pH 7.7 [70], while the same gene in P.
lividus were not at pH 7.7, but up-regulated at pH 7.0 [71]. In coccolithophores, a decreasing calcification between ancient (past 400 years) and actual specimens was shown [72]. However they also discovered a heavily calcified coccolithophore in modern water with low pH, Emiliana huxleyi morphotype evidencing variation in responses of phyla to environmental forcing factors [72]. In the case of marine mollusks, low pH induced mainly a decrease of calcification both at the transcriptomic and the animal level, but some exceptions were reported [73][74][75][76][77]. These data highlight the difficulties to extrapolate the results obtained from one species to the whole phylum.
In our study, we hypothesize that the upregulation of genes encoding proteins involved in ion transport counteract the negative effect of low pH on the calcification process ( Figure 5). If the aragonite saturation state (V arag ) decreased in the extracellular calcifying medium (ECM) due to a decreasing pH in the ECM [78,79], the upregulation of ion transporters would allow to increase ion concentration and to maintain the same calcification efficiency in a less favorable pH conditions ( Figure 5). Upregulation of genes involved in ion transport could indeed help to increase the saturation state of aragonite, which depends on the concentration of Ca 2+ and CO 3 22 ( Figure 5). The increase in expression of genes coding for organic matrix proteins needs to be further investigated. However, we postulate that they may modify nucleation, crystal growth inhibition and orientation, and because of their negative charge help to locally increase V arag by trapping Ca 2+ ions [1,56]. It is also possible that the relative composition of the skeleton (inorganic versus organic matrix) could change under low pH; an increase in the organic matrix proportion could help sustain coral growth ( Figure 5). In Figure 4. Gene expression for key biological functions following exposure to various pH levels for three weeks. The data included q-RT-PCR results for samples exposed to pH 7.8 and 7.2, and RNAseq results for samples exposed to pH 7.4. Quantification was normalized by comparison with results for exposure to pH 8.1 (present seawater pH); the results are presented as a log2 fold change in expression. doi:10.1371/journal.pone.0058652.g004 summary, the concomitant and putative additive effect of the increase in ion transport and deposition of organic matrix proteins may help P. damicornis to cope with low pH exposure ( Figure 5). This hypothesis is supported by recent studies performed on marine calcifying invertebrates at several levels of organization. In larvae of the sea urchin P. lividus, there was no effect of low pH on the calcification rate at the organism level with an upregulation of biomineralization related genes at the transcriptomic level [71]. In the mollusk Mytilus edulis exposed to low pH, gene transcription of a tyrosinase putatively involved in the periostracum formation was strongly upregulated and organic matrix protein coding gene remained expressed at a high level. Both mechanisms were associated to the protection and the maintenance of shell formation [74]. In the red abalone Haliothis rufescens, organic matrix protein coding genes remained expressed at a high level during exposure to low pH [80]. At the organism level examples of calcification increase were reported in several organisms such as Crepidula fornicata, Arbacia punctulata [77]. Finally, a recent work on the coral Stylophora pistillata has shown that between pH 8.1 and pH 7.4, the animal controls sufficiently the pH in the ECM to allow the maintenance of calcification [81]. Even if all these examples strengthen our hypothesis, it should be confirmed in P. damicornis further by direct quantification of skeletal Ca 2 + incorporation, organic matrix skeletal content, and measurements of pH in the ECM.
In tropical corals energy is supplied by prey capture (heterotrophic nutrition) and translocation of photosynthates from symbionts to their host [82]. The upregulation of genes involved in heterotrophic and autotrophic pathways may reflect an increase in coral energetic requirements under low pH stress. In this context, it is known that calcification is an energetically costly process consuming up to 30% of the daily energy budget [1]; the need to maintain calcification under stress conditions could in part explain the increase in energetic needs and energy production ( Figure 5). In this context, the upregulation of lipolysis and betaoxidation metabolic pathways shows that corals may use their fatty acid reserves under low pH conditions. However, these reserves are not inexhaustible and should be replenished, whereas global climate change is predicted to cause a decrease in ocean primary production [5,83,84]. Thus, it is uncertain whether corals can maintain their energy production from heterotrophic and autotrophic resources when primary production is predicted to decrease.
Usually high CO 2 concentration driving seawater pH decrease is known to induce bleaching of photosymbiotic organisms such as reef-building coral and foraminifera [17,22,85], even if the photophysiology of the phototrophic symbiont shows mixed responses to this treatment (see introduction section). Our study is the first to show the upregulation of several photosynthetic related genes at pH,8.1. The lack of photophysiological data in the present study makes it difficult to interpret such upregulations but argue in favor of an increasing photosynthetic capability of the holobiont P. damicornis under pH of 7.8, 7.4 and 7.2. Such a positive effect of high CO 2 concentration on photosymbiotic species just begins to be revealed in the literature. Indeed, the acoel worm Symsagittifera roscoffensis entering in symbiosis with the microalgae Tetraselmis convolutae shows a remarkable symbiosis stability even at a very low pH level (up to pH 6) [86]. In a second example, the sea anemone Anthopleura elegantissima/Symbiodinium muscatinei photosymbiosis, high pCO 2 level (pH 8.08 and pH 7.35) induces a higher rate of photosynthesis and mitotic index of the algae compared to pH 8.1 [87]. Under such conditions the cnidarian host received more of their respiratory carbon from the symbiont that under the actual pH condition [87]. These results taking together confirmed the hypothesis that photosymbiosis could be resistant to high pCO 2 , and that the negative effect of this high pCO 2 could be the results of indirect impact at other levels (hypothesis proposed in [86]).

Uncoordinated downregulation of biological functions under low pH
During exposure to pH 7.4 compared to 8.1, 10% of genes were downregulated and 6% were upregulated. Although more genes were downregulated, the enrichment analysis shows that the upregulation of genes was more readily explained and informative. A total of 32 biological functions were upregulated under low pH conditions (p.0.01), while only 1, the 'cellular macromolecule metabolic process' was downregulated ( Fig. 3; p.0.01) despite a comparable number of contigs with a Gene Ontology term (1696 for the downregulated set, 1495 for the upregulated set). Such a results may reflect that the downregulation measured are un or less coordinated than the upregulation. This kind of uncoordinated downregulation of genes has been reported in several other animal species exposed to stress [88][89][90]. A similar response was also observed in the coral Montastraea faveolata subject to thermal stress [91]. Downregulation could reflect a trade-off mechanism induced by the low pH treatment, providing energy savings necessary for efficient stress responses ( Figure 5). This hypothesis is supported by the results of the q-RT-PCR experiments performed on selected host genes. Indeed, the level of downregulation increased concomitantly with the pH level. At pH 7.2 the downregulation was maximal and, as it has been observed and hypothesized in other species subjected to extreme physiological stress conditions, this may reflect a trade-off mechanism [89,[91][92][93][94].
Post-translational protein modification was one of the biological functions downregulated (GO term enrichment analysis, p,0.05) during low pH treatment, as evidenced by the substantial downregulation of several genes encoding protein kinases (serine/threonine kinase and tyrosine kinase). Protein kinases phosphorylate proteins to modulate their activity, and play key roles in signal transduction pathways [95]. The four protein kinases selected for q-RT-PCR analysis showed similar patterns of expression, with upregulation occurring at pH 7.8 and downregulation occurring at pH 7.4 and 7.2. This variable pattern of transcription could reflect the disturbance or molecular plasticity of signaling in coral cells under stress. Such results are strengthened by a proteomic study of the barnacle Balanus amphitrite exposed to low pH. Indeed, under low pH condition  3 2 transport to the extracellular calcifying medium (ECM) increase through an increasing number of Ca 2+ and HCO 3 2 transporters. The rate of CO 2 conversion to HCO 3 2 increases through an upregulation of carbonic anhydrase (CA). Production and secretion of organic matrix protein (OMP round green shape) also increase. This process may help to maintain a sufficient aragonite saturation state (V arag ) despite a decrease in pH, enabling to maintain the same calcification in less favorable pH conditions. This process is energetically costly and need an increase in energy production (ATP) through an increasing metabolic activity and a general trade off mechanism. Finally, at a pH 7.2, the coral cannot cope anymore with such a low pH and there could be a physiological collapse of the colony. Values of omega aragonite (V) and pH in the ECM at pH 8.1 are from [97]. doi:10.1371/journal.pone.0058652.g005 the phosphoproteome is modified revealing a regulation of kinase activity and/or synthesis [96]. The biological effects of all these downregulation events were difficult to resolve because of the very large number of pathways in which kinases are involved, and this will require further study.

Supporting Information
File S1 Biological function, annotation, BlastX top hit species, BlastX evalue, contig number and primer sequence of the selected candidate genes.

(DOCX)
File S2 Table results of the GO term enrichment analysis.

(DOCX)
File S3 Numerical value corresponding to the qRT-PCR. Quantification was normalized by comparison with results for exposure to pH 8.1 (present seawater pH); the results are presented as a log2 fold change in expression. (DOCX)