Skip to main content

Integrated application of multi-omics provides insights into cold stress responses in pufferfish Takifugu fasciatus



T. fasciatus (Takifugu fasciatus) faces the same problem as most warm water fish: the water temperature falls far below the optimal growth temperature in winter, causing a massive death of T. fasciatus and large economic losses. Understanding of the cold-tolerance mechanisms of this species is still limited. Integrated application of multi-omics research can provide a wealth of information to help us improve our understanding of low-temperature tolerance in fish.


To gain a comprehensive and unbiased molecular understanding of cold-tolerance in T. fasciatus, we characterized mRNA-seq and metabolomics of T. fasciatus livers using Illumina HiSeq 2500 and UHPLC-Q-TOF MS. We identified 2544 up-regulated and 2622 down-regulated genes in the liver of T. fasciatus. A total of 40 differential metabolites were identified, including 9 down-regulated and 31 up-regulated metabolites. In combination with previous studies on proteomics, we have established an mRNA-protein-metabolite interaction network. There are 17 DEMs (differentially-expressed metabolites) and 14 DEGs-DEPs (differentially co-expressed genes and proteins) in the interaction network that are mainly involved in fatty acids metabolism, membrane transport, signal transduction, and DNA damage and defense. We then validated a number of genes in the interaction network by qRT-PCR. Additionally, a number of SNPs (single nucleotide polymorphisms) were revealed through the transcriptome data. These results provide key information for further understanding of the molecular mechanisms of T. fasciatus under cold stress.


The data generated by integrated application of multi-omics can facilitate our understanding of the molecular mechanisms of fish response to low temperature stress. We have not only identified potential genes and SNPs involved in cold tolerance, but also show that some nutrient metabolites may be added to the diet to help the overwintering of T. fasciatus.


Temperature has profound effects on the physiology and behavior of ectotherms, especially teleosts [1]. Within the appropriate temperature range, fish have the capacity to cope with natural temperature changes, such as daily and seasonal changes [2]. However, low temperature can cause irreversible damage [3], especially in farmed species that are unable to escape unfavorable environmental conditions. Low temperature causes alterations in the cellular biological functions, structures, and metabolism, as well as folding, assembly, activity, and stability of proteins [4, 5]. However, we have only poor knowledge about the mechanisms underlying these organisms’ responses to low temperature. Investigating how and to what extent organisms adapt to low temperature at various molecular and physiological levels could provide a better understanding of the molecular mechanisms underlying environmental adaptation.

It has been shown that fish can gradually establish cold adaptive phenotypes through extensive biochemical, metabolic and physiological regulations [6]. Well-defined biochemical and physiological acclimations include producing temperature-specific isozymes [7], altering the content of membrane lipid and the degree of fatty acid unsaturation [8], recruiting different muscle fiber types [9], synthesizing molecular chaperones [10], and altering mitochondrial densities and their properties [11]. However, this amount of information is very limited for a comprehensive understanding of the changes in the organisms. Recently, a rise of omics applications has made it convenient to acquire bioinformation in many studies, especially the screening of key pathways or major genes for the regulation of economically-important traits. Transcriptomics has already been employed to assess responses to cold environment [12, 13]. A comparison of the zebrafish and tilapia transcriptomes revealed that the FoxO signal regulates the low-temperature limit [14]. The comparative studies of transcriptomes and microRNAomes showed that notothenioid lineages were adapted to cold by evolutionary modulation of the multi-functional TGF-β signaling pathway [15]. Metabolomics has also been applied to study the metabolic differences of gilthead sea bream under low-temperature stress, and that study found that the most perturbed metabolic pathways are related to methionine cycle in liver at low water temperature (11 °C) [16]. In the past decade several studies have elaborated on the mechanisms involved in the responses of fish to low temperature. However, the molecular mechanisms of teleost’s response to low-temperature stress have not yet been analyzed using multi-omics.

The biological phenomena are complex and variable, and the regulation of gene expression is complicated. When conducting a single-omics study, the conclusions are often not comprehensive. Therefore, there is a bottleneck in the study of single omics. Multi-omics combines two or more -omics research methods, such as genomics, transcriptomics, proteomics, and/or metabolomics. In teleosts, the combination of transcriptomics and proteomics has revealed some important biological phenomena. For example, it has highlighted the molecular mechanism of dopamine in reproduction of goldfish [17]; revealed that tilapia gut and liver may collaborate immunologically [18]; revealed that dietary vegetable oil can cause indelible injuries to intestines of Atlantic salmon (Salmo salar) [19]; and explored the differential genes and proteins in the red skeletal muscle of rainbow trout (Oncorhynchus mykiss) and the gilthead seabream (Sparus aurata) during sustained swimming [20]. In addition, the combination of proteomics and metabolomics can be used to discriminate fish liver tumors in wild flatfish [21]. Transcriptomics and metabolomics have also been used to identify health impacts of exposure to copper in fish [22]. Unfortunately, the use of multi-omics to study the response mechanisms of fish under low-temperature stress is still close to negligible, which greatly hinders our comprehensive understanding of this biological phenomenon.

As an anadromous and widely distributed species in the South China Sea, the East China Sea, and inland waters in China and Korean Peninsula [23], Takifugu fasciatus (T. fasciatus) (pufferfish) is an important farmed fish of high commercial value [24]. The consumption of T. fasciatus has a long history in the middle and lower reaches of the Yangtze River in China [25]. Furthermore, the pufferfish is also popular in Korea [26] and Japan [27]: pufferfish sashimi is widely known. As a kind of warm-water fish [28], T. fasciatus has an optimum growing temperature between 23 °C and 32 °C [29]. However, in farms as well as rivers, the water temperature in winter is far below that optimal growth temperature, which causes massive deaths of T. fasciatus and leads to large economic losses. Therefore, it is necessary to understand the mechanisms of low-temperature tolerance in T. fasciatus.

In the present study, transcriptomics and metabolomics were used to conduct comparative quantitative omics analysis of T. fasciatus liver after exposure to low-temperature water, and the results were combined with our previous proteomics research to further broaden our perspective on the mechanisms of T. fasciatus responding to low-temperature stress. The results not only provided new and important information about the cold response of T. fasciatus, but also illustrated more completely the relationship between warm-water fish and low-temperature environments.


Transcriptome analysis of T. fasciatus liver under cold stress

In order to identify mRNA expression profiles in T. fasciatus liver under cold stress six cDNA libraries representing the CG (control groups: G1, G2, G3) and EG (experimental group: Ga, Gb, Gc) were constructed with total RNA and subjected to Illumina deep sequencing. In total, 58,246,984, 62,328,458, and 49,318,414 raw reads were obtained from the CG libraries G1, G2, and G3, respectively, and 71,189,264, 51,105,490, and 61,379,788 raw reads were obtained from the libraries Ga, Gb, and Gc, respectively. After quality filtering, approximately 96.46, 97.07, and 96.56% clean reads from CG and 97.01, 96.13, 97.15% clean reads from EG libraries qualified for assembling (Additional file 1: Table S1). These clean reads were mapped to the genome of T. fasciatus; approximately 89.94% (G1), 89.54% (G2), and 89.42% (G3) clean reads from CG and 90.03% (Ga), 90.34% (Gb), 89.91% (Gc) clean reads from EG libraries were located (Additional file 2: Table S2). In all, we identified 20,991 genes, of which 19,924 were able to be matched and annotated by the genome (Additional file 3: Table S3). Subsequently, 2544 up-regulated and 2622 down-regulated genes (differentially expressed genes, DEGs) were identified in T. fasciatus liver, which showed significantly different expression between CG and EG (Additional file 4: Figure S1).

A total of 273 GO clusters that changed significantly (p-value < 0.05) after low temperature stress were identified (Additional file 5: Table S4). The gene ontology (GO) analysis of DEGs produced three major functional categories: biological process, cellular component, and molecular function. Among the molecular function category, “binding” was most commonly represented, followed by “ion binding”. Gene involved in the “intracellular” and “intracellular part” were notably represented in the cellular components. In the category of biological process, a significant proportion of clusters were assigned to “cellular response to DNA damage stimulus” and “macromolecule catabolic process” (Additional file 6: Figure S2). By performing KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses, a total of 12 pathways (Additional file 7: Table S5) that changed significantly (p -value < 0.05) after the low-temperature treatment were identified. Among these pathways, ribosome biogenesis in eukaryotes, glyoxylate and dicarboxylate metabolism, RNA transport, PPAR signaling pathway, and fatty acid biosynthesis, among others, were significantly enriched.

Metabolome analysis of T. fasciatus liver under cold stress

To understand the metabolome changes associated with T. fasciatus under cold stress, an untargeted metabolomics analysis was performed in livers using a UHPLC-Q-TOF MS platform. A total of 9464 metabolic ion peaks were extracted, including 4085 positive ion peaks and 5379 negative ion peaks. In the positive ion mode, a total of 40 differential metabolites were identified (Table 1), including 9 down-regulated and 31 up-regulated metabolites. The PCA (Principal Component Analysis) model obtained by 7-fold cross-validation showed that the QC (Fig. 1a, b) samples in the positive and negative ion modes were closely clustered, indicating that the experiment was reproducible. The OPLS-DA (orthogonal partial least-squares-discriminant analysis) model of the example set was established (Fig. 1c, d); the model evaluation parameters (Positive ion mode: R2Y = 0.82 cum, Q2 = 0.566 cum; Negative ion mode: R2Y = 0.997 cum, Q2 = 0.812 cum) obtained by 7-fold cross-validation indicated that the model is stable and reliable. The permutation test established 200 OPLS-DA models by randomly changing the order of the categorical variables Y to obtain the R2 and Q2 values of the stochastic model (Fig. 1e, f). All Q2 points from left to right are lower than the original blue Q2 points on the right side, indicating that the model is robust and reliable without overfitting. In summary, the instrument analysis system and the test data of this experiment were stable and reliable.

Table 1 Differential metabolites in positive ion mode
Fig. 1

Quality analysis of metabolomics data. a PCA (Principal Component Analysis) score diagram of samples in positive ion mode. b PCA score diagram of samples in negative ion mode. c OPLS-DA (orthogonal partial least-squares-discriminant analysis) score diagram for positive ion mode. d OPLS-DA score diagram for negative ion mode. e OPLS-DA permutation test for positive ion mode. f OPLS-DA permutation test for negative ion mode

Multi-omics identification of key genes and metabolites

In this study, the transcriptome and proteome of T. fasciatus liver under cold stress were integrated by investigating whether changes in the protein levels correlated with changes in the corresponding transcripts. Based on the level of mRNA and protein expression, we classified the genes and proteins involved into five different types (DEGs-DEPs-SameTrend, DEGs-DEPs-Opposite, DEGs-NDEPs, NDEGs-DEPs, NDEPs-NDEGs, Fig. 2a; Additional file 8: Table S6), and the subsequent analysis was based on this classification. A total of 36 DEPs (differentially expressed proteins) were matched with corresponding DEGs, which included 18 up-regulated DEGs-DEPs and 18 down-regulated DEGs-DEPs that were well matched with the tendency of change in abundance of DEGs. On the other hand, 19 DEPs had opposite trends to the DEGs (Additional file 8: Table S6). Investigating the 55 co-expressed genes (DEGs-DEPs) revealed that they were involved in 81 metabolic pathways. Subsequently, combined proteome and metabolome analyses revealed that a total of 41 pathways were involved in DEPs-DEMs (associated DEPs and DEMs).

Fig. 2

Venn diagrams describing the numbers of mRNAs, proteins and pathways that varied between EG (experimental group) and CG (control group). a the number of transcriptome and proteome at the quantitative and differential expression levels. b the number of pathways involved in DEGs-DEPs (differentially co-expressed genes and proteins) or DEPs-DEMs (associated DEPs and DEMs) and their commonalities

Genes produce proteins through complex transcription and translation processes to regulate the metabolism of organisms. Therefore, to obtain an overview of the correlation between the level of transcription, protein and metabolism, the correlation of 81 pathways (related to DEGs-DEPs) and 41 pathways (related to DEPs-DEMs) was analyzed. A total of 20 identical pathways (Fig. 2b) are involved in DEGs, DEPs and DEMs (Table 2). These pathways contained metabolic processes (e.g. pyrimidine metabolism; biosynthesis of unsaturated fatty acids; amino sugar and nucleotide sugar metabolism; valine, leucine and isoleucine degradation; glutathione metabolism; lysosomal degradation), membrane transport (e.g. ABC transporters), signal transduction (e.g. cAMP signaling pathway), organism systems (e.g. bile secretion; vitamin digestion and absorption; retrograde endocannabinoid signaling; serotonergic synapse; regulation of lipolysis in adipocytes; aldosterone synthesis and secretion; longevity regulating pathway) and diseases (e.g. central carbon metabolism in cancers; morphine addiction; alcoholism; Parkinson’s disease; insulin resistance).

Table 2 DEGs, DEPs and DEMs involved in the pathways

There were 17 DEMs (14 up-regulated with fold change > 1, and 3 down-regulated with fold change < 1) and 14 DEGs-DEPs (8 up-regulated, 3 down-regulated, 3 opposite) in these 20 pathways (Table 2). The interaction between them was clustered mainly into two clusters. One of these was related mainly to the transmembrane transport of bile salts (Fig. 3a), and the other to unsaturated fatty acids, vitamins and adenosine (Fig. 3b).

Fig. 3

Interaction networks of the significant differentially co-expressed genes and metabolites analyzed by Cytoscape software (version 3.0.1). a Interaction networks related to the transmembrane transport of bile salts. b Interaction networks related to unsaturated fatty acids, vitamins and adenosine. Red indicates metabolite; green indicates co-expressed gene; blue indicates pathway

Validation of selected genes by qPCR

We performed real-time qPCR analysis on some selected genes to provide additional mRNA transcript levels to validate our RNAseq results. The results of qPCR analysis revealed that most of these genes shared similar expression tendencies (log2fold change; EG vs CG) with those from transcriptomic data (Additional file 9: Table S7). Although some quantitative differences did exist between the two analytical platforms, the common characteristics of the proteomic data and the qRT-PCR analysis supported the reproducibility and reliability of the transcriptomic data.

SNP detection and validation

For extended application of the RNA-Seq data, structural variations were discovered using the assembled transcriptome. Information on all SNPs detected was presented in Additional file 10: Table S8. A total of 23 potential SNPs were identified by detection in 5 co-expressed genes: 3 up-regulated and 2 down-regulated genes (Additional file 11: Table S9). Twenty-three SNPs were tested for experimental validation in fifty T. fasciatus. Among the 23 primer pairs designed for SNP validation, 19 could amplify target sequences. Within these amplified sequences, 13 SNPs were validated, suggesting that approximately 50% of the predicted SNPs were indeed true SNPs. The primers of validated SNPs are listed in Additional file 12: Table S10.


In the present study, low temperature stress was associated with numerous co-expressed DEGs-DEPs and DEMs mainly involved in fatty acid metabolism, oxidative stress, immune system, membrane transport and signal transduction of T. fasciatus (Fig. 3; Table 2). These systems strongly responded to low-temperature stress, and the further exploration of cold stress may deepen our understanding of the molecular mechanisms in fish adaptation (Fig. 4).

Fig. 4

Model of possible molecular mechanisms of T. fasciatus under cold stress

Fatty acids metabolism

The physical properties of phospholipid membranes are heavily dependent upon the saturation of their constituent fatty acids [30]. Maintaining an appropriate balance between saturated and unsaturated fatty acids, under a variable dietary supply, is therefore an essential compositional requirement for all living organisms. Both arachidonic acid (ARA) and (4Z,7Z,10Z,13Z,16Z,19Z)-4,7,10,13,16,19-docosahexaenoic acid (DHA) are important components of membrane phospholipids [31]. The enrichment of ARA and DHA is mainly related to Acyl-CoA desaturase (ACAD) [31], an enzyme involved in maintaining the balance between saturation and unsaturation of fatty acids. Under cold stress, ACAD was up-regulated in the liver of T. fasciatus, indicating that the production of unsaturated fatty acids was strengthened. This is consistent with the findings on carp responses to low-temperature stress [32].

In our metabolomics study, two types of unsaturated fatty acids showed opposite trends (ARA-down, DHA-up), perhaps because there is an antagonistic relationship between ARA and DHA [33, 34]. This has been reported in yellowtail liver, where ARA was incorporated into phosphatidylinositol more effectively than DHA, and the latter inhibited incorporation of ARA into phosphatidylcholine but not phosphatidylinositol [35]. Of course, a more credible reason was that in the combined multi-omics analysis, the metabolism of ARA was affected by other genes and pathways. Compared to the uncertainty around ARA, the up-regulation of DHA was undoubtedly the response of fish enhancing their tolerance of low-temperature environmental stress [36, 37]. In the endocrine system, ARA is regulated by hormone-sensitive lipase (HSL) through two pathways (Fig. 3b). HSL is a key enzyme in the mobilization of fatty acids from intracellular stores [38]. The findings on grass carp [39] and gilthead sea-bream [40] confirmed that dietary ARA can significantly stimulate HSL expression. On the other hand, in the nervous system, much evidence has accumulated indicating that the guanine nucleotide-binding protein subunit beta-4 (G protein) has a positive effect on the metabolism of ARA [41]. Therefore, down-regulation of this G protein would undoubtedly affect the metabolism of ARA. These results indicate that the metabolism of ARA in T. fasciatus subjected to low-temperature stress is regulated by various pathways and genes.

In terms of fatty acid synthesis, Acetoacetyl-CoA synthetase (AACS) is a ketone body-utilizing enzyme responsible for the synthesis of cholesterol and fatty acids from ketone bodies in lipogenic tissues such as liver and adipocytes [42]. In the present study, expression of AACS was down-regulated at the mRNA level, but was up-regulated at the protein level. Complex post-transcriptional mechanisms can be responsible for the inconsistency between mRNA and protein synthesis [43]. Previous research found that leucine supplementation was effective in increasing adiponectin concentration and in reducing total cholesterol concentration in rats [44]. Therefore, the up-regulation of AACS and L-leucine in T. fasciatus liver was associated with regulation of the synthesis of cholesterol and fatty acids in response to low-temperature stress.

Membrane transport and signal transduction

Several lines of evidence clearly indicate that ABC transporters play an important role in the antioxidant mechanism of fish [45]. The bile salt export pump (BSEP) is an important member of the ABC transporter pathway. It is a membrane glycoprotein located on the hepatocyte membrane that mediates ATP-dependent transport of bile salt from the cells to the capillary bile duct [46]. Up-regulation of BSEP can enhance the transport of bile salts to prevent the production of oxygen free radicals and to decrease lipid peroxides and mitochondrial respiratory chain disorders, thereby protecting the function of biofilm and strengthening the immunity of organisms [47]. The blockage of bile salt transport by biliary obstruction induced a variety of liver diseases in humans and mice [48]. Therefore, the antioxidant role of BSEP in fish may be consistent with that in humans and mice.

In the ABC transporter pathway several organic anions are excreted into the bile via a canalicular multispecific organic anion transporter (cMOAT) [49]. The regulation of cMOAT expression is closely related to the liver function of an organism. Compared with the control group, synthesis of cMOAT decreased when the mRNA was unchanged or up-regulated, suggesting a post-transcriptional regulation mechanism. Further analysis found that the metabolic differences caused by low-temperature stress were directly related to the expression of cMOAT and BSEP. Our results show that up-regulation of BSEP was involved in the transport of oligosaccharides (D-Maltose and maltotriose), amino acids (L-lysine and L-leucine) and monosaccharides (D-mannose), whereas up-regulation of cMOAT was involved in the transport of mineral and organic ion (betaine) [50]. Among these up-regulated metabolites, D-mannose can enhance the organism immunity [51]; D-maltose and maltotriose can enhance energy metabolism; and L-leucine is mainly involved in the synthesis of cholesterol and fatty acids. The downregulation of betaine negatively affects fish growth, antioxidant defense and fatty acid synthesis [52]. These differential metabolites reflect the complex transport compensation mechanism of T. fasciatus in the response to low-temperature stress. Not only that, cMOAT and BSEP are related to the metabolism of taurocholate in the bile secretion pathways. Taurocholate supplementation has a potentiating effect on turbot (Scophthalmus maximus) [53], but has a negative impact on the growth of rainbow trout (Oncorhynchus mykiss) [54]. Therefore, it remains uncertain whether taurocholate up-regulation has any effect on the growth performance of T. fasciatus under low-temperature stress. However, taurocholate might play an important role in teleost fish in response to environmental temperature changes because rainbow trout - a cold-water fish - enhanced the metabolism of taurocholate under heat stress [55], whereas T. fasciatus - a warm-water fish - also strengthened the metabolism of taurocholate under cold stress. These results reflect the membrane transport mechanisms in response of T. fasciatus to low-temperature stress.

DNA damage and defense

An extent of damage in an organism is manifestation of its capacity to cope with stress. According to reports in the literature, uridine phosphorylase2 (UPase2) plays critical roles in the metabolism of pyrimidines [56]. UPase2 regulates the balance of uridine concentration in plasma and tissues by catalyzing the dissimilation and phosphorylation of uridine to uracil for nucleic acid recovery and reuse [57]. Therefore, when UPase2 protein is up-regulated, it indicates that the catalytic reaction of uridine is strengthening, which leads to a decrease in uridine. Hence, in previous studies, when low-temperature stress was found to cause DNA damage in T. fasciatus [58], this might have been caused by an imbalance in uridine metabolism.

The immune and antioxidant systems also play pivotal roles in resisting potential damage caused by environmental stress. Vitamin intake is an important part of teleost immunity [59]. Retinol-binding protein (RBP) is a carrier of vitamin A that is mainly synthesized and released in liver [60]. RBP binds retinol and retinal, and as such is thought to play an important role in vitamin A (retinol) uptake, transport and metabolism [61]. Up-regulation of RBP enhances uptake of vitamin A to enhance fish immunity and growth [62]. For example, RBP was up-regulated in Channa striatus under high-temperature stress [63], suggesting that RBP has a positive effect on teleosts facing environmental temperature changes. In the vitamin digestion and absorption pathway, pantotherate (vitamin B5) and nicotinamide (vitamin B3) are up-regulated at low temperature, suggesting they might enhance low-temperature tolerance of T. fasciatus along with vitamin C [58]. Of course, the dosage of vitamins and their mechanism of action still need further study.

Carbohydrate metabolism is also a part of organism immunity. D-mannose is a hexose stimulating liver to secrete mannose-binding lectin (MBL), ultimately strengthening the immune system [64]. N-acetyl-D-glucosamine (GlcNAc) kinase is a hexokinase involved in carbohydrate metabolism. Considering our experimental results, D-mannose can be phosphorylated by GlcNAc kinase to form mannose-6-phosphate [65]. In mammals, mannose has a protective effect on an intestinal mucosal immune barrier in rats with severe acute pancreatitis [66]. On the other hand, acid phosphatase (ACP) is unique among lysosomal enzymes in that it has high concentration of both mannose and complex type sugars chains, and oligosaccharide chains of lysosomal enzymes contain a large proportion of mannose [67]. Therefore, the down-regulation of ACP decreases dephosphorylation of mannose [68]. Similarly, under low-temperature stress, ACP in sea bream (Sparus aurata) blood also decreased [69]. Interestingly, down-regulated ACP and GlcNAc kinase have opposite functions, namely dephosphorylation and phosphorylation, and the corresponding target metabolic substrate is up-regulated D-mannose. However, there is no evidence to support the competitive relationship between ACP and GlcNAc kinase. Obviously, the metabolic mechanism of D-mannose is still unclear under low-temperature stress, and more research is needed.

Oxidative stress is a physiological response of organisms to environmental changes. Detoxification-related protein glutathione S-transferase (GST) is an important indicator of the level of oxidative stress [70]. In our results, GST was up-regulated at both the mRNA and protein levels, indicating enhanced detoxification processes to remove ROS (reactive oxygen species). Correspondingly, the metabolism of glutathione disulfide was also enhanced as revealed by metabolomics. GST prevents damage to the cell membrane and other macromolecules by promoting the attachment of tripeptide GSH to a number of potentially harmful electrophilic substrates [71]. In Danio rerio, deficiency in the glucose transporter slc2al (Solute carrier family 2, facilitated glucose transporter member 1) has been shown to cause a suite of neural defects [72]. Therefore, up-regulation of sla2a1 in T. fasciatus may indicate a need to strengthen the nervous system under low-temperature stress. Not only that, Slc2a1 can also enhance the metabolism of L-leucine and acetylcarnitine in the cancer pathway and the insulin resistance pathway, respectively. However, this conclusion is derived from studies of human diseases and has not been confirmed in teleosts. A possibility needs to be explored that they may have a similar relationship in teleosts because we found that L-leucine could increase growth performance, improve metabolic activities and enhance non-specific immunities in tilapia [73] and black carp juveniles [74], and acetylcarnitine has antioxidant effects in teleosts [75].


The results showed that transcriptome, proteome and metabolome of T. fasciatus liver changed significantly under low-temperature stress. Although more than 5000 genes were differentially expressed, only 55 relevant proteins (DEGs-DEPs) changed. Meanwhile, a total of 40 DEMs were found at the metabolic level. Through interaction analysis, we found that DEGs-DEPs and DEMs mainly regulate immunity, growth and antioxidant capacity by enhancing metabolism of unsaturated fatty acids, transport of bile salts, vitamin uptake, and antioxidant capacity, which ultimately affected cold tolerance of T. fasciatus (Fig. 4). Co-expressed DEGs-DEPs in these biochemical processes can serve as potential indicators for regulation of cold-tolerant traits, and DEMs can be used as feed additives to investigate whether they can enhance cold tolerance of T. fasciatus. These findings provide new insights into the cold stress responses in fish and could promote research into the molecular breeding of pufferfish Takifugu fasciatus.


Experimental design and sampling

The sampling time point of stress was chosen basing on previous studies [29] and the results of our preliminary-experiment (Additional file 13: Table S11). The results indicated that the time point of 24 h under low temperature stress could be used as a breakthrough point in this study.

Taking experiment needs and the suitable culture density of juvenile T. fasciatus into consideration [28], a total of 120 juvenile T. fasciatus (13 ± 1.76 cm in length, 22 ± 2.85 g in weight), supplied from Zhongyang Group Co., Ltd. (Jiangsu Province, China), were randomly transferred to six aquaria with biofiltered water recirculation systems (equipped with cooling and heating functions; volume 200 L; flow rate 5 L/min), and were divided into the control group (G1, G2 and G3; total of 60 individuals; 2.2 ± 0.3 g/L) and the experimental group (Ga, Gb and Gc; total of 60 individuals; 2.2 ± 0.3 g/L). These individuals were acclimated at 26 ± 1 °C for 2 weeks. During the acclimation period, the commercial fish diet (including 42% w/w protein and 8.0% w/w fat, supplied by Zhenjiang Jiaji Feed Co. Ltd.) was given twice a day until 24 h before the experimental treatments (Fig. 5).

Fig. 5

General workflow and summary of the present study

The water temperature of the CG was maintained at 26 ± 1 °C until the sampling was completed. The water temperature of the three experimental tanks (EG) decreased from 26 °C to 19 °C at a rate of 0.85 °C /1 h, and then maintained at 19 °C for 12 h, before decreasing from 19 °C to 12 °C at the same rate. After the EG was maintained at 12 °C for 24 h, both the experimental and control fish were euthanized with MS-222 at 10 mg L− 1, and quickly removed for liver dissection. Samples of EG (Ga, Gb and Gc) and CG (G1, G2 and G3) were done in three biological replicates. Each group (Ga, Gb, Gc, G1, G2 and G3) contains liver tissue from six different individuals, which were used for RNAseq and qPCR analysis. Ten replicates were consisted for both EG and CG respectively for metabolomics experiments. Ten individuals of EG came from Ga (3), Gb (3), and Gc (4). Similarly, 10 individuals from G1 (3), G2 (3) and G3 (4) were used as CG. Tissue samples were transferred to Eppendorf tubes, immediately flash frozen in liquid nitrogen, and then stored at − 80 °C.

RNA isolation

Total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. The total RNA quantity and purity were analysed with a Bioanalyzer 2100 and an RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number > 7.0.

Assembly and functional annotation of the transcriptome sequencing

For constructing six cDNA library, approximately 5 μg of total RNA per sample was used for the RNA sample preparations. The library for sequencing was generated using an Illumina TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA). Transcriptome sequencing was carried out on an Illumina HiSeq 2500 platform that generated approximately 125-bp paired-end (PE) raw reads by LC Sciences (Houston, TX, USA). After removing adaptor sequences, ambiguous ‘N’ nucleotides (with the ratio of ‘N’ greater than 5%) and low quality sequences (with quality score less than 20), the high-quality trimmed sequences were used for further mapping to the T. fasciatus genome (unpublished data) with Hisat 2.0.4.

Differential expression analysis

The expression level of each transcript was measured as the number of clean reads mapped to its sequence. The mapped clean read number was normalized to RPKM (reads per kilobase of transcript per million mapped reads) with RSEM 1.2.3. We used DESeq to determine the FDR threshold. FDR < 0.05 and fold change > 2 were considered to indicate significant expression abundance. All the DEGs identified in this study were used as reference for the enrichment analysis by GO ( and KEGG ( Hypergeometric test was used to identify the overrepresented GO and KEGG pathway terms with a significance level at 0.05, and Beniamini & Hochberg method was used for the correction of the p-values.

Extraction of metabolites

For metabolite extraction, 60 mg of each sample was weighted and placed on ice, then resuspended in 300 μL of 80% (vol/vol) methanol in water and homogenized by using an Ultra-Turrax homogenizer for 30s, with three 5-s pulses per sample, followed by ultrasonication at low temperature for 30 min, 2 times. Samples were kept at − 20 °C for 1 h and then spun at 4 °C for 15 min at 13, 000 rpm. Supernatant (250 μL per sample) was transferred to a new tube and dried down by using a speed vacuum at 30 °C. Dried samples were stored at − 80 °C.

LC-MS/MS analysis (HILIC/MS)

Analyses were performed using an UHPLC (1290 Infinity LC, Agilent Technologies) coupled to a quadrupole time-of-flight (AB Sciex TripleTOF 6600).

For HILIC separation, samples were analyzed using a 2.1 mm × 100 mm ACQUIY UPLC BEH 1.7 μm column (Waters, Ireland). In ESI positive and negative modes, the mobile phase comprised A = 25 mM ammonium acetate and 25 mM ammonium hydroxide in water and B = acetonitrile. The gradients (vol/vol) were 85% B for 1 min, linearly reduced to 65% in 11 min, and then was reduced to 40% in 0.1 min and kept for 4 min, and then increased to 85% in 0.1 min, with a 5 min re-equilibration period employed.

The ESI source conditions were set as follows: Ion Source Gas1 (Gas1) as 60, Ion Source Gas2 (Gas2) as 60, curtain gas (CUR) as 30, source temperature: 600 °C, IonSpray Voltage Floating (ISVF) ±5500 V. In MS only acquisition, the instrument was set to acquire over the m/z range 60–1000 Da, and the accumulation time for TOF-MS scan was set at 0.20s/spectra. In auto MS/MS acquisition, the instrument was set to acquire over the m/z range 25–1000 Da, and the accumulation time for product ion scan was set at 0.05 s/spectra. The product ion scan was acquired using information dependent acquisition (IDA) with a high sensitivity mode selected. The collision energy (CE) was fixed at 35 V with ±15 eV. Declustering potential (DP) was set at ±60 V.

Data processing

The raw MS data (wiff.scan files) were converted to MzXML files using ProteoWizard MSConvert and processed using XCMS for feature detection, retention time correction and alignment. The MS/MS spectra were searched in the in-house Standard MS/MS library. The library contained the MS/MS spectra of approximately 2700 compounds (primarily polar compounds and some nonpolar compounds), which were acquired using standards. The MS/MS spectra matching score was calculated using dot-product algorithm [76] taking fragments and intensities into account. The matching score cutoff was set as 0.8. The MS/MS spectra matching results were confirmed with standards. The data after completeness check, supplement missing value and delete extreme value were normalized between samples and metabolites respectively to ensure comparisons parallel. Finally, all extracted ion peak areas of the data were normalized and the data was subjected to Pareto-scaling processing in the SIMCA-P 14.1 software (Umetrics, Umea, Sweden).

In the extracted ion features, only the variables having more than 50% of the nonzero measurement values in at least one group were kept. For the multivariate statistical analysis, the MetaboAnalyst ( web-based system was used. After the Pareto-scaling, PCA and OPLS-DA were performed. The leave-one-out cross-validation and the response permutation testing were used to evaluate the robustness of the model. The significantly different metabolites were determined based on the combination of a statistically significant threshold of Variable Importance for the Projection (VIP) values obtained from the OPLS-DA model and two-tailed Student’s t test (p value) on the raw data. Metabolites with both multidimensional statistical analysis VIP > 1 and univariate statistical analysis P value < 0.05 were selected as metabolites with significant differences; VIP > 1 and 0.05 < P < 0.1 were used as differential metabolites (DEMs).

Multi-omics analysis

Transcriptome samples used in this study are parallel to the materials researched in previous proteomic study [77]. Samples of EG and CG were done in three biological replicates, each containing liver tissue from six different individuals, and were used to conduct proteomic analysis by iTRAQ technique (Proteome data are available via ProteomeXchange with identifier PXD010955). A total of 160 proteins with fold change > 1.2 and a p -value < 0.05 were considered to be significantly differentially expressed (DEPs). These differential proteins were combined into the differential genes and metabolites obtained in this study to conduct the multi-omics analysis. First, the combined transcriptome and proteome were analyzed to find their co-expressed differential genes (proteins), and related pathways (related to DEGs-DEPs). Secondly, the analyses of proteomes and metabolomes were combined to identify pathways associated with them (related to DEPs-DEMs). Finally, the two result sets were combined to find regulatory pathways involving differentially expressed genes, proteins and metabolites. Moreover, the relationships among significantly different genes, proteins and metabolites were mapped into an interaction network by cytoscape software to provide rich and comprehensive information for revealing the molecular mechanisms of response to low-temperature stress at different levels.

Validation of significant DEGs by RT-qPCR

A total of 11 differentially regulated mRNAs obtained from integrative analysis of multi-omics were verified with quantitative real-time PCR (qRT-PCR). Primers for qRT-PCR were listed in Additional file 14: Table S12. qRT-PCR with β-actin as an internal control was used to explore mRNA expression. qRT-PCR was performed with an SYBR Green Master kit according to the manufacturer’s protocol (Roche, Basel, Switzerland). The experiments were carried out in triplicate with a total volume of 20 μL in an ABI StepOnePlus, containing 10 μL SYBR Green Master, 4 μL cDNA (500 ng), and 3 μL forward and reverse primers (2 μmol/L). The qRT-PCR was programmed at 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s, and 55 °C for 1 min. The expression level was calculated by the 2-CT method and subjected to statistical analysis.

Detection and validation of SNPs in differential genes

By multi-omics analysis we obtained 14 differentially expressed genes, including 8 up-regulated and 3 down-regulated genes (Table 2). SNPs located on these genes have the most potential to be developed into molecular markers for breeding.

Potential SNPs were called using SAMtools software. SAMtools provides various utilities for manipulating alignments in the SAM format, including sorting, merging, indexing, and generating alignments in a preposition format. First, we called SNPs with the SAMtools mpileup utility. We then piped a BCF output file to SAMtools bcftools, which converted the BCF file into a VCF file. Afterwards, we piped the VCF file into with the varFilter-d100 option, which retained SNPs that had read depth higher than 100. In order to validate the accuracy of the SNPs predictions, fin tissue samples from 50 T. fasciatus were collected for this SNP validation. Genomic DNA was extracted from pterygiophore tissue samples, using a cell/tissue genomic DNA extraction kit (centrifugal column type, Generay PBiotech, Shanghai). Primers were designed to amplify the flanking sequence of the selected SNPs using Premier 5.0. The amplified PCR products were sequenced by utilizing an ABI3730 sequencer and the products were analyzed using DNAMAN 8.0.

Availability of data and materials

The dataset generated by deep-sequencing platforms was deposited in the NCBI Sequence Read Archive (SRA, The NCBI Sequence Read Archive accession number is GSE129226.

The mass spectrometry proteomics data supporting the conclusions of this article is available in ProteomeXchange via the PRIDE ( partner repository with the dataset identifier PXD010955.

The metabolite data set supporting the results of this article is included within the article, and can be found in the Supplemental Information (Additional file 15: Table S13).



Acetoacetyl-CoA synthetase


Acyl-CoA desaturase


Acid phosphatase


Arachidonic acid


ATP synthase-coupling factor 6


Bile salt export pump


Canalicular multispecific organic anion transporter


Differentially expressed genes


Differentially co-expressed genes and proteins


Differential metabolites


Differentially expressed proteins


Differentially expressed proteins and genes with opposite trends; DEPs_NDEGs

Differential expressedly proteins with no significant difference in transcription level


differentially expressed proteins and genes with the same trend


Associated differential proteins and differential metabolites


(4Z,7Z,10Z,13Z,16Z,19Z)-4,7,10,13,16,19-docosahexaenoic acid

G proteins:

Guanine nucleotide-binding proteins subunit base 4




Gene Ontology


Glutathione S-transferase


Hormone-sensitive lipase


Kyoto Encyclopedia of Genes and Genomes


Non differentially expressed proteins with significant difference in transcription level


Genes with no significant difference in transcription and protein levels


retinol-binding protein


Reactive oxygen species


Solute carrier family 2, facilitated glucose transporter member 1


Single nucleotide polymorphisms


Ubiquitin-associated protein 1-like


Uridine phosphorylase


Variable Importance for the Projection


  1. 1.

    Ibarz A, Martín-Pérez M, Blasco J, Bellido D, Oliveira ED, Fernández-Borràs J. Gilthead Sea bream liver proteome altered at low temperatures by oxidative stress. Proteomics. 2010;10(5):963–75.

    CAS  PubMed  Google Scholar 

  2. 2.

    Oppedal F, Dempster T, Stien LH. Environmental drivers of Atlantic salmon behaviour in sea-cages: a review. Aquaculture. 2011;311(1–4):1–18.

    Article  Google Scholar 

  3. 3.

    Bennett WA, Judd FW. Comparison of methods for determining low temperature tolerance: experiments with pinfish, Lagodon rhomboides. Copeia. 1992;1992(4):1059–65.

    Article  Google Scholar 

  4. 4.

    Pan J, Tristramnagle S, Kucerka N, Nagle JF. Temperature dependence of structure, bending rigidity, and bilayer interactions of dioleoylphosphatidylcholine bilayers. Biophys J. 2008;94(1):117–24.

    CAS  PubMed  Article  Google Scholar 

  5. 5.

    Simon SA, Advani S, Mcintosh TJ. Temperature dependence of the repulsive pressure between phosphatidylcholine bilayers. Biophys J. 1995;69(4):1473.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  6. 6.

    Yong L, Song G, Yan J, He X, Li Q, Cui Z. Transcriptomic characterization of cold acclimation in larval zebrafish. BMC Genomics. 2013;14(1):612.

    Article  CAS  Google Scholar 

  7. 7.

    Somero GN, Hochachka PW. Biochemical adaptation to the environment. Fish Physiol. 1971;6(1):99–156.

    Google Scholar 

  8. 8.

    Johnston PV, Roots BI. Brain lipid fatty acids and temperature acclimation . Comp Biochem Physiol. 1964;11(3):303–9.

    CAS  PubMed  Article  Google Scholar 

  9. 9.

    Gerlach GF, Turay L, Malik KT, Lida J, Scutt A, Goldspink G. Mechanisms of temperature acclimation in the carp: a molecular biology approach. Am J Physiol. 1990;259(2):237–44.

    Google Scholar 

  10. 10.

    Fader SC, Yu Z, Spotila JR. Seasonal variation in heat shock proteins (hsp 70) in stream fish under natural conditions. Jthermbiol. 1994;19(5):335–41.

    CAS  Google Scholar 

  11. 11.

    St-Pierre J, Charest PM, Guderley H. Relative contribution of quantitative and qualitative changes in mitochondria to metabolic compensation during seasonal acclimation of rainbow trout Oncorhynchus mykiss. J Exp Biol. 1998;201(21):2961–70.

    CAS  Google Scholar 

  12. 12.

    Kassahn KS, Crozier RH, Ward AC, Stone G, Caley MJ. From transcriptome to biological function: environmental stress in an ectothermic vertebrate, the coral reef fish Pomacentrus moluccensis. BMC Genomics. 2007;8(1):358.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  13. 13.

    Mininni AN, Milan M, Ferraresso S, Petochi T, Di MP, Marino G, Livi S, Romualdi C, Bargelloni L, Patarnello T. Liver transcriptome analysis in gilthead sea bream upon exposure to low temperature. BMC Genomics. 2014;15(1):765.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  14. 14.

    Hu P, Liu M, Liu Y, Wang J, Zhang D, Niu H, Jiang S, Wang J, Zhang D, Han B. Transcriptome comparison reveals a genetic network regulating the lower temperature limit in fish. Sci Rep. 2016;6:28952.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  15. 15.

    Xu Q, Cai C, Hu X, Liu Y, Guo Y, Hu P, Chen Z, Peng S, Zhang D, Jiang S. Evolutionary suppression of erythropoiesis via the modulation of TGF-β signalling in an Antarctic icefish. Mol Ecol. 2015;24(18):4664.

    CAS  PubMed  Article  Google Scholar 

  16. 16.

    Melis R, Sanna R, Braca A, Bonaglini E, Cappuccinelli R, Slawski H, Roggio T, Uzzau S, Anedda R. Molecular details on gilthead sea bream (Sparus aurata) sensitivity to low water temperatures from (1) H NMR metabolomics. Comp Biochem Physiol A Mol Integr Physiol. 2017;204:129.

    CAS  PubMed  Article  Google Scholar 

  17. 17.

    Popesku JT, Martyniuk CJ, Denslow ND, Trudeau VL. Rapid Dopaminergic Modulation of the Fish Hypothalamic Transcriptome and Proteome. Plos one. 2010;5(8):e12338.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  18. 18.

    Wu N, Song YL, Wang B, Zhang XY, Zhang XJ, Wang YL, Cheng YY, Chen DD, Xia XQ, Lu YS. Fish gut-liver immunity during homeostasis or inflammation revealed by integrative transcriptome and proteome studies. Sci Rep. 2016;6:36048.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  19. 19.

    Morais S, Silva T, Cordeiro O, Rodrigues P, Guy DR, Bron JE, Taggart JB, Bell JG, Tocher DR. Effects of genotype and dietary fish oil replacement with vegetable oil on the intestinal transcriptome and proteome of Atlantic salmon (Salmo salar). BMC Genomics. 2012;13(1):448.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  20. 20.

    Planas JV, Martínpérez M, Magnoni LJ, Blasco J, Ibarz A, Fernandezborras J, Palstra AP: Transcriptomic and proteomic response of skeletal muscle to swimming-induced exercise in fish: Springer Berlin Heidelberg; 2013.

  21. 21.

    Stentiford GD, Viant MR, Ward DG, Johnson PJ, Martin A, Wenbin W, Cooper HJ, Lyons BP, Feist SW. Liver tumors in wild flatfish: a histopathological, proteomic, and metabolomic study. Omics-a J Integr Biol. 2005;9(9):281–99.

    CAS  Article  Google Scholar 

  22. 22.

    Santos EM, Ball JS, Williams TD, Wu H, Ortega F, Aerle RV, Katsiadaki I, Falciani F, Viant MR, Chipman JK. Identifying health impacts of exposure to copper using transcriptomics and metabolomics in a fish model. Environ Sci Technol. 2010;44(2):820–6.

    CAS  PubMed  Article  Google Scholar 

  23. 23.

    Kato A, Doi H, Nakada T, Sakai H, Hirose S. Takifugu obscurus is a euryhaline fugu species very close to Takifugu rubripes and suitable for studying osmoregulation. BMC Physiol. 2005;5(1):18.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  24. 24.

    Wen X, Wang L, Zhu W, Wang D, Li X, Qian X, Yin S. Three toll-like receptors (TLRs) respond to Aeromonas hydrophila or lipopolysaccharide challenge in pufferfish, Takifugu fasciatus. Aquaculture. 2017;481.

    CAS  Article  Google Scholar 

  25. 25.

    Huang Y, Wang W, Liu Y, Zhang D. The development and Progress of cultured puffer fish in China. J Chin Institute Food Sci Technol. 2018.

  26. 26.

    Hwang SM, Oh KS. Comparisons of food component characteristics of wild and cultured edible pufferfishes in Korea; 2013.

    Google Scholar 

  27. 27.

    Sainen S, Ozawa K, Mineki M, Noguchi T. The sensory characteristics of non-toxic liver of pufferfish (Takifugu rubripes) cultured at aquaria on land by cooking methods. Jpn J Sens Eval. 2009;13(2–2):115–24.

    Google Scholar 

  28. 28.

    Yoo GY, Lee JY. The effect of feeding frequency, water temperature, and stocking density on the growth of river puffer Takifugu obscurus reared in a zero-exchange water system. Fish Aqua Sci. 2016;19(1):1–7.

    Article  CAS  Google Scholar 

  29. 29.

    Cheng CH, Ye CX, Guo ZX, Wang AL. Immune and physiological responses of pufferfish ( Takifugu obscurus ) under cold stress. Fish Shellfish Immunol. 2017;64:137–45.

    CAS  PubMed  Article  Google Scholar 

  30. 30.

    Gennis RB. Introduction: the structure and composition of biomembranes. New York: Springer; 1989.

    Google Scholar 

  31. 31.

    Mutch DM, Grigorov M, Berger A, Fay LB, Roberts MA, Watkins SM, Williamson G, German JB. An integrative metabolism approach identifies stearoyl-CoA desaturase as a target for an arachidonate-enriched diet. FASEB J. 2005;19(6):599–601.

    CAS  PubMed  Article  Google Scholar 

  32. 32.

    Polley SD, Tiku PE, Trueman RT, Caddick MX, Morozov IY, Cossins AR. Differential expression of cold- and diet-specific genes encoding two carp liver delta 9-acyl-CoA desaturase isoforms. Am J Physiol Regul Integr Comp Physiol. 2003;284(1):41–50.

    Article  Google Scholar 

  33. 33.

    Petrik MB, Mcentee MF, Chiu CH, Whelan J. Antagonism of arachidonic acid is linked to the antitumorigenic effect of dietary eicosapentaenoic acid in Apc (min/+) mice. J Nutr. 2000;130(5):1153–8.

    CAS  PubMed  Article  Google Scholar 

  34. 34.

    Watanabe S, Doshi M, Akimoto K, Kiso Y, Hamazaki T. Suppression of platelet-activating factor generation and modulation of arachidonate metabolism by dietary enrichment with (n-9) eicosatrienoic acid or docosahexaenoic acid in mouse peritoneal cells. Prostaglandins Other Lipid Mediat. 2001;66(2):109.

    CAS  PubMed  Article  Google Scholar 

  35. 35.

    Tanaka T, Dai I, Sakamoto M, Takai Y, Morishige JI, Murakami K, Satouchi K. Mechanisms of accumulation of arachidonate in phosphatidylinositol in yellowtail. FEBS J. 2003;270(7):1466–73.

    CAS  Google Scholar 

  36. 36.

    Kanazawa A: Effects of docosahexaenoic acid and phospholipids on stress tolerance of fish. Aquaculture 1997, 155(1–4):129–134.

    CAS  Article  Google Scholar 

  37. 37.

    Watanabe T. Importance of docosahexaenoic acid in marine larval fish. J World Aquacult Soc. 1993;24(2):152–61.

    Article  Google Scholar 

  38. 38.

    Casado ME, Pastor O, Mariscal P, Canfránduque A, Martínezbotas J, Kraemer FB, Lasunción MA, Martínhidalgo A, Busto R. Hormone-sensitive lipase deficiency disturbs the fatty acid composition of mouse testis. Prostaglandins Leukot Essent Fatty Acids. 2013;88(3):227.

    CAS  PubMed  Article  Google Scholar 

  39. 39.

    Tian JJ, Lei CX, Ji H, Kaneko G, Zhou JS, Yu HB, Li Y, Yu EM, Xie J. Comparative analysis of effects of dietary arachidonic acid and EPA on growth, tissue fatty acid composition, antioxidant response and lipid metabolism in juvenile grass carp, Ctenopharyngodon idellus. Br J Nutr. 2017;118(6):411.

    CAS  PubMed  Article  Google Scholar 

  40. 40.

    Alves MD, Rocha F, Martínez-Rodríguez G, Bell G, Morais S, Castanheira F, Bandarra N, Coutinho J, Yúfera M, Conceição LE. Teleost fish larvae adapt to dietary arachidonic acid supply through modulation of the expression of lipid metabolism and stress response genes. Br J Nutr. 2012;108(5):864.

    Article  CAS  Google Scholar 

  41. 41.

    Akiba S, Sato T, Fujii T. Involvement of a guanine-nucleotide-binding protein-mediated mechanism in the enhancement of arachidonic acid liberation by phorbol 12-myristate 13-acetate and Ca2+ in saponin-permeabilized platelets. Biochim Biophys Acta. 1990;1044(3):291–6.

    CAS  PubMed  Article  Google Scholar 

  42. 42.

    Hasegawa S, Yamasaki M, Fukui T. Degradation of acetoacetyl-CoA synthetase, a ketone body-utilizing enzyme, by legumain in the mouse kidney. Biochem Biophys Res Commun. 2014;453(3):631–5.

    CAS  PubMed  Article  Google Scholar 

  43. 43.

    Smillie CL, Sirey T, Ponting CP. Complexities of post-transcriptional regulation and the modeling of ceRNA crosstalk. Crit Rev Biochem Mol Biol. 2018;53(3):231–45.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  44. 44.

    Torres-Leal FL, Fonseca-Alaniz MH, Teodoro GF, Capitani MDD, Vianna D, Pantaleão LC, Matos-Neto EM, Rogero MM, Jr JD, Tirapegui J. Leucine supplementation improves adiponectin and total cholesterol concentrations despite the lack of changes in adiposity or glucose homeostasis in rats previously exposed to a high-fat diet. Nutr. Metab. 2011;8(1):62.

    CAS  Article  Google Scholar 

  45. 45.

    Ferreira M, Costa J, Reishenriques MA. ABC transporters in fish species: a review. Front Physiol. 2014;5:266.

    PubMed  PubMed Central  Article  Google Scholar 

  46. 46.

    And PJM, Stieger B. Bile Salt Transporters. Annu Rev Physiol. 2002;64(64):635–61.

    Google Scholar 

  47. 47.

    Anwer MS. Cellular regulation of hepatic bile acid transport in health and cholestasis. Hepatology. 2004;39(3):581–90.

    CAS  PubMed  Article  Google Scholar 

  48. 48.

    Sohail MI, Schmid D, Wlcek K, Spork M, Szakacs G, Trauner M, Stockner T, Chiba P. Molecular mechanism of taurocholate transport by the bile salt export pump, an ABC-transporter associated with intrahepatic cholestasis. Mol Pharmacol. 2017;92(4):401.

    CAS  PubMed  Article  Google Scholar 

  49. 49.

    Ito K, Suzuki H, Hirohashi T, Kume K, Shimizu T, Sugiyama Y. Molecular cloning of canalicular multispecific organic anion transporter defective in EHBR. Am J Phys. 1997;272(1):16–22.

    Google Scholar 

  50. 50.

    Duan G, Huang ZM, Hai-Long WU, Wang XY, Sun FE. Determination of three kinds of organic acid in betaine by ion-exclusion chromatography. China Surfactant Detergent Cosmetics. 2010;40(6):465–7.

  51. 51.

    Johnson DR, Habeebu SSM, Klaassen CD. Increase in bile flow and biliary excretion of glutathione-derived Sulfhydryls in rats by drug-metabolizing enzyme inducers is mediated by multidrug resistance protein 2. Toxicol Sci. 2002;66(1):16.

    CAS  PubMed  Article  Google Scholar 

  52. 52.

    Adjoumani JY, Wang K, Zhou M, Liu W, Zhang D. Effect of dietary betaine on growth performance, antioxidant capacity and lipid metabolism in blunt snout bream fed a high-fat diet. Fish Physiol Biochem. 2017;43(6):1733–45.

    CAS  PubMed  Article  Google Scholar 

  53. 53.

    Gu M, Bai N, Kortner TM. Taurocholate supplementation attenuates the changes in growth performance, feed utilization, lipid digestion, liver abnormality and sterol metabolism in turbot (Scophthalmus maximus) fed high level of plant protein. Aquaculture. 2016;468:597–604.

    CAS  Article  Google Scholar 

  54. 54.

    Iwashita Y, Suzuki N, Matsunari H, Sugita T, Yamamoto T. Influence of soya saponin, soya lectin, and cholyltaurine supplemented to a casein-based semipurified diet on intestinal morphology and biliary bile status in fingerling rainbow trout Oncorhynchus mykiss. Fish Sci. 2009;75(5):1307–15.

    CAS  Article  Google Scholar 

  55. 55.

    Kemp CJ, Curtis LR. Thermally modulated biliary excretion of [14C] taurocholate in rainbow. Can J Fish Aqua Sci. 2011;44(4):846–51.

    Article  Google Scholar 

  56. 56.

    Yang G, Ziemba A, Turner J, Pizzorno G. Abstract 52: overexpression of uridine phosphorylase 2 (UPP2) in tissues of uridine phosphorylase (UPP1) knock-out mice. Cancer Res. 2010;70(8 Supplement):52.

    Google Scholar 

  57. 57.

    Temmink OH, Bruin MD, Laan AC, Turksma AW, Cricca S, Masterson AJ, Noordhuis P, Peters GJ. The role of thymidine phosphorylase and uridine phosphorylase in (fluoro) pyrimidine metabolism in peripheral blood mononuclear cells. Int J Biochem Cell Biol. 2006;38(10):1759–65.

    CAS  PubMed  Article  Google Scholar 

  58. 58.

    Cheng CH, Liang HY, Luo SW, Wang AL, Ye CX. The protective effects of vitamin C on apoptosis, DNA damage and proteome of pufferfish (Takifugu obscurus) under low temperature stress. J Therm Biol. 2018;71:128.

    CAS  PubMed  Article  Google Scholar 

  59. 59.

    Cerezuela R, Cuesta A, Meseguer J, Esteban MAÁ. Effects of dietary vitamin D administration on innate immune parameters of seabream (Sparus aurata L.). Fish Shellfish Immunol. 2009;26(2):243–8.

    CAS  PubMed  Article  Google Scholar 

  60. 60.

    Smith JE, Borek C, Goodman DWS. Regulation of retinol-binding protein metabolism in cultured rat liver cell lines. Cell. 1978;15(3):865–73.

    CAS  PubMed  Article  Google Scholar 

  61. 61.

    Parmar MB, Shams R, Wright JM. Genomic organization and transcription of the medaka and zebrafish cellular retinol-binding protein ( rbp ) genes. Mar Genomics. 2013;11(5):1–10.

    PubMed  Article  Google Scholar 

  62. 62.

    Guimarães IG, Pezzato LE, Santos VG, Orsi RO, Barros MM. Vitamin a affects haematology, growth and immune response of Nile tilapia (Oreochromis niloticus, L.), but has no protective effect against bacterial challenge or cold-induced stress. Aquac Res. 2016;47(6):2004–18.

    Article  CAS  Google Scholar 

  63. 63.

    Mahanty A, Purohit GK, Banerjee S, Karunakaran D, Mohanty S, Mohanty BP. Proteomic changes in the liver of Channa striatus in response to high temperature stress. Electrophoresis. 2016;37(12):1704–17.

    CAS  PubMed  Article  Google Scholar 

  64. 64.

    Zhang H, Peatman E, Liu H, Niu D, Feng T, Kucuktas H, Waldbieser G, Chen L, Liu Z. Characterization of a mannose-binding lectin from channel catfish (Ictalurus punctatus). Res Vet Sci. 2012;92(3):408–13.

    CAS  PubMed  Article  Google Scholar 

  65. 65.

    Wilson JR, Williams D, Schachter H. The control of glycoprotein synthesis: N-acetylglucosamine linkage to a mannose residue as a signal for the attachment of L-fucose to the asparagine-linked N-acetylglucosamine residue of glycopeptide from alpha1-acid glycoprotein. Biochem Biophys Res Commun. 1976;72(3):909–16.

    CAS  PubMed  Article  Google Scholar 

  66. 66.

    Zhang WF, Li ZT, Fang JJ, Wang GB, Yu Y, Liu ZQ, Wu YN, Zheng SS, Cai L. Effect of mannose on the lung function of rats with acute pancreatitis. J Biol Regul Homeost Agents. 2018;32(3):627.

    CAS  PubMed  Google Scholar 

  67. 67.

    Imai K, Yoshimura T. Endocytosis of lysosomal acid phosphatase; involvement of mannose receptor and effect of lectins. Biochem Mol Biol Int. 1994;33(6):1201.

    CAS  PubMed  Google Scholar 

  68. 68.

    Makrypidi G, Damme M, Müllerloennies S, Trusch M, Schmidt B, Schlüter H, Heeren J, Lübke T, Saftig P, Braulke T. Mannose 6 Dephosphorylation of lysosomal proteins mediated by acid phosphatases Acp2 and Acp5. Mol Cell Biol. 2012;32(4):774–82.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  69. 69.

    Mateus AP, Costa R, Gisbert E, Pis P, Andree KB, Estévez A, Power DM. Thermal imprinting modifies bone homeostasis in cold challenged sea bream (Sparus aurata, L.). J Exp Biol. 2017;220(19):jeb.156174.

    Article  Google Scholar 

  70. 70.

    Adeyemi JA, Atere TG, Deaton LE. Oxidative damage and changes in glutathione S-transferase activity in juvenile African catfish, Clarias gariepinus exposed to cypermethrin and chlorpyrifos. Biokemistri. 2014;25(3):113–7.

    Google Scholar 

  71. 71.

    Armstrong RN. Glutathione S-transferases: reaction mechanism, structure, and function. Chem Res Toxicol. 1991;4(2):131.

    CAS  PubMed  Article  Google Scholar 

  72. 72.

    Jensen PJ, Gitlin JD, Carayannopoulos MO. GLUT1 deficiency links nutrient availability and apoptosis during embryonic development. J Biol Chem. 2006;281(19):13382–7.

    CAS  PubMed  Article  Google Scholar 

  73. 73.

    Ma YM, Yang MJ, Wang S, Li H, Peng XX. Liver functional metabolomics discloses an action of L-leucine against Streptococcus iniae infection in tilapias. Fish Shellfish Immunol. 2015;45(2):414–21.

    CAS  PubMed  Article  Google Scholar 

  74. 74.

    Wu C, Chen L, Lu Z, Gao J, Chu Y, Li L, Wang M, Zhang G, Zhang M, Ye J. The effects of dietary leucine on the growth performances, body composition, metabolic abilities and innate immune responses in black carp Mylopharyngodon piceus. Fish Shellfish Immunol. 2017;67:419–28.

    CAS  Article  Google Scholar 

  75. 75.

    Pancotto L, Mocelin R, Marcon M, Herrmann AP, Piato A. Anxiolytic and anti-stress effects of acute administration of acetyl-L-carnitine in zebrafish. Peerj. 2018;6:e5309.

    PubMed  PubMed Central  Article  Google Scholar 

  76. 76.

    Stein SE, Scott DR. Optimization and testing of mass spectral library search algorithms for compound identification. J Am Soc Mass Spectrom. 1994;5(9):859–66.

    CAS  PubMed  Article  Google Scholar 

  77. 77.

    Wen X, Hu YD, Zhang XY, We XZ, Wang T, Yin SW. iTRAQ-based quantitative proteomic analysis of Takifugu fasciatus liver in response to low-temperature stress. J Proteome. 2019;201:27–36.

    CAS  Article  Google Scholar 

Download references


The authors thank all members of the laboratory for valuable discussions.


This work was supported by the National Key R&D Program of china (No.2018YFD0900301), the National Natural Science Foundation of China (No. 31800436), the Natural Science Foundation of Jiangsu Province of China (No. BK20180728), the National Spark Program of China (2015GA690040), the National Finance Projects of Agro-technical popularization (TG15–003), and Project Foundation of the Academic Program Development of Jiangsu Higher Education Institution (PAPD). Our funding agencies did not play a role in the study design, data collection, analysis, interpretation of the data, or preparation of the manuscript.

Author information




SWY and XW co-conceived this study, and supervised the experiments; XZW and XYZ performed the experiments; YDH and TW conducted the data analysis and created figures and tables; XW wrote the manuscript. All the authors reviewed and approved the manuscript.

Corresponding authors

Correspondence to Tao Wang or Shaowu Yin.

Ethics declarations

Ethics approval and consent to participate

All experiments were performed in accordance with the Guidelines for the Care and Use of Laboratory Animals in China. The method of euthanasia followed the Chinese Laws for Animal Experimentation. This study was approved by the Ethics Committee of Experimental Animals at Nanjing Normal University (grant No. SYXK 2015–0028, Jiangsu). Furthermore, species Takifugu fasciatus is not listed in CITES.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. Sequencing results of six cDNA libraries. (XLSX 10 kb)

Additional file 2:

Table S2. Mapping results between transcriptome and reference genome. (XLSX 10 kb)

Additional file 3:

Table S3. The results of transcriptome annotation referring to genome. (XLSX 10878 kb)

Additional file 4:

Figure S1. Volcano plot of gene expression levels in EG (experimental group) compared with CG (control group). (PNG 122 kb)

Additional file 5:

Table S4. GO analysis of differentially expressed genes. (XLSX 191 kb)

Additional file 6:

Figure S2. GO analysis of significant differentially expressed genes. (PNG 212 kb)

Additional file 7:

Table S5. KEGG analysis of differentially expressed genes. (XLSX 16 kb)

Additional file 8:

Table S6. The classification of integrated analysis results on transcriptome and proteome. (XLSX 308 kb)

Additional file 9:

Table S7. qPCR verification results of transcriptomes. (DOCX 15 kb)

Additional file 10:

Table S8. The list of SNPs. (XLSX 19240 kb)

Additional file 11:

Table S9. Location information of the predicted SNPs. (DOCX 15 kb)

Additional file 12:

Table S10. Primers for detecting SNPs. (DOCX 15 kb)

Additional file 13:

Table S11. The results of preliminary experiment. (DOCX 15 kb)

Additional file 14:

Table S12. Primers for verifying the transcriptome. (DOCX 15 kb)

Additional file 15:

Table S13. The metabolite data set of this study. (XLSX 3761 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Wen, X., Hu, Y., Zhang, X. et al. Integrated application of multi-omics provides insights into cold stress responses in pufferfish Takifugu fasciatus. BMC Genomics 20, 563 (2019).

Download citation


  • Environmental stress
  • Multi-omics
  • Low temperature
  • mRNA-protein-metabolite
  • Takifugu fasciatus