- Email: [email protected]
Comparative transcriptomic analysis of shrimp hemocytes in response to acute hepatopancreas necrosis disease (AHPND) causing Vibrio parahemolyticus infection
Accepted Manuscript Comparative transcriptomic analysis of shrimp hemocytes in response to acute hepatopancreas necrosis disease (AHPND) causing Vibrio parahemolyticus infection Zhihong Zheng, Fan Wang, Jude Juventus Aweya, Ruiwei Li, Defu Yao, Mingqi Zhong, Shengkang Li, Yueling Zhang PII:
S1050-4648(17)30773-8
DOI:
10.1016/j.fsi.2017.12.032
Reference:
YFSIM 5017
To appear in:
Fish and Shellfish Immunology
Received Date: 17 July 2017 Revised Date:
8 December 2017
Accepted Date: 20 December 2017
Please cite this article as: Zheng Z, Wang F, Aweya JJ, Li R, Yao D, Zhong M, Li S, Zhang Y, Comparative transcriptomic analysis of shrimp hemocytes in response to acute hepatopancreas necrosis disease (AHPND) causing Vibrio parahemolyticus infection, Fish and Shellfish Immunology (2018), doi: 10.1016/j.fsi.2017.12.032. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1
ACCEPTED MANUSCRIPT Comparative transcriptomic analysis of shrimp hemocytes in response to
2
acute hepatopancreas necrosis disease (AHPND) causing Vibrio
3
parahemolyticus infection
4
Zhihong Zheng#, Fan Wang#, Jude Juventus Aweya, Ruiwei Li, Defu Yao, Mingqi Zhong,
6
Shengkang Li, Yueling Zhang*
7
RI PT
5
Department of Biology and Guangdong Provincial Key laboratory of Marine Biotechnology,
9
Shantou University, Shantou, Guangdong 515063, China
M AN U
SC
8
10 11
Running title: Shrimp hemocytes response to AHPND V. parahemolyticus infection
12
#
13
*Correspondence to:
These authors contributed equally to this work.
Prof. Yue-ling Zhang, Ph.D. Department of Biology
15
School of Science
16
Shantou University
19 20 21
EP
18
Shantou, Guangdong 515063, China Tel: +86-754-82902580
AC C
17
TE D
14
Fax: +86-754-82902767 E-mail: [email protected]
22 23 24 25 1
ACCEPTED MANUSCRIPT
Abstract
27
The recent emergence of acute hepatopancreas necrosis disease (AHPND) in shrimps has
28
posed a major challenge in the shrimp aquaculture industry. The Pir toxin proteins carried
29
by some strains of Vibrio parahaemolyticus are believed to play essential roles in the
30
pathogenesis of AHPND. However, few studies have so far explored how the host immune
31
system responds to these bacteria. In this study, AHPND V. parahaemolyticus (with Pir)
32
and non-AHPND V. parahaemolyticus (without Pir) were injected into two groups of
33
shrimps, and the hemocytes collected for comparative transcriptomic analyses. A total of
34
1064 differentially expressed genes (DEGs) were identified, of which 910 were
35
up-regulated and 154 were down-regulated. Gene Ontology (GO) annotation and Kyoto
36
Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed that many DEGs
37
were involved in a number of biological processes such as cellular process, metabolic
38
process and single-organism process in the AHPND V. parahaemolyticus injected group
39
than the non-AHPND V. parahaemolyticus injected group. Among these, major metabolic
40
processes such as carbohydrate metabolism, lipid metabolism and amino acid metabolism
41
were further identified as the major responsive gene groups. We observed that genes
42
involved in cell growth and anti-apoptosis including src, iap2, cas2, cytochrome P450, gst
43
and cytochrome c oxidase were strongly activated in the AHPND V. parahaemolyticus
44
group than in the non-AHPND V. parahaemolyticus group. Collectively, our results
45
unveiled that shrimp hemocytes respond to AHPND related strain of Vibrio
46
parahaemolyticus infection at the transcriptional level, which is useful in furthering our
47
understanding of AHPND.
48
Keywords: Litopenaeus vannamei; AHPND; Pir toxin; transcriptome; metabolism;
49
anti-apoptosis
AC C
EP
TE D
M AN U
SC
RI PT
26
50 51 2
ACCEPTED MANUSCRIPT
1. Introduction
53
The Pacific white leg shrimp, Litopenaeus vannamei, is one of the most important shrimp
54
species cultured worldwide and especially in China. According to the annual report of the
55
Food & Agriculture Organization (FAO), a total of 7, 894, 204 metric tons (mt) of shrimp
56
was supplied globally in 2014 with over 4, 303, 204 mt harvested directly from aquaculture
57
facilities [1]. However, the development of the global shrimp industry is saddled with
58
various diseases caused by viruses, bacteria, fungi and protozoa, hence impede the growth
59
of shrimp aquaculture [2].
SC
RI PT
52
60
Acute hepatopancreatic necrosis disease (AHPND), an emerging disease with an initial
62
outbreak in Asia around 2009 as a destructive disease, has caused massive production
63
losses in shrimp aquaculture [3]. The causative agent of AHPND is reported to be a special
64
strain of AHPND Vibrio parahemolyticus, which carries an extrachromosomal plasmid with
65
a region coding toxin genes [4]. These toxin-encoding genes, which are called
66
photorhabdus insect-related (Pir) toxins, are currently not detected in other non-AHPND V.
67
parahemolyticus species. The Pir toxins are binary toxins, designated PirA and PirB, which
68
produce typical symptoms of AHPND, and induce cell death of shrimp hepatopancreas,
69
including massive sloughing of hepatopancreas’s tubule epithelial cells [3, 5-6]. In inducing
70
hepatopancreatic cell death, it is believed that V. parahemolyticus may first attach and
71
colonize the stomach of the shrimp with its peritrichous pili-like structures [7], and then the
72
soluble Pir toxins are released into the hepatopancreas to induce organic necrosis by an
73
unknown mechanism [8]. Following the colonization of the shrimp’s stomach, other
74
opportunistic bacteria can then invade the circulatory system thereby leading to the
75
aggregation of hemocytes in the hepatopancreas [9], hence, triggering the systemic
76
inflammatory response. However, it is still largely unknown about how shrimp immune
77
system especially hemocytes respond to AHPND V. parahaemolyticus infection, which
AC C
EP
TE D
M AN U
61
3
ACCEPTED MANUSCRIPT 78
might be helpful to explain how these toxins kill the shrimp.
79
To address this issue, shrimps were challenged with normal saline, non-AHPND V.
81
parahaemolyticus and AHPND V. parahaemolyticus respectively in this study. Total RNAs
82
extracted from three independent treated samples per group were pooled, followed by
83
RNA-Seq analysis and the responsive genes compared and analyzed using GO term
84
enrichment and KEGG pathway enrichment. Our results provide a different perspective of
85
the impact of AHPND V. parahaemolyticus on hemocytes, in the pathogenesis of AHPND.
SC
RI PT
80
86
2. Materials and methods
88
2.1. Shrimp culture and pathogen challenge
89
Litopenaeus vannamei (weighing about 10-12g) were purchased from a local shrimp farm
90
(Niutianyang, Shantou, Guangdong, China) and cultured in laboratory tanks filled with
91
aerated seawater at 26°C. Shrimps were acclimatized for 2 days before challenge
92
experiments. Two types of V. parahemolyticus were used in this experiment, i.e., AHPND
93
Vibrio parahaemolyticus (PD-2) and non-AHPND Vibrio parahaemolyticus (BCRC12959)
94
strains, which were kind gifts from Dr. Chu Fang Lo of National Cheng Kung University.
95
The AHPND V. parahaemolyticus strain (denoted VA) was identified as the pathogen
96
responsible for AHPND in shrimp, while the non-AHPND V. parahaemolyticus strain
97
(denoted VN) was used as control. Bacteria at the logarithmic growth phase were cultured
98
in an improved Luria-Bertani media (tryptone: 10g/L, yeast extract: 5g/L, NaCl 30g/L)
99
at 30 °C with shaking at 150 rpm until the OD600 of the culture reached 0.6, counted and
100
diluted to the appropriate CFU with 0.65 % normal saline. Shrimps were randomly grouped
101
into an experimental group and two control groups (30 shrimps per group). For the
AC C
EP
TE D
M AN U
87
4
ACCEPTED MANUSCRIPT experimental group, 100 ul of 1.0×105CFU/ml VA was injected via the second abdominal
103
segment of shrimp, while for the controls, 100 ul of 1.0×105CFU/ml VN or 0.65 %
104
normal saline (denoted NS) was injected as above. For mortality analysis, forty shrimps
105
from each group were used and shrimp mortality was recorded at two hours interval till 96
106
hours post injection. Three independent experiments were performed with different batches
107
of shrimps, and the Log-rank (Mantel-Cox) Test method was used to analyze the
108
differences in cumulative mortality between groups. The mortality graphs were plotted
109
using GraphPad Prism (Version 5.01).
SC
RI PT
102
M AN U
110
2.2. RNA extraction and transcriptome sequencing
112
Hemolymph was collected at 24 hours post infection from the pericardial sinus with a
113
sterile needle and syringe into an equal volume of ice-cold anti-coagulant buffer (450 mM
114
NaCl, 10 mM KCl, 10 mM EDTA-Na2, 10 mM HEPES. pH 7.0). The hemocytes were
115
harvested by centrifugation at 800 g for 10 min at 4 °C and then immediately used for total
116
RNA extraction. Total RNA was extracted from hemocytes with mirVana miRNA Isolation
117
Kit (Ambion) according to the manufacturer’s instructions. The concentration and purity of
118
total
119
(Nano-drop Technologies, Wilmington, DE). The total RNA integrity was evaluated using
120
an Agilent 2100 Bioanalyzer system (Agilent Technologies, Palo Alto, California.) with a
121
RNA integrity number (RIN) > 8.5. Only high quality RNA samples were used for
122
constructing the cDNA libraries.
was
quantified
with
a
NanoDrop
2000
spectrophotometer
AC C
RNA
EP
TE D
111
123 124
To prepare the cDNA libraries for sequencing, total RNA from each sample was diluted to
125
200 ng/µL, and then three independent RNA samples for each treatment were pooled into
126
one. Three RNA pools, i.e., VA, VN and NS, were produced and then sequenced at the 5
ACCEPTED MANUSCRIPT Beijing Genomics Institute (BGI, Wuhan, China). Briefly, mRNA with poly (A) were
128
isolated from the total RNA using oligo (dT) beads, purified, fragmented (100 bp~400 bp)
129
with an ultra-sonicator and reverse transcribed into first strand cDNA using random primers.
130
Subsequently, sequencing adapters were connected to the short fragments, and the resultant
131
cDNA libraries used for paired-end RNA-seq. Assembled sequence data from this article
132
were submitted to GenBank under accession numbers PRJNA385392.
RI PT
127
133
2.3. Assembly and functional annotation
135
Raw reads generated using the Illumina Hiseq 2000 platform were filtered to remove
136
adaptors, sequences with unknown nucleotides larger than 5 %, and low quality reads (more
137
than 20 % Q ≤ 10 bases). Following this, the clean reads were used for de novo assembly to
138
produce unigenes using Trinity software (release-20130225) [10]. Functional assignment of
139
unigenes was carried out by Blastx with the protein databases Nr, Swiss-Prot, KEGG and
140
COG (e-value < 0.00001), and by Blastn with the nucleotide database Nt (e-value <
141
0.00001). Next, Blast2GO software was used to obtain the Gene Ontology (GO)
142
annotations [11]. The number of unigenes annotated based on each database was noted.
TE D
M AN U
SC
134
EP
143
2.4. Expression pattern and clustering of differentially expressed genes (DEGs)
145
For each sample, clean reads were mapped to unigenes using Bowtie2 (v2.2.5) [12], and the
146
expression of unigenes calculated using the Fragments Per kb per Million reads (FPKM)
147
method [13]. The read counts of all comparisons were normalized to the aligned FPKM to
148
obtain the relative expression. To ensure high quality of DEGs, False Discovery Rate (FDR)
149
control method was used, with a smaller FDR indicating a larger difference in expression
150
between two samples. Unigenes with FDR ≤0.001 and |log2Ratio|≥1 were chosen as the
151
DEGs [14].
AC C
144
152 6
ACCEPTED MANUSCRIPT 2.5. GO term and KEGG pathway enrichment analysis of DEGs
154
Using the GO term enrichment and KEGG pathway analysis results, the significant GO
155
terms in DEGs were further enriched using hypergeometric test (p≤0.05) [15]. Significantly
156
enriched metabolic pathways or signal transduction pathways in DEGs with Q≤0.05 were
157
chosen for further analysis. Next, WEGO software was used to obtain GO functional
158
classifications and distribution of the genes [16], while the biological functions of unigenes
159
were analyzed using the online KEGG Automatic Annotation Server (KAAS)
160
(http://www.genome.jp/kegg/kaas/) [17].
SC
RI PT
153
161
2.6. Validation of differentially expression genes by qRT-PCR
163
To validate the transcriptome data, 1 µg of the high quality total RNA samples, extracted
164
using the mirVana
165
reverse-transcribed using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara,
166
Dalian, China), following the manufacturer’s protocol. Six unigenes were selected for
167
validation by qRT-PCR using a LightCycler 480 RT-PCR system (Roche Applied Science,
168
Switzerland) with the following program: one cycle at 95 °C for 10 min and 45 cycles of
169
95 °C for 15 s and 60 °C for 30 s. The qPCR was carried out in triplicate for each sample,
170
and five shrimps were analyzed in each group. The fold changes in gene expression were
171
computed using the relative quantification method and normalized to L. vannamei EF1α
172
expression. The gene specific qRT-PCR primers designed with Primer 5 software are listed
173
in Table S1.
(Ambion, Austin,
TX, USA) was
EP
TE D
miRNA Isolation Kit
AC C
174
M AN U
162
175
2.7. Statistical analyses
176
All data are presented as means ± standard deviation (SD). Data Normality test was
177
checked by the Shapiro-Wilk test. One-way analysis of variance (ANOVA) and a multiple
178
comparison (Tukey) test were used to compare the significant differences with the SPSS 7
ACCEPTED MANUSCRIPT 179
13.0 program (SPSS Inc., Chicago, IL, USA). P< 0.05 was considered statistically
180
significant.
181
3. Result
183
3.1. VA treated group showed higher mortality than VN treated group
184
To examine which of these bacteria (VA and VN) causes more mortality in shrimp, the
185
same dose (104cfu) of VA and VN was injected into shrimps, with about 40 shrimps per
186
each treatment group. The results showed a cumulative mortality of about 97.5% in VA at
187
96 hours 4 days post infection (hpi), while only 62.5% cumulative mortality was observed
188
in VN at the same time point (Fig 1).
M AN U
189
SC
RI PT
182
3.2. De novo assembly and annotation of Unigenes
191
Three cDNA libraries representing VA, VN and NS groups were constructed. A total of
192
145,161,854 raw reads were generated from these three groups (Table 1). After data
193
cleaning, 45,658 assembled sequences with an average length of 1324 nt were obtained
194
using Trinity (Table 2). The sequence length (nt) ranges from 300 nt to ≥3000 nt with
195
distribution shown in Fig 2. Next, all unigenes were annotated using Blastx to protein
196
databases (Nr, Nt, Swiss-Prot, KEGG, COG and GO) and Blastn to nucleotide databases
197
(Nt), with 26634 unigenes assigned, which represented about 58.33% of the total (Table 3).
EP
AC C
198
TE D
190
199
In the GO term annotation, a total of 11,457 unigenes from the assembly were assigned GO
200
terms of biological process (23 subcategories), cellular component (16 subcategories) and
201
molecular function (17 subcategories). The major GO term subcategories (i.e., representing
202
top 30 % of unigenes) were cellular process, cell, cell part, metabolic process, binding,
203
single-organism process, catalytic activity，organelle, and biological regulation (Fig 3A).
204
Using the COG database to further understand the protein orthologies of the assembled 8
ACCEPTED MANUSCRIPT unigenes, 11,403 were successfully annotated into 26 specific protein function definitions
206
or orthologous categories (Fig 3B). Among these protein function categories, the largest
207
four were general function prediction only (5753); translation, ribosomal structure and
208
biogenesis (3605); transcription (2684); and replication, recombination and repair (2340).
209
To identify the biological processes by which the annotated unigenes are involved in,
210
KEGG pathway analysis was further carried out. In the KEGG analysis, 18,599 unigenes
211
were annotated into 6 major KEGG pathways, including cellular processes, environmental
212
information processing, genetic information processing, human diseases, metabolism and
213
organismal systems. These annotated unigenes were further divided into 42 level 2
214
subcategories pathways. Except the global map which had no pictorial information, the
215
largest subcategory group, infection diseases (Bacterial), had 2609 annotated genes,
216
followed by signal transduction (2400), immune system (2089) and translation (2060) (Fig
217
3C).
218
Table S2.
SC
M AN U
Apart from these, 299 level 3 KEGG subcategories were annotated and are listed in
TE D
219
RI PT
205
3.3. Classification and analysis of the differentially expressed genes (DEGs)
221
To analyze and characterize the DEGs following the different treatments, a cutoff of
222
FDR (false discovery rate) ≤ 0.05 and log2 Ratio ≥ 1 was used as the threshold for selecting
223
DEGs. Based on this, a large number of genes were observed to be dysregulated in VA
224
group (1988) than in VN group (743) relative to control group (NS), while 1064 genes were
225
differentially expressed between VA group and VN group (Fig. 4A and Table S3-S6). To
226
further analyze the DEGs between the AHPND V. parahaemolyticus and non-AHPND V.
227
parahaemolyticus samples, DEGs were further annotated with GO and KEGG databases. In
228
the GO enrichment analysis, a total of 345 DEGs were successfully enriched in biological
229
process (231), cellular component (154) and molecular function (281) (Fig 4B). For the
AC C
EP
220
9
ACCEPTED MANUSCRIPT KEGG pathway enrichment analysis, 709 of the DEGs were annotated, relative to all
231
unigenes and between VA and VN (Fig 4C, Table S7). When VA was compared to VN, 20
232
subcategories of the KEGG pathway were enriched with the highly significant biological
233
pathways being carbohydrate metabolism, xenobiotics biodegradation and metabolism,
234
infectious diseases (bacterial), lipid metabolism and amino acid metabolism (Fig 4C (i)).
235
On the other hand, 22 subcategories were enriched with lower but significantly changes in
236
signal transduction, infectious diseases (viral), cancers (specific types), immune system and
237
nervous system (Fig 4C (ii)).
SC
RI PT
230
M AN U
238
3.4. Validation of RNA-seq transcriptome data by qRT-PCR
240
In order to substantiate the transcriptome data as well as better understand the response of
241
shrimp hemocytes to AHPND V. parahaemolyticus and non-AHPND V. parahaemolyticus,
242
validation of the transcriptome data using qRT-PCR was performed. As mentioned above,
243
the AHPND V. parahemolyticus can induced massive sloughing of hepatopancreas tubule
244
epithelial cells [6], therefore to further explore the detailed hemocytes response induced by
245
AHPND V. parahemolyticus, DEGs were categorized with KEGG pathway enrichment
246
analyses in Fig. 4C. We found that the genes involved in cell growth and anti-apoptosis
247
were strongly induced in VA group compared with VN group in the transcriptome (Table
248
S8). Thus, six genes, src, iap2, cas2, cytochrome P450, gst and cytochrome c oxidase, were
249
selected for validation using qPCR (Fig 5). The qPCR results showed a significant
250
up-regulation in the expression of the oncogene src as well as the anti-apoptosis genes iap2
251
and cas2 in VA relative to VN and NS. Similarly, a significant up-regulation was observed
252
in the expression of the oxidoreductase gst, cytochrome P450 and cytochrome c oxidase in
253
VA relative to VN or NS. Thus, the qPCR results were generally consistent with the
254 255
transcriptome data.
AC C
EP
TE D
239
10
ACCEPTED MANUSCRIPT
4. Discussion
257
AHPND has caused severe losses in the shrimp industry in the last couple of years [5]. To
258
curtail this loss, several different control or preventive methods including the use of
259
inactivated bacteria antigen, IgY feed additives, phage therapy and nano-bubble technology
260
[18] have been carried out, but with little or no success. This enigma drove scientist to
261
comprehensively explore this disease. In this study, we have mimicked the late stage of
262
AHPND with intramuscularly injection of two strains of V. parahemolyticus, which allow a
263
direct contact between the bacteria and hemocytes. Following this, a comparative
264
transcriptomic analysis was carried out. Our results revealed that shrimp hemocytes
265
probably had a strong anti-apoptosis response against VA infection.
M AN U
SC
RI PT
256
266
As the most important immunocytes in shrimps, hemocytes are the major guardians against
268
pathogen invasion in the circulatory system [19]. The morphology and function of
269
hemocytes is reported to be affected by various bacteria components, as LPS injection has
270
been shown to cause significant decrease in total hemocyte count (THC) in Penaeus
271
monodon [20]. Similarly, the exotoxin from Clostridium difficile has been reported to
272
induce cytoskeletal reorganization in lepidopteran pests [21]. In the present study, our
273
results have shown that VA could induce stronger responses than VN, as indicated by the
274
higher cumulative mortality (Fig 1) as well as the higher number of differentially expressed
275
genes (DEGs) between different V. parahemolyticus treated hemocytes (Fig 4). To further
276
explore the detail response of the hemocytes, DEGs were categorized using KEGG pathway
277
analysis and revealed that anti-apoptosis-related pathways were more strongly activated in
278
VA. Apoptosis is a genetically programmed cell death mechanism, which could be induced
AC C
EP
TE D
267
11
ACCEPTED MANUSCRIPT by a wide range of stimuli, via the caspase cascade [22]. Our results suggest that VA could
280
not only induce hepatopancreas necrosis, but also cause hemocytes apoptosis, which greatly
281
impaired shrimp immune system. Thus, some anti-apoptosis genes like the oncogene src
282
kinase, which plays a key role in cell growth, division, migration, and survival signaling
283
pathways [23], were highly expressed in VA than in VN. Our results reveals how
284
hemocytes resist the stress exerted by AHPND V. parahemolyticus challenge in terms of
285
transcriptional response.
SC
RI PT
279
286
In conclusion, this comparative transcriptomic study depicts how shrimp hemocytes
288
respond to AHPND V. parahemolyticus infection. Future studies will explore how the
289
shrimp immune system can be prompted to respond to AHPND strain of bacteria and to
290
eventually find treatment for AHPND.
TE D
291
M AN U
287
Disclosure statement
293
The authors declare no conflicts of interest.
294
EP
292
Acknowledgments
296
This work was sponsored by National Natural Science Foundation of China (Nos.
297
31672689 & 31372558), Natural Science Foundation of Guangdong Province (No.
298
2017A030311032), Shantou University Scientific Research Foundation for Talents (No.
299
NTF16001) and Guangdong SAIL Foundation for Distinguished Scholars (No. 14600703).
AC C
295
12
ACCEPTED MANUSCRIPT 300 301 302
References [1] S. Interpreters, FAO (Food and Agriculture Organization of the United Nations), International Organization 1 (2014) 121-123. [2] S. Thitamadee, A. Prachumwat, J. Srisala, et al., P. Jaroenlak , P. V. Salachan, K.
304
Sritunyalucksana, et al., Review of current disease threats for cultivated penaeid shrimp
305
in Asia, Aquaculture, 452 (2016) 69-87.
RI PT
303
[3] L. Tran, L. Nunan, R.M. Redman, L.L. Mohney, C.R. Pantoja, K. Fitzsimmons, et al.,
307
Determination of the infectious nature of the agent of acute hepatopancreatic necrosis
308
syndrome affecting penaeid shrimp, Dis. Aquat. Organ. 105(1) (2013) 45-55.
SC
306
[4] Y.T. Yang, I.T. Chen, C.T. Lee, C.Y. Chen, S.S. Lin, L.I. Hor, et al., Draft genome
310
sequences of four strains of Vibrio parahaemolyticus, three of which cause early
311
mortality syndrome/acute hepatopancreatic necrosis disease in shrimp in China and
312
Thailand, Genome Announc. 2 (5) (2014) e00816-14.
314
[5] D.V. Lightner, R.M. Redman, C.R. Pantoja, B.L. Noble, L. Tran, Early mortality
TE D
313
M AN U
309
syndrome affects shrimp in Asia, Global Aquaculture Advocate. 15 (1) (2012) 40. [6] C.T. Lee, I.T. Chen, Y.T. Yang, T.P. Ko, Y.T. Huang, J.Y. Huang, et al., The opportunistic
316
marine pathogen Vibrio parahaemolyticus becomes virulent by acquiring a plasmid that
317
expresses a deadly toxin, Proc Natl Acad Sci U. S. A. 112 (34) (2015) 10798-10803.
EP
315
[7] W. Soonthornchai, S. Chaiyapechara, P. Jarayabhand, K. Soderhall, P. Jiravanichpaisal,
319
Interaction of Vibrio spp. with the inner surface of the digestive tract of Penaeus
320
monodon, PLoS ONE. 10 (8) (2015) e0135783.
AC C
318
321
[8] H.C. Lai, T.H. Ng, M. Ando, C.T. Lee, I.T. Chen, J.C. Chuang, et al., Pathogenesis of
322
acute hepatopancreatic necrosis disease (AHPND) in shrimp, Fish Shellfish Immunol. 47
323
(2) (2015) 1006-1014.
324
[9] W. Soonthornchai, S. Chaiyapechara, S. Klinbunga, W. Thongda, S. Tangphatsornruang,
325
T. Yoocha, et al., Differentially expressed transcripts in stomach of Penaeus monodon in 13
ACCEPTED MANUSCRIPT 326
response to AHPND infection, Dev. Comp. Immunol. 65 (2016) 53-63.
327
[10] M.G. Grabherr, B.J. Haas, M. Yassour, J.Z. Levin, D.A. Thompson, I. Amit, et al.,
328
Full-length transcriptome assembly from RNA-Seq data without a reference genome, Nat.
329
Biotechnol. 29 (7) (2011) 644-652. [11] A. Conesa, S. Gotz, J.M. Garcia-Gomez, J. Terol, M. Talon, M. Robles, Blast2GO: a
331
universal tool for annotation, visualization and analysis in functional genomics research,
332
Bioinformatics. 18 (2005) 3674-3676.
334
[12] B. Langmead, S.L. Salzberg, Fast gapped-read alignment with Bowtie 2, Nat. Methods.
SC
333
RI PT
330
9 (4) (2012) 357-359.
[13] A. Mortazavi, B.A. Williams, K. McCue, L. Schaeffer, B. Wold, Mapping and
336
quantifying mammalian transcriptomes by RNA-Seq, Nat. Methods 5 (7) (2008)
337
621-628.
340 341 342 343
false discovery rate controlling procedures, Bioinformatics. 19 (3) (2003) 368-375.
TE D
339
[14] A. Reiner, D. Yekutieli, Y. Benjamini, Identifying differentially expressed genes using
[15] M.D. Young, M.J. Wakefield, G.K. Smyth, A. Oshlack, Gene ontology analysis for RNA-seq: accounting for selection bias, Genome Biol. 11 (2) (2010) R14. [16] J. Ye, L. Fang, H. Zheng, Y. Zhang, J. Chen, Z. Zhang, et al., WEGO: a web tool for
EP
338
M AN U
335
plotting GO annotations, Nucleic Acids Res. 34 (2) (2006) 293-297. [17] Y. Moriya, M. Itoh, S. Okuda, A.C. Yoshizawa, M. Kanehisa, KAAS: an automatic
345
genome annotation and pathway reconstruction server, Nucleic Acids Res. 35 (2) (2007)
346
182-185.
AC C
344
347
[18] I. Hirono, S. Tinwongger, Y. Nochiri, H. Kondo, Latest research on acute
348
hepatopancreatic necrosis disease (AHPND) of penaeid shrimps, Philippines.(2016)
349
3-10.
350
[19] L. Vazquez, J. Alpuche, G. Maldonado, C. Agundis, A. Pereyra-Morales, E. Zenteno,
351
Review: Immunity mechanisms in crustaceans, Innate Immunol. 15 (3) (2009) 179-188. 14
ACCEPTED MANUSCRIPT 352
[20] J.A. Xian, X.X. Zhang, H. Guo, D.M. Wang, A.L. Wang, Cellular responses of the tiger
353
shrimp Penaeus monodon haemocytes after lipopolysaccharide injection, Fish Shellfish
354
Immunol 54 (2016) 385-390. [21] A. Castagnola, S.P. Stock, Common virulence factors and tissue targets of
356
entomopathogenic bacteria for biological control of Lepidopteran Pests, Insects. 5 (1)
357
(2014) 139-166.
360 361
cell survival: a review, Biochim. Biophys. Acta. 1813 (1) (2011) 238-259.
[23] R.J. Roskoski, Src protein-tyrosine kinase structure, mechanism, and small molecule inhibitors, Pharmacol Res. 94 (2015) 9-25.
362 363 364
EP AC C
367
TE D
365 366
SC
359
[22] L. Portt, G. Norman, C. Clapp, M. Greenwood, M.T. Greenwood, Anti-apoptosis and
M AN U
358
RI PT
355
15
ACCEPTED MANUSCRIPT
Figure legends
369
Fig. 1. Cumulative mortality graphs for VA, VN, and NS injected shrimps. Forty shrimps
370
were injected with VA (AHPND V. parahemolyticus) or VN (non-AHPND V.
371
parahemolyticus), at the same dosage of 104cfu per shrimp or with NS (Sterile normal
372
saline). The mortality tests were carried out in triplicates for each sample. Significant
373
difference is indicated by asterisks, **p < 0.01.
374
Fig. 2. Length distribution of all assembled unigenes.
SC
375
RI PT
368
376
Fig. 3. Functional enrichment of assembled genes. (A) Gene ontology (GO) classifications
378
of non-redundant unigenes. All annotated unigenes were categorized into 3 categories: A:
379
biological process, B: cellular component, C: molecular function. (B) The cluster of
380
orthologous groups (COG) classification of putative proteins. In total, 11,403 unigenes
381
were successfully annotated to 26 categories. (C) The Kyoto Encyclopedia of Genes and
382
Genomes(KEGG) classification of non-redundant unigenes. In total, 18,599 unigenes were
383
assigned to 6 special KEGG pathways, including Cellular Processes (a), Environmental
384
Information Processing (b), Genetic Information Processing (c), Human Diseases (d),
385
Metabolism (e), Organismal Systems (f).
TE D
EP
AC C
386
M AN U
377
387
Fig. 4. Analysis of differentially expressed genes between VN (non-AHPND V.
388
parahemolyticus) and VA (AHPND V. parahemolyticus) treated groups. (A) Comparison of
389
the differentially expressed genes (DEGs) between groups. (B) GO terms functional
390
enrichment analysis for VA and VN. (C) KEGG pathway enrichment analysis and
391
annotation compared between all unigenes and DEGs between VA and VN, (a) KEGG
392
pathway enrichment analysis between VA and VN with percentages higher than all
393
unigenes; (b). KEGG pathway enrichment analysis between VA and VN with percentages 16
ACCEPTED MANUSCRIPT 394
lower than all unigenes.
395
Fig. 5. Validation of differentially expressed genes (DEGs) by qPCR. To validate the
397
RNA-seq data, the relative gene expression levels of 6 selected DEGs were examined by
398
qPCR. EF1a was used as an internal control. Different letters on the graph denotes
399
significant statistical difference (p<0.05).
AC C
EP
TE D
M AN U
SC
RI PT
396
17
ACCEPTED MANUSCRIPT 400 401
Table 1. Summary of RNA-seq data for L. vannamei injected with AHPND and
402
non-AHPND V. parahemolyticus Samples
Total Raw Reads
Total Clean Reads
Q20 Percentage
GC Percentage
47,123,094
46,275,446
97.27%
47.52%
VN
47,941,818
47,076,076
97.71%
47.59%
NS
50,096,942
49,084,770
97.77%
47.31%
Total
145,161,854
142,436,292
97.58%
47.47%
RI PT
VA
Abbreviations: NS，normal saline treated group；VA，AHPND V. parahemolyticus treated group；VN，
404
non-AHPND V. parahemolyticus treated group.
AC C
EP
TE D
M AN U
SC
403
18
ACCEPTED MANUSCRIPT 405
Table 2. Summary of De novo assemble of L. vannamei hemocyte transcriptome Sample
Total Number
Total Length (nt)
Mean Length (nt)
N50
VA
53,879
48,513,536
900
2084
VN
53,401
49,225,356
922
2143
NS
53,910
48,998,100
All
45,658
60,455,885
RI PT
406
2131
1324
2574
AC C
EP
TE D
M AN U
SC
407
909
19
ACCEPTED MANUSCRIPT 408 409
Table 3. Functional annotation of unigenes from L. vannamei hemocytes transcriptome Annotation mode
NR
NT
Swiss-Prot
KEGG
COG
GO
ALL
23,223
14,516
21,041
18,599
11,403
11,457
26,634
50.86%
31.79%
46.08%
40.74%
Number of Unigenes with
Percentage of all assembled unigenes (45,658)
24.97%
AC C
EP
TE D
M AN U
SC
410
RI PT
hits or assignments
20
25.09%
58.33%
ACCEPTED MANUSCRIPT
Fig. 1
AC C
EP
TE D
M AN U
SC
RI PT
411
412 21
ACCEPTED MANUSCRIPT
Fig. 2
AC C
EP
TE D
M AN U
SC
RI PT
413
414 415 22
Fig. 3
417
A.
AC C
EP
TE D
M AN U
SC
416
RI PT
ACCEPTED MANUSCRIPT
418 23
ACCEPTED MANUSCRIPT
B.
AC C
EP
TE D
M AN U
SC
RI PT
419
420
24
ACCEPTED MANUSCRIPT
C.
AC C
EP
TE D
M AN U
SC
RI PT
421
422
25
Fig. 4
424
A.
AC C
EP
TE D
M AN U
SC
423
RI PT
ACCEPTED MANUSCRIPT
425 26
ACCEPTED MANUSCRIPT
B.
427
AC C
EP
TE D
M AN U
SC
RI PT
426
428
27
ACCEPTED MANUSCRIPT
C(i).
430
.
AC C
EP
TE D
M AN U
SC
RI PT
429
28
ACCEPTED MANUSCRIPT
C (ii).
AC C
EP
TE D
M AN U
SC
RI PT
431
432 433 29
ACCEPTED MANUSCRIPT
Fig. 5
AC C
EP
TE D
M AN U
SC
RI PT
434
435 30
ACCEPTED MANUSCRIPT
Research highlights ►
We found 1064 genes differentially expressed between the AHPND and Non-AHPND
V. parahemolyticus infected hemocytes in Litopenaeus vannamei. this, the genes involved in the major metabolism processes, cell growth and
RI PT
► Among
anti-apoptosis were further identified as the major responsive genes in response to the
AC C
EP
TE D
M AN U
SC
AHPND V. parahemolyticus.
S1050-4648(17)30773-8
DOI:
10.1016/j.fsi.2017.12.032
Reference:
YFSIM 5017
To appear in:
Fish and Shellfish Immunology
Received Date: 17 July 2017 Revised Date:
8 December 2017
Accepted Date: 20 December 2017
Please cite this article as: Zheng Z, Wang F, Aweya JJ, Li R, Yao D, Zhong M, Li S, Zhang Y, Comparative transcriptomic analysis of shrimp hemocytes in response to acute hepatopancreas necrosis disease (AHPND) causing Vibrio parahemolyticus infection, Fish and Shellfish Immunology (2018), doi: 10.1016/j.fsi.2017.12.032. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1
ACCEPTED MANUSCRIPT Comparative transcriptomic analysis of shrimp hemocytes in response to
2
acute hepatopancreas necrosis disease (AHPND) causing Vibrio
3
parahemolyticus infection
4
Zhihong Zheng#, Fan Wang#, Jude Juventus Aweya, Ruiwei Li, Defu Yao, Mingqi Zhong,
6
Shengkang Li, Yueling Zhang*
7
RI PT
5
Department of Biology and Guangdong Provincial Key laboratory of Marine Biotechnology,
9
Shantou University, Shantou, Guangdong 515063, China
M AN U
SC
8
10 11
Running title: Shrimp hemocytes response to AHPND V. parahemolyticus infection
12
#
13
*Correspondence to:
These authors contributed equally to this work.
Prof. Yue-ling Zhang, Ph.D. Department of Biology
15
School of Science
16
Shantou University
19 20 21
EP
18
Shantou, Guangdong 515063, China Tel: +86-754-82902580
AC C
17
TE D
14
Fax: +86-754-82902767 E-mail: [email protected]
22 23 24 25 1
ACCEPTED MANUSCRIPT
Abstract
27
The recent emergence of acute hepatopancreas necrosis disease (AHPND) in shrimps has
28
posed a major challenge in the shrimp aquaculture industry. The Pir toxin proteins carried
29
by some strains of Vibrio parahaemolyticus are believed to play essential roles in the
30
pathogenesis of AHPND. However, few studies have so far explored how the host immune
31
system responds to these bacteria. In this study, AHPND V. parahaemolyticus (with Pir)
32
and non-AHPND V. parahaemolyticus (without Pir) were injected into two groups of
33
shrimps, and the hemocytes collected for comparative transcriptomic analyses. A total of
34
1064 differentially expressed genes (DEGs) were identified, of which 910 were
35
up-regulated and 154 were down-regulated. Gene Ontology (GO) annotation and Kyoto
36
Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed that many DEGs
37
were involved in a number of biological processes such as cellular process, metabolic
38
process and single-organism process in the AHPND V. parahaemolyticus injected group
39
than the non-AHPND V. parahaemolyticus injected group. Among these, major metabolic
40
processes such as carbohydrate metabolism, lipid metabolism and amino acid metabolism
41
were further identified as the major responsive gene groups. We observed that genes
42
involved in cell growth and anti-apoptosis including src, iap2, cas2, cytochrome P450, gst
43
and cytochrome c oxidase were strongly activated in the AHPND V. parahaemolyticus
44
group than in the non-AHPND V. parahaemolyticus group. Collectively, our results
45
unveiled that shrimp hemocytes respond to AHPND related strain of Vibrio
46
parahaemolyticus infection at the transcriptional level, which is useful in furthering our
47
understanding of AHPND.
48
Keywords: Litopenaeus vannamei; AHPND; Pir toxin; transcriptome; metabolism;
49
anti-apoptosis
AC C
EP
TE D
M AN U
SC
RI PT
26
50 51 2
ACCEPTED MANUSCRIPT
1. Introduction
53
The Pacific white leg shrimp, Litopenaeus vannamei, is one of the most important shrimp
54
species cultured worldwide and especially in China. According to the annual report of the
55
Food & Agriculture Organization (FAO), a total of 7, 894, 204 metric tons (mt) of shrimp
56
was supplied globally in 2014 with over 4, 303, 204 mt harvested directly from aquaculture
57
facilities [1]. However, the development of the global shrimp industry is saddled with
58
various diseases caused by viruses, bacteria, fungi and protozoa, hence impede the growth
59
of shrimp aquaculture [2].
SC
RI PT
52
60
Acute hepatopancreatic necrosis disease (AHPND), an emerging disease with an initial
62
outbreak in Asia around 2009 as a destructive disease, has caused massive production
63
losses in shrimp aquaculture [3]. The causative agent of AHPND is reported to be a special
64
strain of AHPND Vibrio parahemolyticus, which carries an extrachromosomal plasmid with
65
a region coding toxin genes [4]. These toxin-encoding genes, which are called
66
photorhabdus insect-related (Pir) toxins, are currently not detected in other non-AHPND V.
67
parahemolyticus species. The Pir toxins are binary toxins, designated PirA and PirB, which
68
produce typical symptoms of AHPND, and induce cell death of shrimp hepatopancreas,
69
including massive sloughing of hepatopancreas’s tubule epithelial cells [3, 5-6]. In inducing
70
hepatopancreatic cell death, it is believed that V. parahemolyticus may first attach and
71
colonize the stomach of the shrimp with its peritrichous pili-like structures [7], and then the
72
soluble Pir toxins are released into the hepatopancreas to induce organic necrosis by an
73
unknown mechanism [8]. Following the colonization of the shrimp’s stomach, other
74
opportunistic bacteria can then invade the circulatory system thereby leading to the
75
aggregation of hemocytes in the hepatopancreas [9], hence, triggering the systemic
76
inflammatory response. However, it is still largely unknown about how shrimp immune
77
system especially hemocytes respond to AHPND V. parahaemolyticus infection, which
AC C
EP
TE D
M AN U
61
3
ACCEPTED MANUSCRIPT 78
might be helpful to explain how these toxins kill the shrimp.
79
To address this issue, shrimps were challenged with normal saline, non-AHPND V.
81
parahaemolyticus and AHPND V. parahaemolyticus respectively in this study. Total RNAs
82
extracted from three independent treated samples per group were pooled, followed by
83
RNA-Seq analysis and the responsive genes compared and analyzed using GO term
84
enrichment and KEGG pathway enrichment. Our results provide a different perspective of
85
the impact of AHPND V. parahaemolyticus on hemocytes, in the pathogenesis of AHPND.
SC
RI PT
80
86
2. Materials and methods
88
2.1. Shrimp culture and pathogen challenge
89
Litopenaeus vannamei (weighing about 10-12g) were purchased from a local shrimp farm
90
(Niutianyang, Shantou, Guangdong, China) and cultured in laboratory tanks filled with
91
aerated seawater at 26°C. Shrimps were acclimatized for 2 days before challenge
92
experiments. Two types of V. parahemolyticus were used in this experiment, i.e., AHPND
93
Vibrio parahaemolyticus (PD-2) and non-AHPND Vibrio parahaemolyticus (BCRC12959)
94
strains, which were kind gifts from Dr. Chu Fang Lo of National Cheng Kung University.
95
The AHPND V. parahaemolyticus strain (denoted VA) was identified as the pathogen
96
responsible for AHPND in shrimp, while the non-AHPND V. parahaemolyticus strain
97
(denoted VN) was used as control. Bacteria at the logarithmic growth phase were cultured
98
in an improved Luria-Bertani media (tryptone: 10g/L, yeast extract: 5g/L, NaCl 30g/L)
99
at 30 °C with shaking at 150 rpm until the OD600 of the culture reached 0.6, counted and
100
diluted to the appropriate CFU with 0.65 % normal saline. Shrimps were randomly grouped
101
into an experimental group and two control groups (30 shrimps per group). For the
AC C
EP
TE D
M AN U
87
4
ACCEPTED MANUSCRIPT experimental group, 100 ul of 1.0×105CFU/ml VA was injected via the second abdominal
103
segment of shrimp, while for the controls, 100 ul of 1.0×105CFU/ml VN or 0.65 %
104
normal saline (denoted NS) was injected as above. For mortality analysis, forty shrimps
105
from each group were used and shrimp mortality was recorded at two hours interval till 96
106
hours post injection. Three independent experiments were performed with different batches
107
of shrimps, and the Log-rank (Mantel-Cox) Test method was used to analyze the
108
differences in cumulative mortality between groups. The mortality graphs were plotted
109
using GraphPad Prism (Version 5.01).
SC
RI PT
102
M AN U
110
2.2. RNA extraction and transcriptome sequencing
112
Hemolymph was collected at 24 hours post infection from the pericardial sinus with a
113
sterile needle and syringe into an equal volume of ice-cold anti-coagulant buffer (450 mM
114
NaCl, 10 mM KCl, 10 mM EDTA-Na2, 10 mM HEPES. pH 7.0). The hemocytes were
115
harvested by centrifugation at 800 g for 10 min at 4 °C and then immediately used for total
116
RNA extraction. Total RNA was extracted from hemocytes with mirVana miRNA Isolation
117
Kit (Ambion) according to the manufacturer’s instructions. The concentration and purity of
118
total
119
(Nano-drop Technologies, Wilmington, DE). The total RNA integrity was evaluated using
120
an Agilent 2100 Bioanalyzer system (Agilent Technologies, Palo Alto, California.) with a
121
RNA integrity number (RIN) > 8.5. Only high quality RNA samples were used for
122
constructing the cDNA libraries.
was
quantified
with
a
NanoDrop
2000
spectrophotometer
AC C
RNA
EP
TE D
111
123 124
To prepare the cDNA libraries for sequencing, total RNA from each sample was diluted to
125
200 ng/µL, and then three independent RNA samples for each treatment were pooled into
126
one. Three RNA pools, i.e., VA, VN and NS, were produced and then sequenced at the 5
ACCEPTED MANUSCRIPT Beijing Genomics Institute (BGI, Wuhan, China). Briefly, mRNA with poly (A) were
128
isolated from the total RNA using oligo (dT) beads, purified, fragmented (100 bp~400 bp)
129
with an ultra-sonicator and reverse transcribed into first strand cDNA using random primers.
130
Subsequently, sequencing adapters were connected to the short fragments, and the resultant
131
cDNA libraries used for paired-end RNA-seq. Assembled sequence data from this article
132
were submitted to GenBank under accession numbers PRJNA385392.
RI PT
127
133
2.3. Assembly and functional annotation
135
Raw reads generated using the Illumina Hiseq 2000 platform were filtered to remove
136
adaptors, sequences with unknown nucleotides larger than 5 %, and low quality reads (more
137
than 20 % Q ≤ 10 bases). Following this, the clean reads were used for de novo assembly to
138
produce unigenes using Trinity software (release-20130225) [10]. Functional assignment of
139
unigenes was carried out by Blastx with the protein databases Nr, Swiss-Prot, KEGG and
140
COG (e-value < 0.00001), and by Blastn with the nucleotide database Nt (e-value <
141
0.00001). Next, Blast2GO software was used to obtain the Gene Ontology (GO)
142
annotations [11]. The number of unigenes annotated based on each database was noted.
TE D
M AN U
SC
134
EP
143
2.4. Expression pattern and clustering of differentially expressed genes (DEGs)
145
For each sample, clean reads were mapped to unigenes using Bowtie2 (v2.2.5) [12], and the
146
expression of unigenes calculated using the Fragments Per kb per Million reads (FPKM)
147
method [13]. The read counts of all comparisons were normalized to the aligned FPKM to
148
obtain the relative expression. To ensure high quality of DEGs, False Discovery Rate (FDR)
149
control method was used, with a smaller FDR indicating a larger difference in expression
150
between two samples. Unigenes with FDR ≤0.001 and |log2Ratio|≥1 were chosen as the
151
DEGs [14].
AC C
144
152 6
ACCEPTED MANUSCRIPT 2.5. GO term and KEGG pathway enrichment analysis of DEGs
154
Using the GO term enrichment and KEGG pathway analysis results, the significant GO
155
terms in DEGs were further enriched using hypergeometric test (p≤0.05) [15]. Significantly
156
enriched metabolic pathways or signal transduction pathways in DEGs with Q≤0.05 were
157
chosen for further analysis. Next, WEGO software was used to obtain GO functional
158
classifications and distribution of the genes [16], while the biological functions of unigenes
159
were analyzed using the online KEGG Automatic Annotation Server (KAAS)
160
(http://www.genome.jp/kegg/kaas/) [17].
SC
RI PT
153
161
2.6. Validation of differentially expression genes by qRT-PCR
163
To validate the transcriptome data, 1 µg of the high quality total RNA samples, extracted
164
using the mirVana
165
reverse-transcribed using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara,
166
Dalian, China), following the manufacturer’s protocol. Six unigenes were selected for
167
validation by qRT-PCR using a LightCycler 480 RT-PCR system (Roche Applied Science,
168
Switzerland) with the following program: one cycle at 95 °C for 10 min and 45 cycles of
169
95 °C for 15 s and 60 °C for 30 s. The qPCR was carried out in triplicate for each sample,
170
and five shrimps were analyzed in each group. The fold changes in gene expression were
171
computed using the relative quantification method and normalized to L. vannamei EF1α
172
expression. The gene specific qRT-PCR primers designed with Primer 5 software are listed
173
in Table S1.
(Ambion, Austin,
TX, USA) was
EP
TE D
miRNA Isolation Kit
AC C
174
M AN U
162
175
2.7. Statistical analyses
176
All data are presented as means ± standard deviation (SD). Data Normality test was
177
checked by the Shapiro-Wilk test. One-way analysis of variance (ANOVA) and a multiple
178
comparison (Tukey) test were used to compare the significant differences with the SPSS 7
ACCEPTED MANUSCRIPT 179
13.0 program (SPSS Inc., Chicago, IL, USA). P< 0.05 was considered statistically
180
significant.
181
3. Result
183
3.1. VA treated group showed higher mortality than VN treated group
184
To examine which of these bacteria (VA and VN) causes more mortality in shrimp, the
185
same dose (104cfu) of VA and VN was injected into shrimps, with about 40 shrimps per
186
each treatment group. The results showed a cumulative mortality of about 97.5% in VA at
187
96 hours 4 days post infection (hpi), while only 62.5% cumulative mortality was observed
188
in VN at the same time point (Fig 1).
M AN U
189
SC
RI PT
182
3.2. De novo assembly and annotation of Unigenes
191
Three cDNA libraries representing VA, VN and NS groups were constructed. A total of
192
145,161,854 raw reads were generated from these three groups (Table 1). After data
193
cleaning, 45,658 assembled sequences with an average length of 1324 nt were obtained
194
using Trinity (Table 2). The sequence length (nt) ranges from 300 nt to ≥3000 nt with
195
distribution shown in Fig 2. Next, all unigenes were annotated using Blastx to protein
196
databases (Nr, Nt, Swiss-Prot, KEGG, COG and GO) and Blastn to nucleotide databases
197
(Nt), with 26634 unigenes assigned, which represented about 58.33% of the total (Table 3).
EP
AC C
198
TE D
190
199
In the GO term annotation, a total of 11,457 unigenes from the assembly were assigned GO
200
terms of biological process (23 subcategories), cellular component (16 subcategories) and
201
molecular function (17 subcategories). The major GO term subcategories (i.e., representing
202
top 30 % of unigenes) were cellular process, cell, cell part, metabolic process, binding,
203
single-organism process, catalytic activity，organelle, and biological regulation (Fig 3A).
204
Using the COG database to further understand the protein orthologies of the assembled 8
ACCEPTED MANUSCRIPT unigenes, 11,403 were successfully annotated into 26 specific protein function definitions
206
or orthologous categories (Fig 3B). Among these protein function categories, the largest
207
four were general function prediction only (5753); translation, ribosomal structure and
208
biogenesis (3605); transcription (2684); and replication, recombination and repair (2340).
209
To identify the biological processes by which the annotated unigenes are involved in,
210
KEGG pathway analysis was further carried out. In the KEGG analysis, 18,599 unigenes
211
were annotated into 6 major KEGG pathways, including cellular processes, environmental
212
information processing, genetic information processing, human diseases, metabolism and
213
organismal systems. These annotated unigenes were further divided into 42 level 2
214
subcategories pathways. Except the global map which had no pictorial information, the
215
largest subcategory group, infection diseases (Bacterial), had 2609 annotated genes,
216
followed by signal transduction (2400), immune system (2089) and translation (2060) (Fig
217
3C).
218
Table S2.
SC
M AN U
Apart from these, 299 level 3 KEGG subcategories were annotated and are listed in
TE D
219
RI PT
205
3.3. Classification and analysis of the differentially expressed genes (DEGs)
221
To analyze and characterize the DEGs following the different treatments, a cutoff of
222
FDR (false discovery rate) ≤ 0.05 and log2 Ratio ≥ 1 was used as the threshold for selecting
223
DEGs. Based on this, a large number of genes were observed to be dysregulated in VA
224
group (1988) than in VN group (743) relative to control group (NS), while 1064 genes were
225
differentially expressed between VA group and VN group (Fig. 4A and Table S3-S6). To
226
further analyze the DEGs between the AHPND V. parahaemolyticus and non-AHPND V.
227
parahaemolyticus samples, DEGs were further annotated with GO and KEGG databases. In
228
the GO enrichment analysis, a total of 345 DEGs were successfully enriched in biological
229
process (231), cellular component (154) and molecular function (281) (Fig 4B). For the
AC C
EP
220
9
ACCEPTED MANUSCRIPT KEGG pathway enrichment analysis, 709 of the DEGs were annotated, relative to all
231
unigenes and between VA and VN (Fig 4C, Table S7). When VA was compared to VN, 20
232
subcategories of the KEGG pathway were enriched with the highly significant biological
233
pathways being carbohydrate metabolism, xenobiotics biodegradation and metabolism,
234
infectious diseases (bacterial), lipid metabolism and amino acid metabolism (Fig 4C (i)).
235
On the other hand, 22 subcategories were enriched with lower but significantly changes in
236
signal transduction, infectious diseases (viral), cancers (specific types), immune system and
237
nervous system (Fig 4C (ii)).
SC
RI PT
230
M AN U
238
3.4. Validation of RNA-seq transcriptome data by qRT-PCR
240
In order to substantiate the transcriptome data as well as better understand the response of
241
shrimp hemocytes to AHPND V. parahaemolyticus and non-AHPND V. parahaemolyticus,
242
validation of the transcriptome data using qRT-PCR was performed. As mentioned above,
243
the AHPND V. parahemolyticus can induced massive sloughing of hepatopancreas tubule
244
epithelial cells [6], therefore to further explore the detailed hemocytes response induced by
245
AHPND V. parahemolyticus, DEGs were categorized with KEGG pathway enrichment
246
analyses in Fig. 4C. We found that the genes involved in cell growth and anti-apoptosis
247
were strongly induced in VA group compared with VN group in the transcriptome (Table
248
S8). Thus, six genes, src, iap2, cas2, cytochrome P450, gst and cytochrome c oxidase, were
249
selected for validation using qPCR (Fig 5). The qPCR results showed a significant
250
up-regulation in the expression of the oncogene src as well as the anti-apoptosis genes iap2
251
and cas2 in VA relative to VN and NS. Similarly, a significant up-regulation was observed
252
in the expression of the oxidoreductase gst, cytochrome P450 and cytochrome c oxidase in
253
VA relative to VN or NS. Thus, the qPCR results were generally consistent with the
254 255
transcriptome data.
AC C
EP
TE D
239
10
ACCEPTED MANUSCRIPT
4. Discussion
257
AHPND has caused severe losses in the shrimp industry in the last couple of years [5]. To
258
curtail this loss, several different control or preventive methods including the use of
259
inactivated bacteria antigen, IgY feed additives, phage therapy and nano-bubble technology
260
[18] have been carried out, but with little or no success. This enigma drove scientist to
261
comprehensively explore this disease. In this study, we have mimicked the late stage of
262
AHPND with intramuscularly injection of two strains of V. parahemolyticus, which allow a
263
direct contact between the bacteria and hemocytes. Following this, a comparative
264
transcriptomic analysis was carried out. Our results revealed that shrimp hemocytes
265
probably had a strong anti-apoptosis response against VA infection.
M AN U
SC
RI PT
256
266
As the most important immunocytes in shrimps, hemocytes are the major guardians against
268
pathogen invasion in the circulatory system [19]. The morphology and function of
269
hemocytes is reported to be affected by various bacteria components, as LPS injection has
270
been shown to cause significant decrease in total hemocyte count (THC) in Penaeus
271
monodon [20]. Similarly, the exotoxin from Clostridium difficile has been reported to
272
induce cytoskeletal reorganization in lepidopteran pests [21]. In the present study, our
273
results have shown that VA could induce stronger responses than VN, as indicated by the
274
higher cumulative mortality (Fig 1) as well as the higher number of differentially expressed
275
genes (DEGs) between different V. parahemolyticus treated hemocytes (Fig 4). To further
276
explore the detail response of the hemocytes, DEGs were categorized using KEGG pathway
277
analysis and revealed that anti-apoptosis-related pathways were more strongly activated in
278
VA. Apoptosis is a genetically programmed cell death mechanism, which could be induced
AC C
EP
TE D
267
11
ACCEPTED MANUSCRIPT by a wide range of stimuli, via the caspase cascade [22]. Our results suggest that VA could
280
not only induce hepatopancreas necrosis, but also cause hemocytes apoptosis, which greatly
281
impaired shrimp immune system. Thus, some anti-apoptosis genes like the oncogene src
282
kinase, which plays a key role in cell growth, division, migration, and survival signaling
283
pathways [23], were highly expressed in VA than in VN. Our results reveals how
284
hemocytes resist the stress exerted by AHPND V. parahemolyticus challenge in terms of
285
transcriptional response.
SC
RI PT
279
286
In conclusion, this comparative transcriptomic study depicts how shrimp hemocytes
288
respond to AHPND V. parahemolyticus infection. Future studies will explore how the
289
shrimp immune system can be prompted to respond to AHPND strain of bacteria and to
290
eventually find treatment for AHPND.
TE D
291
M AN U
287
Disclosure statement
293
The authors declare no conflicts of interest.
294
EP
292
Acknowledgments
296
This work was sponsored by National Natural Science Foundation of China (Nos.
297
31672689 & 31372558), Natural Science Foundation of Guangdong Province (No.
298
2017A030311032), Shantou University Scientific Research Foundation for Talents (No.
299
NTF16001) and Guangdong SAIL Foundation for Distinguished Scholars (No. 14600703).
AC C
295
12
ACCEPTED MANUSCRIPT 300 301 302
References [1] S. Interpreters, FAO (Food and Agriculture Organization of the United Nations), International Organization 1 (2014) 121-123. [2] S. Thitamadee, A. Prachumwat, J. Srisala, et al., P. Jaroenlak , P. V. Salachan, K.
304
Sritunyalucksana, et al., Review of current disease threats for cultivated penaeid shrimp
305
in Asia, Aquaculture, 452 (2016) 69-87.
RI PT
303
[3] L. Tran, L. Nunan, R.M. Redman, L.L. Mohney, C.R. Pantoja, K. Fitzsimmons, et al.,
307
Determination of the infectious nature of the agent of acute hepatopancreatic necrosis
308
syndrome affecting penaeid shrimp, Dis. Aquat. Organ. 105(1) (2013) 45-55.
SC
306
[4] Y.T. Yang, I.T. Chen, C.T. Lee, C.Y. Chen, S.S. Lin, L.I. Hor, et al., Draft genome
310
sequences of four strains of Vibrio parahaemolyticus, three of which cause early
311
mortality syndrome/acute hepatopancreatic necrosis disease in shrimp in China and
312
Thailand, Genome Announc. 2 (5) (2014) e00816-14.
314
[5] D.V. Lightner, R.M. Redman, C.R. Pantoja, B.L. Noble, L. Tran, Early mortality
TE D
313
M AN U
309
syndrome affects shrimp in Asia, Global Aquaculture Advocate. 15 (1) (2012) 40. [6] C.T. Lee, I.T. Chen, Y.T. Yang, T.P. Ko, Y.T. Huang, J.Y. Huang, et al., The opportunistic
316
marine pathogen Vibrio parahaemolyticus becomes virulent by acquiring a plasmid that
317
expresses a deadly toxin, Proc Natl Acad Sci U. S. A. 112 (34) (2015) 10798-10803.
EP
315
[7] W. Soonthornchai, S. Chaiyapechara, P. Jarayabhand, K. Soderhall, P. Jiravanichpaisal,
319
Interaction of Vibrio spp. with the inner surface of the digestive tract of Penaeus
320
monodon, PLoS ONE. 10 (8) (2015) e0135783.
AC C
318
321
[8] H.C. Lai, T.H. Ng, M. Ando, C.T. Lee, I.T. Chen, J.C. Chuang, et al., Pathogenesis of
322
acute hepatopancreatic necrosis disease (AHPND) in shrimp, Fish Shellfish Immunol. 47
323
(2) (2015) 1006-1014.
324
[9] W. Soonthornchai, S. Chaiyapechara, S. Klinbunga, W. Thongda, S. Tangphatsornruang,
325
T. Yoocha, et al., Differentially expressed transcripts in stomach of Penaeus monodon in 13
ACCEPTED MANUSCRIPT 326
response to AHPND infection, Dev. Comp. Immunol. 65 (2016) 53-63.
327
[10] M.G. Grabherr, B.J. Haas, M. Yassour, J.Z. Levin, D.A. Thompson, I. Amit, et al.,
328
Full-length transcriptome assembly from RNA-Seq data without a reference genome, Nat.
329
Biotechnol. 29 (7) (2011) 644-652. [11] A. Conesa, S. Gotz, J.M. Garcia-Gomez, J. Terol, M. Talon, M. Robles, Blast2GO: a
331
universal tool for annotation, visualization and analysis in functional genomics research,
332
Bioinformatics. 18 (2005) 3674-3676.
334
[12] B. Langmead, S.L. Salzberg, Fast gapped-read alignment with Bowtie 2, Nat. Methods.
SC
333
RI PT
330
9 (4) (2012) 357-359.
[13] A. Mortazavi, B.A. Williams, K. McCue, L. Schaeffer, B. Wold, Mapping and
336
quantifying mammalian transcriptomes by RNA-Seq, Nat. Methods 5 (7) (2008)
337
621-628.
340 341 342 343
false discovery rate controlling procedures, Bioinformatics. 19 (3) (2003) 368-375.
TE D
339
[14] A. Reiner, D. Yekutieli, Y. Benjamini, Identifying differentially expressed genes using
[15] M.D. Young, M.J. Wakefield, G.K. Smyth, A. Oshlack, Gene ontology analysis for RNA-seq: accounting for selection bias, Genome Biol. 11 (2) (2010) R14. [16] J. Ye, L. Fang, H. Zheng, Y. Zhang, J. Chen, Z. Zhang, et al., WEGO: a web tool for
EP
338
M AN U
335
plotting GO annotations, Nucleic Acids Res. 34 (2) (2006) 293-297. [17] Y. Moriya, M. Itoh, S. Okuda, A.C. Yoshizawa, M. Kanehisa, KAAS: an automatic
345
genome annotation and pathway reconstruction server, Nucleic Acids Res. 35 (2) (2007)
346
182-185.
AC C
344
347
[18] I. Hirono, S. Tinwongger, Y. Nochiri, H. Kondo, Latest research on acute
348
hepatopancreatic necrosis disease (AHPND) of penaeid shrimps, Philippines.(2016)
349
3-10.
350
[19] L. Vazquez, J. Alpuche, G. Maldonado, C. Agundis, A. Pereyra-Morales, E. Zenteno,
351
Review: Immunity mechanisms in crustaceans, Innate Immunol. 15 (3) (2009) 179-188. 14
ACCEPTED MANUSCRIPT 352
[20] J.A. Xian, X.X. Zhang, H. Guo, D.M. Wang, A.L. Wang, Cellular responses of the tiger
353
shrimp Penaeus monodon haemocytes after lipopolysaccharide injection, Fish Shellfish
354
Immunol 54 (2016) 385-390. [21] A. Castagnola, S.P. Stock, Common virulence factors and tissue targets of
356
entomopathogenic bacteria for biological control of Lepidopteran Pests, Insects. 5 (1)
357
(2014) 139-166.
360 361
cell survival: a review, Biochim. Biophys. Acta. 1813 (1) (2011) 238-259.
[23] R.J. Roskoski, Src protein-tyrosine kinase structure, mechanism, and small molecule inhibitors, Pharmacol Res. 94 (2015) 9-25.
362 363 364
EP AC C
367
TE D
365 366
SC
359
[22] L. Portt, G. Norman, C. Clapp, M. Greenwood, M.T. Greenwood, Anti-apoptosis and
M AN U
358
RI PT
355
15
ACCEPTED MANUSCRIPT
Figure legends
369
Fig. 1. Cumulative mortality graphs for VA, VN, and NS injected shrimps. Forty shrimps
370
were injected with VA (AHPND V. parahemolyticus) or VN (non-AHPND V.
371
parahemolyticus), at the same dosage of 104cfu per shrimp or with NS (Sterile normal
372
saline). The mortality tests were carried out in triplicates for each sample. Significant
373
difference is indicated by asterisks, **p < 0.01.
374
Fig. 2. Length distribution of all assembled unigenes.
SC
375
RI PT
368
376
Fig. 3. Functional enrichment of assembled genes. (A) Gene ontology (GO) classifications
378
of non-redundant unigenes. All annotated unigenes were categorized into 3 categories: A:
379
biological process, B: cellular component, C: molecular function. (B) The cluster of
380
orthologous groups (COG) classification of putative proteins. In total, 11,403 unigenes
381
were successfully annotated to 26 categories. (C) The Kyoto Encyclopedia of Genes and
382
Genomes(KEGG) classification of non-redundant unigenes. In total, 18,599 unigenes were
383
assigned to 6 special KEGG pathways, including Cellular Processes (a), Environmental
384
Information Processing (b), Genetic Information Processing (c), Human Diseases (d),
385
Metabolism (e), Organismal Systems (f).
TE D
EP
AC C
386
M AN U
377
387
Fig. 4. Analysis of differentially expressed genes between VN (non-AHPND V.
388
parahemolyticus) and VA (AHPND V. parahemolyticus) treated groups. (A) Comparison of
389
the differentially expressed genes (DEGs) between groups. (B) GO terms functional
390
enrichment analysis for VA and VN. (C) KEGG pathway enrichment analysis and
391
annotation compared between all unigenes and DEGs between VA and VN, (a) KEGG
392
pathway enrichment analysis between VA and VN with percentages higher than all
393
unigenes; (b). KEGG pathway enrichment analysis between VA and VN with percentages 16
ACCEPTED MANUSCRIPT 394
lower than all unigenes.
395
Fig. 5. Validation of differentially expressed genes (DEGs) by qPCR. To validate the
397
RNA-seq data, the relative gene expression levels of 6 selected DEGs were examined by
398
qPCR. EF1a was used as an internal control. Different letters on the graph denotes
399
significant statistical difference (p<0.05).
AC C
EP
TE D
M AN U
SC
RI PT
396
17
ACCEPTED MANUSCRIPT 400 401
Table 1. Summary of RNA-seq data for L. vannamei injected with AHPND and
402
non-AHPND V. parahemolyticus Samples
Total Raw Reads
Total Clean Reads
Q20 Percentage
GC Percentage
47,123,094
46,275,446
97.27%
47.52%
VN
47,941,818
47,076,076
97.71%
47.59%
NS
50,096,942
49,084,770
97.77%
47.31%
Total
145,161,854
142,436,292
97.58%
47.47%
RI PT
VA
Abbreviations: NS，normal saline treated group；VA，AHPND V. parahemolyticus treated group；VN，
404
non-AHPND V. parahemolyticus treated group.
AC C
EP
TE D
M AN U
SC
403
18
ACCEPTED MANUSCRIPT 405
Table 2. Summary of De novo assemble of L. vannamei hemocyte transcriptome Sample
Total Number
Total Length (nt)
Mean Length (nt)
N50
VA
53,879
48,513,536
900
2084
VN
53,401
49,225,356
922
2143
NS
53,910
48,998,100
All
45,658
60,455,885
RI PT
406
2131
1324
2574
AC C
EP
TE D
M AN U
SC
407
909
19
ACCEPTED MANUSCRIPT 408 409
Table 3. Functional annotation of unigenes from L. vannamei hemocytes transcriptome Annotation mode
NR
NT
Swiss-Prot
KEGG
COG
GO
ALL
23,223
14,516
21,041
18,599
11,403
11,457
26,634
50.86%
31.79%
46.08%
40.74%
Number of Unigenes with
Percentage of all assembled unigenes (45,658)
24.97%
AC C
EP
TE D
M AN U
SC
410
RI PT
hits or assignments
20
25.09%
58.33%
ACCEPTED MANUSCRIPT
Fig. 1
AC C
EP
TE D
M AN U
SC
RI PT
411
412 21
ACCEPTED MANUSCRIPT
Fig. 2
AC C
EP
TE D
M AN U
SC
RI PT
413
414 415 22
Fig. 3
417
A.
AC C
EP
TE D
M AN U
SC
416
RI PT
ACCEPTED MANUSCRIPT
418 23
ACCEPTED MANUSCRIPT
B.
AC C
EP
TE D
M AN U
SC
RI PT
419
420
24
ACCEPTED MANUSCRIPT
C.
AC C
EP
TE D
M AN U
SC
RI PT
421
422
25
Fig. 4
424
A.
AC C
EP
TE D
M AN U
SC
423
RI PT
ACCEPTED MANUSCRIPT
425 26
ACCEPTED MANUSCRIPT
B.
427
AC C
EP
TE D
M AN U
SC
RI PT
426
428
27
ACCEPTED MANUSCRIPT
C(i).
430
.
AC C
EP
TE D
M AN U
SC
RI PT
429
28
ACCEPTED MANUSCRIPT
C (ii).
AC C
EP
TE D
M AN U
SC
RI PT
431
432 433 29
ACCEPTED MANUSCRIPT
Fig. 5
AC C
EP
TE D
M AN U
SC
RI PT
434
435 30
ACCEPTED MANUSCRIPT
Research highlights ►
We found 1064 genes differentially expressed between the AHPND and Non-AHPND
V. parahemolyticus infected hemocytes in Litopenaeus vannamei. this, the genes involved in the major metabolism processes, cell growth and
RI PT
► Among
anti-apoptosis were further identified as the major responsive genes in response to the
AC C
EP
TE D
M AN U
SC
AHPND V. parahemolyticus.