Integrated analysis of mRNA-miRNA expression in Tilapia infected with Tilapia lake virus (TiLV) and identifies primarily immuneresponse genes

Journal Pre-proof Integrated analysis of mRNA-miRNA expression in Tilapia infected with Tilapia lake virus (TiLV) and identifies primarily immuneresponse genes Yingying Wang, Qing Wang, Yingying Li, Jiyuan Yin, Yan Ren, Cunbin Shi, Sven M. Bergmann, Xinping Zhu, Weiwei Zeng PII:

S1050-4648(20)30049-8

DOI:

https://doi.org/10.1016/j.fsi.2020.01.041

Reference:

YFSIM 6777

To appear in:

Fish and Shellfish Immunology

Received Date: 18 October 2019 Revised Date:

27 December 2019

Accepted Date: 22 January 2020

Please cite this article as: Wang Y, Wang Q, Li Y, Yin J, Ren Y, Shi C, Bergmann SM, Zhu X, Zeng W, Integrated analysis of mRNA-miRNA expression in Tilapia infected with Tilapia lake virus (TiLV) and identifies primarily immuneresponse genes, Fish and Shellfish Immunology (2020), doi: https:// doi.org/10.1016/j.fsi.2020.01.041. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2020 Published by Elsevier Ltd.

1

Integrated analysis of mRNA-miRNA expression in Tilapia infected with Tilapia

2

Lake virus (TiLV) and identifies primarily immuneresponse genes

3

Yingying Wanga,b, #, Qing Wang a, #, Yingying Li a, Jiyuan Yin a, Yan Ren a, Cunbin Shi a

4

, Sven M. Bergmannc, Xinping Zhu a,b,*, Weiwei Zenga,*

5 6

a

7

Laboratory of Aquatic Animal Immune Technology of Guangdong Province, Pearl

8

River Fisheries Research Institute, Chinese Academy of Fishery Sciences, Guangzhou,

9

Guangdong, the PR China; 510380

Key Laboratory of Fishery Drug Development of Ministry of Agriculture, Key

10

b

11

China

12

c

13

Institute of Infectology, Südufer 10, 17493 Greifswald-Insel Riems, Germany

College of Fisheries and Life Science, Shanghai Ocean University, Shanghai 201306,

Friedrich-Loeffler-Institut (FLI), Federal Research Institute for Animal Health,

14 15

#

16

* To whom correspondence should be addressed:

17

Xinping Zhu, PhD, E-mail: [email protected]

18

Weiwei Zeng, PhD, E-mail: [email protected]

19 20 21 22 23 24 25 26 27 28 29 30

These authors contributed equally to this work.

31

Abstract

32

We investigated differential gene expression in Tilapia infected with the Tilapia

33

Lake virus (TiLV).We used high-throughput sequencing to identify mRNAs and

34

miRNAs involved in TiLV infection progression We identified 25,359 differentially

35

expressed genes that included 863 new genes. We identified 1770, 4142 and 4947

36

differently expressed genes comparing non-infected controls with 24 and 120 h

37

infections and between the infected groups, , respectively. These genes were enriched

38

to 291 GO terms and 62 KEGG pathways and included immune system progress and

39

virion genes. High-throughput miRNA sequencing identified 316 conserved miRNAs,

40

525 known miRNAs and 592 novel miRNAs. Furthermore, 138, 198 and 153

41

differently expressed miRNAs were found between the 3 groups listed above,

42

respectively. Target prediction revealed numerous genes including erythropoietin

43

isoform X2, double-stranded RNA-specific adenosine deaminase isoform X1, bone

44

morphogenetic protein 4 and tapasin-related protein that are involved in immune

45

responsiveness. Moreover, these target genes overlapped with differentially expressed

46

mRNAs obtained from RNA-seq. These target genes were significantly enriched to

47

GO terms and KEGG pathways including immune system progress, virion and Wnt

48

signaling pathways. Expression patterns of differentially expressed mRNA and

49

miRNAs were validated in 20 mRNA and 19 miRNAs by qRT-PCR. We also were

50

able to construct a miRNA-mRNA target network that can further understand the

51

molecular mechanisms on the pathogenesis of TiLV and guide future research in

52

developing effective agents and strategies to combat TiLV infections in Tilapia.

53

Keywords: tilapia, tilapia lake virus (TiLV), RNA-seq, miRNA-seq, integrative

54

interaction

55 56

1. Introduction

57

Tilapia is the second most abundant farmed freshwater fish worldwide (>3.5

58

million tons annually) because of its rapid growth and disease resistance [1-3]. Tilapia

59

are currently a critical protein source in developing countries and major producers are

60

China, Thailand, Indonesia, Egypt, Brazil, Vietnam, Bangladesh, the Philippines,

61

Colombia and Uganda [1]. Overall, there are few significant viral diseases of tilapia

62

and infections generally have limited impacts [4, 5]. However, the novel RNA virus

63

tilapia lake virus (TiLV) has been found to infect these fish causing extensive

64

mortality rates on farms in

65

Sequence analysis of TiLV demonstrated that this is a novel orthomyxo-like virus with

66

a 10-segment negative-sense RNA genome [6]. The first segment shares weak

67

sequence homology to influenza C virues while the other genomic segments show no

68

homology to known viruses [6]. Unfortunately, there are no available vaccines against

69

TiLV although these infections are increasing. Therefore, the identification of TiLV

70

genes that are activated during infections are necessary to reveal the underlying

71

mechanisms used by TiLV for infection.

Israel, Ecuador, Colombia, Egypt and Thailand

[2, 6-8].

72

MicroRNAs (miRNA) are endogenous 21-23 nucleotide (nt) fragments of

73

non-coding RNAs that act as post-transcriptional regulators of target mRNA stability

74

and translation [9, 10]. MiRNAs are wildly expressed in eukaryotic organisms as well

75

as viruses and play critical roles in development, metabolism, immune responses,

76

disease responses and virus infections [11-14]. These miRNAs have also been

77

identified in virus infections of fish. For example, miR-1, miR-30b, miR-150, and

78

miR-184 play important roles in red-spotted grouper nervous necrosis virus infections

79

in the orange-spotted grouper [15]. The miR-214 suppresses replication of the

80

rhabsovirus snakehead vesiculovirus by regulating host glycogen synthase [16]. A

81

recent study has showed that hemolymph exosome exosomal 157 differently

82

expressed miRNAs are related to whitespot syndrome virus infection in red swamp

83

crayfish [17].

84

The anti-viral responses of fish are similar to their mammalian counterparts [18,

85

19]. For example, the grouper STAT1a protein is a key element in response to

86

Iridovirus and Nodavirus infections [20]. The zebrafish fibroblast growth factor

87

receptor (FGFR3) is a negative regulator of innate immune responses by inhibiting

88

TANK-binding kinase 1 (TBK1) during viral infections. In contrast, the black carp

89

TAB1 is an activator of antiviral innate immunity against grass carp reovirus and

90

spring viremia of carp virus infections [21]. The involvement of

miRNAs

in

91

pathogen immune evasion or host responses in fish has been established and miRNA

92

regulatory pathways are beginning to be uncovered [22].

93

In the present study, we used high-throughput RNA sequencing to identify

94

differentially expressed mRNA and miRNAs produced in response to TiLV infection

95

in the spleens of Tilapia. This was performed to s explore mechanisms of TiLV-related

96

immune activation in the host to reveal underlying mechanisms for the prevention of

97

TiLV infections.

98 99

2. Materials and methods

100

2.1 Ethics statement

101

All experimental procedures were approved by China Guangdong Province Science

102

and Technology Department (SYXK (Yue) 2014-0136).

103 104

2.2 RNA isolation and RNA-seq library construction

105

The TiLV-2017A strain used in this study was kindly gifted by Dr. Sven Bergmann

106

from Institute of Infectology, Friedrich-Loffler-Institut (Greifswald, Germany ). Nile

107

tilapia with body lengths of 15 ± 0.5 cm and average weights of 100 g, were obtained

108

from a fish farm in Guangzhou (Guangdong, China). The fish were randomly checked

109

by RT-PCR to verify they were not infected with TiLV as previously described [2].

110

Tilapia were infected TiLV-2017A by intraperitoneal injection with 104 TCID50/0.1

111

mL phosphate buffered saline (PBS) and control fish received PBS alone in the same

112

volume. Spleens were collected, minced and homogenized 24 and 120h post-infection

113

(hpi) from controls(NC) and infected (T) fish and Total RNA was extracted using

114

Trizol reagent according to the manufacturer’s instructions (Invitrogen, Carlsbad, CA,

115

USA). RNA integrity and purity were measured by UV spectroscopy using a

116

NanoDrop ND-1000 instrument (NanoDrop, Wilmington, DE, USA). The sample

117

names were as follows: R-NC-1, R-NC-2, R-NC-3, R-T-24h-1, R-T-24h-2, R-T-24h-3,

118

R-T-120h-1, R-T-120h-2 and R-T-120h-3.

119 120

The bulk of the ribosome RNA (rRNA) was removed from samples using a

121

Ribo-Zero kit (Epicentre, Madison, WI, USA) and mRNA was isolated using oligo dT

122

beads. Enriched mRNA was fragmented and reverse transcribed using random

123

priming using a commercial kit and purified with a QiaQuick PCR extraction kit

124

(Qiagen, Hilden, Germany), The cDNA fragments were end repaired and

125

poly-adenylated then

126

provided by the manufacturer (Illumina, San Diego, CA, USA). Sequencing of cDNA

127

was performed using an Illumina HiSeq 2500 platform at Gene Denovo

128

Biotechnology(Guangzhou, China).

ligated to Illumina sequencing adapters using the protocol

129

Low quality and rRNA mapped raw reads were removed using the alignment tool

130

Bowtie2 (https://www.nature.com/articles/nmeth.1923) and clean reads were mapped

131

to the reference genome (Oreochromis_niloticus l1.0) using TopHat2 V2.0.3.12

132

(http://ccb.jhu.edu/software/tophat/index.shtml) and transcripts were reconstructed

133

using Cufflinks (http://cole-trapnell-lab.github.io/cufflinks/) and TopHat2. The

134

reconstructed transcripts were aligned to the reference genome to identify novel genes

135

or splice variants and divided into twelve categories using by Cuffcompare

136

(http://cole-trapnell-lab.github.io/cufflinks/). The class code “u” indicates transcripts

137

either unknown or located in intergenic spacers were defined as novel genes.

138

Gene abundance was quantified using RSEM (http://deweylab.biostat.wisc.edu/

139

rsem/) and a set of reference transcript sequences were obtained and preprocessed.

140

The Bowtie alignment program was utilized to realign RNA-seq reads to the reference

141

transcripts and these results were used to estimate gene abundance that was then

142

normalized using the fragments per kilobase of transcript per million mapped reads

143

(FPKM) method based on the following formula: FPKM =

10 NL/10

144

EdgeR software (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html)

145

was used to identify differentially expressed genes that possessed fold changes in

146

expression ≥2. Expression patterns of differently expressed mRNAs were determined

147

using expression data for each sample normalized to 0, log2 (v1/v0) and log2 (v2/v0)

148

that was then clustered using short time-series expression miner (STEM) software

149

(http://www.cs.cmu.edu/~jernst/stem). The clustered profiles with p-values ≤0.05

150

were defined as significant profiles. The differently expressed genes with significant

151

profiles were subjected to Gene Ontology (GO) and KEGG pathway enrichment

152

analysis using the GO database (http://www.geneontology.org/). Gene numbers of

153

every term were calculated and significantly enriched GO terms were defined using a

154

hypergeometric test [23]. Pathway enrichment analysis of different expressed mRNAs

155

was performed using the KEGG database (http://www.kegg.jp/kegg/kegg1.html).

156 157

2.3 Small RNA library construction and miRNA identification

158

Small RNAs (18–30 nt) were identified by polyacrylamide gel electrophoresis of

159

RNA samples taken immediately following rRNA removal (see above) and

160

appropriate sized bands were excised from the gels and small RNAs were isolated

161

using Gel Recovery Kit (Tiangen，China). Illumina adapters were added by ligation

162

using the protocol provided by the manufacturer (Illumina) were reverse transcribed

163

to generate a cDNA library and then sequenced as described above.

164

Low-quality raw reads and clean tags that were identified as rRNA, scRNA,

165

snoRNA, snRNA and tRNA using the Rfam database 11.0 (http://rfam.xfam.org/)

166

were removed and enriched clean tags were aligned to the reference genome

167

Oreochromis_niloticus l1.0(see above) and reads that

168

repeat sequences were removed. The clean tags were then used as searchterms in the

169

miRBase database Release 21 (http://www.mirbase.org/) to identify known miRNAs.

170

Unannotated tags were defined as novel miRNA candidates through their genome

171

positions

172

(http://www.tbi.univie.ac.at/RNA).

and

hairpin

structures

using

mapped to exons, introns and

the

software

Mireap_v0.2

173

The expression profiles of miRNAs were calculated and normalized to transcripts

174

per million (TPM) as the following formula: TPM=Actual miRNA counts/Total

175

counts of clean tags × 106 [24]. In addition, the heatmaps of total miRNA including

176

exist miRNA, known miRNA and novel miRNA were drew to display miRNA

177

expression levels in different groups and to cluster miRNAs with similar expression

178

pattern.

179 180

Differentially expressed miRNAs across different groups were identified using the following formula:

181 182 183

miRNAs with a fold change ≥2 and P value <0.05 were considered as significant differently expressed miRNAs.

184 185

The miRNA targets were predicted using RNAhybrid (v2.1.2) +svm_light (v6.01),

186

Miranda (v3.3a) and TargetScan v7.0 together (http://www.targetscan.org/vert_71/)

187

and the miRNA-target network was constructed using Cytoscape v3.6.0

188

(http://manual.cytoscape.org/en/stable/index.html).

189 190 191

All target genes of differentially expressed miRNAs were enriched to GO terms using the Gene Ontology database and KEGG analysis (see above).

192 193

Selected miRNA targets were validated by qRT-PCR analysis of cFDNA

194

generated from mRNA samples using random priming and were quantified using a

195

SYBR Green qPCR SuperMix (Bio-Rad, Hercules, CA, USA) and PCR primers as

196

indicated(Table 1) on an ABI Real Time PCR system (Applied Biosystems, Watham,

197

MA, USA).

198

Expression levels of each gene was normalized to the level of GAPDH expression and

199

groups were compared using the 2−∆∆CT method [25, 26]. P values < 0.05 were

200

considered statistically significant.

Ct values were calculated using software provided with the instruments.

201

The miRNA-Seq data was validated by first generating cDNA copies of miRNAs

202

using the miRcute miRNA First-Strand cDNA Synthesis kit (Tiangen, Beijing, China)

203

and real-time PCR as described above (see Table 1 for PCR primers). Levels of

204

miRNA were normalized to the levels of the splicesomal RNA U6 and compared

205

using the 2−∆∆CT method as described above.

206 207

3. Results

208

3.1 RNA-seq analysis and validation

209

We generated 9 cDNA libraries for virus infected (T) and control (NC) fish

210

including 3 libraries for PBS controls taken at 120 hpi and 6 for virus-infected sample

211

groups taken at 24 and 120 hpi. We obtained an average of 55,828,800 clean reads

212

from each group at a level of 98% (Table S1). The clean reads were mapped to the

213

Tilapia reference genome and >82% could be mapped to the genome (Table S2).

214

Furthermore, we identified 25,359 (87.67%) known genes and 863 (12.33%) novel

215

genes. Interestingly, the quantities of identified genes were greater for the

216

virus-infected groups and these most likely were comprised of TiLV genes (Table S3).

217

However, gene abundance for the T-120 h group was less than for the other 6 groups

218

Fig. S1).

219 220

We used the Pearson correlation coefficient (PCC) to assess the reproducibility of

221

the RNA-Seq results according to the expression levels for each group. The

222

correlations within the 3 groups of 3 replicates ranged from 97.2 to 99.9% indicating

223

only minor differences between replicates (Fig. S2). We identified 1770 differently

224

expressed genes (860 up-regulated and 910 down-regulated) in the T-24 group and

225

4142 genes exhibited different in the R-T-120 group (1757 up-regulated and 2385

226

down-regulated) compared with the NC control group (Fig. 1 and Fig. S3, Tables S4

227

and S5). When we compared the virus-infected groups, 4947 genes were significantly

228

altered (Fig. 1 and Fig. S3, Table S6). These results indicated that different regulated

229

genes responded to TiLV infection during different stages. We sampled 20 of these

230

differentially expressed genes (Table 2) and validated the RNA-Seq results using

231

qRT-PCR. We found that expression levels for these gene were mostly consistent with

232

the RNA-Seq results (Fig.2). We also identified genes that were differentially

233

expressed between 24 and 120 hpi such as GeneID_100703853 (Table S6).

234 235

3.3 Trend analysis for RNA-Seq

236

The differentially expressed genes we found were different when the infected groups

237

were compared one by one to the controls. The differently expressed genes were

238

clustered into eight profiles and most were enriched in profiles 3 and 4 (Fig. 3A and

239

3B, Tables S7 and S8). For example, in profile 3, genes that were upregulated from 0

240

-24 hpi were downregulated at 24-120 hpi (Fig. 3C). In contrast, expression levels of

241

genes that decreased during 0-24 hpi were increased from 24-120 hpi (Fig. 3D).In

242

these groups of differentially expressed genes, we found 291 GO terms that were

243

enriched. These included genes for immune system progress, behavior, biological

244

adhesion, biological regulation, virion, virion part, cell junction, antioxidant activity,

245

binding and catalytic activity (Fig. 4A-C, Table S9). Interestingly, virion and virion

246

part were not involved in GO terms of differently expressed genes between the NC

247

control group and the T-24 h infected group (Fig. 4A). KEGG pathway enrichment

248

analysis revealed that these

249

PPAR signaling, phagosome, starch and sucrose metabolism, drug metabolism,

250

metabolism of xenobiotics by cytochrome P450, Salmonella infection and others (Fig.

251

5A-C, Table S9).

genes were associated with 62 pathways that included

252 253

3.5 miRNA-Seq analysis

254

In our 3 control and 6 experimental groups we identified 13029865 (92.13%),

255

12538729 (92.16%) and 12735309 (91.16%) clean tags for the control (NC) group,

256

14523784 (87.5412%), 13722343 (92.43%) and 15394438 (93.34%) for the infected

257

group T-24h and 13865591 (93.10%), 15165378 (92.50%) and 15445062 (88.36%)

258

for the infected T-120h group (Table S10). Furthermore, a length analysis of clean

259

tags indicated that 92.78 -96.59% of them were distributed between 19 and 24 nt,

260

consistent with the typical length of miRNAs (Table S11).

261 262

3.6 miRNA identification

263

The clean miRNA populations were further sorted by eliminating any rRNA,

264

scRNA, snoRNA, snRNA and tRNA sequences. This resulted in the identification of

265

the following conserved miRNAs for each of the groups: NC (298, 305 and 300), T-24

266

h (302, 310 and 317) and T-120 h (303, 308 and 316) (Table S12). We also identified

267

that these 3 miRNA population groups were contained known miRNA orthologs from

268

other animal and plant species as follows: NC (412, 417 and 374), T-24 h (387, 525

269

and 443) and T-120 h (374, 400 and 437) (Table S13). We also identified novel

270

miRNAs as follows: NC (472, 421 and 409), T-24 h (519, 554 and 555) and T-120h

271

(538, 592 and 482) (Table S14). Interestingly, none of the novel miRNAs were

272

derived from the TiLV genome. The hairpin structures of four novel miRNAs are

273

displayed in Fig. 6A-D.

274

When we examined the populations of differentially expressed miRNAs between

275

virus infected and controls, we found in 138 (91 up-regulated and 47 down-regulated)

276

in T-24 h and 198 in T-120 h (126 up-regulated and 72 down-regulated) (Fig.7A,

277

Tables S15 and 16). A comparision between T-24 h and T-120 h groups, we found 153

278

miRNA levels that were dramatically altered (Fig. 7B, Table S17). These data

279

indicated that there are groups of differentially regulated miRNAs that are expressed

280

during the different stages of TiLV infection. We then selected 19 differently

281

expressed miRNAs for validation by qRT-PCR. These results were consistent with

282

those levels generated using miRNA-Seq (Fig .8A and 8B). The miRNAs exhibited

283

differential expression patterns at different infection stages such as oni-miR-7133 that

284

was downregulated at 24 hpi but upregulated at 120 hpi (Fig. 8A and 8B).

285 286

We performed a trend analysis to identify expression patterns of these

287

differentially expressed miRNAs. Consistent with the results of RNA-Seq, the

288

differently expressed miRNAs were clustered into eight profiles and most of the genes

289

were enriched in profiles 3 and 6 (Fig. 9A and 9B, Tables S18 and S19). In profile 3,

290

expression levels of miRNAs were upregulated over 0-24 hpi while being

291

downregulated from 24-120 hpi. In contrast, expression levels of miRNAs decreased

292

during 0 -24 hpi and increased from 24-120 hpi (Fig. 9C).

293 294 295

3.9 Target prediction of differently expressed miRNAs and functional analysis The biological roles for miRNAs are achieved by regulating target gene

296

expression. We identified immune target genes during TiLV infection that included

297

erythropoietin isoform X2 (epo), double-stranded RNA-specific adenosine deaminase

298

isoform X1 (adar), bone morphogenetic protein 4 (bmp4) and tapasin-related protein

299

(tapbpl). More importantly, these targets overlapped with differentially expressed

300

mRNAs (Table S4 and S5). To further identify related physiological processes and

301

pathways, we utilized GO and KEGG pathway enrichment. Our results paralled that

302

of the sets of differently expressed mRNAs where the GO terms included immune

303

system progress, behavior, biological adhesion, biological regulation, virion, virion

304

part, cell junction, antioxidant activity, binding and catalytic activity were all enriched

305

(Fig. 10A-C, Table S20). In addition, KEGG pathway enrichment analysis revealed

306

targets associated with melanogenesis, ErbB, Wnt, VEGF and MAPK signaling

307

pathways along with lysosome, apoptosis and endocytosis (Fig. 11 A-C, Table S21).

308 309

3.10 Integrated analysis of mRNA and miRNA

310

Our results indicated that the differently expressed mRNA and miRNA populations

311

shared similar expression trends. In addition, the GO terms of differently expressed

312

mRNA were consistent with those for targets of differently expressed miRNAs.

313

Therefore, we performed an integrative analysis to see if predicted miRNA targets

314

were identified as differently expressed mRNAs. We found that in the groups NC,

315

T-24 h and T-120 h possessed 3951, 12,480 and 13,784 mRNA-miRNA pairs that

316

were altered in opposite directions (Supplementary Fig. S3-5, Table S22). This

317

indicated that mRNA expression was negatively correlated to miRNA expression in

318

response to TiLV infection. Twenty-one of these mRNA-miRNA pairs were listed in

319

Table 3, in which miRNAs had been validated by qPCR. We were able to construct a

320

miRNA-mRNA network indicating that most of these miRNAs play roles in TiLV

321

infections (Fig. 12). An integrated GO and KEGG pathway analysis revealed that 208

322

GO terms were enriched and included immune system progress, behavior, biological

323

adhesion, biological regulation, virion, virion part, cell junction, antioxidant activity,

324

binding and catalytic activity (Fig. 13A-C, Table S23). In addition, these mRNAs

325

were grouped into 58 pathways that included PPAR signaling, cell adhesion

326

molecules, phagosome, intestinal immune network for IgA production, starch and

327

sucrose metabolism, drug metabolism, ABC transporters and Salmonella infection

328

(Fig. 14A-C, Table S24).

329 330

4. Discussion

331

In the present study, we identified the expression profile of mRNAs and miRNAs

332

from normal and TiLV-infected tilapia. We found 26,122 differentially expressed

333

genes containing 863 new genes and 1433 miRNAs including 592 novel miRNAs.

334

This is the largest number of miRNAs identified from a single fish study to date. For

335

example, 381 host miRNAs and 9 viral miRNAs were obtained from Japanese

336

flounder (Paralichthys olivaceus) infected with megalocytivirus [27] and only 29

337

miRNAs were identified at the metamorphic stage of Japanese flounder [28]. A recent

338

study identified 289 miRNAs containing 44 novel miRNAs were present in zebrafish

339

vitellogenic ovarian follicular cells [29]. Moreover, 381 known miRNAs and 926

340

novel miRNAs related to the immune response were demonstrated in rainbow trout

341

upon

342

present study identified 863 new genes and 592 novel miRNAs supplements the the

343

tilapia genome [31, 32] and the known miRNA pool in fish [33].

infection withv

Aeromonas salmonicida

subsp. salmonicida [30]. Our

344

We performed this study to discover mechanisms used by TiLV to cause infections.

345

In general, miRNAs have been shown to play important roles in fish upon virus

346

infection by regulating target genes associated with the immune response, apoptosis

347

and virus replication [16, 34-38]. To our knowledge, this is the first study that

348

characterizes the expression profiles of mRNAs and miRNAs in tilapia upon TiLV

349

infection. We identified about 5000 differently expressed mRNAs and most were

350

associated with the immune response including a predicted interleukin-15-like

351

isoform X1 (102081732), four class II histocompatibility antigens (109194343,

352

100698680, 100703853 and 100702693) and two polymeric Ig receptor-like genes

353

(102080724 and 100693666) [39-41] (Table S9). In addition, we also found about 200

354

differently expressed miRNAs regulating genes involved in the immune response

355

through

KEGG

pathway

analysis,

including

102081732,

356

100703853,100703681,100702693, 100698680, 109194343, 102080724, 100693666,

357

epo, adar bmp4 and tapbp. These results indicated that Tilapia mount a robust immune

358

response to TiLV infection and miRNAs are key players that modulate the response to

359

viral infection by regulating immume-releated mRNAs.

360

Further functional enrichment analysis of differently expressed mRNAs and

361

miRNAs demonstrated that processes including immune system progress, behavior,

362

biological adhesion, biological regulation, virion, virion part, cell junction,

363

antioxidant activity, binding, catalytic activity are closed related to TiLV infection.

364

The immune system plays a critical role in antagonizing virus infections in fish [42,

365

43]. In addition, miRNAs play important roles in fish upon virus infection through

366

regulating target genes associated with immune response. For instance, miR-21

367

modulates the inflammatory response through regulating jnk and ccr7 expression in

368

grass carp after infection [34]. MiR-7a suppresses host immune response via

369

PI3K/AKT/GSK3β

370

Edwardsiella tarda infection [35]. Moreover, pathogens infect fish through adhesion

371

to skin mucus and cells and the virion facilitates infection [44, 45]. Antioxidant

372

activity is also involved in modulating immune responsiveness in fish [46, 47].

373

Catalytic activity of protein kinase is also required for the induction of immune

374

response upon virus infection [48]. Thus, these catagories of mRNA and miRNAs

375

converge to regulate the response to TiLV infection. Interestingly, virion and virion

376

part were not involved in GO terms of differently expressed mRNAs between the NC

377

and T-24 h groups suggesting that TiLV may replicate after 24 hpi.

signaling

pathway

in

Paralichthys

olivaceus

following

378

Our KEGG pathway enrichment analysis revealed that differently expressed

379

mRNAs were associated with PPAR signaling, phagosome, starch and sucrose

380

metabolism, drug and P-450 metabolism, Salmonella infection, while targets of

381

differently expressed miRNAs were most enriched for pathways of melanogenesis,

382

ErbB, Wnt, VEGFand MAPK signaling along with , lysosome, apoptosis and

383

endocytosis regulation. Recent studies have indicated that PPAR signaling pathway

384

exerts antiviral activity against fish viral infections [49, 50]. The phagosome pathway

385

contributes to Flavobacterium columnare infection in the top mouth culter (Culter

386

alburnus) [51]. In addition, cytochrome P450 is involved in disease defense responses

387

of channel catfish [52]. Furthermore, miRNAs play important roles in fish upon virus

388

infection by regulating target genes associated with apoptosis [36]. However, no

389

studies have revealed the involvement of the other pathways in fish infection. Notably,

390

ErbB signaling, Salmonella infection, Wnt, VEGF and MAPK signaling and lysosome

391

play important roles in infection [53-58]. This indicates that TiLV-related mRNAs and

392

miRNAs may respond to infection through the above pathways.

393

Interactions between miRNA and target genes are critical for the balance of gene

394

expression in fish following infection [38,59,60]. In this study, about 14,000

395

mRNA-miRNA pairs exhibited opposite expression profiles including the verified

396

miRNAs(oni-miR-192-x,oni-miR-215-x,oni-miR-10856, oni-miR-192, oni-miR-7133,

397

novel-m0060-3p, oni-miR-235-y, oni-miR-92b, oni-miR-541-x, oni-miR-7133-y,

398

novel-m0344-5p, oni-miR-10556-x) and their target genes (109204230, 100702234,

399

109194343, 109202867, 109194454, 100703853, 100703681, 100702693, 100701693,

400

109204294 and100701969). Therefore, miRNAs may modulate host immune

401

response by altering target gene expression following TiLV infection.

402

In summary, we identified a series of TiLV-related mRNA and miRNAs using

403

high-throughput mRNA and miRNA sequencing. We identified 26,122 expressed

404

genes including 863 new genes and 1433 miRNAs containing 592 novel miRNAs. We

405

found ~5000 mRNAs and 200 miRNAs that were differentially expressed after TiLV

406

infection. These were classified into the immune system category and involved high

407

level enrichment for in PPAR, phagosome and starch and sucrose metabolism

408

pathways. The targets of differently expressed miRNAs were classified into the

409

immune system category and were enriched in the melanogenesis, ErbB and Wnt

410

pathways.

411 412

Acknowledgements

413

This work was supported by the Natural Science Foundation of Guangdong Province

414

(NO. 2018A030313757), Central Public ‐ interest Scientific Institution Basal

415

Research Fund, CAFS (NO. 2018SJ ‐YB03), Central Public interest Scientific

416

Institution Basal Research Fund CAFS (NO. 2019ZD0703), the China Agriculture

417

Research System (NO. CARS-45), Central Public-interest Scientific Institution Basal

418

Research Fund, CAFS (NO. 2019XN-001).

419 420 421

Disclosure Statements

422

The authors declare no conflict of interest.

423 424

References

425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452

[1] P. Tattiyapong, W. Dachavichitlead, W. Surachetpong, Experimental infection of Tilapia Lake Virus (TiLV) in Nile tilapia (Oreochromis niloticus) and red tilapia (Oreochromis spp.), Veterinary microbiology 207 (2017) 170-177. [2] J.E. Kembou Tsofack, R. Zamostiano, S. Watted, A. Berkowitz, E. Rosenbluth, N. Mishra, T. Briese, W.I. Lipkin, R.M. Kabuusu, H. Ferguson, J. Del Pozo, A. Eldar, E. Bacharach, Detection of Tilapia Lake Virus in Clinical Samples by Culturing and Nested Reverse Transcription-PCR, Journal of clinical microbiology 55(3) (2017) 759-767. [3] R. Wang, W.X. Wang, Diet-specific trophic transfer of mercury in tilapia (Oreochromis niloticus): Biodynamic perspective, Environ Pollut 234 (2018) 288-296. [4] M. Shlapobersky, M.S. Sinyakov, M. Katzenellenbogen, R. Sarid, J. Don, R.R. Avtalion, Viral encephalitis of tilapia larvae: primary characterization of a novel herpes-like virus, Virology 399(2) (2010) 239-47. [5] L. Bigarre, J. Cabon, M. Baud, M. Heimann, A. Body, F. Lieffrig, J. Castric, Outbreak of betanodavirus infection in tilapia, Oreochromis niloticus (L.), in fresh water, J Fish Dis 32(8) (2009) 667-73. [6] E. Bacharach, N. Mishra, T. Briese, M.C. Zody, J.E. Kembou Tsofack, R. Zamostiano, A. Berkowitz, J. Ng, A. Nitido, A. Corvelo, N.C. Toussaint, S.C. Abel Nielsen, M. Hornig, J. Del Pozo, T. Bloom, H. Ferguson, A. Eldar, W.I. Lipkin, Characterization of a Novel Orthomyxo-like Virus Causing Mass Die-Offs of Tilapia, MBio 7(2) (2016) e00431-16. [7] M. Eyngor, R. Zamostiano, J.E. Kembou Tsofack, A. Berkowitz, H. Bercovier, S. Tinman, M. Lev, A. Hurvitz, M. Galeotti, E. Bacharach, A. Eldar, Identification of a novel RNA virus lethal to tilapia, Journal of clinical microbiology 52(12) (2014) 4137-46. [8] W. Surachetpong, T. Janetanakit, N. Nonthabenjawan, P. Tattiyapong, K. Sirikanchana, A. Amonsin, Outbreaks of Tilapia Lake Virus Infection, Thailand, 2015-2016, Emerg Infect Dis 23(6) (2017) 1031-1033. [9] R.C. Lee, R.L. Feinbaum, V. Ambros, The C. elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14, Cell 75(5) (1993) 843-54. [10] B. Wightman, I. Ha, G. Ruvkun, Posttranscriptional regulation of the heterochronic gene lin-14 by lin-4 mediates temporal pattern formation in C. elegans, Cell 75(5) (1993) 855-62.

453 454 455 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

[11] X. Cai, A. Schafer, S. Lu, J.P. Bilello, R.C. Desrosiers, R. Edwards, N. Raab-Traub, B.R. Cullen, Epstein-Barr virus microRNAs are evolutionarily conserved and differentially expressed, PLoS pathogens 2(3) (2006) e23. [12] Y. Tutar, miRNA and cancer; computational and experimental approaches, Curr Pharm Biotechnol 15(5) (2014) 429. [13] C. Backes, E. Meese, A. Keller, Specific miRNA Disease Biomarkers in Blood, Serum and Plasma: Challenges and Prospects, Mol Diagn Ther 20(6) (2016) 509-518. [14] I.W. Boss, R. Renne, Viral miRNAs and immune evasion, Biochimica et biophysica acta 1809(11-12) (2011) 708-14. [15] W. Chen, L. Yi, S. Feng, L. Zhao, J. Li, M. Zhou, R. Liang, N. Gu, Z. Wu, J. Tu, L. Lin, Characterization of microRNAs in orange-spotted grouper (Epinephelus coioides) fin cells upon red-spotted grouper nervous necrosis virus infection, Fish & shellfish immunology 63 (2017) 228-236. [16] C. Zhang, N. Li, X. Fu, Q. Lin, L. Lin, J. Tu, MiR-214 inhibits snakehead vesiculovirus (SHVV) replication by targeting host GS, Fish & shellfish immunology 84 (2019) 299-303. [17] H. Yang, X. Li, J. Ji, C. Yuan, X. Gao, Y. Zhang, C. Lu, F. Li, X. Zhang, Changes of microRNAs expression profiles from red swamp crayfish (Procambarus clarkia) hemolymph exosomes in response to WSSV infection, Fish & shellfish immunology 84 (2019) 169-177. [18] K. Murata, J.H. Kang, S. Nagashima, T. Matsui, Y. Karino, Y. Yamamoto, T. Atarashi, M. Oohara, M. Uebayashi, H. Sakata, K. Matsubayashi, K. Takahashi, M. Arai, S. Mishiro, M. Sugiyama, M. Mizokami, H. Okamoto, IFN-lambda3 as a host immune response in acute hepatitis E virus infection, Cytokine 125 (2019) 154816. [19] N. Johnson, A.F. Cunningham, A.R. Fooks, The immune response to rabies virus infection and vaccination, Vaccine 28(23) (2010) 3896-901. [20] J. Zhang, X. Huang, S. Ni, J. Liu, Y. Hu, Y. Yang, Y. Yu, L. Zhou, Q. Qin, Y. Huang, Grouper STAT1a is involved in antiviral immune response against iridovirus and nodavirus infection, Fish & shellfish immunology 70 (2017) 351-360. [21] Z. Zou, X. Xie, W. Li, X. Song, Y. Tan, H. Wu, J. Xiao, H. Feng, Black carp TAB1 up-regulates TAK1/IRF7/IFN signaling during the antiviral innate immune activation, Fish & shellfish immunology 89 (2019) 736-744. [22] M.D. Mehta, P.T. Liu, microRNAs in mycobacterial disease: friend or foe?, Front Genet 5 (2014) 231. [23] J. Cao, S. Zhang, A Bayesian extension of the hypergeometric test for functional enrichment analysis, Biometrics 70(1) (2014) 84-94. [24] A. Mortazavi, B.A. Williams, K. McCue, L. Schaeffer, B. Wold, Mapping and quantifying mammalian transcriptomes by RNA-Seq, Nature methods 5(7) (2008) 621-8. [25] S. Sadeghi, Z. Hojati, H. Tabatabaeian, Cooverexpression of EpCAM and c-myc genes in malignant breast tumours, Journal of genetics 96(1) (2017) 109-118. [26] A. Arocho, B. Chen, M. Ladanyi, Q. Pan, Validation of the 2-DeltaDeltaCt calculation as an alternate method of data analysis for quantitative PCR of BCR-ABL P210 transcripts, Diagnostic molecular pathology : the American journal of surgical pathology, part B 15(1) (2006) 56-61. [27] B.C. Zhang, J. Zhang, L. Sun, In-depth profiling and analysis of host and viral microRNAs in Japanese flounder (Paralichthys olivaceus) infected with megalocytivirus reveal involvement of microRNAs in host-virus interaction in teleost fish, BMC genomics 15 (2014) 878. [28] C. Xie, S. Xu, L. Yang, Z. Ke, J. Xing, J. Gai, X. Gong, L. Xu, B. Bao, mRNA/microRNA Profile at the

497 498 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

Metamorphic Stage of Olive Flounder (Paralichthys olivaceus), Comparative and functional genomics 2011 (2011) 256038. [29] Y. Zayed, X. Qi, C. Peng, Identification of Novel MicroRNAs and Characterization of MicroRNA Expression Profiles in Zebrafish Ovarian Follicular Cells, Frontiers in endocrinology 10 (2019) 518. [30] Y. Cao, D. Wang, S. Li, J. Zhao, L. Xu, H. Liu, T. Lu, Z. Mou, A transcriptome analysis focusing on splenic immune-related mciroRNAs of rainbow trout upon Aeromonas salmonicida subsp. salmonicida infection, Fish & shellfish immunology 91 (2019) 350-357. [31] M.A. Conte, W.J. Gammerdinger, K.L. Bartie, D.J. Penman, T.D. Kocher, A high quality assembly of the Nile Tilapia (Oreochromis niloticus) genome reveals the structure of two sex determination regions, BMC genomics 18(1) (2017) 341. [32] J.H. Xia, Z. Bai, Z. Meng, Y. Zhang, L. Wang, F. Liu, W. Jing, Z.Y. Wan, J. Li, H. Lin, G.H. Yue, Signatures of selection in tilapia revealed by whole genome resequencing, Scientific reports 5 (2015) 14168. [33] T.T. Bizuayehu, I. Babiak, MicroRNA in teleost fish, Genome biology and evolution 6(8) (2014) 1911-37. [34] L. Tao, X. Xu, Y. Fang, A. Wang, F. Zhou, Y. Shen, J. Li, miR-21 targets jnk and ccr7 to modulate the inflammatory response of grass carp following bacterial infection, Fish & shellfish immunology 94 (2019) 258-263. [35] Y. Liu, Y. Liu, M. Han, X. Du, X. Liu, Q. Zhang, J. Liu, Edwardsiella tarda-induced miR-7a functions as a suppressor in PI3K/AKT/GSK3beta signaling pathway by targeting insulin receptor substrate-2 (IRS2a and IRS2b) in Paralichthys olivaceus, Fish & shellfish immunology 89 (2019) 477-485. [36] S. Ni, Y. Yan, H. Cui, Y. Yu, Y. Huang, Q. Qin, Fish miR-146a promotes Singapore grouper iridovirus infection by regulating cell apoptosis and NF-kappaB activation, The Journal of general virology 98(6) (2017) 1489-1499. [37] S. Ni, Y. Yu, J. Wei, L. Zhou, S. Wei, Y. Yan, X. Huang, Y. Huang, Q. Qin, MicroRNA-146a promotes red spotted grouper nervous necrosis virus (RGNNV) replication by targeting TRAF6 in orange spotted grouper, Epinephelus coioides, Fish & shellfish immunology 72 (2018) 9-13. [38] R. Andreassen, B. Hoyheim, miRNAs associated with immune response in teleost fish, Developmental and comparative immunology 75 (2017) 77-85. [39] S. Liu, Y. Du, X. Sheng, X. Tang, J. Xing, W. Zhan, Molecular cloning of polymeric immunoglobulin receptor-like (pIgRL) in flounder (Paralichthys olivaceus) and its expression in response to immunization with inactivated Vibrio anguillarum, Fish & shellfish immunology 87 (2019) 524-533. [40] D.J. Wcisel, J.A. Yoder, The confounding complexity of innate immune receptors within and between teleost species, Fish & shellfish immunology 53 (2016) 24-34. [41] G.J. Niu, S. Wang, J.D. Xu, M.C. Yang, J.J. Sun, Z.H. He, X.F. Zhao, J.X. Wang, The polymeric immunoglobulin receptor-like protein from Marsupenaeus japonicus is a receptor for white spot syndrome virus infection, PLoS pathogens 15(2) (2019) e1007558. [42] C.J. Secombes, L.H. Chappell, Fish immune responses to experimental and natural infection with helminth parasites, Annual Review of Fish Diseases 6(96) (1996) 167-177. [43] M.D. Fast, Fish immune responses to parasitic copepod (namely sea lice) infection, Developmental and comparative immunology 43(2) (2014) 300-12. [44] I. Tsilinskii, E.F. Karpova, [Variability of attenuated alphavirus strains in virion aggregate infection], Voprosy virusologii 28(5) (1983) 601-7. [45] D.M. Owen, H. Huang, J. Ye, M. Gale, Jr., Apolipoprotein E on hepatitis C virion facilitates infection

541 542 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

through interaction with low-density lipoprotein receptor, Virology 394(1) (2009) 99-108. [46] W.D. Jiang, K. Hu, Y. Liu, J. Jiang, P. Wu, J. Zhao, Y.A. Zhang, X.Q. Zhou, L. Feng, Dietary myo-inositol modulates immunity through antioxidant activity and the Nrf2 and E2F4/cyclin signalling factors in the head kidney and spleen following infection of juvenile fish with Aeromonas hydrophila, Fish & shellfish immunology 49 (2016) 374-86. [47] M. Reyes-Becerril, V. Sanchez, K. Delgado, K. Guerra, E. Velazquez, F. Ascencio, C. Angulo, Caspase -1, -3, -8 and antioxidant enzyme genes are key molecular effectors following Vibrio parahaemolyticus and Aeromonas veronii infection in fish leukocytes, Immunobiology 223(10) (2018) 562-576. [48] N. Taghavi, C.E. Samuel, Protein kinase PKR catalytic activity is required for the PKR-dependent activation of mitogen-activated protein kinases and amplification of interferon beta induction following virus infection, Virology 427(2) (2012) 208-16. [49] Y. Wang, Y. Yu, Q. Wang, S. Wei, S. Wang, Q. Qin, M. Yang, PPAR-delta of orange-spotted grouper exerts antiviral activity against fish virus and regulates interferon signaling and inflammatory factors, Fish & shellfish immunology 94 (2019) 38-49. [50] S. Mohapatra, T. Chakraborty, M.A. Reza, S. Shimizu, T. Matsubara, K. Ohta, Short-term starvation and realimentation helps stave off Edwardsiella tarda infection in red sea bream (Pagrus major), Comparative biochemistry and physiology. Part B, Biochemistry & molecular biology 206 (2017) 42-53. [51] L. Zhao, J. Tu, Y. Zhang, J. Wang, L. Yang, W. Wang, Z. Wu, Q. Meng, L. Lin, Transcriptomic analysis of the head kidney of Topmouth culter (Culter alburnus) infected with Flavobacterium columnare with an emphasis on phagosome pathway, Fish & shellfish immunology 57 (2016) 413-418. [52] J. Zhang, J. Yao, R. Wang, Y. Zhang, S. Liu, L. Sun, Y. Jiang, J. Feng, N. Liu, D. Nelson, G. Waldbieser, Z. Liu, The cytochrome P450 genes of channel catfish: their involvement in disease defense responses as revealed by meta-analysis of RNA-Seq data sets, Biochimica et biophysica acta 1840(9) (2014) 2813-28. [53] J. Ho, D.L. Moyes, M. Tavassoli, J.R. Naglik, The Role of ErbB Receptors in Infection, Trends in microbiology 25(11) (2017) 942-952. [54] O.H. Pham, S.J. McSorley, Protective host immune responses to Salmonella infection, Future microbiology 10(1) (2015) 101-10. [55] J.L. Smith, S. Jeng, S.K. McWeeney, A.J. Hirsch, A MicroRNA Screen Identifies the Wnt Signaling Pathway as a Regulator of the Interferon Response during Flavivirus Infection, Journal of virology 91(8) (2017). [56] T. Weinkopff, C. Konradt, D.A. Christian, D.E. Discher, C.A. Hunter, P. Scott, Leishmania major Infection-Induced VEGF-A/VEGFR-2 Signaling Promotes Lymphangiogenesis That Controls Disease, Journal of immunology 197(5) (2016) 1823-31. [57] F. Qi, S. Bai, D. Wang, L. Xu, H. Hu, S. Zeng, R. Chai, B. Liu, Macrophages produce IL-33 by activating MAPK signaling pathway during RSV infection, Molecular immunology 87 (2017) 284-292. [58] Y. Miao, G. Li, X. Zhang, H. Xu, S.N. Abraham, A TRP Channel Senses Lysosome Neutralization by Pathogens to Trigger Their Expulsion, Cell 161(6) (2015) 1306-19. [59] X. Liu, J. Tu, J. Yuan, X. Liu, L. Zhao, F.U. Dawar, M.N. Khattak, A.M. Hegazy, N. Chen, V.N. Vakharia, L. Lin, Identification and Characterization of MicroRNAs in Snakehead Fish Cell Line upon Snakehead Fish Vesiculovirus Infection, International journal of molecular sciences 17(2) (2016). [60] N. Wang, R. Wang, R. Wang, Y. Tian, C. Shao, X. Jia, S. Chen, The integrated analysis of RNA-seq and microRNA-seq depicts miRNA-mRNA networks involved in Japanese flounder (Paralichthys olivaceus) albinism, PloS one 12(8) (2017) e0181761.

585 586

587

Figure legends

588

Figure 1. Characteristics of mRNA expression levels between different groups. (A)

589

NC, non-infected control (B) infected T-24h and (C) infected T-120h. Differentially

590

expressed mRNAs are shown in red (up-regulated) or green (down-regulated). (D)

591

Numbers of differentially expressed mRNAs.

592 593

Figure 2. Validation of RNA-Seq data . Expression levels of 20 selected mRNA as

594

determined by qRT-PCR.

595 596

Figure 3. Trend analysis for differently expressed mRNAs. (A) Profiles of 8

597

differently expressed mRNAs. (B) Numbers of differentially expressed mRNAs in

598

each profile. (C) Details of profile 3 in which mRNAs were upregulated at 0-24

599

hpiand downregulated at24 -120 hpi. (D) Details of profile 4 in which mRNAs were

600

decreased at 0-24 hpi and increased at24-120 hpi.

601 602

Figure 4. GO enrichment analysis for differently expressed mRNAs. GO

603

enrichment histograms and GO terms for differently expressed mRNAs in (A) NC,

604

non-infected control (B) infected T-24h and (C) infected T-120h.

605 606

Figure 5. KEGG pathway enrichment analysis for differently expressed mRNAs.

607

The KEGG pathway enrichment scatter plots for differently expressed mRNAs in (A)

608

NC, non-infected control (B) infected T-24h and (C) infected T-120h.

609 610

Figure 6. Hairpin structures of four novel miRNAs. The secondary structures of

611

four novel miRNA identified in this study(A) novel-m0001 (B) novel-m0002 (C)

612

novel-m0005 and (D) novel-m0006.

613 614

Figure 7. Characteristics of miRNA expression levels between different groups.

615

(A) Numbers of differentially expressed miRNAs. (B) Expressed miRNAs (P < 0.05)

616

among the 3 different groups. Colors from green to red stand for z-scores determined

617

from a dimensionality reduction of FPKM value that revealed decreasing miRNA

618

levels for each group.

619 620

Figure 8. Validation of miRNAs by qRT-PCR. (A) Expression of 19 selected

621

miRNAs (P < 0.05) .

622

Colors from green to red stand for z-scores determined from a dimensionality

623

reduction of FPKM value that revealed decreasing miRNA levels for each group.

624

(B) 19 selected miRNAs expression levels determined by qRT-PCR.

625 626

Figure 9. Trend analysis for differently expressed miRNAs. (A) Profiles of 8

627

differently expressed mRNAs (B) Numbers of differentially expressed mRNAs in

628

each profile. (C) Details of profile 3 in which mRNAs were upregulated at 0-24

629

hpiand downregulated at24 -120 hpi. (D) Details of profile 6 in which mRNAs were

630

increased at 0-24 hpi.

631

Figure 10. GO enrichment analysis for target genes regulated by differently

632

expressed miRNAs. The GO enrichment histograms and GO terms for target genes

633

regulated by differently expressed miRNAs in groups (A) NC, non-infected control (B)

634

infected T-24h and (C) infected T-120h.

635 636

Figure 11. KEGG pathway enrichment analysis for target genes regulated by

637

differently expressed miRNAs. KEGG pathway enrichment scatter plots for target

638

genes regulated by differently expressed miRNAs in groups (A) NC, non-infected

639

control (B) infected T-24h and (C) infected T-120h are shown.

640 641

Figure 12. miRNA-mRNA regulatory network between validated miRNAs and

642

negatively regulated target genes. View of miRNA-mRNA regulatory network

643

according to 19 validated miRNAs and their negatively regulated target genes.

644 645

Figure 13. GO enrichment analysis for negatively regulated target genes of

646

differently expressed miRNAs. (A-C) The GO enrichment histograms and GO terms

647

for negatively regulated target genes of differently expressed miRNAs in different

648

groups are shown.

649 650

Figure 14. KEGG pathway enrichment analysis for negatively regulated target

651

genes of differently expressed miRNAs. The KEGG pathway enrichment scatter

652

plots for negatively regulated target genes of differently expressed miRNAs in

653

different groups are shown.

654

Table 1 Primers list Name

Sequence(5`-3`)

GeneID_109194299 F

GATGCCGCTGATTGGTGTTC

GeneID_109194299 R

GCAGAAGCGACCTTTTTGGG

GeneID_109194454 F

CTGTGGAAAACATGGCTGCG

GeneID_109194454 R

ACGCTGCACACGAGTATTGT

GeneID_100701969 F

CAGTCACCATCCTGCCATGT

GeneID_100701969 R

ATGGCGAGCTTGTTCCTCTC

GeneID_102081732 F

GCTCCCGGTGATCCTCTTTC

GeneID_102081732 R

GGCACATGTAGGCGTACTCA

GeneID_100702234 F

TGGATGGTGAAGAGGTTCGC

GeneID_100702234 R

GTCCAGGAGACGTTGACAGG

GeneID_100701693 F

ACATCTGGCCCTGACACAAC

GeneID_100701693 R

AACACAGAGGGTCCAACACC

GeneID_100703853 F

CAACACCTGGATCAGAGCGA

GeneID_100703853 R

GCATGGATGAGTCCCAGTCA

GeneID_100702693 F

ATCCAGGTGTTGGACCCTCT

GeneID_100702693 R

GCTGCACTCGTTTCCTTTGA

GeneID_100701975 F

CTCCTCAGACTCGTCCATGC

GeneID_100701975 R

TTGGGCCTTCCTCCTGTAGT

GeneID_109202867 F

AGAGGAACAAGCTCGCCATC

GeneID_109202867 R

TCGGGCCTTCCTCTTGTAGT

GeneID_100698680 F

TATCCCATGGACGAGGCAGA

GeneID_100698680 R

CCGTCTGAGGATACGGGTTG

GeneID_109204231 F

AGGGACGGAAAGCAGACAAC

GeneID_109204231 F

GCTGCACCAACAGCAATCTT

GeneID_109204230 F

TTGTATACACGGTGCCCAGG

GeneID_109204230 R

CACTGGATGGCCGTTTTTGG

GeneID_109204294 F

CGACATTTACAGCTGCACCG

GeneID_109204294 R

GTTACGACTCCCAGCAGACC

GeneID_109194343 F

CATGTTCATGTGCAGCGTGT

GeneID_109194343 R

GCTGATGTTCTCTCCGGGTT

GeneID_102080724 F

CGCTCACATGAGAGACCCAA

GeneID_102080724 R

TGTCCCCTTTAGCACACCAA

GeneID_100693666 F

GACTCACGGTGGAGAGTTGG

GeneID_100693666 R

GAATAGGGGGAACGGCGAAA

GeneID_100703681 F

TGAGCTGGCTGAGAGATGGA

GeneID_100703681 R

TCCAGACCTGGGTGTGTACT

GeneID_102081732 F

GCTCCCGGTGATCCTCTTTC

GeneID_102081732 R

GGCACATGTAGGCGTACTCA

GeneID_102079294 F

GTCAGACCACACGGCTCTAC

GeneID_102079294 R

AAGCAGGAGAGACGAACACC

miR-192-x-F

ACACTCCAGCTGGGATGACCTATGAATTGACAG

miR-215-x-F

ACACTCCAGCTGGGATGACCTATGAATTGACAG

oni-miR-10856-F

ACACTCCAGCTGGGCTGATCCAGGACCATCTC

oni-miR-192-F

ACACTCCAGCTGGGATGACCTATGAATTGAC

oni-miR-7133-F

ACACTCCAGCTGGGGATGTTGAGTATCAAACT

novel-m0443-3p-F：

ACACTCCAGCTGGGATTTCTGTGAAACTCGGC

novel-m0060-3p-F

ACACTCCAGCTGGGTGTTCAGACTCTCTGTGT

oni-miR-235-y-F

ACACTCCAGCTGGGTATTGCACTTGTCCCGGCC

oni-miR-92b-F

ACACTCCAGCTGGGTATTGCACTCGTCCCGGC

oni-miR-541-x-F

ACACTCCAGCTGGGAAGGGATTCTGATGTTGGTC

oni-miR-7133-y-F

ACACTCCAGCTGGGTAGTTTGATACACAGCACA

novel-m0344-5p-F

ACACTCCAGCTGGGAACTCCAGTTATCAAGGGC

oni-miR-10556-x-F

ACACTCCAGCTGGGGCAAAGAAAAGCGTCTG

oni-miR-329-y-F

ACACTCCAGCTGGGAACACACCCAGCTAACCTT

oni- miR-122-y-F

ACACTCCAGCTGGGAACGCCATTATCACACTA

novel-m0075-3p-F

ACACTCCAGCTGGGAAACAGAAAGTGAGATTGT

novel-m0134-3p-F

ACACTCCAGCTGGGAAACAGAAAGTGAGATTGT

miR-466-x-F

ACACTCCAGCTGGGTGTGTGTGTGTGTGTGTG

oni-miR-10556b-F

ACACTCCAGCTGGGAAAGAAAAGCGTCTGGAC

oni-miRNA-R

CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGAGAAGTCC

U6-F

CATCTACAGGTGCTGCTAAGG

U6-R

TCTTGAGTCTGCAGGTAAGATC

Table 2 MicroRNAs exhibiting differential expression levels in the T-24 and T-120 groups compared to NC control group miRNA

Trend of expression NC vs T-24

NC vs T-120

GeneID_109194299

Down

-

GeneID_109194454

-

Down

GeneID_100701969

Down

Down

GeneID_102081732

-

Down

GeneID_100702234

-

Up

GeneID_100701693

Up

Down

GeneID_100703853

Up

Down

GeneID_100702693

Up

Up

GeneID_100701975

-

Up

GeneID_109202867

-

Up

GeneID_100698680

Up

Down

GeneID_109204231

-

Up

GeneID_109204230

-

Down

GeneID_109204294

-

Down

GeneID_109194343

-

Down

GeneID_102080724

Up

Up

GeneID_100693666

Up

Up

GeneID_100703681

Up

Up

GeneID_102081732

-

Down

GeneID_102079294

-

Up

“-” means no change.

Table 3 mRNA-miRNA pairs changed in the opposite expression trends mRNA

Trend of expression

miRNA

Trend of expression

oni-miR-7133 miR-192-x GeneID_109204230

Down

oni-miR-192

Up

miR-215-x oni-miR-10856 GeneID_109194343

Down

oni-miR-92b

Up

miR-235-y miR-10556x novel-m0344-5p GeneID_109194454

Down

miR-7133-y

Up

oni-miR-10856 oni-miR-7133 miR-192-x GeneID_109204294

Down

oni-miR-192

Up

miR-215-x oni-miR-10856 GeneID_100701969

Down

oni-miR-10856

Up

GeneID_100702234

Up

novel-m0060-3p

Down

GeneID_109202867

Up

miR-541-x

Down

GeneID_100703681

Up

novel-m0060-3p

Down

GeneID_100701693

Up

novel-m0060-3p

Down

High lights 1. In this study, we first identified the mRNA and miRNA expression profiles of TiLV-infected tilapia. 2. We identified approximately 5,000 differently expressed mRNAs, most of which are associated with immune responses. 3. In this study, approximately 14,000 mRNA-miRNA pairs showed opposite expression profiles, including 20 validated miRNAs and their target genes. 4.miRNA can regulate the host's immune response by changing the expression of target genes after TiLV infection.