- Research article
- Open Access
Insights into high-pressure acclimation: comparative transcriptome analysis of sea cucumber Apostichopus japonicus at different hydrostatic pressure exposures
BMC Genomics volume 21, Article number: 68 (2020)
Global climate change is predicted to force the bathymetric migrations of shallow-water marine invertebrates. Hydrostatic pressure is proposed to be one of the major environmental factors limiting the vertical distribution of extant marine invertebrates. However, the high-pressure acclimation mechanisms are not yet fully understood.
In this study, the shallow-water sea cucumber Apostichopus japonicus was incubated at 15 and 25 MPa at 15 °C for 24 h, and subjected to comparative transcriptome analysis. Nine samples were sequenced and assembled into 553,507 unigenes with a N50 length of 1204 bp. Three groups of differentially expressed genes (DEGs) were identified according to their gene expression patterns, including 38 linearly related DEGs whose expression patterns were linearly correlated with hydrostatic pressure, 244 pressure-sensitive DEGs which were up-regulated at both 15 and 25 MPa, and 257 high-pressure-induced DEGs which were up-regulated at 25 MPa but not up-regulated at 15 MPa.
Our results indicated that the genes and biological processes involving high-pressure acclimation are similar to those related to deep-sea adaptation. In addition to representative biological processes involving deep-sea adaptation (such as antioxidation, immune response, genetic information processing, and DNA repair), two biological processes, namely, ubiquitination and endocytosis, which can collaborate with each other and regulate the elimination of misfolded proteins, also responded to high-pressure exposure in our study. The up-regulation of these two processes suggested that high hydrostatic pressure would lead to the increase of misfolded protein synthesis, and this may result in the death of shallow-water sea cucumber under high-pressure exposure.
The ocean is warming because of global climate change, forcing the bathymetric migrations of shallow-water marine invertebrates [1, 2]. As such, the ability of a shallow-water invertebrate to acclimatize to deep-sea environments during its lifetime is vital. The bathymetric migrations of marine fauna are predicted to be constrained by the combined effects of temperature, hydrostatic pressure, and oxygen concentration . Among them, hydrostatic pressure is thought to be the major environmental factor that limits the vertical distribution of extant marine fauna [3, 4]. Many studies have examined the tolerance of shallow-water invertebrates to high hydrostatic pressure and low temperature (reviewed by Brown & Thatje 2014) , indicating that many extant marine benthic invertebrates can tolerate hydrostatic pressure outside their known natural distributions, and a low temperature can impede high-pressure acclimation. Although a few studies focused on DEGs responding to high-pressure exposure [6,7,8], transcriptome analysis was seldom applied to relevant studies, and the molecular mechanisms of shallow-water invertebrates to acclimatize to high-pressure environment is not yet fully understood. This question is important in the present context of climate change and ocean warming.
Most extant deep-sea fauna are accepted to have originated from shallow waters as a consequence of a series of extinction events during the Phanerozoic [9, 10]. The colonization of the deep sea occurs throughout selection and during the slow genetic drift of species that gradually adapt to life in this area , whereas the high-pressure acclimation of shallow-water fauna involve physiological plasticity in response to a simulated immersion in the high-pressure environments. However, both evolutionary adaptation and phenotypic acclimation are essential for adaptation to high pressure . Transcriptome analysis has been applied widely to study the adaptation mechanisms of deep-sea fauna based on the comparisons of congeneric species that have different vertical distribution profiles. Common adaptation patterns have been observed in different taxa of deep-sea living fauna . Many biological processes, including alanine biosynthesis , antioxidation [14, 15], energy metabolism [13, 16], immunity [16, 17], fatty acid metabolism , and genetic information processing , are related to deep-sea adaptation.
Somero (1992) has reviewed the effects of hydrostatic pressure on shallow-water organisms . One of the most sensitive molecular assemblages of hydrostatic pressure is lipid bilayer [11, 19,20,21,22]. High pressure leads to a reduction of membrane fluidity, impeding physiological membrane functions, such as transmission [20, 23], transmembrane transportation, and cell movement [24, 25]. The effects of high hydrostatic pressure and low temperature are similar [26, 27]. Parallel effects can be detected on the basis of membrane composition with an increase in hydrostatic pressure of 100 MPa and a reduction in temperature of 13–21 °C . A high hydrostatic pressure causes the depolymerization of protein structures, whereas a low temperature negatively affects protein activity, and both factors induce an increase in protein chaperoning, thereby decreasing the stabilization of secondary RNA and DNA structures [28, 29]. High pressure can also strengthen hydrogen bonds. Consequently, processes that include DNA replication, transcription, and translation are impeded [30, 31].
The sea cucumber Apostichopus japonicus (phylum: Echinodermata) is a temperate species mainly distributed along the coastal area of eastern Asia . It is also a popular food in China because of its high nutritional and medicinal value. Sea cucumbers of Echinodermata are not only ubiquitous in coastal areas but also widespread at abyssal depth [33, 34]. Since deep-sea species do not obtain new genes, but utilize gene sets homologous to their coastal relatives to adapt to deep-sea environments , we predicted that A. japonicus has the potential to acclimatize to high-pressure environment, and used this species in high-pressure incubations. A pressure vessel was used to perform high-pressure exposure on experimental samples, provide a stable and controllable experimental context, and examine pressure acclimation accurately .
Hydrostatic pressure tolerance of A. japonicus and experimental design
To examine the pressure tolerance of A. japonicus, we incubated 10 individuals at different high-pressure conditions and measured their mortality rate before formal experiments for transcriptome analysis. There were 30% individuals died after 24-h incubation at 35 MPa, but no individual died at 25 MPa. Additionally, eversion was not observed at 25 MPa, which is usually happened when sea cucumbers are stressed. Consequently, 3 pressure conditions were set: 0.1 MPa (atmospheric pressure), 15 MPa (pressure at the depth of 1500 m), and 25 MPa (pressure at the depth of 2500 m). A total of 9 individuals (3 individuals from each experimental group) were high-pressure incubated for transcriptome analysis. The RNA of body wall tissue from each individual was sequenced, and paired reads of these 9 samples were assembled into one tanscriptome.
Sequencing, assembly and annotation
Three experimental groups (P0.1, experimental group incubated at 0.1 MPa; P15, experimental group incubated at 15 MPa; and P25, experimental group incubated at 25 MPa) were used for comparative transcriptome analysis. Each experimental group had three replications. Sequencing qualities are listed in Additional file 2: Table S1. Paired reads from the nine samples were assembled into 553,507 unigenes with a total length of 481,946,001 bp and an N50 length of 1204 bp. BUSCO completeness of the transcriptome were 91.5% (single-copy: 28.4%, duplicated: 63.1%, fragmented: 7.2%, missing: 1.3%). There were 14, 23, and 7% unigenes annotated in the databases of Swiss-Prot, Protein family (Pfam), and Kyoto Encyclopedia of Genes and Genomes (KEGG), respectively.
DEGs involved in high-pressure acclimation
Three combinations, namely, P15 vs. P0.1, P25 vs. P0.1, and P25 vs. P15, were subjected to differential expression analysis by using the DESeq2 R package (v1.22.2) . In this study, up-regulated genes were considered as activated genes in response to high-pressure exposure because only essential processes can be maintained, whereas nonessential processes are reduced outside the optimal range [37,38,39,40]. A total of 598 genes, 1375 genes, and 542 genes were significantly up-regulated in the combinations of P15 vs. P0.1, P25 vs. P0.1, and P25 vs. P15, respectively (Fig. 1). In addition, quantitative real-time reverse transcription-PCR (qPCR) analysis was used to validate the reliability of the RNA-seq results. A total of 14 DEGs were employed for qPCR analysis, and the Pearson correlation coefficients between RNA-seq and qPCR results ranged from 0.81 to 0.99.
Three groups of DEGs comprising 38 linearly related DEGs (LRGs), 244 pressure-sensitive DEGs (PSGs), and 257 high-pressure-induced DEGs (HPGs) (Fig. 1) were identified according to their gene expression patterns. LRGs were up-regulated among the three combinations. PSGs were up-regulated only in P15 vs. P01 and P25 vs. P01. HPGs were up-regulated only in P25 vs. P01 and P25 vs. P15. The expression pattern of LRGs was linearly correlated with hydrostatic pressure (R2 > 0.99, Fig. 2a). The PSGs were significantly up-regulated at 15 MPa and remained at a similar high level at 25 MPa (Fig. 2b). The HPGs were significantly up-regulated at 25 MPa but were not significantly up-regulated at 15 MPa (Fig. 2c).
Swiss-Prot annotation of LRGs, PSGs, and HPGs
The expression patterns of 38 DEGs are linearly related to hydrostatic pressure, and 14 of them are annotated in the Swiss-Prot database (Additional file 3: Table S2). Their functions are mainly involved in homeostasis maintenance (7 genes) and lysosomal activities (3 genes) (Fig. 3a). Four of the seven homeostasis maintenance genes, namely, E3 ubiquitin-protein ligase NEURL1 (NEURL1), E3 ubiquitin-protein ligase RNF14 (RNF14), E3 ubiquitin-protein ligase dbl4 (dbl4), and E3 ubiquitin-protein ligase rbrA (rbrA), are involved in ubiquitination. The three other genes involved in homeostasis maintenance are DnaJ homolog subfamily B member 4 (DNAJB4), cytochrome P450 2 U1 (Cyp2u1), and interleukin-1 receptor-associated kinase 4 (IRAK4). DnaJ, also known as heat shock protein 40, is a molecular chaperone protein regulating the ATPase activity of heat shock protein 70 (HSP70) . Cytochrome P450 proteins (CYPs) are known for their antioxidative functions . The IRAK4 protein is a key regulatory kinase of innate immunity . Three genes, namely, syntaxin-12 (STX12) that regulates protein transport between late endosomes and the trans-Golgi network, TBC1 domain family member 15 (TBC1D15) that promotes fusion events between late endosomes and lysosomes , and zinc finger FYVE domain-containing protein 1 (ZFYVE1) that has been related to vacuolar protein sorting and endosome function, are implicated in lysosomal activities. Two genes, namely, CCAAT/enhancer-binding protein beta (CEBPB) that regulates the glucose homeostasis  and glycogen debranching enzyme (AGL) that facilitates the breakdown of glycogen and serves as glucose storage, participate in energy metabolism . Two genes, namely, ATP-binding cassette sub-family A member 3 (Abca3)  and putative phospholipase B-like 2 (PLBD2), function in lipid metabolism.
A total of 244 genes are PSGs, and 70 of them were annotated in Swiss-Prot database (Additional file 4: Table S3). These 70 genes were grouped into seven different biological processes, namely, homeostasis maintenance (15 genes), signal transduction (15 genes), genetic information processing (12 genes), lysosomal activities (5 genes), membrane related functions (5 genes), lipid metabolism (2 genes), and others (16 genes) (Fig. 3b). Of the 15 genes grouped in homeostasis maintenance, 6 are involved in stress responses, including ankyrin repeat and LEM domain-containing protein 1 (ANKLE1) involving DNA damage response and DNA repair, CREB3 regulatory factor (CREBRF) involving unfolded protein response, and MAP kinase-activated protein kinase 2 (MAPKAPK2) involving cell migration, cell cycle control, DNA damage response, and transcriptional regulation; 6 are implicated in immune response, including histidine triad nucleotide-binding protein 2 (HINT2) involving apoptosis; and 3 participate in ubiquitination. Of the 12 genes grouped in genetic information processing, 7 function in transcription.
A total of 257 genes are HPGs, and 123 of them were annotated in Swiss-Prot database (Additional file 5: Table S4). These genes were grouped into six different biological processes, namely, homeostasis maintenance (23 genes), genetic information processing (22 genes), signal transduction (12 genes), lysosomal activities (11 genes), membrane related functions (7 genes), lipid metabolism (2 genes) and others (46 genes) (Fig. 3c). Of the 23 genes grouped in homeostasis maintenance, 13 are involved in ubiquitination, including ubiquitin-conjugating enzyme E2 R2 (UBE2R2), E3 ubiquitin-protein ligase NEDD4, PELI1, RBBP6, and RNF31; 8 are implicated in stress response, including cytochrome P450 Cyp3a11 and CYP3A6, heat shock 70 protein IV (HSP70IV), AN1-type zinc finger protein 2B (Zfand2b), ankyrin repeat and zinc finger domain-containing protein ANKZF1 and Ankzf1; and 2 participate in immune response. Zfand2b is a recently identified heat shock protein . ANKZF1 and Ankzf1 play a role in the cellular response to hydrogen peroxide. Of the 22 genes grouped in genetic information processing, 12 and 7 are involved in transcription and translation, respectively.
KEGG and Pfam enrichment analysis
The KEGG enrichment analysis of LRGs, PSGs, and HPGs were separately implemented by using the KOBAS software . No significantly enriched KEGG pathway existed in any groups of genes except the pathway of endocytosis in HPGs. A total of 14 genes were annotated in this KEGG pathway. Additionally, KEGG enrichment analysis was applied to 539 genes of the assemblage of LRGs, PSGs, and HPGs. Endocytosis was also the most significantly enriched KEGG pathway (Additional file 1: Figure S1). A total of 17 genes were annotated in this KEGG pathway, and most of them were involved in clathrin-dependent endocytosis (Fig. 4 and Additional file 6: Table S5).
The Pfam enrichment analysis of LRGs, PSGs, and HPGs were separately implemented by using fisher.test function of R software  in LRGs, PSGs, and HPGs. A total of 13, 13 and 20 gene families were significantly enriched in LRGs, PSGs, and HPGs, respectively (Fig. 5).
A total of 13 gene families were significantly enriched in LRGs (Fig. 5a and Additional file 7: Table S6). Three of them, namely, bZIP Maf transcription factor (bZIP Maf), bZIP transcription factor (bZIP 1), and basic region leucine zipper (bZIP 2), are involved in transcription. Two gene families, namely, ring finger domain (zf-RING 2) and zinc-RING finger domain (zf-RING 5), are implicated in the ubiquitination pathway. Two gene families, namely, cytokine-induced anti-apoptosis inhibitor 1apoptosis inhibitor 1 (CIAPIN1) and winged helix-turn-helix transcription repressor (HrcA DNA-bdg), participate in oxidative stress and heat-shock stress response, respectively.
A total of 13 gene families were significantly enriched in PSGs (Fig. 5b and Additional file 8: Table S7). Five of them were involved in transcription (bZIP Maf, bZIP 1, bZIP 2, vestigial family [Vg Tdu], and sterile alpha motif domain [SAM PNT]). Two gene families, namely, ligated ion channel L-glutamate- and glycine-binding site (Lig chan-Glu bd) and ligand-gated ion channel (Lig chan), are implicated in transmembrane ion transportation. The Mus7/MMS22 family (Mus7) participates in DNA damage repair.
A total of 20 gene families were significantly enriched in HPGs (Fig. 5c and Additional file 9: Table S8). Six of them are involved in genetic information related functions. RNA polymerase Rpb1 domain 5 (RNA pol Rpb1 5) catalyzes DNA-dependent RNA polymerization. 50S ribosome-binding GTPase (MMR HSR1) is required for the complete activity of a protein interacting with the 50S ribosome. Rit1 DUSP-like domain (Init tRNA PT) participates in the initiation and elongation of translation. PRP1 splicing factor (PRP1 N) is implicated in mRNA splicing. The regulator of RNA polymerase sigma subunit (Rsd AlgQ) and bZIP Maf function in transcription. Four gene families participate in endocytosis, including ADP ribosylation factor (Arf), Snf7, VHS protein domain (VHS) and coatomer WD associated region (Coatomer WDAD).
The optimum temperature of A. japonicus ranges from 10 °C to 17 °C , and A. japonicus hibernates in winter. The characteristics of A. japonicus in hibernation states were quite different from higher animals, but more closely resembled a semi-dormant state. The shift from normal to hibernation was a chronic process, indicated by the gradual depression of metabolic rate of about 71.7% . The water temperature nearly stays constant at 2 °C below the depth of 2000 m . As such, this species is not likely to survive in the deep-sea environments because of the low temperature. However, the scientific question of this study is how shallow-water invertebrates acclimatize to high-pressure environment, and we suggested that the acclimation mechanisms identified in the species A. japonicus are similar to other sea cucumber species. Thus we did not simulate the same environments as the deep sea in this study, but examined the molecular responses of A. japonicus to high-pressure exposures at 15 °C to prevent variation caused by hibernation, and set hydrostatic pressure as the only variation.
Homeostatic effort is required to maintain internal conditions within their physiological tolerance boundaries outside optimum. Consequently, only essential processes can be maintained, whereas nonessential processes are reduced [37,38,39,40]. Survival under such condition is time limited. Although A. japonicus can survive at 25 MPa for 24 h, whether it can survive at such pressure condition for longer time is currently unclear. New, et al. (2014) found that the acclimation period of shallow-water shrimp Palaemonetes varians to high-pressure condition was 1 week . Thus a long-term high-pressure incubation (1–4 weeks) of A. japonicus can provide information to answer this question. However, since the pressure system in used was isolated, we only incubated A. japonicas for 24 h to avoid the deterioration of water qualities. The 24-h high-pressure incubation in this study is a first approach. Long-term and time-series high-pressure exposures are the future goal to fully address the molecular mechanisms of A. japonicus to acclimatize to high-pressure exposure.
Although LRGs, PSGs, and HPGs have different expression patterns, their up-regulated biological processes are similar. The biological process homeostasis maintenance has the highest proportion in the three groups of DEGs. Additionally, representative biological processes, such as antioxidation, stress response, and immune response, are relevant in many other studies about deep-sea adaptation; similarly, some representative genes, such as HSPs, CYPs, and zinc finger protein, are also involved in deep-sea adaptation [13, 15,16,17,18, 54]. It has been proved that the ability of antioxidation can be beneficial to high pressure adaptation: the bacterium Shewanella piezotolerans mutant OE100, which enhanced antioxidant defense capacity by experimental evolution under H2O2 stress, has better tolerance to high pressure . HSPs were also reported to play important role in the maintenance of protein structure which is highly influenced by high pressure . However, DEGs involved in ubiquitination observed in this study were not identified in most relevant studies about deep-sea adaptation. Three enzymes are involved in ubiquitination, including E1 ubiquitin-activating, E2 ubiquitin-conjugating, and E3 ubiquitin-ligating enzymes. Most DEGs participating in ubiquitination in our results were annotated as E3 ubiquitin ligase of RING domin type. E3 ligases can recognize target substrates and facilitate the transfer of ubiquitin from an E2 ubiquitin-conjugating enzyme to its substrate. The number of ubiquitin transferred to substrate can be multiple. Therefore, these modifications can have diverse effects on the substrate, including proteasome-dependent proteolysis, modulation of protein function, structure, assembly, and localization (reviewed by Deshaies & Joazeiro, 2009 ).
Endocytosis is the most significantly enriched KEGG pathway in this study. Endocytosis in eukaryotic cells is characterized by the continuous and regulated formation of prolific numbers of membrane vesicles at the plasma membrane . In general, these vesicle types result in the delivery of their contents to lysosomes for degradation. Studies on deep-sea mussels have reported that endocytosis is essential for the acquisition of symbionts [16, 18]. As such, this process has been expanded to the mussel genome. Therefore, we assumed that high pressure could accelerate the development of a deep-sea symbiotic system. Additionally, one of the effects of protein ubiquitination is proteasome-dependent proteolysis, which can activate the following endocytosis. Ubiquitination and endocytosis can collaborate with each other and regulate the elimination of misfolded proteins which resulted from high hydrostatic pressure. The significant up-regulation of these two processes suggested that high hydrostatic pressure would lead to the increase of misfolded protein synthesis, and this may be one of the main reasons resulting in the death of shallow-water sea cucumber under high-pressure exposure.
Gene families involving genetic information related functions, especially transcription, were highly enriched in the three groups of DEGs. Since high pressure can strengthen hydrogen bonds and impedes genetic information related processes [30, 31], the up-regulation of these genes can remit the effects of high pressure. Additionally, genes related to this process were also significantly positive selected in deep-sea amphipod Hirondellea gigas . This study suggested that low temperature in deep-sea environments results in the positive selection of these gene families. However, the incubation temperature in our experiments was optimal. We assumed that high pressure also plays an important role in the positive selection of gene families related to genetic information processing. High pressure can cause DNA chain breakage and damage . Thus, high frequencies of DNA repair are needed. The gene family Mus7 and the genes ANKLE1 and MAPKAPK2 that participate in the repair of replication-associated DNA damage were also found significantly up-regulated at high-pressure condition in our study.
Shallow-water sea cucumber A. japonicus could survive 100% under 25 MPa at 15 °C for at least 24 h. However, whether this shallow-water species could survive at this high-pressure condition for more than 24 h or permanently remained unclear. The 24-h high-pressure incubation in this study is a first approach. Long-term and time-series high-pressure exposures are the future goal to fully address high-pressure acclimation mechanisms.
Although LRGs, PSGs, and HPGs had different expression patterns, their up-regulated biological processes are similar. Our results also indicated that genes and biological processes involving high-pressure acclimation were similar to those related to deep-sea adaptation. Representative biological processes, including antioxidation, stress response, genetic information processing, and DNA repair, were found significantly induced by high-pressure exposure in our results. In addition, some representative genes, such as HSPs, CYPs, and zinc finger protein, were also detected at high pressure. Moreover, two important biological processes, namely, ubiquitination and endocytosis, which can collaborate with each other and regulate the elimination of misfolded proteins, were found significantly induced by high-pressure exposure. The significant up-regulation of these two processes suggested that high hydrostatic pressure would lead to the increase of misfolded protein synthesis, and this may result in the death of shallow-water sea cucumber under high-pressure exposure.
Collection, maintenance, and high-pressure incubation of A. japonicus
Juvenile individuals of A. japonicus (length: 5 ± 1 cm, weight: 3.9 ± 0.5 g) were caught from an aquaculture farm in Shandong, China on December 2017. The length was measured in sea water when individuals were under normal activities, and wipe-dry living individuals were used for weight measurement. They were maintained at a closed recirculating aquacultural system (seawater was partially refreshed twice a week). The sea cucumber were reared in aerated filtered seawater (salinity: 34.7–35.3, 1 μm filtered, natural light cycle), and were fed with algae powder twice a week; unconsumed food was removed after 24 h via refreshing seawater. Experimental individuals were starved for 3 d prior to pressurization. The experimental individuals were maintained in the aquacultural system at 15 °C for 2 weeks to acclimatize them to laboratory environments and allow them to recover from collection and handling stress.
The hydrostatic pressure system was set to 15 °C by using circulating water bath at least 6 h prior to each experiment. Three individuals were placed inside the stainless pressure chamber (volume: ~ 20 l, internal diameter 20 cm, internal depth 65 cm)  which is full of filtered seawater at a constant temperature (15 ± 0.8 °C), and maintained for 1 h to allow acclimation and recovery from handling stress. Then, the pressure vessel was pressurized at 1 MPa per minute to experimental pressure by using hydraulic pump. The pressure chamber was sealed and isolated during this time period. After hydrostatic pressure exposure for 24 h, the pressure was released instantaneously, and these 3 high-pressure incubated individuals were removed from the pressure chamber and snap frozen in liquid nitrogen. The maximum time elapsed between departure from experimental pressure and flash freezing is 3 min. The flash-frozen individuals were then stored at − 80 °C for further RNA extraction.
YSI Professional Plus (YSI Inc., USA) and HACH DR 1900 (HACH Company, USA) were used to ascertain the stabilization of seawater quality before and after each experiment. No significant difference was found between the experimental context before pressurization (dissolved oxygen: 5.76 ± 0.3 mg.L− 1, salinity: 35.0 ± 0.3, pH: 8.1 ± 0.1, NO2-N: 0.0055 ± 0.0005 mg.L− 1, NH3-N: 0.015 ± 0.005 mg.L− 1, NO3-N: 0.015 ± 0.005 mg.L− 1) and after pressurization (dissolved oxygen: 5.50 ± 0.4 mg.L− 1, salinity: 35.0 ± 0.3, pH: 8.1 ± 0.1, NO2-N: 0.007 ± 0.0005 mg.L− 1, NH3-N: 0.03 ± 0.005 mg.L− 1, NO3-N: 0.015 ± 0.005 mg.L− 1).
RNA extraction, sequencing, assembly and annotation
Details of this part have been described in our previous study . In brief, approximately 50 mg body wall tissue from each individual was used for RNA extraction because body wall is easy to be separated from the frozen main body and this tissue is also commonly used for sea cucumber RNA extraction . The tissue was dissected before melted, and immediately transferred into 1 ml of QIAzol (from RNeasy Plus Universal Kit) and total RNA was extracted by RNeasy Plus Universal Kit (QIAGEN, UK) according to the manufacturer’s protocol. A total of 1.5 μg RNA per sample was used for the RNA sample preparations. Then, library preparations were sequenced on Illumina Hiseq X platform and 150 bp paired-end reads were generated. Transcriptome assembly was accomplished based on the clean data by using Trinity . At last, transcripts were hierarchical clustering by Corset . BUSCO (v.3.0.2) was used to estimate the completeness of the final assembly with default settings and using the metazoan orthologs. All unigenes were annotated in three databases, namely, Swiss-Prot, Pfam and KEGG, by using diamond (v0.8.22), HMMER (v3.0) and KEGG Automatic Annotation Server (KAAS), respectively.
Differential expression analysis
FPKM and read count were calculated by using RSEM , and FPKM were normalized by using TMM method. Then, differential expression analysis was implemented by using DESeq2 R package (v1.22.2)  to identify DEGs involved in high-pressure acclimation. Only genes with an adjusted p-value < 0.01 and |log2 (fold change)| > 1 were regarded as DEGs. In this study, up-regulated genes were considered as activated genes in response to high-pressure exposure because only essential processes can be maintained, whereas nonessential processes are reduced outside the optimal range [37,38,39,40]. The relationships of DEGs among different combinations (P15 vs. P0.1, P25 vs. P0.1, and P25 vs. P15) were visualized using the VennDiagram R package (v1.6.20) . The normalized FPKM of each gene were scaled to the control treatment (P0.1) and converted to the relative fold change (RFC). Then, log2 (RFC) was used to evaluate differential expression patterns. The expression patterns of LRGs, PSGs, and HPGs were shown on line graphs by using Graphpad Prism 7 (Graphpad software, Inc., USA). Points represented the mean of log2 (RFC) of all genes, and error bars represent standard deviation.
The LRGs, PSGs, and HPGs annotated in Swiss-Prot were grouped in several biological processes. Their normalized FPKM, except DEGs grouped into “others”, were visualized by using the pheatmap R package (v1.0.12) . The distance between rows was calculated with vegan R package (v2.5–5)  by using the bray method. Pfam and KEGG enrichment analyses (based on Fisher’s exact test) were implemented to identify enriched biological processes in each experimental group by using fisher.test function of R software and KOBAS software , respectively. The results of enrichment analysis were visualized by using the ggplot2 R package (v3.2.1) .
RNA-seq validation by qPCR
A total of 14 DEGs were employed for qPCR by StepOnePlus Real-Time PCR System (Applied Biosystems, USA) to validate the RNA-seq results. The primer sequences (Additional file 10: Table S9) were designed by Primer Premier 5.0 software (Premier Biosoft International, Palo Alto, CA, USA). Details of this part have been described in our previous study . In brief, relative standard curve method was used for expression level analysis with cytb and β-actin as internal controls. Six 3-fold serial dilutions were performed on a cDNA template. The melting curve analysis was performed at the end of each PCR to confirm only one PCR product was amplified. At last, Pearson correlation coefficients between RNA-seq and qPCR results were calculated. Pearson correlation coefficients (ranging from 0.81 to 0.99) are supplied in Additional file 10: Table S9, which validated the reliability of the RNA-seq results.
Availability of data and materials
Relevant data supporting this manuscript are contained within the figures of this manuscript or provided in the supplementary material. The clean data of transcriptome were available from National Center for Biotechnology Information Sequence Read Archive database (SRA accession numbers: SRR8269388, SRR8269389, SRR8269390, SRR8269391, SRR8269384, SRR8269385, SRR8269386, SRR8269387, SRR8269383). The transcriptome assembly project has been deposited at DDBJ/ENA/GenBank under the accession GHDI00000000. The version described in this paper is the first version, GHDI01000000. Some R scripts used in this study have been uploaded to GitHub: https://github.com/CHEN-Jiawei-Jason/R_scripts_for_BMC_Genomics.git
Differentially expressed genes
Fragments per kilobase of transcript per million fragments mapped
Linearly related DEGs
Experimental group incubated at atmospheric pressure
Experimental group incubated at 15 MPa
Experimental group incubated at 25 MPa
Quantitative real-time reverse transcription-PCR.
Relative fold change
Pinsky ML, Worm B, Fogarty MJ, Sarmiento JL, Levin SA. Marine taxa track local climate velocities. Science (New York, NY). 2013;341(6151):1239–42.
Brown A, Thatje S. The effects of changing climate on faunal depth distributions determine winners and losers. Glob Chang Biol. 2015;21(1):173–80.
Young CM, Tyler PA, Fenaux L. Potential for deep sea invasion by Mediterranean shallow water echinoids:pressure and temperature as stage-specific dispersal barriers. Mar Ecol Prog Ser. 1997;154:197–209.
Thatje S, Hillenbrand CD, Larter R. On the origin of Antarctic marine benthic community structure. Trends Ecol Evol. 2005;20(10):534–40.
Brown A, Thatje S. Explaining bathymetric diversity patterns in marine benthic invertebrates and demersal fishes: physiological contributions to adaptation of life at depth. Biol Rev Camb Philos Soc. 2014;89(2):406–26.
Morris JP, Thatje S, Ravaux J, Shillito B, Fernando D, Hauton C. Acute combined pressure and temperature exposures on a shallow-water crustacean: novel insights into the stress response and high pressure neurological syndrome. Comp Biochem Physiol A Mol Integr Physiol. 2015;181:9–17.
Cottin D, Brown A, Oliphant A, Mestre NC, Ravaux J, Shillito B, Thatje S. Sustained hydrostatic pressure tolerance of the shallow water shrimp Palaemonetes varians at different temperatures: insights into the colonisation of the deep sea. Comp Biochem Physiol A Mol Integr Physiol. 2012;162(4):357–63.
Brown A, Thatje S, Morris JP, Oliphant A, Morgan EA, Hauton C, Jones DOB, Pond DW. Metabolic costs imposed by hydrostatic pressure constrain bathymetric range in the lithodid crab Lithodes maja. J Exp Biol. 2017;220(Pt 21):3916–26.
Harnik PG, Lotze HK, Anderson SC, Finkel ZV, Finnegan S, Lindberg DR, Liow LH, Lockwood R, Mcclain CR, Mcguire JL. Extinctions in ancient and modern seas. Trends Ecol Evol. 2012;27(11):608–17.
Raup DM, Sepkoski JJ. Mass extinctions in the marine fossil record. Science. 1982;215(4539):1501–3.
Delong EF, Yayanos AA. Adaptation of the membrane lipids of a deep-sea bacterium to changes in hydrostatic pressure. Science. 1985;228(4703):1101–3.
Somero GN. Biochemical ecology of deep-sea animals. Experientia. 1992;48(6):537–43.
Lan Y, Sun J, Tian R, Bartlett DH, Li R, Wong YH, Zhang W, Qiu JW, Xu T, He LS, et al. Molecular adaptation in the world's deepest-living animal: insights from transcriptome sequencing of the hadal amphipod Hirondellea gigas. Mol Ecol. 2017;26(14):3732–43.
Xie Z, Jian H, Jin Z, Xiao X. Enhancing the Adaptability of the Deep-Sea Bacterium Shewanella piezotolerans WP3 to High Pressure and Low Temperature by Experimental Evolution under H2O2 Stress. Appl Environ Microbiol. 2018;84(5):e02342–17.
Zhang J, Sun QL, Luan ZD, Lian C, Sun L. Comparative transcriptome analysis of Rimicaris sp. reveals novel molecular features associated with survival in deep-sea hydrothermal vent. Sci Rep. 2017;7(1):2000.
Sun J, Zhang Y, Xu T, Zhang Y, Mu H, Zhang Y, Lan Y, Fields CJ, Hui JHL, Zhang W, et al. Adaptation to deep-sea chemosynthetic environments as revealed by mussel genomes. Nat Ecol Evol. 2017;1(5):121.
Oliver TA, Garfield DA, Manier MK, Haygood R, Wray GA, Palumbi SR. Whole-genome positive selection and habitat-driven evolution in a shallow and a deep-sea urchin. Genome Biol Evol. 2010;2(10):800–14.
Zheng P, Wang M, Li C, Sun X, Wang X, Sun Y, Sun S. Insights into deep-sea adaptations and host-symbiont interactions: a comparative transcriptome study on Bathymodiolus mussels and their coastal relatives. Mol Ecol. 2017;26(19):5133–48.
Somero GN. Adaptations to high hydrostatic pressure. Annu Rev Physiol. 1992;54:557–77.
Wann KT, Macdonald AG. The effects of pressure on excitable cells. Comp Biochem Physiol Part A Physiol. 1980;66(1):1–12.
Macdonald AG. Hydrostatic pressure as an environmental factor in life processes. Comp Biochem Physiol Part A Physiol. 1997;116(4):291–7.
Winter R, Dzwolak W. Exploring the temperature-pressure configurational landscape of biomolecules: from lipid membranes to proteins. Philos Trans A Math Phys Eng Sci. 2005;363(1827):537–62 discussion 562-533.
Siebenaller JF, Garrett DJ. The effects of the deep-sea environment on transmembrane signaling. Comp Biochem Physiol Part B Biochem Mol Biol. 2002;131(4):675.
Oger PM, Jebbar M. The many ways of coping with pressure. Res Microbiol. 2010;161(10):799–809.
Winter R. Pressure effects on lipid membranes. Berlin Heidelberg: Springer; 2013.
Thatje S, Casburn L, Calcagno JA. Behavioural and respiratory response of the shallow-water hermit crab Pagurus cuanensis to hydrostatic pressure and temperature. J Exp Mar Biol Ecol. 2010;390(1):22–30.
Calligari PA, Calandrini V, Ollivier J, Artero JB, Hartlein M, Johnson M, Kneller GR. Adaptation of extremophilic proteins with temperature and pressure: evidence from initiation factor 6. J Phys Chem B. 2015;119(25):7860–73.
Place SP, Hofmann GE. Constitutive expression of a stress-inducible heat shock protein gene, hsp70, in phylogenetically distant Antarctic fish. Polar Biol. 2005;28(4):261–7.
Schmid B, Klumpp J, Raimann E, Loessner M, Stephan R, Tasara T. Role of cold shock proteins in growth of listeria monocytogenes under cold and osmotic stress conditions. Appl Environ Microbiol. 2009;75(6):1621–7.
Phadtare S, Alsina J, Inouye M. Cold-shock response and cold-shock proteins. Curr Opin Microbiol. 1999;2(2):175–80.
Yayanos AA, Pollard EC. A study of the effects of hydrostatic pressure on macromolecular synthesis in Escherichia coli. Biophys J. 1969;9(12):1464–82.
Han Q, Keesing JK, Liu D. A review of sea cucumber aquaculture, ranching, and stock enhancement in China. Rev Fisheries Sci Aquaculture. 2016;24(4):326–41.
Miller AK, Kerr AM, Paulay G, Reich M, Wilson NG, Carvajal JI, Rouse GW. Molecular phylogeny of extant Holothuroidea (Echinodermata). Mol Phylogenet Evol. 2017;111:110–31.
Li Y, Kong X, Chen J, Liu H, Zhang H. Characteristics of the Copper, Zinc Superoxide Dismutase of a Hadal Sea Cucumber (Paelopatides sp.) from the Mariana Trench. Marine Drugs. 2018;16(5):169.
Chen J, Liu H, Cai S, Zhang H. Comparative transcriptome analysis of Eogammarus possjeticus at different hydrostatic pressure and temperature exposures. Sci Rep. 2019;9(1):3456.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Peck LS, Webb KE, Bailey DM. Extreme sensitivity of biological function to temperature in Antarctic marine species. Funct Ecol. 2004;18(5):625–30.
Feder ME. Temperature Biology of Animals. Science. 1988;241:1108.
Peck LS. Feeding, metabolism and metabolic scope in Antarctic marine ectotherms. In: Society for Experimental Biology Seminar; 1998.
Peck LS, Morley SA, Pörtner HO, Clark MS. Thermal limits of burrowing capacity are linked to oxygen availability and size in the Antarctic clam Laternula elliptica. Oecologia. 2007;154(3):479–84.
Qiu XB, Shao YM, Miao S, Wang L. The diversity of the DnaJ/Hsp40 family, the crucial partners for Hsp70 chaperones. Cell Mol Life Sci. 2006;63(22):2560–70.
Denisov IG, Makris TM, Sligar SG, Schlichting I. Structure and chemistry of cytochrome P450. Chem Rev. 2005;105(6):2253–77.
Suzuki N, Suzuki S, Saito T. IRAKs: key regulatory kinases of innate immunity. Curr Med Chem - Anti-Inflammatory Anti-Allergy Agents. 2005;4(1):13–20.
Peralta ER, Martin BC, Edinger AL. Differential effects of TBC1D15 and mammalian Vps39 on Rab7 activation state, lysosomal morphology, and growth factor dependence. J Biol Chem. 2010;285(22):16814–21.
Croniger CM, Millward C, Yang J, Kawai Y, Arinze IJ, Liu S, Harada-Shiba M, Chakravarty K, Friedman JE, Poli V, et al. Mice with a deletion in the gene for CCAAT/enhancer-binding protein beta have an attenuated response to cAMP and impaired carbohydrate metabolism. J Biol Chem. 2001;276(1):629–38.
Bao Y, Jr TLD, Chen YT. Human glycogen debranching enzyme gene (AGL): complete structural organization and characterization of the 5′ flanking region. Genomics. 1996;38(2):155–65.
Ban N, Matsumura Y, Sakai H, Takanezawa Y, Sasaki M, Arai H, Inagaki N. ABCA3 as a lipid transporter in pulmonary surfactant biogenesis. J Biol Chem. 2007;282(13):9628–34.
Rossi A, Riccio A, Coccia M, Trotta E, La Frazia S, Santoro MG. The proteasome inhibitor bortezomib is a potent inducer of zinc finger AN1-type domain 2a gene expression: role of heat shock factor 1 (HSF1)-heat shock factor 2 (HSF2) heterocomplexes. J Biol Chem. 2014;289(18):12705–15.
Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36(Database issue):D480–4.
R Core Team. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2018.
Bao J, Dong SL, Tian XL, Wang F, Gao QF, Dong YW. Metabolic rates and biochemical compositions of Apostichopus japonicus (Selenka) tissue during periods of inactivity. Chin J Oceanol Limnol. 2010;28(2):218–23.
Pickard GL. Descriptive physical oceanography: an introduction. Oxford: Pergamon; 1990.
New P, Brown A, Oliphant A, Burchell P, Smith A, Thatje S. The effects of temperature and pressure acclimation on the temperature and pressure tolerance of the shallow-water shrimp Palaemonetes varians. Mar Biol. 2014;161(3):697–709.
Zhang Y, Sun J, Chen C, Watanabe HK, Feng D, Zhang Y, Chiu JMY, Qian PY, Qiu JW. Adaptation and evolution of deep-sea scale worms (Annelida: Polynoidae): insights from transcriptome comparison with a shallow-water species. Sci Rep. 2017;7:46205.
Deshaies RJ, Joazeiro CA. RING domain E3 ubiquitin ligases. Annu Rev Biochem. 2009;78:399–434.
Mellman I. Endocytosis and molecular sorting. Annu Rev Cell Dev Biol. 1996;12(1):575–625.
Abram A, Rob VH, Kristof V, Michiels CW. An SOS response induced by high pressure in Escherichia coli. J Bacteriol. 2004;186(18):6133–41.
Yang A, Zhou Z, Pan Y, Jiang J, Ying D, Guan X, Sun H, Shan G, Zhong C. RNA sequencing analysis to capture the transcriptome landscape during skin ulceration syndrome progression in sea cucumber Apostichopus japonicus. BMC Genomics. 2016;17(1):459.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng QD, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–U130.
Davidson NM, Oshlack A. Corset: enabling differential gene expression analysis for de novo assembled transcriptomes. Genome Biol. 2014;15(7):410.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323.
Chen H. VennDiagram: generate high-resolution Venn and Euler plots; 2014.
Kolde R: Pheatmap: Pretty Heatmaps. R package version 1.0.10. (2015).
Oksanen J, Blanchet G, Kindt R, Legendre P, Minchin PR, O’Hara RB, Simpson GL, Solymos P, Stevens M, Wagner H: Vegan: Community Ecology Package. R package version 2.5–2. (2018).
Wickham H. ggplot2: Elegant Graphics for Data Analysis. New York: Springer-Verlag; 2016.
The author J.C. would like to acknowledge Xue Kong, Ruoyu Liu, Simin Xu and Zhiyan Chen for their assistance on this study.
This work was supported by Strategic Priority Research Program of CAS (XDA22040502), The National Key Research and Development Program of China (2016YFC0304905, 2018YFC0309804), and Major scientific and technological projects of Hainan Province (ZDKJ2016012). The funding bodies had no role in the design of the study, the collection, analysis, and interpretation of data, and in writing the manuscript.
Ethics approval and consent to participate
The sea cucumbers (A. japonicus) used here were commercially cultured, and no wild individuals were involved in this study. All experiments were conducted in accordance with the legal requirements of the local government.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1 Figure S1.
Results of KEGG enrichment.
Additional file 2 Table S1.
Quality of sequencing. P0.1: experimental group incubated at atmospheric pressure; P15: experimental group incubated at 15 MPa; P25: experimental group incubated at 25 MPa.
Additional file 3 Table S2.
Swiss-Prot annotation of linearly related DEGs. DEGs: differentially expressed genes.
Additional file 4 Table S3.
Swiss-Prot annotation of pressure-sensitive DEGs. DEGs: differentially expressed genes.
Additional file 5 Table S4.
Swiss-Prot annotation of high-pressure-induced DEGs. DEGs: differentially expressed genes.
Additional file 6 Table S5.
DEGs involved in the pathway clathrin-dependent endocytosis. DEGs: differentially expressed genes.
Additional file 7 Table S6.
Gene family analysis results of linearly related DEGs. DEGs: differentially expressed genes.
Additional file 8 Table S7.
Gene family analysis results of pressure-sensitive DEGs. DEGs: differentially expressed genes.
Additional file 9 Table S8.
Gene family analysis results of high-pressure-induced DEGs. DEGs: differentially expressed genes.
Additional file 10: Table S9.
Information of the primers used in quantitative real-time reverse transcription-PCR (qPCR) analysis and the PCC between RNA-seq and qPCR results (with cytb and ß-actin as internal control). PCC: Pearson correlation coefficients.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Liang, L., Chen, J., Li, Y. et al. Insights into high-pressure acclimation: comparative transcriptome analysis of sea cucumber Apostichopus japonicus at different hydrostatic pressure exposures. BMC Genomics 21, 68 (2020). https://doi.org/10.1186/s12864-020-6480-9
- Hydrostatic pressure
- Differentially expressed gene
- Sea cucumber