Transcriptomic analyses on muscle tissues of Litopenaeus vannamei provide the first profile insight into the response to low temperature stress
Transcriptomic analyses on muscle tissues of Litopenaeus vannamei provide the first profile insight into the response to low temperature stress
Wen Huang 0 1
Chunhua Ren 0 1
Hongmei Li 0 1 2
Da Huo 0 1 2
Yanhong Wang 0 1
Xiao Jiang 0 1
Yushun Tian 0 1 2
Peng Luo 0 1
Ting Chen 0 1
Chaoqun Hu 0 1
0 CAS Key Laboratory of Tropical Marine Bio-resources and Ecology (LMB), South China Sea Institute of Oceanology, Chinese Academy of Sciences , Guangzhou, Guangdong , China , 2 Key Laboratory of Applied Marine Biology of Guangdong Province and Chinese Academy of Sciences (LAMB), South China Sea Institute of Oceanology, Chinese Academy of Sciences , Guangzhou, Guangdong, China, 3 South China Sea Bio-Resource Exploitation and Utilization Collaborative Innovation Center, Guangzhou, Guangdong , China
1 Editor: Sebastian D. Fugmann, Chang Gung University , TAIWAN
2 University of the Chinese Academy of Sciences , Beijing , China
The Pacific white shrimp (Litopenaeus vannamei) is an important cultured crustacean species worldwide. However, little is known about the molecular mechanism of this species involved in the response to cold stress. In this study, four separate RNA-Seq libraries of L. vannamei were generated from 13ÊC stress and control temperature. Total 29,662 of Unigenes and overall of 19,619 annotated genes were obtained. Three comparisons were carried out among the four libraries, in which 72 of the top 20% of differentially-expressed genes were obtained, 15 GO and 5 KEGG temperature-sensitive pathways were fished out. Catalytic activity (GO: 0003824) and Metabolic pathways (ko01100) were the most annotated GO and KEGG pathways in response to cold stress, respectively. In addition, Calcium, MAPK cascade, Transcription factor and Serine/threonine-protein kinase signal pathway were picked out and clustered. Serine/threonine-protein kinase signal pathway might play more important roles in cold adaptation, while other three signal pathway were not widely transcribed. Our results had summarized the differentially-expressed genes and suggested the major important signaling pathways and related genes. These findings provide the first profile insight into the molecular basis of L. vannamei response to cold stress.
Data Availability Statement: All relevant data are
within the paper and its Supporting Information
The Pacific white shrimp, Litopenaeus vannamei, is one of the most important cultural shrimp
species, and its production (by weight) has reached nearly 71% of the total economic penaeid
shrimp production worldwide [
1, 2, 3, 4
]. Naturally, L. vannamei can inhabit estuaries, lagoons
or marine areas, but they are sensitive to low temperatures . With high economic value, the
farming areas of L. vannamei have been exploded from the southern to the northern coast, and
(A201601A03), the National High Technology
Research & Development Program of China (863
Program) (2012AA10A404-4), and the Science and
Technology Service Network Initiative of the
Chinese Academy of Sciences (KFJ-SW-STS-146).
The funders had no role in study design, data
collection and analysis, decision to publish, or
preparation of the manuscript.
to the inland water in China. Therefore, low temperatures have become one of the major
constrained factors to t L. vannamei culture, since that the culturing water temperature of the
Northern China in winter is usually lower than that needed for L. vannamei survival.
Water temperature is one of the most common abiotic stress factors for aquatic ectotherms
]. It influences nearly all biochemical and physiological progresses of aquatic organisms [
and low temperatures may cause severe growth inhibition and mortality [
]. In fishes, the
understanding of biochemical and molecular mechanisms under temperature stress have been
well developed. In species such as rainbow trout (Oncorhynchus mykiss) , channel catfish
(Ictalurus punctatus) [
] and common carp (Cyprinus carpio L.) [
], studies showed that
about 10 percent of the investigated genes were found to be up-regulated under temperature
stress, and they were involved in an extensive gene expression network [
]. Digital gene
expression in Nile tilapia (Oreochromis niloticus) revealed the differentially expressed-genes
when comparing the fish cultured in 10ÊC and 30ÊC [
]. These studies provided useful
insights and promoted the molecular breeding work in fish species. However, in L. vannamei,
the gene expression profile in response to temperature stress has not yet been illustrated.
RNA-Seq technology is a crucial method for quantifying the transcriptional expression in
non-model species or in those organisms that lack whole-genome sequencing data [
Studies on channel catfish (I. punctatus) have suggested that both the de novo synthesis of
cold-induced proteins and the modification of existing proteins were required for adaptation
to low temperatures . Transcriptional regulation research in the coral reef fish Pomacentrus
moluccensis revealed the existence of differentially expressed genes in this species [
Moreover, RNA-Seq has already been applied as a useful method on other stress response studies in
L. vannamei, such as high salinity [
] and low salinity [
17, 18, 19, 20
]. However, the specific
expression profile report on L. vannamei under the stress of cold temperature is still lacked.
Molecular signaling pathways can help to improve the understanding of the activity and
coordination of the cell-expression network in response to extracellular stress. Calcium is a
ubiquitous second messenger in eukaryotic signal transduction cascades [
]. It impacts
nearly every aspect of cellular life  and fulfills its key role in transferring extracellular stress
signals into chromatin and in regulating protein phosphorylation [
]. Structural comparisons
of mammalian Ca2+/calmodulin-dependent protein kinase II (CaMKII) and plant
calciumdependent protein kinases (CDPKs) of calcium signaling showed that their kinase domain
containing the highly conserved subdomains of typical eukaryotic Ser/Thr protein kinases
]. Mitogen-activated protein kinases (MAPKs), regulated by MAPK kinase kinase-MAPK
kinase-MAPK (MKKK-MKK-MAPK) phosphorelay system, are Ser/Thr protein kinases and
are key components for vital signal transduction pathways that converting extracellular stimuli
into a wide range of cellular responses in the organisms from yeast to humans [
Consequently, in order to understand the molecular basis of L. vannamei under cold temperatures,
investigations conducted on differentially expressed genes in signaling were necessary.
Shrimps are poikilothermic species, for which changing of water temperatures can conduct
to each part of their bodies and influence their cellular activities. Muscle compose nearly 60%
body weight of L. vannamei, and it is the main tissue for converting biochemical energy in
]. Furthermore, we found that the muscle tissues of L. vannamei turned whitish
immediately when they were chilled in low temperature water, but regained pellucidity after
being returned to the control temperature, indicating that the muscle tissues of L. vannamei
may process overt physiology activities in response to acute cold stress. In the present study, a
transcriptomic analysis has been carried out on the muscle of L. vannamei. Gene expression
profiles under the stress of low temperature has been characterized by comparing various
libraries from corresponding temperature treatments, including acute cold stress for 2 hours,
endurance cold stress for 48 hours and recovering to normal temperature for 2 hours. The
2 / 21
identified candidate genes will be useful and helpful in revealing a deeper insight into the
molecular adaptation basis under low temperature stress in L. vannamei.
Materials and methods
Shrimp care and experiments were carried out according to the Care and Use of Agricultural
Animals in Agricultural Research and Teaching and approved by the Science and Technology
Bureau of China. Approval from the Department of Wildlife Administration was not required
for the experiments conducted in this paper. All experiments in this paper were performed
with permits obtained from the Government of the People's Republic of China and endorsed
by the Animal Experimentation Ethics Committee of Chinese Academy of Sciences.
Treatments and cold challenges of the shrimp materials
Healthy shrimp with an average weight of 7.09±3.22 g were treated in gradient temperatures
(9, 10, 11, 12, 13, 14, 15, 18, 21, 24, 27, 30 and 33ÊC) for 2 hours, 48 hours and then recovered
at a control temperature (28ÊC) for 2 hours. Survival was evaluated by observing the shrimp's
ability to move spontaneously or after gentle prodding [
]. The temperature of semi-death
group was set as the stress cold temperature.
Thirteen degree was set as the stress cold temperature. Four groups were placed either in
control temperature (CT), 13ÊC for 2 hours (T_2), 13ÊC for 48 hours (T_48) and then, the
T_48 group was placed to recover temperature for 2 hours (RC). Muscle tissues from the
surviving individuals were put in centrifuge tubes and dipped into liquid nitrogen immediately.
Three biological replicates in equal amounts of each group (2 μg of RNA for each individual,
and total 6 μg of RNA for each group) were pooled for subsequent RNA-Seq experiments.
RNA isolation and cDNA library construction
Four separate libraries (CT, T_2, T_48 and RC) were sampled. Total RNA was extracted by
TRIzol Reagent (Invitrogen, USA) according to the manufacturer's instructions. The quality of
the RNA products was verified as described by Ren et al. (2014) [
]; RNase-free DNase I
(Invitrogen, USA) was used to degrade any possible DNA, and oligo (dT)-coated magnetic
beads were mixed with the total RNA to concentrate the ployA tailed mRNA. Fragmentation
was performed by incubation in an NEB Next First Strand Synthesis Reaction Buffer (NEB,
USA), and the second strand was generated with the buffer, dNTPs, RNase H and DNA
polymerase-I. Adapters were ligated to the synthesized cDNA fragments after an end repair step.
RNA-Seq and annotation
The library was constructed using an Illumina HiSeq 4000 sequencing platform located at the
BGI Company (Shenzhen, China; http://www.genomics.cn/index). The clean read data were
deposited to the US National Center for Biotechnology Information (NCBI) Sequence Read
Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra) with the accession No. SRP095377.
After assembling, functional annotation by gene ontology terms (GO, http://www.
geneontology.org) was analyzed using the Blast2GO program [
]. The Kyoto Encyclopedia of
Genes and Genomes Pathway (KEGG; http://www.genome.jp/kegg), the Non-redundant
NCBI collection of nucleotide (NT) and protein (NR) sequence database (ftp://ftp.ncbi.nih.
gov/blast/db/FASTA/), the Swiss-Prot database (http://web.expasy.org/docs/swiss-prot_
guideline.html) and the Cluster of Orthologous Groups (COG) protein database (http://www.
ncbi.nlm.nih.gov/COG/) were used to predict and classify the possible functions [
3 / 21
Quantification of differentially expressed genes
Three comparisons were taken to the four libraries, including CT-VS-T_2 (T_2/CT),
T_2-VS-T_48 (T_48/T_2), and T_48-VS-RC (RC/T_48). The differentially expressed genes
(DEGs) in comparison of T_2/CT represented the transcripts in response to acute cold stress
(placed from control temperature to acute 13ÊC). The DEGs in comparison of T_48/T_2
represented the transcripts in response to endurance cold stress (endured the cold stress for 48
hours). The DEGs in comparison of RC/T_48 represented the transcripts for recover
adaptation (placed from stress cold to control temperature). The clean reads per kilobase per million
(RPKM) value was used to estimate transcript abundance on the basis of eliminating the
influence of different gene lengths and sequencing discrepancies [
]. Therefore, the RPKM values
could be directly used to compare the differences of gene expression among samples. The
algorithm developed by Audic and Claverie (1997) [
] was used to compare the differences in
various gene expression levels. A value of the FDR (false discovery rate) of less than 0.001 and a
value of the Log2 ratio of less than 1 were set as the criteria values.
Identification of major response genes and pathways
The coexisted DEGs in different comparisons of total genes, GO annotated genes and KEGG
annotated genes were filtered. Four important kinds of signaling genes regarding Calcium and
MAPK cascades, Transcription factor and Serine/threonine-protein kinase signaling genes
were searched from the total annotated genes [
]. To understand the expression model of
target genes, hierarchical clustering analysis was performed with the software Heatmap Illustrator
Real-time PCR validation
To verify the expression profiles of genes in our RNA-Seq results, the represented 15 DEGs in
comparison of T_2/CT were selected for validation of the Illumina sequences by real-time
PCR analysis. All primers used were listed in Table 1. The treatments of L. vannamei, sampling
of muscle, isolation of total RNA and treatment of DNase were performed as described above.
The first cDNA strand was synthesized by PrimeScript1 Reverse Transcriptase Kit (Takara,
China) with an oligo (dT) primer, a qPCR analysis was conducted in a RotorGene RG3000
real-time PCR system (Corbett research, Australia), PCR reactions were conducted using a
SYBR Premix Ex TaqTM II (Takara, China), and the reactions were exposed to an initial
denaturation (94ÊC for 3 min) followed by 40 cycles of 94ÊC for 25 s, 60ÊC for 15 s, and 72ÊC for
15 s. The relative transcript abundance was obtained by normalizing with expression of the L.
vannamei β-actin gene based on the 2-ΔΔCT method [
Cold stress temperature
A batch of L. vannamei were treated with different temperatures, and the survival rates were
summarized (Fig 1B). The results showed that L. vannamei could not survive when treated
with 11ÊC or lower temperatures, the survival rates were slightly affected in the water at 14ÊC,
15ÊC or higher temperatures (Fig 1B). The survival rates at 12ÊC and 13ÊC were 9.30% and
47.22%, respectively. We used 13ÊC as the cold stress temperature since a nearly half survival
rate (47.22%) of shrimp was observed at this temperature. The same batch of shrimp from the
same cultural environment was used for the following experiments.
4 / 21
Overview of the libraries
After filtering, the quality metrics of the reads are showed in Table 2: the total raw reads of the
samples were from 87.69 to 90.94 Mb, the number of clean reads per library ranged from 65.23
to 66.02 Mb, and the clean reads ratio ranged from 71.95% to 75.28%. The mean value of clean
reads Q20 and Q30 percentages were 97.41% and 94.34%, respectively. The clean reads were
assembled and clustered into Unigenes, and the quality metrics of Unigenes are presented in
Table 3. The total number of All-Unigenes was 29,662, the total length was 29,632,338 bp and
the mean length was 999 bp. The annotated genes in NR, NT, SwissProt, KEGG, COG, GO
database were 17,103 (57.66%), 11,976 (40.37%), 15,552 (52.43%), 13,936 (46.98%), 8,558
(28.85%) and 2,504 (8.44%), respectively (Fig 2A). Overall 19,619 Unigenes were annotated
(Fig 2A). The main Venn diagram of NR, KEGG, COG and SwissProt annotation are showed
in Fig 2B, the largest interior gene numbers were 8,264. The percentage of annotated species
with NR database are displayed in Fig 2C, the top three annotated species were Zootermopsis
nevadensis (the largest blue part, 12.34%), Daphnia pulex (the second large tenne part, 5.31%)
and Tribolium castaneum (the third large gray part, 2.95%), respectively.
5 / 21
Fig 1. The cold challenge and survival rate of L. vannamei. (A): Treatment processes of the shrimp. (B):
Gradient temperature challenge and survival rates. Control: treated in control temperature of 28℃; T_2:
treated in related temperature for 2 hours; T_48: treated in related temperature for 48 hours; Recover: treated
in control temperature of 28℃ for 2 hours recovering. Red frame indicated the stress temperature.
GO and KEGG classification of the DEGs
GO annotation was conducted for the three comparisons. In the comparison group of T_2/
CT, 870 of 1471 DEGs could be assigned by GO classification, and the equivalent numbers for
other comparison groups were 1147 of 1943 DEGs (T_48/T_2) and 690 of 1477 DEGs (RC/
T_48); the detailed GO classification data (Pathway, PathwayID, and GeneID) are listed in S1
Table. The most annotated GO pathway in each comparison was Catalytic activity (GO:
KEGG annotation was conducted for the three comparisons. In the comparison of T_2/CT,
1202 of 1471 DEGs were assigned to 83 pathways. In the comparison of T_48/T_2, 1556 of
aQ20: the rate of bases whitch quality is greater than 20
b Q30: The rate of bases which quality is greater than 30
6 / 21
a N50: a weighted median statistic that 50% of the TotalLength is contained in Unigenes great than or equal to this value
b N70: a weighted median statistic that 70% of the TotalLength is contained in Unigenes great than or equal to this value
c N90: a weighted median statistic that 90% of the TotalLength is contained in Unigenes great than or equal to this value.
d GC (%): the percentage of G and C bases in all Unigenes
1943 DEGs were assigned to 91 pathways. In the comparison of RC/T_48, 1208 of 1477 DEGs
were assigned to 77 pathways. The details of the KEGG classification (Pathway, PathwayID,
and GeneID) of the three comparisons are listed in S2 Table. The top two annotated KEGG
pathways in all the three comparisons were Metabolic pathways (ko01100) and RNA transport
pathways (ko03013). The RNA transport pathway maps (there were no Metabolic pathway
maps in KEGG database) were traced and downloaded from the KEGG database (S1 Fig), the
annotated DEGs in the three comparisons were summarized (S3 Table), the major discrepant
DEGs in RNA transport pathway were the components of exon-junction complex (EJC). The
annotated up- and down-regulated genes of EJC complex were marked and displayed in Fig 3.
The EJC complex mainly participated in mRNA translation (S1 Fig, Fig 3A), the EJC outer
Fig 2. Annotation of the transcriptome. (A): The summary results with main annotation database; the X
axis were labeled with various items; Y axis represented the proportion genes annotated by the functional
database; 29,662 was the total numbers of Unigenes, only the percentage numbers were marked in other
blue dots. (B): The Venn diagram between NR, KEGG, COG and Swissprot; the annotated gene numbers
were exhibit in the related areas. (C): The pie chart of percentage of annotated species with NR database; the
main species names were displayed under the pie chart, species (with the corresponding color) were
displayed with descending order.
7 / 21
Fig 3. The main components and DEGs in RNA transport pathways. (A): mRNA translation processes. (B): the annotated
different expressed genes in EJC complex of T_2/CT. (C): the annotated different expressed genes in EJC complex of T_48/T_2.
(D): the annotated different expressed genes in EJC complex of RC/T_48. The up-regulated genes were marked with red frame,
down-regulated genes were marked with green frame.
shell and Transiently interacting factor genes were the most two kinds of responding
components in EJC complex (S3 Table, Fig 3B, 3C and 3D).
Represented GO and KEGG classification columnar figures in different comparisons are
showed in Fig 4. GO classes contained three categories, including Biological processes, Cellular
components and Molecular functions. The top two annotated classes of Biological processes in
the comparison of T_2/CT were Cellular processes (97 annotated genes) and Metabolic
processes (96 annotated genes), the major classes of Cellular components were Cell (57 annotated
genes) and Cell parts (57 annotated genes), and the major classes of Molecular functions were
Binding (94 annotated genes) and Catalytic activity (117 annotated genes). The top two
annotated KEGG classes in comparison of T_2/CT were Global and overview maps (107 annotated
genes of Metabolism) and Translation (76 annotated genes of Genetic Information
Processing). The most annotated GO classes in comparison of T_48/T_2 were Catalytic activity (151
annotated genes) in Molecular function, Cell (65 annotated genes) and Cell junction (65
annotated genes) in Cellular component, Metabolic process (128 annotated genes) in Biological
process; the most top two annotated KEGG classes were Global and overview maps (123
annotated genes of Metabolism) and Translation (106 annotated genes of Genetic Information
Processing). The most annotated GO classes in comparison of RC/T_48 were Catalytic activity (96
8 / 21
Fig 4. GO and KEGG classification of the DEGs in each comparison. (A): Comparison of T_2/CT. (B): Comparison of T_48/
T_2. (C): Comparison of RC/T_48. The Y-axis represents the number of DEGs in each item.
9 / 21
annotated genes) in Molecular function, Cell (36 annotated genes) and Cell part (36 annotated
genes) in Cellular component, Cellular process (84 annotated genes) in Biological process; the
most top two KEGG classes were Global and overview maps (90 annotated genes of
Metabolism) and Translation (82 annotated genes of Genetic Information Processing).
Quantification of transcripts
The transcript abundance of each Unigene from the four related libraries is listed in S4 Table.
DEGs were identified through a pairwise comparison between various transcripts by setting a
threshold FDR (false discovery rate) value of 0.001 and a |log2FoldChange| value of 1 based on an
algorithm described by Audic and Claverie (1997) and Ren et al. (2014). From the three
comparisons (T_2/CT, T_48/T_2 and RC/T_48), a large number of DEGs were identified, and the
filtered results (P value < 0.05, FDR
0.001, and estimated absolute |log2FoldChange|
exhibited in S5±S7 Tables. The columnar and Volcano-plot results of the DEGs are showed in
Fig 5. The X-axis of Volcano-plot represented the comparison value of log2FoldChange, the
YFig 5. Column and Volcano-plot of the DEGs in comparisons. The columnar results of the DEGs are listed
and the corresponding Volcano-plot figures are shown. Plots in the Volcano-plot figures represent the DEGs.
The Y-axis in the columnar figures represents the DEG numbers in comparisons. In the Volcano-plot figures,
the X-axis represents the comparison value of log2FoldChange, and the Y-axis represents the value of |-log10FDR|.
The red color refers to up-regulated DEGs, blue color refers to down-regulated DEGs, and the black plot refers
to no significantly different DEGs between the pairs of libraries (FDR value 0.001, |log2FoldChange| value 1).
10 / 21
axis represented value of |-log10FDR|. From the pictures, CT-vs-T_2 (T_2/CT) had 904 up- and
567 down-regulated genes, T_2-vs-T_48 (T_48/T_2) had 634 up- and 1309 down-regulated
genes and T_48-vs-RC (RC/T_48) had 997 up- and 480 down-regulated genes. The
down-regulated transcripts were far more than the up-regulated genes in comparison of T_48/T_2,
indicating that more genes were suppressed under the endurance cold stress.
Identification of major response genes
The Venn diagrams in Fig 6A were used to exhibit the coexisted DEGs in different
comparisons. It showed that 234 genes were coexisting in the three comparisons of Total DEGs, and 35
of coexisted GO annotated DEGs and 48 of coexisted KEGG annotated DEGs were also
Top 20% of differently-expressed coexisted genes in the three comparisons were selected
and taken to Hierarchical clustering analysis, 72 genes were obtained, the detailed GeneID and
other information are showed in S8 Table, the heat map figures are showed in Fig 6B. Take the
Unigene of CL532.Contig1_All as an example, significantly up-regulated (log2FoldChange
value = 9.33) in comparison of T_2/CT, down-regulated (log2FoldChange value = -9.33) in
comparison of T_48/T_2 and up-regulated (log2FoldChange value = 8.60) in comparison of RC/T_48.
The results indicated that this gene was sensitive to acute cold temperature and recover
temperature, and it was down-regulated under the endurance cold stress environment. This gene
might play important role in adaptation of changing temperature but was greatly inhibited in
long time cold stress environment. The Unigene of Unigene3735_All was taken as another
example, for which was significantly down-regulated (Log2FoldChange Value = -7.63) in
comparison of T_2/CT, up-regulated (Log2FoldChange Value = 8.12) in comparison of T_48/T_2 and
down-regulated (log2FoldChange Value = -8.12) in comparison of RC/T_48. The results indicated
this gene was transcriptional inhibited under the environment of changing temperature but
was highly expressed under the environment of endurance cold stress, suggesting that it might
play important role in response to long time cold stress.
Identification of major pathways
Total 35 coexisted GO-annotated DEGs and 48 coexisted KEGG-annotated DEGs were
selected, the detailed GeneID information are showed in S9 and S10 Tables. The major
(Absolute Value of log2FoldChange 5.0) Unigene in GO and KEGG Pathways were filtered, and 15
GO pathways and 5 KEGG pathways were recognized as the temperature-sensitive pathways
in response to cold stress (Table 4). The Unigene of CL596.Contig2_All, described as
Eukaryotic translation initiation factor 5 gene (gi|646705355|), was annotated by Translation
initiation factor 5 pathway (KEGG: K03262) and Cellular macromolecule metabolic (GO: 0044260)
and primary metabolic process (GO: 0044238), implying that this gene might play important
role in several molecular functions in response to cold stress.
Major signaling genes in response to cold stress
Four important kinds of signaling genes regarding Calcium signaling, MAPK cascades,
Transcription factors and Serine/threonine-protein kinase were fished out (S11±S14 Tables) and
clustered into heat map. Generally, the Calcium signaling genes, MAPK cascades genes and
Transcription factor genes were not widely transcribed in response to cold stress (Fig 6C, 6D
and 6E). Four Unigenes of Calcium signaling genes, 4 Unigenes of MAPK cascades genes and
5 Unigenes of Transcription factor genes were obtained when filtered by the top two levels of
log2FoldChange (absolute value of log2FoldChange 6.0 in Calcium signaling genes, absolute value
of log2FoldChange 5.0 in MAPK cascades genes and absolute value of log2FoldChange 3.5 in
11 / 21
Fig 6. Venn diagrams and heat maps of the major coexisted DEGs. (A): The Venn diagrams of Total DEGs, GO annotated
DEGs and KEGG annotated DEGs in the three comparisons; DEG numbers were labeled in the corresponding areas. (B):
12 / 21
Hierarchical clustering analysis of the top 20% of Total coexisting DEGs in the three comparisons. (C): Hierarchical clustering
analysis of the annotated Calcium signaling genes in cold stress. (D): Hierarchical clustering analysis of the annotated MAPK
cascades genes in cold stress. (E): Hierarchical clustering analysis of the annotated Transcription factor genes in cold stress. (F):
Hierarchical clustering analysis of the annotated Serine/threonine-protein kinase genes in cold stress.
Transcription factor genes, Table 5). Whereas, conspicuous patches of
Serine/threonine-protein kinase genes were observed (Fig 6F), 16 Unigenes were obtained when filtered by the top
two levels of log2FoldChange (absolute value of log2FoldChange 5.0). The results suggested that
Serine/threonine-protein kinase might participate in the adaptation of acute or endurance
Real-time PCR validation
The Real-time PCR outcomes were further summarized. The qPCR results of the 15 genes in
the treatments of CT and T_2 are showed in Fig 7A. Grouped comparison results between
aT_2/CT: The value of Log2FoldChange in comparison of T_2/CT
bT_48/T_2: The value of Log2FoldChange in comparison of T_48/T_2
cRC/T_48: The value of Log2FoldChange in comparison of RC/T_48
13 / 21
aT_2/CT: The value of Log2FoldChange in comparison of T_2/CT
b T_48/T_2: The value of Log2FoldChange in comparison of T_48/T_2. Ca: Calcium signaling pathway genes Mapk: Mitogen-activated protein kinases
cascades genes; Ser/Thr: Serine/threonine-protein kinase genes; TFs: Transcription factor genes
qPCR and RNA-Seq were also displayed, and the results showed all the 15 candidate genes in
qPCR verification consisted with the results of the RNA-Seq technology (Fig 7B).
14 / 21
Fig 7. Validation for the 15 transcriptomic DEGs in the comparison of T_2/CT. (A): The Real-time PCR results of the 15 genes in
treatments of CT and T_2; Y axis represented the Log2 value of the related expression. (B): Grouped results of the 15 genes; the blue
columnar mean the Real-time PCR results, red columnar mean the RNA-Seq results; X axis indicated the gene names, Y axis represented
the value of log2FoldChange value in the comparison of T_2/CT, the FoldChange value were obtained by2±ΔΔCT method. Bars in this figure
represented the standard deviation (SD).
Temperature is one of the most important environmental variances for aquatic ectotherms,
and it has been investigated for its key role in affecting overall biological processes in various
]. In larval zebrafish (Danio rerio), cold stress induced-transcripts were generated
from 2 and 48 hours treatments under 16ÊC . In tilapia (O. niloticus), an analysis of gene
expressions was carried out under cold stress by comparing different temperatures, from 30ÊC
followed by a decrement of 1ÊC/Day to 10ÊC [
]. In Manila clam (Ruditapes philippinarum),
a transcriptome study was conducted in response to cold stress with an extreme low
temperature of -1ÊC [
]. Stresses could be classified into mild, chronic or acute, which represented a
15 / 21
dramatic shift in the environmental conditions [
]. In our study, shrimp were treated in
acute cold temperatures for 2 hours, 48 hours and then in a recovery temperature for 2 hours,
and a nearly 50% survival rate (47.22%) was observed in the 13ÊC group. The stress
temperature in our result was similar to that reported by Fan et al. (2013) [
], in which 13ÊC for 24
hours was regarded as a suitable cold stress temperature to investigate the altered proteins and
genes in L. vannamei.
RNA-Seq technology has been demonstrated as a straightforward method to reveal the
general molecular regulation processes in species underlying the response to stressful
]. In the previous studies for L. vannamei, RNA-Seq technology has been employed
to investigate the different gene expression levels involved in various environment factors and
developmental processes, such as a transcriptome analysis under acute ammonia stress [
some digital gene expression analyses in response to low salinity stress [
17, 18, 19, 20, 39
the transcriptional profiling during ontogenesis . Parts of the molecular mechanisms have
already been characterized in response to temperature in the shrimp. The protein content of
CAT, GST, Ferritin and HSP60 have been found to be influenced by temperature stress in L.
]. Lipid peroxidation, catalase activity and glutathione-S-transferase in Palaemon
elegans and Palaemon serratus have been proved as the biomarkers for monitoring aquatic
]. In addition, the decreasing water temperatures has been shown to induce
oxidative stress in L. vannamei [
]. Furthermore, 20 proteomic spots and several genes
have demonstrated to play important roles in temperature regulation of L. vannamei .
However, the profile information of the regulation of gene expression under low temperatures
in L. vannamei was still lacked. In the present study, transcriptomic analyses were conducted
under cold temperature stress on the muscle of L. vannamei. To our knowledge, this is the first
report of transcriptional profiling under acute low temperature stress in L. vannamei.
Transcriptome analyses may provide insights into the molecular basis of cold tolerance in
many species. Transcriptomic characterization under temperature stress in zebrafish (D. rerio)
revealed the transcriptional and post-transcriptional regulation of the cold-acclimated
], the multiple temperature-regulated biological processes and pathways [
], and the
over-expressed biological processes (such as circadian rhythm, energy metabolism, lipid
transport and metabolism) . An RNA-Seq analysis in Melanotaenia duboulayi provided
candidate genes in response to changing temperatures in freshwater fish [
]. The study of
transcriptomic analysis on the marine fish ectoparasitic ciliate Cryptocaryon irritans under a
low temperature revealed the required genes for cell survival and a deeper dormancy state
]. In our present study, a transcriptomic analysis for L. vannamei under acute low
temperature was carried out. In total of the four libraries, 904 up- and 567 down-regulated genes in
comparison of T_2/CT, 634 up- and 1309 down-regulated genes in comparisons of T_48/T_2,
and 997 up- and 480 down-regulated genes in comparisons of RC/T_48 were detected. The
top 20% of coexisted DEGs were fished out. In our result, 15 GO pathways and 5 KEGG
pathways were annotated as temperature-sensitive pathways. Our study promoted the
understanding of profiling gene expression of L. vannamei in response to cold stress.
Acute stress places cells in danger, and rapid adaptation is crucial for maximizing cell
]. To adapt to either biotic or abiotic stress, cellular signaling pathways may
coordinate multiple integrated and regulated layers of different steps of mRNA biogenesis in inner
]. Thus, the important signaling pathways regarding cellular stress response
adaptation have been investigated in many studies. Calcium ions (Ca2+) impact nearly every aspect
of cellular life ; Teets et al. (2013) demonstrated that goldenrod gall fly (Eurosta solidaginis)
tissues used calcium signaling to instantly detect decreases in temperature and to trigger
downstream cold-hardening mechanisms [
]. In the transcriptomic analysis reported by Ren
et al. (2014), the gene members of the calcium signaling, including cation/calcium exchangers,
16 / 21
calcium-binding proteins, calmodulin-like proteins, etc., were found to be involved in the
response to cold stress in Chrysanthemum morifolium [
]. The MAPKs influenced the cellular
signal transduction in organisms from yeast to human. When yeast under osmostress, the
Ste11 MAP3Ks were activated, followed by activation of Pbs2 MAP2K that combined to Hog1
and then phosphorylation of the specific osmostress transcription factors for cell adaption
]. The p38 family in mammalian species contains 4 isoforms (MAPK11-14) with many
overlapping functions in cells, and downstream targets of p38 MAPKs include several kinases
that are involved in the control of gene expression and nuclear proteins, such as transcription
factors (CREB, ATF1, NFκB, STAT1 and STAT3) and regulators of chromatin remodeling
(H3 and HMG14) [
]. In this paper, the expression model for the genes in four signal
pathways were investigated under cold stress. The major differentially-expressed signaling genes
were suggested involve in the adaptation of cold stress. Serine/threonine-protein kinase genes
might play more important roles in response to cold stress, while Calcium, MAPK cascades
and Transcription factor signaling genes were not widely transcribed in response to cold
In conclusion, the overall transcriptional expression and quantification analyses under the
cold stress temperature of 13ÊC were conducted for Litopenaeus vannamei. Top 20% of
differentially expressed genes were filtered. The most annotated and temperature-sensitive of GO
and KEGG pathways were found out. The major differentially-expressed genes in signaling
were suggested. The results showed that Ser/Thr kinase signal pathway might play more
important roles in participating the cold adaptation, while Calcium, MAPK cascades and
Transcription factors signaling genes were not widely transcribed. Our study provided the first
insight into the molecular basis and supplied an important reference for the mechanism of L.
vannamei responded to cold stress.
S1 Fig. The RNA transport pathway map from KEGG database (map03013).
S1 Table. GO classification of DEGs in each comparison.
S2 Table. KEGG pathway classification of the DEGs in each comparison.
S3 Table. Summary of DEGs of the main components of RNA transport pathway.
S4 Table. The transcription level of each unigene derived from the number of relevant
reads of the four libraries. The ªGeneLengthº column gives the length of exon sequence. "-"
represented absence of this unigene in the related group of transcript. RPKM: clean reads per
kilobase per million for each unigene in each libraries, respectively.
S5 Table. Different transcribed genes between the comparison of T_2/CT. a The
ªGeneLengthº column gives the length of exon sequence. b CT: reads per kb per million reads
(RPKM) for each unigene in transcription of Control. c T_2: reads per kb per million reads
(RPKM) for each unigene in transcription of T_2. d Log2FoldChange(T_2/CT): the ratio
between the RPKM in CT to the RPKM in T_2. The criteria applied for assigning significance
were: Pvalue < 0.05, FDR 0.001, and estimated absolute |log2FoldChange| 1. Genes listed
in descending order of absolute |log2FoldChange|. NA: No annotation were found in the
17 / 21
S6 Table. Different transcribed genes between the comparison of T_48/T_2. a The
ªGeneLengthº column gives the length of exon sequence. b T_2: reads per kb per million reads
(RPKM) for each unigene in transcription of T_2. c T_48: reads per kb per million reads
(RPKM) for each unigene in transcription of T_48. d Log2FoldChange(T_48/T_2): the ratio
between the RPKM in T_2 to the RPKM in T_48. The criteria applied for assigning
significance were: Pvalue < 0.05, FDR 0.001, and estimated absolute |log2FoldChange| 1. Genes
listed in descending order of absolute |log2FoldChange|. NA: No annotation were found in the
S7 Table. Different transcribed genes between the comparison of RC/T_48. a The
ªGeneLengthº column gives the length of exon sequence. b T_48: reads per kb per million reads
(RPKM) for each unigene in transcription of T_48. c RC: reads per kb per million reads
(RPKM) for each unigene in transcription of RC. d Log2FoldChange(RC/T_48): the ratio
between the RPKM in T_48 to the RPKM in RC. The criteria applied for assigning significance
were: Pvalue < 0.05, FDR 0.001, and estimated absolute |log2FoldChange| 1. Genes listed
in descending order of absolute |log2FoldChange|. NA: No annotation were found in the
S8 Table. Top 20% of coexisted DEGs in the three comparisons. a T_2/CT: The value of
Log2FoldChange in comparison of T_2/CT; b T_48/T_2: The value of Log2FoldChange in
comparison of T_48/T_2; c RC/T_48: The value of Log2FoldChange in comparison of RC/
S9 Table. The coexisted DEGs in GO annotation. a T_2/CT: The value of Log2FoldChange
in comparison of T_2/CT; b T_48/T_2: The value of Log2FoldChange in comparison of T_48/
T_2; c RC/T_48: The value of Log2FoldChange in comparison of RC/T_48.
S10 Table. The coexisted DEGs of KEGG annotation. a T_2/CT: The value of
Log2FoldChange in comparison of T_2/CT; b T_48/T_2: The value of Log2FoldChange in comparison
of T_48/T_2; c RC/T_48: The value of Log2FoldChange in comparison of RC/T_48.
S11 Table. The differential gene expression of Calcium signaling genes in each
S12 Table. The differential gene expression of MAPK cascades genes in each comparison.
S13 Table. The differential gene expression of Transcription factors (TFs) in each
S14 Table. The differential gene expression of Serine/threonine-protein kinase genes in
18 / 21
The authors are grateful to Lin Huang (Boya College, Sun Yat-Sen University, Guangzhou
510275, China) for helping us with English Copy-editing.
Conceptualization: WH TC CH.
Data curation: WH YW XJ PL TC.
Formal analysis: WH CR TC.
Funding acquisition: WH RC TC CH.
Investigation: WH HL DH YT.
Methodology: WH CR TC.
Project administration: CR CH.
Resources: CR CH.
Software: WH YW XJ PL TC.
Supervision: CR CH.
Validation: WH RC YW XJ PL TC.
Visualization: WH TC.
Writing ± original draft: WH.
Writing ± review & editing: RC TC CH.
19 / 21
20 / 21
1. FAO. Fisheries and Aquaculture Department: The State of World Fisheries and Aquaculture . Rome, Italy: Food and Agriculture Organization of the United Nations; 2012 : 37 .
2. Kim S , Pang Z , Seo H , Cho Y , Samocha T , Jang I . Effect of bioflocs on growth and immune activity of Pacific white shrimp, Litopenaeus vannamei postlarvae . Aquac Res . 2014 ; 45 : 362 ± 371 .
3. Chen T , Ren C , Wang Y , Luo P , Jiang X , Huang W , et al. Molecular cloning, inducible expression and antibacterial analysis of a novel i-type lysozyme (lyz-i2) in Pacific white shrimp, Litopenaeus vannamei . Fish Shellfish Immun . 2016 ; 54 : 197 ± 203 .
4. Ghaffari N , Sanchez-Flores A , Doan R , Garcia-Orozco K , Chen P , Ochoa-Leyva A , et al. Novel transcriptome assembly and improved annotation of the whiteleg shrimp (Litopenaeus vannamei), a dominant crustacean in global seafood mariculture . Sci Rep . 2014 ; 4 : 7081 . https://doi.org/10.1038/ srep07081 PMID: 25420880
5. Long Y , Song G , Yan J , He X , Li Q , Cui Z. Transcriptomic characterization of cold acclimation in larval zebrafish . BMC Genomics . 2013 ; 14 : 612 . https://doi.org/10.1186/ 1471 -2164-14-612 PMID: 24024969
6. Wang Q , Tan X , Jiao S , You F , Zhang P. Analyzing Cold Tolerance Mechanism in Transgenic Zebrafish (Danio rerio) . PLoS ONE . 2014 ; 9 ( 7 ):e102492. https://doi.org/10.1371/journal.pone. 0102492 PMID: 25058652
7. Long Y , Yan J , Song G , Li X , Li X , Li Q , et al. Transcriptional events co-regulated by hypoxia and cold stresses in Zebrafish larvae . BMC Genomics . 2015 ; 16 : 385 . https://doi.org/10.1186/s12864-015-1560- y PMID: 25975375
8. Vornanen M , Hassinen M , Koskinen H , Krasnov A . Steady-state effects of temperature acclimation on the transcriptome of the rainbow trout heart . Am J Physiol Regul Integr Comp Physiol . 2005 ; 289 : R1177 ± 1184 . https://doi.org/10.1152/ajpregu.00157. 2005 PMID: 15932967
9. Ju Z , Dunham RA , Liu Z. Differential gene expression in the brain of channel catfish (Ictalurus punctatus) in response to cold acclimation . Mol Genet Genomics . 2002 ; 268 : 87 ± 95 . https://doi.org/10.1007/ s00438-002 -0727-9 PMID: 12242503
10. Gracey A , Fraser E , Li W , Fang Y , Taylor R , Rogers J , et al. Coping with cold: An integrative, multitissue analysis of the transcriptome of a poikilothermic vertebrate . PNAS . 2004 ; 101 : 16970 ± 16975 . https:// doi.org/10.1073/pnas.0403627101 PMID: 15550548
11. Long Y , Li L , Li Q , He X , Cui Z. Transcriptomic Characterization of Temperature Stress Responses in Larval Zebrafish . PLoS ONE . 2012 ; 7 ( 5 ):e37209. https://doi.org/10.1371/journal.pone. 0037209 PMID: 22666345
12. Yang C , Jiang M , Wen H , Tian J , Liu W , Wu F , et al. Analysis of differential gene expression under lowtemperature stress in Nile tilapia (Oreochromis niloticus) using digital gene expression . Gene . 2015 ; 564 : 134 ± 140 . https://doi.org/10.1016/j.gene. 2015 . 01 .038 PMID: 25617524
13. Hornett E , Wheat C . Quantitative RNA-Seq analysis in non-model species: assessing transcriptome assemblies as a scaffold and the utility of evolutionary divergent genomic reference species . BMC Genomics . 2012 ; 13 ( 1 ): 361 .
14. Ren L , Sun J , Chen S , Gao J , Dong B , Liu Y , et al. A transcriptomic analysis of Chrysanthemum nankingense provides insights into the basis of low temperature tolerance . BMC Genomics . 2014 ; 15 : 844 . https://doi.org/10.1186/ 1471 -2164-15-844 PMID: 25277256
15. Kassahn K , Crozier R , Ward A , Stone G , Caley M. From transcriptome to biological function: environmental stress in an ectothermic vertebrate, the coral reef fish Pomacentrus moluccensis . BMC Genomics . 2007 ; 8 : 358 . https://doi.org/10.1186/ 1471 -2164-8-358 PMID: 17916261
16. Kumlu M , Kumlu M , Turkmen S . Combined effects of temperature and salinity on critical thermal minima of pacific white shrimp Litopenaeus vannamei (Crustacea: Penaeidae) . J Therm Biol . 2010 ; 35 : 302 ± 304 .
17. Gao W , Tan B , Mai K , Chi S , Liu H , Dong X , et al. Profiling of differentially expressed genes in hepatopancreas of white shrimp (Litopenaeus vannamei) exposed to long-term low salinity stress . Aquaculture . 2012 ; 364 ( 365 ): 186 ± 191 .
18. Chen K , Li E , Li T , Xu C , Wang X , Lin H , et al. Transcriptome and Molecular Pathway Analysis of the Hepatopancreas in the Pacific White Shrimp Litopenaeus vannamei under Chronic Low-Salinity Stress . PLoS ONE . 2015 ; 10 ( 7 ):e0131503. https://doi.org/10.1371/journal.pone.0131503 PMID: 26147449 19 . Wang X , Wang S , Li C , Chen K , Qin J , Chen L , et al. Molecular pathway and gene responses of the pacific white shrimp Litopenaeus vannamei to acute low salinity stress . J Shellfish Res . 2015 ; 34 ( 3 ): 1037 ± 1048 .
20. Zhao Q , Pan L , Ren Q , Hu D. Digital gene expression analysis in hemocytes of the white shrimp Litopenaeus vannamei in response to low salinity stress . Fish Shellfish Immun 2015 ; 42 : 400 ± 407 .
21. Berridge M , Lipp P , Bootman M. The versatility and universality of calcium signaling . Nat Rev Mol Cell Biol . 2000 ; 1 ( 1 ): 11 ± 21 . https://doi.org/10.1038/35036035 PMID: 11413485
22. Cheng S , Willmann M , Chen H , Sheen J. Calcium Signaling through Protein Kinases. The Arabidopsis Calcium-Dependent Protein Kinase Gene Family . Plant Physiol . 2002 ; 129 : 469 ± 485 . https://doi.org/ 10.1104/pp.005645 PMID: 12068094
23. Clapham D. Calcium Signaling . Cell . 2007 ; 131 ( 14 ): 1047 ± 1058 .
24. Dodd A , Kudla J , Sanders D. The Language of Calcium Signaling . Annu Rev Plant Biol . 2010 ; 61 : 593 ± 620 . https://doi.org/10.1146/annurev-arplant- 070109 -104628 PMID: 20192754
25. Cuevas B , Abell A , Johnson G . Role of mitogen-activated protein kinase kinase kinases in signal integration . Oncogene . 2007 ; 26 : 3159 ± 3171 . https://doi.org/10.1038/sj.onc. 1210409 PMID: 17496913
26. Cargnello M , Roux P. Activation and Function of the MAPKs and Their Substrates, the MAPK-Activated Protein Kinases . Microbiol Mol Biol R . 2011 ; 75 ( 1 ): 50 ± 83 .
27. Morrison D. MAP Kinase Pathways . CSH Perspect Biol . 2012 ; 4 : a011254 .
28. Teets N , Elnitsky M , Benoit J , Lopez-Martinez G , Denlinger D , Lee R . Rapid cold-hardening in larvae of the Antarctic midge Belgica antarctica: cellular cold-sensing and a role for calcium . Am J Physiol Regul Integr Comp Physiol . 2008 ; 294 : 1938 ± 1946 .
29. Conesa A , GoÈtz S , GarcÂõa-GoÂmez JM , Terol J , TaloÂn M , Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research . Bioinformatics. 2005 ; 21 ( 18 ): 3674 ± 3676 . https://doi.org/10.1093/bioinformatics/bti610 PMID: 16081474
30. Kanehisa M , Araki M , Goto S , Hattori M , Hirakawa M , Itoh M , et al. KEGG for linking genomes to life and the environment . Nucleic Acids Res . 2008 ; 36 ( suppl 1 ):D480± D484 .
31. Mortazavi A , Williams B , McCue K , Schaeffer L , Wold B . Mapping and quantifying mammalian transcriptomes by RNA-Seq . Nat Methods . 2008 ; 5 ( 7 ): 621 ± 628 . https://doi.org/10.1038/nmeth.1226 PMID: 18516045
32. Audic S , Claverie J. The significance of digital gene expression profiles . Genome Res . 1997 ; 7 ( 10 ): 986 ± 995 . PMID: 9331369
33. Livak K , Schmittgen T. Analysis of relative gene expression data using real-time quantitative PCR and the 2− ΔΔCT Method . Methods . 2001 ; 25 ( 4 ): 402 ± 408 . https://doi.org/10.1006/meth. 2001 .1262 PMID: 11846609
34. Abele D , Puntarulo S. Formation of reactive species and induction of antioxidant defence systems in polar and temperate marine invertebrates and fish . Comp Biochem Phys A . 2004 ; 138 ( 4 ): 405 ± 15 .
35. Nie H , Jiang L , Huo Z , Liu L , Yang F , Yan X . Transcriptomic responses to low temperature stress in the Manila clam, Ruditapes philippinarum . Fish Shellfish Immunol . 2016 ; 55 : 358 ± 66 . https://doi.org/10. 1016/j.fsi. 2016 . 06 .008 PMID: 27288255
36. Nadal E , Ammerer G , Posas F . Controlling gene expression in response to stress . Nat Rev Genet . 2011 ; 12 : 833 ± 845 . https://doi.org/10.1038/nrg3055 PMID: 22048664
37. Fan L , Wang A , Wu Y . Comparative proteomic identification of the hemocyte response to cold stress in white shrimp, Litopenaeus vannamei . J Proteomics . 2013 ; 80 : 196 ± 206 . https://doi.org/10.1016/j.jprot. 2012 . 12 .017 PMID: 23396037
38. Lu X , Kong J , Luan S , Dai P , Meng X , Cao B , et al. Transcriptome Analysis of the Hepatopancreas in the Pacific White Shrimp (Litopenaeus vannamei) under Acute Ammonia Stress . PloS ONE . 2016 ; 11 ( 10 ):e0164396. https://doi.org/10.1371/journal.pone. 0164396 PMID: 27760162
39. Hu D , Pan L , Zhao Q , Ren Q . Transcriptomic response to low salinity stress in gills of the Pacific white shrimp, Litopenaeus vannamei . Mar Genom . 2015 ; 24 : 297 ± 304 .
40. Quispe R , Justino E , Vieira F , Jaramillo M , Rosa R , Perazzolo L. Transcriptional profiling of immunerelated genes in Pacific white shrimp (Litopenaeus vannamei) during ontogenesis . Fish Shellfish Immunol . 2016 ; 58 : 103 ± 107 . https://doi.org/10.1016/j.fsi. 2016 . 09 .024 PMID: 27637731
41. Zhou J , Wang L , Xin Y , Wang W , He W , Wang A , et al. Effect of temperature on antioxidant enzyme gene expression and stress protein response in white shrimp, Litopenaeus vannamei . J Therm Biol . 2010 ; 35 : 284 ± 289 .
42. Vinagre C , Madeira D , MendoncËa V , Dias M , Roma J , Diniz M. Effect of temperature in multiple biomarkers of oxidative stress in coastal shrimp . J Therm Biol . 2014 ; 41 : 38 ± 42 . https://doi.org/10.1016/j. jtherbio. 2014 . 02 .005 PMID: 24679970
43. Qiu J , Wang W , Wang, Liu Y , Wang A . Oxidative stress, DNA damage and osmolality in the Pacific white shrimp, Litopenaeus vannamei exposed to acute low temperature stress . Comp Biochem Phys C . 2011 ; 154 : 36 ± 41 .
44. Li B , Xian J , Guo H , Wang A , Miao Y , Ye J , et al. Effect of temperature decrease on hemocyte apoptosis of the white shrimp Litopenaeus vannamei . Aquacult Int . 2014 ; 22 : 761 ± 774 .
45. Smith S , Bernatchez L , Beheregaray L . RNA-Seq analysis reveals extensive transcriptional plasticity to temperature stress in a freshwater fish species . BMC Genomics . 2013 ; 14 : 375 . https://doi.org/10.1186/ 1471 -2164-14-375 PMID: 23738713
46. Yin F , Sun P , Wang J , Gao Q . Transcriptome analysis of dormant tomonts of the marine fish ectoparasitic ciliate Cryptocaryon irritans under low temperature . Parasite Vector . 2016 ; 9 : 280 .
47. Teets N , Yi S , Lee R , Denlinger D. Calcium signaling mediates cold sensing in insect tissues . PNAS . 2013 ; 110 ( 22 ): 9154 ± 9159 . https://doi.org/10.1073/pnas.1306705110 PMID: 23671084