RNA Sequencing Analysis to Capture the Transcriptome Landscape during Tenderization in Sea Cucumber Apostichopus japonicus

Sea cucumber (Apostichopus japonicus) is an economically significant species in China having great commercial value. It is challenging to maintain the textural properties during thermal processing due to the distinctive physiochemical structure of the A. japonicus body wall (AJBW). In this study, the gene expression profiles associated with tenderization in AJBW were determined at 0 h (CON), 1 h (T_1h), and 3 h (T_3h) after treatment at 37 °C using Illumina HiSeq™ 4000 platform. Seven-hundred-and-twenty-one and 806 differentially expressed genes (DEGs) were identified in comparisons of T_1h vs. CON and T_3h vs. CON, respectively. Among these DEGs, we found that two endogenous proteases—72 kDa type IV collagenase and matrix metalloproteinase 16 precursor—were significantly upregulated that could directly affect the tenderness of AJBW. In addition, 92 genes controlled four types of physiological and biochemical processes such as oxidative stress response (3), immune system process (55), apoptosis (4), and reorganization of the cytoskeleton and extracellular matrix (30). Further, the RT-qPCR results confirmed the accuracy of RNA-sequencing analysis. Our results showed the dynamic changes in global gene expression during tenderization and provided a series of candidate genes that contributed to tenderization in AJBW. This can help further studies on the genetics/molecular mechanisms associated with tenderization.


Introduction
The sea cucumber (Apostichopus japonicus) is naturally found along the coasts of China, Japan, Korea, and Russia of Northeast Asia and intensively cultured in East Asian countries [1]. A. japonicus achieved the maximum single-species production value and revenue in northern China (Liaoning and Shandong provinces) with nearly 220,000 tons production [2] and estimated value of about 4 Bn USD in 2015. Traditionally, sea cucumbers are generally processed into dry products, however, in recent years new products, such as instant sea cucumbers, with more nutritional value are being developed. Nevertheless, it is difficult to maintain fine and stable texture during regular thermal processing even for instant sea cucumbers due to the distinctive physiochemical structure of A. japonicus body wall (AJBW). This particular problem has caused substantial economic losses in the seafood industry [3,4].
Tenderness is generally considered as the basic characteristic of meat related to mouth feeling quality. Therefore, tenderization, a common technology in meat processing, was developed to break down the collagens to obtain better palatability. It has been reported that low temperature heating was effective to improve the tenderness of beef and pork [5]. The improved tenderness of meat heated at low temperature may be attributed partially to the action of proteolytic enzymes, causing weakening of fibrils, and solubilization of collagen [6]. Similarly, our previous studies also demonstrated that the hardness and chewiness of AJBW decreased after the low temperature heating time was extended and developed a method using nuclear magnetic resonance and magnetic resonance imaging to assess the proton changes of AJBW [7,8]. Tenderization process mainly involves (1) hydrolysis of muscle tissues by the endogenous proteases, like cathepsins, calpains, and caspases [9][10][11], and (2) assessing physicochemical properties of the muscle protein [12,13]. Previous reports have demonstrated that the release and activation of endogenous enzymes could be induced by apoptosis during meat tenderization [14,15]. However, there are no reports on the potential biological processes of the activation of endogenous proteases in A. japonicus during tenderization process.
The transcriptome is a set of all RNA transcripts and its modifications can exert its effect of protein translation on the phenotype of the organism [16]. Therefore, transcriptome analysis is essential for elucidating the underlying molecular constituents of cells and tissues in various biological processes. In this context, RNA-Seq technology has been widely applied in histological analysis [17], immunology [18], physiology [19], embryonic development, and gene markers [20] in A. japonicus. In addition, it has been used to investigate the effect of heat stress on marine animals such as rainbow trout [21], abalone [22], coral [23], and zooplankton [24]. These published articles found that heat shock could induce oxidative stress, immune response, and protein degradation, thus explaining the mechanism of the heat tolerance.
In sea cucumber, endogenous proteases, like cysteine proteinases [25] and matrix metalloproteinases [26], have been found to degrade the protein in the body wall, while their gene expressions were regulated by some inflammatory cytokines such as tumor necrosis factor (TNF) and interleukin-1 (IL-1) [27]. In this study, the transcriptome of AJBW during tenderization was determined using Illumina Hiseq 4000 and comparative analysis was performed to understand (1) the type of endogenous protease involved and (2) the physiological and biochemical processes that were induced at the mRNA levels of AJBW during tenderization. These new findings would help to understand the tenderization mechanism of AJBW and provide a theoretical basis for the determination of suitable thermal processing conditions.

Transcriptome Sequencing and Quality Control
In this study, nine cDNA libraries were constructed from AJBW following the treatment at different times and sequenced to obtain a comprehensive analysis of tenderization progression. Totals of 48,381,506-58,609,160 raw reads were generated in each library, and the clean reads were 47,514,898-57,471,652 reads with valid ratios of over 97%. Q20 > 99% and Q30 > 91% in each library suggested that the quality of the data was superior (Table 1). In addition, the evenly distributed reads in every position of the genes indicated that the randomness of the breaking was fair. However, we found that the sample correlations were more than 0.6 except in the T_1h_1 samples (the least one was 0.394, Supplementary Figure S1). Thus, the cDNA library of T_1h_1 was excluded to ensure the results were consistent and practical for further analysis of differentially-expressed genes (DEGs).

Reads Mapping to the Reference Genome Dataset
All the clean reads were mapped to the reference genes expressed in the A. japonicus genome [28]. We found that 70.44-82.03% of the reads were perfectly matched to the reference genome in each library. In unique mapped reads, 39.54-44.52% (18,785,024,127) and, in multimapped reads, 28.73-43.71% (16,510,860,562) was matched (Table 2). Gene expression of each library showed a normal distribution (Supplementary Figure S2), indicating consistency in the gene expression of three biological replicates. Meanwhile, most reads matched with exon regions from each library (Supplementary Figure S3).  In total, 52,878 genes were detected in three stages (CON, T_1h, and T_3h) in the AJBW through RNA-seq analysis. Among them, 57.40% (30,350 genes) were mapped to the reference genes with homologous sequences in at least one of the databases (Supplementary Table S1). Moreover, 5415 (17.84%) showed significant matches to Gene Ontology (GO); 3301 (10.88%) to Kyoto Encyclopedia of Genes and Genomes (KEGG); 19,149 (63.09%) to NCBI nonredundant protein sequences (Nr); 12,834 (42.29%) to Swissprot; and 14,875 (49.01%) to Protein family (Pfam). In short, approximately 69.57% (21,115 genes) of the total annotated genes from the AJBW were noted by the main five common databases.

Analysis of DEGs
Compared to the fresh AJBW (CON), 721 and 806 significant DEGs (|log2Ratio| ≥ 1 and p ≤ 0.05) were identified in T_1h and T_3h, respectively ( Figure 1a). Out of these DEGs, 345 upregulated genes and 376 downregulated genes were identified from the comparison of T_1h vs. CON; while 368 upregulated and 438 downregulated were screened from the comparison of T_3h vs. CON. In addition, a strategy was developed to focus on key tenderization genes by comparing the gene expression between T_3h and T_1h. Consequently, only 81 upregulated genes and 94 downregulated genes were screened from the comparison of T_3h vs. T_1h. At the two time points studied, after treatment, the total downregulated genes (750) was slightly more than upregulated genes (631) (Figure 1b,c). There were 146 DEGs consistently upregulated or downregulated in T_1h and T_3h. These upregulated or downregulated DEGs with annotation are presented in Supplementary Table S2. The consistent DEGs contained some important extracellular matrix (ECM)-associated genes (7), such as alpha-2 collagen, 72 kDa type IV collagenase, and matrix metalloproteinase 16 precursor; five cytoskeletal genes, such as myosin heavy chain, troponin I, and titin; and six immune-related genes, such as lipo-polysaccharide (LPS)-induced TNF-alpha factor, fibrinogen-like protein A, complement component 3 C3 (C3), and complement factor B (Bf).
Based on the primary results, the main DEGs that might exhibit important functions when AJBW was tenderized are shown in Table 3. (The total 94 genes with specific functions are shown  in Supplementary Table S3.) The specific genes associated with heat treatment were classified into five groups.

Validation of the Results by RT-qPCR
To confirm the accuracy of the RNA-Seq transcriptome data, 10 consistently upregulated genes were selected for RT-qPCR analysis at three stages of tenderization (CON, T_1h, and T_3h). The results were consistent with the findings achieved by RNA-Seq analysis (Table 4). Single peak in the melting curve was detected, suggesting that all PCR products were specifically amplified. The change scales of the 10 gene expressions were also similar in RT-qPCR analysis when compared with those of RNA-Seq, indicating that the RNA-Seq successfully identified the DEGs.

Gene Ontology Analysis of DEGs
The GO functional analysis, including "cellular component", "molecular function", and "biological process" involved in tenderization, was evaluated. GO terms were assigned to 328 and 487 DEGs for T_1h and T_3h, respectively (Supplementary Table S4). The top 50 terms of each GO analysis are presented in Figure 2. For T_1h vs CON, 58 DEGs could be annotated into 150 biological process, 71 cellular components, and 107 molecular functions. Similarly, between the libraries of T_3h and CON, 105 DEGs could be annotated into 247 biological process, 103 cellular components, and 137 molecular functions. With respect to the two comparisons, the terms "oxidation-reduction process" and "proteolysis" both consisted more enriched DEGs in the "biological process" category, and the terms "cytoplasm" and "extracellular space" consisted more enriched DEGs in the "cellular components" category, while the term "ATP binding" comprised the most enriched DEGs in the "molecular function" category.

Pathway Enrichment Analysis of DEGs
To further understand the biochemical pathways involved in tenderization, enrichment analysis of DEGs was performed as demonstrated in Supplementary Table S5. The top 20 KEGG enrichment results are displayed in Figure 3. The DEGs from the two comparisons were mapped onto 54 and 69 pathways in the KEGG database, respectively, in which three and thirteen pathways were found to be enriched significantly (p < 0.05). In T_1h vs. CON, eleven enriched pathways were related to metabolism and three enriched pathways involving tight junction, focal adhesion, and regulation of actin cytoskeleton were related to cellular processes. Moreover, two immune-related pathways, including the hypoxia inducible factor-1 (HIF-1) signaling pathway and the Ras-proximate-1 (Rap1) signaling pathway, were also involved. In the comparison of T_3h vs CON, nine enriched pathways were related to metabolism, mainly involving nitrogen, glutathione, arginine, and proline metabolism. In addition, the most representative pathways incorporated notch the signaling pathway, antigen processing, and presentation, as well as Fc gamma R-mediated phagocytosis related to immune reaction.

Discussion
In this study, the first transcript analysis of mRNA expression levels was performed using cDNA libraries from AJBW during tenderization process using Illumina Hiseq 4000. Nearly five million valid reads sequenced in each library reached saturation as previously reported [29]. Through the analysis of the DEGs, we acquired a broad understanding of genes involved in the AJBW during tenderization. We found 631 DEGs in T_1h vs. CON and 750 DEGs in T_3h vs. CON. Thousands of variations in the transcription of genes occurred during tenderization. A number of genes related to oxidative stress response, immune response, apoptosis process, cytoskeleton, and ECM varied in terms of their expression levels during tenderization (Table 3). These results will be beneficial in future investigations of the molecular mechanism associated with tenderization in AJBW. Moreover, based on GO (functions) and KEGG (pathways) analyses, the function of DEGs during A. japonicus tenderization was enriched. The analyses showed that multiple GO terms and KEGG pathways were involved in proteolysis, oxidation-reduction process, metabolism, and immune reaction.

Genes of Endogenous Protease
In previous studies, cathepsin L and collagenase have been purified and characterized from the body wall of sea cucumber [30,31]. In addition, the gelatinolytic metalloproteinase (GMP) and cysteine protease were found to degrade the collagen protein [25,32]. In this study, we found that the genes of 72 kDa type IV collagenase and matrix metalloproteinase 16 precursor were significantly upregulated in T_1h and T_3h (Table 3). The 72 kDa type IV collagenase, also known as 72 kDa gelatinase or matrix metalloproteinase-2 (MMP-2), is an ubiquitous metalloproteinase that could degrade collagens (I, IV, V, VII, X, XI, and XIV), gelatin, elastin, fibronectin, aggrecan, and other ECM proteins [27]. Matrix metalloproteinase-16 (MMP-16), also called MT3-MMP, can degrade various components of the ECM, such as collagen type III, gelatin, casein, and fibronectin. The AJBW dermis is a typical catch connective tissue (or called mutable collagenous tissue) that contains a large amount of ECM consisting mainly of collagen fibrils, proteoglycans and microfibrils [33]. Approximately 70% of the total body wall protein was insoluble collagen fibers and the collagen protein belonged to collagen type I formed by (α1) 2 α2 [34]. The noncollagenous protein from the body wall of sea cucumber contains glycoprotein (400 kDa) as the main component exhibiting structural similarity to fibronectins from other vertebrate and invertebrate animals [35]. With respect to the results obtained, we inferred that MMP-2 and MMP-16 were the main endogenous proteases that facilitated the tenderization process by degrading the ECM of the AJBW. Moreover, the processes of oxidative stress, immune response, apoptosis and reorganization of cytoskeleton and ECM were involved during tenderization in AJBW, which might be attributed to activation of MMPs. This is inconsistent with previous studies that also concluded that MMP gene expression may be modulated in both physiological and pathophysiological events [36].

Genes Associated with Oxidative Stress Response
Numerous studies on marine animals, such as coral [37], pacific oyster [38], and marine snail [39], have reported that exposure to thermal stress can induce reactive oxygen species (ROS) accumulation following antioxidant defense mechanism. In our study, oxidant response genes, such as glutathione S-transferase (GST) and glutathione peroxidase (GPx), showed significant upregulated patterns (Table 3), and superoxide dismutase, catalase, and microsomal GST decreased in AJBW during the heating period of 3 h. Similarly, we previously detected two times increment in signal intensity of ROS-derived radicals in T_1h vs CON [40]. Previous studies reported that the transcripts and activities of the antioxidant enzymes could be altered differently by temperature gradients combined with exposure periods [41][42][43]. Moreover, oxidative stress could induce apoptosis [44], autophagy [45], and the immune response [46] by stimulating various types of cascade reactions, besides regulating collagen synthesis and MMPs activities [47]. Our results indicated that the response to oxidative stress induced by heat stress was essential in AJBW during tenderization.

Genes Associated with Immune Response
Innate immune system is the primary immune system for invertebrates due to the lack of acquired or adaptive immune system to defense against the pathogens [48]. As mentioned before, MMPs could be regulated by some cytokines (IL-1β) and transforming growth factor-β1 (TGF-β1) [27], in addition to some specific MMPs that can control chemokine activity [49]. This suggested that immune response could have been stimulated in AJBW during tenderization.
The superfamily of fibrinogen-related proteins (FREPs), such as fibrinogen-like protein, tenascin, ficolin, angiopoietin, and hepassocin, have been reported to play vital roles in innate immune responses [50] and regeneration [51] in invertebrates. Previous studies have also demonstrated that FREP A is widely distributed in the body wall, intestines, longitudinal muscles and respiratory channel of A. japonicus [52]. In our study, we found that the genes of FREP A and fibrinogen C domain-containing protein 1 were upregulated, and the genes of ficolin-2 and tenascin-N were downregulated in AJBW during tenderization. FREPs are supposed to exhibit 'antibody-like' properties in biological mechanisms of recognition and binding of invading pathogens in marine animals like shrimps [53], mussels [54], and marine snails [55]. The roles of FREPs in AJBW during tenderization remain unclear and further investigations are required to reveal their function.
In invertebrates, the complement system is an innate immune response that attacks the surfaces of foreign cells [56]. It consists of three complement pathways such as classical, alternative and lectin pathways. C3 genes have been involved in the three pathways, and Bf genes, that can activate C3, have been involved in the alternative pathway [57]. In our results, the genes of C3, Bf, and Bf-2 were upregulated in T_1h vs. CON and T_3h vs CON. Some previous studies also reported C3-2 could be triggered in A. japonicus by LPS challenge [58] or Vibrio splendidus [59]. However, we could not identify this gene in our transcript libraries. Similarly, an important component of the lectin complement pathway-mannose-binding lectin (MBL)-was also not identified. However, lactose-binding lectins (LBL) and mannan-binding lectin serine protease (MASP) were found to be upregulated in tenderized AJBW (Table 3). It is not known that whether the LBL (just like MBL or ficolin) would combine with MASP to activate C3, thus, further investigation is required. Although it is possible that these complement pathways could be stimulated in AJBW during tenderization that would trigger the terminal pathway to generate membrane attack complex (MAC) affecting the permeabilities of biofilms. Besides, it has been demonstrated that MAC can both induce apoptosis leading to tissue damage [60] and clear of apoptotic cells [56], which suggested that the complement system may be related to apoptosis in AJBW during tenderization.
In A. japonicus, heat shock proteins (HSPs) as a bioindicator of thermal stress, could act as innate immune agents and stimulate the immune system [61]. They are produced by the cells in response to environmental stress factors such as heat shock [62], osmotic stress [63], and ultraviolet light [64]. In this study, four encode HSP70 mRNA were significantly upregulated during the tenderization process. This is in consistence with other reports on A. japonicus [65] and chocolate chip cucumber, Isostichopus badionotus [43], in which respective genes were upregulated due to heat stress. HSP70 had been found to be related to apoptosis [66] and protein degradation [67]. Thus, it can be concluded that the thermal treatment could enhance the expression of HSP70 and was critical to A. japonicus during tenderization.

Genes Associated with Apoptosis Process
Cell apoptosis could be induced by thermal treatment in marine animals [37]. Similarly, MMPs could also affect apoptosis [68]. In our study, based on the GO annotation, two genes encoding apoptosis-related proteins (poly-U-binding factor 60 kDa (PUF60) and TFIIH basal transcription factor complex helicase XPB) were detected and upregulated during tenderization. The first echinoderm PUF60 was identified in gonad, coelomocytes, intestine, respiratory tree, and body wall of sea cucumber Stichopus monotuberculatus. Moreover, the overexpressed PUF60 could induce apoptosis [69]. In addition, two antiapoptosis-related genes (Bax inhibitor-1 and Src family kinase 5) were also identified with downregulated expression during tenderization. Bax belongs to the Bcl-2 family that governs the mitochondrial outer membrane permeabilization and performs proapoptotic function. Bax inhibitor-1 is located in endoplasmic reticulum (ER) membranes and protects the cells from ER stress-induced apoptosis [70]. Src family kinase (belonging to family of non-receptor tyrosine kinases) could inhibit cytokines activation and control the cell apoptosis process [71]. The downregulation of these two genes (Bax inhibitor-1 and Src family kinase 5) further indicated that the original antiapoptosis mechanisms were damaged and irreversible apoptosis occurred in the cells of tenderized AJBW. These results were supported by our previous findings about the DNA damage and enhancement of caspase-3 activity in AJBW during tenderization [40].

Genes Associated with Reorganization of Cytoskeleton and ECM
Thermal treatment can affect the ultrastructure of the cell, which has been found in catfish [72], scallops [73], and corals [74]. Genes encoding cytoskeleton-associated proteins, such as myosin, actin, titin, troponin, tropomodulin, and gelsolin, were upregulated during tenderization (Table 3). Cell volume could be adjusted through manipulation of transporters and cytoskeletal reorganization, which may affect the phagocytic activity [58]. Three genes encoding myosin light chain kinase smooth muscle were upregulated during tenderization, which could have phosphorylated the myosin light chain causing the contractile activity of the smooth muscle [75]. In addition, it has been proved that actin can protect the cells from oxidative stress [76], indicating that the upregulation of actin may build the defense against oxidative damage. The exposure to thermal treatment also induced expression of genes involved in ECM synthesis and remodeling including α collagen gene and MMP-2, a MMP-16 precursor. Besides, integrin β (an important transmembrane receptor in immune systems) was also upregulated ( Table 3). The integrin signaling can regulate ECM remodeling and control subsequent cell behavior and tissue organization [77]. Thus, it can be suggested that ECM reorganization was induced in AJBW during tenderization.

Animals Materials
Live A. japonicus (150 ± 20 g, body weight) were procured from a local aquatic market in Dalian, China. The viscera and heads were removed and each body wall was cut into large pieces and transferred onto the petri plates with three pieces per plate before incubating at 37 • C for 1 h (T_1h) and 3 h (T_3h). This temperature is suitable to activate endogenous enzymes activities. All the body wall chunks, including the incubated and fresh as control samples (CON), were frozen immediately with liquid nitrogen and then stored at −80 • C until RNA extraction.

mRNA Library Construction and Sequencing
Total RNA was extracted from nine samples (three biological replicates each of CON, T_1h and T_3h) using miRNeasy minikit (Qiagen, Duesseldorf, Germany) following the manufacturer's instructions. The quantity and purity of the total RNA were determined using Bioanalyzer 2100 and RNA 6000 Nano Lab Chip Kit (Agilent, Santa Clara, CA, USA) with RIN > 7. Approximately 10 µg of total RNA was used for transcriptome cDNA library construction using mRNA Seq sample preparation kit (Illumina, San Diego, CA, USA). The average insert size for the paired-end libraries was 300 bp (±50 bp). Then we performed the paired-end sequencing on an Illumina Hiseq 4000 (LC Sciences, Houston, TX, USA) following the manufacturer's recommended protocol and generated a total of million paired-end reads of 150 bp length. Prior to assembly, the low quality reads (reads containing sequencing adaptors; reads containing sequencing primer; and nucleotides with q quality score < 20) were removed. The raw sequence data have been submitted to the NCBI Short Read Archive [78] under accession code GSE123077.

RNA-seq Reads Mapping and DEG Testing
We aligned the reads of CON, T_1h and T_3h with A. japonicus reference genome [28] using HISAT package that initially removes a portion of the reads based on quality information accompanying each read and then maps the reads to the reference genome. The mapped reads were assembled using StringTie and merged to reconstruct a comprehensive transcriptome using perl scripts. Later, the expression levels of all transcripts were estimated using StringTie and Ballgown gene expression levels were calculated according to FPKM (Fragments Per Kilobase of exon model per Million mapped reads). The DEGs were selected with log2 (fold change) >1 or log2 (fold change) < −1 with statistical significance (p < 0.05) using R package-Ballgown. The genes were searched using Nr, SwissProt protein database, and Pfam database with E value ≤ 10 −5 . To further identify the functional classifications, all the genes were mapped to the GO database [79] and KEGG resource [80] with a corrected p-value ≤ 0.05 and a Q ≤ 0.05, respectively.

RT-qPCR Validation
To validate the RNA-seq results, 10 DEGs were employed on the qPCR software (version 3.0, qTOWER 2.2, Janalytik Jena, Jena, Germany). Total RNA was isolated from three individuals in each stage (CON, T_1h and T_3h) and was reverse-transcribed into cDNA templates using the PrimeScript™ RT reagent Kit with gDNA Eraser (TaKaRa, Otsu, Japan) according to the manufacturer's instructions. The reaction condition consists of two steps: 37 • C for 15 min and then 85 • C for 5 s. The primer pairs for the selected DEGs are listed in Supplementary Table S6. Amplification of cytochrome b (Cyt b) gene was selected as the reference gene [61]. Each sample run was carried out in triplicate, along with Cyt b. The RT-qPCR amplification was performed in a volume of 20 µL containing 10 µL of 2× TB Green Premix Ex Taq™ II, 2 µL of cDNA template, 1 µL each of forward and reverse primers (10 µM), and 6 µL deionized water using TB Green TM Premix Ex Taq™ II Kit (Tli RNase H Plus, TaKaRa, Otsu, Japan) according to the manufacturer's instructions with slight modification. The PCR conditions were 95 • C for 30 s, 40 cycles (heating at 95 • C for 5 s and annealing at 60 • C for 30 s), 95 • C for 30 s, 60 • C for 30 s, and 95 • C for 15 s. The 2 −∆∆CT method was used to analyze the comparative expression levels. Statistical analysis was performed using SPSS 16.0 software (IBM Corp., Armonk, NY, USA) with the test conducted as a two-tailed Independent Samples t-test and a significance level of p < 0.05.

Conclusions
In this study, we used RNA-seq to identify the various DEGs in the tenderized AJBW and found that the genes of MMP-2 and MMP-16 were significantly upregulated which could efficiently explain the changes in tenderness in AJBW during 37 • C treatment. Moreover, the results suggested that the tenderization is a complex process that involved oxidative stress, immune response, autolysis, and reorganization of the cytoskeleton and ECM. These specific relationships of the regulatory mechanisms and endogenous protease activation warrant further investigations that can help to further improve the thermal processes required to maintain the high quality of A. japonicus.

Conflicts of Interest:
The authors declare no conflicts of interest.