MicroRNA and mRNA interactions coordinate the immune response in non-lethal heat stressed Litopenaeus vannamei against AHPND-causing Vibrio parahaemolyticus

While Vibrio parahaemolyticus (VPAHPND) has been identified as the cause of early mortality syndrome (EMS) or acute hepatopancreatic necrosis disease (AHPND) in shrimp, mechanisms of host response remain unknown. Understanding these processes is important to improve farming practices because this understanding will help to develop methods to enhance shrimp immunity. Pre-treatment of shrimp with 5-minute chronic non-lethal heat stress (NLHS) for 7 days was found to significantly increase Litopenaeus vannamei survival against VPAHPND infection. To elucidate the mechanism involved, mRNA and miRNA expression profiles from the hemocyte of L. vannamei challenged with VPAHPND after NLHS with corresponding control conditions were determined by RNA-Seq. A total of 2,664 mRNAs and 41 miRNAs were differentially expressed after the NLHS treatment and VPAHPND challenge. A miRNA-mRNA regulatory network of differentially expressed miRNAs (DEMs) and differentially expressed genes (DEGs) was subsequently constructed and the interactions of DEMs in regulating the NLHS-induced immune-related pathways were identified. Transcriptomic data revealed that miRNA and mRNA interactions contribute to the modulation of NLHS-induced immune responses, such as the prophenoloxidase-activating system, hemocyte homeostasis, and antimicrobial peptide production, and these responses enhance VPAHPND resistance in L. vannamei.

. Effect of chronic non-lethal heat stress (NLHS) on the survival of VP AHPND -challenged shrimp. Shrimps were stressed under the NLHS condition or reared normally under the NH condition. After 3 days recovery, shrimp were infected with VP AHNPD by immersion. The TSB contained 1.5% NaCl (instead) for the control group. Shrimp survival was measured every 6 h post-infection (hpi) for 53 h. NLHS-VP (○) and NH-VP (■) represent NLHS-and NH-treated shrimps challenged with VP AHPND , respectively. NLHS (▼) and NH (▲) represent NLHSand NH-treated shrimp immersed with TSB containing 1.5% NaCl. The experiment was performed in triplicates, and the survival percentage was calculated as the mean ± 1 standard error (S.E.) at each time point.
All clean reads from mRNA libraries were concatenated and de novo assembled using Trinity, which generated 174,835 putative genes (unigenes) or 205,137 isotigs. The generated reference assembly has an N50 isotig length of 1,074 represented by various isotigs ranging from 201 bp to 22,966 bp (Table S2). The isotigs/transcripts were annotated by searching their sequences using BLAST against transcripts predicted from the available L. vannamei genome in NCBI Genbank, Swiss-Prot, GO, Cluster of Orthologous Groups (COG), and KEGG Pathway databases. A total of 47,401 (or 23.11%) sequences had significant hits (E-value ≤ 10) to the Swiss-Prot database and the majority of these sequences were homologous to Homo sapiens (26.54%), Mus musculus (16.77%), and Drosophila melanogaster (13.55%) genes ( Fig. 2A). BLAST2GO mapped 184,422 level 2 gene ontologies (Fig. 2B), while COGs classified 11,350 sequences into different categories (Fig. 2C). Searching against the KEGG database showed that 33,475 sequences were mapped to a KEGG orthology, but only 20,183 were grouped into the reference pathways. The metabolic pathways, biosynthesis of secondary metabolites and biosynthesis of antibiotics were among the top 20 KEGG pathways (Fig. 2D) represented in the transcriptome assembly. A protein BLAST was also completed using the predicted coding sequences from the Trinotate protocol and these annotations, along with other supplementary information such as the transmembrane regions are shown in Table S3.

Differentially expressed genes (DEGs) in NLHS-treated L. vannamei upon Vp AHpnD challenge.
In this study, we aimed to study the effect of NLHS treatment on transcription in VP AHPND -challenged shrimp. Therefore, only NH-VP-responsive genes and NLHS-VP-responsive genes were identified. Differentially expressed genes (DEGs) in NLHS-treated L. vannamei (upon VP AHPND challenge) were identified by pairwise comparisons among the relevant groups as shown in the volcano plots (Fig. S1). Three biological replicates were used in the experiments and differential gene expression is represented as fold change against a specific group (Table S4; Table S4_ summary_fold-changes_mRNAs.xlsx). Between the VP AHPND challenged groups at 0 hpi (0 NH-VP) and 6 hpi (6 NH-VP), 792 genes were differentially expressed and 272 of these genes were significantly up-regulated (FDR < 0.05) in 6 NH-VP, whereas 520 genes were down-regulated in 6 NH-VP. The gene, identified as L. vannamei Relish small isoform gene (FJ416145), had the highest up-regulation (362-fold) in this particular comparison. Between the 0 NH-VP and VP AHPND challenged group at 24 hpi (24 NH-VP), 676 genes were differentially expressed; 224 and 452 genes were up-and down-regulated, respectively, in the 24 NH-VP group. The Relish small isoform gene was also up-regulated in this group (144-fold), together with P. monodon triosephosphate isomerase gene homolog (7.4-fold). In the NLHS-treated shrimp challenged with VP AHPND at 0 hpi (0 NLHS-VP) vs 6 hpi (6 NLHS-VP) comparison, 522 genes were differentially expressed, and 262 of these were significantly up-regulated in the 6 NLHS-VP group, whereas 260 genes were down-regulated. The gene homolog of lipoprotein aminopeptidase was found to be 5.3-fold www.nature.com/scientificreports www.nature.com/scientificreports/ up-regulated in the 6 NLHS-VP group, 3.9-fold up-regulated in the 6 NH-VP group, and 3.7-fold up-regulated in the 24 NH-VP group. Between the 0 NLHS-VP and 24 NLHS-VP, 272 genes were up-regulated in the 24 NLHS-VP group, whereas 520 genes were down-regulated.

Sequence analysis of shrimp miRNAs.
To supplement the information derived from the transcriptome, we also explored the global miRNA expression to obtain a glimpse of some regulators that are associated with the observed gene expression. To do this, we analyzed the miRNAs expressed in VP AHPND -infected and NLHS-treated shrimp (mir_NLHS-VP) by sequencing them at various sampling times (0, 6, and 24 h post infection; hpi). High-throughput sequencing generated 1,086,629 total raw reads in the 0 mir_NLHS-VP, 879,272 in the 6 mir_ NLHS-VP, and 1,114,328 in the 24 mir_NLHS-VP. The high-quality sequences that passed initial quality filters were 948,089, 771,799, and 956,249 reads, respectively.
Final filtering and analysis generated a majority of the non-redundant sequences 20-22 nucleotides (nt) in length (Fig. 3A). Searching the NCBI nucleotide database revealed that, on average, 25% of the sequences are most likely contaminating RNAs (Fig. 3B). After removal of these contaminating mRNA, rRNA, and tRNA homologs, the final counts of sequences were 78,000, which were mapped to miRBase 22.1 generating 77,415 sequences with hits. The percentages of matched mature miRNA sequences were 93.86%, 93.93%, and 94.19%, respectively, for the 0 mir_ NLHS-VP, 6 mir_NLHS-VP, and 24 mir_NLHS-VP libraries. Sequences with unknown identities and homologs were also listed. Of those, forty-one miRNA homologs were identified from the NLHS-VP experimental group (Table S5).
Under the NLHS-VP condition, Relish gene expression was significantly higher in all experimental groups than the respective controls. The dynamin gene was up-regulated 2-fold at 6 hpi and down-regulated 2-fold at 24 hpi. The lipoprotein receptor gene was up-regulated nearly 2-fold at 6 hpi. Importin7, JHEH-1, DNAJ5, PO1, and PO2 were significantly down-regulated 1.5-to 10-fold at 6 hpi and 24 hpi (Fig. 4). It should be noted, the expression pattern determined from the RNA-Seq data was similar to the expression pattern of selected DEGs determined from the RT-qPCR results.
Meanwhile, stem loop RT-qPCR analysis revealed that the expression levels of all 10 chosen DEMs were altered in shrimp hemocytes following NLHS treatment and VP AHPND challenge by about 1.5-to 8-fold. For  . Validation of RNA-Seq using RT-qPCR. Eight representative genes (relish, lipoprotein receptor, dynamin, importin7, juvenile hormone epoxide hydroxylase 1; JHEH-1, DNAJ5, prophenoloxidase 1; PO1, and prophenoloxidase 2; PO2) were evaluated for their expression in hemocytes of shrimp under the NLHS and NH conditions in response to VP AHPND infection and are referred to as NLHS-VP and NH-VP, respectively. Total RNA from hemocytes of NLHS-VP and NH-VP L. vannamei at 0, 6, and 24 hpi was used for cDNA synthesis. The relative expression levels of eight genes were determined by RT-qPCR and normalized against EF-1α, the internal reference. The relative expression ratio was calculated using the 2 −ΔΔCT method. The experiments were completed using triplicates. The expression level was calculated relative to that of the normal shrimp under the NH condition at 0 h after the VP AHPND challenge. The bar graphs are the data from RT-qPCR presented as means ± standard deviations and the triangles (▲) are data from the RNA-Seq. Asterisks indicate significant difference (P < 0.05) from the respective VP AHPND infected NH shrimp at 0 hpi. (2020) 10:787 | https://doi.org/10.1038/s41598-019-57409-4 www.nature.com/scientificreports www.nature.com/scientificreports/ Correlation of DEMs and DEGs of NLHS-treated shrimp in response to VP AHpnD infection. The DEGs of 3,980 NLHS-VP-responsive genes and 3,141 NH-VP-responsive genes were compared and used to construct a Venn diagram highlighting specific groups of DEGs (Fig. 6A). A grouping of 2,664 DEGs were considered to be the NLHS-VP-responsive genes, while another grouping of 1,825 DEGs were categorized as the NH-VP-responsive genes. Their intersection with 1,316 DEGs were expressed in response to VP AHPND infection whether or not the shrimp were treated with NLHS. In the future, it will be interesting to analyze the NLHSand VP-responsive genes to broaden our knowledge regarding the shrimp immune response against NLHS and VP AHPND infection.
The DEMs of NLHS-VP analyzed in this study and of NH-VP identified in our previous work 13 were analyzed to further identify miRNAs that regulate immune genes of NLHS-treated shrimp infected with VP AHPND . As with the expression profiles of the DEMs, a Venn diagram was also created to highlight specific groupings of DEMs between the libraries of NLHS-VP and NH-VP. Eighteen DEMs were specifically grouped into NLHS-VP-responsive miRNAs and two DEMs were added to a NH-VP-responsive miRNAs group (Fig. 6A).
The DEMs from the NLHS-VP-responsive miRNAs that group only with their corresponding target mRNAs from the sequencing dataset were analyzed using CU-mir (in-house) and RNA-hybrid software. This analysis facilitated the identification and functionalizing of specific miRNA-mRNA interactions, which then served as a clue to the general regulatory mechanisms underlying the immune response of shrimp under the NLHS and VP AHPND infection. 1,833 DEM-DEG pairs with negative correlations were identified and included in a miRNA-target network ( Table S6). Some of the biological functions that might be regulated by the NLHS-VP miRNAs were found to include: "Defense & Homeostasis", "Energy & Metabolism", "Cell cycle & DNA Synthesis/ repair", "Gene expression & Protein synthesis/degradation", "Receptor", "Signaling & Communication", "Transporter", "Hypothetical protein", and "Unknown". Several miRNAs such as lva-miR-7278-5p, lva-miR-6813-5p, lva-miR-745b, lva-miR902l-5p, lva-miR-502b-3p, and lva-miR-2898 had high degrees of connectivity and might play crucial roles in the regulatory network. Meanwhile, genes involved in "Defense & Homeostasis", "Energy & Metabolism", and "Cell cycle & DNA Synthesis/repair" were the most common miRNA targets ( Fig. 6B,C).
In order to further characterize the identified biological pathways, an enrichment analysis specific for immune-related pathways was done for the identified DEGs of miRNA-mRNA pairs. Three pathways that changed significantly (P-value < 0.05) were "hemocyte homeostasis", "prophenoloxidase system", and "AMP production". Given these immune pathways and the information on the canonical members of these pathways (in related studies), we were able to find sequence homologs of these members in our RNA-Seq dataset. We then used publicly available information on shrimp as well as other sequences in our transcriptome data to further annotate these homologs. This sequence information was used to lookup gene expression patterns in the identified pathways using RT-qPCR. Figure 7 shows the expression profiles of some canonical members of pathways identified in this study. It is of note that Fig. 7 presents the expression profile of PO1 and PO2 under the NLHS condition, which is also presented in Fig. 4, which reveals that PO1 and PO2 of the prophenoloxidase system are highly up-regulated, from 1.5-to 8-fold, in the NLHS-VP group. Looking at the members of the IMD pathway, IMD, IKKβ, and Relish, all had higher gene expression in the NLHS-VP group. Meanwhile, Toll1, Toll2, Toll3, MyD88, TRAF6, Pelle, and Dorsal genes in the Toll pathway were not significantly changed in the NLHS-VP group. For the hemocyte homeostasis pathway, transglutaminase and inhibitor of apoptosis protein were analyzed and found to be highly up-regulated while the caspase gene was down-regulated in the NLHS-VP group. Confirming the results of the pathway enrichment analysis, these findings suggest that the biological immune pathways, "hemocyte homeostasis", "prophenoloxidase system", and "AMP production via IMD pathway", might play important roles in the enhancement of shrimp antibacterial immunity against VP AHPND upon modulation of NLHS (Fig. 8). The NLHS trigger could be considered a prior conditioning and preparation mechanism for the shrimp to fight later infections, thereby enhancing bacterial immunity against VP AHPND and improving survival rates.

Discussion
The AHPND is known to be caused by VP AHPND , which accumulates in the stomach of shrimp and secretes PirA/B toxins in the hepatopancreas 14 . The mechanism by which AHPND kills shrimp is currently unclear. Recent data demonstrates that the genes of Toll and IMD pathways, and their downstream antimicrobial peptides (AMPs) are suppressed in the stomach and in hemocytes, but are overexpressed in the hepatopancreas 15 . This implies that while the stomach and hepatopancreas are major VP AHPND targets, the hemocytes, being the major immune organs of shrimp, may also provide informative clues regarding the immune mechanisms against disease 16 .
The previous observation that the non-lethal heat shock enhances the production of heat shock proteins and, subsequently, increases the expression of some immune-related genes resulting in enhanced immunity 4 , is of research interest because of its potential application in developing preventive strategies for diseases in shrimp. For instance, a short-term hyperthermic treatment that is suggested to reduce gill-associated virus replication in P. monodon, may prove to be a simple and effective prophylactic strategy 17 .
In this study, we found that NLHS improved the survival of VP AHPND -infected L. vannamei (Fig. 1), revealing that the NLHS may indeed modulate immune factors to alleviate the mortality caused by AHPND. However, previous studies have shown that the tolerance and survival rate of L. vannamei after VP AHPND infection is not influenced by Hsp70 accumulation or changes in immune-related proteins, such as proPO and hemocyanin 4 . This indicates that the mechanisms of resistance to VP AHPND infection under the NLHS conditions remain unexplored and may (potentially) be explained by some undescribed immune proteins. Thus, we performed RNA-Seq and small RNA-Seq analyses of either NLHS-treated or NH control shrimp infected by VP AHPND to explore the genes, gene networks, and miRNAs that regulate these unknown immune mechanisms. www.nature.com/scientificreports www.nature.com/scientificreports/ Figure 5. Relative expression analysis of miRNAs in response to VP AHPND infection following the NLHS and NH treatments in L. vannamei hemocyte. Total small RNAs from hemocyte of VP AHPND -infected L. vannamei under NH-and NLHS-treated conditions which are NH-VP and NLHS-VP, respectively, were used as templates for specific stem-loop first strand cDNA synthesis. Relative expression levels of 10 miRNAs (lva-miR-7170-5p, lva-miR-2169-3p, lva-miR-184, lva-miR-92b-5p, lva-miR-317, lva-miR-4901, lva-miR-92a-3p, lva-miR-61, lva-miR-2898, and lva-miR-6090) were determined by RT-qPCR and normalized against U6, the internal reference, at 0, 6, and 24 hpi. The bar graphs are data from RT-qPCR presented as means ± standard deviations and triangles (▲) are data from the small RNA-Seq. The results were derived from triplicate experiments. Asterisks indicate significant differences (P < 0.05) from the respective VP AHPND -infected shrimp at 0 hpi. www.nature.com/scientificreports www.nature.com/scientificreports/ Here, sequencing of mRNAs generated approximately 400 million (400 M) reads across 18 libraries for de novo assembly, which is 4-times and 6-times higher than those of recent transcriptome libraries from the hepatopancreas 18 and stomach 19 , respectively. These 18 libraries represent 3 biological replicates with an average of 22 M reads, each with 6 samples ensuring optimum statistical power to detect differentially expressed genes 20 . We also assembled a higher number of 205,137 isotigs with higher 1,074 N 50 isotig length. Regarding the number of unigene sequences with COG and KEGG annotations, this study includes a higher number of hits compared to Rao, et al. 18 However, the percentage of BLASTx hits of our reference assembly to the SWISS-Prot database and the L. vannamei genome-transcribed data was lower (23.11% vs 26.13%) 18 . The percentage (26.89%) of unknown contigs in our assembly may indicate that we mined more novel genes than previously described, a collection of interesting genes related to immune response and bacterial tolerance to characterize in the future.
Small RNA libraries prepared from the NLHS-treated and VP AHPND -infected L. vannamei hemocytes (NLHS-VP) generated 3 million reads across 3 libraries, which revealed 41 DEM homologs, and 27 up-regulated and 14 down-regulated miRNAs. Previously, there were 47 up-regulated and 36 down-regulated miRNA homologs identified in VP AHPND challenged L. vannamei hemocytes 10 . In this study, eight miRNA homologs that included lva-miR-184, lva-let-7, lva-miR-9, lva-miR-305, lva-miR-71, lva-miR-2, lva-miR-274, and lva-miR-317 were found to be DEMs, similar to a previous study from 2018 10 . Predicted target genes of NLHS-VP-responsive miRNA homologs such as caspase, c-type lectin, and Kazal-type serine proteinase inhibitor were similar to those of VP AHPND -responsive miRNAs identified in a previous study 10 , confirming the roles of NLHS-VP-responsive miRNAs in the regulation of immune-related genes in shrimp. The RNA-Seq data was useful in identifying genes associated with bacterial infection and NLHS response and in providing sequence data, which can be used to predict the interactions with miRNAs. Relish, a gene identified in the sequencing dataset as a NH-VP-responsive gene, was found to be highly expressed in both NH-VP and NLHS-VP groups based on RT-qPCR confirmation. This indicates that the RNA-seq dataset can correctly identify candidate genes that are of significant relevance to the experimental treatment. Nevertheless, it should be noted that the observed profiles from the RNA-seq dataset could only be detected during either heat stress or bacterial infection (but not both), as demonstrated by RT-qPCR. For example, the dynamin and lipoprotein receptor www.nature.com/scientificreports www.nature.com/scientificreports/ were predicted in the RNA-Seq dataset to be up-regulated genes in the NLHS-VP group. The RT-qPCR expression analysis confirmed this, but also demonstrated a similar expression profile in the NH-VP group. Likewise, the JHEH-1, importin-7, and DNAJ5 were predicted to be down-regulated genes in the NLHS-VP dataset and were again confirmed by RT-qPCR. However, these genes also showed the same down-regulated expression in the NH-VP group. The RNA-Seq data, therefore, needs validation by RT-qPCR to confirm the mechanisms of immune modulation that are specific for NLHS and for bacterial infection.
Gene expression data from the RNA-seq and RT-qPCR analyses also identified some genes that are up-regulated after bacterial infection when there is a prior NLHS applied to shrimp. A synergistic effect could be seen in the expression of some of these genes, e.g. lipoprotein receptor and dynamin, whose expression were not significantly changed in NH-VP, but then changed in the NLHS-VP group. This finding also supports that there is an immune modulation mechanism by NLHS prior to infection.
Within the interactome, the NLHS-VP-responsive miRNAs were mapped against DEGs. In this study, lva-miR-6813-5p, lva-miR-7278-5p, and lva-miR-745b were found to be down-regulated miRNAs that target the highest number of up-regulated NLHS-VP-responsive genes. Meanwhile, lva-miR-2898 and lva-miR-902l-5p were found to be up-regulated miRNAs that target the highest number of the down-regulated NLHS-VP-responsive genes.
In hemocyte homeostasis, caspase and transglutaminase (TGase) are two important proteins for most crustacean species 22,23 . TGase transcript expression in NLHS-treated shrimp was up-regulated at 24 hpi upon VP AHPND challenge. This expression profile is expected because TGase is known to be involved in hemolymph www.nature.com/scientificreports www.nature.com/scientificreports/ coagulation 25 . Experimental data showed that suppression of TGase results in low hemocyte counts and high bacterial counts in M. japonicas 25 , which supports that TGase is an important immune factor during bacterial infection. Up-regulation of TGase prior to bacterial infection is a good indication of how NLHS modulates the shrimp immune system in order to better survive the infection. At the same time, the down-regulation of caspase and up-regulation of inhibitor of apoptosis transcripts might synergize with the regulation of TGase by suppressing the apoptosis pathway. This suggests that NLHS-treated shrimp can better maintain hemocyte homeostasis during VP AHNPD infection.
Phenoloxidase (PO) is a key enzyme in the proPO-system that triggers the non-enzymatic conversion of phenolic substances to quinones, leading to the production of cytotoxic intermediates and melanin 26 . Also known as the melanization cascade, the prophenoloxidase system is a major innate defense system in invertebrates that involves the melanization of pathogens and damaged tissues 27 . When the system is suppressed by gene silencing in shrimp, it increases susceptibility to bacterial infection 28 . The prophenoloxidase system was identified in this study as putatively regulated by the NLHS-VP-responsive miRNAs. This pathway was examined through the expression of its canonical members, such as PO1, a putative lva-miR-3689c target. Analysis by RT-qPCR showed that NLHS treatment causes down-regulation of lva-miR-3689c and up-regulation of PO1 at 24 h post VP AHPND infection, suggesting that the proPO system is modulated by a lva-miR-3689c/PO1 interaction. This demonstrates that the lva-miR-3689c/PO1 interaction plays a crucial role in bacterial defense that leads to some form of resistance against the VP AHPND infection. Up-regulation of the two LvPO transcripts in L. vannamei (identified in a previous study) indicates an increase in disease resistance to VP AHPND 29 , which is consistent with the results presented here. The last identified pathway putatively regulated by the NLHS-VP is the antimicrobial peptide/protein (AMP) production. Shrimp AMPs are diverse and generally cationic peptides that primarily protect the host against a broad range of microorganisms. Because they also have a secondary role in modulating other immune effectors and different biological pathways, their synthesis and induction is governed by a complex interplay of various factors. Specifically, the Toll and IMD pathways are regarded as key regulating mechanisms in the transcription of AMP genes 21 . The anti-lipopolysaccharide factor (ALF) is one of major shrimp AMPs and its expression is regulated by the Toll and IMD pathways 30,31 . In this study, the transcription of ALF was found to be modulated by NLHS, along with other AMPs such as penaeidin and AMP type 2. These results consistently support AMP's importance in the immune mechanisms governing VP AHPND infection and VP AHPND toxin resistance 32,33 . According to the miRNA/mRNA network, lva-miR-3869c and lva-miR-9000 are predicted regulators of AMP type 2 and ALF AA-K, respectively, indicating the crucial roles of miRNA in modulating AMP gene expression in NLHS-treated shrimp prior to VP AHPND infection. To further characterize the regulatory mechanisms of the Toll and IMD pathways, the expression of the genes related to these pathways were analyzed. Analysis of the IMD pathway showed that transcripts of IMD, IKKε and Relish are highly expressed in the NLHS-VP shrimp, consistent with the expression of downstream genes, such as penaeidin. Previous studies have found that the high expression of IMD after V. anguillarum infection is an important regulatory mechanism in shrimp immunity against Gram-negative bacteria 34 . The IMD expression observed in the current study is consistent with the hypothesis of immune modulation by NLHS.
In conclusion, an integrated analysis of miRNAs and mRNAs can effectively elucidate the mechanisms that govern the tolerance of shrimp to VP AHPND infection after NLHS treatment. The interaction of NLHS-VP miR-NAs to their predicted mRNA targets may have a strong regulatory influence on the 3 immune-related pathways in shrimp, although these interactions will need to be directly tested. By providing new insights regarding the regulatory roles of miRNAs in the biological changes that occur in shrimp during NLHS and bacterial infection, the present study improves our understanding of the mechanisms that underlie NLHS treatment and mediate the immunity of shrimp against VP AHPND infection. Once we understand how NLHS activate the shrimp immune pathway to fight VP AHPND infection, biomolecules that can activate those immune pathways can be further identified and applied in the field. Shrimp samples. Healthy shrimps weighing 2-4 grams were obtained from a local shrimp farm and acclimatized in rearing tanks with an ambient temperature of 30 °C, water salinity of 20 parts per thousand (ppt) and constant aeration before use in experiments. The shrimps were fed with commercial pellets 4 times a day during the course of the experiments. The shrimps were sampled to determine if they were free from VP AHPND and WSSV-infections by PCR prior to experiments using the specific primer as described by Sirikharin et al. 35 and Tummamunkong et al. 36 (Fig. S2). www.nature.com/scientificreports www.nature.com/scientificreports/ After rearing at 30 °C, the chronic non-lethal heat stress (NLHS) was applied to the shrimp. The shrimps were divided into four groups of 10 shrimps each. Two groups were NLHS-treated by placing the shrimp in tanks containing 10 L of sea water at 38 °C for 5 min daily for 7 days, shrimp were given a 3-day recovery period in their respective rearing tanks. The other two groups were reared in the tank at the ambient temperature (30 °C) as control groups of non-heat (NH) treatment. Shrimp were then challenged with VP AHPND by immersion in tanks containing the bacterial ioculum at a final concentration of 1.5 × 10 6 CFU/mL. The uninfected control group was immersed in the TSB containing 1.5% NaCl. The shrimp survival was observed for 53 h. Experiment was completed in triplicates. Statistical analyses of the results were conducted using GraphPad Prism version 6. The infographic outlining all experimental groups is shown in Fig. S3.

Non-lethal heat stress (NLHS
RnA extraction. NLHS and bacterial challenge experiments were performed as described above and approximately 500 µL of hemolymph of VP AHPND -challenged NLHS-treated and VP AHPND -challenged NH control shrimp at 0, 6, and 24 h post infection (hpi) time points were drawn out from the ventral sinus using a sterile syringe pre-loaded with an equal volume of anticoagulant (27 mM sodium citrate, 336 mM NaCl, 115 mM glucose, and 9 mM EDTA, pH 5.6) 37 . Hemocytes were immediately collected by centrifugation at 800 × g for 10 min at 4 °C and kept in liquid nitrogen. The hemocytes from 30 individuals were pooled and extracted for large and small RNAs using the mirVana miRNA Isolation Kit (Ambion, Life Technologies) following the manufacturer's protocol. These experiments were done using triplicates. RNA integrity was evaluated using the Agilent 2100 Bioanalyzer chip, RNA 6000 Pico Kit and Small RNA Kit (Agilent) for large and small RNA preparations, respectively. The RNA concentrations were determined using the Qubit RNA HS Assay Kit on the Qubit 2.0 flourometer (ThermoFisher Scientific).
RnA sequencing (RnA-Seq) and data analysis. Six cDNA libraries that included 0 NLHS-VP, 6 NLHS-VP, 24 NLHS-VP, 0 NH-VP, 6 NH-VP, and 24 NH-VP with three biological replicates were prepared from 4 μg total RNA following the manufacturer's protocol for TruSeq Stranded mRNA LT Sample Prep Kit (Illumina). Eighteen indexed libraries were normalized, pooled and then sequenced with a 1% PhiX spike-in control using the NextSeq 500 Mid Output v2 Sequencing Kit (Illumina) in a NextSeq 500 desktop sequencer (Illumina). Additional adapter trimming and quality control of raw reads was performed using the FastQ Toolkit available through the BaseSpace (Illumina) public app repository. High quality reads were assembled to form a reference assembly in Trinity v2.06 software 38 . Transcript abundance was estimated using RSEM wrapped by scripts included in Trinity. Differentially expressed genes (DEGs) were detected using the edgeR software 39 and R program 40 for each treatment group and checked for sequence quality and correlation (Fig. S4). The DEGs were selected based on FDR < 0.05 and fold change >2. Pairwise comparisons between relevant groups were analyzed using 3 biological replicates and expressed as the fold change against a specific group. Bootstraps and permutation resampling were set to a default (500) as well as the largest genes size (5000). No correction was applied to FDR, but only FDRs lower than 0.05 were considered significant and listed in Table S4. The de novo assembled sequences to the L. vannamei transcripts predicted from L. vannamei genome data available in GenBank (https:// www.ncbi.nlm.nih.gov/genome/?term = vannamei).
Small RnA-Seq and data analysis. The cDNA libraries from small RNA from VP AHPND -infected NLHS-treated shrimp hemocytes at 0, 6, and 24 hpi were constructed following the manufacturer's instruction and the TruSeq Small RNA Library Preparation Kit (Illumina). Three indexed libraries were normalized, pooled, and sequenced with a PhiX control spiked at 7.5% using MiSeq Reagent Kits v2 (Illumina) in a MiSeq sequencer (Illumina). The 5′-, 3′-adapter trimming and quality control of raw reads were performed using tools in a Galaxy instance (https://usegalaxy.org/) 58 . High quality small RNA sequences with lengths shorter than 18 nucleotides, and longer than 24 nucleotides, were removed. Homology search for contaminating RNA, such as mRNA, rRNA, and tRNA was conducted using BLASTn against the NCBI nucleotide and Rfam database. After discarding the contaminating RNA, the remaining sequences were searched against miRBase 22.1 (http://www.mirbase.org/) in order to identify known miRNA homologs. Based on the number of reads from 3 libraries cut off >10, the miRNA homologs were selected for the differentially expressed miRNA (DEM) analysis. The specific procedures were as follows: (1) treatment and control groups were normalized to the same orders of magnitude. Formula: Normalized expression level = miRNA expression level/total expression level of the sample × normalized magnitude; (2) Normalized results were used to calculate the fold change and P-value 61 . Quantitative real-time pcR analysis. Several transcripts from the reference assembly were selected for quantitative real-time PCR analysis (RT-qPCR) to evaluate and confirm the differential expression profiles reported by RNA-seq analysis. The gene specific primers (Table S7) were designed by Primer3 as packaged in miRnA target prediction. The miRNA targets were identified by comparing the miRNA sequences with transcriptome data using CU-Mir software developed by our research group (http://shrimp-irn.org/mirtarget/ index.php) 63 . The software searched for the sequences on mRNA that match (perfectly) or mismatch (by one nucleotide) the seed sequences (2-8 nucleotides from the 5′-end) of miRNA. The percent complementary was calculated from the number of nucleotides that perfectly match the target mRNAs per total length of miRNA sequences. The percent total length complementary cutoff was set at 55%. The RNAhybrid (http://bibiserv.techfak.uni-bielefeld.de/rnahybrid/) was also used to predict genes targeted by the miRNAs with the parameters of free energy <−15.0 kcal/mol 65 . miRnA/mRnA interaction network analysis. In order to define all possible miRNA-mRNA interactions involved in a specific dataset of immune-related genes, only NLHS-VP-responsive DEGs and DEMs with negative correlations were grouped. NLHS-VP-responsive DEMs were used as queries to search for mRNA targets from NLHS-VP-responsive DEGs. Again, the miRNA/mRNA binding were predicted using RNAhybrid with the parameters of free energy <−15 kcal/mol. These target mRNAs were mapped against RNA-Seq data to determine their gene functions. Subsequently, the miRNAs/mRNA pairs involved in a specific dataset of immune-related genes were included in the integrated analysis of the NLHS-VP miRNA/mRNA network. The workflow of data integration approach used to build the shrimp NLHS-VP miRNA/mRNA network is shown in Fig. S5.