Two Independently Comparative Transcriptome Analyses of Hemocytes Provide New Insights into Understanding the Disease-Resistant Characteristics of Shrimp against Vibrio Infection

Simple Summary Identification of disease-resistant genes and their application for disease-resistance breeding is an important way to solve disease problems in shrimp aquaculture. Although some transcriptomics studies have been performed to identify genes related to pathogen infection and host defense, there is still a lack of deep understanding of the host-resistant characteristics against V. parahaemolyticus infection. The present study established a method to obtain hemocytes from shrimp with clear resistant abilities against V. parahaemolyticus infection. We also proposed a strategy to identify key biological processes and genes related to disease resistance based on the analysis of the present transcriptome data. Abstract Vibrio parahaemolyticus carrying plasmid encoding toxins PirA and PirB is one of the causative agents leading to the severe disease of AHPND in shrimp aquaculture. However, there is a lack of deep understanding of the host-resistant characteristics against V. parahaemolyticus infection. Here, we established a method to obtain hemocytes from shrimp with different V. parahaemolyticus-resistant abilities and performed comparative transcriptome analysis on the expression profiles at the background level of hemocytes from shrimp in two independent populations. Principal component analysis and sample clustering results showed that samples from the same population had a closer relationship than that from shrimp with similar disease-resistant abilities. DEGs analysis revealed that the number of DEGs between two populations was much more than that between V. parahaemolyticus-resistant and susceptible shrimp. A total of 31 DEGs and 5 DEGs were identified from the comparison between V. parahaemolyticus-resistant and susceptible shrimp from populations 1 and 2, respectively. DEGs from population 1 were mainly cytoskeleton-related genes, metabolic related genes, and immune related genes. Although there was no DEGs overlap between two comparisons, DEGs from population 2 also included genes related to cytoskeleton and metabolism. The data suggest that these biological processes play important roles in disease resistance, and they could be focused by comprehensive analysis of multiple omics data. A new strategy for screening key biological processes and genes related to disease resistance was proposed based on the present study.


Introduction
Shrimp is one of the most important aquatic animals in aquaculture. However, frequent diseases caused by bacteria, viruses, and parasites always lead to severe economic losses, which greatly hinder the healthy development of shrimp aquaculture. AHPND,

Experimental Animals
Two batches of Litopenaeus vannamei purchased from Rizhao breeding farm in Shandong were used for two independent experiments. Before experiments, the shrimp were cultured in aerated seawater laboratory aquaculture tank at 26 • C for seven days. The shrimp weights were 10 ± 0.5 g in the first batch (population 1) and 19.5 ± 1 g in the second batch (population 2). The hepatopancreas from three individuals was dissected, mixed, and ground. A total of 0.18 g tissue sample was homogenized in 1 mL PBS and then spread on the TCBS solid medium plate. The plate was put in the 30 • C incubator overnight to Biology 2023, 12, 977 3 of 13 confirm whether the shrimp were free of Vibrio pathogens. Three biological replicates were performed for each population.

VIE Fluorescent Labeling and Hemolymph Collection
The VIE distinguished with eight different colors (red, pink, orange, brown, white, purple, green, yellow) was injected into each shrimp in the third segment with different color fluorescence. A total of 80 and 70 individuals were used in two experiments, respectively. Then, the shrimp were cultured in aerated seawater at 25 ± 1 • C for seven days and fed with artificial food twice a day. The seawater was changed every day. After one week, 200 µL hemolymph was exsanguinated from each shrimp with a sterile syringe containing an equal volume of ice-cold anticoagulant (27 mM trisodium citrate, 336 mM sodium chloride, 115 mM glucose, 9 mM Na 2 EDTA·2H 2 O, pH 7.4). The hemocytes were separated by centrifugation at 1000 g for 10 min and stored in −80 • C. Then, all the shrimp were returned to the laboratory tanks.

V. parahaemolyticus Immersion and Hepatopancreas Collection
One week after exsanguination, the shrimp were infected with V. parahaemolyticus AG01, a strain isolated in our lab [11]. During V. parahaemolyticus infection, shrimp were cultured in the same condition while the seawater was not changed for two days. The bacteria strain was cultured in tryptic soy broth with 2% sodium chloride liquid media. The cultured bacteria were further checked by positive PCR amplification of PirA vp and PirB vp genes. Bacterial titer was counted with a hemocytometer under a light microscope. The soak infection doses in two batches of shrimp were set at the lethal concentration 50, 2.5 × 10 6 CFU/mL and 5 × 10 6 CFU/mL, respectively. The mortality was continuously recorded for two days. The hepatopancreas of moribund shrimp and survival shrimp after 48 h post infection was collected and stored at −80 • C. In addition, the hepatopancreas of moribund shrimp was also collected for pathogen detection. The tissue was treated as described in Section 2.1, and the Vibrio strain grown on the TCBS plate was collected as PCR template for 16S rDNA amplification and Sanger sequencing. The Vibrio strain was confirmed as V. parahaemolyticus before further analysis.

DNA Extraction and Bacteria Load Detection
The DNA from the hepatopancreas was extracted using TIANGEN Plant Genomic DNA Extraction Kit according to the manufacturer's instruction. The concentration and purity of DNA were determined using NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA) by coagulation and 1% gel electrophoresis. PirA copy number was estimated by Ependorf MasterCycler EP Realplex (Ependorff, Germany) using TaqMan Assay. A forward primer PirA-F and reverse primer PirA-R were used to amplify a 284 bp fragment of PirA gene, with a TaqMan hydrolysis probe 5 -FAM-CCGCCAGCCATAAATGGCGCACC-BHQ1-3 ). The amplification program consisted of one cycle of contamination digestion at 37 • C for 2 min, one cycle of pre-denaturation at 95 • C for 5 min, and 45 cycles of denaturation at 95 • C for 10 s and extension at 60 • C for 30 s. A standard curve was obtained using serial dilutions of plasmid PirA (full-length ORF of PirA gene of V. parahaemolyticus was cloned into pUC57 vector), which was used to quantify the V. parahaemolyticus copy number. Each assay was carried out in four parallels. Therefore, all shrimp were divided into resistant group and susceptible group (9 hpi, 12 hpi, and the surviving individual) according to the time of death and Vibrio loads.

RNA Extraction and Transcriptome Sequencing
Based on the death and survival data and Vibrio load data, all shrimp were divided into V. parahaemolyticus-resistant and susceptible individuals. The collected hemocytes of V. parahaemolyticus-resistant and susceptible shrimp were mixed, respectively. For population 1, three biological replicates were prepared for resistant and susceptible shrimp. The samples were designated as HcS1-1, HcS1-2, and HcS1-3 for resistant shrimp and  designated as HcD1-1, HcD1-2, and HcD1-3 for susceptible shrimp. For population 2,  four replicates were prepared for resistant and susceptible shrimp. The samples were  designated as HcS2-1, HcS2-2, HcS2-3, and HcS2-4 for resistant shrimp and designated as  HcD2-1, HcD2-2, HcD2-3, and HcD2-4 for susceptible shrimp. The total RNA of hemocytes was extracted by RNAiso Plus (TaKaRa, Kyoto, Japan) according to the manufacturers' protocols. The total RNA was qualified by 1% agarose gel electrophoresis and quantified by NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA). About 1 µg total RNA from each sample was pre-treated with gDNA Eraser to remove genomic DNA according to the manufactures' instruction (TaKaRa, Kyoto, Japan). Fragmentation buffer was used to break the mRNA into short fragments. The first cDNA strand was synthesized using random hexamers, followed by the addition of buffer, dNTPs, RNase H, and DNA polymerase I to synthesize the second cDNA strand. Poly(A) was added to connect to the sequencing adaptor. Finally, the Illumina HiSeq2500 platform was used to sequence the library at Guangzhou Gene Denovo Biotech Co., Ltd. (Guangzhou, China).

Clean Data Mapping and Annotations
In order to ensure the data quality, the clean reads obtained from the preliminary filtering were subjected to more stringent filtering. High-quality clean reads were obtained by eliminating all the reads containing the A base and more than 10% of the reads containing N. The clean reads were mapped to the reference genome of L. vannamei by HISAT2.2.4 [12]. The mapped reads were assembled using StringTie v 1.3.1 [13] in a reference-based approach. The reconstructed transcripts were aligned to the reference genome. Annotation was carried out by blast against NR, KEGG, COG, and Swiss-Prot databases.

Differential Expression and Enrichment Analysis
The DEGs in hemocytes between V. parahaemolyticus-resistant and susceptible shrimp were identified by calculating the FPKM value. DEGs were analyzed using the edgeR package (http://www.rproject.org/, 10 January 2023), with the parameter of FDR < 0.05 and the absolute fold change ≥ 2. The correlation analysis between samples was carried out among samples to ensure the reliability of experimental data.
The GO function and KEGG pathway enrichment analyses of DEGs were performed using the online OmicShare tools (http://www.omicshare.com/tools, 12 January 2023). The Q value < 0.05 was considered statistically significant.

Quantitative Real-Time PCR
The validation was conducted by qRT-PCR using the RNA prepared for RNA-seq. The first strand of cDNA was reverse transcribed from 1 µg RNA with PrimeScript™ RT Reagent Kit with gDNA Eraser (TaKaRa, Kyoto, Japan). The products were diluted by 20-fold with nuclease-free water. Primers for qRT-PCR (Table S1) were designed for six DEGs from population 1 and all five DEGs from population 2 using Primer Premier 5.0 software. Additionally, 18s rRNA was used as an internal control to standardize the expression levels. The qRT-PCR was performed using THUNDERBIRD ® SYBR ® qPCR Mix on the Eppendorf Mastercycler ep realplex (Eppendorf, Germany) in a total volume of 10 µL, containing 5 µL qPCR mix, 1 µL diluted cDNA, 0.3 µL each of forward and reverse primers, and 3.4 µL nuclease-free water. The amplification programs were set as follows: 95 • C for 2 min; 40 cycles of 95 • C for 15 s, annealing temperature for 15 s, and 72 • C for 30 s; and followed by a melting curve. Three technical replicates were set for each sample. The expression analysis was calculated by 2 ∆∆ct and expression profiles of DEGs were shown with log2 fold change values.

The Loads of V. parahaemolyticus in Hepatopancreas of Infected Shrimp
Detection of Vibrio pathogens showed that there was no strain on the TCBS plate ( Figure 1A), suggesting that the shrimp used in experiments were free of the specific pathogens. The accumulative mortality rates for two populations were 71.25% and 48.57%, respectively ( Figure 1B). The standard curve equation of PirA Vp was established as y = −3.2029x + 40.884, where x represented the log value of Vibrio copy number per µL standard plasmid PirA, and y represented the Ct value. The correlation coefficient was 0.9652, which showed a good linear relationship between the concentration of standard substance and the Ct value ( Figure 1C).
showed a good linear relationship between the concentration of standard substance and the Ct value ( Figure 1C).
As shown in Figure 1D, the loads of V. parahaemolyticus in hepatopancreas of shrimp died at 9 hpi and 12 hpi were very low, which were 5.07 and 59.99 copies per ng hepatopancreas DNA. The loads of V. parahaemolyticus were high in hepatopancreas of shrimp died at 15 hpi, 18 hpi, 24 hpi, and 36 hpi, all of which were more than 1000 copies per ng of hepatopancreas DNA. The load of V. parahaemolyticus in hepatopancreas of survival shrimp (still alive at 48 hpi) was 10.57 copies per ng of hepatopancreas DNA. According to the death time and the load of V. parahaemolyticus in hepatopancreas, shrimp that died at 9 hpi and 12 hpi were deemed as susceptible individuals, and those with survival after 48 hpi were deemed as resistant individuals. Hemocytes from V. parahaemolyticus-resistant and susceptible shrimp in two populations were used for independent comparative transcriptome analysis.  As shown in Figure 1D, the loads of V. parahaemolyticus in hepatopancreas of shrimp died at 9 hpi and 12 hpi were very low, which were 5.07 and 59.99 copies per ng hepatopancreas DNA. The loads of V. parahaemolyticus were high in hepatopancreas of shrimp died at 15 hpi, 18 hpi, 24 hpi, and 36 hpi, all of which were more than 1000 copies per ng of hepatopancreas DNA. The load of V. parahaemolyticus in hepatopancreas of survival shrimp (still alive at 48 hpi) was 10.57 copies per ng of hepatopancreas DNA. According to the death time and the load of V. parahaemolyticus in hepatopancreas, shrimp that died at 9 hpi and 12 hpi were deemed as susceptible individuals, and those with survival after 48 hpi were deemed as resistant individuals. Hemocytes from V. parahaemolyticus-resistant and susceptible shrimp in two populations were used for independent comparative transcriptome analysis.

The Correlation of the Transcriptome Data from All Samples
The relationship of all samples from two populations were analyzed by PCA and sample clustering methods. The results showed that samples from the same population had a close relationship (Figure 2A,B). However, the difference between V. parahaemolyticus-resistant and susceptible shrimp was small in each population ( Figure 2C,D). The results suggest that shrimp from the same population, even with different V. para-Biology 2023, 12, 977 6 of 13 haemolyticus-resistant abilities, have more similar transcriptome profiles than those from different populations.
The relationship of all samples from two populations were analyzed by PCA and sample clustering methods. The results showed that samples from the same population had a close relationship (Figure 2A,B). However, the difference between V. parahaemolyticus-resistant and susceptible shrimp was small in each population ( Figure 2C,D). The results suggest that shrimp from the same population, even with different V. parahaemolyticus-resistant abilities, have more similar transcriptome profiles than those from different populations.

DEGs between Two Populations and between V. parahaemolyticus Resistant and Susceptible Shrimp
Corresponding to the sample relationship, much more DEGs were identified between two populations than those between V. parahaemolyticus-resistant and susceptible shrimp in each population. Between two populations, there were a total of 984 DEGs, including 496 up-regulated DEGs and 488 down-regulated DEGs in the hemocytes from population 2 ( Figure 3A). There were 678 DEGs between the V. parahaemolyticus-resistant shrimp from population 1 and 2 ( Figure 3B, HcS1-vs-HcS2). There were 342 DEGs between the V. parahaemolyticus-susceptible shrimp from population 1 and 2 ( Figure 3B, HcD1-vs-HcD2).
A total of 121 DEGs were overlapping between HcS1-vs-HcS2 and HcD1-vs-HcD2 ( Figure 3C). The expression trends of these DEGs showed obvious differences between

DEGs between Two Populations and between V. parahaemolyticus Resistant and Susceptible Shrimp
Corresponding to the sample relationship, much more DEGs were identified between two populations than those between V. parahaemolyticus-resistant and susceptible shrimp in each population. Between two populations, there were a total of 984 DEGs, including 496 up-regulated DEGs and 488 down-regulated DEGs in the hemocytes from population 2 ( Figure 3A). There were 678 DEGs between the V. parahaemolyticus-resistant shrimp from population 1 and 2 ( Figure 3B, HcS1-vs-HcS2). There were 342 DEGs between the V. parahaemolyticus-susceptible shrimp from population 1 and 2 ( Figure 3B, HcD1-vs-HcD2). vs-HcD1 and HcS2-vs-HcD2 ( Figure 3C), suggesting that distinct genes contribute to the disease-resistant ability of shrimp against V. parahaemolyticus infection in two populations.

Enriched GO Items and KEGG Pathways between V. parahaemolyticus-Resistant and Susceptible Shrimp
GO and KEGG enrichment results showed that some DEGs in population 1 were enriched in GO items and KEGG pathways. The enrichment items in cell composition were mainly "actin cytoskeleton", "myosin complex", "cytoskeleton", etc. The enriched items in molecular function were "motor activity", "nucleoside-triphosphatase activity", "pyrophosphatase activity", etc. ( Figure 4A). The enriched KEGG pathways included "viral myocarditis", "hypertrophic cardiomyopathy", "dilated cardiomyopathy" and "thyroid hormone signaling pathway" (Figure 4B). There were no GO item and KEGG pathway enriched for DEGs in population 2 because only five DEGs were identified. A total of 121 DEGs were overlapping between HcS1-vs-HcS2 and HcD1-vs-HcD2 ( Figure 3C). The expression trends of these DEGs showed obvious differences between two populations ( Figure 3D). However, there were no overlapping DEGs between HcS1vs-HcD1 and HcS2-vs-HcD2 ( Figure 3C), suggesting that distinct genes contribute to the disease-resistant ability of shrimp against V. parahaemolyticus infection in two populations.

Enriched GO Items and KEGG Pathways between V. parahaemolyticus-Resistant and Susceptible Shrimp
GO and KEGG enrichment results showed that some DEGs in population 1 were enriched in GO items and KEGG pathways. The enrichment items in cell composition were mainly "actin cytoskeleton", "myosin complex", "cytoskeleton", etc. The enriched items in molecular function were "motor activity", "nucleoside-triphosphatase activity", "pyrophosphatase activity", etc. ( Figure 4A). The enriched KEGG pathways included "viral myocarditis", "hypertrophic cardiomyopathy", "dilated cardiomyopathy" and "thyroid hormone signaling pathway" (Figure 4B). There were no GO item and KEGG pathway enriched for DEGs in population 2 because only five DEGs were identified.

Detailed Analysis of DEGs between V. parahaemolyticus-Resistant and Susceptible Shrimp
Among the 31 DEGs in population 1, 19 had functional annotations ( Table 1). Most of these DEGs exhibited higher expression levels in V. parahaemolyticus-susceptible shrimp. DEGs encoded cytoskeleton-related proteins, including troponin, myosin, septin-4, and actin, as well as immune and signal transduction-related genes, such as immunityassociated nucleotide-binding protein 13-like, E3 ubiquitin-protein ligase TRIM32, and MAM and LDL-receptor class A domain-containing protein 2-like, which had higher expression levels in V. parahaemolyticus susceptible shrimp. Sugar metabolism-related genes including glyceraldehyde-3-phosphate-dehydrogenase and alpha-(1,6)-fucosyltransferase-like were also highly expressed in V. parahaemolyticus-susceptible shrimp.

Detailed Analysis of DEGs between V. parahaemolyticus-Resistant and Susceptible Shrimp
Among the 31 DEGs in population 1, 19 had functional annotations ( Table 1). Most of these DEGs exhibited higher expression levels in V. parahaemolyticus-susceptible shrimp. DEGs encoded cytoskeleton-related proteins, including troponin, myosin, septin-4, and actin, as well as immune and signal transduction-related genes, such as immunity-associated nucleotide-binding protein 13-like, E3 ubiquitin-protein ligase TRIM32, and MAM and LDL-receptor class A domain-containing protein 2-like, which had higher expression levels in V. parahaemolyticus susceptible shrimp. Sugar metabolism-related genes including glyceraldehyde-3-phosphate-dehydrogenase and alpha-(1,6)-fucosyltransferase-like were also highly expressed in V. parahaemolyticus-susceptible shrimp. In population 2, one cytoskeleton-related genes-encoding tubulin α-3 was also highly expressed in V. parahaemolyticus-susceptible shrimp. The other four DEGs, including dihydropyrimidinase-like isoform X3, prohibitin, iroquois-class homeodomain protein IRX-2-like, and diacylglycerol kinase 1, were all highly expressed in V. parahaemolyticus-resistant shrimp ( Table 2). In order to validate the expression of DEGs in the transcriptome, we compared the expression trends from transcriptome data and qRT-PCR data for six DEGs from population 1 and five DEGs from population 2. The results showed that all genes had similar expression changes between transcriptome data and qRT-PCR data ( Figure 5), suggesting that the differential expression analysis in transcriptome data was credible. In population 2, one cytoskeleton-related genes-encoding tubulin α-3 was also highly expressed in V. parahaemolyticus-susceptible shrimp. The other four DEGs, including dihydropyrimidinase-like isoform X3, prohibitin, iroquois-class homeodomain protein IRX-2-like, and diacylglycerol kinase 1, were all highly expressed in V. parahaemolyticus-resistant shrimp ( Table 2). In order to validate the expression of DEGs in the transcriptome, we compared the expression trends from transcriptome data and qRT-PCR data for six DEGs from population 1 and five DEGs from population 2. The results showed that all genes had similar expression changes between transcriptome data and qRT-PCR data ( Figure 5), suggesting that the differential expression analysis in transcriptome data was credible.

Discussion
In order to screen genes in shrimp hemocytes related to V. parahaemolyticus resistance, the present study established a method to collect hemocytes from shrimp with different

Discussion
In order to screen genes in shrimp hemocytes related to V. parahaemolyticus resistance, the present study established a method to collect hemocytes from shrimp with different resistant abilities against V. parahaemolyticus infection. The results showed that shrimp more susceptible to V. parahaemolyticus infection died early with low pathogen load. When shrimp exhibited higher resistant ability against V. parahaemolyticus infection, they could suffer from higher pathogen load and stayed alive for a longer time. Shrimp with the highest resistant ability against V. parahaemolyticus infection could clear the pathogen in hepatopancreas after a longer time (48 h) post infection. Shrimp families with distinct resistances against specific pathogens were usually constructed, while individuals in the same family always exhibited different disease-resistant abilities. The method used in the present study could collect hemocytes from shrimp before pathogen infection with clear pathogen resistance information.
Free of Vibrio species ensured that the transcriptome profiles of shrimp hemocytes could represent the background expression levels of genes that were not influenced by specific pathogens. Immune challenge test showed that two populations had distinct resistant abilities against V. parahaemolyticus infection. This might be attributed to the different genetic backgrounds and animal sizes of the two populations. Nevertheless, different individuals in each population exhibited various disease-resistant abilities. Therefore, two independent transcriptome analyses on hemocytes of shrimp with different V. parahaemolyticus-resistant abilities were carried out to identify common genes related to the biological trait. However, the identified DEGs in two populations shared no overlap. The result was comparable with some transcriptome studies on shrimp hepatopancreas, which also identified rare common genes related to V. parahaemolyticus resistance [9,10]. It might be caused by the different genetic backgrounds of shrimp. In other words, different populations might employ distinct genes to resist V. parahaemolyticus infection. This brings difficulties for the identification of key genes or regulatory processes related to specific pathogens.
Although there was no overlap between DEGs identified from the two populations, many DEGs were reported relevant to immunity or pathogen infection. The cytoskeleton has four main components including actin, microtubules, intermediate filaments, and septins, which all have important functions in cell-autonomous immunity [14]. An E3 ubiquitin-protein ligase TRIM32 gene was crucial under oxidase stress and during Vibrio infection in L. vannamei [15]. In intestine, the MAM and LDL-receptor class A domaincontaining protein could regulate enterohepatic bile acid signaling, while bile acids are vital factors in mucosal immunity and inflammation [16,17]. Glyceraldehyde-3-phosphatedehydrogenase is the key enzyme in the glycolytic pathway, which benefits viral replication during WSSV in shrimp [18]. These genes were differentially expressed in hemocytes between V. parahaemolyticus-resistant and susceptible shrimp in population 1, indicating that they might contribute to the shrimp resistance against V. parahaemolyticus infection.
Even if the DEGs in population 2 were totally different than those in population 1, they were also related to immunity and pathogen infection. In the central nervous system, dihydropyrimidinase-like 3 could regulate the inflammatory response of activated immune cells microglia, which respond to infection and inflammation by producing cytokines and phagocytosing cell debris and pathogens [19]. In Fenneropenaeus chinensis, two prohibiting genes played important roles in hemocytes during WSSV infection [20]. Diacylglycerol kinase 1 is responsible for the ATP-dependent phosphorylation of diacylglycerol to phosphatidic acid, both of which are crucial second messengers involved in many immune processes such as regulating Toll-like receptor-induced cytokine production [21]. Tubulin alpha-3 is the major constituent of microtubules, one component of cytoskeleton. These genes might be responsible for variable resistant abilities of shrimp against V. parahaemolyticus infection in population 2.
Although no common genes related to V. parahaemolyticus resistance were identified in hemocytes from two populations, genes in specific GO item-like cytoskeleton could be found in DEGs from both transcriptome data. It means that the key processes or genes related to V. parahaemolyticus resistance could be obtained by comprehensive analysis of multiple omics data information. Therefore, a strategy was proposed here to approach this purpose ( Figure 6). In this proposal, individuals from the same full sibling family, which have relatively consistent genetic background, will be used to construct V. parahaemolyticusresistant and susceptible hemocytes samples based on the method established in the present study. Differentially expressed genes (Gene Set, GS), as well as genetic makers (GM), will be identified after comparative analysis between V. parahaemolyticus-resistant and susceptible samples from each family. The GS and GM information from each family will be pooled respectively to obtain the GS union and GM union. The GS union will be put for GO and KEGG enrichment analyses to identify key processes and genes related to V. parahaemolyticus resistance. Further association analysis between these identified key processes and genes and GM union will help to obtain useful genetic markers. The information will provide a basis for understanding regulatory mechanisms of shrimp defense against V. parahaemolyticus infection and for disease-resistance breeding. this purpose ( Figure 6). In this proposal, individuals from the same full sibling family, which have relatively consistent genetic background, will be used to construct V. parahaemolyticus-resistant and susceptible hemocytes samples based on the method established in the present study. Differentially expressed genes (Gene Set, GS), as well as genetic makers (GM), will be identified after comparative analysis between V. parahaemolyticus-resistant and susceptible samples from each family. The GS and GM information from each family will be pooled respectively to obtain the GS union and GM union. The GS union will be put for GO and KEGG enrichment analyses to identify key processes and genes related to V. parahaemolyticus resistance. Further association analysis between these identified key processes and genes and GM union will help to obtain useful genetic markers. The information will provide a basis for understanding regulatory mechanisms of shrimp defense against V. parahaemolyticus infection and for disease-resistance breeding. The proposed strategy was built based on the transcriptome data from two independent shrimp populations. The distinct genetic backgrounds among individuals might limit the identification of DEGs between V. parahaemolyticus-resistant and susceptible shrimp. Therefore, we suggested to identify DEGs using shrimp from the same full sibling family in the proposed strategy. A further investigation based on this proposal will be implemented to verify it.

Conclusions
In the present study, we identified genes related to V. parahaemolyticus resistance from hemocytes of two independent shrimp populations. Although no common genes existed in identified DEGs from two populations, there were some DEGs involved in the same biological processes such as cytoskeleton and metabolism. The data suggest that different genes might contribute to the resistant abilities in different shrimp populations, while these genes probably play functions in specific biological processes. Therefore, a strategy The proposed strategy was built based on the transcriptome data from two independent shrimp populations. The distinct genetic backgrounds among individuals might limit the identification of DEGs between V. parahaemolyticus-resistant and susceptible shrimp. Therefore, we suggested to identify DEGs using shrimp from the same full sibling family in the proposed strategy. A further investigation based on this proposal will be implemented to verify it.

Conclusions
In the present study, we identified genes related to V. parahaemolyticus resistance from hemocytes of two independent shrimp populations. Although no common genes existed in identified DEGs from two populations, there were some DEGs involved in the same biological processes such as cytoskeleton and metabolism. The data suggest that different genes might contribute to the resistant abilities in different shrimp populations, while these genes probably play functions in specific biological processes. Therefore, a strategy was proposed based on the present data and will be verified to identify key processes and genes related to disease resistance.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/biology12070977/s1, Table S1: Information of primers used in the present study. Table S2: The identified DEGs from population 1. Table S3: The identified DEGs from population 2.
Author Contributions: Conceptualization, S.L. and F.L.; methodology, K.Z. and W.D.; statistical analysis, S.L. and K.Z.; writing-original draft preparation, S.L. and K.Z.; writing-review and editing, S.L. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement:
This study used shrimp as experimental animals, which are not endangered invertebrates. In addition, there is no genetically modified organism used in the study. According to the national regulation (Fisheries Law of the People s Republic of China), no permission is required to collect the animals and no formal ethics approval is required for this study.

Informed Consent Statement: Not applicable.
Data Availability Statement: The original contributions presented in the study are included in the article/Supplementary Material. The raw data of RNA-Seq have been deposited to NCBI data bank with the accession number PRJNA974583. Further inquiries can be directed to the corresponding authors.

Acknowledgments:
We would like to thank the research team for their preparation of the materials for RNA-seq and technical support. Thanks for the data service provided by the Oceanographic Data Center, Chinese Academy of Sciences (CASODC).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: AHPND acute hepatopancreatic necrosis disease