Epigenome-associated phenotypic acclimatization to ocean acidification in a reef-building coral.

There are increasing concerns that the current rate of climate change might outpace the ability of reef-building corals to adapt to future conditions. Work on model systems has shown that environmentally induced alterations in DNA methylation can lead to phenotypic acclimatization. While DNA methylation has been reported in corals and is thought to associate with phenotypic plasticity, potential mechanisms linked to changes in whole-genome methylation have yet to be elucidated. We show that DNA methylation significantly reduces spurious transcription in the coral Stylophora pistillata. Furthermore, we find that DNA methylation also reduces transcriptional noise by fine-tuning the expression of highly expressed genes. Analysis of DNA methylation patterns of corals subjected to long-term pH stress showed widespread changes in pathways regulating cell cycle and body size. Correspondingly, we found significant increases in cell and polyp sizes that resulted in more porous skeletons, supporting the hypothesis that linear extension rates are maintained under conditions of reduced calcification. These findings suggest an epigenetic component in phenotypic acclimatization that provides corals with an additional mechanism to cope with environmental change.

INTRODUCTION

Over the last century, the anthropogenic production of CO2 has led to warmer (+0.74°C) and more acidic (−0.1 pH) oceans (), resulting in increasingly frequent and severe mass bleaching events worldwide that precipitate global coral reef decline (, ). To mitigate this decline, proposals to augment the stress tolerance of corals through genetic and nongenetic means have been gaining traction (), in the light of environmentally induced alterations in DNA methylation leading to potentially heritable phenotypic acclimatization in model systems (, ).

Our organism of interest, Stylophora pistillata, is a globally distributed scleractinian coral with an available draft genome (). Previous work has demonstrated its plasticity and resilience in the face of high pCO2 (partial pressure of CO2) conditions (, ). Remarkably, this coral remains capable of calcifying in seawater at a significantly reduced pH of 7.2 even when aragonite, the main component of coral skeletons, falls below its saturation point in seawater (that is, Ωaragonite < 1) (). Previous work showed that the linear extension rate of S. pistillata in acidic seawater is not significantly different from that measured under control conditions (pH 8.0), despite the significantly reduced calcification rates (). Here, we aimed to investigate whether epigenetic mechanisms regulate phenotypic acclimatization to long-term pH stress in corals.

RESULTS AND DISCUSSION

S. pistillata cultured in varying pH conditions indicate long-term impact of pH stress on coral colonies

S. pistillata colonies were cultivated in an ongoing aquarium experiment with four experimental conditions: seawater pHs of 7.2, 7.4, 7.8, and 8.0 (control). Conditions in these aquaria were identical, except for the pCO2 and resulting differences in seawater carbonate chemistry of each tank (Table 1). Coral fragments of the same genet were added to the system annually. Fragments used for whole-genome bisulfite sequencing (WGBS) and RNA sequencing (RNA-seq) were added to the system in 2012 and collected in 2014. Similarly, samples for the follow-up laser microdissections, amplicon-specific bisulfite sequencing, and porosity measurements were added in 2013 and collected in 2015. Calyx and cell size measurements were performed on samples added in 2014 and 2015 and collected in 2016 and 2017, respectively. In all cases, sampled fragments were subjected to the respective pH treatments for approximately 2 years.

Table 1

Carbonate chemistry parameters in the four experimental seawater pH treatments.

Parameters of carbonate seawater chemistry were calculated from total scale pH, total alkalinity (TA), temperature, and salinity using the open-access CO2SYS package () using constants from the study of Mehrbach et al. () as refit by Dickson and Millero (). Values shown are means ± SD. pHT, pH in treatment tanks; Ωar, saturation state of aragonite.

Treatment name	pHT	TA(mmol/kgSW)	pCO2(μatm)	HCO3−(μmol/kgSW)	CO32−(μmol/kgSW)	Total carbon(μmol/kgSW)	Ωar
7.2	7.22	2496.62	3513.15	2384.13	46.25	2528.27	0.72
	±0.01	±6.99	±67.20	±4.42	±1.08	±3.62	±0.02
7.4	7.43	2474.16	2109.3	2299.43	71.66	2429.86	1.11
	±0.01	±5.89	±59.63	±0.76	±2.13	±1.23	±0.03
7.8	7.81	2461.77	798.77	2082.22	155.17	2259.64	2.41
	±0.01	±5.97	±14.36	±0.60	±2.75	±1.75	±0.04
8	7.95	2447.48	537.66	1951.14	202.41	2168.53	3.14
	±0.00	±7.16	±4.88	±1.90	±2.25	±4.02	±0.03

Gene body methylation reduces spurious transcription and transcriptional noise

We performed WGBS on DNA extracted from three replicate nubbins per tank and obtained data from 98% of all CpGs in the genome with a per-sample, per-position mean coverage of ~25× (discussion S1 and data S1). The S. pistillata genome is sparsely methylated [1,406,097 base pairs (bp), 7% of all CpGs], similar to that of other invertebrates, for example, 1% in the bee Apis mellifera (), 2% in the wasp Nasonia vitripennis (), and 12% in the termite Zootermopsis nevadensis (). The vast majority of cytosine methylation in S. pistillata occurs in genic rather than intergenic regions (76.4% versus 23.6%; Fig. 1A), similar to other invertebrates (–). There were more methylated positions residing in introns than in exons (801,767 versus 270,806; Fig. 1A), similar to genomes of Z. nevadensis (), but unlike those of other invertebrates (, ). Surprisingly, this remained true even when we accounted for the longer lengths of introns: A higher proportion of CpGs in introns were methylated than in exons (11.3% versus 8.6%; Fig. 1B), unlike Z. nevadensis ().

Fig. 1

Epigenetic landscape of S. pistillata.

(A) More than half of all methylated positions in S. pistillata are located in annotated introns. (B) Introns have proportionally more methylated positions (11.3%) than exons (8.6%) or intergenic regions (3.3%) even when accounting for the different amounts of CpG dinucleotides in the respective regions. (C) Distribution of methylated positions across a standardized gene model with flanking 4-kb regions. Solid lines depict transcription start site (left) and transcription termination site (right). Exons and introns in the plot have normalized lengths that correspond to their respective mean lengths in S. pistillata (from left to right: 373, 1258, 363, 1114, 302; 253, 967, 312, 1037, and 669 bp). arb. units, arbitrary units.

In contrast to vertebrates, methylation in the putative promoter regions of S. pistillata is scarce and does not seem to affect gene expression (see discussion S2); methylated positions are instead predominantly located within gene bodies (Fig. 1C). Gene body methylation has recently been shown, in mice, to be established via cross-talk between the transcriptional machinery and histone modifications. RNA polymerase II–mediated transcription establishes new methylated positions along the gene body through the action of SetD2, H3K36Me3, and Dnmt3b. Methylation, in turn, reduces spurious transcription from cryptic promoters within these gene bodies ().

Analyses of our RNA-seq data strongly suggest that these functions are conserved in corals. We find that gene body methylation increases with gene expression in S. pistillata (Fig. 2, A and B), similar to other corals (, ) and other invertebrates (, ). We observed that methylated genes show significantly lower levels of spurious transcription relative to unmethylated genes (Fig. 2C). Furthermore, consistent with the repressive effect of methylation on expression, our analysis indicates that methylation in gene bodies also reduces transcriptional noise (lower variability of gene expression levels), echoing previous reports in corals and other organisms (Fig. 2D and discussion S3) (, ). These findings also suggest that the biological functions of gene body methylation are likely conserved across Metazoa, if at all present in the organism.

Fig. 2

Effect of gene body methylation on genic expression.

(A) The expression of methylated genes is significantly higher than that of unmethylated genes (P < 10−300, Student’s t test). (B) Expression values for methylated genes are exponentially proportional to methylation level; however, the relationship of expression to methylation density is nonlinear. Methylation density at low levels is exponentially proportional to expression values, but it plateaus at ~40%. (C) Expression levels of the first six exons were calculated as natural logarithms of fold changes (FC) relative to the expression of the first exon. The difference in expression levels is likely driven by the reduction of cryptic transcription initiation in methylated genes. The difference is greater in highly methylated genes (median methylation level >80%). Asterisks represent P values from t tests of methylated (orange) or highly methylated genes (red) against unmethylated genes (peach), and colored accordingly. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. (D) There is a linear relationship between the inverse of the coefficient of variation and the log10-transformed mean expression values from all samples. The coefficient of variation, defined as the SD of measured expression values divided by their mean, is consistently lower in methylated genes than unmethylated genes.

Differentially methylated genes are linked to general growth and stress responses

On the basis of these findings, we sought to further elucidate the potential role of DNA methylation in phenotypic acclimation. Using generalized linear models (GLMs), we identified genes that undergo differential methylation in response to pH treatment (data S2). In general, we observed a significant increase in overall DNA methylation with decreasing pH, echoing similar observations in Pocillopora damicornis (). To validate these changes in methylation, we performed amplicon-specific bisulfite sequencing of selected genes on the original samples, as well as on samples of an independent experiment. These analyses showed high correlation of results obtained from WGBS and amplicon-specific bisulfite sequencing (r2 > 0.8, P < 0.01; fig. S1A) and further confirmed a high degree of reproducibility of DNA methylation changes across independent experiments (P < 0.01; fig. S1B). Analyses on laser-microdissected oral and aboral tissues further highlighted that most of the selected genes displayed strong and consistent tissue-specific methylation patterns, similar to findings in vertebrates (). These patterns were, in some cases, also correlated with known tissue-specific functions or activities of these genes (discussion S4).

On the basis of the skeletal phenotypes reported in previous studies analyzing the effects of decreased pH levels on corals (, –), we initially expected that many of the differentially methylated and differentially expressed genes would be involved in biomineralization. To our surprise, we did not observe this pattern: Fewer than expected biomineralization-related genes had methylation at all—68 of 148 genes were methylated (46%), which is lower than the global rate of 55% (Fisher’s exact test, P < 0.05). In the functional enrichment analysis of differentially methylated genes, we observed few processes linked to biomineralization, for example, “otolith development,” “calcium ion transport,” and “voltage-gated calcium channel complex” for corals grown at pH 7.2, but these processes were not present at all for corals grown at higher pHs. For the 68 biomineralization-related methylated genes, we observed a minor shift in methylation of organic matrix (OM) and extracellular matrix (ECM) components, with some known OM genes upmethylated and/or up-regulated and others downmethylated and/or down-regulated, implying a possible change in the scaffolding structure of the OM and ECM (discussion S5). As a general trend, we observed that Ca2+-binding OM constituents exhibited increased methylation and expression, whereas genes involved in cell-cell adhesion exhibited decreased methylation and expression.

Instead of biomineralization-related pathways, a functional enrichment analysis of differentially methylated genes revealed processes linked to general growth and stress responses (Fig. 3A and data S3). More specifically, we observed that the methylation levels of many genes in the mitogen-activated protein kinase (MAPK) signaling and cell growth pathways changed significantly (Fig. 3B). Remarkably, we found that the methylation levels of genes involved in the negative regulation of JNK (c-Jun N-terminal kinase) and MAPK increased; conversely, the methylation levels of genes that positively regulate the same kinases decreased (Fig. 3C). For instance, the methylation of JIP1 (JNK-interacting protein 1; SpisGene10613), a negative regulator of JNK, increased in both aboral and oral tissues at pH 7.2, while the methylation of TRAF6 (tumor necrosis factor receptor–associated factor 6; SpisGene580), a positive regulator of JNK, decreased in both tissues. The JNK pathway has previously been shown to control cell, organ, and body sizes in response to stress in mice () and Drosophila (). The consistent differential methylation of JNK effectors that we observed in S. pistillata in response to decreasing pH levels therefore suggested a change in cell and body sizes.

Fig. 3

Differential methylation of genes in the MAPK/JNK signaling and cell growth pathways in the coral S. pistillata under pH stress.

(A) Most of the general growth– and stress response–related Gene Ontology (GO) terms are significantly enriched at pH 7.2. GTPase, guanosine triphosphatase. (B) A composite pathway diagram consisting of KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways “ko04010”: “MAPK signaling pathway” and “ko04110”: “cell cycle” was produced in Cytoscape. At pH 7.2, more genes respond with a significant increase (red) than with a decrease (blue) in methylation; unshaded genes are methylated but not differentially regulated. (C) The methylation of genes that negatively regulate MAPK and JNK (red lines) increases, whereas the methylation of genes that positively regulate the same kinases (blue lines) decreases in a reciprocal manner (error bars denote ±1 SE).

Key genes in MAPK and cell growth pathways are differentially expressed

Analysis of differentially expressed genes in the MAPK and cell growth pathways highlighted several key genes that were differentially expressed in a manner that suggests a delay in cell division and an increase in cell size and body growth. For instance, Mos, a kinase that has been shown to prevent the degradation of cyclin B and thereby delay the onset of the anaphase (), was up-regulated. Similarly, GADD45, a gene that binds to cyclin-dependent kinase 1 (CDK1) and prevents its association with cyclin B, was also up-regulated, potentially delaying the onset of the M phase (). In contrast, Ras, which plays a key role in the progression of the G1 phase, was down-regulated (). The kinases and phosphatases that regulate Ras experienced concomitant changes that indicated an overall reduction in Ras activity (fig. S2). For the latter, we observed expression increases in c-JUN, a transcription factor that promotes cellular growth, and Hsp70, previously observed in coral larvae under prolonged pH stress and presumably expressed to assist in the correct folding of cellular proteins (). The differential expression of these genes was subsequently corroborated via a quantitative reverse transcription PCR (RT-qPCR) analysis (fig. S3 and data S4).

Physiological measurements correlate with changes in methylation and expression

In light of our data, we hypothesized that changes in DNA methylation (Fig. 4A) and expression of key genes correlate with an increase in cell size at pH 7.2. We first confirmed that individual cell sizes were significantly larger in corals grown at pH 7.2 relative to the control (total n = 4728 measurements across 14 nubbins; Fig. 4B and data S5A). We then hypothesized that this cell-level phenotypic change would be accompanied by a corresponding increase in polyp size and, ultimately, skeletal structure. Specifically, we posited that larger cells form larger polyps and larger corallite calyxes (the cup-shaped openings in the skeleton that house the polyps). As measuring the size of polyps is difficult because of their expandable/contractible nature, we opted to measure the size of the calyx as a proxy. As predicted, our analyses confirmed that the calyxes in corals grown at pH 7.2 relative to the control were significantly larger (total n = 181 measurements across 10 nubbins; Fig. 4C and data S5B). Finally, we also confirmed that the skeleton was significantly more porous (n = 6 nubbins; Fig. 4, D and E) in samples from the same treatment. All the observed phenotypes were consistent with our hypothesis.

Fig. 4

Effect of lower pH levels on cell growth.

(A) Mean methylation levels were significantly increased at pH 7.2. (B) Cell sizes were significantly larger in nubbins grown in the pH 7.2 tank. (C) Larger cell sizes were associated with skeletal structure that contained larger calyxes. (D) The skeletal porosity was significantly higher at lower pH. (E) Representative longitudinal sections of S. pistillata skeletons under pHs of 7.2 and 8.0. Error bars represent 1 SE. Asterisks denote significance of t test P values. **P < 0.01; ***P < 0.001.

It is important to note that while our results show a strong correlation between changes in DNA methylation and the resulting phenotype, they do not show a sensu stricto causal relationship between changes in the methylation state and the phenotype. However, the high reproducibility of the changes in methylation and the presence of tissue-specific DNA methylation patterns lend support for a function of DNA methylation in phenotypic plasticity. Our results suggest that the observed phenotypic changes under pH stress are mediated through differential methylation and expression of known stress response pathways that control cell proliferation and cell growth (, ). We propose that these cellular phenotypic changes, together with shifts in OM proteins, are among the drivers of morphological changes in the skeleton of this species under seawater acidification observed here and described previously (). These morphological changes toward a more porous skeleton are possibly a means by which S. pistillata can maintain linear extension rates in the face of depressed calcification rates under seawater acidification (, ). Such a trait would be advantageous in the benthic environment where competition for space and light is an important selective pressure ().

CONCLUSION

In conclusion, our results suggest that DNA methylation could offer corals greater ability to buffer the impacts of environmental changes by fine-tuning gene expression in response to changing environmental conditions, thereby providing additional time for genetic adaptation to occur. Better understanding of the mechanisms underlying coral resilience will also provide additional avenues for reef restoration efforts, such as the human-assisted acclimatization of corals in specialized nurseries [“designer reefs” ()], especially if these modifications can be passed down to future generations. These efforts might prove crucial to averting large-scale losses of extant coral reefs in light of recent global declines due to climate change.

MATERIALS AND METHODS

Definition of “methylation density” and “methylation level”

Methylation density refers to the number of CpG positions that were deemed as methylated within a predefined stretch of DNA sequence. Qualitatively, we used “sparse” and “dense” to describe methylation density. Quantitative definitions differ depending on context: In Fig. 1, for instance, we were interested in looking at whether there were more methylated positions in a particular genomic feature (exons, introns, or intergenic regions). In this case, we displayed, and discussed, the abundance of methylation positions per standardized length of the genomic feature with arbitrary units.

However, when we compared methylation densities across genes, we had had to consider the effect of gene length on the number of methylated positions (the longer a gene, the higher the number of methylated positions), and CpG biases present in genes (comparing genes of equal length, the one with more CpGs in it will have more methylated positions). We standardized it by defining methylation density as the fraction of CpG positions in the gene that were deemed methylated, that is, if a gene has five CpG dinucleotides, there are 10 positions that could be methylated. If 8 positions are methylated, methylation density = 8/10 = 0.8. The range of methylation densities, for this definition, is [0, 1].

Methylation level refers to the abundance of methylation at a particular methylated position. Qualitatively, we used “weak” and “strong” to describe methylation level. Quantitatively, if a particular cytosine has a coverage of 10×, with eight Cs and two Ts mapping to it, its methylation level is 8/10 = 0.8. Thus, methylation level also has a range of [0, 1].

Maintenance of S. pistillata colonies

Colonies of the tropical coral S. pistillata were exposed to long-term seawater acidification in an experimental aquarium setup that has been maintained continuously for several years and described in previous investigations (, ). Briefly, coral fragments were kept in aquaria supplied with Mediterranean seawater (exchange rate of 70% hour−1) at a salinity of 38 g liter−1, temperature of 25°C, and irradiance of 170-mmol photons m−2 s−1 on a 12-hour/12-hour photoperiod provided by HQI-10000K metal halide lamps (BLV Nepturion). Carbonate chemistry was manipulated by bubbling CO2 to reduce the control pH (pH 8.0) to the target values of pH 7.8, 7.4, and 7.2. pH and temperature were monitored continuously by electrodes linked to a monitoring system (Enoleo) that controlled CO2 bubbling. Weekly pH measurements were also made using the indicator dye m-cresol purple (Acros Organics; no. 199250050) adapted from the study of Dickson et al. (). Absorbance was measured using a spectrophotometer (UVmc2; Safas). Measurements of TA were made weekly according to protocols described by Dickson et al. (). TA was measured via titration with 0.03 N HCl containing 40.7 g NaCl liter−1 using a Metrohm Titrando 888 Dosimat controlled by Tiamo software to perform automated titrations of 4-ml samples, and alkalinity was calculated using a regression routine based on U.S. Department of Energy guidelines (). For each sample run, certified seawater reference material supplied by the Dickson Laboratory (Scripps Institution of Oceanography, La Jolla, CA) was used to verify acid normality. Parameters of carbonate seawater chemistry were calculated from total scale pH, TA, temperature, and salinity using the open-sourced CO2SYS package () using constants from the study of Mehrbach et al. () as refit by Dickson and Millero () (Table 1).

Coral sampling

Briefly, the experimental setup described previously () has been ongoing since the early 2010s. Coral fragments of the same genet were added to the system periodically (every ~12 months). For our experiments, coral fragments were sampled from the treatment tanks at different time points; all of these fragments, however, were subjected to their respective pH treatments for approximately 2 years.

Identification of methylated CpGs

DNA was extracted from S. pistillata nubbins (triplicates of four conditions) using the nuclei isolation approach to minimize contamination with symbiont DNA, as previously described (). Briefly, cells from S. pistillata from a nubbin of about 3 cm were harvested using a water pick in 50 ml of 0.2 M EDTA solution refrigerated at 4°C. Extracts were passed successively through a 100-μm and a 40-μm cell strainer (Falcon) to eliminate most of the algal symbionts. Extracts were then centrifuged at 2000g for 10 min at 4°C. The supernatant was discarded, and the resulting pellets were homogenized in lysis buffer (G2) of the Qiagen Genomic DNA Isolation Kit (Qiagen). The DNA was extracted following the manufacturer’s instructions using Genomic-tip 100/G (Qiagen). DNA concentration was determined by optical density using an Epoch Microplate Spectrophotometer (BioTek). Contamination with Symbiodinium DNA was assessed via PCR targeting the multicopy gene RuBisCO (GenBank accession no. AY996050). The negligible presence of Symbiodinium DNA in the extracts was further supported by the low mapping rates (<0.2%; data S1) of the bisulfite-converted reads against the (converted) Symbiodinium microadriaticum genome ().

Bisulfite DNA libraries were prepared following a modified version of the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England Biolabs). Methylated TruSeq Illumina adapters (Illumina) were used during the adapter ligation step followed by bisulfite conversion with the EpiTect Bisulfite Kit (Qiagen), with the following cycling conditions (95°C for 5 min, 60°C for 25 min, 95°C for 5 min, 60°C for 85 min, 95°C for 5 min, 60°C for 175 min, and then three cycles of 95°C for 5 min and 60°C for 180 min. Hold at 20°C for ≤5 hours.). The final library was enriched with the KAPA HiFi HotStart Uracil+ ReadyMix PCR Kit (KAPA Biosystems) following the manufacturer’s instructions for bisulfite-converted next-generation sequencing library amplification. Final libraries were quality-checked using the Bioanalyzer DNA 1000 chip (Agilent) and quantified using Qubit 2.0 (Thermo Fisher Scientific) and then pooled in equimolar ratios and sequenced on the HiSeq 2000 platform. Sequencing of the libraries resulted in 1.53 billion read pairs across 12 samples (data S1). The raw sequences were trimmed using cutadapt v1.8 (). Trimmed reads were then mapped to the S. pistillata genome, deduplicated, and scored on a per-position basis for methylated and unmethylated reads using Bismark v0.13 ().

Lambda DNA, which can be spiked-in to estimate the combined nonconversion and mis-sequencing rate during the bisulfite treatments, was not used in our sequencing runs. However, as we observed that the rates of non-CpG methylation (CHG and CHH, where H = non-G base) were at 0.1% in all samples (data S1), the combination of nonconversion and mis-sequencing would be—at worst—0.1%, if we assumed that CHG and CHH methylation does not occur in S. pistillata.

As an additional guard against false positives, three stringent filters were used to ensure that methylated positions were bona fide. First, the probability of methylated positions arising from chance on a per-position basis was modeled using a binomial distribution B(n, p), where n represents the coverage (methylated + unmethylated reads) and p represents the probability of sequencing error (set to 0.01, an order worse than the worst-case error rate discussed above). We kept positions with k methylated reads if p(X ≥ k) < 0.05 (post–false discovery rate correction). Second, methylated positions had to have at least a methylated read in all three biological replicates of at least one condition. Finally, the median coverage of the position across all 12 samples had to be ≥10. These steps ensured that methylated positions were highly replicable and highly covered.

Assignment of genomic context to methylated cytosines

On the basis of the genome annotation of S. pistillata (in the form of a GFF3 file) and the positional coordinates of the methylated cytosines (in a tab-separated file produced by Bismark), a Python script was written to annotate every methylated cytosine with additional genomic context. The script determined whether the methylated position resides in a genic or intergenic region: For the former, an additional check occurs to determine whether it is in an exon or an intron. Subsequently, distances to the 5′ and 3′ end of each genomic feature (gene/intergenic region/exon/intron) were calculated.

Identification of differentially methylated genes

Using the methylation level of genes at pH 8.0 as a control, GLMs were implemented in R to identify genes that were differently methylated at pHs of 7.2, 7.4, and 7.8, respectively. The general formula used waswhere “methylated, non_methylated” was a two-column response variable denoting the number of methylated and nonmethylated reads at a particular position, while “pH” and “position” were predictor variables for the pH of the environment and the genomic coordinate of the position, respectively. Data from individual replicates were entered separately—this approach assigned equal weightage to each replicate, instead of having a disproportionate skew toward the replicate with the highest coverage if the data were pooled. Genes with <5 methylated positions were filtered out to reduce type I errors.

Identification of differentially methylated promoters

As promoter regions were not defined in the S. pistillata genomes, it was assumed to be located in a window of 4 kb upstream of the transcription start site. A GLM similar to the one described in the previous paragraph was used to identify differentially methylated promoters; however, because of the scarcity of methylated positions in these windows, genes with ≥2 methylated positions were retained for this analysis.

Identification of differentially expressed genes

High-quality total RNA was extracted for library creation from the same S. pistillata nubbins used in the WGBS (that is, also triplicates from four conditions). Directional mRNA libraries were produced using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs), as described previously (). A total of 674 million paired-end reads (read length of 101 bp) were retrieved from six lanes on the HiSeq 2000 platform (Illumina).

Trimming was intentionally left out to increase the number of mapped reads and to reduce bias (). The resulting 362 million trimmed reads were mapped to S. pistillata gene models with kallisto v0.42.4 () to produce TPM (transcripts per million reads) values. On the basis of these values, sleuth () was used to identify differentially expressed genes by contrasting all biological replicates of pHs of 7.2, 7.4, and 7.8 against the controls (pH 8.0).

Calculation of exon expression relative to first exon

As kallisto is a pseudomapper, that is, it assigns reads to gene models but not the exact location of where the read maps to the gene model, reads from all 12 replicates were mapped separately using HISAT2 (v2.1.0) () against the S. pistillata genome. Genomic positions that correspond to exonic locations were extracted with a Perl script, and mean coverages were computed to produce RPKM (reads per kilobase per million mapped reads) values on a per-replicate, per-exon basis with a Python script.

Genes with six or more exons were selected for the analysis. Furthermore, as ratio-based computations are skewed by lowly expressed genes, genes with an overall mean RPKM of >0.5 and nonzero expression in the first six exons across all 12 replicates were selected. Ratios were subsequently natural log–transformed to improve normality, resulting in 1141 unmethylated genes, 5750 methylated genes, and 2955 highly methylated genes (genes with median methylation >80%).

Functional enrichment of methylated genes

GO term annotations were obtained from literature (). GO term enrichment analyses were carried out with topGO () with default settings. GO terms with P < 0.05 and occurring ≥5 times in the background set were considered significant. Multiple-testing correction was not applied on the resulting P values as the tests were considered to be nonindependent ().

KEGG Orthology (KO) annotations were merged from results of gene model annotation and KAAS (KEGG Automatic Annotation Server) (www.genome.jp/tools/kaas/), with parameters “GHOSTZ,” “eukaryotes,” and “bi-directional best hit.” On the basis of the KO annotation, KEGG pathway enrichment was carried out using Fisher’s exact test. Pathways with P < 0.05 were considered significant.

Two KEGG pathways stood out in terms of functionality—“ko04010”: “MAPK signaling pathway” and “ko04110”: “cell cycle.” To generate the pathway diagrams, XML representations of the pathways were imported into Cytoscape (). Genes with no putative homologs in S. pistillata were removed from the pathway; the remaining genes were colored according to their respective changes in methylation level at pH 7.2. To substantiate the functionality of differential methylation in these pathways, a Python script was written to illustrate the general hypermethylation of negative regulators of JNK and MAPK pathways and the general hypomethylation of positive regulators of the same pathways.

RT-qPCR verification of key growth genes

There is only one reference amplicon (against β-actin) designed specifically for this organism (). As normalization against a single reference gene is rarely acceptable (), additional control genes were selected from our RNA-seq data. Selected genes were highly abundant (TPMs > 1000) and expressed at roughly equivalent levels in all sequenced libraries.

Primers were designed, whenever possible, to include a large intronic region in addition to an exonic size of ~100 bp. A Python script that interfaced with Primer3 v2.3.6 () was written to optimize primer design (for example, melting temperatures of ~60°C and GC% of 30% to 70%), avoid long runs of a single base and no strong secondary structure.

A preliminary RT-PCR was run to confirm that amplicons produced band sizes corresponding to the amplified exonic region; as expected, none of the primers produced detectable bands that included the intronic regions. Six reverse transcription reactions (on total RNA from all pH 7.2 and 8.0 replicates) were carried out with SuperScript III First-Strand Synthesis SuperMix (Invitrogen) with the supplied oligo-dT primers. The subsequent qPCR was carried out using Platinum SYBR Green qPCR SuperMix-UDG (Invitrogen) in a 7900HT Fast Real-Time PCR System (Applied Biosystems). Both protocols were carried out according to the manufacturer’s instructions. Primer sequences, amplicon sizes, and the analysis of the RT-qPCR results are fully detailed in data S4.

Laser microdissection of S. pistillata samples

Apexes of colonies were prepared as described previously (). Briefly, apexes of S. pistillata were fixed in 3% paraformaldehyde in S22 buffer [450 mM NaCl, 10 mM KCl, 58 mM MgCl2, 10 mM CaCl2, and 100 mM Hepes (pH 7.8)] at 4°C overnight and then decalcified using 0.5 M EDTA in Ca2+-free S22 at 4°C. They were then dehydrated in an ethanol series and embedded in Paraplast. Cross sections (6 μm thick) were cut and mounted on polyester membrane (0.9 μm) frame slides (Leica Microsystems). The Leica AS LMD system, with a pulsed 337-nm ultraviolet laser on an upright microscope, was used for the microdissections. The laser beam can be moved with a software-controlled mirror system that allows selecting target cells and tissues. Target cells can be preselected on a monitor with a freehand drawing tool, and then the computer-controlled mirror moves the laser beam along the preselected path and the target cells are excised from the section. The dissected part then falls into a PCR tube under gravity. The collection by gravity ensures quick and contamination-free processing of the dissected tissue sections.

DNA from eight sections (duplicates of aboral and oral tissues from pH 7.2 and controls) were extracted and subjected to bisulfite conversion using EpiTect Plus Bisulfite Conversion kit following the manufacturer’s instructions (Qiagen).

Validation of differential methylation via amplicon-specific bisulfite sequencing

As dissected tissues contain minute amounts of DNA, we decided to perform amplicon-specific bisulfite sequencing to validate the methylation levels within our amplicons of interest. A nested PCR design was used to generate these amplicons. Outer primers were pooled in the initial PCR run (35 cycles), and the resulting reaction mix was then evenly split into amplicon-specific individual PCRs (35 cycles) with their respective inner primers. As the amplicons were to be sequenced on the Illumina MiSeq platform, the inner primer pairs were designed with additional overhangs to facilitate downstream library creation, per the 16S Metagenomic Sequencing Library Preparation guide (Illumina).

Amplicons were selected within genes that exhibited differential methylation at pH 7.2 relative to control. As bisulfite conversion produces Watson and Crick strands that are no longer reverse complements of each other, primers were designed to amplify the sense strand of the gene product. To optimize primer design, a self-written Python script that interfaced with Primer3 v2.3.6 () was used to select primer pairs that were fully located within nonmethylated regions. This important consideration avoids the need for degenerate primers (one degenerate base is required per methylated CpG), reducing ordering costs and dodging potential PCR amplification biases. In total, 46 amplicons (that is, 184 primers) were designed.

To assess the primers, test nested PCRs were performed on converted S. pistillata total DNA; PCR products were subsequently visualized on an agarose gel. Primers that completely failed to produce amplicons with expected band sizes were discarded (n = 3). The remaining primers were grouped into three batches (of ~14 primers each) to reduce unintended products that might arise from pooling too many outer primers in the initial reaction.

To validate our original findings and to test the accuracy of this approach, the same amplicon-specific bisulfite approach was carried out on DNA from all three samples of pH 7.2 and control that were used to construct the original whole-genome bisulfite libraries. This allowed for direct comparisons of methylation levels assayed via amplicon-specific bisulfite sequencing and WGBS.

With the use of Nextera XT indices (Illumina), libraries were pooled and sequenced on 1.5 MiSeq runs. A total of 14.2 million paired-end reads (read length of 300 bp) were produced. These reads were cleaned and mapped to the S. pistillata genome with the same pipeline used to process reads from WGBS, with the sole exception of skipping the deduplication step (distinct amplicons from the same locus map to the same genomic coordinates and thus are erroneously deemed duplicates).

We used a very conservative coverage threshold for the downstream analyses: Only methylated positions with read coverages greater than 100 in all samples were retained. This filtering step increased the precision of the measured methylation levels and reduced the effect of noise on the results.

Measurement of cell sizes

As described previously (), branches of 2- to 3-cm size were sampled from colonies grown in pH 7.2 and 8.0 treatments and placed in a 7% MgCl2 solution to anesthetize tissues. Oral discs (the apparent portion of the polyp) were then cut from the colony under a microscope using microdissection scissors with 5-mm blades. Cells were then dissociated using a sterile tube pestle and homogenized by repeated pipetting using a 200-μl pipette. Suspended cells (50 μl) were mounted between the slide and coverslip, and 30 random pictures of the surface were taken. Cell sizes on these pictures were subsequently measured using SAISAM software (Microvision Instruments). In all, 4728 measurements were taken across 14 tested fragments (data S5A).

Measurement of corallite calyx size

Branches of similar size (2 cm in length) were sampled from colonies in the pH 7.2 and 8.0 treatments and placed in a 10% NaClO solution for 2 hours to remove tissues. Skeletons were rinsed several times in tap water, followed by ultrapure H2O, and then dried at 40°C for at least 24 hours. Samples (five replicates per pH condition) were coated with gold/palladium and observed at 4 kV with a JEOL 6010LV electron microscope. Diameters of the corallite calyxes were measured using manufacturer-provided SMILE VIEW software (JEOL). A total of 181 measurements were performed across 10 fragments (data S5B).

Measurement of skeletal porosity

Five coral fragments were analyzed per treatment (pH 7.2 and 8.0; data S5C). The noninvasive, high-resolution imaging method of x-ray microcomputed tomography (micro-CT) was used to measure skeletal porosity (), as described. Micro-CT analysis was carried out at the Polyclinique Saint Jean, Cagnes-sur-Mer, France, with a SkyScan 1173 compact micro-CT (SkyScan). A microfocus x-ray tube with a focal spot of 10 mm was used as a source (80 kV and 100 μA). The sample was rotated 360° between the x-ray source and the camera. Rotation steps of 1.5° were used. At each angle, an x-ray exposure was recorded on the distortion-free flat-panel sensor (resolution, 2240 × 2240 pixels). The resulting slice was made of voxels, the three-dimensional equivalent of pixels. Each voxel was assigned a gray value derived from a linear attenuation coefficient that relates to the density of materials being scanned. All specimens were scanned at the same voxel size. The radial projections were reconstructed into a three-dimensional matrix of isotropic voxels ranging from 5 to 10 mm, depending on the exact height of the coral tip. X-ray images were transformed by NRecon software (SkyScan) to reconstruct two-dimensional images for quantitative analysis. From these images, evaluation of the morphometric parameters was performed using CT Analyzer software (SkyScan). A manual grayscale threshold was implemented on the first set of images and then applied to all specimens. For each sample, a digital region of interest was created to extend through 100 μm of skeleton at 7-mm distance from the apex, corresponding to about 15 slices. The percentage of negative space in relation to the skeleton was then determined, providing the measure of porosity.