Virus-derived small RNAs in the penaeid shrimp Fenneropenaeus chinensis during acute infection of the DNA virus WSSV

Abstract

Small interfering RNAs (siRNAs) and microRNAs (miRNAs) are two classes of small RNAs (sRNAs) that are critical for virus-host interplay via the RNA interference (RNAi) pathway. One virus-derived siRNA and numerous miRNAs has been reported for the double-stranded DNA virus white spot syndrome virus (WSSV), however, the expression profiles of these different types of sRNAs have not been assessed. Here, by sequencing the sRNAs and mRNAs of WSSV-infected Chinese shrimp (Fenneropenaeus chinensis), we found that the viral transcripts were universally targeted by WSSV-derived siRNAs, supporting a pivotal role for RNAi in the anti-viral immunity of shrimp. The genesis of WSSV-derived siRNAs was associated with long RNA structures. Moreover, by separating miRNAs from siRNAs, 12 WSSV miRNAs were identified. Investigation of conserved viral miRNA targets in different host species indicated the involvement of viral miRNAs in host immune responses. Collectively, our data provide new insights into the role of the RNAi pathway in the interplay between DNA viruses and crustaceans.

Virus-derived small RNAs in the penaeid shrimp Fenneropenaeus chinensis during acute infection of the DNA virus WSSV

OPEN

R

A

P

Chengzhang Liu, Fuhua Li,, Yumiao Sun, Xiaojun Zhang,, JianboYuan, HuiYang & Jianhai Xiang,

Small interfering RNAs (siRNAs) and microRNAs (miRNAs) are two classes of small RNAs (sRNAs) that are critical for virus-host interplay via the RNA interference (RNAi) pathway. One virus-derived siRNA and numerous miRNAs has been reported for the double-stranded DNA virus white spot syndrome virus Here, by sequencing the sRNAs and mRNAs of WSSV-infected Chinese shrimp (Fenneropenaeus chinensis), we found that the viral transcripts were universally targeted by WSSV-derived siRNAs, supporting a pivotal role for RNAi in the anti-viral immunity of shrimp. The genesis of WSSV-derived indicated the involvement of viral miRNAs in host immune responses. Collectively, our data provide new insights into the role of the RNAi pathway in the interplay between DNA viruses and crustaceans.

The RNA interference (RNAi) pathway degrades mRNAs or inhibits their translation via gene silencing. Silencing occurs through RNA-induced silencing complexes (RISCs) consisting of an Argonaute protein (AGO) and a guiding small RNA (sRNA) that is antisense to the mRNA target1. In eukaryotes ranging from yeasts to humans, RNAi is an evolutionarily conserved mechanism that regulates gene expression and protects against viruses and transposons. Three classes of sRNAs are involved in RNAi: microRNAs (miRNAs) with a length of 1824 nt, small interfering RNAs (siRNAs) of 2123 nt, and PIWI-interacting RNAs (piRNAs) of 2431 nt. These types of sRNAs guide the RNAi process in similar ways but employ dierent AGO proteins2.

The origins of siRNA and miRNA dier. miRNA biogenesis starts in the nucleus, wherein a primary miRNA (pri-miRNA) is transcribed and processed into shorter hairpin RNAs (precursor miRNAs or pre-miRNAs) by Drosha. A pre-miRNA is then exported to the cytoplasm and processed by Dicer to produce miRNA duplexes3.

However, in most animal groups, siRNAs originate in the cytoplasm, wherein endogenous or exogenous dsRNAs are cleaved by Dicer2. In Drosophila, subsequent to the activity of Dicer, the RNAi pathways of sRNA duplexes from two dierent origins intersect: duplexes with central mismatches are sorted for AGO1, whereas highly paired duplexes are loaded into AGO2, thus forming dierent RISC complexes4,5.

Previous studies conducted in shrimp have identied components of the RNAi pathway, such as Drosha, Pasha, TRBP, Dicer1, Dicer2, AGO1, AGO2, and AGO3, revealing an RNAi system similar to that of insects69.

Concomitantly, many shrimp miRNAs have been studied10, with an emphasis on those involved in innate immunity against viruses, such as miR-1001113. To date, no reports have described shrimp piRNAs, perhaps due to the lack of available shrimp genome data.

Virus infection may induce the production of virus-derived small RNAs (vsRNAs) including siRNAs, miRNAs and/or piRNAs in diverse animals14. Various viruses express miRNAs to manipulate their own gene expression or that of the host to facilitate their life cycles. On the other hand, RNAi provides vital anti-viral defense for the host14. In this defense, viral dsRNA is recognized by Dicer2 and cleaved into vsiRNAs (virus-derived small

Key Laboratory of Experimental Marine Biology, Institute of Oceanology, Chinese Academy of Sciences, Qingdao Laboratory for Marine Biology and Biotechnology, Qingdao National Laboratory for Marine Science

SCIENTIFIC REPORTS

1

www.nature.com/scientificreports/

interfering RNAs)15. These vsiRNAs are then loaded into an RISC to guide slicing of the target viral RNAs by an AGO2 present in the complex. These vsiRNAs are characterized by continuous coverage of genomic regions that are distinct from the discrete distribution of miRNAs1618. Anti-viral RNAi is particularly important for plants and animals that lack a protein-based adaptive immune system19. It has been suggested that due to evolutionary pressure from viruses, Dicer in insects has diversied into Dicer1 and Dicer2 to reduce competition for Dicer between endogenous pre-miRNAs and exogenous viral dsRNAs20. Moreover, dsRNA-specic Dicer2 and AGO2 have become the most rapidly evolving Drosophila genes during the evolutionary arms race between Drosophila and rapidly evolving viruses21. vsiRNA was originally thought to be conned to RNA viruses with dsRNA genomes or replication intermediates that are targeted by host Dicer. However, genome-wide siRNA production has recently been discovered in DNA viruses that infect plants and insects17,18,2225. In shrimp, natural

siRNAs have also been reported, but only for the vp28 gene of white spot syndrome virus (WSSV)26. To our knowledge, genome-wide viral siRNA production has not been documented in any crustacean to date.

WSSV is a distinct DNA virus that has both economic and academic signicance. With the capability to infect a wide range of crustaceans, including shrimp, crab, and craysh, WSSV is a major pathogen in shrimp aquaculture and has caused serious losses worldwide since the 1990s27. It diers profoundly from all presently known viruses and is the sole representative of the Nimaviridae family28,29. Its circular double-stranded DNA (dsDNA) genome is among the largest of all animal viruses (305 kb). One of the rst animal DNA viruses shown to produce natural siRNAs26, it has also been reported to express numerous miRNAs30,31. However, for WSSV, natural siRNAs were only discovered for the vp28 gene26, and it has been suggested that small RNA deep sequencing may provide more detailed insights to this issue32. Moreover, although the deep-sequencing technique has the potential to discriminate between miRNAs and siRNAs33, their dierences have not been addressed in WSSV.

Thus, we aimed to determine whether WSSV genes other than vp28 are targeted by shrimp RNAi, and to clarify the proportions of miRNA and siRNA of WSSV. For these purposes, we conducted an RNA-seq study of sRNAs and mRNAs from the Chinese Shrimp (Fenneropenaeus chinensis) during latent infection (LI) and acute infection (AI) with WSSV. To identify the origin of viral sRNAs, the sRNAs were compared to the genome and mRNAs of WSSV. Bioinformatic analyses were then carried out to classify the viral sRNAs (Fig.1).

Results

In pleopod DNA from cephalothorax samples, the WSSV copy number was 119 times higher in AI shrimp than in LI shrimp ((2.931.37)105/ng and (2.460.23)103/ng, respectively, n = 10 each). Consistent with the DNA data, the mRNA diversity and expression also increased dramatically in WSSV, with 111 transcripts detected during the AI stage and none detected during the LI stage34.

Genome-wide viral sRNA production during acute infection. Aer removing sRNAs that aligned with host mRNAs, the LI cephalothorax sample provided 1,664 unique sRNAs that aligned with the WSSV genome. During the AI stage, the expression of WSSV-derived sRNA increased abruptly by 188-fold (Table1).

There were 545 unique WSSV-derived sRNAs found exclusively in LI cephalothoraxes, while 130,518 were found only in the AI stage. Additionally, 1119 viral sRNAs were present in both stages, accounting for 67% of the viral sRNAs in the LI stage but only 0.85% of those in the AI stage. The accumulated expression of these overlapping sRNAs was up-regulated 61-fold during the AI stage. None of the viral sRNAs detected in the LI stage were signicantly down-regulated during the AI stage.

Viral sRNAs of 1835 nt in length were found to be ubiquitously distributed along the WSSV genome (Fig.2C). In the gures, we refer to the viral strand that is transcribed from le to right as the plus strand and to its complement as the minus strand. A distinct sRNA hotspot was identied at 151 kb157.5kb on the minus strand. The strand-specic density of unique sRNAs across the genome provided an average value of 0.22 (Table1). Thus, on average, an unique sRNA was found on either strand at intervals no larger than 5 nt. The majority (82%) of the 305-kb viral genome was encompassed by sRNAs on either the plus (55%) or minus (58%) strand, and the sequencing depth per nucleotide ranged from 01,463 for the plus strand (average = 10) and from 035,587 (average=31) for the minus strand.

WSSV mRNAs were selectively processed to generate viral sRNAs. The ubiquitous production of viral sRNAs has raised the question of whether they are vsiRNAs or just normal mRNA degradation products. This was assessed from three perspectives.

First, the sRNA reads were classied to clarify their major components. The clean reads were annotated by comparison with the Rfam non-coding RNA database and the miRBase miRNA database. The majority of the sRNAs belonged to known sRNA categories (Supplementary Fig. S1). Considering that host piRNAs, endogenous siRNAs, and novel miRNAs cannot be identied based on sequence similarity, the proportion of genuine sRNAs should be even higher than the current ndings. These results indicated that the majority of the total sRNAs were not degraded mRNAs.

Second, to evaluate whether the unidentied short RNAs were degradation products, they were mapped to protein coding sequences (CDSs) from the transcriptome data of the same samples (see Methods). If the mRNAs were degraded, both host and viral coding mRNAs should have produced abundant short RNA fragments. The sRNA distribution frequency in the WSSV CDSs (median = 141.9 reads per kilobase of transcript per million reads (RPKM)) was 53.2 times higher than that in the host CDSs (2.67 RPKM, t-test, P=0; F-test, P=0) (Fig.3). This suggests that sRNAs were generated from viral mRNAs at productive ratios much higher than those from host mRNAs, which means viral sRNAs cannot be attributed to normal degradation of mRNAs. In other words, they originated from selective processing of viral transcripts.

Third, the distribution of sRNA lengths indicated that the WSSV-derived sRNAs were products of Dicer processing. The distinct enrichment of unique sRNAs with a length of 2123 nt is a characteristic of Dicer processing

SCIENTIFIC REPORTS

2

www.nature.com/scientificreports/

Infected shrimp ssue

RNA-seq

RNA-seq

23,360,150 sRNA reads

WSSV genome

64,188,426 Illumina reads

mapping

Trinity

WSSV sRNAs

46,813 unigenes

Non-coding regions

No

con nuously

distributed

Blast

mapping

Yes

46,702 shrimp unigenes

111 WSSV unigenes

srnaloop miRdeep

miRANDA

RT-PCR

2,656 target genes

12 WSSV miRs

ORF nding

mapping

mapping

111 WSSV ORFs

10 pri-miRs

WSSV-derived siRNAs

Figure 1. Roadmap of the bioinformatics analysis.

(Fig.4E,F), in stark contrast to the unique host sRNAs, which predominantly displayed lengths of 2627 nt (Fig.4AC), a very small proportion of which were miRNAs (Supplementary Fig. S1).

In brief, the combined information for the length distribution and classication of sRNAs as well as the level of sRNA production from coding sequences suggested that the majority of WSSV-derived sRNAs were not degraded mRNAs but were, in fact, Dicer products.

Relationship between the distribution of viral sRNAs and mRNAs. Because RNAi provides an anti-viral mechanism by selectively degrading viral mRNAs, which results in the production of viral siRNAs, a correlation between the genomic distribution of viral sRNAs and mRNAs is expected. In fact, regions that lacked mRNA expression usually also lacked sRNA production (Fig.2B,C, Supplementary Fig. S2C,D), indicating that viral sRNAs were derived from mRNAs. A weak correlation was observed between the expression of viral sRNA and mRNA (Supplementary Fig. S3, R= 0.52), suggesting that factors other than the quantity of mRNA also aected the genesis of sRNA. The viral transcripts that matched the most sRNAs were located at ~155 kb in the virus genome, including wsv271 (VP136) and wsv277 (Fig.5). A second sRNA hotspot residing at ~220 kb was represented by wsv360 (VP664). Hotspots for sRNAs of dierent lengths usually coincided in terms of the genomic location. However, at certain locations, their proportions were quite dierent, suggesting the presence of dierent classes of sRNAs (Fig.2C, Supplementary Fig. S2 D).

Genomic characteristics of WSSV-derived sRNAs. Whether WSSV-derived sRNAs were miRNAs, piRNAs or siRNAs was what we wanted to nd out next.

We rst evaluated the sRNAs for piRNA characteristics. piRNAs are generally longer than siRNAs or miRNAs35. Known characteristics of piRNAs include an overlap of 10 nt between the 5 ends of the sense and antisense

SCIENTIFIC REPORTS

3

www.nature.com/scientificreports/

Category Stages Parameter 2123 nta 1835 nt

Count 204 1,664 Densityb 0.00033 0.0027 2123 nt/1835 nt 0.12

Count 34,904 131,637 Density 0.06 0.22 2123 nt/1835 nt 0.27 Acute/latent Fold 171 79

Unique sRNA

Latent

Acute

Latent RPMc 23 194

2123 nt/1835 nt 0.12

Acute RPM 8719 36631

2123 nt/1835 nt 0.24 Acute/latent Fold 385 188

Table 1. WSSV-derived sRNA expression in the cephalothoraxes of F. chinensis at LI and AI stages.

a2123nt represents sRNAs in the size range of 2123 nucleotides, which is the typical length of siRNA. 1835 nt stands for sRNAs in the range of 1835 nucleotides, which covers the length ranges of siRNA, miRNA and piRNA. bThe strand-specic density here was calculated as the number of unique WSSV-derived sRNAs divided by twice the length of the WSSV genome. cRPM: WSSV-derived sRNA reads per million sRNA reads.

Total sRNA

A Coding structures

CDS in the plus strand

CDS in the minus strand

mRNA contigs from the RNA-seq transcriptome

B mRNA reads (log)

0510

Ln(coverage + 1)

C sRNA reads (log)

505

1820nt 2123nt 2425nt 2631nt 3235nt

0 50000 100000 150000 200000 250000 300000

Genome position (nt)

Figure 2. Distribution of sRNAs and mRNAs along the WSSV genome from cephalothoraxes of F. chinensis during acute WSSV infection. (A) Coding structures of the WSSV genome. Arrows represent annotated coding sequences from NCBI. (B) Coverage of the mRNA reads (orientation unknown) on a logarithmic scale. (C) Coverage of sRNA reads on a logarithmic scale. Bars with a coverage less than 0 represent sRNAs in the reverse direction. Red arrows between (B,C) indicate regions where coldspot of mRNA and sRNA match each other.

piRNAs, as well as a strong sequence bias for uridine at position 1 (1 U bias) in piRNAs bound by PIWI and a bias for adenosine at position 10 (10A bias) in piRNAs bound by AGO336. However, the majority of WSSV-derived sRNAs did not display the above characteristics of piRNAs (Supplementary Figs S4 and S5). Moreover, the typical length of piRNAs is 2431 nucleotides, but the WSSV-derived sRNAs were predominantly 21~23nt in length (Fig.4E,F,K,L). This is a feature of Dicer products, indicating that the majority of sRNAs were either miRNAs or siRNAs.

SCIENTIFIC REPORTS

4

www.nature.com/scientificreports/

cephalothorax intestine

A

WSSV CDS

420

B

WSSV CDS

659

1000

618

368

394 1000

577

sRNA production (RPKM)

341

536

100

315

494

289

100

453

263

412

237

371

10

210

330

10

184

289

158

248

132

207

1

106

1

166

80

124

53

83

0.1

27 0.1

42

0.1 1 10 100 1000 10000

Host CDS counts

1

0.1 1 10 100 1000 10000 100000

Host CDS counts

1

mRNA CDS (FPKM)

mRNA CDS (FPKM)

cephalothorax

intestine

C

sRNA production (RPKM)

0100200300400

D

050100200300

Host WSSV

Host WSSV

Figure 3. Comparison of sRNA production from coding sequences (CDSs) of WSSV and host genes in the cephalothoraxes and intestines of F. chinensis during acute infection. For each CDS, the production of sRNAs has been normalized by the length of the CDS. (A,B) Scatter plots in logarithmic style. (C,D) Box plots. The sRNA production in the WSSV CDSs (median=141.9 RPKM) was 53.2 times higher than that in the host CDSs (2.67 RPKM). The dierence were very signicant (t-test, P=0; F-test, P=0).

The WSSV-derived sRNAs were then assessed for miRNA characteristics. Because a hairpin-like RNA structure is required for processing of pre-miRNAs by Dicer1,2, we scanned the non-coding region of the WSSV genome for pre-miRNAs and found 1,255 potential pre-miRNAs. Because miRNAs are generated from the stem regions of pre-miRNAs, the sRNA frequencies were compared between the potential pre-miRNA structures, the non-coding genomic regions, and the whole viral genome (Fig.6). We found that sRNA frequency was signicantly higher in the stems of hairpin structures than in other regions (Wilcoxon rank-sum test, W= 1.0 1010,

P=2.21016), indicating the production of miRNAs. However, potential miRNAs (sRNAs generated from hairpin stems) made up only 3.6% of all the WSSV-derived sRNAs at the AI stage. Taken together, the above information suggested that viral miRNAs were produced during AI, although they only made up a minor portion of the WSSV-derived sRNAs.

Because miRNAs and piRNAs were deemed unlikely to comprise the major component of WSSV-derived sRNAs, the sRNAs was hypothesized to be primarily siRNAs. One distinction between siRNA and miRNA is that siRNAs can be generated from viral ORFs23, while viral miRNAs are usually generated from non-coding regions37,38. We aligned the viral sRNAs to the WSSV ORFs from the transcriptome data, and found that 78.0% of the total WSSV-derived sRNAs originated from protein-coding sequences. This suggested that most were not miRNAs but siRNAs. Our sequencing data is consistent with previously reported vp28 siRNAs26. For the AI stage, 1,817 sRNA reads mapped to the vp28 CDS region (Supplementary Fig. S6), and of these, 22 reads overlapped with the published vp28 probes (Supplementary Text S1). In addition, the WSSV-derived sRNAs showed much higher sequence diversity than the host sRNAs, which were mainly miRNAs. Although the WSSV-derived sRNAs made up only 1.4% of the total AI sRNAs, they accounted for 14% of the unique sRNAs. High sequence diversity is a characteristics of siRNAs, because unlike miRNAs, siRNA genesis is not conned to pre-miRNA structures.

Although DNA viruses lack a dsRNA intermediate that could serve as an siRNA precursor, there may be other sources of dsRNA. Genesis of virus-derived siRNA from overlapping transcripts has been proposed by

SCIENTIFIC REPORTS

5

www.nature.com/scientificreports/

Cephalothorax (latent) Cephalothorax (acute ) Intestine (acute)

A B C

D E F

G H I

J K L

Host Unique

RPM

02000600010000

18 20 22 24 26 28 30 32 34

18 20 22 24 26 28 30 32 34 18 20 22 24 26 28 30 32 34

WSSV Unique

RPM

050010001500

18 20 22 24 26 28 30 32 34

18 20 22 24 26 28 30 32 34 18 20 22 24 26 28 30 32 34

RPM

0100000250000

Host Total

18 20 22 24 26 28 30 32 34

18 20 22 24 26 28 30 32 34 18 20 22 24 26 28 30 32 34

WSSV Total

RPM

01000200030004000

18 20 22 24 26 28 30 32 34

18 20 22 24 26 28 30 32 34 18 20 22 24 26 28 30 32 34

Host both strand WSSV plus strand WSSV minus strand

Figure 4. Length distribution of WSSV and host sRNA reads. Total and unique are abbreviations for total reads and unique reads, respectively.

Bronkhorst et al. in a DNA virus23. This explanation is consistent with the overlapping nature of the ORFs in both strands of the WSSV genome. However, unlike Bronkhorst et al., we did not identify a strong correlation for the siRNA density between the two strands of individual ORFs (Supplementary Fig. S7). This indicated that convergent transcripts were not the major source of WSSVderived siRNAs. siRNAs can also be derived from hairpin structures via a pathway that comprises Dicer2, Argonaute 2 and Loquacious, yielding another potential source of viral siRNA39,40. Because the hairpin structures that generate siRNAs can be much longer than pre-miRNAs, we examined the association between the sRNA frequency and long RNA structures. The RNA structure of the 6.5-kb sRNA hotspot in the minus strand at 151kb157.5kb in the WSSV genome was predicted using RNAfold (Fig.7A). More than 60 percent (3,924) of the nucleotides were located in the double-stranded region. The WSSV-derived sRNAs from the sequencing data were mapped to the structures. The number of nucleotides in the predicted double strands was calculated for each sRNA, and their relationship to the sequencing counts for the same sRNA was plotted (Fig.7B). sRNAs with more paired nucleotides tended to have higher copy numbers. Signicant copy number dierences were identied between sRNAs (average count=12.0) with greater numbers (12) of paired nucleotides and those (average count = 7.5) with lower numbers (11) of paired nucleotides (t=1.2109, t-test with unequal sample sizes and unequal variances; W=1.1109, P=1.6105, Wilcoxon

SCIENTIFIC REPORTS

6

www.nature.com/scientificreports/

Figure 5. Distribution of sRNAs and mRNAs at the hotspot at ~155kb in the WSSV genome in cephalothoraxes of F. chinensis during acute infection. (A) Coding structures. Arrows represent coding sequence annotation from NCBI. (B) Sequencing coverage of mRNA reads (orientation unknown). (C) Sequencing coverage of sRNA reads. Bars with a coverage of less than 0 represent sRNAs in the reverse direction. The lengths of the sRNAs are indicated in dierent colors.

A B

Total sRNA

Whole genome

Noncoding

Hairpin stems

Unique sRNA

Coverage

0.00.51.01.52.02.5

Coverage

0.00.20.40.60.8

20 25 30 35

20 25 30 35

sRNA length (nt)

sRNA length (nt)

Figure 6. sRNA (18-24 nt) expression in dierent regions of the WSSV genome during acute infection in the cephalothoraxes of F. chinensis. (A) Total sRNA coverage. (B) Unique sRNA coverage. An illustration of hairpin structures is shown in the legend, with dierent regions illustrated in the same colors used in (A,B). Total sRNA frequency was signicantly higher in hairpin structures than in other regions (Wilcoxon rank-sum test, W=1.01010, P=2.21016).

rank-sum test). The same results were also found for sRNAs obtained from an independent shrimp intestine sample (Fig.7C). Thus, the production of WSSV-derived sRNAs is related to long hairpin structures.

Briey, the majority of WSSV-derived sRNAs were most likely siRNAs, but their genomic characteristics also suggested the presence of miRNAs as a minor component.

Taking the existence of siRNA into consideration, we rened the procedure for miRNA identication to reduce false-positive results by excluding candidates that overlapped with adjacent

SCIENTIFIC REPORTS

7

www.nature.com/scientificreports/

Figure 7. Relationship between sRNA copy numbers and RNA structures at the sRNA hotspot at ~155kbof WSSV genome during acute infection. (A) Predicted RNA structures. (B,C) Relationship between the total counts of an sRNA and the number of paired nucleotides corresponding to the sRNA. Unique sRNAs (1824 nt) were mapped to the RNA secondary structure in (A). The mapped site for each individual sRNA was then counted for the paired nucleotides. Nucleotides forming double-stranded bonds were counted, while those in bulge structures were excluded. Using this strategy, the number of paired nucleotides for an sRNA was obtained. This number was related to the copy number of the same unique sRNA, which was obtained by sequencing. Signicant dierences were found for the copy numbers between sRNAs (average count=12.0) with higher (12) paired nucleotides and those (average count=7.5) with lower (11) paired nucleotides (t=1.2109, t-

test with unequal sample sizes and unequal variances; W=1.1109, P=1.6105, Wilcoxon rank-sum test).

sRNAs (see classification of viral sRNAs in the Methods). Fourteen WSSV miRNAs were predicted from the sequencing data. To conrm the expression of these viral miRNAs, RT-PCR assays were conducted using stem-loop RT primers specic for detection of mature miRNAs41 (a loop structure in the primer prevents its hybridization with miRNA precursors, Supplementary Dataset S1). The expression of 12 of the 14 candidates was conrmed (Fig.8, Table2). The hairpin structures of their pre-miRNAs were conserved in the 4 WSSV strains (Supplementary Fig. S8; Text S2). Ten of the miRNAs were found to be novel, whereas 2 overlapped with reported WSSV miRNAs: 11 nt of wssv-m5 overlapped with wssv-miR-N131, and 16 nt of wssv-m4 overlapped with wssv-miR-21130.

Expression of the viral miRNA in shrimp stomach and lymphoid organ was examined by RT-PCR. The results showed that the expression of wssv-m1 and m12 were much higher than others (Fig.8A). Expression of most

SCIENTIFIC REPORTS

8

www.nature.com/scientificreports/

Figure 8. WSSV miRNA expression at 24 and 48h post-injection in the stomach and the lymphoid organ of F. chinensis. (A) The viral miRNA expression levels compared with 18S rRNA. (B) Fold change of the miRNAs expression levels compared with their own levels at hour 0. OKA and ST stands for lymphoid organ and stomach, respectively.

Genome Position

Strand Position 5/3

wssv-m1 GUUGUUGUUCAUGUUGAGGGCAUU 152840~152863 5 33576 wssv-m2 AGAAGAGGACUUUGGGGUAGAC 79223~79244 5 31474 wssv-m3 UGACGAUAAGAGGAGGCAGUAG + 236406~236427 5 33547

wssv-m4

GAAAUGGCUGUCCAGAAAUCUGG

wssv-m5 CAGGUCUCUUCUAGGUAGAAG 41679~41699 5 39470 wssv-m6 CUGGAGAGACUAUAGGGGAAUA + 173945~173966 3 33544 wssv-m7 GAAAGGUUGUUGAUGGAGAUG 290022~290042 3 33511 wssv-m8 GCGUCUCUUGAGUCUUGUAGUA 250561~250582 5 37970 wssv-m9 UAAGAGAAUAAGGAGAAUCUGAG + 246949~246971 3 39453 wssv-m10 GUGACUGGAGUGUUGAUUGGGU + 248192~248213 5 / wssv-m11 UAACAUUAUCGUGGAACAGAAA 43233~43254 3 39471 wssv-m12 GGGGUUUCGUUUGUGGUUG 154722~154740 5 33576

Table 2. WSSV miRNAs veried by RT-PCR. aTranscriptome sequences from Li et al.34. The prex Unigene was omitted. The sequences can be found in Text S3 in the Supplemental material.

of the viral miRNAs continued to increase up to 48 hours post injection (hpi) (Fig.8B). The expression levels of wssv-m5 and m6 were signicantly higher in the lymphoid organ than in the stomach at 24h, while at 48h, their tissue expression patterns were the opposite of those observed at 24 h. Clustering of the variations in miRNA expression suggested that wssv-m7, m10 and m11 had the most similar expression patterns.

Viral miRNA target analysis. Twenty-two WSSV genes were found to be potentially regulated by viral miRNAs, among which 14 were targeted by multiple miRNAs, indicating potential co-regulation of these miRNAs (Supplementary Table S1). wssv-m1 had the most viral targets (11), followed by wssv-m6 (7 targets), indicating that these miRNAs might play pivotal roles in viral gene regulation. WSSV miRNAs and their predicted viral target genes were compared among the genomes of 4 dierent WSSV strains from the NCBI. All of the miRNAs and their target sites were conserved among the strains. A comparison of 2 samples (AI cephalothorax and AI intestine) from two independent WSSV outbreaks was also performed. The results showed that all of the miRNA sequences were conserved. Examination of the viral target sequences showed that some target sites were missing due to alternative splicing or lack of expression, but no nucleotide substitutions were found at the target sites.

For the host genes, preliminary prediction using miRanda42 recovered 2,656 mRNAs as potential viral miRNA targets, indicating a broad inuence on host gene expression.

To screen out pivotal target genes in the virus-host interactions, evolutionary data were explored. Considering that WSSV can infect a wide range of hosts, the regulation of a pivotal target gene by WSSV miRNAs should be conserved in dierent host species. Therefore, we searched for shared target genes of WSSV miRNAs in three host species, including two Penaeidae shrimp (F. chinensis and L. vannamei) and a Palaemonidae shrimp

ID Sequence

Pri-miRNA IDa

272861~273501

273826~272883

+ 273523~273848

5

/

SCIENTIFIC REPORTS

9

www.nature.com/scientificreports/

(Exopalaemon carincauda). All three species are susceptible to WSSV infection. The results revealed that 66 target genes were shared by all three hosts (Supplementary Dataset S2; Fig. S9A), among which 22 had binding sites that were evolutionarily conserved (Supplementary Dataset S3; Fig. S9B).

Among the 22 host genes that contained conserved binding sites for WSSV miRNAs, 15 have been reported to be involved in innate immunity, the cell cycle, or other virus-related processes. The PL10A/Belle/DDX3 dead-box helicase subfamily is associated with antiviral signaling and RNAi43, and it is targeted by diverse viruses for inhibition44. p120-catenin regulates a wide range of biological processes, including cell proliferation, inammation and innate immunity45. The 14-3-3 gene participates in several signal transduction pathways, and it has been reported to be involved in anti-viral responses in shrimp46 and humans47. Eukaryotic translation initiation factor 5A (eIF5A) is implicated in cell cycle control and has been associated with WSSV infection48. The predicted regulation will be veried in future studies.

Visualization of the viral miRNA regulation network. To gain a better understanding of the intricate relationships between viral miRNAs and their targets, the interactions were visualized with Cytoscape. Gene annotation and expression data were also visualized (Supplementary Fig. S9). As indicated by the symbol size, it is readily apparent that the most abundant WSSV miRNAs were wssv-m1 and m12, while the most abundant host target gene was MLR (myosin light chain 2). Host genes that were targeted by multiple miRNAs were placed in the middle of the network, which revealed that the most strongly down-regulated mRNAs (blue circles) were targeted by multiple viral miRNAs, indicating possible co-regulation. As indicated by the thick line, the strongest regulation of a host gene was predicted between wssv-m1 and CTND2 (Supplementary Fig. S9A). All of the WSSV genes were up-regulated during acute infection, as indicated in red (Supplementary Fig. S9C).

Discussion

In recent years, virus-derived sRNAs (vsRNAs) have been detected in animals infected with various viruses14,18.

Interestingly, the RNAi mechanisms associated with vsRNAs seems to be quite dierent between insect and mammalian hosts49. Supplementary Table S2 summarizes some recent studies of vsRNAs in animal viruses. In some cases, vsRNAs are ubiquitously and continuously produced, such that they can be assembled into long genomic contigs and even an entire viral genome17. Most of these discoveries were associated with RNA viruses. For DNA viruses, ubiquitous sRNA production has only been reported in insects22,24. The present study expands these ndings to include a crustacean. In previous reports, the classication of vsRNAs was conducted according to their genomic characteristics and experimental evidence. In a few studies, piRNAs were identied in the presence of an abundance of 2430-nt sRNAs displaying strong 5 terminal uridine uridine bias and 10 nt adenine bias. In most of the other vsiRNA reports, viral sRNAs were classied as siRNAs because their continuous coverage of genomic regions was distinct from the discrete distribution of miRNAs14. In the cases of RVFV, VSV and WSSV, both virus-derived miRNAs and siRNAs were generated during the virus-host interaction. It is believed that viral miRNAs were encoded and expressed by the virus to regulate viral or host gene expression, while virus-derived siRNAs were generated via host RNAi response during which viral RNAs were detected and subjected for degradation14,50.

In the present study, adjacent siRNAs were considered during the identication of miRNAs, which has not yet been performed for WSSV. Following the miRbase protocol51, we successfully recovered known miRNAs (wssv-miR-211 and wssv-miR-N1) and siRNAs (for vp28). The biogenesis and functionality of these siRNAs and miRNAs have been experimentally validated previously26,30,31, supporting the reliability of our methods. Nevertheless, the sample types used in the sequencing studies here were limited. Additional sequencing studies should be carried out in the future using multiple tissue types at more infection stages, preferably from a dierent host species.

Consistent with most of the results of genome-wide viral sRNA analyses, the majority of WSSV-derived sRNAs were also suggested to be siRNAs. Canonical miRNAs also contributed to WSSV-derived sRNA production but only to a minor proportion of the sRNAs. Nonetheless, these viral miRNAs were dramatically up-regulated during the AI stage (Fig.8), indicating their relevance to this stage. Although WSSV miRNAs were less abundant than WSSV-derived siRNAs, their amounts might still be present in sufficient quantities for target-regulation. Exploiting conserved miRNA targets from dierent host species enabled us to screen for important candidate targets of WSSV miRNAs, which indicated their roles in evading host innate immune responses such as antiviral signaling, inammation and apoptosis. On the other hand, the autoregulatory WSSV miRNAs may contribute to the control of virus replication. These hypotheses should be veried in future studies.

One noteworthy characteristic of WSSV-derived sRNAs is their dramatic change in expression along the viral genome. The presence of sRNA hotspots has previously been reported for other viruses23,24,52. In this study, the presence of a distinct hotspots at 151.0~157.5 kb in the minus strand of the WSSV genome was conrmed in independent samples (Figs2 and S2). Analysis showed that the production of sRNA from secondary structures of single-stranded RNAs (ssRNAs) might have produced the sRNA hotspot. This is consistent with previous studies which found that viral siRNAs can be derived from long hairpin structures via a pathway that comprises Dicer2, AGO2 and Loquacious39,40.

Our results expand the scope of WSSVderived vsiRNA discovery from vp28 to the level of the whole genome. The ndings of genome-wide siRNAs in a crustacean are consistent with the results obtained for insects. Taken together, these ndings suggest that RNAi-mediated defense against DNA viruses may be common among arthropods. Studies of the functional role of siRNAs from individual WSSV genes should be carried out in the future, using mutant animals23 or by knocking down Dicer2 and/or AGO226, as well as siRNA inhibitors26.

In contrast to viral infection of insect vectors, which show considerable resistance53, WSSV infection of the penaeid shrimp is virulent and usually lethal. This indicates that despite the degradation of viral mRNAs by host

SCIENTIFIC REPORTS

10

www.nature.com/scientificreports/

RNAi, which results in a large amount of viral siRNA production, WSSV can still overwhelm the host defenses, supporting the presence of viral countermeasures to the host RNAi system. The substantial co-expression of WSSV mRNAs with sRNAs during acute infection supports this speculation (Fig.2). It has been proposed that a rapidly replicating virus can saturate the RNAi machinery54. However, if a virus can replicate so rapidly that viral RNA saturates the RNAi system, host RNAi should have already been futile against the virus prior to saturation. An alternative theory is that the virus uses junk RNAs as a decoy to divert the RNAi machinery so that the core part of the viral genome or transcripts achieve relative protection. This hypothesis is supported by the nding that hotspot vsiRNAs possess poor antiviral activity compared with coldspot vsiRNAs52. From a dierent perspective, Jayachandran et al. (2012) argued that the degradation of viral transcripts by host RNAi may be advantageous to the virus by preventing its over-replication, which may lead to the premature death of host cells24.

Previous studies have also identied various families of Viral Suppressors of RNA Silencing (VSR) proteins19.

One VSR class, called siRNA sequesters binds to siRNAs and prevents their incorporation into RISC, which may protect viral mRNAs and promote the accumulation of viral siRNAs. This provides a potential explanation for the co-existence of WSSV mRNA and siRNA observed in the present study. We searched for homologs to siRNA sequesters in WSSV proteins, but we only identied sequences with low similarity. This was not unexpected because VSRs with similar functions can be strikingly diverse as a consequence of their independent evolution19.

The discovery of genome-wide vsiRNAs in shrimp only raises questions regarding the complex mechanisms underlying their genesis and their roles in the context of RNAi against a DNA virus.

Materials and Methods

The ow chart outlines the applied bioinformatic analysis technique (Fig.1).

Preparation of shrimp and RNA samples. Cephalothorax samples of WSSV-carrying F. chinensis collected during the LI and AI stages of infection, prepared as described in our previous transcriptome study34, were used for the transcriptome and sRNA sequencing. Briey, WSSV-carrying shrimp obtained from a local farm were cultured in tanks, and individuals in the LI and AI stages were classied according to WSSV copy numbers55.

The cephalothoraxes of 10 LI shrimps and another 10 AI individuals were sampled, frozen, and preserved in liquid nitrogen until DNA or RNA extraction. An independent RNA sample was also obtained from the intestines of AI F. chinensis collected at another local farm at a dierent time, following the same procedures described above.

The shrimp used for the time series and organ-specic miRNA RT-PCR have been described in our previous study55. In short, F. chinensis received an intramuscular injection of 800 copies of WSSV, and tissues were dissected at dierent hours post injection (hpi) and immediately preserved individually in liquid nitrogen. Stomachs and lymphoid (Oka) organs from 3 shrimp at 0, 24 and 48 hpi were pooled for RNA extraction.

Total RNA was extracted using Unizol reagent (Sangon, China), treated with DNase I, veried by electrophoresis in a 1% agarose gel and quantied with a NanoDrop 1000 spectrophotometer (NanoDrop, USA).

Transcriptome sequencing. Transcriptome sequencing was performed as previously described34. Briey, Polyadenylated (polyA+) RNA was puried using Sera-mag oligo (dT) beads, fragmented, reverse-transcribed using random hexamers, end-repaired, and adaptor-ligated. Products of 300500 bp were separated, amplied, and pair-end sequenced using a Genome Analyzer II (Illumina, San Diego, USA) at the Beijing Genomics Institute (BGI), Shenzhen, China.

sRNA sequencing. Ten micrograms of total RNA was used to construct the sRNA libraries with the Small RNA Sample Prep Kit (Illumina, San Diego, USA). Briey, sRNA fragments with an approximate length of 1835 nt were isolated and puried following 15% denaturing polyacrylamide gel electrophoresis (PAGE). Subsequently, 3 and 5 RNA adaptors were ligated to the sRNAs using T4 RNA ligase. RT-PCR was performed to reverse transcribe the adaptor-ligated sRNAs into cDNAs, which were further amplied for 15 PCR cycles. The products were puried using 10% PAGE and checked for size, purity, and concentration with an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, USA). The sRNA libraries were sequenced using an Illumina/Solexa 1G Genome Analyzer (Illumina, San Diego, USA) at the BGI. The LI and AI sRNA sequencing data have been submitted to the SRA database at the NCBI and are now available under accession numbers SRX807107 and SRX806880, respectively.

Classification of host sRNAs. The clean sRNA reads were annotated by comparison with the Rfam non-coding RNA database (http://xfam.org/) and the miRBase miRNA database (http://mirbase.org/), using Blastn. The majority of total sRNAs were successfully classied into known sRNA categories, which made up 67%, 60% and 70% of the total sRNAs from the LI cephalothorax, AI cephalothorax, and AI intestine, respectively (Supplementary Fig. S1).

Viral sRNAs were identied by comparison with WSSV genome and transcriptome sequences.

Four WSSV genomes were downloaded from NCBI (GI: 19481591, 426202315, 417410589 and 58866698). Aer removing adaptor or low-quality sequences, sRNAs with a length of 1835 nt were mapped to the WSSV genome using bowtie2 (http://bowtie-bio.sourceforge.net) and SAMtools (http://samtools.sourceforge.net/) to separate WSSV and host sRNAs. Because the mapping results were very similar between the 4 versions, the most recently updated genome, GI417410589 was chosen for subsequent analyses. The genomic distribution of mapped vsiRNAs was visualized using R (https://www.r-project.org/).

Transcriptome sequences from the same sample were annotated as previously described. Transcripts matching plants or microbes other than WSSV were excluded from the analysis. Ribosomal RNA transcripts were also removed. sRNAs were also mapped to annotated transcripts using bowtie2 and SAMtools. The sRNA expression

SCIENTIFIC REPORTS

11

www.nature.com/scientificreports/

of each transcript was calculated as RPKM and then integrated according to the NR annotation to summarize the expression of transcripts from the same gene.

Comparison of sRNA production from host and WSSV CDSs. To evaluate whether the unclassied short RNAs (Supplementary Fig. S1) were degradation products, they were mapped to protein coding sequences (CDSs) from the transcriptome sequencing data obtained from the same samples34. The CDSs were obtained by aligning mRNAs to known proteins in the NCBI NR database, using Blastx. For AI cephalothorax, 120 WSSV CDSs and 17,240 host CDSs (FPKM0.1) were compared for unique sRNAs mapping to the sense strand of each CDS. A considerable proportion (27.6%) of the host CDSs showed no sRNA production; by contrast, 100% of the WSSV CDSs displayed sRNA production. sRNA distribution frequency in the CDSs was calculated in RPKM. Dierence in sRNA frequencies between host CDSs and WSSV CDSs were compared using t-test and F-test.

Genomic characterization of virus-derived sRNAs. We refer to the viral strand that is transcribed from le to right as the plus strand and to its complement as the minus strand.

At rst, the viral sRNAs were evaluated for miRNA characteristics. Sequences from the coding and non-coding regions of the WSSV genome were obtained according to the gene structure annotation (GI417410589, NCBI), using in-house Perl (https://www.perl.org/) scripts. pre-miRNA-like hairpin structures along the non-coding region of the WSSV genome were predicted with the srnaloop algorithm (http://arep.med.harvard.edu/miRNA/). The stability of the structures was evaluated with RNAfold (http://www.tbi.univie.ac.at/RNA/) based on calculations of minimum free energy (MFE<25kcal/mole). Then the sequences of both arms and the loop regions of the structures were retrieved with Perl scripts. The sRNA frequencies in the hairpin arms, loops, and other genomic regions were compared, for which the regions were trimmed by 24 nt to exclude sRNAs that stretched across the borders. For sRNAs in the range of 1824 nt (the length of miRNAs), the dierences were checked using the Wilcoxon rank-sum test and t-test, using R.

In order to assess the viral sRNAs for piRNA characteristics, the nucleotide preferences of sRNAs and the distance between the 5 ends of the sense and antisense sRNAs were visualized using R.

To verify the association between the sRNA frequency and RNA secondary structures in long sequence ranges, the RNA secondary structures of the sRNA hotspot at 151157.5 kb in the minus strand of the WSSV genome were predicted with RNAfold (MFE<25kcal/mole).

Stacks of multiple oset sRNA reads (2123 nt) covering a continuous region were classied as possible siRNAs33 if the dierence in their 5 start positions was >5 nt. Next, two strategies, the miRDeep2 and the srnaloop methods, were combined to identify viral miRNAs among the remaining sRNAs. Initially, the miRDeep2 pipeline analysis was performed according to Friedlander et al.56.

Briey, all of the clean reads were mapped to the genomes, and potential pre-miRNA structures mapped by stacking 1824-nt sRNAs were excised from the genome and assessed in terms of their energetic stability by RNAfold to identify pre-miRNA hairpins. If a read fell within a pre-miRNA-like structure resembling a pair of Dicer-processed miRNAs and its miRNA-star, a score was given to reect the likelihood of the pre-miRNA being genuine. To recover miRNAs that might have been overlooked by miRDeep2, another procedure was also adopted. This procedure was initiated by ltering for pre-miRNA-like structures along the non-coding region of the WSSV genome using the srnaloop algorithm and by retrieving sequences in both hairpin arms. A stack of 1824 nt sRNA reads (copy number 2) residing in a hairpin arm was considered a possible miRNA if there were no adjacent overlapping sRNA stacks. Repetitive sequences were removed from the miRNA candidates. To distinguish novel miRNAs from known miRNAs, the predicted miRNA was compared with all registered miRNAs in miRBase and with previously reported WSSV miRNAs that were not registered in miRBase30,57. Finally, the remaining 2123-nt sRNAs that aligned with the WSSV genome were considered to be potential viral siRNAs.

Quantitative real-time PCR. Mature miRNAs were rst reverse-transcribed with a looped RT primer41 (Supplementary Dataset S1, LC-Sciences, Houston, USA) complementary to the 3 end of the miRNA, and the subsequent cDNA was quantied by real-time PCR using Platinum SYBR Green qPCR SuperMix (Thermo Fisher, Waltham, USA) with the ABI PRISM 7900HT system (Thermo Fisher, Waltham, USA)41. The expression levels

were normalized to 18S rRNA and evaluated using the comparative CT method. Statistical analyses were performed using paired Students t-tests.

Viral miRNA target analysis. Unlike plant miRNAs, which bind to the coding regions of their target mRNAs, animal miRNA binding sites are usually conned to the 3 UTR58. Binding sites in 3 UTRs have been reported to show high regulatory efficacy, while those in ORFs and 5 UTRs have only marginal or no efficacy58.

Thus, we focused on 3 UTRs for miRNA target prediction. Two strategies were used to predict WSSV miRNA targets, namely the miRanda strategy42 and the Targetscan strategy58. The miRanda strategy was used to assess the potential scope of the inuence of viral miRNAs on host genes, whereas the Targetscan strategy was performed to identify functionally important target genes that are conserved among dierent host species or virus strains.

For the miRanda strategy, mRNA data obtained from the same samples as the miRNA data were used as the target background, as described previously34. The coding structures of the mRNA transcripts were analyzed using a combination of 2 algorithms. First, the transcripts were compared to the sequences deposited in the NCBI protein database (Nr) to identify orthologous coding sequences, using Blastx (E=1e5). Second, for transcripts with no known homology, ESTscan (http://estscan.sourceforge.net/) was used to predict the coding structure. The 3 UTR sequences were retrieved and evaluated for miRNA binding sites using miRanda42, which requires strict alignment in the seed region (miRNA oset positions 28).

SCIENTIFIC REPORTS

12

www.nature.com/scientificreports/

For the Targetscan strategy, the 3 UTR data from three host species and four WSSV strains were used to predict conserved target genes. For the host genes, in addition to the F. chinensis data, the transcriptome sequences of Litopenaeus vannamei and Exopalaemon carincauda assembled in our previous study were used59,60. Homologous transcripts were identified by blast alignment (E = 1e5) between unigenes of the three hosts, requiring bi-directional best hit (BBH). Coding structures were analyzed for 46,813, 66,815, and 68,693 transcriptome unigenes from F. chinensis, L. vannamei, and E. carincauda, resulting in 6,257, 6,031, and 9,322 3 UTRs, respectively. Homologous 3 UTR sequences were identied in 602 genes, and were searched for miRNA targets. For viral genes, homologous 3 UTRs were extracted from four WSSV genomes in the NCBI database (GI: 19481591, 417410589, 426202315, and 58866698). mRNAs without NCBI (Nr) annotation were excluded. 3 UTRs longer than 25 nt were retrieved and aligned using MUSCLE (http://www.drive5.com/muscle/). Conserved 3 UTR regions were searched for miRNA binding sites using Targetscan58.

Visualization of the viral miRNA regulation network. A viral miRNA regulatory network containing 12 WSSV miRNAs and their target genes was constructed. The expression levels and fold-changes of the miRNAs and their target genes during AI in shrimp cephalothoraxes were integrated into the network, which was visualized with Cytoscape (http://cytoscape.org/). The miRNAs and their target genes are represented by dierent symbols. Fold changes in the expression levels of miRNAs and targets during acute WSSV infection are indicated by the color gradient. The absolute expression level of an miRNA (RPKM) or a target gene (FPKM, fragments per kilobase of transcript per million mapped reads) is indicated by the size of the symbol.

References

1. Fire, A. et al. Potent and specic genetic interference by double-stranded RNA in Caenorhabditis elegans. Nature 391, 806811, doi: 10.1038/35888 (1998).2. Kim, V. N., Han, J. & Siomi, M. C. Biogenesis of small RNAs in animals. Nature Reviews Molecular Cell Biology 10, 126139, doi: 10.1038/nrm2632 (2009).3. Tang, G. L. siRNA and miRNA: an insight into RISCs. Trends in Biochemical Sciences 30, 106114, doi: 10.1016/j.tibs.2004.12.007 (2005).

4. Forstemann, K., Horwich, M. D., Wee, L., Tomari, Y. & Zamore, P. D. Drosophila microRNAs are sorted into functionally distinct argonaute complexes aer production by Dicer-1. Cell 130, 287297, doi: 10.1016/j.cell.2007.05.056 (2007).

5. Tomari, Y., Du, T. & Zamore, P. D. Sorting of Drosophila small silencing RNAs. Cell 130, 299308, doi: 10.1016/j.cell.2007.05.057 (2007).

6. Su, J. et al. A key gene of the RNA interference pathway in the black tiger shrimp, Penaeus monodon: Identication and functional characterisation of Dicer-1. Fish & Shellsh Immunology 24, 223233, doi: 10.1016/j.fsi.2007.11.006 (2008).

7. Chen, Y. H. et al. Isolation and characterization of cDNAs encoding Ars2 and Pasha homologues, two components of the RNA interference pathway in Litopenaeus vannamei. Fish & Shellsh Immunology 32, 373380, doi: 10.1016/j.fsi.2011.11.032 (2012).

8. Wang, S., Chen, A. J., Shi, L. J., Zhao, X. F. & Wang, J. X. TRBP and eIF6 Homologue in Marsupenaeus japonicus Play Crucial Roles in Antiviral Response. PLos One 7, doi: 10.1371/journal.pone.0030057 (2012).

9. Unajak, S., Boonsaeng, V. & Jitrapakdee, S. Isolation and characterization of cDNA encoding Argonaute, a component of RNA silencing in shrimp (Penaeus monodon). Comparative Biochemistry and Physiology B-Biochemistry & Molecular Biology 145, 179187, doi: 10.1016/j.cbpb.2006.07.002 (2006).

10. Ruan, L. W. et al. Isolation and identication of novel microRNAs from Marsupenaeus japonicus. Fish & Shellsh Immunology 31, 334340, doi: 10.1016/j.fsi.2011.05.023 (2011).

11. Huang, T., Xu, D. & Zhang, X. Characterization of host microRNAs that respond to DNA virus infection in a crustacean. Bmc Genomics 13, doi: 10.1186/1471-2164-13-159 (2012).

12. Yang, L., Yang, G. & Zhang, X. The miR-100-mediated pathway regulates apoptosis against virus infection in shrimp. Fish & Shellsh Immunology 40, 146153, doi: 10.1016/j.fsi.2014.06.019 (2014).

13. Yang, G., Yang, L., Zhao, Z., Wang, J. & Zhang, X. Signature miRNAs Involved in the Innate Immunity of Invertebrates. PLos One 7, doi: 10.1371/journal.pone.0039015 (2012).

14. Ding, S. W. RNA-based antiviral immunity. Nature Reviews Immunology 10, 632644, doi: 10.1038/nri2824 (2010).15. Ding, S. W. & Lu, R. Virus-derived siRNAs and piRNAs in immunity and pathogenesis. Curr. Opin. Virol. 1, 533544, doi: 10.1016/j. coviro.2011.10.028 (2011).

16. Donaire, L. et al. Deep-sequencing of plant viral small RNAs reveals eective and widespread targeting of viral genomes. Virology 392, 203214, doi: 10.1016/j.virol.2009.07.005 (2009).

17. Kreuze, J. F. et al. Complete viral genome sequence and discovery of novel viruses by deep sequencing of small RNAs: A generic method for diagnosis, discovery and sequencing of viruses. Virology 388, 17, doi: 10.1016/j.virol.2009.03.024 (2009).

18. Wu, Q. F. et al. Virus discovery by deep sequencing and assembly of virus-derived small silencing RNAs. Proceedings of the National Academy of Sciences of the United States of America 107, 16061611, doi: 10.1073/pnas.0911353107 (2010).

19. Ding, S. W. & Voinnet, O. Antiviral immunity directed by small RNAs. Cell 130, 413426, doi: 10.1016/j.cell.2007.07.039 (2007).20. Ghildiyal, M. & Zamore, P. D. Small silencing RNAs: an expanding universe. Nature Reviews Genetics 10, 94108, doi: 10.1038/ nrg2504 (2009).

21. Obbard, D. J., Jiggins, F. M., Halligan, D. L. & Little, T. J. Natural selection drives extremely rapid evolution in antiviral RNAi genes. Current Biology 16, 580585, doi: 10.1016/j.cub.2006.01.065 (2006).

22. Ma, M. J. et al. Discovery of DNA Viruses in Wild-Caught Mosquitoes Using Small RNA High throughput Sequencing. PLos One 6, 7, doi: 10.1371/journal.pone.0024758 (2011).

23. Bronkhorst, A. W. et al. The DNA virus Invertebrate iridescent virus 6 is a target of the Drosophila RNAi machinery. Proceedings of the National Academy of Sciences of the United States of America 109, E3604E3613, doi: 10.1073/pnas.1207213109 (2012).

24. Jayachandran, B., Hussain, M. & Asgari, S. RNA Interference as a Cellular Defense Mechanism against the DNA Virus Baculovirus. Journal of Virology 86, 1372913734, doi: 10.1128/jvi.02041-12 (2012).

25. Sabin, L. R. et al. Dicer-2 Processes Diverse Viral RNA Species. PLos One 8, 16, doi: 10.1371/journal.pone.0055458 (2013).26. Huang, T. Z. & Zhang, X. B. Host defense against DNA virus infection in shrimp is mediated by the siRNA pathway. Eur. J. Immunol. 43, 137146, doi: 10.1002/eji.201242806 (2013).

27. Chou, H. Y., Huang, C. Y., Wang, C. H., Chiang, H. C. & Lo, C. F. Pathogenicity of a baculovirus infection causing white spot syndrome in cultured penaeid shrimp in Taiwan. Diseases of Aquatic Organisms 23, 165173, doi: 10.3354/dao023165 (1995).

28. Yang, F. et al. Complete genome sequence of the shrimp white spot bacilliform virus. Journal of Virology 75, 1181111820, doi: 10.1128/jvi.75.23.11811-11820.2001 (2001).29. van Hulten, M. C. W. et al. The white spot syndrome virus DNA genome sequence. Virology 286, 722, doi: 10.1006/viro.2001.1002 (2001).

SCIENTIFIC REPORTS

13

www.nature.com/scientificreports/

30. He, Y. D. & Zhang, X. B. Comprehensive characterization of viral miRNAs involved in white spot syndrome virus (WSSV) infection. Rna Biology 9, 10191029, doi: 10.4161/rna.20741 (2012).

31. Huang, T. Z., Cui, Y. L. & Zhang, X. B. Involvement of Viral MicroRNA in the Regulation of Antiviral Apoptosis in Shrimp. Journal of Virology 88, 25442554, doi: 10.1128/jvi.03575-13 (2014).

32. Bronkhorst, A. W., Miesen, P. & van Rij, R. P. Small RNAs tackle large viruses RNA interference-based antiviral defense against DNA viruses in insects. Fly 7, 216223, doi: 10.4161/y.25708 (2013).

33. Kozomara, A. & Griffiths-Jones, S. miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Research 39, D152D157, doi: 10.1093/nar/gkq1027 (2011).

34. Li, S. H., Zhang, X. J., Sun, Z., Li, F. H. & Xiang, J. H. Transcriptome Analysis on Chinese Shrimp Fenneropenaeus chinensis during WSSV Acute Infection. PLos One 8, 15, doi: 10.1371/journal.pone.0058627 (2013).

35. Siomi, M. C., Sato, K., Pezic, D. & Aravin, A. A. PIWI-interacting small RNAs: the vanguard of genome defence. Nature Reviews Molecular Cell Biology 12, 246258, doi: 10.1038/nrm3089 (2011).

36. Brennecke, J. et al. Discrete small RNA-generating loci as master regulators of transposon activity in Drosophila. Cell 128, 10891103, doi: 10.1016/j.cell.2007.01.043 (2007).

37. Umbach, J. L. et al. MicroRNAs expressed by herpes simplex virus 1 during latent infection regulate viral mRNAs. Nature 454, 780 U108, doi: 10.1038/nature07103 (2008).

38. Zhu, Y., Haecker, I., Yang, Y., Gao, S.-J. & Renne, R. gamma-Herpesvirus-encoded miRNAs and their roles in viral biology and pathogenesis. Curr. Opin. Virol. 3, 266275, doi: 10.1016/j.coviro.2013.05.013 (2013).

39. Okamura, K. et al. The Drosophila hairpin RNA pathway generates endogenous short interfering RNAs. Nature 453, 803U808, doi: 10.1038/nature07015 (2008).40. Kawamura, Y. et al. Drosophila endogenous small RNAs bind to Argonaute 2 in somatic cells. Nature 453, 793U795, doi: 10.1038/ nature06938 (2008).

41. Chen, C. F. et al. Real-time quantication of microRNAs by stem-loop RT-PCR. Nucleic Acids Research 33, doi: 10.1093/nar/gni178 (2005).

42. Enright, A. J. et al. MicroRNA targets in Drosophila. Genome Biology 5, 14 (2004).43. Pek, J. W. & Kai, T. DEAD-box RNA helicase Belle/DDX3 and the RNA interference pathway promote mitotic chromosome segregation. Proceedings of the National Academy of Sciences of the United States of America 108, 1200712012, doi: 10.1073/ pnas.1106245108 (2011).

44. Schroeder, M. Human DEAD-box protein 3 has multiple functions in gene regulation and cell cycle control and is a prime target for viral manipulation. Biochemical Pharmacology 79, 297306, doi: 10.1016/j.bcp.2009.08.032 (2010).

45. Perez-Moreno, M. et al. p120-Catenin mediates inflammatory responses in the skin. Cell 124, 631644, doi: 10.1016/j. cell.2005.11.043 (2006).

46. Wanna, W., Thipwong, J., Mahakaew, W. & Phongdara, A. Identication and expression analysis of two splice variants of the 14-3-3 epsilon from Litopenaeus vannamei during WSSV infections. Molecular Biology Reports 39, 54875493, doi: 10.1007/s11033-011-1351-8 (2012).

47. Aoki, H., Hayashi, J., Moriyama, M., Arakawa, Y. & Hino, O. Hepatitis C virus core protein interacts with 14-3-3 protein and activates the kinase Raf-1. Journal of Virology 74, 17361741, doi: 10.1128/jvi.74.4.1736-1741.2000 (2000).

48. Phongdara, A., Laoong-u-thai, Y. & Wanna, W. Cloning of eIF5A from shrimp Penaeus monodon, a highly expressed protein involved in the survival of WSSV-infected shrimp. Aquaculture 265, 1626, doi: 10.1016/j.aquaculture.2007.02.005 (2007).

49. Schirtzinger, E. E. et al. Repertoire of virus-derived small RNAs produced by mosquito and mammalian cells in response to dengue virus infection. Virology 476, 5460, doi: 10.1016/j.virol.2014.11.019 (2015).

50. Umbach, J. L. & Cullen, B. R. The role of RNAi and microRNAs in animal virus replication and antiviral immunity. Genes & Development 23, 11511164, doi: 10.1101/gad.1793309 (2009).

51. Ambros, V. et al. A uniform system for microRNA annotation. Rna-a Publication of the Rna Society 9, 277279, doi: 10.1261/ rna.2183803 (2003).

52. Siu, R. W. C. et al. Antiviral RNA Interference Responses Induced by Semliki Forest Virus Infection of Mosquito Cells: Characterization, Origin, and Frequency-Dependent Functions of Virus-derived Small Interfering RNAs. Journal of Virology 85, 29072917, doi: 10.1128/jvi.02052-10 (2011).

53. Fablet, M. Host Control of Insect Endogenous Retroviruses: Small RNA Silencing and Immune Response. Viruses-Basel 6, 44474464, doi: 10.3390/v6114447 (2014).

54. Geiss, B. J., Pierson, T. C. & Diamond, M. S. Actively replicating West Nile virus is resistant to cytoplasmic delivery of siRNA. Virology Journal 2, 13, doi: 10.1186/1743-422x-2-53 (2005).

55. Sun, Y. M., Li, F. H. & Xiang, J. H. Analysis on the dynamic changes of the amount of WSSV in Chinese shrimp Fenneropenaeus chinensis during infection. Aquaculture 376, 124132, doi: 10.1016/j.aquaculture.2012.11.014 (2013).

56. Friedlander, M. R., Mackowiak, S. D., Li, N., Chen, W. & Rajewsky, N. miRDeep2 accurately identies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Research 40, 3752, doi: 10.1093/nar/gkr688 (2012).

57. Huang, Y. J. et al. The expression of interleukin-32 is activated by human cytomegalovirus infection and down regulated by hcmvmiR-UL112-1. Virology Journal 10, doi: 5110.1186/1743-422x-10-51 (2013).

58. Grimson, A. et al. MicroRNA targeting specicity in mammals: Determinants beyond seed pairing. Molecular Cell 27, 91105, doi: 10.1016/j.molcel.2007.06.017 (2007).59. Wei, J., Zhang, X., Yu, Y., Li, F. & Xiang, J. RNA-Seq reveals the dynamic and diverse features of digestive enzymes during early development of Pacic white shrimp Litopenaeus vannamei. Comparative Biochemistry and Physiology D-Genomics & Proteomics 11, 3744, doi: 10.1016/j.cbd.2014.07.001 (2014).

60. Sun, Y. et al. Transcriptome Analysis of the Initial Stage of Acute WSSV Infection Caused by Temperature Change. PLos One 9, doi: 10.1371/journal.pone.0090732 (2014).

Acknowledgements

This work was supported by the Major State Basic Research Development Program of China (No. 2012CB114403), the Hi-Tech Research and Development Program of China (No. 2012AA092205 and No. 2012AA10A405) and the Scientic and Technological Innovation Project of Qingdao National Laboratory for Marine Science and Technology (No.2015ASKJ02). We thank Shihao Li for the transcriptome data from the same samples. We thank Jiquan Zhang for providing advices. We thank High Performance Computing Center of IOCAS for providing the computational resource.

Author Contributions

J.X. and F.L. formulated the idea of the paper and supervised the research. C.L., Y.S. and H.Y. performed the experiments and analysis. C.L. and J.Y. prepared the gures and tables. C.L. wrote the paper. F.L., X.Z. and J.X. revised the manuscript.

SCIENTIFIC REPORTS

14

www.nature.com/scientificreports/

Supplementary information accompanies this paper at http://www.nature.com/srep

Competing nancial interests: The authors declare no competing nancial interests.

How to cite this article: Liu, C. et al. Virus-derived small RNAs in the penaeid shrimp Fenneropenaeus chinensis during acute infection of the DNA virus WSSV. Sci. Rep. 6, 28678; doi: 10.1038/srep28678 (2016).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

SCIENTIFIC REPORTS

15