Genome-wide DNA methylation and transcription analysis reveal the potential epigenetic mechanism of heat stress response in the sea cucumber Apostichopus japonicus

DNA methylation is an important epigenetic modification that regulates many biological processes. The sea cucumber Apostichopus japonicus often suffers from heat stress that affects its growth and leads to significant economic losses. In this study, the mRNA expression patterns and DNA methylation characteristics in the body wall of A. japonicus under heat stress were analyzed by whole-genome bisulfite sequencing (WGBS) and transcriptome sequencing (RNA-seq). We found that CpG was the main DNA methylation type, and heat stress caused a significant increase in the overall methylation level and methylation rate, especially in the intergenic region of the A. japonicus genome. In total, 1,409 differentially expressed genes (DEGs) and 17,927 differentially methylated genes (DMGs) were obtained by RNA-seq and WGBS, respectively. Association analysis between DNA methylation and transcription identified 569 negatively correlated genes in both DMGs and DEGs, which indicated that DNA methylation affects on transcriptional regulation in response to heat stress. These negatively correlated genes were significantly enriched in pathways related to energy metabolism and immunoregulation, such as the thyroid hormone signaling pathway, renin secretion, notch signaling pathway and microRNAs in cancer. In addition, potential key genes, including heat shock protein (hsp70), calcium-activated chloride channel regulator 1(clca1), and tenascin R (tnr), were obtained and their expression and methylation were preliminarily verified. The results provide a new perspective for epigenetic and transcriptomic studies of A. japonicus response to heat stress, and provide a reference for breeding sea cucumbers resistant to high temperatures.

DNA methylation is an important epigenetic modification that regulates many biological processes. The sea cucumber Apostichopus japonicus often suffers from heat stress that affects its growth and leads to significant economic losses. In this study, the mRNA expression patterns and DNA methylation characteristics in the body wall of A. japonicus under heat stress were analyzed by wholegenome bisulfite sequencing (WGBS) and transcriptome sequencing (RNA-seq). We found that CpG was the main DNA methylation type, and heat stress caused a significant increase in the overall methylation level and methylation rate, especially in the intergenic region of the A. japonicus genome. In total, 1,409 differentially expressed genes (DEGs) and 17,927 differentially methylated genes (DMGs) were obtained by RNA-seq and WGBS, respectively. Association analysis between DNA methylation and transcription identified 569 negatively correlated genes in both DMGs and DEGs, which indicated that DNA methylation affects on transcriptional regulation in response to heat stress. These negatively correlated genes were significantly enriched in pathways related to energy metabolism and immunoregulation, such as the thyroid hormone signaling pathway, renin secretion, notch signaling pathway and microRNAs in cancer. In addition, potential key genes, including heat shock protein (hsp70), calcium-activated chloride channel regulator 1(clca1), and tenascin R (tnr), were obtained and their expression and methylation were preliminarily verified. The results provide a new perspective for epigenetic and transcriptomic studies of A. japonicus response to heat stress, and provide a reference for breeding sea cucumbers resistant to high temperatures.

Introduction
Sea cucumber (Apostichopus japonicus) is a traditional health product with high nutritional and economic value, which has been an important aquaculture species and farmed in many Asian countries, including China, Japan, and Korea (Cao, 2014). As a typical temperate species, the growth of A. japonicus is often hindered by extremely high or low temperatures, but especially high temperatures (Li et al., 2019). Changes in temperature usually reduce the activity and food intake of sea cucumbers, resulting in reduced immunity, which can lead to disease and death (Huo et al., 2021;Li et al., 2021;Huo et al., 2022). Therefore, revealing the molecular mechanism of their sensitivity to temperature and pressure can expand our understanding of how marine animals adapt to environmental challenges.
With the development of high-throughput sequencing technology, RNA-seq has been widely used in aquatic animals. The transcriptome map of the body wall and intestinal tissue of A. japonicus was first constructed in 2011, which provided data for studying the regeneration mechanism of A. japonicus (Sun et al., 2011). Subsequent studies constructed the transcriptome maps of sea cucumbers at different development stages and during summer dormancy stage, and numerous genes related to development, growth, metabolism and immunity were identified (Du et al., 2012;Zhao et al., 2014;Zhou et al., 2016). At present, studies on gene expression patterns of A. japonicus under heat stress have mainly focused on muscle, intestine, and respiratory tree tissues (Xu et al., 2018;Li et al., 2019;Chen, 2020). However, research on the body wall, which is the target organ of skin ulceration syndrome (SUS) is relatively rare.
Epigenetics is a kind of hereditary variation that does not change DNA sequence, and its regulatory mechanism mainly includes DNA methylation modification, histone modification, chromosome remodeling, and non-coding RNA regulation. These modifications interact and jointly regulate genome function (Kulis and Esteller, 2010;Crabtree et al., 2020) and these epigenetic modifications are easily induced by the environment (Feil and Fraga, 2012). DNA methylation is an important part of epigenetics research and one of the modification methods of DNA. It participates in the regulation of gene expression, cell growth and development, and response to adversity and stress by affecting the spatial conformation, stability, and interactions with proteins of nucleic acids (Entrambasaguas et al., 2021).
In aquatic animals, environmental factors influence gene expression by affecting DNA methylation status, so that the biological activities and phenotypes of organisms can be changed to adapt to the environment . Studies have found that DNA methylation mediates the biological process and then affects sex changes in Dicentrarchus labrax (Kuhl et al., 2011). A study showed that the total methylation rate of genomic DNA in the scallop (Patinopecten yessoensis) was affected by acute temperature stress (Wu et al., 2016). Li et al. (2018) found that after salt stress treatment, the exon methylation level of the growth-related gene igf1 was negatively correlated with its liver expression level in Cynoglossus semilaevis. Studies have also found that DNA methylation plays an important role in regulating the response to adversity in sea cucumbers. Methylation-sensitive amplified polymorphism (MSAP) analysis revealed that methylation levels significantly increased in sea cucumbers with SUS (Gao et al., 2017). Li et al. (2017) found that the expression of the apparent modification-related genes (dnmt1, hdac3, and mll5) of sea cucumbers significantly changed at high temperatures. Zhao et al. (2015) found that sea cucumber methylation levels during summer dormancy were higher than during non-summer dormancy. Although some progress has been achieved in the study of transcription pattern and DNA methylation characteristics of sea cucumber under heat stress, the effect of epigenetic modification on the transcription process is still unclear.
In this study, whole-genome bisulfite sequencing (WGBS) and transcriptome sequencing were used to gain insights into the DNA methylation characteristics and transcription pattern of A. japonicus under heat stress. We analyzed the correlation between gene expression and DNA methylation, and obtained potential key pathways and functional genes. These findings provide the first insight into the epigenetic mechanism of heat stress response in the body wall tissue of A. japonicus.

Experimental design and tissue collection
The experimental sea cucumbers were from Qingdao Ruizi Group Co., Ltd., they had an average weight of 50 ± 2g and were healthy, active, and free from disease or infection. Sea cucumbers were randomly divided into the control (CS) group at 17°C and the heat stress (HS) group at 32°C (the median lethal temperature of sea cucumber according to a previous study by our lab; Zhang et al., 2022), with 50 sea cucumbers in each group. The initial water temperature of the HS group was 17°C, and the temperature was increased at a rate of 1°C per day until the water temperature reached 32°C. On the fifth day, when the temperature reached 32°C, three healthy and three SUS sea cucumbers were randomly selected from the CS and HS groups, respectively. The body wall tissues were sectioned, placed in 2 ml EP tubes and put into liquid nitrogen for quick freezing. All samples were stored at -80°C for transcriptome and genome-wide DNA methylation sequencing analyses.

RNA-seq library construction, sequencing, and data analysis
Total RNA was isolated using TRIzol reagent (Qiagen, Carlsbad, CA, USA) following the manufacturer's procedure. The extracted RNA samples from the CS and HS groups were analyzed using a LabChip GX Touch (Perkin Elmer, CA, USA), and samples with RIN number > 7 were selected for library construction. Finally, the cDNA library was sequenced with PE150 on the Illumina Novaseq 6000 by LC-Bio (Hangzhou, China).
Cutadapt 4.2 (Kechin et al., 2017) was used to remove adaptor reads. High-quality sequences were mapped to the sea cucumber reference genome (assembled by our lab, unpublished data) using HISAT2 (Pertea et al., 2016). The transcriptomes from the six samples were combined and a complete transcriptome was reconstructed using gffcompare (Pertea and Pertea, 2020). The expression levels were estimated for all transcripts using StringTie (Pertea et al., 2016) and analyzed for mRNA expression levels by calculating of Fragments per Kilobase Million (FPKM). The differentially expressed genes (DEGs) were performed with FDR < 0.005 and fold change > 2 by DESeq2 (Love et al., 2014).

DNA methylation library construction, sequencing, and data analysis
Total DNA was isolated using the QIAamp Fast DNA Tissue Kit (Qiagen, 51404). The A260/280 ratio was read by a spectrophotometer and DNA was available when the A260/280 ratio was in the range of 1.8 to 2.0. DNA samples were subjected to bisulfite conversion after sonication using the Bisulfite Conversion Kit (Swift, MI, USA) to construct the WGBS library. The library was then sequenced on an Illumina Hiseq 4000 platform by LC-Bio (Hangzhou, China).
Cutadapt and perl scripts were used to remove the adapter contamination, low quality, and undetermined bases. Reads that passed quality control were mapped to the reference genome of sea cucumbers using Bismark (Krueger and Andrews, 2011). After alignment, reads were further deduplicated using SAMtools (Li et al., 2009). Plink (V1.90b6.21) software was used for principal component analysis (PCA). PerScript and MethPipe (Song et al., 2013) were used to locate each cytosine in the alignment sequence. The DNA methylation level was determined by the ratio of the number of reads supporting C (methylation) to the total number of reads (methylated and unmethylated). DSS software was used to calculate the differentially methylated regions (DMRs) with 1,000 bp slide windows, 500 bp overlap, P < 0.05 (from TSS to TES) (Akalin et al., 2012). In this study, genes were considered to be differentially methylated genes (DMGs) when at least one DMR was located within 3000 bp upstream of the start codon.

Association analysis of RNA-seq and WGBS
Pearson's method was used to test the relative expression data of DEGs and DMRs, and to determine the correlation between gene expression and DNA methylation. The selected genes were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment anslyses to characterize their function and response pathways.
Methylation sites were verified by BSP. The DNA of six cucumbers was modified with bisulfite using an EZ DNA Methylation-Gold ™ Kit (Takara, Japan). The modified DNA samples were subjected to BSP amplification using the Takara EpiTaqHS Kit (Takara, Japan). The six genes chosen for validation included hsp70, snrnp25 (small nuclear ribonucleoprotein U11/U12 subunit 25), gpatch2 (G patch domain containing 2), BSL78_04364, BSL78_07932 and BSL78_12691. The amplification cycling was performed as follows: 3 min at 98°C and then 10 s at 98°C, 30 s at 55°C, 30 s at 72°C, and 7 min at 72°C for 35 cycles. After the recombinant plasmid was constructed, positive clones were selected for sequencing analysis The qRT-PCR and BSP primers are listed in Table S1. The 2 -△△CT method was used to analyze the comparative mRNA expression levels. One-way analysis of variance conducted with SPSS 19.0 software (SPSS Inc., Chicago, IL, USA) was used for statistical analysis. The level of significance was set to P < 0.05. All data are presented as mean ± SD (standard deviation of the mean) (N = 3).

Transcriptome analysis
In total, 73.33G raw reads were generated from six RNA-seq samples. All raw reads were submitted to the SRA database under the accession number PRJNA901272. In this study, the Q20 ranged from 99.83% to 99.90%. An average of 75.03% of clean sequences could be mapped to the reference genome of A. japonicus, of which 40.14% to 45.70% had unique alignment and 24.99% to 26.75% had multiple alignments to the genome (Table 1).
Data analysis showed that the general distributions of FPKM values for six samples were similar; the CS group and the HS group had similar gene expression levels ( Figure S1). In total, 29,290 genes were detected from transcriptome sequencing and 1,409 DEGs were selected, including 528 up-regulated DEGs and 881 down-regulated DEGs (Figure 1). The information for the 1,409 DEGs is shown in Table S2.
The up-regulated DEGs were significantly enriched in 340 GO terms, and the down-regulated DEGs were significantly enriched in 308 GO terms (P < 0.05). In the GO terms of up-regulated DEGs, the following terms were frequently enriched in protein folding, response to heat and extracellular exosomes and stress-related terms ( Figure 2A). The down-regulated DEGs were mainly enriched in terms related to protein decomposition and metabolism, such as calcium ion binding, cellular protein catabolism, and hydrolase activity ( Figure 2B). KEGG enrichment showed that the up-regulated DEGs were significantly enriched in 18 KEGG signaling pathways, including lysine degradation, Toll and IMD signaling pathway, and protein processing in endoplasmic reticulum ( Figure 2C). The down-regulated DEGs were significantly enriched in 24 KEGG signaling pathways, including notch signaling pathway, ubiquitin mediated proteolysis, and MAPK signaling pathways ( Figure 2D).

DNA methylation characteristics
The DNA libraries of six samples were constructed and sequenced on the Hiseq 4000 platform. All raw reads were submitted to the SRA database under the accession number PRJNA901271. First, sequencing adapters and low-quality data were removed, and the average value of clean data was 24.04 Gb for each library. In this study, the Q20 value varied from 96.78% to 97.70%, and the genome mapping rate varied from 41.24% to 42.59% for each sample. The conversion of cytosine under the bisulfite treatment was between 99.35% and 99.66%, with an average of 99.55%, which indicated high efficiency of bisulfite treatment (Table 2). PCA analysis showed that both CS and HS groups had effective repliactions ( Figure S2).
The average DNA methylation levels of C sites and different types (CpG, CHG, and CHH) in the body wall tissues of sea cucumbers under heat stress are shown in Table 1. The whole genome methylation level of the 23 chromosomes of A. japonicas is also shown in a Circos plot ( Figure 3A). The percentage of methylated C sites to the total number of C sites was 3.94 ± 0.17% in the HS samples and 3.70 ± 0.03% in the CS samples, which indicated a significant increase in the percentage of methylated C sites in the genome under heat stress (P < 0.05). Statistical results of the average methylation level of different types (CpG, CHG, and CHH) showed that the proportions of mCpG, mCHG, and mCHH in the CS group were 24.14%, 0.56%, and 0.47%, respectively. While the mCpG mCHG and mCHH in HS group accounted for 24.22%, 0.98% and 0.76%, respectively. For the sea cucumber genome, the DNA methylation ratio of CpG sites was significantly higher than that of the CHG and CHH sites (Table 2). Overall, DNA methylation of the sea cucumber was mainly observed in CpG sites (84%) and rarely in CHG (13.6%) and CHH sites (3.8%). CpG DNA methylation levels varied to a greater extent across functional regions, whereas CHH and CHG methylation levels varied to a lesser extent (Table S3). The promoter and exon regions had the highest proportion of CpG methylation levels ( Figure 3B).
The DSS method was used to identify genome-wide DNA methylation changes between the CS and HS groups. In total, 748,347 DMRs were identified in the sea cucumber genome. Among them, 48,392 DMRs were in promoter, 69,090 DMRs were in exons, 129,779 DMRs were in introns, 481,072 DMRs were in the intergenic region, and 20,014 DMRs were in the CGI.CG island region. The results showed that the percentage of methylated Heatmap of DEGs between CS and HS with FDR < 0.005 and fold change > 2. genes in the intergenic region of sea cucumbers was the highest, accounting for 64.02%, whereas the percentage of methylated genes in the CGI.CG island region was the lowest, accounting for 2.72% ( Figure 3C). Overall, 17,927 DMGs were identified in DMRs. The GO enrichment analysis showed that DMGs were significantly enriched in substance transport-related terms, such as protein binding, actin binding, and calcium ion binding ( Figure 4A). In addition, the KEGG enrichment analysis showed that these DMGs were significantly enriched in RNA transport, MAPK signaling pathway, ubiquitin mediated proteolysis, and basal transcription factors ( Figure 4B).  Note: mC: Percentage of methylated C cites in the whole genome; mCpG: Percentage of methylated C cites in CpG context region to the total number of C sites; mCHG: Percentage of methylated C cites in CHG context region to the total number of C sites; mCHH: Percentage of methylated C cites in CHH context region to the total number of C sites. Data were at least three independent biological replicates. Asterisk (*) indicates significant differences between treatments (P <0.05).

Association analysis between DMRs and DEGs
To investigate the effect of DNA methylation on transcription, 808 shared genes were identified in both DMGs and DEGs ( Figure 5A). The results of association analysis of DMRs and DEGs are shown in Table S4. In total, 569 negatively correlated genes were identified, which accounted for 70.42% of shared genes, including 293 hyper-methylated and down-regulated genes, and 276 hypo-methylated and up-regulated genes. In addition, 223 negatively correlated genes with DMRs located in the promoter region and 276 negatively correlated genes with DMRs located in the gene body region were identified. The results showed that more genes had a negative correlation between transcription patterns and DNA methylation characteristics ( Figure 5B). The gene expression and methylation levels of these genes in the CS and HS groups are shown in Figure 5C. Enrichment analysis of these negatively correlated genes indicated that genes were mainly enriched in GO terms such as transmembrane transport, integrated component of membrane, and protein binding ( Figure 5D). KEGG enrichment results showed that genes were mainly enriched in the thyroid hormone signaling pathway, renin secretion, cGMP-PKG signaling pathway, notch signaling pathway, antigen processing and presentation, intestinal immune network for IgA production, and microRNAs in cancer. Some key functional genes were screened, including hsp70, clca1, tnr, bag3, traf7, and fut4 ( Figure 5E).

Validation analysis of key genes
The results showed that the expression levels of tnr, clca1, hsp70, act-3, bag3 and traf7 were up-regulated in sea cucumber under heat stress, and the expression levels of srpr, wdr20 and fut4 were downregulated in sea cucumber under heat stress. The qRT-PCR results showed that the expression of DEGs was consistent with the transcriptome sequencing results ( Figure 6A). In this study, six DMGs were randomly selected for BSP, and the WGBS results were verified. According to a standard regression analysis, the results were consistent between BSP and WGBS ( Figure 6B).

Discussion
RNA-seq has been widely used in studies of gene expression patterns of A. japonicus under heat stress on muscle, intestine, and respiratory tree tissues (Xu et al., 2018;Li et al., 2019;Chen, 2020). Li et al. (2019) found that 46 DEGs were involved in the innate and adaptive immunity of A. japonicus muscle tissue under heat stress, which contained many genes of heat shock family. Interestingly, genes of heat shock families such as hsp26, hsp70, hsc70 and hsp90 were selected in intestine tissue under heat stress (Xu et al., 2018). Chen (2020) also revealed that DEGs in the intestine were significantly enriched in immune-related GO terms after heat stress of A. japonicus, such as toll-like receptor 3, carboxyterminal kinesin 2-like, and complement component C3. Studies have also shown that A. japonicus can respond to heat stress by regulating the expression of immune-related genes with some noncoding RNAs such as miRNA, circRNA and lncRNA (Huo et al., 2020;Huo et al., 2021;Chang et al., 2022). In the present study, the down-regulated DEGs were significantly enriched in immunerelated pathways such as notch signaling pathway, ubiquitin mediated proteolysis, and MAPK signaling pathway. It was reported that these pathways can regulate the production of cytokines and participate in immune response (Meurette and Mehlen, 2018;Holdsworth et al., 2020;Liu et al., 2021). These results indicated that heat stress may induce a significant decrease in the immune regulation ability of sea cucumbers and provided an insight into the occurrence of SUS of A. japonicus under heat stress.
Research on aquatic organisms has shown that DNA methylation can be altered in response to temperature shift, diet composition, ploidy, and heterosis (Zhuo et al., 2019;Han et al., 2021;Yuan et al., 2021). Our results showed that methylation levels in the body wall of sea cucumbers increased after heat stress. Interestingly, a similar methylation pattern was also observed in the Yesso scallop (Patinopecten yessoensis), Pacific oyster (Crassostrea gigas) and Pacific abalone (Haliotis discus hannai), whose DNA methylation levels increased under heat stress (Huang et al., 2021;Wang et al., 2021;Yuan et al., 2021). In addition, previous studies showed that the DNA methylation levels in body wall with SUS were markedly higher than those of healthy body wall in A. japonicus (Gao et al., 2017;Sun et al., 2020). These results demonstrate that aquatic animals generally respond to environmental stress by increasing the DNA methylation level. However, a different result was obtained in a previous study in which the DNA methylation level of the A. japonicus digestive tractdecreased after heat stress at 32°C (Wen et al., 2021). We speculate that this divergence may be due to the difference in tissue type and selection of reference genome. Our results indicated that the body wall DNA methylated cytosines of sea cucumbers under heat stress were mainly located on CpG dinucleotides. Previous studies have also shown that DNA methylation in the body wall and intestinal tissues of A. japonicus mainly occurred at the CpG sites during Vibrio splendens infestation, aestivation, and heat stress Yang et al., 2020;Wen et al., 2021). In addition, the DNA methylation pattern in A. japonicus was consistent with research results on the Pacific abalone, Yesso scallop, Pacific oyster (Huang et al., 2021;Yuan et al., 2021;Venkataraman et al., 2022), and other metazoans (Lechner et al., 2013). Although CpG was the main type of DNA methylation, there were some differences in the characteristics of CpG in different regions of the genome. The DNA methylation level of the exon region in A. japonicus intestine during aestivation  was lower than that in our study. We speculate that the difference may be related to temperature selection, because temperature can affect the methylation level of exons (Iamjan et al., 2021). The DNA methylation characteristics of different tissues after heat stress might be different, and WGBS data are needed for more tissue types in the future. The relationship between methylation and gene expression was explored by analyzing the correlation between the expression level of DEGs and the methylation characteristics of DMRs. As a result, 808 shared genes, including 569 negative regulatory genes and 239 positive regulatory genes, were screened by DNA methylation and transcription association analysis; this indicated that DNA methylation mainly regulates gene expression through negativeregulation. This was consistent with previous studies on DNA methylation of Pacific abalone and European sea bass (Anastasiadi et al., 2017;Huang et al., 2021). Both negative and positive correlations between methylation levels and gene expression levels have also been found in the body wall and intestine of A. japonicus (Chen, 2020;Sun et al., 2020).The effect of methylation on gene expression often depends on the location in the genome. Generally, methylation of the gene body promotes gene expression, whereas methylation of the promoter inhibits gene expression (Palou-Maŕquez et al., 2021). However, our study showed that most of the negative regulatory genes were located in the gene body, and there were few negative regulatory genes in the promoter. The same result was found in a study on intestinal DNA methylation in A. japonicus under heat stress (Chen, 2020). This condition may be caused by re-methylation of CpG islands within the gene or inhibition of antisense transcription factors (Lou et al., 2014).
KEGG enrichment of 569 negative regulatory genes showed that they were mainly enriched in energy metabolism pathways, such as thyroid hormone signaling pathway, renin secretion and cGMP-PKG signaling pathway, as well as notch signaling pathway, antigen processing and presentation, intestinal immune network for IgA production, and microRNAs in cancer, which are related to immunoregulation. The cGMP-PKG signaling pathway can generally mediate the expression of cytokine PDE5 and affect the proliferation and apoptosis of cancer cells . The thyroid hormone signaling pathway is an essential signaling pathways for regulating growth, development and energy metabolism, and the downstream molecules can participate in the development and metamorphosis of marine invertebrates (Gilleron et al., 2006). The notch signaling pathway, as a classical immunerelated signaling pathway, has a wide and close relationship with the occurrence and development of tumors. The signal transduction between adjacent cells through the Notch receptor can regulate cell differentiation, proliferation, and apoptosis (Sajadimajd et al., 2022). Therefore, these pathways may be important adaptation strategies of sea cucumbers under heat stress.
The DNA methylation levels and expression levels of negative regulatory genes in these pathways were randomly verified by qRT-PCR and BSP. These negatively regulated genes are ideal targets in the future to reveal the epigenetic regulation of heat stress response in A. japonicus. For example, hsp70 can reduce the effects of toxic substance accumulation during stress (Radons, 2016). Furthermore, hsp70 may play a role in both innate and acquired immune systems (Treweek et al., 2015). Studies have found that the expression level of hsp70 significantly increased significantly with the increase of heat stress time. The rapid and persistent response of hsp70 indicates its critical role in the heat stress response of A. japonicus (Xu et al., 2015;Xu et al., 2016). The up-regulated expression of hsp70 in sea cucumbers under heat stress has also been detected in our previous study (Chang et al., 2022). Additionally, clca1 was found to be a key gene in IL-13 (Proinflammatory interleukin-13) dependent mucosal metaplasia (Alevy et al., 2012). In Cynoglossus semilaevis, clca1 acts as an activator of hepatic macrophages and regulates the expression of cytokines to inhibit pathogenic infection (Li and Sun, 2016). Tnr plays an important role in the development of the nervous system, the regeneration of tissues and organs after trauma and it can regulate the morphology of extracellular matrix (Meloty-Kapella et al., 2006). Another study showed that tnr could act as an important ECM component involved in neurodevelopmental regeneration and participates in the intestinal regeneration of A. japonicus (Yamazaki et al., 2020). We speculated that A. japonicus could repair tissue damage during heat stress and resist environmental changes by up-regulating the expression of tnr.

Conclusion
This study was the first attempt to investigate the gene transcription pattern and genome-wide DNA methylation characteristics in the body wall tissue of sea cucumbers in response to heat stress. We found significant changes in overall DNA methylation levels of sea cucumbers under heat stress, and obtained a list of DEGs and DMGs. In particular, many negatively correlated genes in both DMGs and DEGs were identified, which indicated that DNA methylation affects transcriptional regulation during heat stress. Further analysis showed that these negatively correlated genes were significantly enriched in pathways related to energy metabolism and immunoregulation. These results provide important targets for further revealing the response or adaptation mechanism of sea cucumbers to heat stress, and deepened our understanding of the epigenetic regulatory mechanism of heat stress response in sea cucumbers.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA901271 https://www.ncbi.nlm.nih.gov/, PRJNA901272.

Author contributions
ML and YW conceived and designed the experiments. MC, JG, XR, BL, XL, JW, ZZ, YY and CW performed the experiments. MC and JG analyzed the data. MC and JG wrote the paper. All authors read and approved the manuscript. All authors contributed to the article and approved the submitted version.