Comparative transcriptome analysis provides comprehensive insights into the heat stress response of Marsupenaeus japonicus

Aquaculture 502 (2019) 338–346

Contents lists available at ScienceDirect

Aquaculture journal homepage: www.elsevier.com/locate/aquaculture

Comparative transcriptome analysis provides comprehensive insights into the heat stress response of Marsupenaeus japonicus

T

⁎

Jinbin Zhengb, Jiawen Caoc, Yong Maoa,b, , Yongquan Sub, Jun Wangb a

State Key Laboratory of Marine Environmental Science, Xiamen University, Xiamen 361102, China College of Ocean and Earth Sciences, Xiamen University, Xiamen 361102, China c College of the Environment and Ecology, Xiamen University, Xiamen 361102, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Marsupenaeus japonicus Heat stress Comparative transcriptome Immune molecules Heat shock protein Antioxidant enzyme Metabolism

The kuruma shrimp Marsupenaeus japonicus is an economically important crustacean that is widely distributed throughout the Indo-West Pacific. However, M. japonicus exhibits relatively poor thermotolerance properties and frequently suffers from mass mortality during the hot summer months in aquaculture systems. Herein, we performed a comparative transcriptome analysis of M. japonicus in response to heat stress, focusing on immune responses, heat shock proteins, antioxidant systems and metabolism alterations under heat stress. The results revealed that M. japonicus generated adaptive responses to maintain physiological homeostasis under heat stress via the transcriptional up-regulation of heat shock proteins and antioxidant enzymes, enhancement of fatty acid metabolism. However, acute and prolonged exposure to heat stress resulted in the continually decreased expression of immune molecules and led to the repression of glycometabolism, the tricarboxylic acid cycle and the respiratory chain. This study provides a systematic overview of the M. japonicus heat stress response and elucidates potential molecular mechanisms underlying this phenomenon and summer mortality syndrome in this species, which is beneficial for proposing effective strategies to address summer mortality syndrome and promote kuruma shrimp production.

1. Introduction Water temperature is one of the most crucial ecological factors in aquatic ecosystems that modulates many aspects of aquatic organisms, including their metabolism, growth, development, reproduction and immune response (Azaza et al., 2008; Bermudes-Lizárraga et al., 2017; Cheng et al., 2004; Munday et al., 2008; Pankhurst and Munday, 2011; Wang and Chen, 2006; Xu et al., 2017). In the context of global warming (Hoegh-Guldberg and Bruno, 2010; Poloczanska et al., 2013), intense temperature fluctuations and increased water temperature extremes have posed serious threats to aquatic organisms in aquaculture systems and caused tremendous losses for aquaculture industries (Chen et al., 2016; Liang et al., 2014; Shao et al., 2015; Lang et al., 2009; Rebl et al., 2013). The kuruma shrimp Marsupenaeus japonicus, which is widely distributed throughout the Indo-West Pacific (Tsoi et al., 2007; Tsoi et al., 2005), is one of the main shrimp species with a high aquaculture interest and is mainly cultured in China (FAO, 2018). However, its farming industry has been limited by mass mortality during warm summer months, which causes massive production loss and deeply

⁎

obstructs the development of the M. japonicus aquaculture industry. Consequently, the annual production of M. japonicus is only approximately 50,000 tons, which is far below that of other penaeid shrimp such as Litopenaeus vannamei and Penaeus monodon (FAO, 2018). To date, no effective method exists for preventing M. japonicus mass mortality during hot summer months. Better understanding the molecular mechanisms underlying the heat stress response of M. japonicus will provide a theoretical foundation for elucidating effective strategies to address the summer mortality syndrome and promote kuruma shrimp production. Nevertheless, few studies have focused on illustrating the heat stress response of M. japonicus, and the molecular mechanisms underlying this phenomenon are not fully elucidated. With the rapid development of next-generation sequencing (NGS) technologies, transcriptome analysis has become an attractive, lowcost, highly accurate method for studying non-model organisms that lack genomic information. Recently, an increasing number of transcriptome analyses have been performed in a variety of commercially important aquatic organisms using RNA-Seq technology (Lim et al., 2016; Lv et al., 2014; Nguyen et al., 2016; Shiel et al., 2015; Xing et al., 2018; Xiong et al., 2014; Xue et al., 2013). However, until recently,

Corresponding author at: College of Ocean and Earth Sciences, Xiamen University, South Xiangan Road, 361102 Xiamen, China. E-mail address: [email protected] (Y. Mao).

https://doi.org/10.1016/j.aquaculture.2018.11.023 Received 24 March 2018; Received in revised form 1 November 2018; Accepted 10 November 2018 Available online 20 December 2018 0044-8486/ © 2018 Elsevier B.V. All rights reserved.

Aquaculture 502 (2019) 338–346

J. Zheng et al.

2.4. Bioinformatics analysis of the transcriptome data

relatively few transcriptome analyses have been conducted in M. japonicus, and no information on the transcriptome analysis of M. japonicus in response to heat stress is available. In this context, we performed a comparative transcriptome analysis of M. japonicus in response to heat stress, focusing on the immune response, heat shock proteins (HSPs), antioxidant systems and metabolism alterations under heat stress. This study provides an insight into the gene expression profiles of M. japonicus in response to heat stress and the potential molecular mechanisms underlying the summer mortality syndrome in M. japonicus, which is beneficial for constructing optimal strategies to improve its cultivability during the summer. Furthermore, the tremendous amount of data obtained from RNA-Seq provide valuable information and candidate genes for the genetic improvement of M. japonicus.

2.4.1. Quality control and de novo assembly Raw data were firstly processed through in-house perl scripts. In this step, reads shorter than 100 bp, reads without adaptors and low quality reads were removed. Then, adaptors and barcodes were removed from the reads, and reads shorter than 25 bp after trimming and contaminant removal were discarded. The screened high-quality sequences (clean reads) were de novo assembled using Trinity version trinityrnaseq_r20131110 as described for de novo transcriptome assembly without a reference genome (Grabherr et al., 2011). The clean reads were mapped to the assembled transcriptome using Bowtie2. 2.4.2. Identification of differentially expressed genes and clustering analysis Differentially expressed genes (DEGs) were obtained by the comparison of fragments per kilobase of transcript per million mapped reads (FPKM) values between the unstressed and heat stressed samples (Mortazavi et al., 2008). Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by edgeR program package through one scaling normalized factor. Differential expression analysis between two samples was performed using the DEGseq R package, and the p-value was adjusted using the q-value. Transcripts with fold changes (FCs) ≥ 2 and q-values < 0.005 were considered DEGs.

2. Materials and methods 2.1. Animals All M. japonicus prawns (4.76 ± 1.22 g in body weight) were obtained from an aquaculture farm in Dongshang (Zhangzhou, Fujian, China) and acclimated in environmentally controlled flat-bottomed rectangular tanks (70 cm × 50 cm × 40 cm) with aerated seawater at 28 °C and a salinity of 28‰ prior to experimentation. The seawater was renewed daily, and the prawns were fed twice daily with commercial pellets (Fuxing, China).

2.4.3. Functional annotation For functional annotation, all transcripts were searched against the databases NCBI non-redundant (nr) protein (ftp://ftp.ncbi.nih.gov/ blast/db/), Pfam (http://pfam.xfam.org/), Swiss-Prot (http://www. uniprot.org/), Gene Ontology (GO) (http://www.geneontology.org/), Clusters of Orthologous Groups of proteins (COG) (http://www.ncbi. nlm.nih.gov/COG/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/) using BLASTX programs with an E-value threshold of 10−5.

2.2. Heat stress and samples The heat shock experiment was conducted as previously described with slight modifications (Rungrassamee et al., 2010). A group of 50 shrimp was transferred from the acclimation tank (28 °C) to a tank containing seawater preheated and maintained at 32 °C. The hepatopancreases and gills of five individuals were collected from unstressed shrimp (control) and shrimp subjected to heat stress for 3 h (3 h HS) and 96 h (96 h HS). All samples were immediately frozen at −80 °C until the extraction of total RNA for analysis.

2.5. Validation of mRNA-Seq data by quantitative real-time PCR Several interested genes, including HSPs, heat shock transcription factor 1 (HSF1), immune molecules, antioxidant enzymes and metabolism-related genes, were selected for validation by quantitative realtime PCR (qRT-PCR), using the same samples used for transcriptome sequencing. QRT-PCR was performed on the Applied Biosystems QuantStudio 6 Flex Real-time PCR System (Applied Biosystems, USA) with SYBR® Premix DimerEraser™ (TaKaRa, Japan) according to the manufacturer's instructions. Five biological and three technical replicates were performed, and elongation factor 1- α (EF1-α) served as the reference for internal standardization. The specific primers for qRTPCR are listed in Table 1. The PCR cycling conditions were as follows: 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s, 60 °C for 30 s and 72 °C for 30 s. Melting curve analysis was performed to evaluate the specificity of amplification. The slopes and regression curves of the standard curves were calculated to determine the PCR efficiency (E) with 7 tenfold template dilutions. The E-value was calculated according to the eq. E = 10(−1/slope) − 1, and the 2−ΔΔCt method was used to estimate the relative mRNA expression levels as previously described (Livak and Schmittgen, 2001). Statistical significance was determined using one-way analysis of variance (ANOVA) followed by multiple comparison testing with the least significant distance (LSD) t-test using SPSS 17.0 software. P values < .05 were considered to represent statistically significant differences.

2.3. Library construction and sequencing 2.3.1. Total RNA extraction and qualification Total RNA was extracted using RNAiso Plus (TaKaRa, Dalian China) according to the manufacturer's instructions. RNA degradation and contamination were assessed on 1% agarose gels. The NanoPhotometer® spectrophotometer (IMPLEN, CA, USA), Qubit®2.0 Fluorometer (Life Technologies, CA, USA) and Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) were used to detect the purity, concentration and integrity of the RNA samples. The cut-off thresholds for the RNA integrity number (RIN) values and concentrations were 7.0 and 300 ng/μl, respectively.

2.3.2. Library preparation and sequencing Six libraries were constructed by pooling equal amounts of total RNA (300 ng per sample) extracted from the hepatopancreases and gills of five biological samples at different heat stress stages, unstressed-Hp, 3 h HS-Hp, 96 h HS-Hp, unstressed-G, 3 h HS-G and 96 h HS-G. Transcriptome sequencing libraries were generated using the NEBNext®Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) according to the manufacturer's recommendation and quantitative realtime PCR (qRT-PCR) was performed for library quality control. The qualified libraries were sequenced on an Illumina HiSeq X-ten platform and 150 bp paired-end reads were generated.

3. Results 3.1. Sequencing and de novo assembly of the transcriptome After the removal of low-quality reads, 24,017,119–28,936,095 339

Aquaculture 502 (2019) 338–346

J. Zheng et al.

Table 1 PCR primers for the validation of RNA-Seq data by qRT-PCR.

Table 2 Sequencing and assembly statistics of the transcriptome data.

Name

Primer sequence (5′-3′)

Target

Primer efficiency

hsp10-RT-F hsp10-RT-R hsp60-RT-F hsp60-RT-R hsp70-RT-F hsp70-RT-R hsp90-RT-F hsp90-RT-R HSF1-RT-F HSF1-RT-R ALF6-RT-F ALF6-RT-R ALF3-RT-F ALF3-RT-R PPAE-RT-F PPAE-RT-R SPI-RT-F SPI-RT-R CytMnSODRT-F CytMnSODRT-R MitMnSODRT-F MitMnODRT-R CZSOD2-RTF CZSOD2-RTR CZSOD4-RTF CZSOD4-RTR Gpx-RT-F Gpx-RT-R Gpx3-RT-F Gpx3-RT-R CAT-RT-F CAT-RT-R CS-RT-F CS-RT-R IDH-RT-F IDH-RT-R SDH-RT-F SDH-RT-R F-ATP-α-RT-F F-ATP-α-RTR F-ATP-β-RT-F F- ATP-β-RTR PEPCK-RT-F PEPCK-RT-R GP-RT-F GP-RT-R SCD-RT-F SCD-RT-R FAS-RT-F FAS-RT-R LACS-RT-F LACS-RT-R EF1-α-F

TGAGTTGAAAATGGCTGGTGC GTCCTGGCTCCCTCTCCTACA GGATGTCGCAGAACGATGTCTAG CGAGGTACTTCGATGGAATGATTAA TCATCAACGAGAGCACGAA TCCGCAGTTTCCTTCATCT CTCTATGATCGGTCAGTTCGG CGGCGTTCCTCAAGGTATTCT ATGGCAAGTTTTGTGAGGCAG GAAAGTGTGGGTGGGCGAA ATACTTCCTCGGCACTGTCAC TATGGTTTTGCTTCCTCCTCA TGGAGGAACGAGAAGACCGA ACTGGCTGCGTGTGCTGG GACCCGTTTCGGAGAGTG CGTTGGGAGCCAGTGTGA AACTTCGCTCAGAACCCG ACAGCCAAGGCAGTCGTC CAGTGTCGGAGTGAAAGGC

qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR

98.48% 98.48% 99.45% 99.45% 104.91% 104.91% 99.18% 99.18% 99.78% 99.78% 99.32% 99.32% 91.91% 91.91% 94.28% 94.28% 101.92% 101.92% 93.588%

CTGAAGGTAGTAGGCGTGC

qRT-PCR

93.588%

TCAGGCTGGGGATGGTTGG

qRT-PCR

94.972%

CATAATCAGGTCGCACATTCTTG

qRT-PCR

94.972%

GAGCCGACACTTCCGACATC

qRT-PCR

91.972%

TGCGTGAGGGTGAGCGTTC

qRT-PCR

91.972%

AGGCTCCATCGTCATCTCACA

qRT-PCR

95.591%

CTTCCATCTGTCACCGCTCG

qRT-PCR

95.591%

GAAACGGCTTTGTTCCTGACTT CGGGTGGTAGAAGAGTTGGC AGAGCCTGGAGCCACTGCC GTTGTCCCCGTTGACCTGC GTGAGATTGGCAAACGAACT GATGAAGAAGATAGGGGTGT TGGGCTGGAAACTTCAGTCACA GAGACATTACCGCCTTCGTGG GCAAGGTCGGAATCAAATGTG GAATGGGCTCCCTGAAAACG TTCCACCCTACGGGCATTTAT ACGGCTGACGACATCTCTTGA GAAACTGATAATGTCGGTGTCGT TTTCCTTTGAGTGGTTTTGACAG

qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR

95.398% 95.398% 93.672% 93.672% 94.89% 94.89% 100.707% 100.707% 100.437% 100.437% 101.025% 101.025% 97.378% 97.378%

TGAATCCAAAGTGAAATACCCGA CGACTCCTGCTCCACCGAAG

qRT-PCR qRT-PCR

101.698% 101.698%

CGGTGGGAACTCTCTCTTGG CCTGTGGGTTTGTGATGCC TGGTATCACTCCTCGTCGTTG GCAAGGGGTTTTAGTTTGGTTAG GACTGACGCAGACCCACACA CACCAGCGGCAGGTAGAACT AACCTACCGCCACATTCCACT CATAGAGCACATCACATACGAGCA GAGTGCTGAGTTGGGACGATG CCCTTTGGAGGACCAGTTGTG GGAACTGGAGGCAGGACC

96.724% 96.724% 97.928% 97.928% 96.850% 96.850% 96.409% 96.409% 102.668% 102.668% 101.18%

EF1-α-R

AGCCACCGTTTGCTTCAT

qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR qRT-PCR Internal control Internal control

Library

Clean reads

Clean reads ratio

Mapped reads

Mapped ratio

Unstressed-Hp 3 h HS-Hp 96 h HS-Hp Unstressed-G 3 h HS-G 96 h HS-G

28,936,095 26,260,553 26,770,939 24,017,119 25,721,806 27,975,940

86.97% 86.84% 86.76% 85.64% 85.83% 87.36%

23,587,503 22,050,373 21,841,815 19,485,132 20,872,322 23,210,704

81.52% 83.97% 81.59% 81.13% 81.15% 82.97%

775 bp. After mapping the clean reads to the reference transcriptome, 19,485,132–23,587,503 mapped reads were identified in the six samples (Table 2). 3.2. DEG identification To identify genes potentially involved in the heat stress response of M. japonicus, DEGs at different heat stress stages were identified by comparing the relative gene expression abundances. A total of 301, 443, 1136 and 1228 DEGs were detected in the hepatopancreas and gill at 3 h HS and 96 h HS, respectively (Table 3). More genes were downregulated than up-regulated at 3 h HS in both the hepatopancreas and gill, while the inverse result was observed at 96 h HS (Table 3). Moreover, the number of DEGs in the gill was much greater than that in the hepatopancreas (Table 3). 3.3. Functional annotation and pathway assignment of the DEGs The DEGs at different stages of heat stress showed similar distributions of GO categories. In the molecular function category, most DEGs were determined to be involved in catalytic activity, structural molecule activity and binding (Fig. 2). Among the DEGs categorized as a biological process, metabolic processes, cellular processes, single-organism processes and responses to stimuli were the most represented (Fig. 2). Additionally, most of the cellular component-related DEGs were associated with cells, organelles and cell parts (Fig. 2). To identify the biological pathways that are activated or depressed during the response of M. japonicus to heat stress, the DEGs were mapped to referential canonical pathways in the KEGG database. A total of 186 KEGG pathways were mapped (Supplementary 1), among which a batch of nutrient metabolism pathways involved in glycometabolism, lipid metabolism and amino acid metabolism were identified in the hepatopancreas (Supplementary 1). Additionally, two crucial energy metabolism- related pathways, the tricarboxylic acid cycle (TCA cycle) and respiratory chain oxidative phosphorylation pathway were enriched in the gill (Supplementary 1). 3.4. Verification of transcriptomic data by qRT-PCR The relative mRNA expression levels of several interested genes were measured by qRT-PCR. Although the exact FCs in the genes at several data points varied between the RNA-Seq and qRT-PCR analyses, the differential expression tendencies detected by the two methods showed good consistency (Fig. 3). 4. Discussion

101.18%

M. japonicus is one of the most economically important shrimp species worldwide (Tsoi et al., 2005, 2007). However, it frequently suffers from mass mortality during warm summer months, which results in enormous economic losses and deeply hinders the development of the M. japonicus aquaculture industry. Herein, a comparative transcriptome analysis was conducted to elucidate the molecular mechanisms underlying the heat stress response and summer mortality syndrome in M. japonicus. Illumina sequencing yielded a total of 95,763

clean reads were generated from the six libraries via transcriptome sequencing (Table 2). De novo assembly was performed using the merged reads of all the samples, generating a total of 95,763 transcripts as the M. japonicus reference transcriptome. The N50 length of the assembled transcripts was 1369 bp, and the average contig length was 340

Aquaculture 502 (2019) 338–346

J. Zheng et al.

Table 3 Number of DEGs between libraries. DEG Set

All DEGs

Up-regulated

Down-regulated

Annotated

Unstressed-Hp vs 3 h HS-Hp Unstressed-Hp vs 96 h HS-Hp Unstressed-G vs 3 h HS-G Unstressed-G vs 96 h HS-G

301 443 1136 1228

137 240 189 730

164 203 947 498

160 252 452 542

Fig. 1. Effects of heat stress on the expression of heat shock proteins, antioxidant enzymes and immunerelated genes in the hepatopancreas and gill. The green arrow indicates down-regulation, and the red arrow indicates up-regulation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

valuable information for follow-up studies on the specific biological and metabolic processes involved in the heat stress response of M. japonicus. Shrimp lack an adaptive immune system and rely on their innate cellular and humoral immune responses to combat invading pathogens. The cellular immune reactions include phagocytosis, nodule formation and encapsulation, whereas the humoral responses include the prophenoloxidase (proPO) system, the clotting cascade and the synthesis and release of a wide array of antimicrobial peptides (AMPs) (Li and Xiang, 2013; Tassanakajon et al., 2013). Heat stress has been reported to reduce the immune ability and increase the susceptibility of prawns to pathogens (Cheng et al., 2005; Wang et al., 2006). However, very little is known about the molecular mechanisms by which heat stress influences the shrimp immune system. In this study, the transcription of numerous immune molecules, e.g., crustin, α-2 macroglobulin homolog, ALF, PPAE, prophenoloxidase b, C-type lectin, E-type lectin, peroxinectin and Spätzle, were significantly decreased in either the hepatopancreas or gill of M. japonicus under heat stress (Fig. 1), revealing that heat stress can obviously suppress shrimp immune capacity. Previous studies have revealed that heat stress inhibits the immune capability of the Pacific whiteleg shrimp Litopenaeus vannamei for a short period, as the shrimp establish an immune adaptability after a period of immune regulation (Lu et al., 2007). However, in the present study, the transcriptional levels of a batch of M. japonicus immune

transcripts with an average length of 775 bp and an N50 of 1369 bp. The percentage of clean reads mapped to the reference transcriptome was > 81.13% in the six libraries (Table 2), suggesting a good assembly quality of the transcriptome. Numerous DEGs were identified at different stages of heat stress in either the hepatopancreas or gill of M. japonicus (Table 3), which provided a stepping-stone for further mining thermotolerance-related genes as well as valuable genomic resources for breeding programs aimed at improving the thermoresistance performance of M. japonicus. To understand the possible functions of the DEGs, all DEGs were searched against the databases NCBI nr protein, Pfam, Swiss-Prot, GO, COG and KEGG. However, approximately half of the DEGs were not annotated (Table 3), which may be due to the limited information available on the gene backgrounds of crustaceans and arthropods (Li et al., 2012). Catalytic activity, structural molecule activity, binding, metabolic processes, cellular processes, single-organism processes, responses to stimuli, cells, organelles and cell parts were the highly represented GO terms observed in the DEGs. Enrichment in these GO categories indicated that cellular activities as well as, metabolic and cellular processes in M. japonicus are impacted by heat stress. A total of 186 KEGG pathways were found to be altered under heat stress, implying that a broad spectrum of biological processes is involved in the heat stress response of M. japonicus. These annotations provided 341

Aquaculture 502 (2019) 338–346

J. Zheng et al.

Fig. 2. Distribution of GO second-level functional annotations. (A), GO classifications of DEGs between the unstressed-Hp and 3 h HS-Hp groups. (B), GO classifications of DEGs between the unstressed-Hp and 96 h HS-Hp groups. (C), GO classifications of DEGs between the unstressed-G and 3 h HS-G groups. (D), GO classifications of DEGs between the unstressed-G and 96 h HS-G groups.

the hepatopancreas or gill of M. japonicus under heat stress (Fig. 1). The up-regulation of heat HSPs and antioxidant enzymes is generally considered as an adaptive process by which aquatic organisms cope with heat stress, which is possibly associated with relieving oxidative injury and providing protection for proteins. Heat stress has obvious effects on the metabolic characteristics of M. japonicus. In this study, the transcription levels of several key ratelimiting enzymes in glycometabolism pathways, such as phosphenolpyruvate carboxykinase (PEPCK), α-glucosidase and glycogen phosphorylase (GP), were significantly down-regulated in the hepatopancreas of M. japonicus under heat stress. However, a batch of enzymes involved in fatty acid metabolism, e.g., stearoyl-CoA desaturase (SCD), fatty acid synthase (FAS) and long-chain fatty acid acyl-CoA synthetase (LACS), were up-regulated under heat stress. Above all, these results provide a global picture of a shift in energy usage from glucose to lipids. PEPCK is a key enzyme in the gluconeogenesis pathway, which catalyzes oxalacetate to phosphoenolpyruvate (PEP), generally considered the first committed step in gluconeogenesis, and thus strongly affects intracellular glucose biosynthesis (Hanson and Reshef, 1997).α-glucosidase, which catalyzes the hydrolysis of glycosidic bonds in oligosaccharides or glycoconjugates, is fundamental for the degradation of diet polysaccharides to supply the monosaccharide unit, principally starch (Arellano-Carbajal and Olmos-Soto, 2002; Chiba, 1997). Glycogen phosphorylase is a key enzyme involved in glycogen catabolism that catalyzes the initial and key control step of glycogen degradation (Roach et al., 2011; Zhao et al., 2014). The down-regulated

molecules exhibited a tendencies to decrease as the heat stress time prolonged. The impacts of heat stress on the aquatic organism immune response may be species-specific, and the continuously decreased immune molecule expression levels during the course of heat stress might contribute to the mass mortality of M. japonicus during hot summer months. Previous studies have demonstrated that heat stress can disturb the balance between the production and elimination of reactive oxygen species (ROS) and may further facilitate the formation of ROS (Lushchak, 2011). Excessive ROS induced by heat stress can cause oxidative injury to cellular macromolecules, such as DNA, proteins and lipids (Lesser, 2006). To protect against the toxicity caused by ROS and maintain proper cellular functions, organisms have developed effective and efficient antioxidant defense systems, including superoxide dismutase (SOD), catalase (CAT), glutathione peroxidase (Gpx) and glutathione S-transferase (GST), to eliminate excessive ROS and maintain a ROS balance (Jacob, 1995; Limón-Pacheco and Gonsebatt, 2009). In addition, HSPs, also known as molecular chaperones, help to sustain cell homeostasis by refolding denatured proteins, degrading unstable or misfolded proteins, preventing proteins from aggregation and restoring the function of proteins denatured by heat stress (Feder and Hofmann, 1999; Hendrick and Hartl, 1993; Kregel, 2002). In the present study, the transcription of several HSPs (hsp10, hsp60, hsp70, hsp90) as well as those of antioxidant enzymes (copper/zinc SOD isoform 2 (CZSOD2), CZSOD4, cytosolic MnSOD (CytMnSOD),mitochondrial MnSOD (MitMnSOD), CAT, Gpx and GST) were significantly increased in either 342

Aquaculture 502 (2019) 338–346

J. Zheng et al.

Fig. 3. Comparison of the expression profiles of nine genes determined by RNA-Seq and qRT-PCR in the hepatopancreas (A) and gill (B). Each bars represents the expression fold change in a gene compared to that in the control group. Bars with asterisks indicate values (mean ± SD) that were significantly different at p < .05 (n = 5, ANOVA) from those in the control group.

343

Aquaculture 502 (2019) 338–346

J. Zheng et al.

Fig. 4. Impact of heat stress on the TCA cycle and ETC pathways in the gill. The green arrow indicates down-regulation of the transcription level. Abbreviations of enzymes involved in the TCA cycle are as follows: CS, citrate synthase; ACO, aconitase; IDH, isocitrate dehydrogenase; α-KGDH, α-ketoglutarate dehydrogenase; SCoAS, succinyl-CoA synthetase; SDH, succinate dehydrogenase; FUM, fumarase; and MDH, malate dehydrogenase. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. Proposed scenario involving the heat stress response of M. japonicus. The red arrow indicates activation, and the green, flat arrow indicates repression. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

2004; Ntambi, 1995). Because unsaturated fatty acids are key components of cellular membranes, SCD plays an essential role in the regulation of membrane composition and fluidity (Castro et al., 2011). The increased SCD transcriptional levels are considered an important adaptative mechanism of aquatic organisms in response to temperature

transcription of these key enzymes implies that glycometabolism is potentially suppressed in the hepatopancreas of M. japonicus under heat stress. SCD is a key pace-making enzyme in the biosynthesis of monounsaturated fatty acids that catalyzes the initial step in long-chain unsaturated fatty acids formation (Kim and Ntambi, 1999; Ntambi et al., 344

Aquaculture 502 (2019) 338–346

J. Zheng et al.

fluctuations (Hsieh et al., 2007; Hsieh and Kuo, 2005). FAS is a multifunctional enzyme that catalyzes the synthesis of long-chain fatty acids, which are essential constituents of membrane lipids and important substrates for cellular energy metabolism (Smith et al., 2003; Wakil, 1989). LACS, catalyzes the synthesis of acyl-CoA, and plays important roles in both lipid synthesis and the degradation of fatty acids via β-oxidation (Coleman et al., 2002; Suzuki et al., 1990). The up-regulation of these rate-limiting enzyme genes implies an enhancement of fatty acid metabolism, which is regarded as a potential energy budget strategy in response to heat stress (Thorne et al., 2010). The temperature adaptation of aquatic organisms relies on their glycometabolism and mitochondrial metabolism capacities, which includes the respiratory chain and TCA cycle (Anestis et al., 2007). The TCA cycle and respiratory chain oxidative phosphorylation are central biochemical pathways in cellular energy metabolism (Fernie et al., 2004; Oexle et al., 1999). The TCA cycle, an iconic metabolism pathway in aerobic organisms, is a cyclic flux that is responsible for the total oxidation of acetyl-CoA and provides reduced molecules (NADH and FADH2) for the mitochondrial respiratory chain to fuel ATP formation (Briere et al., 2006; Nunes-Nesi et al., 2013; Sweetlove et al., 2010). The mitochondrial respiratory chain, also referred to as the electron transport chain (ETC), has a pivotal function in supplying the cell with energy by synthesis of ATP via oxidative phosphorylation (Fernie et al., 2004). Generally, aquatic organisms are forced to generate more energy to maintain physiological homeostasis under heat stress (Wang et al., 2006), however, in the present study, the transcriptional levels of several critical enzymes and components in the TCA cycle and ETC were down-regulated in the M. japonicus gill under heat stress (Fig. 4). The transcriptional levels of citrate synthase (CS), isocitrate dehydrogenase (IDH), succinate dehydrogenase (SDH) and ATP synthase were significantly down-regulated to the level that can neither be detected by RNA-Seq nor qRT-PCR at 96 h HS (Fig. 4). These results suggest that heat stress may suppress the energy metabolism of M. japonicus via a transcriptional mechanism. The mismatch between energy supply and demand resulted in a shortage of energy that was available for normal energy-requiring processes as well as for sustaining physiological homeostasis under heat stress, thus potentially leading to mass mortality during hot summer months.

Acknowledgments This research was supported by the Project of China Agriculture Research System (Grant No. CARS-48), the Project of Fujian Provincial Department of Science and Technology (Grant No. 2016NZ0001-4), the Marine Economy Innovation and Development Project of Xiamen (Grant No. 16CZY009SF05). Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.aquaculture.2018.11.023. References Anestis, A., Lazou, A., Pörtner, H.O., Michaelidis, B., 2007. Behavioral, metabolic, and molecular stress responses of marine bivalve Mytilus galloprovincialis during long-term acclimation at increasing ambient temperature. Am. J. Physiol. Regul. Integr. Comp. Physiol. 293, R911–R921. Arellano-Carbajal, F., Olmos-Soto, J., 2002. Thermostable α-1, 4-and α-1, 6-glucosidase enzymes from Bacillus sp. isolated from a marine environment. World J. Microb. Biot. 18, 791–795. Azaza, M.S., Dhraief, M.N., Kraiem, M.M., 2008. Effects of water temperature on growth and sex ratio of juvenile Nile tilapia Oreochromis niloticus (Linnaeus) reared in geothermal waters in southern Tunisia. J. Therm. Biol. 33, 98–105. Bermudes-Lizárraga, J., Nieves-Soto, M., Medina-Jasso, A., Pina-Valdez, P., 2017. Effect of temperature and salinity on larval survival and development of Litopenaeus vannamei. Rev. MVZ Cordoba 22, 5844–5853. Briere, J.J., Favier, J., Gimenez-Roqueplo, A.P., Rustin, P., 2006. Tricarboxylic acid cycle dysfunction as a cause of human diseases and tumor formation. Am. J. Physiol. Cell Phys. 291 (6), C1114–C1120. Castro, L.F.C., Wilson, J.M., Gonçalves, O., Galante-Oliveira, S., Rocha, E., Cunha, I., 2011. The evolutionary history of the stearoyl-CoA desaturase gene family in vertebrates. BMC Evol. Biol. 11, 132. Chen, N., Luo, X., Gu, Y., Han, G., Dong, Y., You, W., Ke, C., 2016. Assessment of the thermal tolerance of abalone based on cardiac performance in Haliotis discus hannai, H. gigantea and their interspecific hybrid. Aquaculture 465, 258–264. Cheng, W., Hsiao, I.S., Hsu, C.H., Chen, J.C., 2004. Change in water temperature on the immune response of Taiwan abalone Haliotis diversicolor supertexta and its susceptibility to Vibrio parahaemolyticus. Fish Shellfish Immunol. 17, 235–243. Cheng, W., Wang, L.U., Chen, J.C., 2005. Effect of water temperature on the immune response of white shrimp Litopenaeus vannamei to Vibrio alginolyticus. Aquaculture 250, 592–601. Chiba, S., 1997. Molecular mechanism in α-glucosidase and glucoamylase. Biosci. Biotechnol. Biochem. 61, 1233–1239. Coleman, R.A., Lewin, T.M., Van Horn, C.G., Gonzalez-Baró, M.R., 2002. Do long-chain acyl-CoA synthetases regulate fatty acid entry into synthetic versus degradative pathways. J. Nutr. 132, 2123–2126. FAO, 2018. Food and Agriculture Organization of the United Nations Global Production Statistics. http://www.fao.org/fishery/statistics/global-aquaculture-production/ query/en, Accessed date: 30 March 2018. Feder, M.E., Hofmann, G.E., 1999. Heat-shock proteins, molecular chaperones, and the stress response: evolutionary and ecological physiology. Annu. Rev. Physiol. 61, 243–282. Fernie, A.R., Carrari, F., Sweetlove, L.J., 2004. Respiratory metabolism: glycolysis, the TCA cycle and mitochondrial electron transport. Curr. Opin. Plant Biol. 7, 254–261. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q.D., Chen, Z.H., Mauceli, E., Hacohen, N., Gnirke, A., Rhind, N., Palma, F.D., Birren, B.W., Nusbaum, C., Lindblad-Toh, K., Friedman, N., Regev, A., 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. Hanson, R.W., Reshef, L., 1997. Regulation of phosphoenolpyruvate carboxykinase (GTP) gene expression. Annu. Rev. Biochem. 66, 581–611. Hendrick, J.P., Hartl, F.U., 1993. Molecular chaperone functions of heat-shock proteins. Annu. Rev. Biochem. 62, 349–384. Hoegh-Guldberg, O., Bruno, J.F., 2010. The impact of climate change on the world's marine ecosystems. Science 328, 1523–1528. Hsieh, S.L., Kuo, C.M., 2005. Stearoyl-CoA desaturase expression and fatty acid composition in milkfish (Chanos chanos) and grass carp (Ctenopharyngodon idella) during cold acclimation. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 141, 95–101. Hsieh, S.L., Hu, C.Y., Hsu, Y.T., Hsieh, T.J., 2007. Influence of dietary lipids on the fatty acid composition and stearoyl-CoA desaturase expression in hybrid tilapia (Oreochromis niloticus × O. aureus) under cold shock. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 147, 438–444. Jacob, R.A., 1995. The integrated antioxidant system. Nutr. Res. 15, 755–766. Kim, Y.C., Ntambi, J.M., 1999. Regulation of stearoyl-CoA desaturase genes: role in cellular metabolism and preadipocyte differentiation. Biochem. Biophys. Res. Commun. 266, 1–4. Kregel, K.C., 2002. Invited review: heat shock proteins: modifying factors in physiological stress responses and acquired thermotolerance. J. Appl. Physiol. 92, 2177–2186. Lang, R.P., Bayne, C.J., Camara, M.D., Cunningham, C., Jenny, M.J., Langdon, C.J., 2009.

5. Conclusion In conclusion, based on comparative transcriptome analysis, we proposed a comprehensive scenario involving the heat stress response of M. japonicus (Fig. 5). M. japonicus can generate adaptative responses to maintain physiological homeostasis under heat stress via the transcriptional up-regulation of heat shock proteins and antioxidant enzymes, and enhancement of fatty acid metabolism. Nevertheless, acute and prolonged exposure to heat stress resulted in continual repression of the transcription of numerous immune molecules transcription level and led to the transcriptional inhibition of key enzymes and components in the glycometabolism and cellular energy production pathways. The potential dysfunction of the immune response and energy metabolism under heat stress might contribute to the mass mortality of M. japonicus during hot summer months. Considering that a large variety of biological processes are affected by heat stress, multiple aspects (such as strengthening immunity and antioxidant ability, improving dietary fatty acid and carbohydrate composition, etc.) should be considered when constructing strategies to prevent the mass mortality of M. japonicus during hot summer months.

Conflicts of interest The authors declare that they have no conflicts of interest.

345

Aquaculture 502 (2019) 338–346

J. Zheng et al.

Roach, P.J., Skurat, A.V., Harris, R.A., 2011. Regulation of glycogen metabolism. Compr. Physiol. 19, 609–647. Rungrassamee, W., Leelatanawit, R., Jiravanichpaisal, P., Klinbunga, S., Karoonuthaisiri, N., 2010. Expression and distribution of three heat shock protein genes under heat shock stress and under exposure to Vibrio harveyi in Penaeus monodon. Dev. Comp. Immunol. 34, 1082–1089. Shao, Y., Li, C., Chen, X., Zhang, P., Li, Y., Li, T., Jiang, J., 2015. Metabolomic responses of sea cucumber Apostichopus japonicus to thermal stresses. Aquaculture 435, 390–397. Shiel, B.P., Hall, N.E., Cooke, I.R., Robinson, N.A., Strugnell, J.M., 2015. De novo characterisation of the greenlip abalone transcriptome (Haliotis laevigata) with a focus on the heat shock protein 70 (HSP70) family. Mar. Biotechnol. 17, 23–32. Smith, S., Witkowski, A., Joshi, A.K., 2003. Structural and functional organization of the animal fatty acid synthase. Prog. Lipid Res. 42, 289–317. Suzuki, H., Kawarabayasi, Y., Kondo, J., Abe, T., Nishikawa, K., Kimura, S., Hashimoto, T., Yamamoto, T., 1990. Structure and regulation of rat long-chain acyl-CoA synthetase. J. Biol. Chem. 265, 8681–8685. Sweetlove, L.J., Beard, K.F.M., Nunes-Nesi, A., Fernie, A.R., Ratcliffe, R.G., 2010. Not just a circle: flux modes in the plant TCA cycle. Trends Plant Sci. 15, 462–470. Tassanakajon, A., Somboonwiwat, K., Supungul, P., Tang, S., 2013. Discovery of immune molecules and their crucial functions in shrimp immunity. Fish Shellfish Immunol. 34, 954–967. Thorne, M.A.S., Burns, G., Fraser, K.P.P., Hillyard, G., Clark, M.S., 2010. Transcription profiling of acute temperature stress in the Antarctic plunderfish Harpagifer antarcticus. Mar. Genomics 3, 35–44. Tsoi, K.H., Wang, Z.Y., Chu, K.H., 2005. Genetic divergence between two morphologically similar varieties of the kuruma shrimp Penaeus japonicus. Mar. Biol. 147, 367–379. Tsoi, K.H., Chan, T.Y., Chu, K.H., 2007. Molecular population structure of the kuruma shrimp Penaeus japonicus species complex in western Pacific. Mar. Biol. 150, 1345–1364. Wakil, S.J., 1989. Fatty acid synthase, a proficient multifunctional enzyme. Biochemistry 28, 4523–4530. Wang, F.I., Chen, J.C., 2006. The immune response of tiger shrimp Penaeus monodon and its susceptibility to Photobacterium damselae subsp. damselae under temperature stress. Aquaculture 25, 34–41. Wang, W.N., Wang, A.L., Liu, Y., Xiu, J., Liu, Z.B., Sun, R.Y., 2006. Effects of temperature on growth, adenosine phosphates, ATPase and cellular defense response of juvenile shrimp Macrobrachium nipponense. Aquaculture 256, 624–630. Xing, L., Sun, L., Liu, S., Li, X., Zhang, L., Yang, H., 2018. Transcriptome analysis provides insights into the mechanism of albinism during different pigmentation stages of the albino sea cucumber Apostichopus japonicus. Aquaculture 486, 148–160. Xiong, Y., Huang, J., Li, X., Zhou, L., Dong, F., Ye, H., Gan, L., 2014. Deep sequencing of the tilapia (Oreochromis niloticus) liver transcriptome response to dietary protein to starch ratio. Aquaculture 433, 299–306. Xu, D., Zhou, S., Yang, H., 2017. Carbohydrate and amino acids metabolic response to heat stress in the intestine of the sea cucumber Apostichopus japonicus. Aquac. Res. 48, 5883–5891. Xue, S., Liu, Y., Zhang, Y., Sun, Y., Geng, X., Sun, J., 2013. Sequencing and de novo analysis of the hemocytes transcriptome in Litopenaeus vannamei response to white spot syndrome virus infection. PLoS One 8, e76718. Zhao, N., Hou, M., Wang, T., Chen, Y.F., Lv, Y., Li, Z.R., Zhang, R., Xin, W.T., Zou, X.Y., Hou, L., 2014. Cloning and expression patterns of the brine shrimp (Artemia sinica) glycogen phosphorylase (GPase) gene during development and in response to temperature stress. Mol. Biol. Rep. 41, 9–18.

Transcriptome profiling of selectively bred Pacific oyster Crassostrea gigas families that differ in tolerance of heat shock. Mar. Biotechnol. 11, 650–668. Lesser, M.P., 2006. Oxidative stress in marine environments: biochemistry and physiological ecology. Annu. Rev. Physiol. 68, 253–278. Li, F.H., Xiang, J.H., 2013. Recent advances in researches on the innate immunity of shrimp in China. Dev. Comp. Immunol. 39, 11–26. Li, C.Z., Weng, S.P., Chen, Y.G., Yu, X.Q., Lü, L., Zhang, H.Q., He, J.G., Xu, X.P., 2012. Analysis of Litopenaeus vannamei transcriptome using the next-generation DNA sequencing technique. PLoS One 7, e47442. Liang, S., Luo, X., You, W., Luo, L., Ke, C., 2014. The role of hybridization in improving the immune response and thermal tolerance of abalone. Fish Shellfish Immunol. 39, 69–77. Lim, H.J., Kim, B.M., Hwang, I.J., Lee, J.S., Choi, I.Y., Kim, Y.J., Rhee, J.S., 2016. Thermal stress induces a distinct transcriptome profile in the Pacific oyster Crassostrea gigas. Comp. Biochem. Physiol. D Genomics and Proteomics 19, 62–70. Limón-Pacheco, J., Gonsebatt, M.E., 2009. The role of antioxidants and antioxidant-related enzymes in protective responses to environmentally induced oxidative stress. Mutat. Res. Genet. Toxicol. Environ. 674, 137–147. Livak, K..J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using realtime quantitative PCR and the 2-ΔΔCt method. Methods 25, 402–408. Lu, Q.P., Fang, B., Jiang, L.X., Liu, J., 2007. The effect of temperature on selected immune parameters of the white shrimp, Litopenaeus vannamei. J. World Aquacult. Soc. 38, 326–332. Lushchak, V.I., 2011. Environmentally induced oxidative stress in aquatic animals. Aquat. Toxicol. 101, 13–30. Lv, J., Liu, P., Gao, B., Wang, Y., Wang, Z., Chen, L., Li, J., 2014. Transcriptome analysis of the Portunus trituberculatus: de novo assembly, growth-related gene identification and marker discovery. PLoS One 9, e94055. Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628. Munday, P.L., Kingsford, M.J., O'callaghan, M., Donelson, J.M., 2008. Elevated temperature restricts growth potential of the coral reef fish Acanthochromis polyacanthus. Coral Reefs 27, 927–931. Nguyen, C., Nguyen, T.G., Van Nguyen, L., Pham, H.Q., Nguyen, T.H., Pham, H.T., Nguyen, H.T., Ha, T.T., Dau, T.H., Vu, H.T., Nguyen, D.D., Nguyen, N.T.T., Nguyen, N.H., Quyen, D.V., Chu, H.H., Dinh, K.D., 2016. De novo assembly and transcriptome characterization of major growth-related genes in various tissues of Penaeus monodon. Aquaculture 464, 545–553. Ntambi, J.M., 1995. The regulation of stearoyl-CoA desaturase (SCD). Prog. Lipid Res. 34, 139–150. Ntambi, J.M., Miyazaki, M., Dobrzyn, A., 2004. Regulation of stearoyl-CoA desaturase expression. Lipids 39, 1061–1065. Nunes-Nesi, A., Araújo, W.L., Obata, T., Fernie, A.R., 2013. Regulation of the mitochondrial tricarboxylic acid cycle. Curr. Opin. Plant Biol. 16, 335–343. Oexle, H., Gnaiger, E., Weiss, G., 1999. Iron-dependent changes in cellular energy metabolism: influence on citric acid cycle and oxidative phosphorylation. BBA Bioenergetics 1413, 99–107. Pankhurst, N.W., Munday, P.L., 2011. Effects of climate change on fish reproduction and early life history stages. Mar. Freshw. Res. 62, 1015–1026. Poloczanska, E.S., Brown, C.J., Sydeman, W.J., Burrows, M., 2013. Global imprint of climate change on marine life. Nat. Clim. Chang. 3, 919–925. Rebl, A., Verleih, M., Köbis, J.M., Kühn, C., Wimmers, K., Köllner, B., Goldammer, T., 2013. Transcriptome profiling of gill tissue in regionally bred and globally farmed rainbow trout strains reveals different strategies for coping with thermal stress. Mar. Biotechnol. 15, 445–460.

346