Genome-Wide DNA Methylation Signatures of Sea Cucumber Apostichopus japonicus during Environmental Induced Aestivation

Organisms respond to severe environmental changes by entering into hypometabolic states, minimizing their metabolic rates, suspending development and reproduction, and surviving critical ecological changes. They come back to an active lifestyle once the environmental conditions are conducive. Marine invertebrates live in the aquatic environment and adapt to environmental changes in their whole life. Sea cucumbers and sponges are only two recently known types of marine organisms that aestivate in response to temperature change. Sea cucumber has become an excellent model organism for studies of environmentally-induced aestivation by marine invertebrates. DNA methylation, the most widely considered epigenetic marks, has been reported to contribute to phenotypic plasticity in response to environmental stress in aquatic organisms. Most of methylation-related enzymes, including DNA methyltransferases, Methyl-CpG binding domain proteins, and DNA demethylases, were up-regulated during aestivation. We conducted high-resolution whole-genome bisulfite sequencing of the intestine from sea cucumber at non-aestivation and deep-aestivation stages. Further DNA methylation profile analysis was also conducted across the distinct genomic features and entire transcriptional units. A different elevation in methylation level at internal exons was observed with clear demarcation of intron/exon boundaries during transcriptional unit scanning. The lowest methylation level occurs in the first exons, followed by the last exons and the internal exons. A significant increase in non-CpG methylation (CHG and CHH) was observed within the intron and mRNA regions in aestivation groups. A total of 1393 genes were annotated within hypermethylated DMRs (differentially methylated regions), and 749 genes were annotated within hypomethylated DMRs. Differentially methylated genes were enriched in the mRNA surveillance pathway, metabolic pathway, and RNA transport. Then, 24 hypermethylated genes and 15 hypomethylated genes were Retrovirus-related Pol polyprotein from transposon (RPPT) genes. This study provides further understanding of epigenetic control on environmental induced hypometabolism in aquatic organisms.


Introduction
Changing environmental conditions is a challenge for virtually all organisms. It may be caused by predictable rhythms (e.g., daily, tidal, seasonal) or episodically (e.g., a sudden environmental change Quantitative RT-PCR was performed on an SYBR Green real-time PCR assay (SYBR PrimeScriptTM RT-PCR Kit II, TaKaRA) with an Eppendorf Mastercycler ep realplex (Eppendorf, Hamburg, Germany). Then, 40 S ribosomal protein S18, ß-tubulin, and NADH dehydrogenase were used as internal controls, as described previously [26]. In our study of reference gene stability [27], the expression of these reference genes was stable during aestivation in the intestine tissue of A. japonicus. Thermal cycling was as follows: 95 • C for 5 s, 40 cycles at 95 • C for 10 s, 60 • C for 20 s and 72 • C for 30 s. The 2 −∆∆CT method was used to analyze the expression levels, and the target genes' Cts were normalized using the geometric mean of these three internal control genes. The relative expression levels were shown as mean ± s.d., and the statistical significance was set at p < 0.05.

Genomic DNA Extraction and Whole-Genome Bisulfite Sequencing
Genomic DNA was extracted from intestine tissues (~25 mg) of A. japonicus using DNeasy Blood and Tissue Kit (Cat No. 69504, QIAGEN) following the manufacturer's protocol. Genomic DNA of intestine tissues from control and aestivation groups was sent to BGI (BGI Tech Solutions Co., Ltd., Shenzhen, China) for whole-genome bisulfite sequencing. Genomic DNA was firstly fragmented by sonication using a Bioruptor to a mean size of approximately 250 bp, followed by the blunt-ending, dA addition to 3'-end, and adaptor ligation. Ligated DNA was bisulfite converted using the EZ DNA Methylation-Gold kit (ZYMO, Irvine, CA, USA). Further, fragments with a length of 100-300 bp were selected and purified by the QIAquick Gel Extraction kit (QIAGEN). The products were then amplified by PCR. At last, we performed 100 bp paired-end sequencing on HiSeq 2000 Sequencing System. Sequence data have been deposited at the NCBI BioProject database under accession PRJNA643989.

Identification of Differentially Methylated Regions
The raw reads were filtered to remove the adaptors, contamination, and low-quality data (unknown bases more than 10%; the ratio of bases whose quality was less than 20). The clean bisulfite reads were mapped to the sea cucumber reference genome [28] using BSMAP [29], and the duplicated reads were removed. The BSMAP script was BSMAP -a filename_1.clean.fq.gz -b filename_2.clean.fq.gz -o filename.sam -d ref.fa -u -v 8 -z 64 -p 4 -n 0 -w 20 -s 16 -f 10 -L 125. The sam files were converted to bam files using scripts (samtools view -S -b -o filename.bam filename.sam; samtools sort -m 2000000000 filename.bam filename.sort.bam; samtools index filename.sort.bam). The mapping rate and bisulfite conversion rate of each sample were calculated. The CpG methylation clustering and Principal Component Analysis were generated by MethylKit [30]. The scripts we used were clusterSamples(meth, dist = "correlation", method = "ward", plot = TRUE) and PCASamples(meth). The methylation level was determined by dividing the number of reads covering each mC by the total reads [31], which equal the mC/C ratio at each reference cytosine [32]. The formula we used was Rm averge = Nm all Nm all +Nnm all . Nm represents the reads number of mC, while Nnm represents the reads number of non-methylation reads. Differentially methylated regions (DMRs) were identified by windows that contained at least 5 CpG (CHG or CHH) sites with a 2-fold change in methylation level and Fisher's test p ≤ 0.05. The heat maps were built to present the distribution of methylation levels in distinct genomic features [33]. Canonical DNA methylation profiles of the entire transcriptional units were divided into different functional elements [33] to study the methylation level changes.

Gene Ontology and Pathway Enrichment of DMRs
GO enrichment analysis provided GO terms that were significantly enriched in the list of DMR-related genes, and filtered DMR-related genes that correspond to specific molecular function, cellular component, and biological processes. The DMR-related genes were mapped to the GO database (http://www.geneontology.org/). The significant enriched GO terms were identified using the hypergeometric test, and the calculated p-value goes through Bonferroni Correction, taking corrected p ≤ 0.05 as a threshold. KEGG pathway enrichment analysis further identified significantly enriched metabolic pathways or signal transduction pathways and helped understand the biological functions of the DMR-related genes [34]. The calculation is the same as that in the GO analysis.

Expression of DNA Methylation Related Enzymes during Aestivation
Quantitative RT-PCR was used to specifically investigate the gene expression levels of DNA methyltransferases (DNMT1, DNMT3a, DNMT3b), Methyl-CpG binding domain proteins (MBD2/3, MBD4, MBD5, MBD6), and DNA demethylases (TET, TDG). Despite DNMT3b, all target genes were up-regulated during aestivation ( Figure 1). More specifically, DNMT1 and DNMT3a were highly expressed in aestivation groups with a fold change of around six and four times. DNMT3a showed no differences between control groups and aestivation groups. All MBD proteins showed higher expression levels in aestivation groups, with a fold change of around three times (MBD2/3), twice (MBD4), ten times (MBD5), and six times (MBD6). TET and TDG were both expressed with a 1.5 fold change in aestivation groups.
Genes 2020, 11, x FOR PEER REVIEW 5 of 16 hypergeometric test, and the calculated p-value goes through Bonferroni Correction, taking corrected p ≤ 0.05 as a threshold. KEGG pathway enrichment analysis further identified significantly enriched metabolic pathways or signal transduction pathways and helped understand the biological functions of the DMR-related genes [34]. The calculation is the same as that in the GO analysis.

Expression of DNA Methylation Related Enzymes during Aestivation
Quantitative RT-PCR was used to specifically investigate the gene expression levels of DNA methyltransferases (DNMT1, DNMT3a, DNMT3b), Methyl-CpG binding domain proteins (MBD2/3, MBD4, MBD5, MBD6), and DNA demethylases (TET, TDG). Despite DNMT3b, all target genes were up-regulated during aestivation ( Figure 1). More specifically, DNMT1 and DNMT3a were highly expressed in aestivation groups with a fold change of around six and four times. DNMT3a showed no differences between control groups and aestivation groups. All MBD proteins showed higher expression levels in aestivation groups, with a fold change of around three times (MBD2/3), twice (MBD4), ten times (MBD5), and six times (MBD6). TET and TDG were both expressed with a 1.5 fold change in aestivation groups.

DNA Methylation of Sea Cucumber during Aestivation
We conducted whole-genome bisulfite sequencing of sea cucumber intestine from aestivation groups and control groups. After quality control of filtering, a total of 1.39 billion clean reads were generated, consisting of 235.2 million, 221.9 million, and 252.4 million reads for each aestivation sample and 214.7 million, 223.2 million, and 240.5 million reads for each non-aestivation sample ( Table 2). The bisulfite conversion rate (%) of all sequencing libraries ranges from 99.82% to 99.90%. After read alignment, clean reads were mapped to the reference genome of sea cucumber with mapping rates ranging from 62.38% to 63.45% (Table 2).

DNA Methylation of Sea Cucumber during Aestivation
We conducted whole-genome bisulfite sequencing of sea cucumber intestine from aestivation groups and control groups. After quality control of filtering, a total of 1.39 billion clean reads were generated, consisting of 235.2 million, 221.9 million, and 252.4 million reads for each aestivation sample and 214.7 million, 223.2 million, and 240.5 million reads for each non-aestivation sample ( Table 2). The bisulfite conversion rate (%) of all sequencing libraries ranges from 99.82% to 99.90%. After read alignment, clean reads were mapped to the reference genome of sea cucumber with mapping rates ranging from 62.38% to 63.45% (Table 2).

Proportion of Methylation Contexts and mCs Distribution across Genomic Features
Hierarchical clustering analysis revealed that the samples were clustered according to their groups (aestivation groups and control groups) (Supplementary Material Figure S1). PCA analysis showed that the first principal component could clearly divide samples into aestivation groups and control groups (Supplementary Material Figure S2). The methylation level distribution of A. japonicus genome showed the methylome's overall characteristics (Supplementary Material Figure S3). The methylation levels of approximately 35% of all mCG were hypermethylated (methylation level >90%). However, only 10% of mCHH and mCHG were hypermethylated (methylation level >90%) compared with mCG. The methylated Cs mostly occur in the form of mCG, followed by mCHH and mCHG. The methylation level distribution of mC and mCG were alike. The methylated Cs mostly occur in the form of mCG; approximately 94% of all detected mCs (Table 3). Only 1.2% and 4% of detected mCs were mCHG and mCHH ( Table 3). The proportion of mCG was significantly lower in the aestivation group than that of the control group (average proportion in control: 94.841%, average proportion in aestivation: 94.689%, p = 0.029). The proportion of non-CpG methylation, mCHG showed no differences between control and aestivation groups (control: 1.202%, aestivation: 1.23%, p = 0.105), whereas proportion of mCHH was significantly higher in the aestivation group than that in the control group (control: 3.956%, aestivation: 4.082%, p = 0.019). The heatmap presents the methylation landscape in different genomic features (whole genome, CGI, downstream 2 kb, upstream 2 kb, mRNA) ( Figure 2 and Supplementary Material Figures S4-S7). CpG islands contained the highest numbers of CpG sites (approximately 10-20 CpG sites in a 200 bp window) compared with other genomic features. About 30% of CpG sites in CpG islands were hypermethylated (methylation levels >90%) in the heatmap (Figure 2 and Supplementary Material Figure S4-S7). The other genomic features (downstream 2 kb, upstream 2 kb, mRNA) generally contained 0-10 CpG sites in the 200 bp window. A higher proportion of CpG sites within mRNA and upstream 2 kb were hypomethylated (methylation levels < 10%) than CpG sites within downstream 2 kb. The DNA methylation patterns across the entire transcriptional units (upstream, first exon, first intron, internal exon, internal exon, last exon, downstream) at the whole-genome level can help study the changes of methylation levels in distinct functional elements (Figure 3 and Supplementary Material Figure S8). Methylation differences between CG and non-CpG methylation (CHG and CHH) are visible (Figure 3 and Supplementary Material Figure S8), as methylation levels of CG are higher than those of CHG and CHH across the entire transcriptional units. Another feature is a distinct elevation in methylation level at internal exons with clear demarcation of intron/exon boundaries during transcriptional unit scanning. The lowest methylation level occurs in the first exons, followed by the last exons and the internal exons.  The average methylation levels of the whole genome, CDS, intron, and mRNA were listed in Table 4. The average methylation levels of CG, CHG, and CHH at the whole genome levels were 27.55%, 0.33%, and 0.30% in the aestivation group, and 26.53%, 0.27%, and 0.23% in the control group. The average methylation levels of CG, CHG, and CHH within CDS were 43.57%, 0.42%, and 0.50% in the aestivation group, and 44.48%, 0.39%, and 0.52% in the control group. The average methylation than other genomic features [35]. The average methylation of CG showed no differences between control and aestivation groups within CDS (p = 0.453), intron (p = 0.92), and mRNA (p = 0.33). The methylation level of CHG showed a significantly increased within intron and mRNA in the aestivation group compared with the control group (intron: p = 0.012; mRNA: p = 0.004). The methylation level of CHH showed significantly increased within intron and mRNA in the aestivation group compared with the control group (intron: p = 0.008; mRNA: p = 0.008).  The average methylation levels of the whole genome, CDS, intron, and mRNA were listed in Table 4. The average methylation levels of CG, CHG, and CHH at the whole genome levels were 27.55%, 0.33%, and 0.30% in the aestivation group, and 26.53%, 0.27%, and 0.23% in the control group. The average methylation levels of CG, CHG, and CHH within CDS were 43.57%, 0.42%, and 0.50% in the aestivation group, and 44.48%, 0.39%, and 0.52% in the control group. The average methylation levels of CG, CHG, and CHH within intron were 19.26%, 0.35%, and 0.29% in the aestivation group, and 19.32%, 0.28%, and 0.24% in the control group. Within mRNA, the average methylation levels of CG, CHG, and CHH were 26.07%, 0.34%, and 0.29% in the aestivation group, and 25.35%, 0.27%, and 0.25% in the control group. The average methylation levels within CDS were the highest, suggesting that cytosine methylation primarily occurs mainly in genes. This result was consistent with other marine invertebrates, such as Pacific oysters. In oyster, CDS showed the highest methylation levels than other genomic features [35]. The average methylation of CG showed no differences between control and aestivation groups within CDS (p = 0.453), intron (p = 0.92), and mRNA (p = 0.33). The methylation level of CHG showed a significantly increased within intron and mRNA in the aestivation group compared with the control group (intron: p = 0.012; mRNA: p = 0.004). The methylation level of CHH showed significantly increased within intron and mRNA in the aestivation group compared with the control group (intron: p = 0.008; mRNA: p = 0.008).

Identification and Enrichment Analysis of Differential Methylated Regions
A total of 37,378 DMRs were identified among all detected sea cucumber samples. A total of 1393 genes were in hypermethylated DMRs, and 749 genes were in hypomethylated DMRs (Supplementary Material Tables S1 and S2). When focusing on promotors, 774 and 465 genes were associated with hypermethylated DMRs and hypomethylated DMRs (Supplementary Material Tables S3 and S4).
GO enrichment analysis of DMR-related genes provided significantly enriched GO terms corresponding to specific molecular function, cellular component, and biological process (Figure 4). The top enriched GO terms in the biological process are cellular process, metabolic process, and single-organism processes. The over-represented GO terms in the cellular component are cell, cell part, and organelle. In terms of molecular function, the top enriched GO terms are binding, catalytic activity, and transporter activity. KEGG pathway enrichment analysis indicated that the annotated genes within hypermethylated DMRs were enriched in metabolic pathways and the mRNA surveillance pathway ( Figure 5). In contrast, the annotated genes within hypomethylated DMRs were enriched in the mRNA surveillance pathway and RNA transport ( Figure 5). Genes 2020, 11, x FOR PEER REVIEW 10 of 16   A total of 44 hypermethylated genes and 29 hypomethylated genes were involved in the mRNA surveillance pathway (Supplementary Material Tables S3 and S4). Among genes associated with mRNA surveillance, many genes were RPPT genes (Retrovirus-related Pol polyprotein from transposon). A total of 24 hypermethylated and 15 hypomethylated RPPT genes were found involved in mRNA surveillance. Enriched differentially methylated genes associated with the mRNA surveillance pathway included Retrotransposable element Tf2 155 kDa protein (type 1 and type 3), Transposon Ty3-I Gag-Pol polyprotein, and Pro-Pol polyprotein. Interestingly, many other members of RPPT genes were also found to be differentially methylated. A total of 37 RPPT genes were in the hypermethylated gene list and 26 RPPT genes in the hypomethylated gene list. Promotors of 42 RPPT genes and 29 RPPT genes were associated with hypermethylated DMRs and hypomethylated DMRs.

Discussion
An increasing number of studies reveal that epigenetic controls contribute to the differential gene expression during environmentally induced hypometabolism and seasonal change [12,[36][37][38][39]. In A. japonicus, the transcriptional changes of methylation associated enzymes DNMT1 and MBD2/3 were up-regulated during aestivation [40,41]. In our study, DNA methylation modifiers, including DNA methyltransferases (DNMT1, DNMT3a), methyl-CpG-binding domain proteins (MBD2/3, MBD4, MBD5, MBD6), DNA demethylases (TET, TDG) were all up-regulated in aestivating groups (Figure 1). In thirteen-lined ground squirrels, the altered expression of MBDs and MECP2 (methyl CpG binding protein 2) expression levels accompanies different stages of hibernation [42]. In read-eared slider turtle, an overall increase in DNMT protein expression in liver and muscle was observed during anoxia-induced hypometabolism [43]. The significant changes in methylation-related enzymes indicate that DNA methylation machinery contributes to regulatory response to environmental induced hypometabolism.
In aquatic species, there is still limited research suggesting that DNA methylation has regulatory functions in genes involved in environmental-induced aestivation. In the present study, we reported a genome-wide examination of DNA methylation in the intestine of sea cucumber A. japonicus during thermal stress-induced hypometabolism. The mapping rates of clean reads approximately ranged from 62.38% to 63.45% (Table 2), which is above the average value reported in other species that bisulfite-converted DNA typically has a mapping efficiency closer to 30-50% [44]. However, the WGBS results obtained from this study are also consistent because bisulfite conversion reduces the complexity of the genomic sequence and reduces the ability of most computational programs to align sequences onto the reference genome.
Further DNA methylation profile analysis was also conducted in our present study across the distinct genomic features and entire transcriptional units. Interestingly, a significant increase in non-CpG methylation (CHG and CHH) was observed within the intron and mRNA regions in aestivation groups. This result suggested that non-CpG methylation may be more vital to regulate aestivation in A. japonicus than CpG methylation. A whole-genome methylation landscape of constant heat-stressed pigs showed similar proportions of methylation at CpG sites but showed significantly decreased percentages at non-CpG sites with control pigs [45]. Non-CpG methylation is controlled through DNMT3a [46] and DNMT3b [47]. DNMT3a was highly up-regulated in the aestivation groups, potentially contributing to higher non-CpG methylation levels. A distinct elevation in methylation level at internal exons was observed with clear demarcation of intron/exon boundaries during transcriptional unit scanning. The lowest methylation level occurs in the first exons, followed by the last exons and the internal exons. This result was consistent with previous findings in other marine invertebrates [48], suggesting that the first exons appear more prone to the mutagenic effects. In contrast, internal exons are more influenced by the regulatory impact of DNA methylation and may play a role in constructing alternative splicing.
To better understand the functional classification of DMRs during aestivation, GO enrichment analysis was conducted to compare aestivation groups and control groups. The top enriched GO terms of DMRs identified in the present study were mainly associated with the term "cellular process" and "metabolic process" which suggests that DNA methylation may have an essential role in regulating global transcription during aestivation-induced hypometabolism. The enriched pathway analysis played a crucial role in identifying the pathways most strongly affected by DNA methylation. The results indicated that the annotated genes within hypermethylated DMRs were enriched in metabolic pathways and the mRNA surveillance pathway. As aestivation led to a low metabolic rate, energy redistribution, and defense mechanisms enhancement [19], the expression patterns of many genes associated with metabolic pathways might be altered. The annotated genes within hypomethylated DMRs were enriched in the mRNA surveillance pathway and RNA transport. The mRNA surveillance pathway is a mechanism that ensures the quality of mRNAs and detects translational errors and degrades abnormal mRNAs [49]. The mRNA surveillance pathway is also required for oxidative stress tolerance, which inhibits the generation of abnormal proteins from translation error, NSD (non-stop decay) [50]. Another recent study revealed that an mRNA surveillance system contributed to E. coli cold shock adaptation program [51]. Therefore, the mRNA surveillance pathway might be one of the essential adaptive strategies through DNA methylation modification for aestivating sea cucumbers (A. japonicus).
KEGG enrichment analysis of the genes with significant methylation difference indicates that these genes were mostly enriched in the mRNA surveillance pathway. Interestingly, mRNA surveillance pathway was mainly composed of retrovirus-related Pol polyprotein from transposon (RPPT) genes, one of the most significantly expanded gene families in A. japonicus genome that predominately accumulated in LG01 [28]. A total of 24 hypermethylated and 15 hypomethylated retrovirus-related pol polyprotein from transposon (RPPT) genes were found involved in mRNA surveillance, which might be potentially important regulators for aestivation. These genes are members of retrotransposons, which regulate transcription levels of their adjacent genes by genes coming under the control of their promoter [52,53].
The regulation of retrovirus-related transposon has been observed in many other species in response to environmentally induced stress. The expression of RPPT 412 was up-regulated when snails suffered thermal stress [54]. RPPT genes, including RPPT 17.6, RPPT 297, and RPPT 412, were significantly up-regulated in response to early heat stress in corals [55]. Enriched differentially methylated genes associated with the mRNA surveillance pathway also included transposon Ty3-I Gag-Pol polyprotein and Pro-Pol polyprotein. Gag-Pol polyprotein genes were detected and regulated when exposed to heat and cold stresses in manila clam [56]. These studies all support our current findings.
Retrotransposons are characterized by typical long terminal repeats (LTRs), which can amplify themselves and insert them into the genome at target sites [57]. The presence of active promoters within retrotransposons, because of reduced DNA methylation, can affect the standard transcription, resulting in reduced mRNA levels or translation capacity, or disrupted transcription initiation [58]. Additionally, retrotransposons could alter the genetic information and restructure genomes through chromosomal rearrangements [59,60]. The transposition rates of retrotransposons could contribute to genome evolution [61]. Retrovirus-related Pol polyprotein from transposon genes are among the significantly expanded gene families in sea cucumber [28]. Adaptive expansion of the RPPT gene family in the A. japonicus genome might be a more efficient regulatory mechanism on global gene expression during aestivation through differential methylation.

Conclusions
Our findings provide a comprehensive, detailed picture of DNA methylation patterns in A. japonicus during environmental-induced aestivation. Most of the methylation-related enzymes were up-regulated during aestivation. Whole methylome analysis suggested that DNA methylation may have an essential role in regulating global transcription during environmental induced hypometabolism. Further DNA methylation profile analysis was also conducted across the distinct genomic features and entire transcriptional units. A significant increase in non-CpG methylation (CHG and CHH) was observed within the intron and mRNA regions in aestivation groups. A total of 1393 genes were in hypermethylated DMRs, and 749 genes were in hypomethylated DMRs. Differentially methylated genes were enriched in the mRNA surveillance pathway, metabolic pathway, and RNA transport. Many retrovirus-related transposon genes modified by DNA methylation might be potentially essential regulators for aestivation. This study provides further understanding of epigenetic control on environmental induced hypometabolism in aquatic organisms. However, explicit confirmation of molecular mechanisms requires future intensive research, which may eventually aid in understanding the complexity of marine species' temperature-induced aestivation mechanisms.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/9/1020/s1, Figure S1: Cluster analysis of similarities among methylation profiles of aestivation groups and control groups; Figure S2: Sample structure identified by principal component analysis (PCA) with the first 2 principal components; Figure S3: Methylation level distribution of mC; Figure S4: Heat maps show CpG density patterns of A2 sample; Figure S5: Heat maps show CpG density patterns of A3 sample; Figure S6: Heat maps show CpG density patterns of C2 sample; Figure S7: Heat maps show CpG density patterns of C3 sample; Figure S8: Canonical DNA methylation profiles of the entire transcriptional units; Table S1: Gene list annotated within hypermethylated DMRs; Table S2: Gene list annotated within hypomethylated DMRs; Table S3: Hypermethylated gene list involved in the mRNA surveillance pathway; Table S4: Hypomethylated gene list involved in the mRNA surveillance pathway.