Heat and Hypoxia Exposure Mediates Circadian Rhythms Response via Methylation Modification in Apostichopus japonicas

As global warming progresses, heat and hypoxia are gradually becoming important factors threatening the survival, reproduction, and development of marine organisms. To determine the effect of heat and hypoxia on Apostichopus japonicus, whole genome methylation of the respiratory tree was determined under heat, hypoxia, and heat-hypoxia conditions [designed as heat stress treatment (HT), hypoxia treatment (LO), and heat-hypoxia combined treatment (HL) groups]. The number of differentially methylated regions (DMRs) under three treatments was determined based on the Venn diagram. The network of the DMRs associated with promoters that were co-existed under the three conditions showed that circadian rhythm was involved based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. Circadian rhythm-related genes, CRY1a, CRY1b, CLC, and TIM, decreased in LO and HL groups, while CRY1a, CRY1b, and BMAL1 increased in the HT group. Bisulfite sequencing PCR (BSP) showed that the methylation levels of CpG island regions in the promoters of CRY1a and CRY1b were upregulated in HT, LO, and HL groups, leading to the decreased promoter activity of CRY1a and CRY1b. RNAi of CRY1a and CRY1b led to increased enzyme activities of two energy-related enzymes, pyruvate kinase (PK) catalyzing the rate-limiting step in glycolysis, and ATPase hydrolyzing ATP to ADP, which were also increased under the three tested conditions. Thus, it was concluded that A. japonicus may respond to the heat, hypoxia, and heat-hypoxia stresses via the DNA methylation of heat, hypoxia, and heat-hypoxia stresses via the DNA methylation of CpG islands of circadian rhythm-related genes, which increased the activity of energy-related enzymes.


INTRODUCTION
Heat stress is a negative factor affecting the production and health of cultured animals, endangering the development of animal husbandry, thus causing serious economic losses and its impacts have been deteriorated due to global warming (St-Pierre et al., 2003;Quinteiro-Filho et al., 2010). In practice, heat stress tends to accompany hypoxia stress commonly, and the increased temperature is always following decreased dissolved oxygen (DO) (Parthasarathy et al., 1992). Therefore, heat and hypoxia are two related stressors that can threaten the cultured animals alone or in synergy (McBryan et al., 2013;Crozier and Hutchings, 2014;Jenny et al., 2016). Beyond a specific critical temperature, the organisms can expand their passive caloric range by increasing the anaerobic metabolic capacity, inducing molecular protective mechanisms and minimizing metabolic costs through metabolic inhibition (Pörtner et al., 2017). When organisms are continuously exposed to heat and hypoxia, the organisms exhibit many physiological and molecular responses to cope with these stresses, one of which is heat acclimation (Carter et al., 2005;Hung et al., 2005;Magalhães et al., 2010).
Heat acclimation is regulated by epigenetic factors from previous work, such as DNA methylation and histone methylation (Horowitz, 2016). DNA methylation regulates cellular response by transferring methyl to the C5 position of cytosine to form 5-methylcytosine (m5C) with no change in sequence. DNA sequence can be modified and reversed by intracellular and extracellular stimuli such as heat, hypoxia, pollutants, and nutrients (Baylin, 2006;Byun and Baccarelli, 2014). Previous researches have shown that the DNA methylation profiles changed under heat stress (Horowitz, 2016). DNA methyltransferases (DNMTs), such as DNMT1 and DNMT3, have been proven in response to heat stress before the formation of heat acclimation, indicating a relationship between methylation modification and heat stress (Dai et al., 2017(Dai et al., , 2018. The pieces of evidence have shown that changes of CpG methylation profile in HSP70 distal promoter region in the chicken and high mitochondrial calcium content in Rattus norvegicus may respond to heat acclimation (Assayag et al., 2012;Kisliouk et al., 2017). Data on mitochondrial DNA methylation in human diseases suggest that DNA methylation may play a potential role under the harmful environmental stress (Iacobazzi et al., 2013). Epigenetic-related metabolites produced by mitochondria play an essential role in processes of epigenetic and nuclear transcription, such as histone modification, chromatin remodeling, and nucleosome localization (Shaughnessy et al., 2014). In general, DNA methylation can help organisms to adapt to environmental changes through different adjustment methods.
Apostichopus japonicus are rich in nutrition and medicinal value and distributed across the northwest coast of the Pacific Ocean (Yang et al., 2005;Yuan et al., 2006). The echinoderms are more sensitive to environmental stress than other marine organisms, such as ascidians, mollusks, and anthozoans (Riedel et al., 2012). The heat and hypoxia stresses in summer far exceeded the optimal conditions for the growth and survival of A. japonicus, resulting in high mortality (Wang et al., 2015;Xu et al., 2016). To see the response of A. japonicus through DNA methylation, the whole-genome DNA methylation of the respiratory tree of A. japonicus under heat, hypoxia, and heathypoxia stresses were determined. Circadian rhythms drive 24h cycles in the physiology and behavior of organisms, and it enables organisms to synchronize their internal biology with their external environment to maximize fitness and survival (Yerushalmi and Green, 2009;Martinez-Bakker and Helm, 2015). Many studies have shown that changes in abiotic environmental conditions, such as heat, hypoxia, and toxic chemicals can affect circadian rhythms (Pittendrigh and Caldarola, 1973;Lim et al., 2006;Claudel et al., 2007;Svendsen et al., 2014). In the present study, the methylation changes of promoters of the circadian rhythm-related genes under the heat, hypoxia, and combined stresses, and their probable function were further explored.

Animals
The A. japonicus in the experiment were collected from the Shandong Oriental Ocean Technology Co., Ltd., of Yantai, China, and the A. japonicus is 9.5-11 cm in length and 3.5-4.5 cm in width with a wet weight of 100 ± 5 g. After being weighed, A. japonicus were kept at 16 ± 0.4 • C with 5 mg/L DO for a week. We used 5 mg/L DO because it is the level of DO in seawater with aeration at 26 • C, which was coincided with the fact that DO in the living aquarium of the A. japonicus in summer is maintained at 5-6 mg/L. During acclimation and experiment, the remaining feed was removed daily. For the stress conditions, 26 • C was used for heat stress and 2 mg/L DO was used for hypoxia stress, for the fact that the A. japonicus experienced the limits of temperature at 26 • C and 2 mg/L DO in summer in its important local living environment of northern Yellow sea and Bohai sea, as well as the definition of hypoxia by the committee on Environment and Natural Resources at the National Science and Technology Council in 2000 (Huang et al., 2012;Liu et al., 2014;Huo et al., 2019). Then, A. japonicus were divided into four groups. The standard control (NC) group was maintained in seawater with 5 mg/L DO at 16 • C. The hypoxia treatment (LO) group was maintained in seawater with 2 mg/L DO at 16 • C. The heat stress treatment (HT) group was maintained with 5 mg/L DO at 26 • C. The heat and hypoxia combined treatment (HL) group was maintained with 2 mg/L DO at 26 • C.
Under the experimental conditions, the environmental temperature was gradually increased from 16 to 26 • C at the rate of approximately 2 • C/h using a 100 W heating rod. The environmental DO level was reduced from 5 to 2 mg/L at the rate of approximately 1 mg/L/h to simulate environmental stress via air-N 2 gas flow adjustment system, which fills the seawater with oxygen and nitrogen continuously with real-time DO probe detection (Cao et al., 2019). The initial time is defined as the time when the stress group reaches the specified stress condition.

Histology Observation
The respiratory tree, as the unique respiratory organ of A. japonicus, plays an important role in respiration and metabolism (Huo et al., 2017). Thus, in this experiment, respiratory tree was chosen to study the mechanism under heat and hypoxia stress. The respiratory tree of A. japonicus from each group (NC, HT, LO, and HL) were collected, which were put into paraformaldehyde solution, then in 70% ethanol. The samples were dehydrated in ethanol at a graded series of 70, 75, 85, 95, and 100%, and then embedded in paraffin after rinsing with xylene. The tissue was cut into 7 mm sections and then stained with H&E (Xu et al., 2015). The tissue sections were observed with light microscopy (Nikon Corporation, Japan).

DNA Extraction and Preparation
After exposure to heat, hypoxia, and heat-hypoxia stress for 48 h, the respiratory trees of six A. japonicus from each group (NC, HT, LO, and HL) were taken and mixed in liquid nitrogen. Genomic DNA was isolated from A. japonicus respiratory tree using a TIANamp Genomic DNA kit (Tiangen Biotech, Beijing, China). DNA quality was detected by 1% agarose gel electrophoresis.

Library Preparation and Quantification
A total of 5.2 µg of genomic DNA was fragmented to 200-250 bp by sonication using a Bioruptor (Diagenode, Belgium), then the blunt-ending, dA addition to 3 ′ -end, and adaptor ligation. Bisulfite converted from ligated DNA was treated with EZ DNA Methylation-Gold Kit (ZYMO Research, CA, USA) which made the template DNA Illumina HighSeq4000 (BGI, Shenzhen, China) was used for the sequencing of ligated DNA. Sequence data have been deposited at the NCBI BioProject database under accession PRJNA748843.

Methylomic Data Analysis
After sequencing, the raw reads were analyzed by BGI programs by removing contamination, adaptor sequences, and low-quality reads to get the clean data. The clean data were mapped to the genome of A. japonicus by whole genome bisulfite sequence mapping program (BSMAP), followed by merging the mapping results and removing the duplication reads (Xi and Li, 2009). The methylation levels in the NC, HT, LO, and HL groups were detected by covering each mC by the total reads covering that cytosine.

Differentially Methylated Region Analysis
The differentially methylated regions (DMRs) were identified by comparison of the treatment groups and the NC group in A. japonicus using windows that contained at least five CpG (CHG or CHH) sites with a two-fold change in methylation level and Fisher's test P ≤ 0.05. The degree of difference of a methyl-cytosine (mCG, mCHG, and mCHH) was calculated while comparing the methylation level of DMRs in different groups by CIRCOS (Krzywinski et al., 2009).
The DMR-related genes were mapped to gene ontogeny (GO) terms in the database for the GO enrichment analysis (http:// www.geneontology.org/) to find the significantly enriched GO terms by the hypergeometric test. The Kyoto Encyclopedia of Genes and Genomes (KEGG) identifies significantly enriched metabolic pathways or signal transduction pathways, and the calculating formula is the same as that in GO analysis (Kanehisa et al., 2008).

Bisulfite Treatment of Genomic DNA
Genomic DNA of each group (NC, HT, LO, and HL) was treated with an EZ DNA methylation Gold TM kit (ZYMO Research, California, USA), and the unmethylated cytosine (C) in the sequence was transformed into uracil (U). In the process of PCR amplification, all uracil (U) was transformed into thymine (T), and the bisulfite sequencing PCR (BSP) products were sequenced to determine whether the CpG site was methylated. All the primers used for BSP are shown in Table 1.

Methylation Modification and Double Luciferase Experiment
To study the effect of methylation level in promoter region on the promoter activity, the CpG free vector was used as described previously (Klug and Rehli, 2006). All the primers used for plasmid construction are shown in Table 1. The sequence of all methylation sites of putative cryptochrome-1a (CRY1a) was screened by BSP and was further used to construct the recombinant plasmid CpGCRY1apro-1. The sequence of all the methylation sites of putative cryptochrome-1b (CRY1b) and the sequence of 84 methylation sites were screened by BSP and were further used to construct the recombinant plasmid CpGCRY1bpro-1 and CpGCRY1bpro-2, respectively. CpGCRY1apro-1, CpGCRY1bpro-1, and CpGCRY1bpro-2 were separately transfected into EPC cells and tested for the promoter activities using the dual luciferase test to the protocol of the CpG Methyltransferase (M.SssI) kit (Thermo Fisher Scientific, MA, USA).

Real-Time Quantitative PCR
The respiratory trees of A. japonicus from each group (NC, HT, LO, and HL) were taken with three biological replicates after exposure to stress for 0, 6, and 18 h. After that, RNA was extracted to detect the expression of circadian According to the TaKaRa M-MLV RTase cDNA synthesis kit, 1 µg of RNA was used for cDNA synthesis (TaKaRa Bio, Dalian, China). The β-actin of A. japonicus was used as the internal standard. All the primers used for realtime quantitative PCR (qRT-PCR) are shown in Table 1. The baseline was set automatically by the software to maintain consistency. The expression levels of the circadian rhythmrelated genes were determined by the 2 − CT method (Livak and Schmittgen, 2001).

Enzyme Activity Measurement
Considering the fact that pyruvate kinase (PK) is the ratelimiting enzyme catalyzing the final step in glycolysis, converting phosphoenolpyruvate (PEP) to pyruvate while phosphorylating ADP to produce ATP and ATPase is a kind of enzyme that catalyzes the hydrolysis of ATP to ADP, the enzyme activities of PK and ATPase were measured because both of which could reflect the energy metabolism of the organism (Lunt and Vander Heiden, 2011). ATPase and PK enzyme activities were measured in the A. japonicus collected from groups of NC, HT, LO, and HL. After 24 h of treatment under heat, hypoxia, and heat-hypoxia stress, the respiratory tree was quickly placed into liquid nitrogen, and the cell lysate was obtained after grinding. The BCA Protein Assay Kit detected the protein concentration. PK activity and ATPase activity were detected using the pyruvate kinase assay kit and the ATPase assay kit (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). To see if the CRY1a and CRY1b gene were involved, dsRNA CRY1a and dsRNA CRY1b were applied. The dsRNA was obtained using the MEGAscript TM T7 Transcription Kit (Thermo Fisher Scientific, MA, USA) by primers whose production was about 500 bp, with the T7 promoter sequence added in the forward primer. The designed fragment did not overlap with the quantitative verification fragment. In the study, 50 µl CRY1a dsRNA or CRY1b dsRNA were injected into A. japonicus using a 1 ml syringe. The expressions of CRY1a and CRY1b after interference were performed to verify its success interference rate. The A. japonicus were divided into the CRY1a interference group and CRY1b interference group. After interference, the enzyme activity of ATPase and PK was measured after 24 h of treatment under heat, hypoxia, and heathypoxia stresses.

Statistical Analysis
Statistical analysis was performed using the two-tailed Student's t-test. One-way ANOVA was applied to determine the differences between the control and the experimental groups. Any significant difference relative to the control for each time point was indicated using * P < 0.05, * * P < 0.01, and * * * P < 0.001.

Changes in the Morphology of the Respiratory Tree Under Standard Control (NC), Heat Stress Treatment (HT), Hypoxia Treatment (LO), and Heat-Hypoxia Combined Treatment (HL) Groups
When the A. japonicus was reared in seawater in the groups of heat stress treatment (HT), hypoxia treatment (LO), and heat-hypoxia combined treatment (HL), the morphology of the respiratory tree changed dramatically. Using light microscopy, the structural changes in the respiratory tree were observed ( Figure 1A). The size of the muscular layer marked the overall outline of the respiratory tree. In the standard control (NC) group, the thickness of connective tissue was high, and the dendritic lining epithelium was well-distributed and well-organized. The atrophy of connective tissue and dendritic lining epithelium was observed in all three groups. The thickness of connective tissue decreased, and the degeneration of connective tissue was the most obvious in HL group among all the three treatment groups. The dendritic lining epithelium and the dendritic lining epithelium, distributing on the inner side of connective tissue, continued to be shrinking in all three treatment groups.

Genome-wide DNA Methylation Patterns of A. japonicus
To further understand the molecular mechanism after the heat, hypoxia, and heat-hypoxia stresses, the genome-wide DNA methylation patterns of A. japonicus were investigated. Approximately 4% of all genomic C sites were methylated. Methylation in the DNA of A. japonicus was found to exist in three sequence sites: CG, CHG, and CHH (where H is A, C, or T), and 25.522% CG, 0.585% CHG, and 0.597% CHH was methylated in the mapping reads. A higher rate of methylated sites was present in the structure gene regions of internal exon and internal introns. Among all the regions, the average level of methylation in promoter regions was the lowest (Figure 1B). The DNA methylation in the other three treatment groups (HT, LO, and HL) showed the same characterization in the control sample.

DMR Analysis
DMRs were detected to characterize the differences in wholegenome DNA methylation levels among NC and HT, LO, and HL groups, respectively. A total of 43,929, 39,820, and 43,901 DMRs were separately identified in the three treatment groups compared with the NC group, which were mainly related to the methylation of CG sites ( Table 2). DMRs could be divided into two contents, one kind located inside the structural gene regions (named DMRs associated with genes), and another kind located inside the promoter region (named DMRs associated with promoters). Approximately 19% of DMRs were identified as DMRs associated with the gene, and 8% of DMRs were identified as DMRs associated with the promoter (Figure 1C). The DMRs associated with genes and DMRs associated with promoters were analyzed by GO and KEGG. The GO analysis showed that DMRs associated with genes or promoters were significantly enriched in the biological process of cellular process, metabolic process, cellular component of membrane, membrane part and molecular function of binding, and catalytic activity in all the three groups of HT, LO, and HL (Figure 2). Kyoto Encyclopedia of Genes and Genomes analysis showed that these DMRs associated with promoters were enriched in a notch signaling pathway and circadian rhythm in the groups of HT, LO, and HL ( Figure 3A). The DMRs associated with genes were significantly enriched in ubiquitin-mediated proteolysis,  lysosome, and other glycan degradation in the groups of HT and HL, while significantly enriched in apoptosis and fatty acid metabolism in the groups of LO (Figure 3B).
Through the analysis of the Venn diagram of the DMRs in the groups of NC, HT, LO, and HL, 3,872 DMRs associated with genes and 1,019 DMRs associated with promoters were identified in all of the three treatments groups simultaneously ( Figure 1D). Then, further analysis of the network of signal pathways of the 3,872 DMRs associated with genes and 1,019 DMRs associated with promoters were performed through KEGG enrichment (Figures 4, 5). The network of DMRs associated with genes was enriched in a big category, such as endocrine resistance, Th1 and Th2 cell differentiation, thyroid hormone signaling pathway, and dorso-ventral axis formation, which may take functions with notch signaling pathway. The network of DMRs associated with promoters was mainly enriched into two categories. One category consisted of the notch signaling pathway and the mRNA surveillance pathway. The other category consisted of the circadian rhythm, the calcium signaling pathway, the oxytocin signaling pathway, and the apelin signaling pathway. The circadian rhythm was first enriched in KEGG in all the treatment groups of DMRs associated with promoters from the results above. The circadian rhythm was further enriched in the network of signal pathways via the KEGG enrichment of DMRs in all of the three treatments groups simultaneously selected from the Venn diagram. In addition, circadian rhythm took the functions as a center and regulated peripheral signaling pathways in the network of DMRs associated with promoters. So, circadian rhythm was chosen for further study.

Relative Expression of Circadian Rhythm-Related Genes Under Heat and Hypoxia
Five circadian rhythms-related genes were identified from the genome of A. japonicus: CRY1a, CRY1b, Putative aryl hydrocarbon receptor nuclear translocator like 1 (BMAL1), Clock (CLC), and Timeless (TIM). The expression of circadian rhythm-related genes was studied after heat, hypoxia, and heathypoxia stress for 0, 6, and 18 h (Figure 6A). At the HT group, BMAL1 increased significantly and peaked at 6 h (6.88fold of the control, P < 0.01). The expressions of CRY1a and CRY1b reached peaked at 18 h (2.56-and 2.63-fold, respectively, P < 0.05). There was no significant difference in CLC and TIM. In the LO group, the expression of CRY1a decreased continuously and reached the lowest level at 18 h (0.55-fold of the control, P < 0.05). However, the expression of CRY1b, CLC, and TIM decreased at 18 h, while BMAL1 had no changes. The expressions of CRY1a, CRY1b, BMAL1, CLC, and TIM were significantly decreased in the group of HL (P < 0.05).

Methylation of the Promoters of Rhythmic Genes
Among the five circadian rhythms-related genes, CRY1a, CRY1b, and CLC, BMAL1 possessed hypermethylated CpG island regions when identified through the whole genome DNA methylomics. There are four high CpG islands as follows: CRY1a promoter region (site 299378-site 299529), CRY1b promoter region (site 52083-site 52112), CLC promoter region (site 3142-site 3218), and BMAL1 exon 1 region (site 280973-site 281011). The methylation sites of CRY1a, CRY1b, BMAL1, and CLC in these regions were identified under heat, hypoxia, and heat-hypoxia stresses ( Figure 7A). The methylation level of CRY1a CpG island region showed significant differences under heat, hypoxia, and heat-hypoxia stresses compared with the group of NC. The probability of methylation was 30% in the HT group and about 10% in the LO group, and 10% HL group compared with the control group. The probability of methylation of CRY1b CpG island region did not change, except for the methylation of site 84 significantly increased in the group of HT, LO, and HL, with about 30% probability of occurrence. The methylation level of BMAL1 and CLC did not change significantly compared with the NC group. The results of the dual-luciferase assay showed that when CpGCRY1apro-1 or CpGCRY1bpro-1 was methylated, the promoter activity containing the CpG island region decreased significantly (P < 0.05, Figure 7B). In addition, the activity of FIGURE 3 | (A) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DMRs associated promoters of HT, LO, and HL groups compared with NC group. (B) KEGG pathway analysis of DMRs associated genes of HT, LO, and HL group compared with NC group. The abscissa represents the proportion of the number of DMRs belonging to a specific KEGG pathway to the whole number of differentially expressed genes, and the ordinate represents the specific KEGG pathway.

Function of CRY1a and CRY1b Under Heat, Hypoxia, and Heat-Hypoxia Stress
The response of aquatic species to heat and hypoxia is to enter a state of low metabolism, reduce metabolic rate, suspend development and reproduction, and survive in the critical ecological changes. The activities of two enzymes, PK and ATPase, which dominated the energy levels, were detected under the three tested conditions (Figure 6B). Compared with the enzyme activity of the normal A. japonicus in the NC group (39.4 U/mgprot), the PK activity in the HT group and HL group increased significantly (P < 0.05), but there was no significant difference in the LO group (P > 0.05). Compared with the ATPase activity of A. japonicus in the NC group (102.3 U/mgprot), the ATPase activity increased significantly in the HT group (141.3 U/mgprot, P < 0.001). In comparison, it decreased in the LO group (66.1 U/mgprot, P > 0.05), and there was no significant difference in the ATPase activity in the HL group.
To see if CRY1a and CRY1b were involved in the energy generation, CRY1a and CRY1b were interfered in vivo. The PK activity in the interference groups was significantly increased under all three conditions (P < 0.001). In the CRY1a interference group, the highest PK activity was found in the HT group (84.5 U/mgprot), followed by the HL and LO groups. In the CRY1b interference group, the highest PK activity was found in the HL group (122.4 U/mgprot), followed by the HT group and LO group. The PK activity in the dsRNA CRY1b group was higher than that in the dsRNA CRY1a group. Compared with the control group, the ATPase activity of the interference groups was also significantly increased (P < 0.001). Among the CRY1a interference groups, ATPase activity was highest in the HT group (883.6 U/mgprot). Among the CRY1b interference groups, ATPase activity was highest in the HL group (874.3 U/mgprot) ( Figure 6B).

DISCUSSION
From the results of the respiratory tree became smaller and shrunk with small diffuse debris suffering from heat and hypoxia stress, it suggested that the respiratory tree was sensitive to heat and hypoxia, and can well-reflect the changes of environmental stress on the respiratory tree (Huo et al., 2018). After analyzing various methylation patterns (CG, CHG, and CHH) of wholegenome DNA in the respiratory tree of A. japonicus, most DNA methylation occurred at CpG sites, which was consistent with the methylation happened in mammals (Lister et al., 2011;Tomizawa et al., 2011;Ziller et al., 2011). It was shown that 19% of DMRs associated with genes and 8% of DMRs associated with promoters might significantly influence the gene expression. DMRs localized in the internal sequence of the gene, such as internal exon and internal introns, which may affect gene expression via changes in chromatin structure or transcription efficiency (Lorincz et al., 2004;Klose and Bird, 2006;Suzuki and Bird, 2008). DMRs located in promoter regions which were detected in this study, the activity of CpG regions in promoters of CRY1a and CRY1b decreased modified by in vitro methylation. There are two explanations for the mechanism of DNA methylation maintaining promoter FIGURE 6 | (A) The relative expression of genes of circadian rhythm under heat stress, hypoxia stress, and heat-hypoxia stress. (B) The relative expression of CRY1a and CRY1b after interference and relative activity of PK and ATPase under heat stress, hypoxia stress, and heat-hypoxia stress for 24 h and the interference groups of dsRNA CRY1a and dsRNA CRY1b under heat stress, hypoxia stress, and heat-hypoxia stress for 24 h.
FIGURE 7 | (A) Methylation sites of CRY1a, CRY1b, BMAL1, and CLC under heat stress, hypoxia stress, and heat-hypoxia stress. The histogram represents the proportion of methylation and unmethylation, black represents methylation and white represents unmethylation. Dot graph: each row represents a single clone, 10 replicates. Each dot represents a CpG dinucleotide, the white dot represents unmethylated CpG, and the black dot represents methylated CpG. (B) Promoter activity of CRY1a, CRY1b, CpGCRY1apro-1, and CpGCRY1bpro-1 selected all methylated sequences from (A), and CpGCRY1bpro-2 selected the methylation sequence at site 84 from (A) (the sequence was screened by repeated experiments).
silencing. First, methylation changes chromatin structure and inhibits transcription factors and co-transcription factors, thereby reducing the gene expression. Second, DNA methylation can recruit methyl binding proteins (MSPs) to inhibit chromatin state through interaction-maintained promoter silencing (Campanero et al., 2000;Joulie et al., 2010;Deaton and Bird, 2011). The results of GO and KEGG of DNA methylation status of promoter and gene regions showed that notch signaling pathway, Th1-Th2 pathway, apelin pathways, and circadian rhythm were enriched, which were all regulated by circadian rhythm (Arjona and Sarkar, 2006;Henley et al., 2009;Cai et al., 2018).
Furthermore, Li discovered Klf2 and Egr1 exerted the effects through a clock gene-controlled process through genome sequence analysis under heat and hypoxia stress . It indicated that circadian rhythm might act as a manager to regulate numerous signaling pathways in response to heat and hypoxia in A. japonicus. Although global DNA methylation has been analyzed in the respiratory tree of A. japonicus by Guo (Guo et al., 2013) and Zhao (Zhao et al., 2015), only the differences in overall DNA methylation levels between the respiratory tree and other tissues were discussed without defined pathways analysis. This study is the first to systematically explore the DNA methylation under heat, hypoxia, and simultaneous heat-hypoxia stresses.
In this study, the expression of CRY1a, CRY1b, BMAL1, CLC, and TIM in HT, LO, and HL groups suggested that circadian rhythm may play an important role in A. japonicus response to environmental stressors. The expressions of CRY1a and CRY1b increased in the group of HT, which were consistent with the study that heat upregulated the expression of CRY (Ji et al., 2018) but were inversely proportional to the results of DNA methylation inhibiting the activity of CpG regions in our experiment. In comparison, decreased expression of CRY1a and CRY1b in the groups of LO and HL were consistent with results of DNA methylation inhibiting the activity of CpG regions in our experiment. Studies have shown no correlation between promoter DNA methylation and transcriptional activity of some genes (Grimm et al., 2013;Vucic et al., 2014;Yoo et al., 2015). Additionally, we only selected 200 bp CpG island regions for promoter activity detection, which cannot represent the whole promoter activity.
Moreover, promoter DNA methylation can also be involved in regulating the expression of other distal genes (Bell et al., 2011;Yoo et al., 2015). Then we compared the stress conditions of the HT group and groups LO and HL. The cause of this phenomenon may be hypoxia stress. The study showed that hypoxia stress increased HIF-1 α, which decreased CRY content to some extent (Dimova et al., 2019). The circadian rhythms may be destroyed in extremely harsh environments, such as heathypoxia. Therefore, DNA methylation is not a single negative regulatory relationship with gene expression, and we need to further study the mechanism of methylation modification on gene expression of CRY1a and CRY1b.
To study the relationship between circadian rhythm and energy metabolism, PK and ATPase enzyme activities in groups of NC, HT, LO, and HL were detected, which were consistent with the expressions of CRY1a and CRY1b. The expression level and enzyme activity increased in the HT group while decreasing in LO and HL groups. In addition, the activity of PK and ATPase increased significantly after interfering with CRY1a and CRY1b, which suggests that CRY1a and CRY1b may inhibit PK activity and ATPase to some extent. As lower basal metabolic rate in heat acclimation organisms, circadian rhythm may contribute to the formation of heat acclimation by regulating the metabolizing enzyme activity. Studies have shown that heat stress delays the circadian rhythm of subsequent activities for several hours under constant temperature and dark conditions, compared with flies without heat stress (Sidote et al., 1998). It is suggested that heat may cause the rhythm delay of CRY1a and CRY1b, leading to dynamically change in activities of PK and ATPase to satisfy the energy demand of the body or a result of metabolic compensation (Lemos et al., 2003).
In general, this study provided a comprehensive analysis of DNA methylation in A. japonicus, and circadian rhythm was determined to be involved in response to the heat, hypoxia, and heat-hypoxia stresses via DNA methylation. Therefore, the results of this study might provide new clues for deciphering the response of A. japonicus to the global warming signals via the circadian rhythm and methylation modification.

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/, PRJNA748843.