PacBio Isoform Sequencing and Illumina RNA Sequencing Provide Novel Insights on Responses to Acute Heat Stress in Apostichopus japonicus Coelomocytes

Introduction

Climate change has been predicted to cause a significant increase in global sea surface temperatures (Collins et al., 2013) and alter amplitudes of diurnal temperature fluctuations (Easterling et al., 2000; Meehl et al., 2000). Over the past several decades, rising sea water temperature has impinged upon all biochemical and physiological processes of marine organisms (Somero, 2010). Heat waves and seasonal variation linked with global warming are causing frequent fluctuations in water temperature in aquatic environments, causing further challenges for marine organisms, particularly those in shallow water locations (Frölicher et al., 2018; Oliver et al., 2018).

Sea cucumbers play an important role in marine ecosystem functionality through nutrient recycling, bioturbation, benthic productivity, and meiofauna population dynamics (Schneider et al., 2011; MacTavish et al., 2012; Oh et al., 2017). To date, more than 1,250 species of sea cucumbers have been identified in benthic areas of both shallow and deep seas world-wide and approximately 20 of them are known to be edible (Jo and Park, 2016). The sea cucumber Apostichopus japonicus, has a high nutritional value and important pharmacological properties. The species has long been popular as a food and folk medicine in Asian markets and is widely cultured in China as one of the most important economic species (Chang et al., 2009; Bordbar et al., 2011). As a poikilothermic marine animal and a temperate species, A. japonicus is very sensitive to temperature variation and shows a limited ability to adapt to environmental temperature changes (Xu et al., 2016). Over a long evolutionary process, sea cucumbers have acquired a high temperature-induced special physiological behavior, estivation, to aid survival and reproduction (Yang et al., 2005). Since the late summer of 2013, continuous extreme high temperatures (up to 30°C) occurred frequently during the summer in Shandong and Hebei provinces of China and have resulted in mass mortality of sea cucumbers cultured in semi-open regions, seriously restricting the sustainable development of sea cucumber aquaculture (Liu et al., 2016).

In previous studies, heat stress (HS) responses of A. japonicus were investigated and it was demonstrated that acute and/or long-term exposure to HS results in extensive changes in gene transcription (Meng et al., 2009; Zhao et al., 2011; Xu et al., 2014, 2016, 2018; Li and Xu, 2018; Huo et al., 2019; Xu and Wang, 2019). Although transcriptome analyses have greatly broadened our knowledge of the HS response in sea cucumbers, the available data are still not enough because of the limited annotations and different sampling times and tissues reported in different studies. Most previous transcriptome analyses of the HS response in sea cucumbers were performed using intestine, respiratory tree, or muscle tissues (Xu et al., 2018; Huo et al., 2019; Li et al., 2019). However, information on the HS responses of coelomocytes at the transcriptional level are still limited. Coelomocytes are suspended in coelomic fluid and function similarly to the hemocytes of other invertebrates. Previous studies have reported that immune function in sea cucumbers is mainly executed by coelomocytes, and that humoral immunity is based on a variety of macromolecules in the coelomic cavity that are secreted by coelomocytes (Wang et al., 2009; Ramírez-Gómez et al., 2010; Xue et al., 2015).

In addition to transcriptional changes, post-transcriptional modifications generated by alternative splicing (AS) also play important roles in the HS response (Fujikake et al., 2005; Liu et al., 2013; Tian et al., 2019). HS-induced AS has been suggested to be a common feature among heat shock factor (hsf) transcripts, being reported from Drosophila melanogaster (Fujikake et al., 2005), zebrafish (Danio rerio) (Råbergh et al., 2000), and catfish (Ictalurus punctatus) (Tan et al., 2019). However, to date, there have been no studies of the potential role of AS in the HS response of sea cucumbers. Therefore, the strategy of applying “Next-generation” sequencing (NGS) technology and single-molecular real-time (SMRT) sequencing combined with genome reference to analyze AS under HS was worth exploring.

Currently, SMRT sequencing has been confirmed as a valuable tool for the improvement and identification of alternative isoforms produced by AS with higher accuracy by optimizing gene structures of genome annotation. In the present study, we aimed to investigate transcriptional and post-transcriptional changes generated by AS in response to HS in sea cucumber coelomocytes. Both Pacific Biosciences (PacBio) SMRT and Illumina-based sequencing were performed to construct transcriptomes from control and 48 h heat-shocked sea cucumbers. Our results help to elucidate the underlying mechanisms of the acute HS response in sea cucumber coelomocytes and provide important clues for further study of the molecular physiology of A. japonicus.

Materials and Methods

Sample Collection for PacBio Sequel Sequencing

One hundred and twenty A. japonicus adults (100–120 g) were collected from an aquaculture farm in Weihai (Shandong, China) in April, 2019. The animals were acclimated in a seawater aquarium system at 15°C and 30 salinity for 2 weeks and fed with an artificial food compound (30% gulfweed powder, 20% murina powder, and 50% sea mud) once a day. All animals were fasted for 2 days prior to experiments.

After acclimation, animals were randomly divided into control (CO, 15°C) or experimental groups (48 h HS, 30°C) with three parallel tanks for each group. Dissolved oxygen (DO) was maintained at 8 mg/L in both groups. Following exposure, coelomocytes were collected immediately and frozen rapidly in liquid nitrogen from all groups. All samples were stored at −80°C for subsequent analysis. All animal protocols for care, experimentation, and sacrifice had the approval of the Experimental Animal Ethics Committee of the Ocean University of China.

RNA Extraction

Total RNA was isolated from coelomocytes using Trizol RNA isolation reagent (TaKaRa, Japan) following the manufacturer’s instructions and digested with RNase-free DNase I (TaKaRa, Japan) to remove genomic DNA. RNA degradation and contamination were monitored using 1% agarose gels. The purity and concentration of RNA was measured using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Rockland, DE, United States) with OD 260/280 reading. Accurate quantification and integrity of RNA was further determined using a Qubit RNA Detection Kit with a Qubit 2.0 Fluorometer (Invitrogen, CA, United States) and Agilent 2100 Bioanalyzer (Agilent Technologies, CA, United States).

PacBio Library Construction and Sequencing

The total RNA samples from control (n = 3) and 48 h heat stressed (n = 3) A. japonicus coelomocytes were pooled and sequenced using the Pacific Biosciences Isoform sequencing (Iso-Seq) platform. The mRNA with polyA tails was enriched by use of Oligo (dT) magnetic beads and reverse transcribed into cDNA using a Clontech SMARTer PCR cDNA Synthesis Kit (TaKaRa, Japan). Then optimized PCR reactions were used to generate double-stranded cDNA, followed by size selection (fragments larger than 4 kb) using a Blue Pippin TM Size-Selection System (Sage Science, MA, United States) as described by PacBio (PN 100-092-800-03). The generated cDNA was then re-amplified using PCR to obtain sufficient cDNA. Full-length cDNA damage/terminal repair, ligation of the SMRT dumbbell-type adapters and SMRTbell template preparation were then carried out. The SMRTbell template was bound with sequencing primer and polymerase and sequenced on the PacBio sequel platform according to the manufacturer’s protocol. The Iso-seq datasets generated are available from the NCBI Sequence Read Archive database (SRA) with accession number PRJNA785124.

SMRT Sequencing Data Processing

The raw SMRT sequencing reads were processed using PacBio SMRTlink v5.0 software. After removing the adapter and low-quality data, the circular consensus sequence (CCS) reads were generated from subreads.bam files and further classified as full-length (FL) and non-full-length (NFL) reads based on following parameters: whether the sequence contains both 5′ primer, 3′ primer, and poly A tail signals. Subsequently, the full-length non-chimeric (FLNC) reads were clustered by Iterative Clustering for Error Correction (ICE) software to generate the cluster consensus isoforms. To obtain FL polished consensus sequences, Quiver software was used to finally obtain the FL polished high quality consensus sequences (accuracy ≧ 99%), which were applied for gene structure optimization of the A. japonicus reference genome (NCBI BioProject ID PRJNA757529).

Open Reading Frame Prediction and Gene Structure Optimization

The open reading frames (ORFs) were analyzed using ANGEL (Shimizu et al., 2006) software for the isoform sequences to obtain the coding sequences (CDS), protein sequences, and UTR sequences. To optimize gene structure of the A. japonicus reference genome, each isoform from FL polished high-quality consensus sequences was compared with the existing gene models of A. japonicus genome annotation. Known isoforms were 5′UTR or 3′UTR extended or shortened. Other isoforms were further classified into span isoforms, exonic overlap isoforms and intron isoforms based on their exon structure (splicing junctions).

RNA-seq Library Construction and Sequencing

Total RNA from nine individuals per group (CO and HS) was pooled into three biological replicates (three individuals each) and a total of six libraries were constructed. The Illumina libraries were generated using NEBNext^® Ultra™ RNA Library Prep Kit (E7530L) (NEB, United States) following manufacturer’s recommendations, and index codes were added to attribute sequences to each sample. The library preparations were sequenced on an Illumina Novaseq6000 platform and paired-end reads were generated. Then, the raw reads were processed using the NGS QC Toolkitv2.3.3 (Patel and Jain, 2012). The first five bases from the 5′ end of the read were trimmed, and low quality bases > 20% or ambiguous bases > 1% were removed. Cleaned reads were obtained for the following analyses. The RNA-seq datasets generated are available from the NCBI SRA with accession number PRJNA687597.

Differentially Expressed Genes, Enrichment Analysis, and Gene Set Enrichment Analysis

The clean reads from the above RNA-Seq data were mapped to the optimized reference genome using HISAT2.2.4 with “-rna-strandness RF” and other parameters as default (Kim et al., 2015). The mapped reads of each sample were assembled using StringTie v1.3.1 in a reference-based approach and then counted and subsequently normalized to fragments per kilobase of transcript per million (FPKM) values to quantify expression abundance. The DESeq R package (1.10.1) (Love et al., 2014) was used to compare the differential expression between CO and HS groups. The P-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate (FDR). An FDR < 0.05 and absolute fold change≧1 were set as the threshold for significantly differentially expressed genes (DEGs).

Gene ontology (GO) analysis of the DEGs was implemented by the GOseq R packages (V1.10.0) (Young et al., 2010) based on the Wallenius non-central hypergeometric distribution. KOBAS software version 3.0.0 (Mao et al., 2005) was used to test the statistical enrichment of DEGs in kyoto encyclopedia of genes and genomes (KEGG) pathways. A q-value < 0.05 was considered as a significant change.

Gene set enrichment analysis (GSEA) was performed using software GSEA and MSigDB (Subramanian et al., 2005) to identify whether a set of genes in specific pathways showed significant differences in two groups. Briefly, we input gene expression matrix and rank genes by the SinaltoNoise normalization method. | NES| > 1, NOM p-val < 0.05 and FDR q-val < 0.25 were considered as significant.

Alternative Splicing Analysis

Alternative Splicing analysis was based on the optimized reference genome. Identification and quantification of AS events including skipped exon (SE), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), mutually exclusive exon (MXE), and retained intron (RI) were carried out using rMATS (replicate multivariate analysis of transcript splicing) (Shen et al., 2014) (version 4.0.1)¹ with default parameters. Choosing an FDR of <0.05, significant splicing events were determined.

Validation Experiments

Quantitative real-time PCR (qPCR) analysis was used to verify the differentially expressed transcripts. Total RNA was isolated from coelomocytes of A. japonicus from the CO and HS groups as described as above. cDNA was synthesized using the PrimeScript™ RT reagent kit with gDNA Eraser (TaKaRa, Japan, Cat. No. RR047A) following the manufacturer’s instructions. All transcript-specific primers for qPCR are listed in Supplementary Table 1. SYBR^® Premix Ex Taq (Tli RNaseH Plus) was used for qPCR (TaKaRa, Japan, Cat. No. RR420A). Each PCR reaction system consisted of 10.0 μl SYBR^® Premix Ex Taq (Tli RNaseH Plus), 0.4 μl of each forward and reverse primer (10 μM), 0.4 μl ROX Reference Dye, 2 μl cDNA, and 6.8 μl dH₂O in a total volume of 20 μl. Then, qPCR was performed on a StepOnePlus™ Real-Time PCR machine (Applied Biosystems, CA, United States). The standard procedure for amplification was as follows: 95°C for 30 s and 40 cycles of 95°C for 5 s, 57°C for 30 s. β-actin and β-tubulin were used as control genes for internal standardization as previous reported (Zhao et al., 2014; Xu et al., 2018; Huo et al., 2019). The 2^–ΔΔCT method was applied to analyze transcript expression levels. One-way ANOVA was conducted followed by Duncan’s multiple tests with SPSS 19.0 software. Significant difference was set as P ≦ 0.05.

To verify the accuracy of AS events, total RNA was reverse-transcribed into cDNA using a PrimeScript RT reagent kit (TaKaRa, Japan, Cat. No. RR047A) following the manufacturer’s instructions. Transcript-specific primers were designed to span the predicted splicing events using Primer 6 software (Supplementary Table 1). PCR conditions were 95°C for 2 min and 35 cycles of 95°C for 15 s, Tm for 15 s, 72°C for 3 min, and 72°C for 5–10 min that depended on the predicted product size. PCR products were monitored on 1% agarose gels stained by GelStain and verified by Sanger sequencing.

Results

Apostichopus japonicus Transcriptome Sequencing Using PacBio Iso-Seq and Reference Genome Optimization

A total of 8,949,487 subreads with an average length of 2,443 were generated. After clustering, 479,677 CCSs were obtained, of which 263,473 full length non-chimeric (FLNC) sequences were further classified containing 5′ primer, 3′ primer, and poly (A) tails (Table 1). To further improve the sequence quality, the NFL reads were used to polish the above cluster consensus isoforms using Quiver software to obtain FL polished high quality consensus sequences (accuracy≧99%). The length distribution of the subreads, CCSs, FLNC reads and polished consensus reads are shown in Supplementary Figure 1. A total of 20,162 polished high quality consensus sequences were applied for gene structure optimization of the A. japonicus reference genome, and 380 genes in the reference genome were optimized with detailed information shown in Supplementary Table 2.

TABLE 1

Table 1. Summary of PacBio sequencing data in A. japonicus.

Differential Alternative Splicing Events Analysis in Response to Heat Stress

In the present study, DAS events were identified by comparing control and HS groups. A total of 1,224 and 1,251 DAS events (Figure 1) were identified using the splice junction counts (JC only) and reads on target and junction counts (JCEC) as the input for rMATS. Skipping exon (SE) dominated in both count methods (61.68 and 61.23%, respectively). Detailed information of DAS events is shown in Supplementary Table 3.

FIGURE 1

Figure 1. Analysis of alternative splicing (AS) events after HS. (A) The number of differentially expressed AS events after 48 h heat shock (48 h_HS) as compared with control (CO) sea cucumbers using the splice junction counts (JC only) as the input for rMATS. (B) The number of differentially expressed AS events in 48 h_HS vs. CO using reads on target and junction counts (JCEC) as the input for rMATS. SE, skipped exon; A5SS, alternative 5′ splice site; A3SS, alternative 3′splice site; MXE, mutually exclusive exon, and RI, retained intron.

Differentially Expressed Genes Analysis in Response to Heat Stress

A total of 1375 DEGs were identified between the control and 48 h HS groups. As shown in Figure 2, 671 genes were up-regulated and 704 genes were down-regulated; FDR was <0.05 and | log2fold-change| ≧1). The expression of DEGs with FDR < 0.01 and | log2fold-change| ≧2.5 from the six samples were clustered into two distinct groups by hierarchical clustering (Supplementary Figure 2). Additional details of the DEGs are presented in Supplementary Table 4.

FIGURE 2

Figure 2. Volcano scatter plots of differentially expressed genes (48 h_HS vs. CO). Red points represent significantly up-regulated genes (false discovery rate, FDR < 0.05 and | log₂fold-change| ≧1). Blue points represent significantly down-regulated genes (FDR < 0.05 and | log₂fold-change| ≧1).

GO, Kyoto Encyclopedia of Genes and Genomes Enrichment Analysis of DEGs and GSEA

Gene ontology enrichment analysis showed that 543 DEGs were significantly enriched in 83 GO categories (Supplementary Table 5) (q-value < 0.05) including 127 DEGs assigned to “response to stress” (BP, GO:0006950), 77 DEGs assigned to “immune response and immune system process” (BP, GO:0006955 and GO: 0002376, respectively), and 23 DEGS assigned to “protein folding” (BP, GO:0006457). The top 20 GO terms are shown in Figure 3.

FIGURE 3

Figure 3. Circular map of the top 20 gene ontology (GO) classifications of differentially expressed genes (DEGs) between 48 h_HS and CO groups. The first (outer) circle shows different Ontology: Biological Process (yellow) and Cellular Component (green). The second circle (orange) shows the background gene numbers enriched in each GO term, the higher the number, the longer the bar; the depth of color represents the q-value, the smaller of the q-value, the deeper the color. The third circle represents the number of up-regulated (dark purple) and down-regulated (light purple) DEGs enriched in each GO term. The fourth (inner, white) circle shows the rich factor value of each GO term (the number of DEGs enriched in this GO term divided by all the numbers of background genes enriched in this GO term). In the fourth circle, each grid of the background grid line represents 0.1.

KEGG enrichment analysis showed that 9 pathways were significantly enriched (q-value < 0.05) after 48 h of HS including “Protein processing in endoplasmic reticulum,” “Prostate cancer,” “PI3K-Akt signaling pathway,” “Antigen processing and presentation,” “Parathyroid hormone synthesis, secretion, and action,” “Spliceosome,” “Systemic lupus erythematosus,” “Thyroid hormone synthesis,” and “Herpes simplex infection” (Figure 4 and Supplementary Table 6), most of which were further classified into four classes: (a) folding, sorting, and degradation, (b) immune and infectious diseases, (c) signal transduction, and (d) transcription. Among them, all 41 genes connected with Protein processing in the endoplasmic reticulum were significantly up-regulated. The top 20 KEGG pathways are shown in Figure 4.

FIGURE 4

Figure 4. Circular map of top 20 KEGG pathways enriched in the differentially expressed genes (DEGs) between 48 h_HS and CO groups. The first (outer) circle: different colors represent different KEGG A-classes: genetic information processing, environmental information processing, cellular processes, organismal systems, and human diseases. The second circle: the number shows the background gene numbers enriched in this pathway, the higher the number, the longer of the bar; depth of color represents the q-value, the smaller the q-value, the deeper of the color. The third circle: dark purple bar represents the number of up-regulated DEGs enriched in this pathway, and the light purple bar represents the number of down-regulated DEGs enriched in this pathway. The fourth (inner) circle shows the rich factor value of each pathway (the number of DEGs enriched in each pathway divided by all the numbers of background genes enriched in this pathway). In the fourth circle, each grid of the background grid line represents 0.1.

The results of GSEA further showed that two significantly enriched pathways were confirmed to be activated in response to HS: “Protein processing in endoplasmic reticulum” and “Antigen processing and presentation” (Figure 5).

FIGURE 5

Figure 5. GSEA applied to validate the significantly enriched pathways.

Validation of RNA-seq Results and Representative DAS Gene

To validate the accuracy of the transcriptome quantification, 12 up-regulated and 7 down-regulated genes (all from the enriched KEGG pathways) were selected and quantified using qRT-PCR. The results of quantitative real-time PCR showed that all up- and down-regulated genes showed the same trends as the high-throughput sequencing data (Supplementary Figure 3). Hsf1 was selected as the representative DAS gene, being the only heat shock transcription factor gene that was significantly alternatively spliced under HS. The differentially spliced exon was located at exon 10 (Figure 6) resulting in two hsf1 transcript isoforms: Ajhsf1-1 and Ajhsf1-2. In the control group of sea cucumbers, the inclusion level of exon 10 was 86.3% but, by contrast, the exon inclusion level was only 7.3% in sea cucumbers under HS (Figure 6 and Supplementary Table 3).

FIGURE 6

Figure 6. The alignment of exon 10 and flanking regions in Ajhsf1-1 and Ajhsf1-2 isoforms. Exon 10 is spliced out in the Ajhsf1-1 isoform, whereas it is retained in the Ajhsf1-2 isoform. The boxes represent the exons and the line represents the introns. Bar graph showing the relative mRNA levels for the two isoforms in CO and HS groups.

Discussion

Seasonal fluctuations in water temperature are unavoidable in semi-open aquaculture areas and, with global warming, it is predicted that suboptimal high temperatures may become more frequent and severe in coastal areas (Frölicher et al., 2018; Oliver et al., 2018). In recent years, extreme weather caused by high temperatures (up to 30°C) has been considered to be an added factor for sea cucumber mortality syndrome in the summer (Liu et al., 2016). The HS response of sea cucumbers is a rapid and complicated process, potentially involving numerous genes and with differential responses to HS by different tissues (Xu et al., 2018; Huo et al., 2019; Li et al., 2019). Coelomocytes, the major component of the non-specific defense system of sea cucumbers, are involved in a series of immune responses (Xue et al., 2015). It is believed that the molecular basis of the HS response in coelomocytes is more complicated than those in intestine, muscle, or respiratory tree. The current study presents the first analysis of transcriptional and post-transcriptional changes generated by AS of genes in response to HS (48 h at 30°C) in sea cucumber coelomocytes.

Overall, a total of 1,375 DEGs identified in response to HS were significantly enriched in nine pathways including protein processing in endoplasmic reticulum, PI3K-AKT signaling pathway, antigen processing and presentation, and the spliceosome. Two of these were confirmed by GSEA to be activated by HS: (a) protein processing in endoplasmic reticulum and (b) antigen processing and presentation.

Heat Shock Response and Endoplasmic Reticulum Stress

The HSP family has been widely confirmed to play an important role as “assistants” that support proper folding, assembly, or disassembly of proteins and protein complexes and promote the proper refolding of abnormal or damaged proteins (Parsell and Lindquist, 1993; Velichko et al., 2013). Synthesis of heat shock proteins under HS conditions is one of the most important molecular responses to stress and is a highly conserved mechanism across phylogeny that helps to maintain cell homeostasis and protect cells and tissues from structural damage in most organisms, including aquatic ectotherms (Roberts et al., 2010). Expression of HSPs in animals exposed to HS has been well characterized in a wide variety of aquatic animals (Meng et al., 2009; Zhao et al., 2011; Xu et al., 2014, 2016). As observed in the present study, significant overexpression of several HSPs was detected in response to HS in sea cucumber coelomocytes. These included the inducible form and its constitutive counterpart of the HSP90 family members, hsp90ab1, hsp90b1, and hsp90aa1 (alias hsp90-alpha). The HSP90 protein family is essential to the folding and stabilization of cellular proteins and is implicated in a wide variety of cell processes as a response to stress (Jackson, 2013). These findings are in line with previous HS studies in aquatic animals including sea cucumbers (Liu et al., 2013; Xu et al., 2018; Huo et al., 2019; Li et al., 2019; Swirplies et al., 2019). Interestingly, sea cucumbers exposed to a moderate high temperature (26°C) showed no changes in hsp90 expression after 24–72 h HS (Huo et al., 2019), which differs from our results, but suggests that the expression of hsp90 may vary in different tissues or has a trigger temperature somewhere in between 26 and 30°C.

We also identified two cofactors, CDC37 and AHSA1, that were significantly up-regulated after 48 h HS. CDC37 is an accessory factor for HSP90, directing HSP90 to protein kinases as substrates and can also function as a molecular chaperone itself, interacting with other HSP90 co-chaperones to support protein folding (MacLean and Picard, 2003). This current finding is similar to the results of a HS study in rainbow trout (Huang et al., 2018). AHSA1 is also a HSP90 chaperone and stimulates the ATPase activity of HSP90 (Shao et al., 2016). It has been reported that HSP90 and AHSA1 are up-regulated in parallel under HS in sea cucumber intestine (Xu et al., 2018), which is consistent with our present study and suggests the potential cooperation of these two genes in response to HS.

In addition, seven members of the HSP40 family were identified and all were up-regulated in our present study, including DNAJA1, DNAJA2, DNAJB4, DNAJB6, DNAJB9, DNAJC3, and DNAJC7. These are all involved in protein translation, folding, unfolding, translocation, and degradation, primarily by stimulating the ATPase activity of HSP70 chaperones (Qiu et al., 2006). HSP70 is usually considered to be the central actor of the HS response, and elevated hsp70 mRNA expression upon HS has been reported in numerous aquatic animals including sea cucumber (Roberts et al., 2010; Li et al., 2019). Thus, in the present study, inducible forms of hsp70 were likely involved in cellular protein folding and remodeling processes within the coelomocytes of stress-exposed sea cucumbers. Strong expression of these HSP genes induced by HS suggest that these family members contribute to “thermotolerance” for sea cucumber cells by providing temporary resistance to stress until potential longer-term adaptive mechanisms can be acquired. HSP70 also cooperates with other cellular chaperones including HSP90 and HSP40 and small heat shock proteins (e.g., HSP26) together constituting a dynamic and functionally versatile chaperone network for protein folding, unfolding, regulation, targeting, aggregation, and disaggregation, in response to HS (Rosenzweig et al., 2019).

In addition, we also identified the most significantly enriched KEGG pathway in the coelomocytes of sea cucumbers in response to HS. This was “protein processing in the endoplasmic reticulum (ER)” and all 41 genes belonging to this pathway were up-regulated in comparison with controls, which is consistent with observations in other aquatic animals and other tissues in sea cucumber in response to HS (Guo et al., 2016; Qian and Xue, 2016; Li et al., 2017; Huang et al., 2018; Xu et al., 2018). It is well known that the ER has strict quality-control systems that ensure the correct folding of newly synthesized proteins and contributes importantly to maintaining cell homeostasis. An accumulation of unfolded or misfolded proteins in the ER induces the unfolded protein response (UPR), which is an evolutionarily conserved stress-detecting system used to adapt to changing environmental or stress conditions and, thereby, re-establish normal ER function (Raghunath et al., 2018). Our previous study reported the key ER chaperone response induced by HS in intestinal cells of sea cucumbers (Wang et al., 2021). A marked signature of the UPR stress response is the coordinated transcriptional up-regulation of resident ER chaperones such as the 78 kDa glucose-regulated protein (Grp78), heat shock protein 90 kDa beta (Grp94), and PDI family members (Lee, 2001; Ma and Hendershot, 2002; Dunlop et al., 2018). The overexpression of ER chaperones in our present study suggested that they may help to correct the folding and assembly of nascent proteins and prevent the transport of proteins or protein subunits that are folded incorrectly in sea cucumbers in response to HS. A second characteristic of the UPR is inhibition of protein synthesis under stress conditions (Prostko et al., 1992; Brostrom et al., 1996; Ma and Hendershot, 2002). In the present study, two eukaryotic translation initiation factor 2-alpha kinase (EIF2AK) proteins (EIF2AK1 and EIF2AK3) were significantly overexpressed in our HS group. EIF2AKs share a common physiological role of phosphorylating eIF2α to inhibit protein synthesis at the level of translation initiation in response to various stress conditions, including oxidative stress, heme deficiency, osmotic shock, and heat shock (Krishna and Kumar, 2018). This indicates another conserved regulatory mechanism present in sea cucumbers and activated in response to HS.

Through a process called ER quality control, proteins that do not mature properly are retained in the ER and eventually targeted for ER-associated degradation (ERAD), which is the third major pathway of the UPR and disposes of misfolded proteins from the ER lumen in response to stress (Brodsky et al., 1999; Nishitoh, 2012; Dunlop et al., 2018). In A. japonicus, we also observed a significant up-regulation of endoplasmin in response to HS, an important protein of the ERAD process (Christianson et al., 2008). ERAD also involves chaperones of the HSP family that act in concert with ubiquitin ligases (Goldberg, 2003).

Hence, in response to HS, the expression of genes from the HSP70, HSP90, and HSP40 families were upregulated in sea cucumber coelomocytes, likely to promote degradation of misfolded proteins and to maintain cell homeostasis. These results provide important evidence of an evolutionarily ancient role for the ER-induced UPR under HS in echinoderms. Further research will be required to confirm this in the future. Additionally, the hypoxia-induced gene, hyou1, was also upregulated by HS and is involved in “protein processing in the ER” pathway. This suggests that thermal and hypoxia tolerance may positively related with each other (Pörtner and Peck, 2010), as has previously been reported for heat-stressed rainbow trout (Huang et al., 2018).

Immune Response

Several previous studies have focused on immune responses within different high temperature stress and/or duration conditions in sea cucumber tissues (intestine and muscle) (Huang et al., 2018; Huo et al., 2019; Li et al., 2019). Here we reported that three enriched pathways were linked to the immune system in response to HS in sea cucumber coelomocytes: “Antigen processing and presentation,” “Systemic lupus erythematosus,” and “Herpes simplex infection.” GSEA confirmed the activated pathway “antigen processing and presentation” in response to HS, which was also reported in muscle tissues of sea cucumbers in response to HS (Li et al., 2019). The GO/pathway analysis also indicated that several HSP transcripts were associated with this pathway (e.g., hsp70, hsp90α) in response to HS. It has been reported that HSPs in the cytosol such as HSP70 and HSP90, are involved in antigen presentation and cross-presentation by delivering chaperones, antigenic peptides to MHC-I molecules (Tsan and Gao, 2009). Our results suggest a potential role of HSPs in the immune system of sea cucumber coelomocytes under HS. Similar results have also been reported in some fish species such as rainbow trout Oncorhynchus mykiss (Huang et al., 2018) and Atlantic salmon Salmo salar (Beemelmanns et al., 2021).

The complement system is a major effector of the immune response and triggers many immunoregulatory mechanisms (Carroll, 2008). It plays a fundamental role in innate immune responses in invertebrates, and echinoderms have been confirmed to have a simplified complement system (Smith et al., 2010). In the current study, the complement component C3-2 (AjC3-2) was significantly down-regulated after 48 h HS (30°C), which is contrary to observations in sea cucumber intestine after 48 h HS at 26°C, suggesting that 30°C stress may be a physiological upper limit temperature for sea cucumbers and lead to impaired immune function of coelomocytes under HS.

Signal Transduction

When ectotherms are exposed to HS, mitochondrial respiration rate is increased and this can result in oxidative stress and cellular damage (Heise et al., 2006; Cheng et al., 2015; Banh et al., 2016). In our study, heat stressed sea cucumbers appeared to have an enhanced oxidative stress response. Members of the forkhead box, class O (FoxO) transcription factor family are important regulators of the cellular stress response and promote antioxidant defenses (Klotz et al., 2015). The vascular endothelial growth factor receptor (VEGFR) has been reported to be activated under oxidative stress and trigger downstream cell survival signaling cascades (Wang et al., 2004). The cyclic AMP-responsive element-binding protein 3-like protein 4 (CREB3L4) is a transcription factor that responds to oxidative stress signals (Pregi et al., 2017). These three up-regulated genes that are potentially involved in oxidative stress responses were identified in sea cucumber coelomocytes after HS, suggesting their potentially important functions in the regulation of signaling to activate responses that can deal with injury caused by HS.

Post-transcriptional Regulation

The spliceosome is a large ribonucleoprotein complex involved in splicing the coding regions (exons) of precursor messenger RNAs (pre-mRNAs) with the concomitant excision of intervening non-coding introns. In sea cucumber coelomocytes, 15 genes were significantly enriched in the spliceosome pathway after HS. Similar results were found in other aquatic animals (Qian and Xue, 2016; Huang et al., 2018). Serine-arginine-rich (SR) and heterogeneous nuclear ribonucleoprotein (hnRNP) proteins are best known as the major splicing factors involved in AS regulation (Lee and Rio, 2015). In our study, three hnRNP genes (hnRNPA1-like, hnRNPA3, and hnRNPK-like), and three SR genes (srsf2, srsf3-like, and srsf7) were identified as responding to HS. Most of these two classes of genes varied oppositely in sea cucumber coelomocytes in response to HS. Indeed, previous studies have observed that these two classes of genes often have antagonistic effects on splice-site usage. Exonic splicing enhancers (ESEs) are often bound by SR proteins whereas exonic splicing silencers (ESSs) are typically bound by hnRNP proteins (Matera and Wang, 2014). Interestingly, a negative correlation between hnRNPK expression (a significant decrease) and hsp70 expression (a significant increase) was also observed in response to HS in the present study. The hnRNPK protein has also been reported to play a role in the heat shock response of cells and control expression of hsp70 mRNAs by inhibiting heat shock-induced transcriptional activity (Kim et al., 2017). Our present findings suggest that these splicing factors may organize a regulatory network of various groups of genes through AS as well as splicing inhibition to respond appropriately to a particular stress.

Our differential splicing analysis further indicated that AS may be an important potential mechanism used by sea cucumbers to deal with heat shock. AS has been reported to play a key role in enhancing regulatory capacities and proteomic complexities in eukaryotes (Ben-Dov et al., 2008) and allows organisms to rapidly modulate protein functions in response to physiological changes. Accumulating evidence indicates that AS is a key player in the response to environmental stresses including HS (Fujikake et al., 2005; Liu et al., 2013). In the DAS list, we identified three hnRNP genes (hnRNPA1-like, hnRNPA3, and hnRNPK-like) and four SR genes (srsf2, srsf3-like, srsf4-like, and srsf7). Previous studies have reported that temperature stress can dramatically alter the AS of pre-mRNAs of several SR and hnRNP genes in catfish (Tan et al., 2019).

Increasing the expression of heat shock proteins (HSPs) is an evolutionarily conserved strategy for eukaryotes to deal with HS. This process is mediated by a master regulator: the heat shock transcription factor, HSF1. HSF1 exists as a monomer in unstressed cells but under stress conditions forms an active trimer that can bind to promoter regions of heat shock genes (Anckar and Sistonen, 2011; Hentze et al., 2016). Tan et al. (2019) reported that the inclusion level of exon 9a in catfish hsfs1 decreased significantly from 100 to 2.7% after heat shock. In accordance, we also found that the inclusion level of exon 10 in Ajhsf1 was significantly reduced from 86.3 to 7.3% after HS. Research on mammals demonstrated that HSF1 isoform ratios are a major factor in the regulation of expression of heat shock protein genes and exon splicing. The HSF1 isoform may be produced to facilitate the activation of the homotrimer and promote the expression of HSPs after HS (Neueder et al., 2014).

Conclusion

The present study focused on transcriptional and post-transcriptional (AS) regulatory levels, combining PacBio Iso-seq and RNA-seq to reveal acute responses to HS in sea cucumber coelomocytes. Strong numbers of up-regulated DEGs were enriched in the protein processing function of the endoplasmic reticulum, providing key evidence of the involvement of an evolutionarily conserved signal transduction mechanism and stress-detecting system in response to HS in echinoderms. Further analysis identified an abundance of RNA-splicing related genes that displayed AS and confirmed that the inclusion level of splicing out an exon in Ajhsf1 was significantly decreased after HS, suggesting that AS is an important potential mechanism used by sea cucumbers in response to HS. Although satisfactory answers to these findings obtained from high-throughput screening will require further functional confirmation, our study improves understanding and opens up new perspectives on the regulatory mechanisms used by echinoderms to respond to environmental stress.

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: NCBI [accession nos: PRJNA785124, PRJNA757529, and PRJNA687597].

Ethics Statement

This animal study was reviewed and approved by the Experimental Animal Ethics Committee of the Ocean University of China.

Author Contributions

MC originally conceived the study. YW, YY, and XC conducted the experiments. MC, YW, and YY contributed to the acquisition of data. MC, YW, YY, and KS contributed to the interpretation of the data. MC wrote the first draft of the manuscript. KS revised the manuscript. All authors contributed to manuscript read and approved the submitted version.

Funding

This work was supported by National Key Research and Development Program of China (2019YFD0902103) and National Natural Science Foundation of China (31972767).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.815109/full#supplementary-material

Supplementary Figure 1 | The length distribution of the subreads, CCSs, FLNC reads and polished consensus reads.

Supplementary Figure 2 | Hierarchical clustering of DEGs with FDR (false discovery rate) <0.01 and | log2fold-change| ≧2.5 among the six sample libraries.

Supplementary Figure 3 | Confirmation of selective DEGs by qPCR method.

Supplementary Table 1 | Primers designed for alternative splicing and qPCR validation.

Supplementary Table 2 | Detailed information of optimized genes in the reference genome based on PacBio Iso-Seq.

Supplementary Table 3 | Detailed information of differential alternative splicing (DAS) events.

Supplementary Table 4 | Detailed information of DEGs with false discovery rate (FDR) <0.01 and | log2fold-change| ≧2.5.

Supplementary Table 5 | Detailed information on significantly enriched GO terms.

Supplementary Table 6 | Detailed information on significantly enriched pathways.

Footnotes

1. ^ http://rnaseq-mats.sourceforge.net/index.html

References

Anckar, J., and Sistonen, L. (2011). Regulation of HSF1 function in the heat stress response: implications in aging and disease. Annu. Rev. Biochem. 80, 1089–1115. doi: 10.1146/annurev-biochem-060809-095203

PubMed Abstract | CrossRef Full Text | Google Scholar

Banh, S., Wiens, L., Sotiri, E., and Treberg, J. R. (2016). Mitochondrial reactive oxygen species production by fish muscle mitochondria: potential role in acute heat-induced oxidative stress. Comp. Biochem. Physiol. Part B Biochem. Mol. Biol. 191, 99–107. doi: 10.1016/j.cbpb.2015.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Beemelmanns, A., Zanuzzo, F. S., Xue, X., Sandrelli, R. M., Rise, M. L., and Gamperl, A. K. (2021). The transcriptomic responses of Atlantic salmon (Salmo salar) to high temperature stress alone, and in combination with moderate hypoxia. BMC Genomics 22:261. doi: 10.1186/s12864-021-07464-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ben-Dov, C., Hartmann, B., Lundgren, J., and Valcárcel, J. (2008). Genome-wide analysis of alternative pre-mRNA splicing. J. Biol. Chem. 283, 1229–1233. doi: 10.1074/jbc.R700033200

PubMed Abstract | CrossRef Full Text | Google Scholar

Bordbar, S., Anwar, F., and Saari, N. (2011). High-value components and bioactives from sea cucumbers for functional foods—A review. Mar. Drugs 9, 1761–1805. doi: 10.3390/md9101761

PubMed Abstract | CrossRef Full Text | Google Scholar

Brodsky, J. L., Werner, E. D., Dubas, M. E., Goeckeler, J. L., Kruse, K. B., and McCracken, A. A. (1999). The requirement for molecular chaperones during endoplasmic reticulum-associated protein degradation demonstrates that protein export and import are mechanistically distinct. J. Biol. Chem. 274, 3453–3460. doi: 10.1074/jbc.274.6.3453

PubMed Abstract | CrossRef Full Text | Google Scholar

Brostrom, C. O., Prostko, C. R., Kaufman, R. J., and Brostrom, M. A. (1996). Inhibition of translational initiation by activators of the glucose-regulated stress protein and heat shock protein stress response systems. J. Biol. Chem. 271, 24995–25002. doi: 10.1074/jbc.271.40.24995

PubMed Abstract | CrossRef Full Text | Google Scholar

Carroll, M. C. (2008). Complement and humoral immunity. Vaccine 26, I28–I33. doi: 10.1016/j.vaccine.2008.11.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, Y., Feng, Z., Yu, J., and Ding, J. (2009). Genetic variability analysis in five populations of the sea cucumber Stichopus (Apostichopus) japonicus from China, Russia, South Korea and Japan as revealed by microsatellite markers. Mar. Ecol. 30, 455–461. doi: 10.1111/j.1439-0485.2009.00292.x

CrossRef Full Text | Google Scholar

Cheng, C., Yang, F., Liao, S., Miao, Y., Ye, C., Wang, A., et al. (2015). High temperature induces apoptosis and oxidative stress in pufferfish (Takifugu obscurus) blood cells. J. Therm. Biol. 53, 172–179. doi: 10.1016/j.jtherbio.2015.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Christianson, J. C., Shaler, T. A., Tyler, R. E., and Kopito, R. R. (2008). OS-9 and GRP94 deliver mutant α1-antitrypsin to the Hrd1-SEL1L ubiquitin ligase complex for ERAD. Nat. Cell Biol. 10, 272–282. doi: 10.1038/ncb1689

PubMed Abstract | CrossRef Full Text | Google Scholar

Collins, W. J., Fry, M. M., Yu, H., Fuglestvedt, J. S., Shindell, D. T., and West, J. J. (2013). Global and regional temperature-change potentials for near-term climate forcers. Atmos. Chem. Phys. 13, 2471–2485. doi: 10.5194/acp-13-2471-2013

CrossRef Full Text | Google Scholar

Dunlop, R. A., Powell, J. T., Metcalf, J. S., Guillemin, G. J., and Cox, P. A. (2018). L-serine-mediated neuroprotection includes the upregulation of the ER stress chaperone Protein Disulfide Isomerase (PDI). Neurotox. Res. 33, 113–122. doi: 10.1007/s12640-017-9817-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Easterling, D. R., Meehl, G. A., Parmesan, C., Changnon, S. A., Karl, T. R., and Mearns, L. O. (2000). Climate extremes: observations, modeling, and impacts. Science 289, 2068–2074. doi: 10.1126/science.289.5487.2068

PubMed Abstract | CrossRef Full Text | Google Scholar

Frölicher, T. L., Fischer, E. M., and Gruber, N. (2018). Marine heatwaves under global warming. Nature 560, 360–364. doi: 10.1038/s41586-018-0383-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujikake, N., Nagai, Y., Popiel, H. A., Kano, H., Yamaguchi, M., and Toda, T. (2005). Alternative splicing regulates the transcriptional activity of Drosophila heat shock transcription factor in response to heat/cold stress. FEBS Lett. 579, 3842–3848. doi: 10.1016/j.febslet.2005.05.074

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldberg, A. L. (2003). Protein degradation and protection against misfolded or damaged proteins. Nature 426, 895–899. doi: 10.1038/nature02263

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, L., Wang, Y., Liang, S., Lin, G., Chen, S., and Yang, G. (2016). Tissue-overlapping response of half-smooth tongue sole (Cynoglossus semilaevis) to thermostressing based on transcriptome profiles. Gene 586, 97–104. doi: 10.1016/j.gene.2016.04.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Heise, K., Puntarulo, S., Nikinmaa, M., Abele, D., and Pörtner, H. O. (2006). Oxidative stress during stressful heat exposure and recovery in the North Sea eelpout Zoarces viviparus L. J. Exp. Biol. 209, 353–363. doi: 10.1242/jeb.01977

PubMed Abstract | CrossRef Full Text | Google Scholar

Hentze, N., Breton, L. L., Wiesner, J., Kempf, G., and Mayer, M. P. (2016). Molecular mechanism of thermosensory function of human heat shock transcription factor Hsf1. eLife 5:e11576. doi: 10.7554/eLife.11576

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, J., Li, Y., Liu, Z., Kang, Y., and Wang, J. (2018). Transcriptomic responses to heat stress in rainbow trout Oncorhynchus mykiss head kidney. Fish Shellfish Immunol. 82, 32–40. doi: 10.1016/j.fsi.2018.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Huo, D., Sun, L., Zhang, L., Yang, H., Liu, S., Sun, J., et al. (2019). Time course analysis of immunity-related gene expression in the sea cucumber Apostichopus japonicus during exposure to thermal and hypoxic stress. Fish Shellfish Immunol. 95, 383–390. doi: 10.1016/j.fsi.2019.09.073

PubMed Abstract | CrossRef Full Text | Google Scholar

Jackson, S. E. (2013). Hsp90: structure and function. Top. Curr. Chem. 328, 155–240. doi: 10.1007/128-2012-356

PubMed Abstract | CrossRef Full Text | Google Scholar

Jo, J., and Park, C. (2016). Phylogenetic analysis of the three color variations of the sea cucumber Apostichopus japonicus. J. Aquacult. Res. Dev 7:3. doi: 10.4172/2155-9546.1000418

CrossRef Full Text | Google Scholar

Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, H. J., Lee, J. J., Cho, J. H., Jeong, J., Park, A. Y., Kang, W., et al. (2017). Heterogeneous nuclear ribonucleoprotein K inhibits heat shock-induced transcriptional activity of heat shock factor 1. J. Biol. Chem. 292, 12801–12812. doi: 10.1074/jbc.M117.774992

PubMed Abstract | CrossRef Full Text | Google Scholar

Klotz, L. O., Sánchez-Ramos, C., Prieto-Arroyo, I., Urbánek, P., Steinbrenner, H., and Monsalve, M. (2015). Redox regulation of FoxO transcription factors. Redox Biol. 6, 51–72. doi: 10.1016/j.redox.2015.06.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Krishna, K. H., and Kumar, M. S. (2018). Molecular evolution and functional divergence of eukaryotic translation initiation factor 2-alpha kinases. PLoS One 13:e0194335. doi: 10.1371/journal.pone.0194335

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, A. S. (2001). The glucose-regulated proteins: stress induction and clinical applications. Trends Biochem. Sci. 26, 504–510. doi: 10.1016/S0968-0004(01)01908-9

CrossRef Full Text | Google Scholar

Lee, Y., and Rio, D. C. (2015). Mechanisms and regulation of alternative pre-mRNA splicing. Annu. Rev. Biochem. 84, 291–323. doi: 10.1146/annurev-biochem-060614-034316

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., Fang, H., and Xu, D. (2019). Effect of seasonal high temperature on the immune response in Apostichopus japonicus by transcriptome analysis. Fish Shellfish Immunol. 92, 765–771. doi: 10.1016/j.fsi.2019.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., and Xu, D. (2018). Understanding microRNAs regulation in heat shock response in the sea cucumber Apostichopus japonicus. Fish Shellfish Immunol. 81, 214–220. doi: 10.1016/j.fsi.2018.07.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Huang, J., Liu, Z., Zhou, Y., Xia, B., Wang, Y., et al. (2017). Transcriptome analysis provides insights into hepatic responses to moderate heat stress in the rainbow trout (Oncorhynchus mykiss). Gene 619, 1–9. doi: 10.1016/j.gene.2017.03.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Sun, N., Liu, M., Liu, J., Du, B., Wang, X., et al. (2013). An autoregulatory loop controlling Arabidopsis HsfA2 expression: role of heat shock-induced alternative splicing. Plant Physiol. 162, 512–521. doi: 10.1104/pp.112.205864

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Zhou, Y., Ru, X., Zhang, M., Cao, X., and Yang, H. (2016). Differences in immune function and metabolites between aestivating and non-aestivating Apostichopus japonicus. Aquaculture 459, 36–42. doi: 10.1016/j.aquaculture.2016.03.029

CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Y., and Hendershot, L. M. (2002). The mammalian endoplasmic reticulum as a sensor for cellular stress. Cell Stress Chaperones 7, 222–229. doi: 10.1379/1466-1268(2002)007<0222:tmeraa>2.0.co;2

PubMed Abstract | CrossRef Full Text | Google Scholar

MacLean, M., and Picard, D. (2003). Cdc37 goes beyond Hsp90 and kinases. Cell Stress Chaperones 8, 114–119. doi: 10.1379/1466-1268(2003)008<0114:cgbhak>2.0.co;2

PubMed Abstract | CrossRef Full Text | Google Scholar

MacTavish, T., Stenton-Dozey, J., Vopel, K., and Savage, C. (2012). Deposit-feeding sea cucumbers enhance mineralization and nutrient cycling in organically-enriched coastal sediments. PLoS One 7:e50031. doi: 10.1371/journal.pone.0050031

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, X., Cai, T., Olyarchuk, J. G., and Wei, L. (2005). Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 21, 3787–3793. doi: 10.1093/bioinformatics/bti430

PubMed Abstract | CrossRef Full Text | Google Scholar

Matera, A. G., and Wang, Z. (2014). A day in the life of the spliceosome. Nat. Rev. Mol. Cell Biol. 15, 108–121. doi: 10.1038/nrm3742

PubMed Abstract | CrossRef Full Text | Google Scholar

Meehl, G. A., Zwiers, F., Evans, J., Knutson, T., Mearns, L., and Whetton, P. (2000). Trends in extreme weather and climate events: issues related to modeling extremes in projections of future climate change. Bull. Am. Meteorol. Soc 81, 427–436. doi: 10.1175/1520-0477(2000)081<0427:TIEWAC>2.3.CO;2

CrossRef Full Text | Google Scholar

Meng, X., Ji, T., Dong, Y., Wang, Q., and Dong, S. (2009). Thermal resistance in sea cucumbers (Apostichopus japonicus) with differing thermal history: the role of Hsp70. Aquaculture 294, 314–318. doi: 10.1016/j.aquaculture.2009.06.015

CrossRef Full Text | Google Scholar

Neueder, A., Achilli, F., Moussaoui, S., and Bates, G. P. (2014). Novel isoforms of heat shock transcription factor 1, HSF1γα and HSF1γβ, regulate chaperone protein gene transcription. J. Biol. Chem. 289, 19894–19906. doi: 10.1074/jbc.M114.570739

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishitoh, H. (2012). CHOP is a multifunctional transcription factor in the ER stress response. J. Biochem. 151, 217–219. doi: 10.1093/jb/mvr143

PubMed Abstract | CrossRef Full Text | Google Scholar

Oh, G. W., Ko, S. C., Lee, D. H., Heo, S. J., and Jung, W. K. (2017). Biological activities and biomedical potential of sea cucumber (Stichopus japonicus): a review. Fish. Aquat. Sci. 20, 1–17. doi: 10.1186/s41240-017-0071-y

CrossRef Full Text | Google Scholar

Oliver, E. C. J., Donat, M. G., Burrows, M. T., Moore, P. J., Smale, D. A., Alexander, L. V., et al. (2018). Longer and more frequent marine heatwaves over the past century. Nat. Commun. 9:1324. doi: 10.1038/s41467-018-03732-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Parsell, D. A., and Lindquist, S. (1993). The function of heat-shock proteins in stress tolerance: degradation and reactivation of damaged proteins. Annu. Rev. Genet. 27, 437–496. doi: 10.1146/annurev.ge.27.120193.002253

PubMed Abstract | CrossRef Full Text | Google Scholar

Patel, R. K., and Jain, M. (2012). NGS QC toolkit: a toolkit for quality control of next generation sequencing data. PLoS One 7:e30619. doi: 10.1371/journal.pone.0030619

PubMed Abstract | CrossRef Full Text | Google Scholar

Pörtner, H. O., and Peck, M. A. (2010). Climate change effects on fishes and fisheries: towards a cause-and-effect understanding. J. Fish Biol. 77, 1745–1779. doi: 10.1111/j.1095-8649.2010.02783.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pregi, N., Belluscio, L. M., Berardino, B. G., Castillo, D. S., and Cánepa, E. T. (2017). Oxidative stress-induced CREB upregulation promotes DNA damage repair prior to neuronal cell death protection. Mol. Cell. Biochem. 425, 9–24. doi: 10.1007/s11010-016-2858-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Prostko, C. R., Brostrom, M. A., Malara, E. M., and Brostrom, C. O. (1992). Phosphorylation of eukaryotic initiation factor (eIF) 2α and inhibition of eIF-2B in GH3 pituitary cells by perturbants of early protein processing that induce GRP78. J. Biol. Chem. 267, 16751–16754. doi: 10.1016/s0021-9258(18)41842-x

CrossRef Full Text | Google Scholar

Qian, B., and Xue, L. (2016). Liver transcriptome sequencing and de novo annotation of the large yellow croaker (Larimichthy crocea) under heat and cold stress. Mar. Genomics 25, 95–102. doi: 10.1016/j.margen.2015.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, X., Shao, Y., Miao, S., and Wang, L. (2006). The diversity of the DnaJ/Hsp40 family, the crucial partners for Hsp70 chaperones. Cell. Mol. Life Sci. 63, 2560–2570. doi: 10.1007/s00018-006-6192-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Råbergh, C. M. I., Airaksinen, S., Soitamo, A., Björklund, H. V., Johansson, T., Nikinmaa, M., et al. (2000). Tissue-specific expression of zebrafish (Danio rerio) heat shock factor 1 mRNAs in response to heat stress. J. Exp. Biol. 203, 1817–1824. doi: 10.1242/jeb.203.12.1817

CrossRef Full Text | Google Scholar

Raghunath, A., Panneerselvam, L., Sundarraj, K., and Perumal, E. (2018). Heat Shock Proteins and Endoplasmic Reticulum Stress. Cham: Springer, 39–78. doi: 10.1007/978-3-319-90725-3_3

CrossRef Full Text | Google Scholar

Ramírez-Gómez, F., Aponte-Rivera, F., Méndez-Castaner, L., and García-Arrarás, J. E. (2010). Changes in holothurian coelomocyte populations following immune stimulation with different molecular patterns. Fish Shellfish Immunol. 29, 175–185. doi: 10.1016/j.fsi.2010.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, R. J., Agius, C., Saliba, C., Bossier, P., and Sung, Y. Y. (2010). Heat shock proteins (chaperones) in fish and shellfish and their potential role in relation to fish health: a review. J. Fish Dis. 33, 789–801. doi: 10.1111/j.1365-2761.2010.01183.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenzweig, R., Nillegoda, N. B., Mayer, M. P., and Bukau, B. (2019). The Hsp70 chaperone network. Nat. Rev. Mol. Cell Biol. 20, 665–680. doi: 10.1038/s41580-019-0133-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Schneider, K., Silverman, J., Woolsey, E., Eriksson, H., Byrne, M., and Caldeira, K. (2011). Potential influence of sea cucumbers on coral reef CaCO3 budget: a case study at one tree reef. J. Geophys. Res. 116:G04032. doi: 10.1029/2011jg001755

CrossRef Full Text | Google Scholar

Shao, J., Wang, L., Zhong, C., Qi, R., and Li, Y. (2016). AHSA1 regulates proliferation, apoptosis, migration, and invasion of osteosarcoma. Biomed. Pharmacother. 77, 45–51. doi: 10.1016/j.biopha.2015.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, S., Park, J. W., Lu, Z., Lin, L., Henry, M. D., Wu, Y., et al. (2014). rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc. Natl. Acad. Sci. U.S.A. 111, E5593–E5601. doi: 10.1073/pnas.1419161111

PubMed Abstract | CrossRef Full Text | Google Scholar

Shimizu, K., Adachi, J., and Muraoka, Y. (2006). ANGLE: a sequencing errors resistant program for predicting protein coding regions in unfinished cDNA. J. Bioinform. Comput. Biol. 4, 649–664. doi: 10.1142/s0219720006002260

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, L. C., Ghosh, J., Buckley, K. M., Clow, L. A., Dheilly, N. M., Haug, T., et al. (2010). Echinoderm immunity. Adv. Exp. Med. Biol. 708, 260–301. doi: 10.1007/978-1-4419-8059-5_14

CrossRef Full Text | Google Scholar

Somero, G. N. (2010). The physiology of climate change: how potentials for acclimatization and genetic adaptation will determine ‘winners’ and ‘losers’. J. Exp. Biol. 213, 912–920. doi: 10.1242/jeb.037473

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A. 102, 15545–15550. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Swirplies, F., Wuertz, S., Baßmann, B., Orban, A., Schäfer, N., Brunner, R. M., et al. (2019). Identification of molecular stress indicators in pikeperch Sander lucioperca correlating with rising water temperatures. Aquaculture 501, 260–271. doi: 10.1016/j.aquaculture.2018.11.043

CrossRef Full Text | Google Scholar

Tan, S., Wang, W., Tian, C., Niu, D., Zhou, T., Jin, Y., et al. (2019). Heat stress induced alternative splicing in catfish as determined by transcriptome analysis. Comp. Biochem. Physiol. Part D Genomics Proteomics 29, 166–172. doi: 10.1016/j.cbd.2018.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, Y., Wen, H., Qi, X., Zhang, X., Liu, S., Li, B., et al. (2019). Characterization of full-length transcriptome sequences and splice variants of Lateolabrax maculatus by single-molecule long-read sequencing and their involvement in salinity regulation. Front. Genet. 10:1126. doi: 10.3389/fgene.2019.01126

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsan, M. F., and Gao, B. (2009). Heat shock proteins and immune system. J. Leukocyte Biol. 85, 905–910. doi: 10.1189/jlb.0109005

PubMed Abstract | CrossRef Full Text | Google Scholar

Velichko, A. K., Markova, E. N., Petrova, N. V., Razin, S. V., and Kantidze, O. L. (2013). Mechanisms of heat shock response in mammals. Cell. Mol. Life Sci. 70, 4229–4241. doi: 10.1007/s00018-013-1348-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J. F., Zhang, X., and Groopman, J. E. (2004). Activation of vascular endothelial growth factor receptor-3 and its downstream signaling promote cell survival under oxidative stress. J. Biol. Chem. 279, 27088–27097. doi: 10.1074/jbc.M314015200

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, T., Sun, Y., Jin, L., Xu, Y., Wang, L., Ren, T., et al. (2009). Enhancement of non-specific immune response in sea cucumber (Apostichopus japonicus) by Astragalus membranaceus and its polysaccharides. Fish Shellfish Immunol. 27, 757–762. doi: 10.1016/j.fsi.2009.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Zheng, Y., Chen, M., and Storey, K. B. (2021) Ultrastructural variation and key ER chaperones response induced by heat stress in intestinal cells of sea cucumber Apostichopus japonicus. J. Oceanol. Limnol. 39, 317–328. doi: 10.1007/s00343-020-9265-8

CrossRef Full Text | Google Scholar

Xu, D., Sun, L., Liu, S., Zhang, L., Ru, X., Zhao, Y., et al. (2014). Molecular cloning of heat shock protein 10 (Hsp10) and 60 (Hsp60) cDNAs and their expression analysis under thermal stress in the sea cucumber Apostichopus japonicus. Comp. Biochem. Physiol. Part B Biochem. Mol. Biol. 171, 49–57. doi: 10.1016/j.cbpb.2014.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, D., Sun, L., Liu, S., Zhang, L., and Yang, H. (2016). Molecular cloning of hsf1 and hsbp1 cDNAs, and the expression of hsf1, hsbp1 and hsp70 under heat stress in the sea cucumber Apostichopus japonicus. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 198, 1–9. doi: 10.1016/j.cbpb.2016.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, D., and Wang, X. (2019). Lysine acetylation is an important post-translational modification that modulates heat shock response in the sea cucumber Apostichopus japonicus. Int. J. Mol. Sci. 20:4423. doi: 10.3390/ijms20184423

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, D., Zhou, S., and Sun, L. (2018). RNA-seq based transcriptional analysis reveals dynamic genes expression profiles and immune-associated regulation under heat stress in Apostichopus japonicus. Fish Shellfish Immunol. 78, 169–176. doi: 10.1016/j.fsi.2018.04.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, Z., Li, H., Wang, X., Li, X., Liu, Y., Sun, J., et al. (2015). A review of the immune molecules in the sea cucumber. Fish Shellfish Immunol. 44, 1–11. doi: 10.1016/j.fsi.2015.01.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, H., Yuan, X., Zhou, Y., Mao, Y., Zhang, T., and Liu, Y. (2005). Effects of body size and water temperature on food consumption and growth in the sea cucumber Apostichopus japonicus (Selenka) with special reference to aestivation. Aquac. Res. 36, 1085–1092. doi: 10.1111/j.1365-2109.2005.01325.x

CrossRef Full Text | Google Scholar

Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11:R14. doi: 10.1186/gb-2010-11-2-r14

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, H., Yang, H., Zhao, H., Chen, M., and Wang, T. (2011). The molecular characterization and expression of heat shock protein 90 (Hsp90) and 26 (Hsp26) cDNAs in sea cucumber (Apostichopus japonicus). Cell Stress Chaperones 16, 481–493. doi: 10.1007/s12192-011-0260-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Chen, M., Wang, T., Sun, L., Xu, D., and Yang, H. (2014). Selection of reference genes for qRT-PCR analysis of gene expression in sea cucumber Apostichopus japonicus during aestivation. Chin. J. Oceanol. Limn. 32, 1248–1256.

Google Scholar