What drives phenotypic divergence among coral clonemates of Acropora palmata?

Evolutionary rescue of populations depends on their ability to produce phenotypic variation that is heritable and adaptive. DNA mutations are the best understood mechanisms to create phenotypic variation, but other, less well-studied mechanisms exist. Marine benthic foundation species provide opportunities to study these mechanisms because many are dominated by isogenic stands produced through asexual reproduction. For example, Caribbean acroporid corals are long lived and reproduce asexually via breakage of branches. Fragmentation is often the dominant mode of local population maintenance. Thus, large genets with many ramets (colonies) are common. Here, we observed phenotypic variation in stress responses within genets following the coral bleaching events in 2014 and 2015 caused by high water temperatures. This was not due to genetic variation in their symbiotic dinoflagellates (Symbiodinium "fitti") because each genet of this coral species typically harbours a single strain of S. "fitti". Characterization of the microbiome via 16S tag sequencing correlated the abundance of only two microbiome members (Tepidiphilus, Endozoicomonas) with a bleaching response. Epigenetic changes were significantly correlated with the host's genetic background, the location of the sampled polyps within the colonies (e.g., branch vs. base of colony), and differences in the colonies' condition during the bleaching event. We conclude that long-term microenvironmental differences led to changes in the way the ramets methylated their genomes, contributing to the differential bleaching response. However, most of the variation in differential bleaching response among clonemates of Acropora palmata remains unexplained. This research provides novel data and hypotheses to help understand intragenet variability in stress phenotypes of sessile marine species.

INTRODUCTION

Acclimatization is a nongenetic process by which an individual can increase its stress tolerance after exposure to a stressor, such as temperature anomalies (see Palumbi, Barshis, Traylor‐Knowles, & Bay, 2014 for a coral example). Acclimatization can lead to phenotypic variability in stress response among clonemates. Nonmutation‐based mechanisms resulting in phenotypic variability in isogenic lines include stochastic gene expression, errors in protein synthesis, protein promiscuity and epigenetic modifications (reviewed by Payne & Wagner, 2019). Of these, only epigenetic mutations have been studied as a mechanism for acclimatization in marine foundation fauna (reviewed by Eirin‐Lopez & Putnam, 2019). Interestingly, nonmutation‐based changes to phenotypes can prolong survival of a genet until such phenotypes become permanent (Yanagida et al., 2015), although these processes are not well understood. Understanding of all the mechanisms that produce phenotypic variability is essential to estimate the evolvability of threatened marine species (Payne & Wagner, 2019).

The large‐scale bleaching event during 2014 and 2015 within the Florida Keys provided an unprecedented opportunity to understand the role of acclimatization and phenotypic variability in framework‐building foundation species of shallow Caribbean coral reefs. Because reef‐building corals harbour intracellular algal symbionts (family Symbiodiniaceae), discerning the relative contribution of host and symbiont to holobiont acclimatization can be difficult. However, the Caribbean elkhorn coral, Acropora palmata, has an uncomplicated symbiosis: it associates with just one symbiont species (Symbiodinium “fitti”) and most colonies (~70%) harbour only one strain of S. “fitti” over space and time (Baums, Devlin‐Durante, & LaJeunesse, 2014) similar to other cnidarian–Symbiodiniaceae associations (e.g., Andras, Kirk, & Drew Harvell, 2011; Thornhill, Xiang, Fitt, & Santos, 2009). Thus, A. palmata is an excellent model to discern the coral host's acclimatization response and phenotypic variability to heat stress.

Like many reef‐building corals, A. palmata frequently reproduces via fragmentation (an asexual process), sometimes forming large, monoclonal stands (Baums, Miller, & Hellberg, 2006; Foster, Baums, & Mumby, 2007; Pinzón, Reyes‐Bonilla, Baums, & LaJeunesse, 2012; Williams, Miller, & Baums, 2014). These colonies represent iterations of the same host–symbiont combination (i.e., they are isogenic replicates), experiencing similar environmental conditions. Initial surveys of A. palmata during the 2014 and 2015 bleaching events documented a range of bleaching responses. This response varied between reefs but also within single, monoclonal stands of A. palmata (see Figure 2). Thus, coral clonemates exhibited different bleaching susceptibilities despite data showing that they share identical (clonal) S. “fitti” symbiont communities, begging the question as to what mechanisms account for such phenotypic variability. Thus, we explored whether differences in the microbiome other than dinoflagellates (e.g., prokaryotes), and/or micro‐environmental differences, such as shading or exposure to water movement, induced epigenetic changes in the host genome that could explain differences in bleaching susceptibilities among ramets of the same genet. The answer may inform our understanding of how reefs might survive climate change and has implications when choosing genets for restoration.

Figure 2

Heterogeneity in bleaching phenotype was observed among ramets (colonies) of the same Acropora palmata genet. (a) Genet at Molasses Reef. (b) Ramets of genet P1002 at Sand Island reef photographed during peak bleaching in September of 2014 and again in February 2015. Genet P1002 harbored a single strain of S.  ‘fitti’ in all ramets. Images suggested that the heterogeneity in bleaching was due to micro‐environmental differences that may have altered the underlying epigenome [Colour figure can be viewed at http://wileyonlinelibrary.com]

Corals associate with a diverse set of prokaryotes (Pollock et al., 2018). The microbiome plays important roles in the functioning of the holobiont, including coral nutrition, element cycling and disease responses (Peixoto, Rosado, Leite, Rosado, & Bourne, 2017). This community can shift in response to stressors such as heat, low pH or disease (Muller, Bartles, & Baums, 2018; Thurber et al., 2009). We are beginning to understand the specific function of some members of the microbiome such as an endosymbiotic cyanobacterium that fixes nitrogen in Montastraea cavernosa (Lesser et al., 2007) and the disease‐causing Vibrio corallyticus (Kimes et al., 2012). Other species, such as Endozoicomonas sp., show patterns of association with coral hosts that suggest an important role, but that role is not yet well understood (Neave et al., 2017; Pollock et al., 2018).

DNA methylation is the most highly studied mechanism in epigenetics and is often used to elucidate how a phenotype is modified without altering the genetic code. DNA methylation occurs at the cytosine bases of eukaryotic DNA, which are converted to 5‐methylcytosine by DNA methyltransferase (DNMT) enzymes. Additional epigenetic changes include histone modifications, chromatin remodelling and gene regulatory mechanisms involving small noncoding RNAs (Danchin et al., 2011). Epigenetic modifications can rapidly produce new phenotypes in response to a change in the environment without mutations in the underlying genetic sequence (Finnegan, 2002; Richards, 2008).

There are limited data on methylation mechanisms in invertebrates. Current evidence has shown that invertebrate genomes are far less methylated than vertebrate genomes (Gavery & Roberts, 2010, 2013; Lyko et al., 2010; Olson & Roberts, 2014; Rivière, 2014; Suzuki, Kerr, Sousa, & Bird, 2007). Additionally, DNA methylation is predominately found in gene bodies in which highly expressed housekeeping genes are hypermethylated and regulated and/or inducible genes are hypomethylated (Dimond & Roberts, 2016; Elango, Hunt, Goodisman, & Soojin, 2009; Gavery & Roberts, 2010; Hunt, Brisson, Soojin, & Goodisman, 2010; Sarda, Zeng, Hunt, & Soojin, 2012). Hypermethylation of those genes essential for biological function is thought to imply that they are “protected” from plasticity in transcriptional opportunities. Such plasticity would be inherently lethal in housekeeping genes (Roberts & Gavery, 2012), and thereby gene body methylation (GBM) in corals is correlated with stable and active transcription (Dixon, Liao, Bay, & Matz, 2018; Liew, Zoccola, et al., 2018). Methylated genes in the anemone Aiptasia show a significant reduction of spurious transcription and transcriptional noise (Li et al., 2018). In contrast, the inducible gene's limited methylation may facilitate, albeit passively, specific transcriptional opportunities including increasing sequence mutations, access to alternative transcription start sites and exon skipping (Roberts & Gavery, 2012). There is evidence for a direct relationship between DNA methylation and phenotypic plasticity, as seen in the determination of caste structure in both honeybees and ants (Kucharski, Maleszka, Foret, & Maleszka, 2008), maternally inherited epigenetic patterns that influence the expression of the agouti gene in mice (Wolff, Kodell, Moore, & Cooney, 1998), and the differential methylation of the human NR3C1 gene in newborns depending on prenatal maternal mood (Oberlander et al., 2008). Scleractinian corals can display strong differences in their DNA methylation in response to stress, demonstrating that de novo DNA methylation may be a driving mechanism for phenotypic plasticity in acclimatization (Putnam, Davidson, & Gates, 2016).

Here, we sampled six different genets from four different reef sites in Florida (Figure 1) and sampled each genet six times (either across the same colony six times, or different ramets of the same genet for a total of six samples). Samples were taken from the upward‐facing side of branches or the bases/trunks of colonies 6 weeks after colonies had experienced a thermal stress event that resulted in differential bleaching. We determined how many S. “fitti” strains were present in each sample and sequenced the 16S gene to characterize variability in the prokaryotic community. We then applied a reduced sequencing technique sensitive to the methylation status of cytosine called MethylRad (Wang et al., 2015) to identify sites that were differentially methylated between reefs, genets, position within the colony and peak bleaching status.

Figure 1

Acropora palmata samples were obtained from four reef sites within the Florida Reef tract, Grecian Rocks, Sand Island, Molasses and Key Largo Dry Rocks. Distance between Grecian Rocks and Key Largo Dry Rocks is about 20 km. Maps created in Google Earth Pro [Colour figure can be viewed at http://wileyonlinelibrary.com]

METHODS

Sample collection

During routine monitoring of established Acropora palmata study sites (Williams et al., 2014) in September 2014 severe bleaching was observed. Colony condition, including a ranking of bleaching severity, was recorded (Williams, Miller, Bright, Pausch, & Valdivia,2017) and detailed photos were taken of each colony. Bleaching severity ranks ranged from 0 indicating normal coloured tissue to 5 for severe bleaching indicating intact tissue that was lacking any colour (appeared white). These ranks were applied based on the overall bleaching severity of the colony as a whole in the field, and ranks were applied by a single observer (D.E.W.). Bleaching severity was monitored in October and again in November, and routine quarterly monitoring was resumed in 2015 and continued through the 2015 bleaching event.

During November 2014 surveys, multiple tissue samples were collected from several precise locations on identified colonies. The tissue collection locations were selected based on their bleaching condition, which varied within and between colonies. Collection sites included four reef sites within about 20 km on the Florida reef tract (Figure 1). We sampled six different genets from these four reef sites in Florida and sampled each genet six times (either across the same colony six times, or different ramets of the same genet for a total of six samples; Table 1). Visual assessment of the tissue condition of the sample spot was scored as normal colour (NC) or slightly pale (SP), very pale (VP) or bleached (BL). The location of the sample within the colony was classified as either exposed branch (high irradiance) or nonbranch (base, lower irradiance). Tissue was collected using a small chisel or snippers and the tissue sample was immediately preserved underwater (~30 s after collection) using RNAlater (Sigma Aldrich).

Table 1

Acropora palmata colonies included in the MethylRAD analysis

Reef	Coral genet	Database ID	Ramet	Field ID	Peak bleaching	Location of within‐colony sampling	Symbiont genet
Grecian Rocks	P1034	12,648	1	GR2‐t844A	BL	N	F109
		12,649	1	GR2‐t844B	BL	B	F109
		12,650	1	GR2‐t844C	NC	B	F109
		12,652	1	GR2‐t844E	SP	N	F109
		12,653	2	GR2‐t846A	BL	B	F447
		12,659	3	GR2‐t847B	NC	N	F109
Grecian Rocks	P2656	12,664	1	GR3‐t398A	BL	B	F476
		12,667	1	GR3‐t398D	NC	N	F476
		12,669	1	GR3‐t398F	NC	N	F476
		12,670	1	GR3‐t398G	BL	B	F476
		12,676	2	GR3‐t873B	SP	B	F477
		12,681	2	GR3‐t873G	BL	B	F477
Key Largo DR	P2598	12,709	1	KL4‐t852A	BL	B	F478
		12,710	1	KL4‐t852B	NC	B	F478
		12,711	1	KL4‐t852C	NC	B	F478
		12,712	1	KL4‐t852D	VP	B	F478
		12,723	2	KL4‐t858A	BL	B	F478
		12,724	2	KL4‐t858B	NC	B	F478
Key Largo DR	P2138	12,715	1	KL4‐t853A	VP	B	F242
		12,717	1	KL4‐t853C	VP	B	F242
		12,718	1	KL4‐t853D	NC	B	F242
		12,719	1	KL4‐t853E	NC	B	F242
		12,720	1	KL4‐t853F	NC	B	F242
		12,721	1	KL4‐t853G	BL	B	F242
Molasses	P2151	12,782	1	ML3‐t915D	NC	N	F448
		12,784	1	ML3‐t915F	NC	N	F448
		12,785	2	ML3‐t955A	VP	B	F448
		12,787	2	ML3‐t955C	VP	B	F448
		12,789	2	ML3‐t955E	NC	B	F448
		12,793	3	ML3‐x01C	BL	B	F448
Sand Island	P1002	12,808	1	SI1‐t959C	BL	B	F449
		12,810	1	SI1‐t959E	NC	B	F449
		12,811	1	SI1‐t959F	VP	N	F449
		12,812	1	SI1‐t959G	NC	N	F449
		12,813	1	SI1‐t959H	BL	N	F449
		12,814	1	SI1‐t959I	NC	B	F449

Samples from four reefs in Florida contained six different coral genets and eight Symbiodinium fitti genets. Phenotypic data included bleaching condition (BL = bleached, VP = very pale, NC = normal colour, SP = slightly pale), and sampling location within a colony (Branch = B, not a branch = N). Key Largo DR = Key Largo Dry Rocks.

Symbiodinium microsatellite analysis

Symbiodinium “fitti” DNA was amplified with primers given by Pinzón, Devlin‐Durante, Weber, Baums, and LaJeunesse (2011) (n = 10 loci). Three additional loci for S. “fitti” were developed de novo (Table S1). Alleles were fluorescently visualized and sized with internal standards on a PRISM 3100 Genetic Analyzer (Applied Biosystems). Briefly, 25–50 ng of template DNA was added to polymerase chain reactions (PCRs) containing 1 × Standard Taq Buffer (New England Biolabs), 2.5 mm MgCl (New England Biolabs), 0.5 mg/ml bovine serum albumin (New England Biolabs) and 0.75 U of Taq (0.325 U for primers 31, 32 and 41; New England Biolabs). Each primer was added at 200 nm to reactions involving loci A3Sym_01, 18, 27 and 28, whereas a primer concentration of 93 nm was added to reactions for primers 31, 32 and 41. Concentrations of 50 nm of the tailed forward primer (see Table S1), 150 nm of the reverse primer and 75 nm of the dye‐labelled T‐oligonucleotide were used for amplifying loci A3Sym_01, 02, 03, 07, 08, 09 and 48. All loci were amplified using the following thermal cycle profile: 94°C for 2 min (one cycle); 94°C for 15 s, primer‐specific annealing temperature (Table S1) for 15 s and 72°C for 30 s (31 cycles); and 72°C for 30 min on a Mastercycler gradient thermal cycler (Eppendorf). Table S2 provides details of microsatellite genotype allele calls.

16S library preparation and bioinformatics

The hypervariable V1 and V2 regions of the 16S rRNA gene were amplified using a general bacterial primer pair with Illumina CS1 and CS2 Fluidigm adapters 27F 5′‐TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGAGAGTTTGATCMTGGCTCAG‐3′ and 355R 5′‐GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGCTGCCTCCCGTAGGAGT‐3′ (Rodriguez‐Lanetty, Granados‐Cifuentes, Barberan, Bellantuono, & Bastidas, 2013). PCRs were performed using 2 U Gotaq (Promega), 1 × Gotaq buffer (Promega), 0.25 mm dNTPs (Bioline), 2.5 mm MgCl2 and 0.25 µm of each primer. Initial denaturation was performed at 95°C for 5 min followed by 30 cycles of 95°C for 30 s, 51°C for 1 min and 72°C for 1 min. Final extension was performed at 72°C for 7 min. PCR products were checked on a 1% agarose gel for a band near 450 bp. Libraries were prepared at the DNA services facility at the University of Illinois and sequenced on a MISEQ2500 platform. Genets P1034, P2138, P1002 and P2151 were run on one chip, and genets P2656 and P2598 were run together on the second Miseq chip.

Paired‐end sample data with quality scores were imported and then reads were joined using vsearch join‐pairs and quality filtered using quality‐filter q‐score‐joined within qiime2 (Bolyen et al., 2018; Caporaso et al., 2010). Chimeras were removed and sequences were denoised using deblur with a p‐trim‐length parameter of 230, which truncates the read at this position. Illumina runs were analysed separately in qiime2 and combined after the deblur analysis. The Naive Bayes classifier was trained on only the region of the target sequences that was sequenced and taxonomy was assigned using silva 132 reference sequences, clustered at 99% sequence similarity. The output was a table containing operational taxonomic units (OTUs). OTUs were normalized using the qiime2 plugin q2‐perc‐norm (Gibbons, Duvallet, & Alm, 2018), a model‐free approach that corrects batch effects in case‐control microbiome studies, with condition of bleaching (normal colour/slightly pale (Not bleached) versus bleached/very pale (Bleached)) as the case control. Correction for batch effects was necessary because genet was confounded in Miseq runs.

Differences in microbiome community composition among genet, reef or bleaching condition were tested using permutational multivariate analysis of variance (PERMANOVA). Bray–Curtis dissimilarity matrices were created using the vegdist() function (method = “bray”) and PERMANOVA tests (999 permutations) were done using the adonis() function in the vegan package (Dixon, 2009). The R vignette variancepartition (Hoffman & Schadt, 2016) was then used to fit a linear mixed model for each order/genus of bacteria and partition the total variance into the fraction attributable to each phenotype in the design, including genet, location of sampling (branch vs. not branch), and the condition of the sample (bleached vs. not bleached, Table S3A).

Methylrad library preparation and bioinformatics

Coral tissue samples were snap frozen in liquid nitrogen and kept at −80°C before extraction using the DNeasy Blood & Tissue Kit (Qiagen) with the following modifications. Time of incubation in extraction buffer was increased to 16–20 hr. Further, only the second elution was retained for library production as this fraction contained the high‐molecular‐weight (HMW) DNA. A total of 100 ng of HMW genomic DNA was digested with FspEI (NEB) at 37°C for 4 hr. Five microlitres of digested DNA was run on an acrylamide gel to verify digestion. Ligation was performed in a 30‐µl reaction with 20 µl of digested DNA, 0.2 µm each of two adaptors, 1 mm ATP and 800 U of T4 DNA ligase (NEB). Ligation was carried out at 4°C overnight. Adaptor and primer sequences are provided in Table S4. Ligation products were amplified in 40‐µl reactions containing 15 µl of ligated product, 0.2 mm of each primer (p1 and p2), 0.3 mm dNTP, 1 × Phusion HF buffer and 0.4 U Phusion high‐fidelity DNA polymerase (NEB). PCR was performed in an Eppendorf Mastercycler with 22 cycles of 98°C for 5 s, 60°C for 20 s, 72°C for 10 s and a final extension of 5 min at 72°C. The target DNA band (~100 bp) was excised from a 10% TBE polyacrylamide gel, and the DNA was diffused from the gel in 50 µl of nuclease‐free water for 12 hr at 4°C. Barcodes were added in a second PCR amplification (total volume 20 µl), which contained 7 µl of gel‐extracted PCR product, 0.2 mm of each primer (p3 and index primer), 0.3 mm dNTP, 1 × Phusion HF buffer and 0.4 U Phusion high‐fidelity DNA polymerase (NEB). Seven PCR cycles with the same profile outlined above were performed. PCR products were purified using a QIAquick minElute PCR purification kit (Qiagen) and then with Ampure XP beads (Agencourt). Libraries were run on an Illumina HiSeq 2500 sequencer using 50 nt single read sequencing.

Raw reads were converted from Illumina 1.8+ to FastQSanger quality scores using fastq groomer (Blankenberg et al., 2010). Illumina adaptors were removed with trim galore! (F. Krueger, Babraham Institute) and trim (Galaxy Version 0.0.1) was used to clip 2 bp off both ends. A reference‐based approach was used by extracting FspE1 methyl sites from an A. palmata genome assembly (version 1.0; Kitchen et al., 2019) using a custom Perl script (provided by Shi Wang) and reads were mapped against these reference sites using bowtie2 (default settings; Langmead & Salzberg, 2012). samtools idxstats (Li et al., 2009) was used to tabulate mapping statistics into counts of mapped reads per reference site for each sample.

Differential DNA methylation analysis was performed in edger (Robinson, McCarthy, & Smyth, 2010). Within the edger vignette, data were converted to reads per million and filtered for only those methylation sites that have at least 10 reads per million in at least six samples. Normalization factors were then calculated to correct for the different compositions of the samples. The effective library sizes are then the product of the actual library sizes and these factors. After library size normalization the data were subjected to nonparametric batch normalization with the combat package for R software (http://www.bu.edu/jlab/wp-assets/ComBat/Abstract.html). Common, trended and tagwise dispersions were estimated and negative binomial generalized linear models were fitted using both the likelihood ratio (LR) test and the more stringent quasilikelihood (QL) F‐test within the edger vignette. E‐values were false discovery rate (FDR)‐corrected.

Methylated sites were matched against the A. palmata genome annotation. Methylated sites were analysed with a principal components analysis within the vegan package in R and visualized using ggbiplot (Wickham, 2016). The R vignette variancepartition (Hoffman & Schadt, 2016) was then used to fit a linear mixed model for each methylation site and partition the total variance into the fraction attributable to each phenotype in the design, including genet, location of sampling (branch vs. not branch), and the condition of the sample (Bleached vs. Not bleached) plus the residual variation (Table S3B).

Differentially methylated sites were tested for Gene Ontology (GO) term functional enrichment with the R/Bioconductor package topgo (Alexa & Rahnenfuhrer, 2018) using the default “weight01” Alexa algorithm with the recommended cutoff of p < 0.05.

RESULTS

Photographic survey results

Prior to the 2014 bleaching event, all sampled colonies had been monitored since at least 2010 and had not bleached until 2014. Following the 2014 bleaching seven of the 12 sampled ramets lost all or nearly all live tissue. Five ramets retained enough live tissue to allow us to evaluate their bleaching condition during the 2015 bleaching event. All five were observed to regain visually normal tissue pigmentation prior to late September 2015 when all five again showed visual signs of bleaching, but generally less severe than in 2014. In many cases the pattern of bleached versus nonbleached areas on a colony was very similar to that observed in 2014 (Figure 2). However, 54% (n = 37) of all colonies present in survey plots (including samples that died) had greater bleaching classifications in 2015, and only 10% (n = 7) of colonies scored one or more classifications of less severely bleached in 2015 than in 2014, indicating that many colonies were more severely impacted by the bleaching of 2015 (Table 2). Survey results thus provided little evidence of short‐term acclimatization.

Table 2

The severity of bleaching was ranked for each surveyed Acropora palmata colony and compared between September 2014 and 2015

(a) Number of observations	2015
	Normal colour	Slightly pale	Pale	Very pale	Bleached	Dead
2014
Normal colour	10	11				3
Slightly pale	2	4	1	2		3
Pale	1	2	5	6		2
Very pale	1		1	5	3	3
Bleached						3

Symbiodinium strain diversity was not correlated with bleaching response

Coral samples were genotyped at 13 microsatellite loci specific for Symbiodinium “fitti” (Pinzón et al., 2011) to determine whether S. “fitti” strain diversity was correlated with bleaching response. Each of the six samples taken per genet was dominated by one S. “fitti” strain with three exceptions (Table 1). Ramet 2 Branch A of Genet P1034 (Grecian Rocks) was dominated by S. “fitti” strain F447 rather than strain F109 like the other five samples from this genet. Ramet 2 Branches B and G of Genet P2656 (Grecian Rocks) contained S. “fitti” strain F477 rather than strain F476 like the other four samples from this genet. These three within‐genet S. “fitti” strain differences were not correlated with the peak bleaching event condition. These findings allow us to discount genetic diversity of S. “fitti”, at least at the strain level, as an explanation for differential bleaching response. However, variable densities of symbiotic dinoflagellates within the host tissues may have influenced bleaching susceptibility (Kemp, Hernandez‐Pech, Iglesias‐Prieto, Fitt, & Schmidt, 2014; Stimson, Sakai, & Sembali, 2002). We do not have prebleaching measurements of symbiont density in tissues, but this may have contributed to variable bleaching within a colony.

Only two species of bacteria differed among coral microbiomes after peak bleaching

Prokaryotic members of the coral microbiome play important roles in coral nutrition, element cycling and disease responses (Peixoto et al., 2017). We thus speculated that the differential bleaching response within and between colonies might be attributable to differences in the microbiome composition. Investigated variables included polyp sampling location (branch vs. not branch), sample peak bleaching condition (bleached vs. not bleached), and host genetic background (genet). The largest partition, not including residuals, was to host genetic background followed by within‐colony sampling location (branch vs. not branch) and then bleaching condition, at the taxonomic rank of both bacterial order and genus (Table 3, Figure 3). The top three bacterial taxa that had the highest variation attributable to host genotype were Methylobacterium (σ 2 = 60%), Rhodobacteraceae (σ 2 = 52%) and Alteromonas (σ 2 = 48.7%). The top three bacterial taxa that had the highest variation attributable to location of the samples (branch vs. not branch) were Cyclobacteriaceae (σ 2 = 41.6%), Tepidiphilus (σ 2 = 25.9%) and Amoebophilus (σ 2 = 23.1%). There were only two bacterial genera in which variance attributable to the bleaching condition of the coral exceeded 5%: Tepidiphilus (σ 2 = 12.1%) and Endozoicomonas (σ 2 = 6.8%). Examining the sub‐OTU sequence output from qiime2 focusing on Tepidiphilus and Endozoicomonas revealed only one strain for both.

Table 3

PERMANOVA results based on Bray–Curtis dissimilarities using abundance data for microbiome community structure at the level of order (a) or genus (b)

Factors	df	Sum Sq	Pseudo‐F	R 2	p
(a) Order
Genet	5	0.69756	2.8573	0.3226	0.001
Residuals	30	1.46477		0.6774
Total	35	2.16232		1.0000
Reef	3	0.36385	2.158	0.16827	1
Residuals	32	1.79848		0.83173
Total	35	2.16232		1.00000
Location	1	0.14355	2.4177	0.06639	0.005
Residuals	34	2.01877		0.93361
Total	35	2.16232		1.00000
Condition	1	0.07035	1.1434	0.03254	0.125
Residuals	34	2.09197		0.96746
Total	35	2.16232		1.00000
(b) Genus
Genet	5	0.70612	3.0791	0.33914	0.001
Residuals	30	1.37596		0.66086
Total	35	2.08208		1.00000
Reef	3	0.3812	2.3906	0.18309	1
Residuals	32	1.7009		0.81691
Total	35	2.0821		1.00000
Location	1	0.11408	1.9709	0.05479	0.024
Residuals	34	1.96800		0.94521
Total	35	2.08208		1.00000
Condition	1	0.05864	0.9854	0.02816	0.169
Residuals	34	2.02344		0.97184
Total	35	2.08208		1.00000

For the Reef, Sampling Location and Bleaching condition comparisons data were blocked using genet designation. df, degrees of freedom; Sum Sq, sum of squares; Pseudo‐F, F value by permutation; bold type indicates statistical significance at p < 0.05, p‐values based on 999 permutations.

Figure 3

Acropora palmata genets are characterized by between and within genet variation in their microbiomes. (a) Principal coordinate analysis of the normalized 16S tag sequence data generated using the package vegan and ggbiplot in R at the taxonomic rank of order. (b) Violin plot showing the contribution of Genet, Location sampled (Branch, or not Branch), the Condition of sample (Bleached or not bleached) and the residuals to variation in the 16S tag normalized microbiome data at the taxonomic rank of order. Genet contributed most of the explained variation, followed by location and condition. Plot shows the kernel probability density of the data at different values, box indicates the interquartile range and horizontal bar indicates median. Generated using the variancepartition package in R [Colour figure can be viewed at http://wileyonlinelibrary.com]

The methylome differs among coral host genets

Image analyses of bleached and unbleached portions of the colonies suggested that micro‐environmental differences such as shading or exposure to water movement might have contributed to differential bleaching responses (Figure 2). Thus, we investigated an epigenetic mechanism of acclimatization, DNA methylation, to address whether different portions of a colony or genet had acclimatized to micro‐environmental differences.

A total of 28,797 sites were analysed for differential methylation, and 64.4% (n = 18,538) of these sites are within predicted genes in the Acropora palmata genome (version 1). Variation attributable to different phenotypic characteristics was determined, and the largest partition, not including residuals, was to host genetic background followed by within‐colony polyp sampling location, i.e., whether a branch location or not (Figure 4b), and lastly whether the polyps were previously bleached or not bleached.

Figure 4

Acropora palmata genets are characterized by between and within genet variation in genome methylation. (a) Principal coordinate analysis of the normalized coral methylome data generated using the package vegan and ggbiplot in R. (b) Violin plot showing the contribution of Genet, Location sampled (Branch, or not Branch), the Condition of sample (bleached or not bleached) and the residuals to variation in the normalized coral methylome data. Genet contributed most of the explained variation, followed by location and condition. Generated using the variancepartition package in R [Colour figure can be viewed at http://wileyonlinelibrary.com]

To determine the biological function of genes in which methylation was more abundant, a GO‐enrichment analysis was completed comparing all methyl sites (n = 6,848 had GO term annotations out of 28,797 sites) to all annotated genes in the A. palmata genome (n = 25,102 had GO term annotation). Significant biological functions included RNA binding, zinc ion binding, transition metal binding, sequence‐specific DNA binding and oxidoreductase activity (Table S5).

Acropora palmata genets showed large and consistent differences in methylation patterns, and clustering analyses clearly grouped samples belonging to the same coral genet together, but also showed some variation of methylation within genets (Figure 4a). On average 9.5% (n = 2,741, QL F‐test) of sites were significant for differential methylation between genets and, of these sites, 57.8% were within predicted genes (n = 1,584). In a GO enrichment of the differentially methylated sites between genets, the top five biological processes were the negative regulation of molecular function, the negative regulation of hydrolase activity, the modification by symbiont of host morphology or physiology, the negative regulation of proteolysis, and the regulation of DNA‐binding transcription factor activity (top 20 nodes, Biological processes: Table S6, Molecular function: Table S7).

Reef location, peak bleaching and intracolony polyp location

Next, we determined whether there were methylation changes attributable to the reef environment. On average 5% (n = 1,511, QL F‐test) of sites were differentially methylated among reefs. Unfortunately, we could not block for genetic background in the QL F‐test analysis because not all reefs were represented by multiple genets. However, a Venn diagram analysis comparing which CCGG sites are differentially methylated between reefs versus those differentially methylated between genets showed there was a 40.1% (n = 1,218) overlap between the factors. Only 9.7% (n = 293) of the differentially methylated sites were exclusively found in the reef comparison compared to 50.2% (n = 1,523) found exclusively in the genet comparison (Figure S1). A GO enrichment for biological processes for these differentially methylated sites exclusively found in comparisons between reefs but not between genets (n = 293) showed that the top five processes were in autophagy and the regulation thereof, amine transport including both regulation and positive regulation, and the regulation of cellular catabolic processes (top 20 nodes, Biological processes: Table S8, Molecular function: Table S9).

Methylation variation within genets was also attributable to whether tissues had recently bleached. In this comparison, we included 18 bleached and 18 nonbleached samples and blocked for genet. One methylation site (Segkk362_pilon‐108100) was significantly different (p < 0.01, FDR‐corrected) when using a LR test, but not when using the more stringent QL F‐test. This methylation site lies within a genomic region that does not have a gene prediction or annotation in A. palmata. The closest gene annotation, at a little over 10,000 bp away, is Transposon TX1 uncharacterized 149‐kDa protein as found in the scleractinian coral Stylophora pistillata.

Next, we looked at whether samples collected from branches versus those collected from the base or trunk of the colony (i.e., not branch) had conserved methylation differences shared among all genets. In this comparison (n = 26 branch, n = 10 not branch) when blocking for genet, there was one methylation site (Segkk1116_pilon‐156192) that was differently methylated (p < 0.001, FDR‐corrected) when using the more stringent QL F‐test and 84 methylation sites when using the LR test (Table S10). Segkk1116_pilon‐156192 is within the gene Galactosylceramide sulfotransferase. Of the 84 differentially methylated sites found in the LR test, 61.9% (n = 52) are within predicted genes. A GO‐enrichment analysis of these differentially methylated sites, using topgo in R, revealed a significant enrichment for the biological processes (Table S11) of AMP biosynthetic and metabolic processes, regulation of fibroblast apoptotic processes, excretion, protein–chromophore linkage, the regulation of viral genome replication, viral life cycle, and viral processes, the regulation of alternative mRNA splicing via spliceosome, immune effector processes, the negative regulation of apoptotic processes and programmed cell death (for enrichment in molecular function see Table S12).

DISCUSSION

Evolutionary rescue of populations depends on a species’ ability to produce phenotypic variation that is heritable and adaptive. DNA mutations are the best understood mechanisms to create phenotypic variation, but other, less well‐studied mechanisms exist. Environmental conditions frequently change over the long lifespan of some reef‐building coral genets (Devlin‐Durante, Miller, Caribbean Acropora Research Group, Precht, & Baums, 2016). Phenotypic plasticity in response to environmental variation is a common trait in corals; however, the mechanisms by which corals achieve this plasticity are not well understood. This can be partly attributed to the difficulty of separating host response from those of algal symbionts and other members of the microbiome. Acropora palmata colonies showed variable responses to increased water temperatures, yet the identity of the associated symbiont strain did not explain this plasticity. In addition, photographic survey results showed little evidence of short‐term acclimatization. Interestingly, the abundance of two species of prokaryotes, Tepidiphilus and Endozoicomonas, differed among microbiome samples depending on their bleaching response and this supports findings that Endozoicomonas plays an important role in the functioning of the coral holobiont (Pollock et al., 2018). Ultradeep sequencing of host and symbionts may reveal the occurrence of somatic mutations in hosts (and/or symbionts) correlated with the bleaching phenotype (Van Oppen, Souter, Howells, Heyward, & Berkelmans, 2011) but nonmutation‐based mechanisms may also play a role (Goldsmith & Tawfik, 2009; Payne & Wagner, 2019), including detection‐based epigenetic modifications (sensu Shea, Pen, & Uller, 2011).

Here we, for the first time, applied a new technique to assay methylation status in wild corals and have found that coral genets show large and consistent differences in the way they methylate their DNA even when growing on the same reef, consistent with the hypothesis that some portion of their methylation patterns was inherited (Dimond, Gamblewood, & Roberts, 2017; Liew, Howells, et al., 2018). Furthermore, there was significant variation in methylation state that was correlated with the location of the polyps within colonies and the bleaching response of the colonies. Changes in epigenetic marks based on the detection of environmental cues provides an avenue to effect phenotypic plasticity (Shea et al., 2011). Detailed data on systematic variation in temperature, flow and light regimes depending on polyp location within a colony and colony morphology are accumulating (Edmunds & Burgess, 2018; Ong, King, Kaandorp, Mullins, & Caley, 2017; Stocking, Rippe, & Reidenbach, 2016). Remarkably, there was some correspondence between genes previously shown to be differentially expressed, and those that were differentially methylated, indicating that methylation differences may translate into gene expression differences. These results are novel because they reveal a potential pathway by which these long‐lived corals can modify their phenotypes in response to the environment. Because these corals can produce large monoclonal stands, modifications of genome methylation produce mosaics of phenotypes despite low genotypic diversity. Whole genome bisulphite sequencing analysis should be performed to further investigate these patterns. Future work on this and other marine foundation species may significantly advance our understanding of the mechanism determining evolvability of threatened species.

Microbial assemblages and symbiotic dinoflagellates are minimally correlated with bleaching condition

Corals live in an intimate symbiosis with algae in the family Symbiodiniaceae (LaJeunesse et al., 2018). While some coral species can host several species of algae, the elkhorn coral, A. palmata, usually hosts Symbiodinium “fitti” and furthermore ramets of the same coral genet often retain the same symbiont strain (Baums et al., 2014). This was also the case here. Genets differed in the S. “fitti” strain they hosted but within each genet, there was remarkable homogeneity. These findings allow us to discount genetic diversity of S. “fitti”, at least at the strain level, as an explanation for differential bleaching responses. However, variable densities of symbiotic dinoflagellates within the host tissues may have influenced bleaching susceptibility (Kemp et al., 2014; Stimson et al., 2002). We do not have prebleaching measurements of symbiont density in tissues, but this may have contributed to variable bleaching within a colony.

Genetic diversity of S. “fitti” below the strain level may also influence the bleaching response. A. palmata colonies host large populations of S. “fitti” that are the result of countless cell divisions. Somatic mutations may accumulate and lead to functional variation among cells of the same S. “fitti” strain, as defined by a shared multilocus microsatellite genotype profile. Experimental evolution of Symbiodiniaceae cultures has recently been shown to result in functional divergence (Chakravarti & van Oppen, 2018). Ultradeep sequencing of S. “fitti” isolated from replicate samples representing large A. palmata genets would probably detect such somatic mutations (Wang et al., 2019), but assigning function to these mutations in symbionts associated with wild coral colonies remains a daunting challenge.

Prokaryotic members of the coral microbiome play important roles in coral nutrition, element cycling and disease responses (Peixoto et al., 2017). We thus speculated that the differential bleaching response within and between colonies might be attributable to differences in the microbiome species composition. However, there were only two bacterial genera in which variance attributable to the bleaching condition of the coral exceeded 5%: Tepidiphilus (σ 2 = 12.1%) and the coral symbiont Endozoicomonas (σ 2 = 6.8%). Each taxon was represented by one strain only. There are three proposed functions of Endozoicomonas in the host coral, including maintaining the structure of the host microbiome, nutrient acquisition and provision, and a role in host health and/or disease (Neave et al., 2017). Evidence for a role in coral health was demonstrated through the altered abundances of Endozoicomonas in multiple scleractinian corals in relation to seawater pH (Webster et al., 2016), sedimentation and wastewater runoff (Ziegler et al., 2016), and the occurrence of diseased lesions (Meyer, Paul, & Teplitski, 2014). Differences in the microbiome that pre‐dated or immediately followed the bleaching event were obscured, because we were only able to sample tissues once they had recovered from bleaching. Previous studies looking at stress events in corals found a shift of the coral holobiont to a more potentially pathogenic state with disease‐associated bacteria and fungi (Thurber et al., 2009), including the fungus Ascomycota (Thurber et al., 2009) and the bacterial taxa Vibrionales (Thurber et al., 2009; Zaneveld et al., 2016) and Oscillatoriales (Zaneveld et al., 2016). A shift to a pathogenic state was not observed in our colonies, indicating that whatever impacts bleaching may have had on the microbiome, those effects were difficult to detect 6 weeks post‐stress using a standard 16S microbiome analysis approach. Future metagenome sequencing may reveal genetic diversity among prokaryotes that shared the same 16S sequencing tag and help to explain some of the unaccounted‐for variance among ramets of the same genet.

Variation in methylation patterns by genet

It is important to distinguish inheritance of epigenetic marks in the soma versus the germline (Shea et al., 2011). During cell proliferation, somatic tissues inherit epigenetic marks from progenitor cells to, for example, give an epithelial cell its identity. Similarly, polyps within a colony and colonies belonging to the same genet share epigenetic signatures via somatic inheritance. However, when comparing genome methylation patterns between A. palmata genets, we found large differences even when the genets lived on the same reef (Figure 2). This suggests that at least some portion of genome methylation was inherited, otherwise a shared environment post‐fertilization should lead to shared methylation patterns among genets. Evidence for inherited germline methylation patterns in corals is accumulating (Dimond et al., 2017). Widespread depletion of CpG dinucleotides was observed in Acropora millepora and is a signature for historical germline DNA methylation (Dixon, Bay, & Matz, 2014). A recent study in Platygyra daedalea for the first time demonstrated intergenerational inheritance of DNA methylation patterns in corals, from parent to sperm, and evidence for maternal and paternal effects in larvae from reciprocal crosses (Liew, Howells, et al., 2018). In most animals with early germline segregation, epigenetic marks are reset during meiosis, and the mechanisms of inheritance of epigenetic marks in coral embryos are unknown (reviewed by Eirin‐Lopez & Putnam, 2019).

Symbiont–host interactions may also influence host genome methylation patterns. Even genets that grew near each other on the same reef hosted a different strain of S. “fitti” while ramets of the same genet usually shared an S. “fitti” strain. Although previously undocumented, S. “fitti” strains may differentially alter the host methylome. Supporting evidence for this comes from genet by genet comparison in a GO enrichment analysis, which showed that the category “modification by symbiont of host morphology or physiology and modulation by symbiont of host cellular processes” was enriched. Differentially methylated genes in this category included the Homeodomain‐interaction protein kinase 2, eIF‐2‐alpha kinase GCN2 (three separate methylation sites within this gene), Gag‐Pol polyprotein (two separate methylation sites within this gene) and TNFAIP3‐interating protein 1. Genotype/genotype interactions between Symbiodiniaceae strains and host genets and their effects on host methylomes deserve further study (reviewed by Parkinson & Baums, 2014).

Methylation patterns vary between locations within the colony

Some of the variance in genome methylation within genets was attributable to long‐term microenvironmental conditions between polyp locations within a colony rather than eukaryotic or bacterial symbiont community composition. Complex skeletal morphologies and varying tissue layer thickness create a variety of intracolonial light microniches, but in general the top of a branch will experience significantly higher solar irradiance than the base or trunk of a colony (Kaniewska et al., 2011; Wangpraseurt, Larkum, Ralph, & Kühl, 2012; Warner & Berry‐Lowe, 2006), and therefore polyps on the tops of branches have an increased need to avoid the damaging effects of excess light energy and the resulting oxidative stress. In A. globiceps, Symbiodiniaceae densities were consistent between internal and external branches but varied with depths (greater densities at lower depths) (Ladrière et al., 2013). The host coral can regulate Symbiodiniaceae densities through nutrient limitation or through digesting or expelling the excess symbiotic algae to maintain relatively low and consistent densities (Dunn, Bythell, Tissier, Burnett, & Thomason, 2002; Falkowski, Dubinsky, Muscatine, & McCloskey, 1993; Muscatine et al., 1998) as one option to reduce oxidative stress in areas or times of higher light exposure (Fitt, McFarland, Warner, & Chilcoat, 2000). The various mechanisms employed by the host, including epigenetic changes as suggested here, help the colony avoid bleaching even in the high‐light exposed polyp locations under nonstressful temperature conditions.

Abnormally high temperatures accompanied by high irradiance can cause a breakdown in the coral–symbiotic algae symbiosis, resulting in expulsion of the algae, a process referred to as bleaching. We observed a higher incidence of bleaching in samples from branch regions (58%) versus those collected from the base (30%) in this study. In addition, there is also strong evidence for a division of labour between coral branch tips and bases in their gene expression (Hemond, Kaluziak, & Vollmer, 2014). Interestingly, galactosylceramide sulfotransferase was significantly differentially methylated between polyps sampled from branches versus other locations within colonies such as the base or trunk (Table S10). This result was obtained with both likelihood ratio T‐test and QL F‐test statistics. Galactosylceramide sulfotransferase is involved in sphingolipid metabolism and was also differentially expressed by colony position, being upregulated in branch tips in A. palmata and A. cervicornis (Hemond et al., 2014). Sphingolipid metabolism may be involved in the regulation of algal symbionts. In anemones, the sphingosine rheostat can regulate the balance between stability and dysfunction in the cnidarian–dinoflagellate partnership (Detournay & Weis, 2011). However, under stressful temperatures, parts of the colonies that were exposed to high irradiation were unable to avoid bleaching even with the epigenetic modifications driven by the internal light gradient.

Symbiodiniaceae in shallow corals must dissipate four times more light energy than is needed for photosynthesis on a bright summer day (Gorbunov, Kolber, Lesser, & Falkowski, 2001). This excess light energy absorbed by chlorophyll can be dissipated through heat loss, re‐emitted as fluorescence, or decayed via the chlorophyll triplet state that produces reactive oxygen species as a byproduct. Here, we identified differentially methylated sites by polyp location that were overrepresented in the GO categories of AMP biosynthetic and metabolic processes. The site is in the gene adenylosuccinate lyase, which catalyses two key steps in AMP synthesis. cAMP induces gene transcription through activation of cAMP‐dependent protein kinase (PKA) and subsequently activation of transcription factors including CREB (cAMP response element binding proteins)/ATF transcription factor family members such as CREM and ATF1 via phosphorylation by PKA. In a differential gene expression analysis of A. palmata fragments kept in complete darkness for 9 days compared to controls, two annotated genes were identified, cAMP‐responsive element modulator and cyclic AMP‐dependent transcription factor ATF‐4 (DeSalvo, Estrada, Sunagawa, & Medina, 2012). The ATF‐4 transcription factor responds to oxidative stress and amino acid starvation (Harding et al., 2003).

The coral skeleton serves as an efficient light‐capturing device and colony and polyp morphology determine the light levels experienced by the intracellular symbionts (Enríquez, Méndez, & Iglesias‐Prieto, 2005; Swain et al., 2018). Symbiodiniaceae can also maximize light absorption and utilization by increasing photosynthetic pigments and photosynthetic efficiency in corals acclimatized to low light (Falkowski & Dubinsky, 1981). Among the enriched GO terms between branch and nonbranch polyp locations was the category protein‐chromophore linkage. The differentially methylated gene was in cyrptochrome‐1. Cryptochromes are flavoproteins that are sensitive to blue light. They regulate the circadian clock in plants and animals. Eight core circadian genes have been identified: Casein kinase 1e (CK1e), Cryptochrome1 (Cry1), Cryptochrome2 (Cry2), Period1 (Per1), Period2 (Per2), Period3 (Per3), Clock and BMAL1 (brain and muscle ARNT‐like protein, Arntl, MOP3). Cryptochrome1 and 2 have been previously reported to display diurnal patterns of transcription in corals, with higher expression found in the light phase than in the dark (Hoadley, Szmant, & Pyott, 2011; Levy et al., 2007, 2011). Cry1 in A. millepora is not under control of an endogenous clock whereas Cry2 is. Both have higher expression during the day (Brady, Snyder, & Vize, 2011). Circadian clock genes affect a large number of downstream processes and thus serve as important nodes in transcriptional networks (Dunlap, 1999) and in the regulation of post‐translational modifications (Gallego & Virshup, 2007; Staiger & Koster, 2011). Differential methylation of these genes may thus be an effective means to alter the transcription of several downstream pathways in response to differential light levels within colonies. Future research related to the transcription and differential methylation of these circadian clock genes in corals is warranted.

We unexpectedly found that the GO terms for the regulation of viral genome replication, viral life cycle and viral processes were enriched in the comparison between branch and nonbranch locations. We are not aware of any data indicating that viral load differs within colonies. Because this GO term was enriched across genets and reefs, we would expect viral loads to differ systematically between branch and nonbranch locations and this hypothesis deserves future testing. Interestingly, the bacterial communities did differ between the tips and the bases of colonies, suggesting that the branches may harbour a specialized microbiome. There is contradicting previous evidence with respect to within‐colony variation of the prokaryotic community in corals. In a previous study on A. palmata, no detectable community‐level differences were found among the prokaryotic microbiota of the uppermost, underside and base of A. palmata (R 2 = 0.20, p = 0.51) (Kemp et al., 2015). In contrast, considerable within‐colony variation of bacterial assemblages was found in Orbicella annularis between the tops and the sides (Daniels et al., 2011). O. annularis also harbours several species of Symbiodiniaceae, and hence further research is required to understand what factors drive within‐coral‐colony diversity of the prokaryotic community.

Methylation patterns vary with bleaching history

Methylation variation within genets was also attributable, to some extent, to whether tissues had recently bleached. The one significant methylation site was located within a genomic region that does not have a gene prediction or annotation. The closest gene annotation, at a little over 10,000 bp away, is Transposon TX1 uncharacterized 149‐kDa protein as found in Stylophora pistillata. By the alteration of splicing and polyadenylation patterns or through functioning as enhancers or promoters, transposable elements can exercise control over neighbouring genes (Girard & Freeling, 1999). Transposable elements are significantly differentially expressed in response to heat stress in corals (DeSalvo et al., 2010; Traylor‐Knowles, Rose, Sheets, & Palumbi, 2017) and in plants (Ito et al., 2011; Pecinka et al., 2010). Overall, less of the variation in methylation was explained by previous bleaching, suggesting that methylation changes may be effective in changing coral transcription in response to longer term differences in the light environment (e.g., between tips and bases) rather than more acute temperature stressors. Future research is required to correlate gene expression and methylation over a range of stress exposures and stress severity.

High fragmentation rates and acute stress events necessitate that A. palmata polyps acclimatize to changes in environmental conditions. Our data suggest that acclimatization may be partially achieved via differential methylation. We suggest that differential genome methylation may be one of the mechanisms by which corals achieve their remarkable phenotypic plasticity in their natural environment. In a transplant experiment in A. millepora, GBM changed subtly, but much less than transcription (Dixon et al., 2018). Dixon et al. also found strong associations between gene body methylation and fitness, although gene body methylation was not directly correlated to transcription, resulting in the authors’ questioning what mechanism connects gene body methylation to phenotype and fitness (Dixon et al., 2018). Yet, methylation is known to affect transcription factor binding both negatively and positively, and thus alter transcriptional regulation bidirectionally, making correlation to gene expression complicated (Yin et al., 2017).

A significant amount of the intragenet variation in phenotypic stress response observed here remains to be explained. Similarly, large residual variances in methylation that could not be attributed to any of the studied factors were observed in Porites porites (Dimond et al., 2017). Epimutations are stochastic events that result in random additions and losses of epigenetic marks (reviewed by Johannes & Schmitz, 2019), some of which are heritable. Distinguishing stochastic versus directed changes in epigenetic marks in corals will require careful experimentation. Plant researchers have made use of mutation accumulation lines grown under controlled conditions to make the distinction (e.g., Becker et al., 2011). Further, ultradeep sequencing of host and symbionts may reveal somatic mutations correlated with the bleaching phenotype (Van Oppen et al., 2011), although mechanisms other than detection‐based DNA methylation changes may also play a role (Goldsmith & Tawfik, 2009; Payne & Wagner, 2019), including selection‐based changes in epigenetic marks (Shea et al., 2011). Future work on this and other marine foundation species may significantly advance our understanding of these mechanisms in determining the evolvability of threatened species.

AUTHOR CONTRIBUTIONS

I.B., M.D. and D.W. designed the study and obtained funding. D.K., M.D., S.V and D.W. performed research and analysed the data. M.D. and I.B. wrote the paper. I.B. supervised the study. All authors edited the paper.

Click here for additional data file.