The experimental setup consisted of two continuous-flow 20L tanks (control and treatment) connected to one of the (360 L) aquarium systems available at the Molecular Geo- and Palaeobiology Laboratory of the Department of Earth and Environmental Sciences, Palaeontology and Geobiology, LMU München (Supplementary Figure 1). Continuous water flow between the experimental tanks and the main 360 L tank was achieved by pumping 1.0–1.6 L/min of water from the main tank to the bottom of the experimental tanks and allowing an equal amount of water from the surface of the tank to return to the main water tank. To simulate ocean acidification conditions, a pH-electrode LE4099 Mettler Toledo connected to a pH computer was used to control the injection of CO₂ into the water of the treatment tank to be kept pH at 7.6 (treatment tank). We used a water pump set to a flow rate of 300 L/h to keep the water circulating in each experimental tank and allow for better mix the CO₂ injected into the treatment tank. Water returning from the treatment tank was equilibrated to pH ~8.0 by letting it flow through a 25 cm column containing limestone grains with diameters ranging from 2 to 5 mm and simultaneously injecting air (300 to 500 L/h). After this, the treated water was mixed with the remaining water in the main (360 L) tank. Both tanks were illuminated using a Mitras LX 6200-HV lamp simulating tropical light cycles; illuminance ranged between 19 and 22 kilolux at the seawater surface of both tanks. During the experiment, water pH in each experimental tank was controlled every day and was recorded every hour from the second day onwards with a pH electrode LE4099 Mettler Toledo and a pH-meter PCE-PHD 1. The pH in the control aquarium fluctuated between 7.90 and 8.26 (mean = 8.08, sd = 0.11) according to the respiration/photosynthesis cycles of coral and symbionts. Because the lowest pH value reached under control condition was close to the worst-case scenario (pH = 7.8) predicted for the end of the century (Meehl et al., 2007), we set the pH in the experimental tank to a lower value (~7.6), which stayed relatively constant (mean = 7.64, sd = 0.02).

Coral nubbins (n = 20) were obtained from branches of two adult M. digitata colonies from our coral culture stock. The colonies were originally donated by Prof. Christian Laforsch. The individual explants were cultured for over 2 months in the main (360 L) aquarium under the same light (19–22 kilolux) and temperature (average ~26°C) conditions used for the experiment. The coral nubbins were not fed either previous to or during the experiment. Approximately 1 month before starting the pH experiment, the nubbins were transferred to the control aquarium of the experimental setup for acclimation under the same light, temperature and feeding conditions as mentioned above. Afterwards, half of the corals were moved to the treatment tank and sampling was conducted after 3 and 42 days of exposure. Five corals were sampled from each tank at each sampling time.

To determine the phylogenetic affinities of the Symbiodinium sp. inhabiting in the colonies of M. digitata used for all experiments, a partial sequence (664 bp) of the ITS region was amplified using primers and PCR conditions detailed in Baillie et al. (2000). ITS sequences from Symbiodinium clades A to G were downloaded from GenBank (accession numbers in Supplementary Figure 2) and aligned with MUSCLE (Edgar, 2004). A maximum likelihood tree was then inferred using PhyML (Guindon et al., 2010) with a GTR substitution model, 100 bootstrap replicates, and optimized nucleotide equilibrium frequency, invariable sites and across-site rate variation. Both analyses were done in SeaView v4.6.2 (Gouy et al., 2009).

At each sampling point, the nubbins were flash-frozen in liquid nitrogen, broken in smaller pieces with mortar and pestle (always in liquid nitrogen) and stored at −80°C until further processing. For RNA extraction, samples were homogenized in TRIzol Reagent (Thermo Fisher Scientific) with a SilentCrusher M and a Dispersion Tool 8 F (Heidolph) and total RNA was extracted using a slightly modified version of the Chomczynski method (Chomczynski and Sacchi, 1987). RNA quality and concentration were controlled with a NanoDrop ND-1000 spectrometer and with a BioAnalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). An additional cleaning step with Agencourt RNAClean XP Magnetic Beads was done if required. RNA concentration and quality requirements for library preparation and sequencing were achieved for all samples obtained after the first sampling time, for the second sampling time only 6 samples (3 control + 3 treatment) met quality requirements and were further used. Library preparation and sequencing (50 bp PE; HiSeq2500) was done at the GeneCore of the European Molecular Biology Laboratory (EMBL) in Heidelberg, Germany. The reads generated for both control and treatment samples of the second sampling time were strand-specific and further used for de novo transcriptome assembly. The reads generated were uploaded to the Short Read Archive of the European Nucleotide Archive under the project accession PRJEB21531 (http://www.ebi.ac.uk/ena/data/view/PRJEB21531).

Illumina read pairs were quality controlled with FastQC v0.63 (Andrews, 2010), low quality bases and reads were removed with Trimmomatic v0.32 (Bolger et al., 2014). The surviving pairs were further processed with the program Filter Illumina v0.40 of the Agalma/Biolite suite (Dunn et al., 2013) to ensure the removal of adapters and low quality reads. Reads from putative bacterial contaminants were filtered out by aligning the surviving pairs against all bacterial genomes on RefSeq as of July 2015 (ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria) and discarding mapping pairs.

De novo transcriptome assembly was done with Trinity v2.0.6 using only strand-specific sequences (i.e., the surviving read pairs from the second sampling time from both conditions) (Grabherr et al., 2011). The assembled contigs from coral and dinoflagellates were separated with Psytrans (https://github.com/sylvainforet/psytrans) using the Acropora digitifera v0.9 (Shinzato et al., 2011) and Symbiodinium minutum Clade B1 v1.0 (Shoguchi et al., 2013) predicted proteins as references. The minicircle sequences of Symbiodinium sp. Chloroplast (Barbrook et al., 2014) were recovered manually from the contigs assigned to the coral. Coding sequences (CDS) from the assembled transcripts were predicted with TransDecoder v2.0.1 (https://transdecoder.github.io). Only transcripts with predicted proteins were used as the reference transcriptome for both coral and symbionts. Completeness of the transcriptomes was assessed by search of the sequences against the Core Eukaryotic Genes Dataset (CEGMA) (Parra et al., 2007; Francis et al., 2013). Data processing and transcriptome assembly were mostly executed on the Molecular Geo- and Palaeobiology Lab's Galaxy server (Goecks et al., 2010). Differentially expressed genes (as defined in the next section; DEGs hereafter) were annotated by BLASTN (E = 10⁻¹⁰) against the Non-Redundant (nr) database of the NCBI and by BLASTX (E = 10⁻¹⁰) against the SwissProt database (both downloaded in February 2015). The best SwissProt hit of each DEG was used to retrieve its associated Gene Ontology terms through the QuickGO online tool (https://www.ebi.ac.uk/QuickGO). In addition, predicted proteins of the DEGs were annotated with the protein identifiers of the A. digitifera v0.9 and S. minutum Clade B1 v1.2 by BLASTP (E = 10⁻¹⁰). A Galaxy history containing some of the steps previously mentioned as well as transcriptome annotations can be found in https://github.com/PalMuc/Montipora_digitata_resources.

For each sampling time, the remaining reads from each sample were mapped with Bowtie2 (Langmead et al., 2009) to the reference transcriptomes of M. digitata and Symbiodinium sp. Read counts from isoforms groups were added to derive a matrix of counts per Trinity-component. Differences between conditions and colonies were assessed with the adonis function of the vegan R library (Oksanen et al., 2013). As in other ocean acidification studies (Davies et al., 2016; Kurman et al., 2017), the sample-dependence effect caused by one single tank per condition could not be corrected (Riebesell et al., 2011; Cornwall and Hurd, 2016; but see Oksanen, 2001 and Schank and Koehnle, 2009 for a discussion on pseudoreplication). DEGs between conditions were detected using the package DESeq2 v1.8.1 (Love et al., 2014). DEGs were selected by their significant change in expression as assessed by a Wald test; P-values were adjusted with the Benjamini-Hochberg correction (Benjamini and Hochberg, 1995). Gene Ontology (GO) terms enrichment was assessed with topGO v1.0 (Alexa and Rahnenfuhrer, 2010) in R v3.2.1. Protein domains of the DEG products were identified using PfamScan and PFAM A release 30.0 (both downloaded from ftp://ftp.ebi.ac.uk/pub/databases/Pfam). The counts and the R scripts used to analyze the data can be found in the project repository (https://github.com/PalMuc/Montipora_digitata_resources).

The M. digitata meta-transcriptome was assembled using over 190 millions of strand-specific read pairs. The transcriptome assembly yielded 179,298 sequences, of which 123,710 (69%) were assigned to the coral host and 55,588 (31%) to the symbiont. These transcript sets resulted in 41,852 and 34,057 predicted CDS for the coral and symbionts, respectively. The N50 length was 1,101 bp for the coral and 744 bp for Symbiodinium sp. G+C content of both reference transcriptomes were 43.9 and 54.5%, consistent with previous reports for the Order Scleractinia and Symbiodinium spp., respectively (Bayer et al., 2012; Shoguchi et al., 2013; Sun et al., 2013; Shinzato et al., 2014). Summary quality statistics for both transcriptomes are summarized in Table 1. Most of the Core Eukarytic Genes (CEGs) were found in the transcriptomes, 89.5% in M. digitata and 79.4% in Symbiodinium sp. (Supplementary Table 1). In addition, the level of fragmentation of the assembled transcriptomes appears to be low as judged by the low number of CEGs that are longer than the queries.

Only 5,541 (13.2%) out of the 41,852 genes contained in the transcriptome of M. digitata transcriptome could not be annotated with any of the databases. A large portion (48.6%) of transcripts had hits against sequences from all three databases used and about 15% of the total transcripts hit exclusively sequences of the A. digitifera genome (Figure 1). On the other hand, most of the sequences in the transcriptome of Symbiodinium sp. were annotated with genes from the S. minutum genome (85.5% of all transcripts), and almost half of them (48.4%) had hits exclusively against this genome. Similarly to the coral transcriptome, only a small fraction of the transcripts (12.9%) lacked any annotation (Figure 1).

No significant differences in global expression were found with the adonis exploratory test between conditions and between colonies for neither the short (p_condition = 0.40, p_colony = 0.19) nor the longer (p_condition = 0.73, p_colony = 0.76) exposures of the coral to low pH. These findings are also supported by the variance in expression of the coral genes at both sampling times (Figure 2) and by the variance-corrected counts of mapped reads per gene per sample (Supplementary Figure 3).

Still, the DeSeq2 gene expression analysis yielded 18 genes with significant (p_adj < 0.05) differences in expression level due to exposure to low pH after 3 days (model design: ~ colony + condition, Figure 3). The log₂-fold changes for these genes ranged between −0.70 and 0.73 (Supplementary Table 2). Six up- and twelve down-regulated genes were found. From the six upregulated genes, two (including a metalloprotease) contain zinc-binding domains and others have RNA/DNA binding domains. Among the downregulated genes, proteins similar to Noelin, to an αB-crystallin and to a “platelet-derived” growth factor were found. A protein with no annotation but containing a domain of the tumor necrosis factor receptor superfamily (TNFR) was also downregulated. GO terms related to organ development, cell differentiation and cell chemotaxis were enriched in this set of genes (Supplementary Table 2). Upregulated genes were rich in GO terms related to the negative regulation of calcium ion-dependent exocytosis, regulation of extracellular matrix assembly, collagen catalysis, changes in bone morphogenesis, negative regulation of cell-substrate adhesion, transmembrane transport of zinc and response to antigenic stimulus. Additionally, the Notch-signaling pathway appeared to be positively regulated, contrary to previous findings (Kaniewska et al., 2012). No DEGs were found in the coral after 42 days of exposure to low pH.

The exploratory analyses of gene expression of Symbiodinium sp. suggest no change at any sampling time (Figure 4). The adonis permutation test resulted in no significant differences between conditions and colony after the first (p_condition = 0.62, p_colony = 0.78) or the second (p_condition = 0.88, p_colony = 0.99) exposure period. No DEGs were found at any time.

The present study used RNA-Seq to assess the coupled response of an adult stony coral (M. digitata) and its microalgal symbionts to low pH. Compared to previous investigations (Kaniewska et al., 2012; Moya et al., 2012, 2015), our results suggest a mild gene expression response of the adult coral holobiont to a 3-days exposure to acidified water in terms of number of DEGs and log₂-fold change. Still, the log₂-fold changes we found in this study are comparable to those observed in the resilient coral Siderastrea siderea under laboratory-induced ocean acidification conditions (Davies et al., 2016). This suggests resilience of M. digitata to low pH.

The few DEGs found after short-term exposure of M. digitata to low pH (7.6) have been shown to be involved in stress response in corals. Zinc-metalloproteases, for instance, seem to be relevant for heat tolerance in these organisms and might play a role in regulation of apoptosis and cell repair (Barshis et al., 2013). Other stress related genes were found among the significantly downregulated genes, like the αB-crystallin-like protein, an HSP20 family member. αB-crystallin seems to play a role in stress response in other corals (e.g., Orbicella annularis and O. faveolata), but its exact function is yet unknown (Downs et al., 2002). A TNFR domain-containing protein is another example, and changes in expression levels of proteins carrying this type of domain display have been found in other corals under stress conditions (Barshis et al., 2013; Seneca and Palumbi, 2015; Yuan et al., 2017). Another down-regulated gene codes for a protein containing a platelet-derived growth factor domain, though the regulation of this kind of protein upon exposure to low pH has not been previously reported in corals. However, it is already known that homologs of human bone-morphogenesis proteins (BMP2/4) participate in skeleton production of marine calcifiers, including corals (Green et al., 2013). In fact, proteins carrying this domain are currently used for bone regeneration in humans, and its use, in combination with coral skeleton, has been proposed to treat bone defects (Parizi et al., 2012).

According to previous investigations, coral larvae can acclimate to acidified water within a few days (Moya et al., 2015). Such a rapid acclimation by the adult coral holobiont would thus explain the lack of change in gene expression after 42 days of acidified water exposure in our experiments, a time by which acclimation would likely have been already completed. Our findings in the long-term response are similar to those of a previous investigation on the response of another adult scleractinian (Acropora millepora) to a 37-days exposure to low pH (elevated pCO₂) that reported absence of change in the expression of target genes associated with calcification and metabolism (Rocker et al., 2015). In agreement with this, previous studies suggest that temperature stress has a greater impact in gene expression of corals than pH (Mayfield et al., 2014; Davies et al., 2016), although resilience to stress conditions varies from one species to another (Loya et al., 2001), and even among populations of the same coral species occurring in habitats with different local conditions (Poli et al., 2017).

Physiological response to stress conditions does not only depend on coral host resilience but also on the interaction with its symbionts (Hoadley et al., 2015). The Symbiodinium sp. studied here, which likely belongs to clade C and is specific to M. digitata (Supplementary Figure 2), did not seem to be affected by the decrease in pH. The more stable environment provided by the coral tissue, which acts as some kind of protection for the algae may serve as an explanation for this observation (Banaszak and Trench, 1996). Nonetheless, the study of the ex hospite transcriptomic response to heat stress of Symbiodinium in the same clade revealed distinct levels of heat-stress tolerance among different populations of this symbiont (Levin et al., 2016). These findings indicate that the tolerance level of different Symbiodinium types and populations to environmental stressors may contribute significantly to the resilience of particular coral holobiont assemblages (Yuyama et al., 2012; Putnam et al., 2017). Still, differences in response to stress between Symbiodinium spp. in hospite and isolated from the host remain largely unexplored, as does the contribution of the resilience of the symbiotic partner to the response of the holobiont.

An alternative explanation for the lack of change in overall gene expression of the symbionts in this study might lie on the constitutive expression of most genes observed in dinoflagellates (Moustafa et al., 2010) and on the reliance on post-transcriptional regulation of gene expression in these organisms (Baumgarten et al., 2013; Liew et al., 2017). However, little is known about the mechanisms of gene expression regulation in Symbiodinium and this topic needs to be more actively investigated given its potential role in holobiont resilience.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer SF and handling Editor declared their shared affiliation.