Phenotypic plasticity can buffer organisms against short-term environmental fluctuations. For example, previous exposure to increased temperatures can increase thermal tolerance in many species. Previous studies have found that this acclimation response is associated with changes in the magnitude of the gene expression response (hereafter, “transcriptional response modulation”). However, mechanisms mediating this gene expression response and, ultimately, phenotypic plasticity remain largely unknown. Epigenetic modifications are good candidates for modulating transcriptional response, as they broadly correlate with gene expression. Here, we investigate changes in DNA methylation as a possible mechanism controlling shifts in gene expression plasticity and thermal acclimation in the reef-building coral Acropora nana. We find that gene expression response to acute stress is altered in corals acclimated to different temperatures, with many genes exhibiting a dampened response to heat stress in corals pre-conditioned to higher temperatures. At the same time, we observe shifts in methylation during both acclimation (11 days) and acute heat stress (24 hours). We observed that the acute heat stress results in shifts in gene-level methylation and elicits an acute transcriptional response in distinct gene sets. Further, acclimation-induced shifts in gene expression plasticity and differential methylation largely occur in separate sets of genes. We found no overall correlation between the magnitude of differential methylation and the change in gene expression plasticity. We do find a small but statistically significant number of genes exhibiting dampened expression response and shifts in methylation (14 genes), which could be candidates for further inquiry. Overall, our results suggest transcriptional response modulation occurs independently of methylation changes induced by thermal acclimation.  

The samples used in this study were taken from a previously published acclimation experiment, and details on the experimental setup can be found in that paper (Fig 1A; Bay & Palumbi, 2015). Briefly, we took whole small colonies of Acropora nana (Studer, 1878)** from the reef on Ofu Island, American Samoa, and placed them in outdoor aquaria for acclimation. Six colonies were placed in each tank, with two tanks per acclimation treatment, totaling 12 colonies per acclimation treatment. Due to logistical constraints, genotype was not replicated across acclimation treatments. The original study had an ambient (29°C), elevated (31°C), and variable (29-33°C daily) acclimation treatment to mimic in situ thermal fluctuation. We only used samples from the two stable treatments (29°C and 31°C) for our purposes. This 2°C temperature difference was sufficient to induce thermal acclimation and is within the typical thermal range. At different time points throughout the experiment, we sampled branches from each colony and subjected them to an acute heat stress assay to test thermal tolerance. For the heat stress assay, two branches were sampled from each colony; one was held at an ambient control temperature (29°C) for 24 hours while the other was subjected to heat stress consisting of a three-hour ramp to 34°C followed by five hours at 34°C, then a decline to 29°C for approximately one hour. The heated samples were incubated at 29°C overnight for the remainder of the 24-hour assay. At 6:00 AM the following morning, branches were cut in half and preserved in 95% ethanol and RNA stabilizing solution (70g Ammonium Sulfate /100ml solution, 10 mM EDTA, 25 mM Sodium Citrate, 5.4 pH) for downstream analysis. To assess thermal tolerance after heat stress, we measured chlorophyll concentration as a proxy for bleaching (Ritchie, 2008). Figure 1B, redrawn from Bay & Palumbi (2015), shows the increase in thermal tolerance (reflected in higher chlorophyll content of heat-stressed samples) for 31°C acclimated corals than 29°C acclimated corals after 11 days. In this study, we examine samples in the acclimation treatments for either 0 days (i.e., collected the previous day and held at 29°C overnight) or 11 days. Note that this means the Day 0 samples in the 31°C tank never actually experienced 31°C acclimation. 

RNASeq

Six samples per treatment group, totaling 48 samples, were extracted and sequenced for RNASeq. Sample preservation, library prep, cDNA sequencing, and read trimming are described in Bay & Palumbi 2015. We aligned trimmed reads to the Acropora millepora genome using STAR aligner software (v. 2.7.0e) (Dobin et al. 2013; Fuller et al. 2020). Filtering parameters that optimized alignment were: --outFilterScoreMinOverLread 0, --outFilterMatchNminOverLread 0, --outFilterMatchNmin 0, --outFilterMismatchNmax 4. We used HTSeq v. 0.9.1 to count all reads that mapped to genes, including reads that map to multiple genes(--nonunique all), using the annotated gene models provided with the Acropora millepora reference genome (Anders et al., 2015; Fuller et al. 2020). We expect that a congener alignment may result in some reads failing to map due to the divergence between Acropora millepora and Acropora nana, thus, these reads are excluded from the remainder of the analyses.

To quantify the symbiont composition within each sample, we aligned RNASeq reads to Symbiodinium goreaui (GenBank accession number: AF333515) and Durusdinium trenchii (GenBank accession number: LC718590) ITS2 sequences using the STAR aligner software (v. 2.7.0e) (Dobin et al. 2013). We used the STAR option --outFilterMultimapNmax 1 to ensure that STAR reports only the best alignment for each read. We counted the number of reads that uniquely aligned to each ITS2 sequence in each sample.

Whole-genome bisulfite sequencing

We extracted DNA from 6 samples per treatment group, totaling 48 samples, using the Qiagen DNeasy Blood & Tissue Kit (Cat. No. 69504). Genomic DNA was sent to Novogene (Sacramento), and Methyl-MaxiSeq libraries were prepared from 300 ng of genomic DNA digested with 2 units of Zymo Research’s dsDNA Shearase™ Plus (Cat. No. E2018-50). The fragments produced were end-blunted, 3’-terminal-A extended, then purified using the Zymo Research DNA Clean & Concentrator™ kit (Cat. No. D4003). The A-tailed fragments were ligated to pre-annealed adapters containing 5’-methylcytosine instead of cytosine. Bisulfite treatment of the fragments was done using the EZ DNA Methylation–Lightning kit (Zymo Research, Cat. No. D5030). PCR was performed with Illumina TruSeq indices, and the size and concentration of the fragments were confirmed on the Agilent 2200 TapeStation. Before sequencing, samples were spiked with Illumina's PhiX Control library. Sequencing was performed using the Novaseq 6000 platform with Paired-End 150 (PE150) reads, aiming for a target coverage of 30x based on the 500 Mb Acropora millepora estimated genome size. The average raw read coverage achieved was 23x.

To analyze WGBS reads, we used Trim Galore! Version 0.6.3 to filter out reads that were less than 20 nt long along with their read-mate (Krueger 2019). We clipped the 5’ ends of the reads to remove possible methylation bias (--clip_R1 10 --clip_R2 10). We aligned trimmed reads to the Acropora millepora genome using bwa-meth with default settings (Fuller et al. 2020; Pedersen et al. 2014). We used SAMtools view (Version 1.9) to filter by excluding unmapped reads, mapped reads with unmapped mates, and/or reads that failed platform quality thresholds (-F 524) (Li et al. 2009). We kept reads with a minimum mapping quality of 2 (-q 2) (Li et al. 2009). To filter out PCR duplicates, we used Picard MarkDuplicates (Version 2.20.2) with default settings (Picard). Finally, to extract methylation calls, we used MethylDackel (https://github.com/dpryan79/MethylDackel). This output summarizes the strand-specific frequency of cytosines and thymines (a proxy for unmethylated cytosine). In downstream statistical analyses, these frequencies are converted to a percentage of DNA methylation at each CpG dinucleotide. We did not mask possible C to T substitutions between the reference genome and Acropora nana samples as collating these transitions were not feasible with our datasets. The minimum depth of methylation calls was set to 10, and the maximum variant fraction was set to 0.75 to exclude possible non-cytosine alleles at reference CpG sites. 

Differential gene expression and heat response gene identification

We analyzed differential gene expression using R Package DESeq2 v. 1.36.0 (Love et al. 2014, R Core Team 2022). First, we filtered for genes with a depth of at least 10 reads in all samples and normalized data with the variance stabilizing transformation. We performed a PCA to investigate the relationship between acclimation and heat stress treatment groups. We identified heat response genes by comparing the heat-stressed samples to their corresponding control counterparts that had not undergone thermal acclimation (29°C  at Day 0 and Day 11, and 31°C at Day 0). Since thermal acclimation can lead to the downregulation of a subset of stress response genes, transcriptional data from thermally acclimated samples were excluded from identifying the heat response genes. Heat response genes are defined by a |log₂^^ fold-change| > 2 and Benjamini-Hochberg false discovery rate (BH p_adj) < 0.01. We classified a total set of heat response genes from the union of the significant differential expression within any of the three sample sets: 29°C (Day 0 and Day 11) and 31°C (Day 0). A more conservative set of “core” heat stress genes were identified by the intersection of the significant heat stress genes across all three contrasts. We analyzed Gene Ontology (GO) terms with the R package, topGO v. 2.50.0 (Alexa and Rahnenfuhrer 2022). We compared the list of “core” heat response genes against a background of all genes that passed our quality filters. Enriched GO terms were identified with the classic Fisher’s test with a p-value < 0.01 and at least 10 transcripts within each category.

Quantifying transcriptional response modulation

Previously, we found no change in gene expression during acclimation but rather a change in the magnitude of gene expression response to heat stress in corals pre-conditioned at higher temperatures (Bay and Palumbi 2015). Here, we quantify this “transcriptional response modulation” using R Package DESeq2 v. 1.36.0 (Love et al. 2014, R Core Team 2022). DESeq2 estimates the interaction term coefficient for each gene and uses the Wald test to test if the interaction term is statistically different from zero (Love et al. 2014). The predictors of the coefficient terms are the levels of the experimental design (Love et al. 2014). In our case, the experimental conditions are acclimation temperature, heat stress treatment, and the interaction between these conditions.To avoid the confounding factor of Day, we used heat stress treatment and control samples from the acclimation temperature control, 29°C (Day 11), and the acclimation treatment, 31°C (Day 11), for comparison to capture the effect of thermal acclimation. We identified two categories of transcriptional response modulation: amplified and dampened expression. Transcripts with amplified expression had a higher magnitude of expression response induced by heat stress in samples acclimated to 31°C than those acclimated to 29°C (Wald statistic > 0 and BH p_adj < 0.1). Transcripts with dampened expression had a smaller magnitude of expression response (Wald statistic < 0 and BH p_adj < 0.1) (Love et al. 2014). To test if amplified and dampened genes have different mean expressions, we log₁₀-transformed each transcript's average normalized count values and performed a two-sample t-test. We analyzed Gene Ontology (GO) terms with the R package, topGO v. 2.50.0 (Alexa and Rahnenfuhrer 2022). We compared lists of amplified and dampened genes against a background of all genes that passed our quality filters. Enriched GO terms were identified with the classic Fisher’s test with a p-value < 0.01 and at least 10 transcripts within each category. 

Effects of treatment on DNA methylation

We used methylKit to analyze CpG methylation calls in R (Akalin et al. 2012, R Core Team 2022). We maintained CpG sites in the ‘methylRawDB’ if there was a coverage count minimum of at least 10, had a maximum cut-off in the 99.9 percentile of read counts covering the site, and was covered in all the samples. We performed principal coordinate analysis (PCoA) using the Bray-Curtis dissimilarity of DNA methylation sites to visualize the relationship between acclimation and heat stress treatment groups. We recalculated the PCoA after removing the following four outliers: 29°C acclimation treatment (Day 11) heat stress sample replicate 1, 29°C acclimation treatment (Day 11) heat stress control replicate 1, 31°C acclimation treatment (Day 11), heat stress treatment replicate 1 and 31°C acclimation treatment (Day 11), heat stress treatment replicate 3.

We assessed the distribution of methylation across CpG sites and genes in our data set. We used methylKit to summarize percent methylation across CpG sites and genes (Akalin et al. 2012). CpG sites were included in the ‘methylRawDB’ if there was a coverage count minimum of at least 10 and had a maximum cut-off in the 99.9 percentile of read counts covering the site. For the CpG site methylation distribution, we used the methylKit ‘unite’ function to collate the methylation data and calculate the methylation percentage where the site was covered in at least 3 replicates per treatment from distinct colonies. For the gene-level methylation distribution, we created a GRanges object by importing annotated A. millepora gene models furnished with the reference genome (Lawrence et al. 2009; Fuller et al. 2020). Then we used the methylKit ‘unite’ function to collate the methylation data, integrate the methylation data with the GRanges object, and calculate the methylation percentage where the gene was covered in at least 3 replicates per treatment.

Previous studies have found that environmental shifts can alter global methylation levels (Putnam et al. 2016; Metzger and Schulte 2017). We conducted distinct tests to examine heat stress and acclimation induced effects within our dataset. We combined methylation data across all samples with the methylKit ‘unite’ function (Lawrence et al. 2009; Akalin et al. 2012; Fuller et al. 2020). We calculated the average percent of methylation at each CpG site covered in at least 3 replicates per treatment across all treatment groups. We excluded missing data from the average methylation percent calculation. We tested the effect of heat stress and acclimation treatments on global CpG methylation percent with an ANOVA.

We tested whether DNA methylation percent varies between genomic features and if there were feature-specific changes in methylation following acclimation. First, we were interested in the methylation percentage within and between genomic features. We created a GRanges object by importing annotated A. millepora exons and introns (Lawrence et al. 2009; Fuller et al. 2020). We created a GRanges object by importing BED files that define predicted coordinates of promoters, transcription start sites (TSS), long interspersed nuclear element (LINE) repeats, short interspersed nuclear element (SINE) repeats, and rolling circle (RC) repeats in the A. millepora genome (https://github.com/Groves-Dixon-Matz-laboratory/benchmarking_coral_methylation/tree/master/windowStats) (Dixon and Matz 2021). We integrated the feature-specific GRanges objects with methylation data across all samples and calculated the methylation percentage where the region was covered in at least 3 replicates per treatment with the methylKit ‘unite’ function  (Akalin et al. 2012). We calculated the average percent of methylation at each promoter, TSS, exon, intron, LINE repeat, SINE repeat, and RC repeat across all samples. We excluded missing data from the average methylation percent calculation. We used an ANOVA followed by a Tukey’s HSD test to test whether DNA methylation varies across genomic features. 

Next, we investigated gene-level differential methylation due to heat stress. We created a GRanges object by importing annotated A. millepora genes and then integrated this object with methylation data across all samples with the methylKit ‘unite’ function (Lawrence et al. 2009; Akalin et al. 2012; Fuller et al. 2020). We set ‘min.per.group’ to 3 and performed gene-level differential methylation analysis using methylKit’s ‘calculateDiffMeth’ function to determine genes with different methylation between heat-stressed and control groups (Akalin et al. 2012). The heat assay and control groups were formed by pooling all heat-stress and control samples, respectively, across all acclimation treatment groups: 29°C (Day 0), 31°C (Day 0), 29°C (Day 11), and 31°C (Day 11). We set the minimum percentage change threshold to 25% difference between heat-stressed and control groups. The results were corrected for false discovery rate using a q-value threshold of 0.05. As we did not account for C-to-T substitutions in our data, these substitutions may artificially increase unmethylated calls at CpG sites, potentially resulting in decreased variance in methylation. After calculating the change in DNA methylation between heat stress treatments, we analyzed Gene Ontology (GO) terms with the R package, topGO v. 2.50.0 (Alexa and Rahnenfuhrer 2022). We compared differentially methylated genes against a background of all genes that passed our quality filters. Enriched GO terms were identified with the classic Fisher’s test with a p-value < 0.01 and at least 10 genes within each category.

We hypothesize that thermal acclimation may lead to changes in DNA methylation, which forms the primary focus of our study. We investigated differences in methylation due to acclimation at the CpG site and gene levels. To capture the effect of thermal acclimation, we formed the acclimation control group by pooling 29°C (Day 11) heat stress treatment and control samples, and we formed the acclimation treatment group by pooling 31°C (Day 11) heat stress treatment and control samples. We used the same filtering criteria as when assaying differences in methylation due to heat stress. After calculating the change in DNA methylation between acclimation and controls, we analyzed GO terms as explained above.

Relationship between gene expression and DNA methylation 

First, we tested the relationship between gene expression and average methylation percent within each coding region by performing a linear regression between the log₁₀ transformed means of expression counts and the average percent methylation of genes in all samples. To test this relationship, we used the linear regression function in base R where the base mean of expression for each gene was the response variable and the average percent of gene-level methylation was the predictor variable (R Core Team 2022). We used ANOVA to determine effect strength and significance.

Then, we tested the relationship between heat stress gene expression plasticity and average methylation percent within each coding region. We used R Package DESeq2 v. 1.36.0 to calculate the log₂ fold-change between all heat-stressed and control samples and to perform log₂ fold-change shrinkage with the ‘ashr’ method to account for the strong variation in log₂ fold-change associated regions of low read counts (Love et al. 2014; Stephens 2017). All heat-stressed and control samples were collated with gene-level methylation percent data. We then performed a linear regression between the gene-level posterior standard deviation of shrunken log₂ fold-changes and percent methylation. We used the linear regression function in base R (R Core Team 2022). The standard deviation was the dependent variable, and the percent methylation was the independent variable. We used ANOVA to determine effect strength and significance.

Relationship between transcriptional response modulation and change in gene level methylation

Our central hypothesis was that shifts in DNA methylation during acclimation result in transcriptional response modulation. In other words, genes with differential methylation between acclimation treatments would be those whose expression response to acute heat stress was either dampened or amplified in 31°C acclimated samples compared to those acclimated at 29°C. We merged acclimation-induced differential methylation data with the analysis of transcriptional response modulation (see above). We collated the measure of transcriptional response modulation–the interaction term coefficient and difference in methylation percent for each gene. We tested the relationship between transcriptional response modulation and change in the gene-level percent DNA methylation using the linear regression function in base R, where for each gene, the measure of transcriptional response modulation was the dependent variable and methylation difference was the independent variable (R Core Team 2022). We used ANOVA to determine effect strength and significance. Finally, we tested the difference in average percent methylation of amplified and dampened genes using a two-sample unpaired wilcoxon test.