Transcriptome analysis reveals differentially expressed genes for improving survival of dormant Apostichopus japonicus in response to high temperature

Background: Aestivation is one of the strategies used by sea cucumbers (Apostichopus japonicas) in order to improve survival in response to the high-temperature and droughty conditions. Previous studies have carried out to investigate the immune or physiological alterations at the aestivation stage. However, it lacks information on the relationship between immunity and physiology. Herein, transcriptome sequencing was used to study gene expression during the aestivation stage. The results of this study provide a comprehensive understanding of the molecular mechanisms that protect sea cucumbers from the high-temperature condition, which favors improving survival in cultured sea cucumbers. Results: The transcriptome analysis of dormant (aestivation) and revival sea cucumbers generated 2,368 differentially expressed genes (down-regulation: 927; up-regulation: 1,441) and 39,081 unchanged genes. Basing on Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses, the down-regulated genes of dormant group were identified to be involved in DNA replication, RNA metabolic process, protein modification and biosynthesis, macromolecule metabolism, and cellular metabolism, which resulted in the inhibition of motility, skeletal development, neural activity, cell proliferation and development of A. japonicas. In contrast, the up-regulated genes were associated with fatty acid metabolism, carbohydrate hydrolysis, and phagocytosis. Protein-protein interaction network analysis further revealed that carbohydrate hydrolysis promoted the phagocytosis activity in the dormant group. Furthermore, the expression pattern of all tested genes in qRT-PCR analysis fitted well with those in RNA-Seq, with the exception of FASNL, which was unchanged in the qRT-PCR but up-regulated in RNA-seq. Conclusions: During the dormant stage, sea cucumbers decreased DNA replication, transcription and translation to achieve a hypometabolic state beneficial for reducing energy consumption. On the contrary, fatty acid metabolism and carbohydrate hydrolysis were increased for energy supply. Moreover, high levels of carbohydrate hydrolysis promoted phagocytosis, which is a crucial innate immune response to infection by pathogens. These results provided new insight into potential molecular mechanisms that enable the sea cucumbers to respond to high temperatures.

As shown in many studies, high-throughput RNA-Sequencing (RNA-Seq) technology is an excellent approach to study the transcriptional levels of large numbers of genes in tissues or organisms. The results of sequencing data can provide the "whole picture" of the gene expression profile in the organisms [18][19][20][21]. Obtained genes were annotated upon Protein family (Pfam), Superfamily, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database for characterizing their functions at the molecular and cellular levels [22][23][24][25]. In order to understand the potential molecular mechanism relating to the aestivation stage of sea cucumber, differentially expressed genes (DEGs) were analyzed to elucidate the functional categories based on the GO and KEGG enrichment analyses [22,23]. Moreover, Venn diagram analysis was used to study the special expression of genes in each group, and then performed GO and KEGG enrichment to analyze functional categories [5,22]. Therefore, the database of Search Tool for the Retrieval of Interacting Genes (STRING) was used to study functional interaction of proteins (including stable complexes, metabolic pathways and a bewildering array of direct and indirect regulatory interactions), which was also usually used to analyze the transcriptomic profile [19,26]. Therefore, our study employed transcriptome sequencing to investigate the alteration of physiological metabolism and immune in the aestivation stage, and then analyze the potential associations between them. The results of this study provide new insight into the potential molecular mechanisms of physiological metabolisms and immune responses, which are crucial in improving the survival of sea cucumber in aquaculture.

Transcriptome quantification
A total of 474,876,912 raw reads were acquired after RNA-seq containing 457,223,622 clean reads (96.28% of the raw reads) after removing adaptor sequences, reads containing poly-N and low-quality reads (Table 1). Moreover, 301,336,131 reads (65,90% of the clean reads) were mapped to the reference genome of A. japonicus which was composed of 279,721,413 (61.18% of the clean reads) uniquely mapped reads and 21,614,718 multiply mapped reads (4.73% of the clean reads). The average Q20 and Q30 of the read data was 97.43% and 92.90%, respectively. The QC content was accounted for 41.18% of the sequencing data.

Differential expression genes (DEGs) characterization
A total of 2,368 DEGs (including 927 down-regulated and 1,441 up-regulated genes) and 39,081 unchanged genes were acquired in the dormant group compared to the revival group (Fig. 1). In the dormant group, GO analysis revealed that the down-and up-regulated genes were annotated to subcategories under three primary categories (i.e. biological process, cellular component, and molecular function) (Fig. 2). The down-regulated genes were assigned to 41 sub-categories, which were mainly classified into DNA replication, RNA metabolic process, protein modification and biosynthesis, macromolecule metabolism, and cellular metabolism ( Fig. 2A). By contrast, the up-regulated genes were assigned to 85 sub-categories, such as peptidase activity, acting on L-amino acid peptides, peptidase activity, proteolysis, endopeptidase activity, metallopeptidase activity, serine-type endopeptidase activity, serine-type peptidase activity, and serine hydrolase activity (all belongs to protein degradation category). Moreover, the carbohydrate metabolic process was also significantly enriched in the dormant group (Fig. 2B).
To better understand the function of DEGs in the dormant group, an analysis of the KEGG enrichment using was performed (Fig. 3). The results showed that the down-regulated genes that related to the gene expression process (i.e. spliceosome, RNA transport, ribosome biogenesis in eukaryotes and ribosome) and proteasome metabolism pathways were significantly enriched (Fig. 3A). The upregulated genes were identified to associate with the pathways relating to phagocytosis (lysosome, endocytosis, and phagosome), lipid metabolism (fatty acid metabolism, steroid biosynthesis, ether lipid metabolism, and glycerophospholipid metabolism), and amino acid metabolism (glycine, serine and threonine metabolism, arginine and proline metabolism, tyrosine metabolism, and histidine metabolism) (Fig. 3B). Otherwise, the metabolisms relating to carbohydrates (including starch and sucrose metabolisms) were also enhanced in the dormant group (Fig. 3B).

Expression of genes between dormant and revival groups
In the Venn diagram analysis, a number of 1,850 or 405 genes were only expressed in the dormant or revival groups, respectively. All of these genes were further analyzed basing on the KEGG database.
The results revealed that the pathways involving in the carbohydrate metabolisms (i.e. starch and sucrose metabolism and other glycan degradation), fatty acid-related metabolism (fatty acid metabolism, biosynthesis of unsaturated fatty acids, and fatty acid elongation) and phagocytosis (phagosome, lysosome, and endocytosis) were significantly enriched in the dormant group compared to the revival group.

Protein-protein interaction (PPI) analysis of up-and down-regulated genes
The STRING database was used to analyze the interactions among different expression genes and the results on the characterization of genes in the PPI diagram, which were shown in Table S1A. The top core down-regulated genes were classified into the transcription and translation, such as methyltransferase NSUN6, ribosomal protein uS7, translation initiation factor eIF4A II, ribosome production factor 1, RNA cytidine acetyltransferase, U3 small nucleolar ribonucleoprotein protein IMP4, pre-rRNA-processing protein TSR1-like, Eukaryotic translation initiation factor 3 subunit B, GMP synthase, and eukaryotic translation initiation factor 6 ( In contrast, the interactions among the up-regulated genes were found to mainly relate to the carbohydrate hydrolysis (including maltase-glucoamylase, hexokinase-1, pancreatic alpha-amylase, alpha-amylase B, alpha-amylase) and the promotion of phagosome (V-type H + -transporting ATPase subunit a) (Fig. 6B). These metabolic pathways that up-regulated genes in the PPI analysis were also consistent with KEGG enrichment of DEGs and Venn diagram analysis. The detailed characterization of related genes was shown in Table S1.

Pathway analysis
The pathways found in the dormant group were elucidated in detail in Fig. 7. Moreover, the functional characterization and gene expression (basing on TPM value) were shown in Table 2

qRT-PCR analysis
Eighteen genes relating to the pathways in Fig. 7 were used to validate using qRT-PCR (Fig. 8). The results revealed that with the exception of FASNL, which was slightly down-regulated, all tested genes related to phagocytosis, fatty acid metabolism, and carbohydrate hydrolysis were increased in the dormant group. The relative expression of genes demonstrated that 17/18 (94%) genes showed a similar trend compared to those in the RNA-seq. These results indicated that the data obtained from RNA-seq is highly reliable. Basing on the RNA-seq, AMYB was detected in the dormant group but not in the revival group. The TPM value ratio (dormant group to revival group) was infinite, which is hard to exhibit in graphics. Thus, the biggest value of ten TPM value was used to replace infinity for showing high expression of AMYB in the dormant group (Fig. 8).

Discussion
The increase in summer temperature facilitates the disease outbreaks caused by microorganisms, which threaten the survival of marine organisms [12,13]. The higher temperature in lower latitudes usually causes worse effects on the animals [15]. To protect against pathogen infections, sea cucumbers had evolved an aestivation behavior in response to high temperature [4,5,20]. As previously comparable reports, most studies of A. japonicas were carried out in northern China [4,5,17,28]. However, the changes of metabolisms and immune responses of dormant A. japonicas adapting to climate warming were scarcely reported. In this study, the transcriptome analysis was performed to investigate the gene expression of A. japonicas in warmer southern China (Fujian province), which is important in exploring the adaption mechanism of A. japonicas in response to the higher temperature.
Totally, 2,368 DEGs (down-regulation: 927 and up-regulation: 1,441) were found in the dormant group compared with the revival group. The enrichment analysis of down-regulated genes revealed that 41 GO terms and 20 KEGG pathways were identified to be associated with the DNA, RNA, protein, and cellular metabolisms in the dormant group. These results implied that the physiological functions were decreased in the dormant group, and furthermore leading to the inhibition of growth and development of A. japonicas during the aestivation period. These results are similar to the previous studies that organisms reduced metabolisms of DNA, RNA, and protein during the dormant period, which further affected the growth and development [4,29,30]. In contrast, the up-regulated genes related to protein degradation, phagocytosis, carbohydrate metabolism, and fatty acid metabolism were identified through the GO and KEGG analyses. Basing on the KEGG database, long-chain fatty acids and saccharides are able to produce acetyl-CoA, which promotes the synthesis of ATP through the TCA cycle. The previous studies on the sea cucumbers cultured in northern China also acquired similar results that carbohydrate catabolism and fatty acid degradation for energy supplying were increased during the aestivation [4,20,30,31]. These results indicated that the dormant animals need to maintain energy supplying via the increasing consumption of stored carbohydrate and fatty acids as the food shortage [11]. Moreover, phagocytosis-related pathways (i.e. endocytosis, lysosome, and phagosome) were up-regulated in dormant sea cucumbers, which may play an important role in the innate immunity as a primary defense in the invertebrates [32][33][34][35][36][37]. Among these phagocytosisrelated genes, C-type lectin genes (macrophage mannose receptor 1-MRC1) was significantly upregulated (3.3-fold) in the dormant group. Previously, the up-regulation of the lysosome and C-type lectin were observed in the dormant sea cucumbers [4]. In the Venn diagram analysis, 2,255 specially expressed genes were obtained. The KEGG enrichment of these genes revealed that the phagocytosis (i.e. endocytosis, lysosome, and phagosome), as well as fatty acid and carbohydrate metabolism were altered in the dormant group, compared to that in the revival group. Therefore, the hypometabolism mainly occurred, but the fatty acid metabolism, carbohydrate metabolism, and phagocytosis activity were increased during the aestivation stage in sea cucumbers.
In the subsequent protein-protein interaction (PPI) analysis, in the dormant sea cucumbers, the core down-regulated genes (methyltransferase NSUN6, ribosomal protein, translation initiation factor, ribosome production factor 1, RNA cytidine acetyltransferase, U3 small nucleolar ribonucleoprotein protein, pre-rRNA-processing protein, eukaryotic translation initiation factor 3 subunit B, GMP synthase, and eukaryotic translation initiation factor 6) were assigned to the transcriptional and translational regulations [38]. Previously, the decreased levels of transcriptional and translational gene expression were also observed during the dormant stage of northern sea cucumber and other animals [4,30,39]. Once the animals enter the aestivation or hibernation states, slow movements decrease their opportunity to acquire adequate food. The food shortages cause energy deficiency for the body; the animals reduce the levels of metabolisms and energy consumption, leading to the inhibition of most of the gene expression and protein translation activities [40]. Furthermore, these down-regulated genes regulate the motility, skeletal development, neural activity, cell proliferation and development of A. japonicas according to Uniprot database analysis [41]. Our findings are consistent with those of other studies and suggested that the sea cucumber reduced the expression of transcriptional and translational genes relating to the hypometabolism (e.g., decreased motility, stopped skeletal development, declined neural activity, and reduced cell proliferation and development) [4,39,40]. This mechanism is accordant to reducing energy consumption in response to the absence of energy intake during the aestivation of sea cucumber at high temperature [2,39].
Basing on STRING database, the HK1 interact with the ubiquitin C variant (UBCV), then improved the phagosome activity via increasing the expression of V-type H+-transporting ATPase subunit A (VATPA), which was ATP-dependent proton pumps responsible for acidification of phagosome [46].  [47]. MRC1 is a typical pattern recognition receptor, which ligates with the mannose-rich glycoconjugates of pathogens [48,49]. The up-regulation of MRC1 indicated that the immune cells improved the recognition of pathogens in dormant A. japonicas, promoting the phagocytosis latterly. The increased expression of ARPC5A, ARPC5B, and PSD promoted the phagocytosis via forming vesicles for trafficking pathogen molecular to the phagosome, a crucial organelle in scavenging pathogens [16,[50][51][52]. Otherwise, the up-regulation of lysosomal GNS and GNSL accelerated the hydrolysis of heparan sulfate, which is important in responses to viral infection.
The heparan sulfate was also found to be a crucial material in triggering of autoimmune in the immune cells. These alterations effectively limited pathogen infections [53][54][55][56][57]. Moreover, the previous studies carried out in sea cucumbers also acquired a similar result that carbohydrate metabolism, lysozyme, and phagocytosis were up-regulated during the aestivation stage [4]. Thus, in the dormant A. japonicas, the high levels of carbohydrate metabolism promoted the phagocytosis activity, which plays a crucial role in scavenging pathogens and improving survival at seasonal high temperatures.
Furthermore, the TPM value was used to characterize the pathway alterations [58]. In the dormant A.
japonicas, high expression of PAMY, AMYB, AMY, and MGAM promoted starch hydrolysis to produce glucose, which is transferred into the TAC cycle via glycolysis pathway [42][43][44][45]. An up-regulation of ACADVL and ACADL in the dormant A. japonicas degraded long-chain fatty acids to acetyl-CoA, which was an essential material involved in the TCA cycle that is used for ATP synthesis [59,60]. A similar finding was observed in a previous study that starch and fatty acids hydrolase were highly expressed during the deep-aestivation of sea cucumbers [4]. These results suggested that the degradation of fatty acids and starch was the main route for providing energy in the dormant sea cucumbers. In addition, the HK1 interacts with VATPA, which then accelerates with H + and they together enter the phagosome, attending to improve phagocytic activities [46]. Thus, the high expression of HK1 indicates its roles in improving the immune system of sea cucumbers during the aestivation period.
Moreover, the phagocytosis-related genes (ARPC5A, ARPC5B, PSD, GNS, GNSL, and MRC1) were upregulated in the dormant group. This result is consistent with the previous results that an upregulation of carbohydrate hydrolase not only catalyzes starch degradation but also promotes phagocytosis activity. Optimal carbohydrate diets have significantly stimulatory effects on phagocytosis [61,62]. Therefore, our results indicated that the dormant A. japonicas consumed carbohydrate and fatty acids for energy supply. High levels of carbohydrate metabolism improved the expression of phagocytosis-related genes beneficial for resistance to pathogens during the aestivation of A. japonicas. However, the above-mentioned results were only based on transcriptome analysis of the RNA level, it is needed to investigate in detail in the future to validate the findings in this current study.

Conclusions
In this study, the transcriptome analysis was used to explore potential molecular mechanisms in the dormant A. japonicas. Compared with the revival group (RG), a total of 2,368 DEGs were identified in the dormant group (DG). Basing on GO and KEGG analysis, the down-regulated genes in DG were mainly annotated to DNA, RNA, and protein metabolic process. PPI analysis further revealed that the down-regulation of DNA, RNA, and protein metabolic process reduced the motility, skeletal development, neural activity, cell proliferation, and development of A. japonicas. In contrast, the upregulated genes were functionally classified into fatty acid and carbohydrate metabolisms. In the dormant A. japonicas, a high level of carbohydrate hydrolysis promoted the phagocytosis, which was a primary innate immune response in the invertebrates. Thus, the improvement of phagocytosis facilitates the A. japonicas increased their resistant ability against the pathogen invasion during the aestivation stage. The findings of this study provided new insights into the potential mechanisms that A. japonicas are able to change their physiological response in order to respond to the high temperature and pathogen invasion at the molecular level. The sea cucumber samples were ground in a ceramic pestle. The total RNA was extracted using Trizol Kit (Invitrogen, USA) following the manufacturer's instructions. The quality and quantity of total RNA were also estimated using 1% gel electrophoresis, Nano Photometer spectrophotometer (IMPLEN, CA, USA), and Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). The obtained total RNA was used to synthesize the cDNA libraries using the PrimeScript RT reagent Kit with gDNA Eraser (Takara, Dalian, China) following the manufacturer's instructions. The cDNA was used as the templates for qRT-PCR and RNA sequencing (RNA-Seq) transcriptome.

Illumina RNA sequencing
At least 1 μg of total mRNA was used to construct cDNA libraries. The purified mRNA from total RNA was cut into fragments by using the NEBNext ® Ultra TM RNA Library Prep Kit (Illumina, USA) for constructing the cDNA library. Put simply, oligo (dT) magnetic beads enriched mRNA with polyA, firstly. NEB fragmentation buffer was used to randomly fragment the mRNA into small fragments, which was then used as templates for synthesizing the first-strand cDNA using random hexamer primer and M-MuLV Reverse Transcriptase. Subsequently, DNA Polymerase I and RNase H were used to produce the second-strand cDNA with dNTPs. Double-strand cDNA was purified accepted end repair, A-base addition, and ligation of the Illumina-indexed adaptors. Finally, 150~200 bp cDNA was used for the PCR amplification and the PCR products were purified using AMPure XP beads (Beckman Coulter, Beverly, USA), which was used to generate cDNA libraries. The cDNA library quality was further assessed using the Agilent Bioanalyzer 2100 system. The concentration of cDNA over 2 nM was used to sequence by an Illumina Hiseq TM 2500 platform.

Sequence assembly and gene functional annotation
Quality control was performed removing adaptor sequences, reads containing poly-N and low-quality reads, then calculated the Q20 and Q30 of clean reads. HISAT2 v2.0.5 software was used to further aligned and obtained the clean reads to reference the genome of A. japonicus that the accession number is PRJNA354676 [27]. The novel transcript assembly was performed using StringTie (version GO enrichment of DEGs was performed with clusterProfiler R package and corrected the bias of gene length. For better understanding the functions of the genes in a biological system, KEGG enrichment was also carried out with clusterProfiler R package, which well-explained the genes at molecular levels. Moreover, the specific expression of genes in dormant and revival groups was analyzed and exhibited in a Venn diagram. Specific genes were also performed with GO and KEGG enrichment.

Protein-protein interaction (PPI) analysis
Protein-protein interaction analysis of DEGs was carried out using the STRING database (http://stringdb.org/), which investigated the direct and indirect associations among different proteins. The networks showing up-and down-regulated genes were drawn using Cytoscape software based on the results of STRING analysis.

Quantitative real-time PCR (qRT-PCR) analysis of metabolism-related genes
To validate the results of gene expression obtained from the transcriptomic sequencing, quantitative real-time PCR (qRT-PCR) was performed to detect the gene expression using the SYBR ® Premix Ex Taq TM II Kit (Takara, Dalian, China) in LightCycler ® 480 (Roche, USA). Total RNA from A. japonicus was extracted using Total RNA kits Ⅱ (Omega, USA) following the manufacturer's instructions. cDNA was synthesized using PrimeScript RT reagent Kit with gDNA Eraser (Takara, Dalian, China) following the manufacturer's instructions. The obtained cDNA was used as templates for qRT-PCR. In this study, 18 genes relating to carbohydrate metabolism, fatty acid metabolism, TCA circle, and phagocytosis pathways were used in the qRT-PCR. A 20-μL reaction mixture contained 10 μL of SYBR ® Premix Ex Taq TM II, 2 μL of cDNA, 0.8 μL (10 mM) of each (forward and reverse) primer (Table 3) Table 3. Sequences of primers used in this study.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.
Additional files 1.doc