Transcriptomic analysis of Pacific white shrimp (Litopenaeus vannamei, Boone 1931) in response to acute hepatopancreatic necrosis disease caused by Vibrio parahaemolyticus

Acute hepatopancreatic necrosis disease (AHPND), caused by marine bacteria Vibrio Parahaemolyticus, is a huge problem in shrimp farms. The V. parahaemolyticus infecting material is contained in a plasmid which encodes for the lethal toxins PirABVp, whose primary target tissue is the hepatopancreas, causing sloughing of epithelial cells, necrosis, and massive hemocyte infiltration. To get a better understanding of the hepatopancreas response during AHPND, juvenile shrimp Litopenaeus vannamei were infected by immersion with V. parahaemolyticus. We performed transcriptomic mRNA sequencing of infected shrimp hepatopancreas, at 24 hours post-infection, to identify novel differentially expressed genes a total of 174,098 transcripts were examined of which 915 transcripts were found differentially expressed after comparative transcriptomic analysis: 442 up-regulated and 473 down-regulated transcripts. Gene Ontology term enrichment analysis for up-regulated transcripts includes metabolic process, regulation of programmed cell death, carbohydrate metabolic process, and biological adhesion, whereas for down-regulated transcripts include, microtubule-based process, cell activation, and chitin metabolic process. The analysis of protein- protein network between up and down-regulated genes indicates that the first gene interactions are connected to oxidation-processes and sarcomere organization. Additionally, protein-protein networks analysis identified 20-top highly connected hub nodes. Based on their immunological or metabolic function, ten candidate transcripts were selected to measure their mRNA relative expression levels in AHPND infected shrimp hepatopancreas by RT-qPCR. Our results indicate a close connection between the immune and metabolism systems during AHPND infection. Our RNA-Seq and RT-qPCR data provide the possible immunological and physiological scenario as well as the molecular pathways that take place in the shrimp hepatopancreas in response to an infectious disease.

Introduction L. vannamei hepatopancreas during AHPND infection, to gain a better understanding of AHPND shrimp pathogenesis and to find new genomic approaches to address this aquaculture problem.

Ethics statement
All shrimp used for experimental work was handled in accordance with the Official Mexican Standard protocols (NOM-062- ZOO-1999). Post larval Litopenaeus vannamei specimens were obtained from Laboratorio de Larvas Granmar S.A. de C.V. and transported to aquarium tanks at CIBNOR with all permits issued by the federal agency CONAPESCA. Juvenile rearing was done following the hatchery's standard procedures until the shrimps weighed 12 g. L. vannamei is not considered as an endangered or protected species.

Biological material
Healthy L. vannamei weighing about 12.26 ± 0.022 g were maintained in 1000 L PVC tanks filled with filtered (1 μm) and aerated seawater, salinity of 35 PSU, at 25˚C. Shrimp were acclimatized and fed ad libitum daily for one week, with commercial feed (Camaronina, Purina, Mexico City). V. parahaemolyticus IPNGS16 was isolated and characterized from dying shrimp obtained from farms in Guasave, Sinaloa, Mexico [16].
PCR cycling conditions were as follows: initial denaturalization 94˚C 2 min, 35 cycles of 95 C 30 s, 60˚C 25 s and 72˚C 30 s, and a final extension of 72˚C 5 min. PCR products were separated in a 1.2% agarose gel (1X TAE) and stained with Red Gel 1X (Biotium, Fremont, CA). PCR products were sequenced in both directions by Sanger method [19]. The molecular weight of the PirB and PHP region amplified is 501 bp and 272 bp (S1 Fig) of plasmid pVPA3-1 (KM067908.1).

Median lethal concentration (LC 50 )
Shrimps were transferred to 40 L plastic tanks containing 20 L of filtered (1 μm) and aerated seawater, and salinity of 35 PSU. Ten organisms were placed in each tank and acclimatized for three days at 28˚C.
To calculate the median lethal concentration (LC 50 ), five bacterial cell concentrations were prepared; each concentration was tested in triplicate.
Briefly, V. parahaemolyticus stored in 30% glycerol at -80˚C, was reactivated in TSB+ medium and incubated for 16 h at 28˚C. Bacterial culture was then centrifuged at 2,600 × g for 20 min at 4˚C; bacteria pellets were washed and resuspended in saline solution (2.5% NaCl). The bacterial solution was adjusted spectrophotometrically to an optical density of 1 at 600 nm (Bio-Rad, Smart-Spec 3000) [16].
Different concentrations of bacteria suspensions (350, 750, 1000, 1500, and 3000 CFU/mL) were applied to the shrimp culture tanks (20 L total volume) and incubated for 24 h. Shrimp mortality was recorded every 6 h by visual inspection in each tank, dead shrimps were removed from the tanks. Shrimp was grown in optimal conditions, however, no cleaning of the tanks was made during the challenging period and temperature was maintained at 28˚C to promote V. parahaemolyticus infection. LC 50 was determined at 24 hours using Probit analysis in IBM SPSS v23 (IBM, Armonk, NY) [20].

Experimental design
Before bacterial challenge, experimental organisms were examined under the microscope to determinate molt stage based on the degree of setae development in uropods [21].
Organisms in the intermolt stage (C1) were placed in 20 L plastic tanks at 28˚C and 35 PSU under continuous aeration and left to acclimatize for three days. A total of 18 tanks (9 control and 9 experimental) were set up with ten organisms per tank. Experimental infection was carried out using a bacterial concentration of 660 CFU/mL (LC 50 at 24 h) per tank; experimental organisms were exposed to the bacteria for at least 72 h.

Hepatopancreas RNA-seq
Total RNA from 24 h samples (control and infected groups) was extracted with Tri-reagent (Sigma-Aldrich, St. Louis, MO) according to manufacturer's protocol. RNA samples were quantified using Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA), the RNA integrity was checked with 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) only samples with an RNA integrity number (RIN) > 8.5 were used. Samples for each treatment (control and infected) were pooled into one tube each and sent for sequencing to GENEWIZ, LLC. (South Plainfield, NJ). mRNA was isolated from total RNA using Oligo d(T) beads and fragmented for 15 min at 94˚C. First and second cDNA strands were synthesized and the 3'-ends were repaired and adenylated. Subsequently, sequencing adapters were ligated to cDNA fragments to prepare a cDNA library.
RNA sequencing library was prepared using the NEBNext Ultra RNA Library Prep Kit (New England Biolabs, Ipswich, MA) from Illumina following manufacturer's recommendations. The universal adapter was ligated to cDNA fragments, followed by index addition and library enrichment with limited PCR cycle. Sequencing libraries were validated using a DNA Chip on the Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA), and quantified using Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA) as well as by quantitative PCR (Invitrogen, Carlsbad, CA) as well as by quantitative PCR (Applied Biosystems, Carlsbad, CA, USA).
The sequencing libraries were then multiplexed and clustered onto a flow cell; and loaded in the illumina HiSeq 2500 instrument (San Diego, CA) according to manufacturer's instructions. Samples were sequenced using a 2X100bp paired-end configuration. The HiSeq Control Software conducted the image analysis and base calling. Raw sequence data generated from Illumina HiSeq 2500 was converted into fastq files and de-multiplexed using illumina bcl2fastq v1.8.4.

Transcript quantification and differential expression analysis
Transcript abundance was estimated using Trinity v2.4.0 downstream analysis pipeline [25]. Transcript abundance was estimated with Bowtie2 v2.3.4 and RSEM using both libraries (control and infected) [26,37]. The transcript expression matrices were built in a database, and the expression patterns were calculated using the Transcript per Kilobase Million (TPM) method [38,39]. The read counting was normalized to obtained relative expression. We used the False Discovery Rate (FDR) method to ensure the statistical significance of differentially expressed transcripts (DETs). Transcript clusters of our DETs of interest were generated using FDR�0.01 and log2Ratio�2 [40]. The resulting DETs clusters were annotated as mentioned above.

GO term and KEGG pathway enrichment analysis of DETs
To select genes of interest from our DETs, a Gene Ontology (GO) enrichment analysis was conducted with 'GOseq' based on the hypergeometric test (P � 0.01) value [41]. Additionally, a pathways enrichment analysis was performed using the KEGG Automatic Annotation Server (KAAS) [42].

Network analysis and protein-protein interaction of DETs
Functional interactions of DETs were predicted using the GeneMANIA algorithm, a tool from the Cytoscape v3.6.0 software [43,44]. The parameters used were based on Gene Ontology (GO) under the term "biological process," and Drosophila melanogaster was used as a source species. Three different networks were constructed depending on their gene expression levels: up-regulated transcripts, down-regulated transcripts and a combination both. The networks were constructed based on different interaction levels: co-expression of genes, physical and genetic interactions, pathway participation, protein co-localization, protein domain similarity, and predicted protein interactions.

Identification of hub genes and DETs clustering
The hub genes are represented as nodes with larger interactions into the network and were identified by calculating the node degree distribution [45,46]. The network that combined the up and down-regulated transcripts was used to perform a community analysis using cluster-Maker and the greedy algorithm [47][48][49]. The resulting clusters were analyzed by enrichment analysis using the DAVID tool [50].

Gene selection and primer design
We selected ten candidate transcripts based on our differential expression analysis, these genes were selected because they are known to participate in the immune and metabolic response against AHPND. Primers for qPCR were designed using the primer3plus software, the list of primers are shown in Table 1. Primers at [25 nM] were synthesized at T4OLIGO (Irapuato, Guanajuato, Mexico) and purified by desalting method (DST) [51].

Relative gene expression analysis
Total RNA extraction from hepatopancreas was conducted as mentioned above. RNA quantity and purity were assessed measuring the A 260nm /A 280nm and A 260nm /A 230nm ratio respectively using Nanodrop 1000 (Thermo Scientific, Chicago, IL). Additionally, RNA quality and 18S and 28S rRNA band integrity were verified by electrophoresis on a 1.2% Bleach gel [1X TAE] [52]. Afterward, 1 μg of total RNA was treated with 1 u of DNAse I amplification grade (Sigma Aldrich, St. Louis, MO). For first-strand synthesis we used ImPromp II reverse transcription system (Promega, Madison, WI) following manufacturer's instructions using 0.5 μg of oligo d (T) for mRNA capture, cDNA was stored at -80˚C until use.
Afterward, 1 μg of total RNA was treated with 1 μL of DNAse I (1 u/μL) amplification grade (Sigma Aldrich, St. Louis, MO). For first-strand synthesis we used ImPromp II reverse transcription system (Promega, Madison, WI) following manufacturer's instructions using 0.5 μg of oligo d(T) for mRNA capture, cDNA was stored at -80˚C until use.
Amplification efficiencies (E) for all primers used in this experiment were conducted by slope calculation of 4-fold serial dilutions starting with 150 ng/μL of cDNA on a threshold value of 0.034 and calculated automatically using the Rotor-Gene 6000 v1.7 software. Based on the 2 -ΔΔCq calculation, the qPCR data was analyzed using a one-way ANOVA; significant statistical differences were obtained with Holm-Sidak multiple comparisons test (α = 0.05). Final statistical analyses were performed in SigmaPlot v11 (Systat Software Inc., Germany). Statistically significant differences were set at p < 0.05 [53].

Characterization of V. parahaemolyticus IPNGS16 strain
To confirm that we were using the correct bacterial strain for our infection experiments, an end-point PCR was performed to screen for the PHP and PirB gene sequences from two different V. parahaemolyticus strains: the IPNGS16 strain which tested positive for the two gene sequences and the CAIM170, a non-AHPND causing strain which is negative for PHP and PirB sequences (S1 Fig).

Experimental infection with V. parahaemolyticus IPNGS16 and lethal concentration 50 (LC 50 )
To determine the LC 50 , first, we tested a range of lethal and sub-lethal bacteria concentrations going from 3,000 to 350 CFU/mL. We observed that using 3,000 CFU/mL, 60% of the  Fig 1A). Using the values obtained in this experiment, a Probit analysis was run to calculate the LC 50 . The LC 50 value obtained was equivalent to 660.95 CFU/mL at 24 hpi ( Fig 1B).

Histopathology of hepatopancreas
Early stages of AHPND in infected shrimps are characterized by shedding of epithelial cells in the hepatopancreas, followed by necrosis of epithelial cells, and hemocytic infiltration at later stages [6]. To confirm that our specimens were indeed developing AHPND due to V. parahaemolyticus infection, we selected 5 infected and 5 control organisms at 0, 12, 24, 48, and 72 hours to perform histological sections, and hematoxylin and eosin staining. The hepatopancreas epithelium is distinguished for having well defined and clear tubule structures (Tub), a star shape lumen (Lum), and B-cells (HpB) (Fig 2A).
Our observations indicate that by 12 hpi, there is a sloughing (Slo) of epithelial cells from the tubule (Fig 2B). At 24 hpi, epithelial tissue shows a clear presence of pyknotic cells, a characteristic feature for necrosis (Nec) (Fig 2C). By 48 hpi, the tubule lumen (Lum) has lost its distinctive star-shape ( Fig 2D) and, by 72 hpi, structural atrophy (Atr) and epithelial cell detachment is seen in B, R, and F cells, which are structural components of the hepatopancreas tubule (Fig 2E), and granular hemocytes (HemG) and hyaline hemocytes (HemH) were observed (Fig 2F).

De novo assembly and functional annotation
Two cDNA libraries were constructed from non-infected and infected shrimp hepatopancreas, the total raw data obtained from both libraries were 166,174,430 paired-end reads. Transcriptomic analysis of Litopenaeus vannamei in response to AHPND caused by V. parahaemolyticus After sequence trimming, a total of 148,571,949 paired-end reads were retrieved from both libraries. Trimmed paired-end reads obtained from the raw data were deposited in the Short Read Archive (SRA) database of The Nacional Center for Biotechnology Information (NCBI) with the following SRA accession codes SRR7986779 (non-infected), and SRR7986780 (infected). According to MIGS standards, the summary of this project is shown in S1 Table [54].
After the de Novo assembly was completed, we found 174,098 transcripts with an N 50 value of 1,538 bp, and an average transcript mean size of 984.07 bp (Table 2), of which 86,388 transcripts (49.6%) were 450bp in length, followed by 34,272 transcripts (19.8%) around 900 bp and 12,375 transcripts were over 3 Kb (Fig 3B). All transcripts obtained were annotated with BLASTx using SwissProt (96,988 hits) and UniProtKB INV (199,139 hits). For BLASTp we only used SwissProt (41,807 hits). Sequence identity percentage and e-value distribution histograms are shown in Fig 3C and 3D, respectively.
Gene ontology (GO) annotations were assigned to 47,886 transcripts (Table 2), the Swis-sProt database was only used to assign annotations. The three main ontologies (Fig 4A) represented are: 41.0% biological process (BP), 32.86% cellular components (CC), and 26.24% molecular functions (MF). The COG database was used to generate orthologous protein groups from the assembled transcripts. Functional categories were identified (Fig 4B), being the translational ribosomal structure and biogenesis categories with the highest number of associated transcripts (14.75%), followed by general function prediction only (12.54%), amino acid transport and metabolism (7.05%), energy production and conversion (6.64%), and carbohydrate transport and metabolism (5.74%). The transcripts obtained were annotated using the KEGG database to determinate the relationship between transcript function and their specific metabolic pathway (Fig 4C). A total of 8,083 transcripts were classified in metabolic pathways; within this category, 1,191 transcripts were in carbohydrate metabolism and 1,279 transcripts related to the family of enzymes ( Fig 4C). Interestingly, the category that contains the highest number of transcripts (12,943) is the processing of genetic information, which includes genes for zinc finger protein, transcription factors, and basic loop-helix, a second group of 4,049 transcripts were classified into the folding, sorting and protein (Fig 4C). Transcriptomic analysis of Litopenaeus vannamei in response to AHPND caused by V. parahaemolyticus

Determination of differentially expressed transcripts (DETs) in hepatopancreas
Based on the FDR�0.01 and Log2Ratio�2 threshold, a total of 915 transcripts were identified as DETs, from which two clusters were obtained, one containing 442 up-regulated transcripts and a second cluster with 473 down-regulated transcripts (Fig 5). Both clusters were annotated with BLASTx and BLASTp using the SwissProt and UniRef90 databases, in S2 and S3 Tables show a list of the up and down-regulated transcripts. To further analyze the functional role of these transcripts, an enrichment analysis was performed to obtain enriched and depleted GO terms in the up and down-regulated transcripts and grouped into three main categories: biological process, molecular function and cellular components (Fig 6A and 6B). Up-regulated transcripts show the following enriched terms for biological process: metabolic process (7,015  Transcriptomic analysis of Litopenaeus vannamei in response to AHPND caused by V. parahaemolyticus

Differentially expressed transcripts at protein-protein analysis
The DETs found in the assembled hepatopancreas transcriptome were analyzed in GeneMA-NIA using Drosophila melanogaster as source species to predict protein-protein interactions between up and down-regulated transcripts. We found 11 protein-protein interactions predicted by the GeneMANIA algorithm. Based on their biological functions, we constructed three regulatory networks, the first network (Network 1) represents up-regulated genes (S2A Fig), it was constructed with 203 predicted proteins, the main enriched terms are chitin metabolic process, which represents 3.94% of the proteins in the network, homophilic cell adhesion via plasma membrane adhesion molecules (2.46%), and oxidation-reduction process (7.88%) see in S4 Table. Network 2 represents down-regulated transcripts (S2B Fig), 230 predicted proteins were used to construct this network, the enriched terms are axonogenesis (3.083%), sarcomere organization (2.20%), and oxidation-reduction process (7.48%) see in S5 Table. Lastly, we constructed a third network combining the up and down-regulated transcripts (Network 3), the enriched terms associated to this network are oxidation-reduction process (7.63%), sarcomere organization (1.84%), and imaginal disc-derived wing morphogenesis (4.73%) (S6 Table).

Identification of hub transcripts
One of our most significant findings come from Network 3, it contains 20 hub genes that represent the maximum number of interactions within any of the three networks ( Table 3). One of the main hub genes identified in Network 3 was thrombospondin, the GO terms assigned to this gene are: wound repair, cell adhesion and immune system response. Interestingly, when thrombospondin is up-regulated, the enriched GO terms are associated with anatomical structure morphogenesis, localization, multi-organism process, cell development, and stress response. Whereas down-regulated thrombospondin enriched GO terms are associated with intercellular transport, male germ-line stem cell population maintenance, protein N-linked glycosylation, neuron projection morphogenesis, and positive regulation of JAK-STAT cascade (Fig 7).

Gene selection, relative gene expression analysis and transcriptome validation
In order to validate candidate transcripts, we selected ten transcripts based on the following criteria: the GO enrichment analysis data, the fold-change value obtained for each transcript and their known function (S4 and S5 Tables). The main enriched GO terms were in oxidation- reduction processes, carbohydrate derivative metabolic processes, apoptosis regulation and biological adhesion (Fig 6A), the selected genes are shown in Table 1 [55]. The RNA-seq expression results were validated by RT-qPCR using the original extracted RNA from both control and infected hepatopancreas samples from different time points after infection (S3 Fig). Based on their relative expression, DETs were grouped into two categories, those related to the immune system (Fig 8) and metabolic processes (Fig 9). LvSOD (Fig 8A) increased its expression at 6, 24, 48 hpi in the infected group compare to controls. However, a significant decreased in expression is observed at 12 hpi. LvPRDX showed a gradual increased expression from 12 to 48 h in the infected group (Fig 8B). Relative expression of the LvCLO gene, which encodes for a coagulation protein, increased significantly only at 3 and 24 hpi (Fig 8C). The LvHYC gene that encodes for hemocyanin increased its expression by 150-fold at 24 hpi and decreases by 48 h, although there are not significant statistical differences between the two time points (Fig 8D). Expression of LvTHBS that encodes for thrombospondin had a slight increase in its expression at 3, and 48 hpi but, a decrease in expression at 6, and 24 hpi ( Fig  8E). Lastly, the cellular apoptosis susceptibility gene LvCAS, showed increased expression at 3, 12, 24 and 48 hpi when compare to the control group, reaching its highest expression (15 fold) at 24 hpi (Fig 8F).
The relative expression of genes connected to metabolism is shown in Fig 9. The phosphoenolpyruvate carboxykinase (LvPEPCK) gene starts to increase its expression by 6 and 12 hpi and by 24 hpi its expression increased 20-fold in comparison to the control group (Fig 9A). The chymotrypsin gene (LvCTR) has small increments in its expression at 3, 6 and 24hpi; however, the only statistically significant difference in increase expression is at 24 hpi when compared to the control group (Fig 9B). Gene expression of hypoxia-up-regulated gene (LvHYO) increased at 6 and 24 hpi; its expression was decreased at 12 and 48 hpi (Fig 9C). Finally, the LvCTSH expression, a gene which encodes for the Cathepsin L protein, showed a significant decrease at 24 hpi (Fig 9D), the rest of the time points were not statistically significantly difference compared to the controls.

Discussion
Acute hepatopancreatic necrosis disease (AHPND) is an emergent disease caused mainly by V. parahaemolyticus, which induces 90-100% mortality in shrimp cultures worldwide [56]. In this study, we have identified a significant number of DETs and started to understand the molecular pathways that take place during AHPND. Before the experimental infection trials, the V. parahaemolyticus strain was characterized by PCR assay by amplifying the PirB and PHP genes. A LC 50 was standardized and set at 660 CFU/mL, this concentration represents 50% mortality at 24 hpi, and these data were used for our subsequent RNA-seq experiment. Juvenile L. vannamei were infected with V. parahaemolyticus IPNGS16 (AHPND positive) by immersion [57]. Following histopathological results at 12, 24, 48 and 72 hpi a sequential progression of hepatopancreas lesions was observed in AHPND infected tissue. Necrosis, epithelial cells sloughing, atrophy, and massive hemocyte infiltration signs were more evident from 48 hpi onwards (Fig 2D), our histopathological results agree with previous reports [4,6,57].
In our differential expression analysis, a total of 915 transcripts were detected between the infected and non-infected control groups. The number of up-regulated genes was considerable more (473) than the down-regulated genes (442); probably due to biological changes and inhibition of basal cellular functions caused by V. parahaemolyticus infection [58]. Our knowledge of the AHPND pathogenesis indicates that V. parahaemolyticus produces toxins that cause cell sloughing [59], interestingly we found two immune-and metabolic-related gene clusters in our KEGG analysis. The immune-related gene cluster contains antimicrobial peptides (AMPs), oxidative stress proteinases and inhibitory proteinases, pattern recognition, coagulation, and signaling pathway proteins, along with other immune-related proteins (S1 and S2 Tables). In the metabolic-related cluster, we found genes related to carbohydrate metabolism, beta-oxidation, lipid metabolic processes, and fatty acid biosynthesis (S1 and S2 Tables). This differential transcript analysis provides us with the basic information about the genes, their expression and the immune and metabolic processes that take place during AHPND [60].
The L. vannamei innate immune system contains cellular and humoral components that work individually and synergistically to protect the organism integrity during infection. Activation of the shrimp immune response depends on diverse mechanisms that include the antioxidant system, hemolymph, wound repair and apoptosis-related proteins [60].
All living organisms produce the reactive oxygen species (ROS) during normal aerobic metabolism. Stress conditions like pH, temperature, hypoxia and microbial infections can shortage oxygen and increase ROS, resulting in oxidative stress within the cells [61]. Increased levels of ROS causes macromolecules damage which in turn affect membranes and enzymes structure and nucleic acids integrity [62]. The endogenous cellular defense mechanisms to stabilize increased levels of ROS include antioxidant molecules and enzymes (catalase, glutathione peroxidase, and superoxide dismutase) [63]. Our GO enriched analysis (P < 0.01) show 442 up-regulated transcripts which were grouped in different categories including oxidationreduction processes, biological process (BP) and oxidoreductase activity for molecular function (MF).
In this work, we identify two oxidative stress-related transcripts, peroxiredoxin (up-regulated) and superoxide dismutase (down-regulated). Peroxiredoxin belongs to a group of antioxidant enzymes that protect the cell from ROS by reduction of a wide range of cellular peroxides [64]. Our gene expression analysis results show that LvPRDX expression increases significantly after 12 hpi (Fig 8B). This result agrees with a previous report where peroxiredoxin 5 expression levels increase during V. angillarum infection in E. carinacauda [65]. We also found that the antioxidant superoxide dismutase (LvSOD) gene was down-regulated in our RNA-seq data (S3 Fig), however in our qPCR results LvSOD expression was increased at 6, 24, and 48 hpi (Fig 8A). Although we are not certain of why LvSOD shows a complete difference in expression, we know that V. parahaemolyticus and white spot syndrome virus (WSSV) infections increase mRNA levels of these genes, and ROS play an essential role in L. vannamei defense against these pathogens [66].
Hemocyanin, an abundant protein present in the hemolymph of arthropods and mollusk [67]. Its primarily function is to transport and store molecular oxygen, along with other multiple physiological functions like protein storage, osmoregulation, molt cycle (post-molt (A) and pre-molt (D)), exoskeleton formations, and antimicrobial activity as well as a non-specific immune molecule [68,69]. Our RNA-seq data shows that multiple isoforms of LvHCY are present in both up and down-regulated transcript sets. Gene expression analysis of LvHCY shows a significant increase expression after 24 hpi compare to the control. Previous studies have shown that Hemocyanin is overexpressed during AHPND and PirAB challenge. Furthermore, Boonchuen and colleagues demonstrated the toxin-neutralizing activity of hemocyanin against V. parahaemolyticus infection [70].
Among of the first immune responses that have the organisms to either fight against pathogens or to prevent loss of hemolymph during tissue damage is clotting.
Our results indicate that hemolymph cottable protein (LvCLO) expression increases significantly at 3 hpi studies have shown that disrupting the ProPO system, a coagulation-related cascade, using miRNA silencing increases bacterial load in the host [71]. Moreover, genomic studies in P. monodon showed up-regulation of coagulation-related genes during stomach infection caused by V. parahaemolyticus [12]. These findings show activation of the coagulation system is a fundamental step during bacterial infection [72].
Thrombospondins (THBSs) are a family of extracellular calcium-binding glycoproteins, in vertebrates THBSs are involved in cell-extracellular matrix interactions, angiogenesis, synaptogenesis and organization of connective tissue [73]. However, very little is known about the function of THBSs in invertebrates, for example in penaeid shrimp, it has been suggested that THBSs are part of the defense response against microbial infection, in Fennerpenaeus chinenesis challenged with V. angillarum and Staphylococcus aureus, TSP is up-regulated in hemocytes, heart, intestine, stomach, and hepatopancreas [74]. Moreover, in P. monodon THBSs are implicated in the immune function of the lymphoid organ against two shrimp pathogens, V. harveyi and the white spot syndrome virus (WSSV) [75] and, RNA-Seq analysis in P. monodon stomach showed up-regulation of three THBSs during AHPND [12]. Our RNA-seq data show up-regulation of thrombospondin, and protein-protein networks show THBS as part of the top-20 of the hub genes in the network 3. Furthermore, qPCR gene expression analysis shows an increased expression of LvTHBS at 3 and 48 hpi suggesting that LvTHBS may be involved in the defense mechanism during AHPND.
One of the most frequent proteins up-regulated during bacterial infections is the cellular apoptosis susceptibility protein (CAS), CAS is a nuclear transport factor associated tumor necrosis factor (TNF)-mediated apoptosis, studies have suggested that in healthy cells, CAS acts as a switch of proliferation or apoptosis [76], and it has been linked to invasion and metastasis of cancer cells [77]. In other animals, CAS is up-regulated in non-specific cytotoxic cells in response to apoptosis [78]. Crustacean CAS in F. chinenesis has been related to ovary development and suggested to function as a nuclear protein. Moreover, CAS upregulation has been shown in Erichor sinensis hemocytes after Spiroplasma eriocheris challenge [79,80]. Furthermore, hemocytes response during AHPND has shown activation of the apoptosis gene CAS2 in L. vannamei [81]. Our RNA-seq and gene expression analysis shows up-regulation of LvCAS, suggesting that the apoptosis process is an essential mechanism response to AHPND.
In addition to the immune-related cluster, we also looked for metabolism-related DETs in our RNA-Seq data. A connection between the immune system and metabolism (immunometabolism) has not been demonstrated in vertebrates yet; this is mainly due to the highly specialized organ structures and their complex functions in both systems. However, invertebrates possess hepatopancreas, the evolutionary forerunner of liver and pancreas, a much more simpler organ that integrates both, metabolic and immune functions, which is an ideal model to study immunometabolism [82].
One of the main enzymes involved in gluconeogenesis is the phosphoenolpyruvate carboxykinase (PEPCK), in hypoxic conditions, PEPCK expression increases to keep glucose homeostasis within the cells [83]. Vibrio infections can decrease oxygen uptake during pathogenesis and induce hypoxia [84]. Our RNA-seq results and gene expression analysis show up-regulation of PEPCK expression by 17-fold at 24 hpi. The LvPEPCK overexpression is probably related to 1) the hypoxia conditions and lack of oxygen uptake due to the Vibrio infection and 2) the difficulty in maintaining glucose levels, during AHPND shrimp cannot feed and use other intermediaries metabolites to biosynthesize de Novo glucose [85].
Cellular oxygen deprivation alters metabolic pathways resulting in a hypoxia response, to prevent further damage due to low oxygen levels, cells protect themselves with a group of chaperone proteins called oxygen-regulated proteins [86]. Interestingly, the hypoxia gene LvHYO was up-regulated at 6 and 24 h in the infected group (Fig 9C). Over-expression of hypoxia genes in bacterial pathogenesis has been related to an increase in bacterial load during early AHPND stages causing hepatopancreatic septicemia, probably due to the high levels of PirAB toxins in the cells [5,11].
Chymotrypsin (CTR) is a serine-protease enzyme that participates in immune reactions and other physiological functions [87]. The synthesis of CTR is located in hepatopancreas tubules, and mRNA levels of CTR increase after infection with V. anguillarum in F. chinensis. Interestingly, our RNA-seq data and qPCR analyzes also show increased LvCTR transcript levels at 24 hpi suggesting that LvCTR expression is important not only to maintain the normal cellular functions but also in the immune response against Vibrio infections, although additional tests will be needed to determine the exact function of CTR during AHPND in shrimp [88].
Cathepsins (CTSHs) are proteases that maintain cellular turnover, acting as scavenger proteins and food digestive. Cathepsin L in shrimp is associated with rapid cell differentiation of F and B-cells in hepatopancreas [89,90]. Our data show down-regulation of LvCTSH at 24 hpi; however, there were no significant differences in expression in any of the time point tested, suggesting that LvCTSH might not play an important role during AHPND. However, as mention above, we cannot rule out the importance of LvCTSH expression in keeping the cell differentiation in the tubule structure especially since our histopathological analysis shows epithelial cells sloughing at 12 hpi.
As mentioned, several studies have looked at the pathogen-host response in shrimp using either different shrimp species or pathogens [5,12,15,75,81]. However, very few studies have investigated the shrimp immunological response at the transcriptome level. For example, Ge et al. (2017) which used the ridgetail white prawn (E. carinicauda) and V. parahaemolyticus to challenge the specimens by intramuscular injection found up-regulation of the hypoxia upregulated gene (HCY) and chymotrypsin (CTR), and down-regulation of the Cathepsin L protein (CTSH) gene [15], these data agree with our results, we also found that these three genes have the same gene expression dynamics after bacteria challenge. Furthermore, we both found very similar enriched GO terms including, antimicrobial peptides, lectins, serine-protease cascade, oxidative-stress, pattern recognition proteins, coagulation, and heat-shock proteins (S2 and S3 Tables). A second study performed by Soonthornchai and colleagues (2016) showed that peroxiredoxin (PRDX), the cottable protein (CLO) gene and thrombospondin (THBS) were up-regulated in the black tiger shrimp (P. monodon) stomach after V. parahaemolyticus infection, which are exactly the same up-regulated genes found in our L. vannamei hepatopancreas samples after infection [12].
In addition to identifying genes involved in the host immune response, like in both previous studies, we also looked for genes that participate in metabolic pathways. Bacterial infections tend to modify their host metabolism via homeostasis disruption [91]. Based upon our analysis, we identified a new set of differentially expressed metabolic genes, among them: the antioxidant superoxide dismutase (SOD), hemocyanin (HCY), the cellular apoptosis susceptibility protein (CAS) gene, phosphoenolpyruvate carboxykinase (PEPCK), chymotrypsin (CTR), and the hypoxia up-regulated (HYO) gene. Interestingly, none of these genes have previously been reported as being differentially expressed during pathogen infection. Determining the function of these genes during the metabolic response will provide us with a deeper understanding of how the metabolic and immune system interact and work towards specific pathogens.

Conclusions
In conclusion, we have constructed a de Novo transcriptome of L. vannamei hepatopancreas infected with V. parahaemolyticus (AHPND) at 24 h using the Hi-Seq 2500 Illumina platform. We found a total of 915 differentially expressed transcripts, from which 442 were up and 473 down-regulated, up-and down-regulated transcripts were found enriched in different categories: metabolic process, oxidation-reduction process, regulation of programmed cell death, biological adhesion mainly located in the extracellular matrix. Protein-protein network analysis showed interactions between up and down-regulated transcripts within the same category, for instance in the cellular process group, response to stress genes were turned on whereas intracellular transport genes were turned off. This protein network analysis also allowed us to select transcripts with the highest number of interactions between the up-and-down-regulated clusters, a clear example is the LvTHBS gene, which due to its interaction and position in the network was defined as a hub gene.
Based on our enriched GO terms results and previous transcriptomic published data, we selected ten candidate genes associated to digestion, metabolism and the immune system to perform qPCR. Gene expression results of the candidate genes showed that each gene has a specific expression dynamics at the different times points analyzed, which agrees with our RNA-seq differential expression data. It will be necessary to perform proteomics analysis to study in depth how the metabolic and immunologic regulation works in shrimp during the AHPND since this knowledge can be the basis for designing therapeutic strategies. Our GO data, protein-protein networks and gene expression analyses strongly suggest a close interaction between the hepatopancreas metabolism and immune system, indicating that both systems respond collectively against a Vibrio parahaemolyticus infection.