Recently hatched juveniles of A. polyacanthus were collected from February to March 2009 from the reef slope at Heron Island (23°27′ S, 151°57′ E) on the Great Barrier Reef (Figure 1A). At this location the average water temperature at 6–8 m depth below sea level was 27°C in summer and 21.8°C in winter (from 1999 to 2008 Australian Institute of Marine Science¹). Fish were transported to James Cook University’s aquarium facility where 120 juvenile fish (<3 months old) were randomly divided into three seasonally-adjusted temperature treatments modeled off the previous 10 years for the collection location: a present-day temperature (control: +0°C) and two within-generational elevated temperature treatments (+1.5 and +3°C above present-day). Fish were maintained in 60 L replicate tanks for each temperature treatment until maturity (11–14 months of age), as previously described (Donelson and Munday, 2012) (Figure 1B). Sampling of fish for tissue was completed during the Austral summer of 2009–2010 when water temperature for the three treatments were +0: 27°C, +1.5: 28.5°C, and +3: 30°C.

Liver was the target tissue because it is a major metabolic organ. Liver is an important long-term storage site for excess energy as glycogen, allowing glucose release to the blood and provision of fuel to the body as required (Wasserman, 2009). Liver also synthesizes bile and heme for fat absorption and oxygen delivery, respectively (Campbell, 2006). Adult fish from each treatment were sampled and whole livers immediately dissected, snap-frozen in liquid nitrogen, and then stored at −80°C until further processing. To avoid sex-related epigenome variations we used only liver tissue from male fish.

Liver genomic DNA (gDNA) of 8 fish (3, 2, and 3 fish from control, within-generational +1.5°C, and within-generational +3°C groups, respectively) was extracted by homogenizing 100 mg of frozen liver tissue in 700 μL CTAB (hexadecyltrimethylammonium bromide) buffer. After adding proteinase K, the homogenate was incubated overnight at 55°C and then RNase treated for 10 min at room temperature. A standard chloroform/isoamyl alcohol extraction and ethanol precipitation method (Sambrook and Russell, 2001) was used. DNA quality and quantity were checked using NanoDrop Spectrophotometer absorbance readings (Invitrogen, Mulgrave, VIC, Australia) and on a 0.8% agarose gel.

gDNA was bisulfite-converted using Methyl-MaxiSeq^TM Kit and sequenced using llumina HiSeq1500/2500 system by Zymo Research (Irvine, CA, United States) as described previously (Ryu et al., 2018). Briefly, 500 ng of genomic DNA was treated with dsDNA Shearase^TM Plus (Zymo Research), end-blunted, 3′ A-extended, and purified by DNA Clean & Concentrator^TM –5 kit (Zymo Research). DNA fragments were then ligated with adapters and filled in. EZ DNA Methylation–Lightning kit (Zymo Research) was used for bisulfite treatment and DNA fragments were sequenced on Illumina HiSeq1500/2500. All procedures were conducted according to the manufacturer’s instructions. 1.26 ∼ 1.37 × 10⁹ reads of 101 bp were generated for each sample (Supplementary Table 1). The methylomic reads for this study were generated together with those from our previous work (Ryu et al., 2018). The quality of these methylome datasets were validated using targeted bisulfite sequencing assays in the same study.

Draft genome assembly and gene models for A. polyacanthus from Heron Island were generated as described in our previous work (Veilleux et al., 2018). The assembled genome consisted of 16,609 genomic scaffolds ranging from 500 ∼ 3,215,819 bp with N50 of 527,731 bp. Base coverage was 129X on average and assembly completeness measured by CEGMA (Parra et al., 2007) was 99.19%. 23,464 gene models were identified based on ab initio and A. polyacanthus transcriptome-based gene prediction by Maker2 pipeline (Holt and Yandell, 2011).

Adapters and low quality bases in reads from the whole genome bisulfite sequencing (WGBS) were trimmed by Trim Galore v0.4 (Krueger, 2015). Trimmed reads were mapped to A. polyacanthus genome sequences from Heron Island using Bismark v0.13.1 (Krueger and Andrews, 2011) with the ‘–non-_directional’ option. Mapping efficiency was 67.17 ∼ 71.23% per sample. The bisulfite conversion rate is often used as a metric to assess the efficiency of bisulfite treatment and is calculated as numT/(numT +numC) × 100, where numC and numT represent the number of times that cytosines and thymines are observed, respectively at CHG and CHH sites. We used methylKit v.1.2.0 (Akalin et al., 2012) to identify differentially methylated regions (DMRs) as follows. Per base coverage of methylomic reads on the genome was 43 ∼53 × calculated by ‘getCoverageStats’ function (Supplementary Table 2). The ‘methRead’ function scanned the methylome mapping information and kept cytosines that were covered by a minimum of 10 reads. Highly covered bases (>99.9%) were excluded to remove PCR bias by ‘filterByCoverage’ function. Variances across entire samples were normalized by ‘normalizeCoverage’ functions. The genome was broadly categorized as eight functional units (CpG island, CpG shore, promoter, 5′ untranslated region (UTR), exon, introns, 3′ UTR, and repeats), following our previous study (Ryu et al., 2018). Methylation status of individual cytosines in each genomic region were summed per sample by the ‘regionCounts’ function. Statistical significance of differential methylation between samples was calculated using the Chi-squared test of the ‘calculateDiffMeth’ function. Multiple testing correction was performed by the Benjamini & Hochberg (BH) method of the ‘p.adjust’ function from p-values from all pairwise comparisons. DMRs for three methylation contexts (CpG, CHH, CHG; H = A, C, T) were identified using the two thresholds: adjusted p-value ≤ 0.05 and ≥ 33.3% methylation difference between two treatments. The adjacent genes to each DMR on the same scaffold were assigned using the ‘bedtools closest’ v2.23 (Quinlan, 2014).

The ‘heatmap.2’ function of the ‘gplots’ package (Warnes et al., 2016) was used to draw heatmaps. For multidimensional scaling analysis, dissimilarity between samples was calculated as ‘1 – Spearman correlation coefficient’ and confidence ellipsoid of each sample was obtained using the ‘bootmds’ function of the ‘smacof’ package (De Leeuw and Mair, 2011) with the setting of alpha = 0.05 and 1,000 bootstrap replications. Statistical significance of separation between experimental groups was assessed by ind.ctest function of the ‘ICSNP’ R package (Nordhausen et al., 2015). p-Values were adjusted for multiple testing correction using the BH method.

Homologs of A. polyacanthus proteins were identified by a BLASTP search against the NCBI non-redundant protein database with e-value of 10^–4 as a threshold (Altschul et al., 1990). Gene Ontology (GO) terms of homologous proteins were then assigned to A. polyacanthus proteins using the NCBI gene2go.gz file² and custom script. GO terms were considered when associated with more than two genes. Only Biological Process among the three GO categories (Biological Process, Molecular Function, and Cellular Component) were analyzed to focus on the biological outcome of each gene rather than a specific molecular activity (Molecular Function) or a location relative to cellular structures (Cellular Component) of the gene (Gene Ontology Consortium, 2008). Statistical significance of GO enrichment in DMR genes from all groups or between experimental groups were calculated by cumulative hypergeometric test and adjusted by the BH method using the custom R script respectively.

To allow comparison of DNA methylation profiles for fish exposed to elevated temperature within- and between-generations, we used the methylome dataset from our previous study (Ryu et al., 2018). This dataset was generated from F2 A. polyacanthus from the Palm Island population (central Great Barrier Reef, Australia; 18°37′ S, 146°30′ E) that had been reared at control and elevated temperatures for two generations (Figure 1A). Briefly, F1 offspring (<3 months old) of wild-caught adult fish were split between three temperature treatments: a present-day control for the collection location +0°C (mean winter: 23.2°C and summer: 28.5°C) as well as +1.5 and +3°C above present-day (seasonally cycling temperature treatments). Fish were reared in treatment conditions in 40–60 L aquaria until 2 years of age when fish were mature. At this time pairs of fish from unrelated families were randomly mated (Donelson et al., 2012 for more details). F2 offspring from the three F1 parental treatments were divided between two F2 developmental conditions, +0 and +3°C immediately after hatching and remained in treatment until 2 years of age. This produced four F2 treatment groups: control (F1: +0°C, F2: +0°C), within-generational (F1: +0°C, F2: +3°C), step (F1: +1.5°C, F2: +3°C), and transgenerational (F1: +3°C, F2: +3°C) (Figure 1B).

To minimize the effect of sex on the downstream analysis, two males and two females were selected for each Palm Island population treatment. DNA extraction from adult liver and WGBS was conducted as written above for Heron Island fish. The methylomic reads were aligned to the genome assembly from Palm Island fish. Each sample has 0.95 ∼ 1.36 × 10⁹ reads of 101 bp and yielded 34 ∼ 62 × coverage depth after alignment to the genome sequences. The same computational methods used to identify DMRs in Heron Island populations were used except for one parameter: to remove the influence of unbalanced sex ratio in this population, the information of males and females was included as the ‘covariates’ parameter in the ‘calculateDiffMeth’ function of methylKit (Akalin et al., 2012). DMRs were identified using the same thresholds (adjusted p-value ≤ 0.05 and ≥ 33.3% methylation difference between two treatments), which resulted in 445 CpG DMR genes.

Orthology was calculated between Heron and Palm Island populations’ gene models using Inparanoid (O’brien et al., 2005). This resulted in 17,423 ortholog groups with 17,569 and 17,758 genes from each population, respectively.

We extracted methylomes from liver samples that were collected from Heron Island populations of A. polyacanthus exposed to control (+0°C) and increased temperatures (+1.5 and +3.0°C) from early life (Figure 1B). This resulted in 841.6 ∼ 884.8 million reads that were uniquely mapped, yielding from 67.2 to 71.2% mapping efficiency (Supplementary Table 2). Methylated cytosines ranged from 61.1 to 63.4% in the CpG context and from 0.77 to 3.68 in CHG and CHH contexts, across samples (Supplementary Table 2). Average bisulfite non-conversion rates were from 0.57 to 2.91% per sample, which is comparable to our previous study (0.84 to 3.46%) (Ryu et al., 2018). Base coverage of methylome reads were ≥ 43 × for the three methylation contexts for all Heron Island samples (Supplementary Table 2).

Differential methylation in CpG, CHH, and CHG contexts was calculated between pairs of treatments and 480, 88, and 3 DMRs were identified, respectively (Table 1). As there were so few CHG DMRs this context was not analyzed further. CpG DMRs were mostly located on exons (22.24%), repeats (21.72%), and introns (20.81%) whilst CHH DMRs were mostly located in repeats (42.66%), introns (16.78%), and CpG shores (13.99%) (Figure 2C and Supplementary Figure 1C, Tables 4, 5).

To determine if within-generational exposure to increased temperatures in the Heron Island population affected global DNA methylation patterns, we performed a multidimensional scaling analysis using the methylation level of DMRs. Individual fish samples were clustered by CpG DMRs and were separated by temperature treatment (Figures 2A,B). Statistical tests also showed clear separation of control samples from the elevated temperature treatments (Supplementary Table 3). These patterns were less evident for the CHH context (Supplementary Figures 1A,B and Table 3), which is similar to our previous observation of methylation patterns in control, within-generational, step, and transgenerational in Palm Island fish (Ryu et al., 2018).

We identified 372 CpG DMR genes (nearest genes to any given CpG DMRs, Supplementary Table 5) in the Heron Island population. GO enrichment analysis of these 372 genes showed 20 enriched terms (adjusted p-value < 0.05) (Supplementary Table 6) mostly related to animal development (e.g., ‘cardiac neural crest cell development involved in heart development,’ ‘central nervous system projection neuron axonogenesis,’ and ‘retina layer formation’) and physiological homeostasis (e.g., ‘cellular chloride ion homeostasis,’ ‘insulin catabolic process,’ and ‘hypotonic response’).

Further investigation of enriched GO terms identified by comparing treatments in a pairwise manner also showed similar but slightly different categories that were treatment-dependent. Most enriched GO terms were related to pathways involved in developmental processes and cardiovascular functions, and the majority of these were found when comparing control to within-generational +1.5°C treated fish (Figure 3). GO terms related to diverse developmental processes, such as ‘bone morphogenesis,’ ‘multicellular organism growth,’ ‘type B pancreatic cell development,’ and ‘neural retina development’ were the most prevalent functions. The enriched terms for cardiovascular physiology included ‘cardiac muscle hypertrophy,’ ‘cardiac myofibril assembly,’ and ‘regulation of heart rate. Nutrient control and heat stress functions were also noteworthy. ‘Insulin catabolic process’ was the most significantly enriched GO term in hypomethylated DMR genes from both elevated temperature groups (both +1.5 and +3°C) compared to control. ‘Glucose homeostasis’ was another enriched term in the comparison between control and within-generational +1.5°C. The GO term ‘regulation of cellular response to heat’ was enriched among DMR genes that were hypomethylated in within-generational +1.5°C fish compared to control.

This analysis also revealed that functions differed depending on the magnitude of temperature change. The greatest number of enriched GO terms were found when the fish reared in elevated temperatures were compared to control fish, but not when the elevated temperature fish (+1.5 and +3°C groups) were compared to each other (Figure 3). For example, the largest number of enriched terms were found among hypomethylated DMR genes in within-generational +1.5°C fish compared to control fish (26 terms), while only two terms were enriched in within-generational +3°C fish compared to within-generational +1.5°C fish. This suggests that major change of physiological functions occurs when a temperature change is experienced during early life, but the magnitude of change, either +1.5 and +3°C, is relatively unimportant.

To compare epigenetic mechanisms between within- and between-generational plasticity, we performed the same GO enrichment analysis on the DNA methylome dataset of the F2 population from Palm Island (Ryu et al., 2018). First, we identified 445 CpG DMR genes by pairwise comparison of WGBS data from the four Palm Island treatments (control, within-generational +3°C, step, and transgenerational) (Supplementary Table 7). Enrichment analysis on all 445 genes in the F2 Palm Island population identified only three significant GO terms (adjusted p-value < 0.05) (Supplementary Table 8), while we identified 20 different GO terms among the 372 CpG DMR genes in the Heron Island population (Supplementary Table 6). Enrichment analysis per treatment pair also showed different patterns with regard to the Heron Island population (Supplementary Figure 2). Five or fewer GO terms were enriched across treatment pairs and there was no overlap with the enriched terms from the Heron Island population.

Because there was no overlap of GO terms between two populations, we examined if there was similarity between within- and between-generational responses at the individual gene level. Among 372 and 445 DMR genes from Heron and Palm Island populations, respectively, only 12 genes were differentially methylated in both populations (Figure 4). These include genes involved in regulation of cell migration (slit2), cell volume homeostasis (znf706), cytoskeleton organization (dtnb), DNA replication (rfc4), immune response (btn1a1), nutrient metabolism (grpr, ren, spon2, and srebf1), sodium channel gating (unc80), and transcriptional regulation (rax2). Interestingly, all 12 genes that were differentially methylated in at least one of the Heron Island within-generational treatments were also differentially methylated in the Palm Island transgenerational or step treatments compared to control or within-generational treatments.

Fish from Heron Island phenotypically adjusted aerobic scope by reducing RMR without a significant change in MMR (Donelson and Munday, 2012). In contrast, the F2 Palm population did not achieve within-generational plasticity, but was able to do so when both parents and offspring were exposed to increased temperatures, either in the step treatment by decreasing RMR or in the transgenerational treatment by increasing MMR and decreasing RMR (Donelson et al., 2012; Bernal et al., 2018). This implies that there could be different physiological mechanisms and underlying genomic responses for within- and between-generational plasticity of aerobic scope, as well as between populations, following exposure to predicted future temperatures. Indeed, we found that each population might employ a different epigenetic toolkit in order to plastically adjust to warmer conditions. Only 12 DMR genes were common among 372 and 445 DMR genes from Heron and Palm Island fish, respectively. However, these 12 genes might be linked to the phenotypic reductions in RMR of both populations. Independent GO analysis for each population also suggested distinct functions may be employed for within- and between-generational plasticity.

These 12 genes shared by the two populations are involved in various regulatory functions for cell migration (slit2) (Zeng et al., 2018), cell volume homeostasis (znf706) (Dossena et al., 2011), cytoskeleton organization (dtnb) (Veroni et al., 2007), DNA replication (rfc4) (Kim and Brill, 2001), immune response (btn1a1) (Smith et al., 2010), nutrient metabolism (grpr, ren, spon2, and srebf1) (Persson et al., 2002; Ferre and Foufelle, 2010; Zhu et al., 2014; Silva et al., 2017), sodium channel gating (unc80) (Lear et al., 2013), transcriptional regulation (rax2) (Wang et al., 2004), and one had unknown function (dytn).

The overlapping differential methylation in genes associated with nutrient delivery and metabolism is not unexpected since these are known to be critical physiological attributes for survival of ectotherms in warming conditions (Portner and Knust, 2007; Portner and Farrell, 2008). The gene ren was differentially methylated in both populations (hypomethylated in Heron Island +1.5°C fish vs. controls, hypermethylated in Heron Island +3°C fish vs. +1.5°C fish, and hypomethylated in Palm Island transgenerational fish. vs. +3°C fish) and is particularly interesting as the encoded protein, RENIN, is the master hormone mediating the volume of extracellular fluid and blood pressure (Persson, 2003). As a component of the renin-angiotensin system, this protein catalyzes angiotensinogen released from the liver, which is eventually converted into angiotensin II, the vasoactive peptide (Persson, 2003). Renin-angiotensin system is also transcriptionally regulated in liver to control glucose and lipid level through hepatic insulin signaling (Silva et al., 2017). We speculate that epigenetic control on this protein may be a key mechanism to improve blood flow and nutrient homeostasis in order to plastically adjust to warmer conditions in A. polyacanthus. The grpr (gastrin-releasing peptide receptor) gene was hypermethylated in Heron Island +3°C fish vs. +1.5°C fish, hypermethylated in Palm Island transgenerational fish vs. +3°C fish, and hypermethylated in Palm Island transgenerational fish vs. step fish. Its encoded product, GRPR, mediates insulin secretion in pancreatic islet cells (Persson et al., 2002) and regulates glucose metabolism by modulating components of the aerobic glycolysis pathway (Rellinger et al., 2015). A recent mice model study suggested that GRPR may play an important role in hepatic ischemia injury (a status characterized by an insufficient supply of oxygen and nutrients), although the underlying molecular mechanism is yet unclear (Guo et al., 2019). SPON2 (spondin-2), whose encoding gene was hypomethylated in Heron Island +3°C fish vs. +1.5°C fish and also hypomethylated in Palm Island transgenerational fish vs. control, controls hepatic lipid metabolism and insulin resistance in rodents (Zhu et al., 2014). Furthermore, SPON2 is known to protect from cardiovascular diseases (Yan et al., 2011) and function in initiation of innate immunity (He et al., 2004) in rodent studies. SREBF1 (sterol regulatory element-binding protein 1), whose encoding gene was differentially methylated in both populations (hypermethylated in Heron Island +1.5°C fish vs. controls, hypomethylated in Heron Island +3°C fish vs. +1.5°C fish, hypomethylated in Palm Island step fish vs. control, and hypomethylated in Palm Island transgenerational fish vs. control), is the master regulator of the lipogenesis in liver in response to environmental signals such as nutrient demand and mediating the transcription of glycolytic and lipogenic genes in insulin- and sterol-dependent manner (Ferre and Foufelle, 2010; Shao and Espenshade, 2012). Epigenetic control on these regulatory genes involved in cell metabolism might be critical for the strict maintenance of homeostatic status at high temperatures.

There is insufficient data available to understand the likely role in thermal plasticity for the other genes that were differentially methylated in both populations (slit2, btn1a1, dtnb, dytn, rax2, rfc4, unc80, and znf706). The slit2 (slit guidance ligand 2) gene was differentially methylated in treatments where temperature was increased +1.5°C within-generations (Heron Island +1.5°C fish vs. controls and Palm Island step fish vs. control). In mammals, the encoded protein, SLIT2, modulates the migration of hepatic stellate cells, the major cell type involved in fibrogenic activity in the liver in chronically abnormal states such as hepatic fibrosis caused by metabolic disorder (Chang et al., 2015; Zeng et al., 2018). The remaining seven genes are not well-studied in relation to stress responses and/or their function in the liver. Although the role of these genes needs to be investigated further in fish liver following within- and between-generation exposure to increased temperatures, our results suggest that they may be useful as epigenetic markers for thermal plasticity.

The epigenetic patterns of within- and between-generational plasticity for the Heron and Palm Island populations of A. polyacanthus could be partially caused by different inheritance mechanisms for DNA methylation. Compared to the epigenomes for within-generational plasticity, which are only transmitted through somatic cell divisions, the epigenomes for between-generational plasticity include transmission through germ and somatic cell division. In mammals, the latter processes include epigenetic reprogramming events during gametogenesis and early embryogenesis, in which DNA methylation marks on the genome are erased and remodeled (Morgan et al., 2005; Messerschmidt et al., 2014). Although not extensively investigated in fish, epigenetic reprogramming differs by fish species. For example, zebrafish did not show global demethylation and recovery, while medaka showed similar reprogramming patterns to mammals (Jiang et al., 2013; Ortega-Recalde et al., 2019; Wang and Bhandari, 2019a, b). If complicated reprogramming event occurs in A. polyacanthus, modifications of epigenetic information during reprogramming (e.g., incomplete erasure and recovery) are possible for Palm Island fish (Kearns et al., 2000; Jacob and Moley, 2005). Furthermore, the exposure of Palm Island fish to increased temperatures for longer time frames compared to Heron Island fish (two vs. one generations, respectively) may increase the chances that the genome could be subjected to modification of DNA methylomes, such as replication-dependent passive demethylation or epigenetic drift that can dilute epigenetic memory (Guo et al., 2014; Messerschmidt et al., 2014; Ciccarone et al., 2018). In future studies, it would be of interest to delve more deeply into the cellular mechanisms of methylation inheritance and reprogramming in A. polyacanthus and other species that show differing capacities for within- and between-generational plasticity.