White spot syndrome virus impact on the expression of immune genes and gut microbiome of black tiger shrimp Penaeus monodon

The gut microbiome plays an essential role in the immune system of invertebrates and vertebrates. Pre and pro-biotics could enhance the shrimp immune system by increasing the phenoloxidase (PO), prophenoloxidase (ProPO), and superoxide dismutase activities. During viral infection, the host immune system alteration could influence the gut microbiome composition and probably lead to other pathogenic infections. Since the JAK/STAT pathway is involved in white spot syndrome virus (WSSV) infection, we investigated the intestine immune genes of STAT-silenced shrimp. During WSSV infection, expression levels of PmVago1, PmDoral, and PmSpätzle in PmSTAT-silenced shrimp were higher than normal. In addition, the transcription levels of antimicrobial peptides, including crustinPm1, crustinPm7, and PmPEN3, were higher in WSSV-challenged PmSTAT-silenced shrimp than the WSSV-infected normal shrimp. Meanwhile, PmSTAT silencing suppressed PmProPO1, PmProPO2, and PmPPAE1 expressions during WSSV infection. The microbiota from four shrimp tested groups (control group, WSSV-infected, PmSTAT-silenced, and PmSTAT-silenced infected by WSSV) was significantly different, with decreasing richness and diversity due to WSSV infection. The relative abundance of Bacteroidetes, Actinobacteria, and Planctomycetes was reduced in WSSV-challenged shrimp. However, at the species level, P. damselae, a pathogen to human and marine animals, significantly increased in WSSV-challenged shrimp. In constrast, Shewanella algae, a shrimp probiotic, was decreased in WSSV groups. In addition, the microbiota structure between control and PmSTAT-silenced shrimp was significantly different, suggesting the importance of STAT to maintain the homeostasis interaction with the microbiota.


Effect of PmSTAT silencing on immune relate genes
The total RNA of P. monodon intestines was extracted and subjected to qRT-PCR to investigate the expression of immune-related genes, including JAK/STAT, Toll, IMD, cytokine, phenol oxidase, and antimicrobial peptide. The transcription levels of Vago, an IFN-like antiviral cytokine, responded to PmSTAT silencing and WSSV challenge (Fig. 1B). The Vago1 was strongly up-regulated in PmSTAT dsRNA, PmSTAT dsRNA + WSSV and WSSV challenged group. It is worth noting that the expression level of Vago1 in PmSTAT dsRNA + WSSV shrimp was significantly higher than that in PmSTAT dsRNA and WSSV-challenged group. In contrast, Vago4 was up-regulated in normal shrimp infected by WSSV but down-regulated in PmSTAT-silenced shrimp infected by WSSV. As shown in Fig. 1B, PmSTAT silencing significantly enhanced Vago5 transcription level, but the WSSV challenge reduced Vago5 transcripts in both normal and PmSTAT-depleted shrimp. www.nature.com/scientificreports/ In PmSTAT-silenced shrimp, the expression levels of the JAK/STAT genes, including PmDOME and PmJAK, were decreased. In constrast, their expression levels were increased upon WSSV infection (Fig. 1C). Regarding the Toll pathway, silencing of PmSTAT enhanced the PmSpätzle expression level more than two-fold, compared to normal shrimp but showed fewer effects on the transcription levels of MyD88, Castus, and Dorsal. The expression level of PmSpätzle in PmSTAT-silenced shrimp was dropped upon WSSV infection (Fig. 1C). Meanwhile, the transcription levels of Relish, a transcription factor in the IMD pathway, in normal and PmSTAT-silenced shrimp slightly increased in response to WSSV infection (Fig. 1C). In this work, the expression level of the inhibitor of kappa B kinase, which plays an essential role in the IKK-NF-ƘB signaling cascade, was also investigated. The expression level of PmIKKβ was increased in response to PmSTAT depletion and WSSV infection (Fig. 1C). PmIKKε1 transcript was increased in both non-infected and WSSV-infected PmSTAT-silenced shrimp, compared to the normal shrimp, while PmIKKε2 transcript was decreased.

Shrimp intestinal microbiome
The dominant phyla in all samples were Proteobacteria, followed by Bacteroidetes, Actinobacteria, and Planctomycetes ( Fig. 2A). However, the relative abundance of phyla Bacteroidetes, Actinobacteria, and Planctomycetes were reduced in WSSV and PmSTAT silenced intestine. Interestingly, Proteobacteria was significantly increased in WSSV-challenged shrimp (P < 0.01). The rest of the phyla, including Bacteroidetes, Actinobacteria, Planctomycetes, Firmicutes, Verrucomicrobia, GN02 (Gracilibacteria), Chloroflexi, and TM (Saccharibacteria), were significantly reduced in WSSV challenged and also in the PmSTAT dsRNA + WSSV intestine. The highest relative abundance was Photobacterium at the genus level, followed by Vibrio and Shewanella (Fig. 2B). The Photobacterium was increased in both WSSV-challenged shrimps. Meanwhile, the relative abundance of Vibrio was similar in all conditions, and P. damselae was significantly increased in both WSSV-challenged shrimps (Fig. 2C). The Linear discriminant analysis effect size (LEfSe) analysis showed that the most LDA scores in each condition, including PBS, PmSTAT dsRNA, WSSV and PmSTAT dsRNA + WSSV was genera Cohesibacter, Shewanella, Roseivirga and Marivita, respectively (Fig. S2). Moreover, potential shrimp probiotics were identified in our samples following the previously reported method 45 . Overall, the relative abundance considering all potential probiotics was not significantly different between the four conditions ( Supplementary Information, Fig. S4). The potential probiotics with abundance in our samples at the genus level were Bdellovibrio, Phaeobacter, Pseudoalteromonas, Rhodobacter, Shewanella, Streptomyces, and Vibrio. While at the species level, we found Phaeobacter sp. DCSW07, Shewanella algae, and Vibrio hepatarius ( Supplementary Information, Fig. S5). Interestingly, when we analyzed the abundance of each of the potential probiotic species separately, we found that Shewanella algae significantly increased in PmSTAT dsRNA compared to all other three conditions. On the contrary, S. algae were depleted in the WSSV group (Fig. 2D); however, the differences were insignificant.
The alpha diversity metrics were calculated to investigate the difference between richness (Chao1 and the observed species index) and diversity (Phylogenetic diversity (PD) and Shannon diversity index) among groups ( Table 1, Rarefaction curves are shown in Fig. S6). The unchallenged groups (PBS and PmSTAT silenced shrimp) showed higher richness and diversity. Contrarily, the challenged groups (PmSTAT silencing with WSSV and WSSV infected) showed lower richness and diversity. Overall, species richness and diversity were reduced in WSSV challenged shrimps and different from the PBS and PmSTAT silencing shrimps which were similar. Table 1. Alpha diversity metrics for all the treatments. Species richness was measured using Chao2 and observed species index. Species diversity was calculated using Phylogenetic diversity (PD) and the Shannon diversity index. Different superscript lowercase letters in the same column indicate significant differences among conditions (P < 0.01).

Chao1
Observed www.nature.com/scientificreports/ The differences in intestinal microbial communities among groups were analyzed by beta diversity. Principal coordinate analysis (PCoA) based on unweighted UniFrac distances exhibited that all samples formed four significantly separated clusters (ANOSIM p = 0.001) (Fig. 3A). The distances in the PCoA are shown in Fig. 3A. Thus, we quantitatively determined the centroids for each group of samples and then calculated the distance between the centroids for all the treatments (Table 2). With this analysis, we found several behaviors: First, the most extended separations of the distances between the control (PBS) and all the treatments suggested that the microbiota of all conditions differed significantly from that of the control. Second, the distance of PmSTAT dsRNA was more significant Vs. WSSV (0.537) and Vs. PmSTAT dsRNA Vs PmSTAT dsRNA + WSSV (0.495), suggesting that the PmSTAT dsRNA microbiota significantly differed from the WSSV and PmSTAT dsRNA + WSSV. Third, Table 2. Pairwise distance between the centroids for all the treatments.

PBS
PmSTAT dsRNA WSSV PmSTAT dsRNA + WSSV   www.nature.com/scientificreports/ the smaller distance between all comparisons was between WSSV Vs. PmSTAT dsRNA + WSSV (0.334), suggesting that the microbiota between those two conditions was more similar than the other group comparisons. This also suggested that the most substantial effect on the microbiota composition caused by WSSV infection, rather than PmSTAT dsRNA depletion. Similar clusters were also observed in the UPGMA tree of unweighted UniFrac distances, in which all WSSVinfected shrimps formed one hand of the three, in which the WSSV-infected shrimps clustered separately from the PmSTAT with WSSV (Fig. 3B). Contrarily, the non-infected shrimps created another hand (Fig. 3B), in which the PBS shrimps were separated from the PmSTAT-silenced shrimps. These findings agreed with the positions observed between the four groups in the unweighted PCoA (Fig. 3A).

Discussion
Gut microorganisms are essential in host functions, including development, nutrition, immunity, and disease resistance. However, host-pathogen interaction is still unclear. Therefore, we investigated the intestinal bacterial and immune-related transcription profile. Shrimp intestines were collected from WSSV unchallenged and challenged. In addition, we also investigated those profiles in suppressed JAK/STAT, an antiviral pathway in shrimps.
Changes in bacteria composition during WSSV infection may influence the expression of shrimp immunity. In mosquitoes, Vago function as an IFN-like antiviral cytokine. The transcription levels of LvVago4 and LvVago5 were up-regulated in WSSV-challenged hemocyte species 38 . Unlike the hemocyte, the transcription level of PmVago4 was slightly up-regulated, and PmVago5 was dramatically down-regulated in WSSV challenged intestine. However, PmVago5 was induced in PmSTAT dsRNA shrimp (Fig. 1B). Meanwhile, PmVago1 was strongly up-regulated by either PmSTAT dsRNA or WSSV challenged. It has been reported that LvTCF, a main downstream effector of Wnt signaling, regulates LvVago1. During WSSV infection in L. vannamei hemocytes, the transcription levels of LvTCF were increased. However, WSSV produced WSV083 to promote the degradation of LvTCF via the ubiquitin-proteasome pathway, which suppressed transcription levels of LvVago1 46 . Furthermore, LvIRF, an interferon regulatory factor, regulates the transcription levels of LvVago4 and LvVago5. In L. vannamei hemocytes, the expression of LvVago4 and LvVago5 were inhibited in LvIRF-silenced shrimp. Moreover, cumulative mortality and WSSV copy numbers were increased in LvVago4-or LvVago5-silenced shrimp 38 . Unlike L. vannamei hemocytes, PmSTAT silencing promoted transcription levels of PmVago1 and PmVago5 in the P. monodon intestine (Fig. 1B). Meanwhile, WSSV infection strongly induced the expression of PmVago1 and slightly induced the expression of PmVago5 (Fig. 1B). Moreover, suppression of PmSTAT reduced WSSV copy numbers and promoted the transcription levels of PmVago1 (Fig. 1A and B). This suggests that PmVago1 might be the immunity frontline against WSSV in the intestine, unlike hemocytes. Figure 1C showed that PmSTAT silencing decreased the transcription levels of PmDOME and PmJAK. When PmSTAT-silenced shrimp were challenged by WSSV, the expression of PmDOME and PmJAK was increased, suggesting that WSSV infection altered the JAK/STAT pathway. Since the Dome functionally reduced the WSSV copy number (e.g., replication), presumably, it should be increased when STAT is inhibited. However, the effect of down-regulation of PmDome in the PmSTAT-silenced intestine could suggest novel functions that need more studies.
The expression of PmDorsal, a transcription factor in the Toll pathway, and PmSpätzle were promoted in the PmSTAT dsRNA group (Fig. 1C). Two WSSV microRNA named WSSV-miRNA-N13 and WSSV-miRNA-N23 were identified 47 . They suppressed the expression of MjDorsal, resulting in the inhibition of MjALF expression. In this study, PmSTAT dsRNA promoted expression of PmDorsal in WSSV-challenged shrimp (Fig. 1C). Previous study has demonstrated that crustinPm1 is controlled through the Toll signaling pathway while crustinPm7 is mediated via both Toll and Imd pathways 48 . Meanwhile, Lvpenaeidin3a is regulated through the Toll pathway 49 . The transcription levels of crustinPm1, crustinPm7, and Pmpenaeidin3 were increased in PmSTAT silenced shrimp, while in WSSV-challenged, they were dramatically decreased (Fig. 1D). The ProPO system is one of the important shrimp immunities 50 . P. monodon hemocyte is the primary cell that express the proteins of this system, while the intestine could not sense transcription in PCR 51,52 . The PmProPO1, PmProPO2, and PmPPAE1 were expressed in low levels (Fig. 1E) and dramatically suppressed in WSSV-challenged shrimp. Interestingly, PmPPAE2 was strongly promoted in WSSV-challenged shrimp. The ProPO and PPAE are important to PO activity. Single or double ProPO silencing decreased PO activity. However, there were no significant differences in PO activity between single or double ProPO silencing 52 . Similar results were shown in single or double PPAE silencing. Moreover, single PPAE silencing promoted the expression of other PPAE 53 . The strongly enhancing PmPPAE2 might result from dramatically suppressed PmPPAE1 by WSSV. However, the expression of PmP-PAE2 was not enhanced in PmSTAT silencing during WSSV infection. Thus, the PmSTAT may be necessary to transcript PmPPAE2.
The antimicrobial peptides, including crustinPm1, crustimPm7, and PEN3, were strongly promoted in PmSTAT-silenced and PmSTAT-silenced + WSSV groups and significantly suppressed in WSSV-infected group (Fig. 1D). A previous study suggested that LvPEN3 could inhibit WSSV virion internalization into hemocytes 54 . Furthermore, FmPEN3 could hinder the growth of Micrococcus lysodeikticus 55 . The stimulation of the host immune system by RNAi could be beneficial to prevent animals from pathogen invasion.
The interaction between bacteria inside the community is complex. For example, the genome sequencing of S. algae, isolated in France, revealed bacteriocins and antimicrobial peptides 69 . The S. algae was increased, while P. damselae subsp. damselae was decreased in PmSTAT-silenced shrimp ( Fig. 2C and 2D). The antimicrobial peptides in the Toll pathway, including crustinPm1, crustimPm7, and PEN3, were strongly promoted in PmSTATsilenced groups and significantly suppressed in the WSSV infected group (Fig. 1D). Thus, these antimicrobial peptides could reduce the population of P. damselae subsp. damselae in PmSTAT-deprived shrimp, giving an advantage in defense against WSSV. These antimicrobial peptides, however, did not kill or inhibit the probiotic S. algae (Fig. 2D).
The WSSV-challenged shrimp (PmSTAT dsRNA + WSSV and WSSV) revealed low richness and diversity (Table 1). Likewise, the bacterial diversity was reduced in shrimp with white feces syndrome (WFS) challenged with V. harveyi 7,70 . Moreover, the intestines of healthy shrimps also have lower diversity than diseased shrimps with AHPND 60 . Notably, no changes in microbial diversity in shrimps infected with WSSV and cotton shrimplike disease were reported 10,71 . Thus the association of greater diversity with better host health in shrimps is, until now, under discussion 8 , and more studies are necessary to understand the role of the microbial ecosystem in health and disease in aquatic organisms. The UPGMA clustering and PCoA analysis using the unweighted UniFrac distances revealed clusters separating the PBS and WSSV-challenged shrimps (Fig. 3). A similar separation of the microbiota from the control by WSSV infection was also observed in L. vannamei 10 . Interestingly the microbiota changes were minimal when we compared the PmSTAT silencing and STAT silencing + WSSV ( Fig. 2 and Table 2). This suggests that changes in the microbiota by the PmSTAT silencing were independent of the WSSV infection.
In a murine model, respiratory viral infections induce secondary bacterial infections. The antiviral immune responses induced by influenza are associated with changes in the microbiota in the respiratory and gastrointestinal tract 72 . The inflammation of the tissues promotes the secretion of type I and II interferons (IFNs), increasing Proteobacteria and Bacteroidetes in the gut 73,74 . Meanwhile, the depletion of type I IFN-alpha/beta receptor in mice improved clearance of secondary Streptococcus pneumoniae infection during influenza infection 75 . However, virus infection also promotes probiotic bacteria. For example, Lactobacillales were enriched in HIV patients 76 . Moreover, bacteria diversity in HIV patients was higher than in the seronegative group 77 , suggesting that the host immune response shifts the host microbiome. Understanding the interaction between host-pathogenmicrobiome could improve strategies to prevent diseases.
Taken together, the shift of immune and immune-related genes could change the bacteria composition in the shrimp intestine, demonstrating the importance of host-microbiota interactions in understanding the diseases. Interestingly, only silencing the PmSTAT caused drastic effects on the microbiota structure compared to the PBS (Fig. 3A), suggesting the importance of the host gene expression to maintain the homeostasis interaction with the microbiota. However, there are limited examples of this relationship in aquatic animals. For instance, in kuruma shrimp (Marsupenaeus japonicus) the silencing of CTL33 expression led directly to intestinal dysbiosis, tissue damage, and shrimp death 78 , suggesting a fine regulation through the evolution of the intestinal microbiota homeostasis in invertebrates. In this regard, it has also been observed that scallop antimicrobial peptides and proteins are implicated in maintaining microbial homeostasis and are critical molecules in orchestrating hostmicrobiota interactions 79 . Our work opens the possibility of using gene silencing to understand the relationship between shrimp microbiota and the host in the absence and diseases.
Methods. This study was conducted under the ethical principles and guidelines according to the animal use protocol 1923021 approved by Chulalongkorn University Animal Care and Use Committee (CU-ACUC). The biosafety guidelines were reviewed and approved by the Institutional Biosafety Committee of Chulalongkorn University (SC-CU-IBC-004-2018). This study was carried out in compliance with the ARRIVE guidelines.

PmSTAT silencing shrimp
The black tiger shrimp, P. monodon, of average 3.12 ± 0.17 g bodyweight, were obtained from Charoen Pokphand Farm in Chanthaburi Province, Thailand. Shrimp were acclimated in 120L tanks at ambient temperature and maintained in aerated water with a salinity of 20 ppt for at least one week before beginning the experiments. P. monodon was injected twice with PmSTAR dsRNA (10 µg per g shrimp) at 0 and 16 h after the first injection. Shrimp's intestines were collected at 24 h post-second injection. Total RNA was extracted using Tissue Total RNA mini kit (Favorgen, Taiwan). The first-strand cDNA was synthesized using the first-strand cDNA Synthesis Kit (Thermo Fisher, USA) according to manufacturer's protocol. The knockdown efficiency was measured by qRT-PCR using EF-1α as internal control and calculated by the 2 −∆∆CT method 80 .

Sampling
Shrimp were divided into four groups, PBS-injected, PmSTAT-silenced, WSSV-challenged, and PmSTAT-silenced combined with WSSV-challenged. Shrimp were cultured in 20L tanks (1 tank per group) at ambient temperature and maintained in aerated water with a salinity of 20 ppt. Each group contained nine shrimp. In the PmSTAT silenced groups, shrimp were injected with PmSTAT dsRNA 10 µg per gram shrimp body weight, and then after 16 h, they were injected with the same amount of dsRNA. Meanwhile, PBS was injected into the shrimp control set. Approximately 6 × 10 6 viral copies of WSSV (quantification described in Fernando et al. 81 ) were injected 24 h after the 2nd PmSTAT dsRNA injection. Shrimp's intestines were collected at 24 h post-WSSV challenge and immediately processed DNA and RNA extraction. Each group contained nine shrimps, and each replicate was pooled from 3 shrimp intestines.

Total RNA and quantitative real-time PCR
Total RNA was isolated by the Tissue Total RNA mini kit (Favorgen, Taiwan) according to the manufacturer's protocol. An equal amount of total RNA (500 ng) from each sample was used for cDNA synthesis using Firststrand cDNA Synthesis Kit (Thermo Fisher). Quantitative real-time RT-PCR using specific primer pair (Supplementary Table 4) for PmSTAT and WSSV immediate-early gene 1 (IE1) was employed to confirm the PmSTAT silencing efficiency, WSSV-challenged and immune-related gene response during WSSV infection by qRT-PCR. Real-time RT-PCR was carried out using an equal amount of cDNAs (2 µl of tenfold diluted cDNA) in the iCycler iQTM Real-Time detection system and the Luna® Universal qPCR Master Mix (Bio-Rad, USA). The qRT-PCR conditions were 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s and 55 °C for 10 s. The qRT-PCR was done in triplicate. The fold difference of mRNA transcription was calculated by 2 −∆∆CT method. Statistical analysis was carried out using the IBM SPSS Statistics 22 program with one-way ANOVA followed by a post hoc test (Tukey's). The result differences were considered significant at P < 0.05.

Sample preparation for microbiome study
Total DNA from the intestines was isolated for microbiota study. To this end, the intestines of nine shrimps have been collected from each treatment group, pooling three intestines as a single sample. The DNA of whole intestines was extracted using the DNA by Quick-DNA Fecal/Soil Microbe Miniprep Kit (Zymo Research, USA). The Qubit™ dsDNA HS (Invitrogen, USA) assay was performed to quantify the DNA concentration. The DNA integrity was confirmed using 1% agarose gel electrophoresis. The V3-V4 regions were amplified from genomic DNA using the following conditions: 95 °C for 3 min (denaturation), followed by 25 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, and a final extension at 72 °C for 5 min. The PCR product was analyzed using 2% agarose gel by electrophoresis and purified using Ampure XP beads (Bechman Coulter). All purified amplicon samples were mixed in equal concentrations and sequenced on the Illumina Miseq sequencing platform using 250 base pair (bp) paired-end (PE) chemistry.

Microbiome analysis
The results from sequencing were filtered and analyzed by Quantitative Insights Into Microbial Ecology (QIIME (version 1.8), http:// qiime. org/ index. html) as previously reported (Caporaso et al., 2010). The filtered sequences were assigned to the same operational taxonomic units (OTUs) with 97% similarity to the Green Genes reference sequence collection. Therefore, the taxonomic information was assigned to the Green Genes reference database. The OTUs at family, genera, and species levels were subjected to a LEfSe analysis to obtain the significantly different taxonomies among treatments with a significant level (alpha) of 0.05 and the LDA threshold > 2. For all data, statistical significance was set at p < 0.05. Alpha diversity index, including Chao1, Observed species, Phylogenetic diversity (PD), and Shannon diversity index, were calculated via QIIME to estimate diversity within the samples. The significant difference in alpha diversity index between groups was carried out using the IBM SPSS Statistics 22 program with one-way ANOVA followed by a post hoc test (Tukey's). The result differences were considered significant at P < 0.01. Unweight UniFrac distances were calculated using the QIIME program and visualized www.nature.com/scientificreports/ using PCoA analysis to estimate the beta diversity. In addition, UPGMA (Unweighted Pair Group Method with Arithmetic mean) was analyzed by QIIME using the beta diversity distance matrix to compare between conditions. Finally, Anosim was performed to measure the differences in the bacterial profiles between groups. The resulting P-value indicated a significant difference between the two groups (P-value < 0.01). The distance between centroids of the PCoA (Fig. 3A) was determined with R using the "usedist" package and the matrix obtained of the beta diversity analysis with the unweighted UniFrac distances. To find the presence of probiotics in our microbiome samples, we followed the analysis used by Ochoa-Romo et al. 45 . First, we obtained a new BIOM file from the joined sequences and clustered it with a 97% identity level against the Silva132 database. This database was used because it contains the reference of most probiotic bacteria. Afterward, the relative abundance of 70 probiotics taxa was extracted from the BIOM file using the summa-rize_taxa_throught_plots.py and an in-house script. Lastly, the Wilcoxon test was performed to determine the significant differential abundance of probiotics among tested groups.

Data availability
All data generated or analysed during this study are included in this published article and its supplementary information file. www.nature.com/scientificreports/