Comparative transcriptome analysis of Sinonovacula constricta in gills and hepatopancreas in response to Vibrio parahaemolyticus infection

Accepted Manuscript Comparative transcriptome analysis of Sinonovacula constricta in gills and hepatopancreas in response to Vibrio parahaemolyticus infection Xuelin Zhao, Xuemei Duan, Zhenhui Wang, Weiwei Zhang, Ye Li, Chunhua Jin, Jinbo Xiong, Chenghua Li PII:

S1050-4648(17)30363-7

DOI:

10.1016/j.fsi.2017.06.040

Reference:

YFSIM 4661

To appear in:

Fish and Shellfish Immunology

Received Date: 11 April 2017 Revised Date:

14 June 2017

Accepted Date: 15 June 2017

Please cite this article as: Zhao X, Duan X, Wang Z, Zhang W, Li Y, Jin C, Xiong J, Li C, Comparative transcriptome analysis of Sinonovacula constricta in gills and hepatopancreas in response to Vibrio parahaemolyticus infection, Fish and Shellfish Immunology (2017), doi: 10.1016/j.fsi.2017.06.040. 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.

ACCEPTED MANUSCRIPT 1

Comparative transcriptome analysis of Sinonovacula constricta

2

in gills and hepatopancreas in response to Vibrio

3

parahaemolyticus infection

RI PT

4

Xuelin Zhao, Xuemei Duan, Zhenhui Wang, Weiwei Zhang, Ye Li, Chunhua Jin,

6

Jinbo Xiong, Chenghua Li*,

7

School of Marine Sciences, Ningbo University, Ningbo 315211, PR China

SC

5

M AN U

8 9 10 11

TE D

12

* Corresponding author:

14

Chenghua Li

15

818 Fenghua Road,

16

Ningbo University,

17

Ningbo, Zhejiang Province 315211, P. R. China

18

Tel: 86-574-87608368，Fax: 86-574-87608368,

19

Email: [email protected]

AC C

EP

13

20

1

ACCEPTED MANUSCRIPT 1

Abstract

3

The razor clam Sinonovacula constricta is an important economic species in China.

4

However, bacterial pathogenic diseases limits S. constricta farming industry for

5

large-scale production. In this study, de novo transcriptome sequencing was

6

performed on S. constricta gills and hepatopancreas under Vibrio parahaemolyticus

7

challenge for 12 h and 48 h, respectively. Transcripts assembly constructed 18,330

8

sequences, each of which was 500 bp long and functionally annotated, and 1,781 and

9

490 transcripts were differentially expressed in the gills and hepatopancreas,

10

respectively. Host immune factors that respond to Vibrio infection were then

11

identified. These factors included up-regulated transcripts with function in non-self

12

recognition, signal transduction, immune effectors and anti-apoptosis. The

13

comparison between the differentially expressed transcripts of the gills and

14

hepatopancreas indicated that immune responses had tissue specificity. As an

15

important external barrier between the environment and the clam, ATP-binding

16

cassette transporters and other ion transporters contribute to immune response in gills,

17

while, transcripts in complement system, such as complement 1 q protein,

18

IgGFc-binding protein, and low affinity immunoglobulin epsilon Fc receptor, were

19

more active in hepatopancreas and often not expressed in gill tissues. Eleven genes

20

were selected to be validated by qRT-PCR and the expressions were consistent with

21

the results of RNA-seq. Our study is the first attempt to identify molecular features in

22

different tissues of S. constricta in response to V. parahaemolyticus infection. These

23

findings improved our understanding of bivalve immunity and defense mechanisms

24

and revealed more potential immune-related genes.

SC

M AN U

TE D

EP

AC C

25

RI PT

2

26

Keywords: Sinonovacula constricta, Vibrio parahaemolyticus, Immune response,

27

RNA-seq

2

ACCEPTED MANUSCRIPT 28

1. Introduction Sinonovacula constricta (Lamarck 1818), a bivalve mollusk of the family

30

Solenidae, is a benthic clam distributed in intertidal zones and estuarine regions in

31

China, Japan and Korea. As one of the four major shellfish species cultivated in China,

32

S. constricta has been cultured for more than 500 years through natural means

33

particularly in Fujian and Zhejiang Provinces. The yield of cultured S. constricta is

34

more than 700,000 tons in 2014, which accounted for 30% of the total mudflat

35

shellfish production in China [1]. However, large-scale S. constricta production

36

through high-density culturing techniques is threatened by degradation of germplasm

37

resources, water pollution and pathogens [2].

SC

RI PT

29

Repeated episodes of mortality caused by bacterial infections occurring during

39

culture of bivalve mollusks reduce production and cause high economic losses [3].

40

Similar to other bivalves, S. constricta normally accumulates rich bacterial microbiota

41

owing to its filter-feeding habit. Vibrio sp. is a Gram-negative halophilic bacterium

42

found abundantly in marine and estuarine environments [4, 5] and is the main causal

43

agents of disease affecting all life stages of bivalves [6-10]. Particular, Vibrio

44

parahaemolyticus cause the mortality of S. constricta [11].

TE D

M AN U

38

S. constricta defends against pathogen invasion mainly through its innate

46

immunity because it lacks an adaptive immune system. Recently, studies have been

47

performed to analyze the molecular mechanisms of bivalve pathogen interactions and

48

were able to discover bivalve immune system components, such as Toll-like receptors

49

(TLRs), apoptosis pathways, and antimicrobial proteins [12-14]. In S. constricta, only

50

a few immune-relevant genes, such as ferritin [15], cathepsin [16, 17], have been

51

separately cloned and characterized. Moreover, only 43 putative immune genes were

52

identified through sequencing and conducting bioinformatics analysis of ESTs [2].

53

Thus, our knowledge of the immune system of S. constricta and different signaling

54

pathways implicated in its immune response remains incomplete.

AC C

EP

45

55

Next generation sequencing (NGS) provides large sequencing datasets in a fast

56

and cost-effective manner, even to non-model organisms. In bivalves, the genome of

57

Crassostrea gigas [18] and draft genome of Pinctada fucata [19] have been 3

ACCEPTED MANUSCRIPT constructed through this technology. Besides, NGS can provide insights into the

59

mechanisms of various biological process in bivalves, including environmental

60

stresses [20, 21], immunity [22, 23], and sex differentiation [24, 25]. Given that S.

61

constricta is an important economic bivalve species, understanding its biology and

62

genome resources is valuable. That is, understanding the molecular mechanisms by

63

which S. constricta adapts to variable environments is necessary for its breeding and

64

conservation. However, despite their importance, studies on these mechanisms are

65

scarce. To date, transcriptome sequencing for S. constricta has only focused on its

66

larval metamorphosis [26] and cadmium detoxification [27].

SC

RI PT

58

Our study aims to characterize the gene regulation features of S. constricta during

68

V. parahaemolyticus infection by profiling the transcripts of the gills and

69

hepatopancreas tissues in infected and non-infected S. constricta at two time points.

70

This study allows for the identification of potential immune genes involved in the

71

interactions between clams and Vibrio sp., as well as provides a theoretical basis for

72

the selective breeding of clam strains resistant to Vibrio sp..

73

2. Materials and Methods

74

2.1. Sample collection

TE D

M AN U

67

The procedure of sample collection was showed in Fig.1. In brief, 36 healthy S.

76

constricta were obtained from a commercial shellfish farm (Ninghai, Zhejiang, China)

77

in October 2014. The samples had an average body weight of 10 g, average body

78

length of 5.5 cm and average shell width of 2.25 cm. The samples were randomly

79

distributed in six tanks and maintained in 10 L of open-circuit-filtered seawater with a

80

salinity level of 20 at 16 °C with aeration. After accommodation for 3 days, three

81

tanks were exposed to V. parahaemolyticus with the final concentration of 107

82

CFU/mL, while the other tanks were used as controls. A sample was randomly

83

selected from each tank at 12 h and another at 48 h. These time points were selected

84

according to the previous work on immune related genes from clam responding to

85

pathogen appearing changes at 12 h and showing significant gene expression at 48 h

AC C

EP

75

4

ACCEPTED MANUSCRIPT [28]. Afterward, the gills and hepatopancreas of the selected samples were dissected

87

and immediately stored with 1 mL of RNAlater (Life technology, Waltham, MA, USA)

88

for later use.

89

2.2. RNA extraction and quantification

RI PT

86

Tissue samples (gills and hepatopancreas) from six individuals with the same

91

treatments were homogeneously mixed to 100 mg. Each mixed sample from eight

92

groups was lysed in 1 mL of TRIzol reagent (Invitrogen, Carlsbad, CA, USA) for

93

RNA extraction in accordance with the manufacturer’s instructions. The total RNA

94

concentration was then measured on a NanoDrop spectrophotometer (Thermo

95

Scientific, Waltham, MA, USA) and quality-assessed using by an Agilent 2100

96

BioAnalyzer (Agilent Technologies, Santa Clara, CA, USA). The total RNA was

97

subsequently treated with Dnase I (Ambion, Thermo Scientific, Waltham, MA, USA)

98

following the manufacturer’s protocol. MicroPoly(A)-Purist™ Kit (Ambion, Thermo

99

Scientific, Waltham, MA, USA) was used to enrich the Poly (A) mRNA from each

M AN U

SC

90

total RNA sample according to manufacturer’s instructions.

101

2.3. Next generation sequencing of transcriptome

TE D

100

A single-end fragment library was constructed following the SOLiD Total

103

RNA-seq Kit protocol (Life Technologies, PN4452437). The mRNA was fragmented

104

by RNase III and purified using a RiboMinus concentration module (Invitrogen,

105

Carlsbad, CA, USA). A hybridization master mix (SOLiD Total RNA-Seq Kit) was

106

then used to link the RNA fragment with the adaptors before reverse transcription.

107

The purified cDNA was size-selected through electrophoresis on a Novex TBE-Urea

108

Gel (Invitrogen, Carlsbad, CA, USA) at 180 V for 20 min, as described by Wang et al.

109

[29]. The cDNA fragments with 150-250 nt were precisely excised and used as

110

amplification template. Approximately 100 ng of cDNA was amplified using a SOLiD

111

Total RNA-Seq Kit. Emulsion PCR and bead enrichment were performed using a

112

SOLiD EZ BeadTM system (Life Technologies, Waltham, MA, USA). The enriched

AC C

EP

102

5

ACCEPTED MANUSCRIPT beads for each sample were then deposited on the sequencing slide. Finally, the

114

libraries were sequenced by the SOLiD 4 sequencing platform and color-space reads

115

were obtained.

116

2.4. Transcriptome analysis

RI PT

113

Data were transformed into sequence data after the initial images ware obtained.

118

The raw reads were trimmed by removing the adapter sequences and ambiguous or

119

low-quality reads (i.e., the proportion of low-quality bases with Q < 5 more than 59%)

120

using cutadapt [30] and Trimmomatic [31], to obtain clean data. The de novo

121

transcriptome assembly was performed using the Trinity program, a short read

122

assembler [32]. The transcripts with more than 500 bp were selected as reference

123

transcriptome. The predicted protein-coding sequences were searched against

124

Swiss-Prot

125

(ftp://ftp.ncbi.nih.gov/blast/db/), Pfam (http://pfam.xfam.org/) by using Blastx with an

126

E-value < 1e-5. Gene names were assigned to each protein sequence on the basis of

127

the best BLAST hit from all BLAST results.

128

2.5. Identification of differentially expressed genes

M AN U

SC

117

NCBI

non-redundant

(Nr)

protein

TE D

(http://www.uniprot.org/),

The clean reads from the two tissues obtained at 12 h and 48 h were mapped back

130

to the transcriptome assembly by using Tophat2 [33] and Bowtie software [34] with

131

default parameters. Differential expression analysis between the control and treated

132

groups of the gills and hepatopancreas at different time points were estimated using

133

Cufflinks and Cuffdiff [35]. The threshold for defining significant differentially

134

expressed (DE) transcripts between two different conditions was set as adjusted

135

p-value (q-value) smaller than 0.05 and absolute log2(foldchange) values greater than

136

1. Only DE transcripts with annotation were considered as candidates of interest for

137

further analysis.

138

2.6. GO and KEGG analysis

AC C

EP

129

6

ACCEPTED MANUSCRIPT 139

For functional annotation, Blast2GO [36] was used to search reference

140

transcriptome against the Gene ontology (GO) terms through sequence similarity blast.

141

The predicted proteins were scanned using the online Kyoto Encyclopedia of Genes

142

and

143

http://www.genome.jp/kegg/kaas/) using single-directional best-hit method (SBH) to

144

obtain KO terms and KEGG pathways. The significantly over-enriched GO and

145

KEGG pathways of the DE genes were calculated using a hypergeometric distribution

146

algorithm and using GOstat [37] and GSEABase packages [38] with a p-value cutoff

147

of 0.05.

148

2.7 Real-time quantitative PCR

(KEGG)

Automatic

Annotation

Server

(KAAS,

M AN U

SC

RI PT

Genomes

For the validation of the RNA-seq expression level, quantitative real-time PCR

150

(qRT-PCR) was conducted for relative quantification of mRNA using SYBR Green

151

real-time PCR kit (Takara) with ABI 7500. The qRT-PCR primers were designed with

152

Primer Premier 5 software, primers with amplification efficiency 95% - 105% by

153

means of standard curve drawn using dilution series. The 20 µl PCR mixture

154

contained 10µl of SYBR Premix Ex Taq II (TaKaRa), 0.8 µl of 10 µM forward primer,

155

0.8 µl of 10µM reverse primer, 0.4 µl of ROX reference dye II, 2 µl of cDNA template

156

and 6 µl of DEPC treated water. The qRT-PCR cycling conditions were 95 °C for 10 s,

157

60 °C for 20 s, and 72 °C for 30 s. At the end of the PCR cycles, melting curve

158

analyses were performed, and β-actin was used as an internal control. Three samples

159

in each group were used to reduce experimental error. The primers were listed in

160

Table S1.

EP

AC C

161

TE D

149

All data were analyzed using 2-∆∆Ct method and the values obtained represented

162

the n-fold difference relative to the control (untreated samples). The data were

163

presented as the relative expression levels (means ± SD, n = 3), Student’s t test was

164

performed to determine the differences between the treated and control groups.

165

Significant differences between the treated and corresponding control groups at each

166

time point were indicated as *P < 0.05. The error bars in the graphs represent the

167

standard error of the mean. 7

ACCEPTED MANUSCRIPT 168

3. Results and discussion

169

3.1. Library sequencing and de novo transcriptome assembly In filter-feeding invertebrates, gills provide first defense in mollusk immunity,

171

and has abundant circulating hemocytes in its filaments [39]. Alternatively, the

172

hepatopancreas is an important organ involved in immune response and oxidative and

173

heat stress of molluscan and crustacean [40]. The present study aims to identify the

174

molecular features of S. constricta under V. parahaemolyticus infection and compare

175

the immune-related genes involved in gills and hepatopancreas at different time points.

176

Owing to the absence of S. constricta genome information, we constructed a reference

177

transcriptome by performing the high-throughput sequencing and de novo assembly. A

178

total of 7,205,607 and 10,301,319 raw reads in gills were obtained from the control

179

groups at 12 and 48 h, respectively, and 13,765,736 and 11,586,049 raw reads

180

obtained from the Vibrio-treated groups at 12 and 48 h, respectively. In

181

hepatopancreas, a total of 7,311,534 and 8,406,059 raw reads were generated from the

182

control groups at 12 and 48 h, respectively, while 10,578,725 and 8,369,232 raw reads

183

were obtained from the treated groups at 12 and 48 h, respectively (Table 1). The

184

short read sequences were deposited into the short read archive of NCBI under the

185

accession number SRP102404. After the low-quality reads were trimmed and filtered,

186

nearly half of the raw data, as clean reads, was used for the de novo assembly of the

187

reference transcriptome on the basis of all sequenced RNA libraries to maximize

188

transcripts diversity. A total of 127,005 contigs ranging from 201 to 7,620 bp were

189

generated using the Trinity program. The size distribution of all the de novo

190

assembled transcripts is shown in Fig. 2A. The size range of 200-500 was dominated

191

in all the assembled transcripts, and the length distribution was the same with others

192

clams [41]. One of reasons that many transcripts were not assembled into full-length

193

sequences is because of the insufficient sequencing coverage. All eight libraries were

194

remapped to the 127,005 transcripts, and about 70% of the reads were remapped. The

195

remapped data was used for further different expression analysis.

AC C

EP

TE D

M AN U

SC

RI PT

170

8

ACCEPTED MANUSCRIPT 196

3.2. Functional annotation The unigenes with length ≥ 500 were annotated on the basis of the identity of the

198

translation frame, and 18,330 transcripts were predicted as putative protein-coding

199

sequences. The putative protein-coding transcripts were searched using Blastx against

200

NCBI nr database and Swissprot database (E-value < 1e-5) (Table 2). As expected, the

201

top three species that had the most similarity to S. constricta sequences were mollusks,

202

such as the Pacific oyster (Crassostrea gigas, 7,605), owl limpet (Lottia gigantean,

203

3,358), and the California sea slug (Aplysia californica, 2,088), which have available

204

genomes (Fig. 2B). Nearly half of the S. constricta transcripts were not annotated, and

205

the results of the annotation were similar to those of the annotation performed in

206

previous transcriptomic studies on bivalves [24, 40-42]. These results may be

207

attributed to the transcripts spanning untranslated mRNA regions, or transcripts

208

containing only nonconserved protein domains [41].

M AN U

SC

RI PT

197

GO assignments were widely used to classify gene functions in terms of

210

biological process, molecular function and cellular component. On the basis of

211

sequence similarity (E-value of 1e-5), 6,173 transcripts were assigned to at least one

212

GO term (Table 2). As summarized in Fig. 3, a total of 1,620, 383 and 881 terms in

213

the three main categories were annotated, respectively. The most dominant terms

214

presented in the three main categories were “cellular process”, “metabolic process”,

215

“catalytic activity”, “binding”, “cell part” and “macromolecular complex” (Fig. 3).

216

Some transcripts were clustered into the immune-related subcategories of response to

217

stimulus (597), immune system process (40), and biological adhesion (47). Such

218

transcripts may be involved in S. constricta defense and resistance toward V.

219

parahaemolyticus infection. In addition, the sequences were searched through KEGG

220

database for systematic analysis [43]. A total of 9,314 transcripts were mapped to 262

221

pathways (Fig. 4). Of the 262 predicted pathways, signal transduction was the

222

predominant group, including 3,372 unigenes. Moreover, a total of 1,232 transcripts

223

were grouped into the immune system subcategory and divided into 19

224

immune-related pathways. Of these subcategories, NOD-like receptor signaling

AC C

EP

TE D

209

9

ACCEPTED MANUSCRIPT pathway included the largest number of transcripts (110), followed by leukocyte

226

transendothelial migration (90). These immune-related pathways can help us elucidate

227

the gene functions and interactions involved in response to Vibrio infection.

228

3.3. Identification of differentially expressed transcripts

RI PT

225

The result of transcriptome assembly was used as a reference for the analysis of

230

global gene expression in two tissues (gills and hepatopancreas) at two time-points

231

(12 h and 48 h) in response to V. parahaemolyticus infection. Overall, 1,781 DE

232

unigenes were detected from the pair-wise comparisons (|log2(foldchange)| > 2,

233

q-value < 0.05) among the gill samples (Fig. 5A), and 490 DE unigenes among the

234

hepatopancreas samples (Fig. 5B). The DE unigenes in each tissue were the union of

235

the DE unigenes between the control group and the treated group at 12 h, those

236

between the control group and the treated group at 48 h, and those between the treated

237

group at 12 h and the treated group at 48 h. These DE unigenes had functional

238

annotations and were further examined for their putative functions involved in

239

immune response during Vibrio infection. Annotated transcripts were subsequently

240

examined, and assigned GO terms and KO terms as references for enrichment

241

analysis. Significant enriched GO terms and KEGG pathways were identified via the

242

Fisher’s exact test (P < 0.05)

243

3.4. Immune response of gill towards vibrio challenge

EP

TE D

M AN U

SC

229

The gill is the first barrier in the body of S. constricta and contains many

245

filaments that increase its superficial area to exchange materials between environment

246

and the body. According to the GO enrichment results of DE transcripts in gills (Fig.

247

5C, Table S2), GO terms related to cell motility were enriched significantly. These

248

GO terms included microtubule-based movement (GO: 0007018, FDR = 1.7E-13),

249

cilium or flagellum-dependent cell motility (GO: 0001539, FDR = 3.02E-09), and

250

locomotion (GO: 0040011, FDR = 6.29E-04), which corresponded to the

251

characteristic of the gill tissues. In addition, obsolete ATP catabolic process (GO:

252

0006200, FDR = 4.26E-12) and single-organism cellular process (GO: 0044763, FDR

AC C

244

10

ACCEPTED MANUSCRIPT = 0.01) were enriched as well, indicating that cellular processes correlated to energy

254

supplied by ATP were active during bacterial infection. The three major enriched

255

KEGG pathways included Spliceosome (KO03040), RNA transport (KO03013), and

256

RNA degradation (KO03018), all of which are related to gene transcription (Table 3).

257

Notably, phagosome (KO04145) was enriched and directly correlated to immune

258

response.

RI PT

253

Innate immune responses provide unique host defenses against invading

260

pathogens in invertebrate. Pattern-recognition receptors (PRRs), the recognition basis

261

of innate immune system, can recognize microbial components known as

262

pathogen-associated molecular patterns (PAMPs), which are special structures

263

expressed on the cell surfaces of pathogens [44]. The hemocytes in gill filaments can

264

sense stimuli through an array of cell surface receptors, as well as trigger immune

265

signaling pathways that result in specific immune responses [45], such as

266

phagocytosis, reactive oxygen species (ROS) production, and secretion of immune

267

effectors and cytokines [46, 47]. DE transcripts involved in immune response listed in

268

Table 4 were identified, indicating that comprehensive host-pathogen interactions

269

existed in gills. In addition, the gill, as the gateway of bivalves, contains many

270

transmembrane proteins and ion channels on its cells. Organic transporters and ion

271

transporters were identified in DE transcripts, which may contribute to immune

272

responses.

EP

TE D

M AN U

SC

259

Among the PRRs in the DE transcripts, C-type lectins (CTLs) and scavenger

274

receptors (SRs) were the two major receptor categories in response to bacterial

275

challenge. The macrophage mannose receptor 1 (MMR1) is an important phagocytic

276

receptor mediating the binding and ingestion of microorganisms with surface

277

mannose residues and soluble mannose-containing glycoproteins [48]. The

278

up-regulated expression of MMR1 with the time suggested that active hemocyte

279

phagocytosis induced by Vibrio enhanced with the prolonged infection time. Perlucin

280

and sialic acid binding lectin (SABL) are all the members of the lectin family and can

281

play a role in non-self-antigen recognition to trigger immune response through glycan

282

recognition [49, 50]. Perlucin is present in Manila clams and hard clams during

AC C

273

11

ACCEPTED MANUSCRIPT microbial infection in previous studies [41, 49]. SRs share common function of

284

recognizing oxidized or acetylated low-density lipoproteins (LDL) [51]. In our results,

285

LDL receptor-related protein 6 (LRP6) expression was up-regulated at 12 h, but the

286

expression levels of SRs only showed up-expressed comparing 12 h and 48 h

287

treatments. These results indicated that LRP6 has a quick response to Vibrio infection.

288

The expression levels of integrins showed up-expressed at 48 h compared with 12 h

289

treatments, similar to those of SRs. Invertebrate integrins are responsible for cellular

290

adhesive processes correlated to several immune responses [52]. The expression of

291

SRs and integrins may suggest that these receptors did not specifically identify of

292

Vibrio, but played roles in immune response.

SC

RI PT

283

Cells have evolved many mechanisms to modulate the signaling pathways

294

involved in many regulators, such as phosphatases and kinases, ubiquitin-related

295

proteins, and transcription factors, to balance activation and suppression of

296

PRR-mediated innate immune responses [53]. During S. constricta response to V.

297

parahaemolyticus, a variety of transcripts encoding kinases and phosphatases were

298

up-regulated (Table S2), suggesting the involvement of MAPK and other

299

kinase-mediated cascades in the regulation of signaling pathways to mediate immune

300

response. In S. constricta, the genes in Rho signaling pathways are up-regulated

301

during an immune response. Rho signaling pathways are reported to regulate the

302

apoptosis process [54]. Calmodulin, as regulator of calcium-regulated pathways, is

303

down-regulated in the gills of S. constricta. This observation is consistent with the

304

result of the previous study, which reported that the down regulation of calmodulin

305

suppresses calcium-regulated pathways under pathogen infection [41]. In our results,

306

E3 ubiquitin-protein ligase was up-regulated compared with 12 h and 48 h treatments,

307

suggesting the activation of autophagy pathway for the degradation of intracellular

308

pathogens [55]. The ubiquitin ligases in the proteasome-ubiquitin pathway are

309

regulated by the presence of bacterial colonization [56]. It is not surprising that the

310

molecule elements of ubiquitin pathway were up-regulated at 48 h than those at 12 h,

311

owing to the bacterial reproduction. These above signaling pathways and signaling

312

molecules

AC C

EP

TE D

M AN U

293

were as

receptor-ligand

interactions

in

regulating the cellular 12

ACCEPTED MANUSCRIPT 313

behaviors/activities in innate immune system of bivalves [57]. Heat shock proteins (HSPs) are highly generated when induced by stress and play

315

important roles in protein folding, the protection of proteins from denaturation or

316

aggregation [58]. HSPs respond to Vibrio infection in clam [59] and shrimp [60]. In

317

the present study, HSPs were highly expressed in S. constricta subjected to V.

318

parahaemolyticus challenge. This finding is in line with those of previous reports.

319

Hypoxia-inducible factors (HIFs) are genes that play a role in adaptation to hypoxia in

320

animals [61]. Host inflammatory cells creates a localized low oxygen environment

321

during inflammation or infection and thus promote the up-expression of HIFs [13].

322

Increased HIF expression was expected in our data, as HIFs respond to hypoxia stress

323

under pathogen infection.

M AN U

SC

RI PT

314

As immune effectors, genes of complement system, antioxidant defense system,

325

lysozyme and related cytokines are differentially expressed during Vibrio infection in

326

this study. High expression levels of these genes were also observed in other bivalves

327

after Vibrio challenge [13, 41, 62]. Stress conditions, including bacterial infection lead

328

to ROS accumulation [63]. Differentially expressed oxidases in this study suggested

329

the need for the host to timely balance out excessive ROS and other toxic

330

intermediates produced under bacterial infection. Genes related to complement

331

systems, including a-macroglobulin complement component family protein,

332

complement component C3 and MAC/perforin- and kringle-domain-containing

333

protein, and lysozyme were up-regulated at 48 h treatment (Table 4) and thus play a

334

pivotal role in killing or clearing pathogens [64, 65].

EP

AC C

335

TE D

324

Apoptosis is an essential host defense process against bacterium, as it removes

336

damaged and infected cells without causing inflammatory destructions around the

337

dying cell [66, 67]. In our results, baculoviral IAP repeat-containing proteins, which

338

are apoptosis suppressors, directly inhibited the caspase activity [68], were

339

up-expression. Interferon regulatory factor (IRF) family encodes transcription factors

340

that play important roles in immune defense, stress responses, reproduction and

341

development [69]. In our results, interferon-induced helicase C domain-containing

342

protein 1 and IRF2 were up-regulated in S. constricta during Vibrio infection. IRF2 is 13

ACCEPTED MANUSCRIPT involved in the immune response to LPS and polyinosinic-polycytidylic acid in pearl

344

oyster and in the regulation of apoptosis [70], suggesting it is a multifunctional

345

transcription factor responding to bacterial challenges [71]. They may share similar

346

function in clams by preventing gills from death under Vibrio infection. This

347

hypothesis is supported by the up-regulation of integrins, which were shown to

348

protect cells from apoptosis and induce anti-apoptosis pathways [41].

RI PT

343

Interestingly, genes encoding transporters including organic or ion transporters,

350

are differentially expressed significantly. ATP binding cassette (ABC) transporters can

351

pump and export the modified xenobiotics out of the cell [72]. In mussels, ABCB and

352

ABCC can control the entry of nutrients and xenobiotics from the external

353

environment and mainly act as a first line of cellular defense [73]. ABCA1 and

354

ABCG1 negatively regulate signaling by TLRs to inhibit inflammatory response [74].

355

In our results, most ABC transporters showed under-regulated expression, suggesting

356

that they contribute to the negative regulation of immune response and reduce the

357

import and export of substrates in S. constricta. To date, there was little report that

358

found ABC transporters play a role in bivalves’ immune response. Furthermore,

359

several metal ion transporters, including plasma membrane calcium ATPase and

360

copper-transporting ATPase 1 (Table 4), were up-expressed in gills in response to

361

bacteria. Meanwhile, calcium acts as a second messenger in numerous cell types, and

362

can be induced to influx from the extracellular space into the immune cells for Ca2+

363

signaling [75]. Copper-transporting ATPase and other ion transporters contribute to

364

host defense by controlling the supply of essential micronutrients in the vicinity of

365

infection sites to reduce parasite survival [76].

366

3.5. Immune response of hepatopancreas towards vibrio challenge

AC C

EP

TE D

M AN U

SC

349

367

A total of 1,781 DE transcripts were obtained in the transcriptomes of gills, but

368

only nearly a quarter of this value was obtained in the transcriptomes of

369

hepatopancreas (490). This results may be attributed to tissue specificity [77]. The

370

hepatopancreas was considered as an important organ involved in immune response

371

and oxidative and heat stress of molluscan and crustacean. However, studies that 14

ACCEPTED MANUSCRIPT compare the difference between tissues at transcriptomic level in invertebrate are few.

373

Our results indicated that the gill may have more important roles than hepatopancreas

374

in immune response against bacteria. The GO enrichment analysis (Table S3) results

375

showed that regulation of cyclin-dependent protein serine/threonine kinase activity

376

(GO: 0000079), regulation of cyclin-dependent protein kinase activity (GO: 1904029),

377

and cytoskeleton-dependent cytokinesis (GO: 0061640) were the first three terms,

378

which were related to cell cycle regulation. In addition, no KEGG pathway was

379

enriched. According to the DE transcripts associated with immune response, a suite of

380

transcripts up-expressed in gills was not observed in hepatopancreas (Table 5), such

381

transcripts included genes in immune signaling pathways, transporters, PRRs and

382

immune effectors.

M AN U

SC

RI PT

372

Meanwhile, C-type lectin was the only up-regulated category of PRRs observed

384

in hepatopancreas. The up-expression of immune effectors, such as complement

385

systems, enzymes related to oxidoreduction and lysozymes, were observed in

386

hepatopancreas and as well as observed in gills. However, more DE transcripts that

387

contribute to the complement system were found in hepatopancreas than in gills,

388

including complement 1 q protein (C1q), IgGFc-binding protein, and low affinity

389

immunoglobulin epsilon Fc receptor (FCER) and so on. The presence of such

390

transcripts possibly contributed to the killing of pathogen in hepatopancreas. In

391

addition, arginine kinase (AK), a phosphagen kinase in the invertebrate energy

392

metabolism, plays a defensive role against viral and bacterial infection [78, 79], and

393

observed in the hepatopancreas under Vibrio challenge.

394

3.6. Validation of differential expressed genes by qRT-PCR

EP

AC C

395

TE D

383

Six and five unigenes related to immune response were selected in gills and

396

hepatopancreas respectively. These unigenes were selected from the DE transcripts in

397

gills and hepatopancreas. These DE transcripts have functions in immune response or

398

are rarely reported in other bivalves. Their expression levels are listed in Tables 4 and

399

5. These transcripts were validated on the basis of their differential expression profiles

400

under Vibrio infection at two time-points. qRT-PCR was performed for the validation. 15

ACCEPTED MANUSCRIPT The correlation between the expression values of qRT-PCR and RNA-seq are showed

402

in Fig. 6. These transcripts included cytochrome c oxidase subunit III (CYC), MAPK,

403

LRP6, MMR1, multidrug resistance protein 1A-like (ABCB1) and multidrug

404

resistance-associated protein 5 (ABCC5) in gills and AR, complement C1q tumor

405

necrosis factor-related protein 4 (C1qtnf4), suppressor of tumorigenicity protein 14

406

(ST14), C1q and FCER in hepatopancreas. As shown in Fig. 6, the expression trends

407

of selective genes by qRT-PCR were basically consistent to RNA-seq analysis results

408

(Tables 4 and 5). In addition, we verified that ABCC5 were only expressed in gills

409

and ST14, C1q and FCER only expressed in the hepatopancreas, indicating that

410

different immune responses occur in gills and hepatopancreas.

411

4. Conclusions

M AN U

SC

RI PT

401

In the present study, we constructed and sequenced the gill and hepatopancreas

413

transcriptomes of S. constricta under V. parahaemolyticus infection at two time-points

414

to understand the mechanisms of immune response in different tissues. The

415

differential expression analysis revealed significant differences in the gene expression.

416

In particular, 1,781 DE transcripts were obtained in the gills and 490 DE transcripts in

417

the hepatopancreas. By comparing the DE transcripts related to immune response of

418

gills with those of hepatopancreas, we found that gill tissues had active responses in

419

pathogen recognition, signal pathways and inhibition of apoptosis. In hepatopancreas,

420

transcripts in the complement systems were up-expressed and several transcripts were

421

expressed specifically. In addition, several ABC transporters and ion transporters

422

might be the specific immune-related markers in gills. Overall, this study provides an

423

opportunity to explore different immune defense mechanisms against V.

424

parahaemolyticus in different tissues of S. constricta and thus may lay the foundation

425

for the selective breeding of disease-resistance S. constricta.

AC C

EP

TE D

412

426 427 428 16

ACCEPTED MANUSCRIPT 429

Acknowledgments This work was financially supported by Zhejiang Major Program of Science and

431

Technology (2016C02055-9, 2015C32004), Natural Science Foundation of Ningbo

432

(2015C10009,2015C10062), Start Research Fund projects for Xuelin Zhao, and the

433

K.C. Wong Magna Fund in Ningbo University.

435

SC

434

RI PT

430

Supplementary data

437

Table S1. Primer information of qRT-PCR.

438

Table S2. GO enrichment analysis of the differentially expressed genes in gills.

439

Table S3. GO enrichment analysis of the differentially expressed genes in

440

hepatopancreas.

444 445 446 447 448 449

TE D

443

EP

442

AC C

441

M AN U

436

450 451 452 453 454 17

ACCEPTED MANUSCRIPT References

456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498

[1] X. Yuan, W. Zhao, China Fishery Statistical Yearbook, Fisheries Agency of China Agriculture Ministry, (2015). [2] B. Feng, L. Dong, D. Niu, S. Meng, B. Zhang, D. Liu, et al, Identification of immune genes of the Agamaki clam (Sinonovacula constricta) by sequencing and bioinformatic analysis of ESTs, Mar. Biotechnol. 12(3) (2010) 282-291. [3] R. Beaz-Hidalgo, S. Balboa, J.L. Romalde, M.J. Figueras, Diversity and pathogenecity of Vibrio species in cultured bivalve molluscs, Environ. Microbiol. Rep. 2(1) (2010) 34-43. [4] A. Ramesh, B.G. Loganathan, K. Venkateswaran, Ecological dynamics of marine luminous bacteria, J. Basic Microbiol. 30(9) (1990) 689-703. [5] F.L. Thompson, T. Iida, J. Swings, Biodiversity of vibrios, Microbiol. Mol. Biol. Rev. 68(3) (2004) 403-431. [6] P.C. Liu, Y.C. Chen, C.Y. Huang, K.K. Lee, Virulence of Vibrio parahaemolyticus isolated from cultured small abalone, Haliotis diversicolor supertexta, with withering syndrome, Lett. Appl. Microbiol. 31(6) (2000) 433-437. [7] A.B.Â.n. Casandra, L.Â.r.P. Marcial Leonardo, S.B. Ricardo, Effect of Vibrio alginolyticus on larval survival of the blue mussel Mytilus galloprovincialis, Dis. Aquat. Organ. 59(2) (2004) 119-123. [8] M.E. Robyn, S.F. Carolyn, A.E. Ralph, P.H. Russell, Pathogenicity testing of shellfish hatchery bacterial isolates on Pacific oyster Crassostrea gigas larvae, Dis. Aquat. Organ. 58(2-3) (2004) 223-230. [9] C. Paillard, S. Gausson, J.L. Nicolas, J.P. le Pennec, D. Haras, Molecular identification of Vibrio tapetis, the causative agent of the brown ring disease of Ruditapes philippinarum, Aquaculture 253(1–4) (2006) 25-38. [10] J. Meng, L. Zhang, B. Huang, L. Li, G. Zhang, Comparative analysis of oyster (Crassostrea gigas) immune responses under challenge by different Vibrio strains and conditions, Molluscan Research 35(1) (2015) 1-11. [11] Q. Qiu, L. Xie, X. Fu, S. Yang. Isolation and identification of Vibrio Parahaemolyticus from Sinonovacula constricta. Anhui Agri. Sci. Bull. 2010 16:46-7. [12] G. Zhang, L. Li, J. Meng, H. Qi, T. Qu, F. Xu, L. Zhang, Molecular basis for adaptation of oysters to stressful marine intertidal environments, Annual Review of Animal Biosciences 4(1) (2016) 357-381. [13] N.G. Ertl, W.A. O’Connor, A. Papanicolaou, A.N. Wiegand, A. Elizur, Transcriptome analysis of the Sydney rock oyster, Saccostrea glomerata: insights into molluscan immunity, PLoS ONE 11(6) (2016) e0156649. [14] A. Tassanakajon, K. Somboonwiwat, P. Amparyup, Sequence diversity and evolution of antimicrobial peptides in invertebrates, Dev. Comp. Immunol. 48(2) (2015) 324-341. [15] C. Li, H. Li, X. Su, T. Li, Identification and characterization of a clam ferritin from Sinonovacula constricta, Fish Shellfish Immunol. 30(4–5) (2011) 1147-1151. [16] D. Niu, K. Jin, L. Wang, B. Feng, J. Li, Molecular characterization and expression analysis of four cathepsin L genes in the razor clam, Sinonovacula constricta, Fish Shellfish Immunol. 35(2) (2013) 581-588. [17] D. Niu, K. Jin, L. Wang, F. Sun, J. Li, Identification of cathepsin B in the razor clam Sinonovacula constricta and its role in innate immune responses, Dev. Comp. Immunol. 41(1)

AC C

EP

TE D

M AN U

SC

RI PT

455

18

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

(2013) 94-99. [18] G. Zhang, X. Fang, X. Guo, L. Li, R. Luo, F. Xu, et al, The oyster genome reveals stress adaptation and complexity of shell formation, Nature 490(7418) (2012) 49-54. [19] T. Takeuchi, T. Kawashima, R. Koyanagi, F. Gyoja, M. Tanaka, T. Ikuta, et al, Draft genome of the pearl oyster Pinctada fucata: a platform for understanding bivalve biology, DNA Res. (2012). [20] X. Zhao, H. Yu, L. Kong, S. Liu, Q. Li, Comparative transcriptome analysis of two oysters, Crassostrea gigas and Crassostrea hongkongensis provides insights into adaptation to hypo-osmotic conditions, PLoS ONE 9(11) (2014) e111915. [21] H.J. Lim, B.M. Kim, I.J. Hwang, J.S. Lee, I.Y. Choi, Y.J. Kim, et al, Thermal stress induces a distinct transcriptome profile in the Pacific oyster Crassostrea gigas, Comp. Biochem. Physiol. D 19 (2016) 62-70. [22] E.E.R. Philipp, L. Kraemer, F. Melzner, A.J. Poustka, S. Thieme, U. Findeisen, et al, Massively parallel RNA sequencing identifies a complex immune gene repertoire in the lophotrochozoan Mytilus edulis, PLoS ONE 7(3) (2012) e33091. [23] L. Zhang, L. Li, Y. Zhu, G. Zhang, X. Guo, Transcriptome analysis reveals a rich gene set related to innate immunity in the eastern oyster (Crassostrea virginica), Mar. Biotechnol. 16(1) (2014) 17-33. [24] V. Teaniniuraitemoana, A. Huvet, P. Levy, C. Klopp, E. Lhuillier, N. Gaertner-Mazouni, et al, Gonad transcriptome analysis of pearl oyster Pinctada margaritifera: identification of potential sex differentiation and sex determining genes, BMC Genomics 15(1) (2014) 1-20. [25] N. Zhang, F. Xu, X. Guo, Genomic analysis of the Pacific oyster (Crassostrea gigas) reveals possible conservation of vertebrate sex determination in a mollusc, G3: Genes|Genomes|Genetics 4(11) (2014) 2207-2217. [26] D. Niu, F. Wang, S. Xie, F. Sun, Z. Wang, M. Peng, J. Li, Developmental transcriptome analysis and identification of genes involved in larval metamorphosis of the razor clam, Sinonovacula constricta, Mar. Biotechnol. 18(2) (2016) 168-175. [27] Z. Wang, Y. Shao, C. Li, W. Zhang, X. Duan, X. Zhao, et al, RNA-seq analysis revealed ROS-mediated related genes involved in cadmium detoxification in the razor clam Sinonovacula constricta, Fish Shellfish Immunol. 57 (2016) 350-361. [28] R. Moreira, A. Romero, M. Milan, L. Bargelloni, B. Novoa, A. Figueras. Gene expression profile analysis of Manila clam (Ruditapes philippinarum) hemocytes after a Vibrio alginolyticus or Perkinsus olseni challenge using an immune-enriched oligo-microarray. (2015) 21-27. [29] L. Wang, X. Sun, Z. Zhou, T. Zhang, Q. Yi, R. Liu, et al, The promotion of cytoskeleton integration and redox in the haemocyte of shrimp Litopenaeus vannamei after the successive stimulation of recombinant VP28, Dev. Comp. Immunol. 45(1) (2014) 123-132. [30] M. Martin. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17 (2014) 10-2. [31] A.M. Bolger, M. Lohse, B. Usadel. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. (2014) btu170. [32] M.G. Grabherr, B.J. Haas, M. Yassour, J.Z. Levin, D.A. Thompson, I. Amit, et al, Full-length transcriptome assembly from RNA-Seq data without a reference genome, Nat Biotech. 29(7) (2011) 644-652.

AC C

499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542

19

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

[33] D. Kim, G. Pertea, C. Trapnell, H. Pimentel, R. Kelley, S.L. Salzberg, TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions, Genome Biol. 14(4) (2013) R36. [34] B. Langmead, Aligning short sequencing reads with Bowtie, Current protocols in bioinformatics (2010) 11.7. 1-11.7. 14. [35] C. Trapnell, A. Roberts, L. Goff, G. Pertea, D. Kim, D.R. Kelley, H. Pimentel, S.L. Salzberg, J.L. Rinn, L. Pachter, Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks, Nat. Protoc. 7(3) (2012) 562-578. [36] A. Conesa, S. Götz, J.M. García-Gómez, J. Terol, M. Talón, M. Robles, Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research, Bioinformatics 21(18) (2005) 3674-3676. [37] T. Beißbarth, T.P. Speed, GOstat: find statistically overrepresented Gene Ontologies within a group of genes, Bioinformatics 20(9) (2004) 1464-1465. [38] M. Morgan, S. Falcon, R. Gentleman, GSEABase: gene set enrichment data structures and methods, R package version 1.4.0. [39] J.K. Seo, J.M. Crawford, K.L. Stone, E.J. Noga, Purification of a novel arthropod defensin from the American oyster, Crassostrea virginica, Biochem. Bioph. Res. Co. 338(4) (2005) 1998-2004. [40] Y. Ren, J. Xue, H. Yang, B. Pan, W. Bu, Transcriptome analysis of Ruditapes philippinarum hepatopancreas provides insights into immune signaling pathways under Vibrio anguillarum infection, Fish Shellfish Immunol. 64 (2017) 14-23. [41] K. Wang, C. del Castillo, E. Corre, E.P. Espinosa, B. Allam, Clam focal and systemic immune responses to QPX infection revealed by RNA-seq technology, BMC Genomics 17(1) (2016) 146. [42] R.B. Leite, M. Milan, A. Coppe, S. Bortoluzzi, A. dos Anjos, R. Reinhardt, et al, mRNA-Seq and microarray development for the Grooved carpet shell clam, Ruditapes decussatus: a functional approach to unravel host-parasite interaction, BMC Genomics 14(1) (2013) 741. [43] M. Kanehisa, S. Goto, KEGG: kyoto encyclopedia of genes and genomes, Nucleic Acids Res. 28(1) (2000) 27-30. [44] S. Akira, S. Uematsu, O. Takeuchi, Pathogen recognition and innate immunity, Cell 124(4) (2006) 783-801. [45] O. Takeuchi, S. Akira, Pattern recognition receptors and inflammation, Cell 140(6) (2010) 805-820. [46] P. Soudant, F.-L.E. Chu, A. Volety, Host–parasite interactions: Marine bivalve molluscs and protozoan parasites, Perkinsus species, J. Invertebr. Pathol. 114(2) (2013) 196-216. [47] R. Hatanaka, Y. Sekine, T. Hayakawa, K. Takeda, H. Ichijo, Signaling pathways in invertebrate immune and stress response, ISJ 6 (2009) 32-43. [48] M. Stein, S. Keshav, N. Harris, S. Gordon, Interleukin 4 potently enhances murine macrophage mannose receptor activity: a marker of alternative immunologic macrophage activation, J. Exp. Med. 176(1) (1992) 287-292. [49] R. Moreira, M. Milan, P. Balseiro, A. Romero, M. Babbucci, A. Figueras, L. Bargelloni, B. Novoa, Gene expression profile analysis of Manila clam (Ruditapes philippinarum) hemocytes after a Vibrio alginolyticus challenge using an immune-enriched oligo-microarray, BMC Genomics 15(1) (2014) 267.

AC C

543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586

20

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

[50] C. Li, S. Yu, J. Zhao, X. Su, T. Li, Cloning and characterization of a sialic acid binding lectins (SABL) from Manila clam Venerupis philippinarum, Fish Shellfish Immunol. 30(4) (2011) 1202-1206. [51] S. Gordon, Pattern recognition receptors: doubling up for the innate immune response, Cell 111(7) (2002) 927-930. [52] K. Terahara, K.G. Takahashi, A. Nakamura, M. Osada, M. Yoda, T. Hiroi, et al, Differences in integrin-dependent phagocytosis among three hemocyte subpopulations of the Pacific oyster “Crassostrea gigas”, Dev. Comp. Immunol. 30(8) (2006) 667-683. [53] X. Cao, Self-regulation and cross-regulation of pattern-recognition receptor signalling in health and disease, Nat. Rev. Immunol. 16(1) (2016) 35-50. [54] F. Cencetti, C. Bernacchioni, F. Tonelli, E. Roberts, C. Donati, P. Bruni, TGFβ1 evokes myoblast apoptotic response via a novel signaling pathway involving S1P4 transactivation upstream of Rho-kinase-2 activation, FASEB J. 27(11) (2013) 4532-4546. [55] V. Rogov, V. Dötsch, T. Johansen, V. Kirkin, Interactions between autophagy receptors and ubiquitin-like proteins form the molecular basis for selective autophagy, Mol. Cell 53(2) (2014) 167-178. [56] C.K. Chun, J.V. Troll, I. Koroleva, B. Brown, L. Manzella, E. Snir, et al, Effects of colonization, luminescence, and autoinducer on host transcription during development of the squid-vibrio association. PNAS 105(32) (2008) 11323-11328. [57] J.E. Humphries, T.P. Yoshino. Cellular receptors and signal transduction in molluscan hemocytes: connections with the innate immune system of vertebrates. Integr. Comp. Biol. 43 (2003) 305-12. [58] I. Horváth, G. Multhoff, A. Sonnleitner, L. Vígh, Membrane-associated stress proteins: more than simply chaperones, BBA-Biomembranes 1778(7) (2008) 1653-1664. [59] Y. Bao, Q. Wang, H. Liu, Z. Lin, A small HSP gene of bloody clam (Tegillarca granosa) involved in the immune response against Vibrio parahaemolyticus and lipopolysaccharide, Fish Shellfish Immunol. 30(2) (2011) 729-733. [60] W. Rungrassamee, R. Leelatanawit, P. Jiravanichpaisal, S. Klinbunga, N. Karoonuthaisiri, Expression and distribution of three heat shock protein genes under heat shock stress and under exposure to Vibrio harveyi in Penaeus monodon, Dev. Comp. Immunol. 34(10) (2010) 1082-1089. [61] H. Piontkivska, J.S. Chung, A.V. Ivanina, E.P. Sokolov, S. Techa, I.M. Sokolova, Molecular characterization and mRNA expression of two key enzymes of hypoxia-sensing pathways in eastern oysters Crassostrea virginica (Gmelin): hypoxia-inducible factor α (HIF-α) and HIF-prolyl hydroxylase (PHD), Comp. Biochem. Physiol. D 6(2) (2011) 103-114. [62] J. De Lorgeril, R. Zenagui, R.D. Rosa, D. Piquemal, E. Bachère, Whole transcriptome profiling of successful immune response to Vibrio infections in the oyster Crassostrea gigas by digital gene expression analysis, PLoS One 6(8) (2011) e23142. [63] U. Bandyopadhyay, D. Das, R.K. Banerjee, Reactive oxygen species: oxidative damage and pathogenesis, CURRENT SCIENCE-BANGALORE- 77 (1999) 658-666. [64] H. Rus, C. Cudrici, F. Niculescu, The role of the complement system in innate immunity, Immunol. Res. 33(2) (2005) 103-112. [65] E. Bachère, Anti-infectious immune effectors in marine invertebrates: potential tools for disease control in larviculture, Aquaculture 227(1) (2003) 427-438.

AC C

587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630

21

ACCEPTED MANUSCRIPT

665 666

EP

TE D

M AN U

SC

RI PT

[66] K. Terahara, K.G. Takahashi, Mechanisms and immunological roles of apoptosis in molluscs, Curr. Pharm. Des. 14(2) (2008) 131-137. [67] T. Kiss, Apoptosis and its functional significance in molluscs, Apoptosis 15(3) (2010) 313-321. [68] L. Zhang, L. Li, G. Zhang, Gene discovery, comparative analysis and expression profile reveal the complexity of the Crassostrea gigas apoptosis system, Dev. Comp. Immunol. 35(5) (2011) 603-610. [69] J. Nehyba, R. Hrdličková, H.R. Bose, Dynamic evolution of immune system regulators: the history of the interferon regulatory factor family, Mol. Biol. Evol. 26(11) (2009) 2539-2550. [70] X.-D. Huang, W.-G. Liu, Q. Wang, M. Zhao, S.-Z. Wu, Y.-Y. Guan, Y. Shi, M.-X. He, Molecular characterization of interferon regulatory factor 2 (IRF-2) homolog in pearl oyster Pinctada fucata, Fish Shellfish Immunol. 34(5) (2013) 1279-1286. [71] T. Tamura, H. Yanai, D. Savitsky, T. Taniguchi. The IRF family transcription factors in immunity and oncogenesis. Annu. Rev. Immunol. 26 (2008) 535-584. [72] D. Schlenk, D.R. Buhler, Xenobiotic biotransformation in the Pacific oyster (Crassostrea gigas), Comp. Biochem. Phy. C 94(2) (1989) 469-475. [73] T. Luckenbach, D. Epel, ABCB-and ABCC-type transporters confer multixenobiotic resistance and form an environment-tissue barrier in bivalve gills, Am. J. Physiol-Reg. I. 294(6) (2008) R1919-R1929. [74] L. Lai, K.M. Azzam, W.C. Lin, P. Rai, J.M. Lowe, K.A. Gabor, et al, MicroRNA-33 regulates the innate immune response via ATP binding cassette transporter-mediated remodeling of membrane microdomains, J. Biol. Chem. 291(37) (2016) 19651-19660. [75] M. Vig, J.P. Kinet, Calcium signaling in immune cells, Nat. Immunol. 10(1) (2009) 21-27. [76] C.C. Chang, M.-S. Yeh, H.K. Lin, W. Cheng, The effect of Vibrio alginolyticus infection on caspase-3 expression and activity in white shrimp Litopenaeus vannamei, Fish Shellfish Immunol. 25(5) (2008) 672-678. [77] M.I. Mohd-Shamsudin, Y. Kang, Z. Lili, T.T. Tan, Q.B. Kwong, H. Liu, et al, In-depth tanscriptomic analysis on giant freshwater prawns, PloS One 8(5) (2013) e60839. [78] J. Arockiaraj, P. Vanaraja, S. Easwvaran, A. Singh, T. Alinejaid, R.Y. Othman, et al, Gene profiling and characterization of arginine kinase-1 (MrAK-1) from freshwater giant prawn (Macrobrachium rosenbergii), Fish Shellfish Immunol. 31(1) (2011) 81-89. [79] R. Rao, Y.B. Zhu, T. Alinejad, S. Tiruvayipati, K.L. Thong, J. Wang, et al, RNA-seq analysis of Macrobrachium rosenbergii hepatopancreas in response to Vibrio parahaemolyticus infection, Gut Pathogens 7(1) (2015) 6.

AC C

631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664

22

ACCEPTED MANUSCRIPT 667

Figure legends

669

Fig. 1. Diagram of sample collection.

670

Fig. 2. Summary of S. constricta de novo assembled transcriptome. A: Assembled

671

transcripts size distribution; B: Distribution of the top 10 species with most

672

homologues to S. constricta. Transcripts were searched using Blastx against NCBI nr

673

database with a cutoff value of E < 1e-05.

674

Fig. 3. GO classification of the assembled transcripts. All the annotated transcripts

675

were grouped into different functional sub-categories: cellular component (A),

676

molecular function (B) and biological process (C).

677

Fig. 4. KEGG annotation analysis of the assembled transcripts.

678

Fig. 5. Grouping and GO enrichment analysis of differentially expressed transcripts.

679

A: The grouping of the gene lists of G12-N, G48-N and G48-12; B: The grouping of

680

the gene lists of H12-N, H48-N and H48-12; C: GO enrichment analysis result (the

681

top three GO terms in subcategories) of the differentially expressed transcripts in gills;

682

D: GO enrichment analysis result (the top three GO terms in subcategories) of the

683

differentially expressed transcripts in hepatopancreas. G: gill tissues; H:

684

hepatopancreas tissues; 12-N: the differentially expressed transcripts between the

685

treated group and control group for 12 h; 48-N: the differentially expressed transcripts

686

between the treated group and control group for 48 h; 48-12: the differentially

687

expressed transcripts between the treated groups for 12 h and 48 h. Target set: the

688

proportion of the differentially expressed transcripts with a certain GO terms in all the

689

differentially expressed transcripts. Reference set: all the transcripts with a certain GO

690

terms in all the transcripts in the transcriptome of S. constricta.

691

Fig. 6. Validation of relative expression levels of eleven transcripts by qRT-PCR

692

compared with RNA-seq. A: Relative expression levels of six genes at 12 h in gills; B:

693

Relative expression levels of six genes at 48 h in gills; C: Relative expression levels

694

of five genes at 12 h in hepatopancreas; D: Relative expression levels of five genes at

695

48 h in hepatopancreas. T-RNA-seq: the foldchange between the treated groups with

696

the control group by RNA-seq.

AC C

EP

TE D

M AN U

SC

RI PT

668

23

ACCEPTED MANUSCRIPT

Table 1 Summary statistics of the transcriptome sequencing for S. constricta from control groups and treated groups. Gill controlled (48h)

Gill treated (48h)

Hepatopancreas controlled (12h)

Hepatopancreas treated (12h)

Hepatopancreas controlled (48h)

Hepatopancreas treated (48h)

7,205,607 3,291,760 45.68% 1,662,752 64.07%

13,765,736 5,670,928 41.20% 3138388 68.84%

10,301,319 4,346,584 42.19% 2,336,054 65.29%

11,586,049 4,996,778 43.13% 2,822,260 64.34%

7,311,534 3,618,758 49.49% 1,793,179 67.31%

10,578,725 5,612,089 53.05% 2,802,107 68.42%

8,406,059 3,476,846 41.36% 1,713,598 68.84%

8,369,232 3,313,028 39.59% 1,671,985 70.13%

EP

TE D

M AN U

SC

RI PT

Gill treated (12h)

AC C

Raw sequencing reads Clean reads Clean reads (%) Mapped reads Mapped reads (%)

Gill controlled (12h)

24

ACCEPTED MANUSCRIPT

Table 2 S. constricta transcriptome assembly and annotation overview Annotation results

Unigene (≥ 500 bp)

Transcripts (bp) Transcripts (≥ 500 bp) Smallest transcripts (bp) Largest transcripts (bp)

127,005 33,134 201 7,620

Nr anntotation Swissprot annotation Pfam hits GO annotation KEGG annotation All annotation

17,770 8,964 12,971 6,173 9,314 18,330

Percentage (%)

RI PT

Counts

53.6% 27.1% 39.1% 18.6% 28.1% 55.3%

AC C

EP

TE D

M AN U

SC

Assembly results

25

ACCEPTED MANUSCRIPT

Pathway ID

DGEs number

Spliceosome RNA transport RNA degradation Indole alkaloid biosynthesis Betalain biosynthesis AGE-RAGE signaling pathway in diabetic complications Amino sugar and nucleotide sugar metabolism Phagosome

ko03040 ko03013 ko03018 ko00901 ko00965 ko04933 ko00520 ko04145

29 26 15 3 3 8 11 14

Unigenes number

p-value

202 187 91 6 6 38 70 99

0.001989 0.005333 0.007165 0.009462 0.009462 0.011207 0.027972 0.032671

AC C

EP

TE D

M AN U

SC

Term Name

RI PT

Table 3 KEGG enrichment of differentially expressed genes in gills.

26

ACCEPTED MANUSCRIPT Table 4 List of the annotated transcripts related to immune response that were differentially expressed in gills during S. constricta response against V. parahaemolyticus. * stands for q-value < 0.05. “inf” designates an infinite fold change as the expression of that transcript in control group was equal to 0.

-1.48111 -0.7494

TE D

EP

AC C

G48-12 log2 (foldchange)

0.494531

3.34852*

-2.80316*

2.79151

1.94559 3.7503* -1.60412

2.68072* 1.66962 2.5683*

0.013308 -5.78764* 3.01341*

-0.83129 -2.51699

0.686374 -0.14313

2.86534* 4.05517*

2.12828*

-0.16484

-1.48204

-1.8568 -2.16016 -2.04973 -1.79389 -2.34073 1.28502

0.286374 -0.07756 0.714786 -0.18744 0.330866 1.59124*

2.79969* 2.9486* 3.1057* 3.26335* 3.63814* 0.999252

-1.63685

1.96685*

3.92935*

-1.04766

-1.41654*

1.2481

-2.76049 -1.22563 -4.12269 0.109073

0.6592 0.869291 0.116669 -1.69507*

3.28146* 2.9451* 5.03946* -1.84771

-2.38841 -0.8979 -2.86835

0.709285 0.362458 1.1089

3.5321* 2.48136* 4.8454*

M AN U

Pathogen recognition receptors (PRRS) Scavenger receptor UN073535 Lysyl oxidase-like protein 2 PREDICTED: deleted in malignant brain UN091254 tumors 1 protein-like C-type lectin UN098373 Macrophage mannose receptor 1 UN097586 Perlucin UN000273 sialic acid binding lectin Integrin UN079222 integrin UN071400 integrin beta subunit Others Low-density lipoprotein receptor-related UN017133 protein 6 Immune signaling and cell communication Phosphatase and kinases UN104438 Cadherin-23 UN016246 CD63 antigen UN081454 Contactin UN017857 Murinoglobulin-2 UN015111 Protein roadkill UN013131 Calcium-dependent protein kinase 31 Receptor-type tyrosine-protein phosphatase UN049721 kappa Receptor-type tyrosine-protein phosphatase UN018054 mu Rho signaling UN085181 Rho GTPase-activating protein 18 UN015554 Rho guanine nucleotide exchange factor 12 UN089657 Rho-associated protein kinase 2 UN089937 Rho-related GTP-binding protein RhoJ Ubiquitin pathway UN001965 E3 ubiquitin-protein ligase HUWE1 UN017946 E3 ubiquitin-protein ligase KCMF1 UN012639 E3 ubiquitin-protein ligase MIB2

G48-N log2 (foldchange)

RI PT

G12-N log2 (foldchange)

Annotation

SC

Transcripts ID

27

ACCEPTED MANUSCRIPT E3 ubiquitin-protein ligase TRIM33 E3 ubiquitin-protein ligase UBR3

-2.7268 -1.91084

0.481416 -0.08133

3.63163* 2.82568*

UN067183

E3 ubiquitin-protein ligase UBR4

-4.86687*

-0.15681

4.94897*

-0.74831 -1.2872 -1.70863 -2.35737 -3.42822*

-1.47661* 0.089658 1.50803* 0.265099 -1.40623

0.0775322 2.42812* 3.03804* 3.83944* -0.634327

-1.78173*

0.0855118

0.2355

3.22746*

-0.0862

SC

-2.34944

-1.13652

0.093206

2.93237*

-2.59383*

1.52445*

3.6964*

-1.23792 -3.21518 -1.97511

0.609134 1.62564* 0.056238

2.54006* 4.74071* 2.99195*

-2.30726

-0.01434

3.45841*

-2.36973 -1.95638 1.38678 -1.30874 -1.66183 -0.21364

1.181 0.679289 1.40829* -0.05542 0.596332 0.399995

4.3888* 3.26636* 0.060597 2.44138* 2.39578* 3.56662*

-1.58853

0.268785

3.21292*

0.996063

3.32094*

1.24598

-0.75808

2.40412*

4.07654*

-2.8195 -1.32916

1.60474 -1.76382*

3.71033* 1.57569

1.30456*

-0.087

-1.25885

-0.35167 0.718335

0.69988 -1.52295*

3.62753* -0.428219

AC C

EP

TE D

M AN U

UN102510 Tripartite motif-containing protein 2 UN089051 Tripartite motif-containing protein 3 UN017506 Tripartite motif-containing protein 56 UN018075 Ubiquitin-associated protein 2 UN016798 TPA: TPA_exp: ubiquitin protein ligase Calcium mediated signal transduction UN004065 calmodulin EF-hand calcium-binding UN074485 domain-containing protein 6 Other pathways Cyclic AMP-responsive element-binding UN008342 protein 3-like protein 2 MAP kinase-interacting UN008369 serine/threonine-protein kinase 1 UN007312 Elongation factor 2 kinase UN006658 Fibropellin-1 UN003343 Regulator of G-protein signaling 12 Endogenous ligands ATP-dependent Clp protease ATP-binding UN074297 subunit clpX UN079956 Heat shock 70 kDa protein 12A UN084396 Heat shock 70 kDa protein 12B UN096277 heat shock protein 23 UN100733 heat shock protein 60 UN085758 hypoxia inducible factor-1alpha UN011236 Hypoxia up-regulated protein 1 Immune effectors Complement system A-macroglobulin complement component UN077936 family protein UN076835 complement component C3 MAC/perforin- and UN010982 kringle-domains-containing protein, partial Antioxidant defense system UN093254 alkaline phosphatase UN102716 CYP450 family 4 cytochrome c oxidase subunit III UN005779 (mitochondrion) UN004935 cytochrome P450 30 UN101645 Cytochrome P450 4F8

RI PT

UN043376 UN017775

28

ACCEPTED MANUSCRIPT

0.503565 -0.12665

-0.08032 -0.09565 0.789628 -1.11522 -1.41611* -0.95406 -6.13664*

TE D

EP

AC C

2.58059 4.72402* 2.7237* 2.45386 -0.775475 -3.06872* -1.89625

RI PT

-3.17502* -4.19092 -0.6881 -0.96258 -0.4584 -0.5667 -2.26372

1.9124* 2.06651*

2.26378 2.79833*

0.063186

2.15935*

0.135977

3.86596*

-2.82491 -0.86246 -1.7619

0.69988 0.217734 -0.44622

3.62753* 2.51241* 2.67862*

-1.73234

2.39913*

5.54074*

-0.52047 -3.72931

0.57583 1.41106

3.24974* 4.7692*

-1.09063

0.189473

2.57393*

-1.5441

1.92852*

3.39403*

-3.96407

0.286376

5.22009*

-2.57634

1.13505

3.74422*

0.087465

1.09889

2.87254*

-2.71751

0.389393

4.05065*

-1.77974

0.723605

3.37858*

-0.23675

-3.70118*

-1.8575

-1.97959 -1.6077 -2.29667 -3.26741*

-0.22839 -1.47007* 0.469854 -1.96834*

2.73829* 1.28518 3.04169* 2.60611

-3.00087*

-1.63206*

2.06721

-0.56656

SC

-2.44646

M AN U

UN005037 D-aspartate oxidase UN061548 Dual oxidase UN001760 Dual oxidase maturation factor 1 UN105827 Eosinophil peroxidase UN004707 Glutathione S-transferase 1, isoform D UN095064 Glutathione S-transferase S1 UN100860 Inositol oxygenase Lysozyme UN004600 cathepsin L1 UN016610 Lysozyme 1 cytokines and cytokine receptors UN008514 Interleukin-17 receptor D Nuclear factor interleukin-3-regulated UN082790 protein Apoptosis UN011652 Apoptosis inhibitor 5 UN088780 Apoptosis-stimulating of p53 protein 1 UN087969 Baculoviral IAP repeat-containing protein 3 Baculoviral IAP repeat-containing protein UN013315 7-A UN000691 caspase-8 UN071978 Programmed cell death 6-interacting protein Interferon-induced helicase C UN091287 domain-containing protein 1 UN003916 interferon regulatory factor 2 Transporter Organic transporter UN008393 ABCB/p-glycoprotein-like protein ATP-binding cassette sub-family A member UN005829 1 ATP-binding cassette sub-family F member UN061986 2 ATP-binding cassette sub-family F member UN072505 3 ATP-binding cassette sub-family G member UN014383 1 ATP-binding cassette sub-family G member UN011332 2-like UN015387 Excitatory amino acid transporter 1 UN102363 Excitatory amino acid transporter 3 UN074622 Multidrug resistance-associated protein 1 UN105119 Multidrug resistance-associated protein 5 PREDICTED: multidrug UN105114 resistance-associated protein 5-like

29

ACCEPTED MANUSCRIPT

UN017426 UN067989 UN007582 UN054000 UN075821 UN080469

-2.64586* -0.6557

-0.993896 2.56277*

-3.12328*

0.96278

4.11871*

-1.56755

-0.2375

2.81388*

-0.02925

2.49901*

0.485314

3.61284*

-1.84081*

0.179116

-0.07745 2.59837* -0.20205 0.061556 -1.59835*

3.67279* 2.477 3.19955* 2.98494* 1.28904

-1.44299 -2.40095* -0.94492 -3.1031 -inf -2.03081 -0.84477 -0.17145

M AN U

UN105255

-0.11306 -2.01088

TE D

UN016503

3.94641*

EP

UN007514

Calcium-transporting ATPase sarcoplasmic/endoplasmic reticulum type Copper-transporting ATPase 1 Electroneutral sodium bicarbonate exchanger 1 Na+/K+ ATPase alpha subunit Neuronal acetylcholine receptor subunit beta-3 plasma membrane calcium ATPase Prestin Sodium-dependent multivitamin transporter Solute carrier family 12 member 4 Zinc transporter ZIP14

0.063784

AC C

UN082922

-2.73159*

RI PT

UN092226 UN017522 Ion transporter

PREDICTED: multidrug resistance protein 1A-like Solute carrier family 46 member 3 taurine transporter

SC

UN096702

30

ACCEPTED MANUSCRIPT Table 5 List of the annotated transcripts related to immune response that were differentially expressed in hepatopancreas during S. constricta response against V. parahaemolyticus. * stands for q-value < 0.05. “inf” designates an infinite fold change as the expression of that transcript in control group was equal to 0.

-1.8803 2.46709* 1.65598*

TE D

EP

AC C

L48-12 log2 (foldchange)

2.45467* 0.93062 0.068284

3.45372* 0.846791 -2.02542

-2.60428* 3.29646* -0.47797 -1.40754

0.615177 2.542* 0.90262 1.59157

1.96314 1.57482 2.20917* 2.81372*

0.004137 5.5823* 5.43788* 1.7603*

3.96865* 3.43111* 3.07173* 0.125869

2.94247* 0.587731 0.909273 0.519807

-2.33337*

-0.26071

2.62051*

3.7813*

0.182257

-3.24503*

-1.22884 -1.0759

1.89877* 3.13853*

1.92252* 2.67929*

1.02177

-5.53199*

-6.79595*

-3.6489*

0.827154

3.45969

0.169786

-1.04247

-2.7739*

1.81325*

0.58746

-0.36271

-1.75443*

-1.35834*

-0.15087

-5.69484*

-inf

-inf

1.62142* 2.64179*

0.031936 2.42962

-0.71651 -1.18563

M AN U

Pathogen recognition receptors (PRRS) C-type lectin UN018015 C-type lectin UN089291 C-type mannose receptor 2 UN102691 sialic acid-binding lectin Immune signaling and cell communication Ubiquitin pathway UN006914 E3 ubiquitin-protein ligase DZIP3 UN000736 Tripartite motif-containing protein 2 UN006861 Tripartite motif-containing protein 47 UN016216 Ubiquitin-protein ligase E3C Endogenous ligands UN081663 70 kDa heat shock protein UN093164 Heat shock 70 kDa protein 12A UN015766 Heat shock 70 kDa protein 12B UN100733 heat shock protein 60 Immune effectors Complement system UN016190 complement 1 q protein Complement C1q tumor necrosis factor-related UN102695 protein 4 UN091243 IgGFc-binding protein UN096803 Kappa-type opioid receptor Low affinity immunoglobulin epsilon Fc UN080616 receptor MAC/perforin- and kringle-domains-containing UN010982 protein, partial UN007456 macrophage expressed protein putative C1q domain containing protein UN082251 MgC1q41 putative C1q domain containing protein UN098877 MgC1q75 UN032053 Trace amine-associated receptor 5 Antioxidant defence system UN004935 cytochrome P450 30 UN105819 D-erythrulose reductase

L48-N log2 (foldchange)

RI PT

L12-N log2 (foldchange)

Annotation

SC

Transcripts ID

31

ACCEPTED MANUSCRIPT Dimethylaniline monooxygenase [N-oxide-forming] 5 Dual oxidase maturation factor 1 extracellular copper/zinc superoxide dismutase Peroxisomal acyl-coenzyme A oxidase 1

0.120379

-0.39919

-0.5439 2.17199 0.512006

1.58017 -inf 3.51487*

2.39274* -inf* 3.08098*

1.60941* -1.42946*

0.134547 1.2616*

-0.86367 1.76039*

1.26316* 1.62704

1.56609* 3.42769*

0.401711

1.30835*

-0.59964 -2.91249* -1.68878*

AC C

EP

TE D

M AN U

SC

UN001760 UN037283 UN008185 Lysozyme UN096611 cathepsin L1 UN088075 i-type lysozyme Apoptosis UN008178 Bcl-2 like 1 protein UN017526 Suppressor of tumorigenicity protein 14 other immune genes UN094400 arginine kinase

inf*

RI PT

UN083817

32

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

Highlights: De novo transcriptome sequencing was performed for S. constricta under Vibrio parahaemolyticus challenge for 12 h and 48 h in gills and hepatopancreas,

RI PT

respectively. A total of 1,781 and 490 transcripts were identified as differentially expressed in gills and hepatopancreas.

M AN U

indicated the tissue specific immune response.

SC

The differentially expressed transcripts between gills and hepatopancreas

AC C

EP

TE D

Some differential expressed genes were further validated by qRT-PCR