Comparison of Gene Expression Between Resistant and Susceptible Families Against VPAHPND and Identification of Biomarkers Used for Resistance Evaluation in Litopenaeus vannamei

Acute hepatopancreatic necrosis disease (AHPND) has caused a heavy loss to shrimp aquaculture since its outbreak. Vibrio parahaemolyticus (VPAHPND) is regarded as one of the main pathogens that caused AHPND in the Pacific white shrimp Litopenaeus vannamei. In order to learn more about the mechanism of resistance to AHPND, the resistant and susceptible shrimp families were obtained through genetic breeding, and comparative transcriptome approach was used to analyze the gene expression patterns between resistant and susceptible families. A total of 95 families were subjected to VPAHPND challenge test, and significant variations in the resistance of these families were observed. Three pairs of resistant and susceptible families were selected for transcriptome sequencing. A total of 489 differentially expressed genes (DEGs) that presented in at least two pairwise comparisons were screened, including 196 DEGs highly expressed in the susceptible families and 293 DEGs in the resistant families. Among these DEGs, 16 genes demonstrated significant difference in all three pairwise comparisons. Gene set enrichment analysis (GSEA) of all 27,331 expressed genes indicated that some energy metabolism processes were enriched in the resistant families, while signal transduction and immune system were enriched in the susceptible families. A total of 32 DEGs were further confirmed in the offspring of the detected families, among which 19 genes were successfully verified. The identified genes in this study will be useful for clarifying the genetic mechanism of shrimp resistance against Vibrio and will further provide molecular markers for evaluating the disease resistance of shrimp in the breeding program.


INTRODUCTION
Litopenaeus vannamei is a commercially important aquaculture species, making up about 85% of total shrimp production in China (Qin et al., 2018). However, the shrimp aquaculture industry is continuously affected by the outbreak of viral and bacterial diseases, which have caused mass mortality and considerable economic losses. Vibrio parahaemolyticus carrying the PirA and PirB toxin genes in its plasmid (VP AHPND ) is one of the most destructive pathogens in shrimp aquaculture, causing acute hepatopancreatic necrosis disease (AHPND) or early mortality syndrome (EMS) in L. vannamei (Lee et al., 2015). VP AHPND also caused AHPND in Penaeus monodon and Exopalaemon carinicauda (Soonthornchai et al., 2016;Ge et al., 2018). Therefore, prevention and control of AHPND are urgently needed in shrimp aquaculture.
Genetic selective breeding of disease resistance broodstock is a feasible and sustainable approach for the disease control. It has been proved to be efficient in controlling Taura syndrome virus (TSV). A disease-resistant line of L. vannamei against TSV had been established, and 18.4% increase in survival rate against TSV infection was obtained after one generation selection (Argue et al., 2002). For the selection of White spot syndrome virus (WSSV) resistance in L. vannamei, it showed that the average survival rates of generations G2 to G5 were 5.57%, 7.78%, 9.52%, and 13.79%, respectively (Huang et al., 2012). After three successive generation selection, the survival rates of E. carinicauda to VP AHPND increased from 26.67% to 36.67% (Ge et al., 2018).
With the development of molecular biology, genomics approach offers a new possibility for accelerating the genetic selection process . Identification of the major genes associated with disease resistance is the first step for marker-assisted selection (MAS) or gene-assisted selection (GAS) (Arora et al., 2019). So far, several genes associated with disease resistance have been reported in aquatic animals. Polymorphism of LvALF and TRAF6 was reported to be associated with the resistance to WSSV in shrimp (Wang et al., 2011;Liu et al., 2014b;Liu et al., 2014a;Yu et al., 2017;Zhang Q. et al., 2019). A major quantitative trait locus (QTL) for the resistance of Atlantic salmon (Salmo salar) against infectious pancreatic necrosis (IPN) was discovered, which has been already applied in marker-assisted breeding of IPN-resistant fish (Houston et al., 2008;Moen et al., 2009;Woldemariam et al., 2020). However, knowledge about the resistance to AHPND in shrimp is still poorly understood. A comparison of individuals or families with significant phenotype difference by transcriptome sequencing is an efficient way for screening the trait-associated genes. Based on the transcriptome data, myosin, myosin heavy chain, and chitinase were proved to be related to growth performance in L. vannamei (Santos et al., 2018). Transcriptome comparison between the families with high growth rate and low growth rate also illustrated that the genes related to cuticle, chitin, and muscle proteins were upregulated exclusively in higher growth families (Santos et al., 2021). Besides, the gene profiles of the Vibrioresistant and Vibrio-susceptible Meretrix petechialis families were analyzed, and several genes such as Big-Def, CTL9, and Bax were identified as candidate resistance-associated genes (Jiang et al., 2017).
In our previous work, we have carried out systematic family selection for the resistance trait to VP AHPND in L. vannamei. A total of 95 families were subjected to VP AHPND challenge test, and the results showed that significant resistance variations existed between different families. In the present study, we selected three AHPND resistant and three susceptible families for transcriptome sequencing and explored the differentially expressed genes (DEGs) in resistant and susceptible families. Several genes related to the resistance of shrimp against AHPND were identified. These genes might be developed as effective molecular markers for evaluating the disease resistance of shrimp, which could facilitate molecular marker-assisted breeding of shrimp.

Selection Resistant and Susceptible Families Against Acute Hepatopancreatic Necrosis Disease
The full-sib families of L. vannamei were produced and stocked separately in Hainan Grand Suntop Ocean Breeding Co., Ltd, Wenchang, China. In order to identify the resistance of AHPND, the shrimp families were challenged by VP AHPND each year, and the families with high survival rate were mated to generate the next-generation families. In 2018, 95 full-sib families with an average body weight of 2.10 g were selected for evaluation of resistance. For the challenge experiment, VP AHPND was prepared according to the method described by Zhang Q. et al. (2019). About 100 healthy shrimp from each family were subjected to VP AHPND immersion infection. Before the experiment, the shrimp were kept in the aquarium at a temperature of 26°C ± 1°C for 3 days to acclimate them to the culture conditions. Then, the VP AHPND were added to the aquarium to make the concentration of VP AHPND as 5 × 10 6 CFU/ml. Then, the dead shrimp were checked every 2 h. The mortality of each family was recorded for 72 h. The survival rates of the tested families are presented in Figure 1.
For each family, the cephalothoraxes of nine individuals were collected, and three individuals were mixed together as one sample, so each family contained three samples. All cephalothorax samples of the six families were rapidly frozen in liquid nitrogen and stored at −80°C for further transcriptome sequencing.

Total RNA Extraction, Library Construction, and Sequencing
Total RNA of the cephalothorax was extracted using RNAiso Plus (Takara, Japan) according to the manufacturer's instructions. RNA purity and concentration were evaluated by electrophoresis on 1% agarose gel and quantified by NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Subsequently, mRNA was enriched by Oligo (dT) beads. Then the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse-transcribed into the first-strand cDNA with random primers. Secondstrand cDNA was synthesized by DNA polymerase I, RNase H, dNTP, and buffer. Then the cDNA fragments were purified with QiaQuick PCR extraction kit, end repaired, poly (A) added, and ligated to Illumina sequencing adapters. The suitable fragments were selected by agarose gel electrophoresis. After PCR amplification, the libraries were sequenced using Illumina HiSeqTM 2,500 by Genedenovo Biotechnology Co., Ltd (Guangzhou, China).

Bioinformatics Analysis
Clean reads were obtained by removing reads containing adapters or more than 10% of unknown nucleotides (N), low-quality reads, and rRNA from the raw data. Then the reads of each sample were mapped to reference genome  by TopHat2 (version 2.0.3.12) (Kim et al., 2013). The reconstruction of transcripts was carried out with software Cufflinks (Trapnell et al., 2012). All reconstructed transcripts were aligned to reference genome, and novel genes were aligned to databases including National Center for Biotechnology Information (NCBI) nonredundant (Nr) database, Swiss-Prot database, and Kyoto Encyclopedia of Genes and Genomes (KEGG) database to obtain protein functional annotation. DEGs between susceptible and resistant libraries were analyzed using the edgeR package (http://www.rproject.org/), with a fold change ≥2 and a false discovery rate (FDR) <0.05. DEGs were then subjected to enrichment analysis of Gene Ontology (GO) functions and KEGG pathways, which were performed using the OmicShare tools (http://www.omicshare.com/tools).

Gene Set Enrichment Analysis
As a complement to the differential expression analyses, gene set enrichment analysis (GSEA) for all 27,331 expressed genes was performed by GSEA software v4.10 (https://www.gsea-msigdb. org/gsea/downloads.jsp) (Subramanian et al., 2005). The submitted gene list was ranked by gene expression value using Signal2Noise method. GSEA was performed with default

Transcriptome Sequencing Data
An overview of sequencing and assembly of the L. vannamei transcriptome is shown in Table 2. A total of 823,161,050 raw reads were obtained, in which 416,197,128 reads are for the resistant families and 406,963,922 reads for the susceptible families. The raw sequencing data were uploaded to the NCBI with the accession numbers SRR15533118-SRR15533135. After low-quality reads were filtered out, a total of 96.55% and 96.85% reads were retained for the resistant and susceptible groups, respectively. After the reads were mapped to the reference genome, a total of 27,331 annotated genes were obtained, of which 2,344 (9.86%) were newly annotated genes.

Identification of Differentially Expressed Genes Between Resistant and Susceptible Families
To identify DEGs involved in VP AHPND resistance, we used FPKM value for comparing the expression levels between the resistant families and susceptible families. A total of 392 unigenes showed differential expression patterns between families VR4013 and VS3868, in which 287 unigenes were highly expressed in VS3868 and 105 unigenes were highly expressed in VR4013 (Figure 2). A total of 1,378 unigenes were differentially expressed between VR3837 and VS3879, including 773 unigenes highly expressed in VS3879 and 605 unigenes highly expressed in VR3837. A total of 2,183 unigenes were differentially expressed between VR4027 and VS3880, including 564 unigenes highly expressed in VS3880 and 1,619 unigenes highly expressed in VR4027. Six DEGs selected from each comparison group were validated by RT-qPCR. The results showed that all of them were consistent with the transcriptome data (Supplementary Figure S1). In order to analyze the function of DEGs between resistant and susceptible families, GO functional enrichment analysis was performed. Interestingly, we found that GO terms of DEGs were very similar among three pairs of families, VR4013-vs.-VS3868, VR3837-vs.-VS3879, and VR4027-vs.-VS3880, which is shown in Figure 3A. In the biological process category, single- organism process was the most enriched subclasses, followed by cellular process. In the cellular component category, cell and cell part were the two most enriched subclasses. While in the molecular function category, catalytic activity and binding were the two most enriched subclasses. KEGG pathway enrichment analysis of three pairwise comparisons is  presented in Figure 3B. In VR4013-vs.-VS3868, "Serotonergic synapse," "Arachidonic acid metabolism," and "PI3K-Akt signaling pathway" were the three most enriched pathways, and "Arachidonic acid metabolism" was also included in top 20 of enrichment pathways in VR3837-vs.-VS3879 and VR4027vs.-VS3880 (p < 0.05).

Differentially Expressed Genes Shared in More Than Two Pairwise Comparisons
Venn diagrams for DEGs of the three pairwise comparisons are shown in Figure 4. A total of 489 DEGs were shared in at least two comparisons, including 196 DEGs highly expressed in the susceptible families and 293 DEGs highly expressed in resistant   Table S1). A total of 16 DEGs were shared in three pairwise comparisons, among which the trophoblast glycoprotein, SE-cephalotoxin-like, peroxidasin-like protein, phosphoenolpyruvate carboxykinase, cytochrome P450, and serpin B6-like isoform X2 genes were identified (Table 3). In order to reduce the false positives, we considered these 489 DEGs as the candidate genes associated with the resistance of shrimp against to VP AHPND . These 489 DEGs were involved in various GO classifications. For the biological process-related genes, most were involved in "single-organism process," "cellular process," and "metabolic process." Most of the cellular component-related genes were associated with "cell," "cell part," and "organelle." And "catalytic activity" and "binding" in the molecular function ontology were the major enriched terms ( Figure 5). The results were consistent with GO functional enrichment analysis in Figure 3A.

Differentially Expressed Genes with High Expression in Resistant Families
There were 293 DEGs highly expressed in resistant families. A total of nine genes annotated as myosin were highly expressed, and the other genes were involved in energy metabolism such as trypsin, chymotrypsinogen A (ChyA), pancreatic lipase (PL), serine protease (SP), aminopeptidase, and phospholipase (Table 4). There were two genes encoding glyceraldehyde 3-phosphate dehydrogenase (GAPDH) and triosephosphate isomerase (TPI), which were involved in the glycolysis pathway (Eanes, 2011).

Differentially Expressed Genes with High Expression in Susceptible Families
A total of 196 DEGs highly expressed in the susceptible families were discovered. According to annotation and function, some immune-related genes were observed to be upregulated in the susceptible families ( Table 5). A total of five genes related to the prophenoloxidase (proPO) system, including C-type lectin, SP, and SP inhibitors (serpin and α-2-macroglobulin) showed high expression in the susceptible families. Several genes are related to metabolic process, including DBH-like monooxygenase protein 1, cytochrome b5, and cytochrome P450. Moreover, six genes annotated as β-arrestin were also highly expressed in the VP AHPND -susceptible families.

Gene Set Enrichment Analysis
By GSEA of GO gene set, there were a total of 265 significantly enriched gene sets between resistant and susceptible groups (nominal p-value < 0.05 and FDR < 0.25). The top 10 gene sets enriched in resistant families were all involved in energy metabolism process, while most gene sets enriched in the susceptible families were related to signal transduction and immunity, such as G protein-coupled receptor activity, G protein-coupled receptor signaling pathway, and pattern recognition receptor signaling pathway ( Table 6). Figure 6 shows GSEA of KEGG subset of canonical pathways. The most significantly enriched pathway in resistant families was oxidative phosphorylation (Figure 6A), which was also related to energy metabolism, supporting the results of GO gene sets. The most significantly enriched pathway in the susceptible families was JAK/STAT signaling pathway ( Figure 6B), including protein inhibitor of activated STAT, cytokine-inducible SH2containing protein, and tyrosine-protein phosphatase nonreceptor type 11 isoform X2.

Verification of Candidate Genes in the Descendant Families
A total of 32 genes were validated in descendant families. For the 16 DEGs shared by three pairwise comparisons, seven genes were verified successfully in the hepatopancreas, and eight genes were verified successfully in the stomach in their offspring ( Table 3, p < 0.05). For the other 16 DEGs shared by two pairwise comparisons, three genes were successfully verified in hepatopancreas (p < 0.05) ( Figure 7A), and six genes were verified successfully in the stomach (p < 0.05) ( Figure 7B). In summary, a total of 19 genes were successfully verified in the hepatopancreas or stomach, including five genes verified successfully in two tissues.  Note. EP, enrichment in phenotype, gene sets enriched in nine resistant (R) samples or nine susceptible (S) samples; Size, number of genes in the gene set; NES, normalized enrichment score; NOM p-value, nominal p-value, the statistical significance of the enrichment score; FDR, false discovery rate; GSEA, Gene set enrichment analysis; GO, Gene Ontology.
Frontiers in Genetics | www.frontiersin.org November 2021 | Volume 12 | Article 772442 9 DISCUSSIONS It is generally considered that the VP AHPND is an opportunistic pathogen. It is pathogenic to cultured shrimp at high concentration or under environmental deterioration condition (Hong et al., 2016). The physiology and health state of shrimp are closely related to disease resistance. It was estimated that the heritability of the resistance against VP AHPND in L. vannamei was near-to-moderate, which indicated that the resistance of shrimp against VP AHPND was hereditable . In this  . These data are expressed as the mean ± SD relative to the reference gene (18S rRNA). A statistically significant difference is indicated with * (p < 0.05) and ** (p < 0.01). Successfully verified genes are marked in red. Unsuccessfully verified genes are marked in black or without any mark.
Frontiers in Genetics | www.frontiersin.org November 2021 | Volume 12 | Article 772442 10 study, the resistant and susceptible families of L. vannamei were obtained through continuous challenge test, which provided reliable genetic material for analysis on the disease resistance of shrimp (Grimholt, et al., 2003;Wang et al., 2014). Understanding the gene expression character of diseaseresistant families and the susceptible families could provide useful information for genetic dissection of disease resistance of shrimp.
Recently, an increasing number of transcriptome studies related to AHPND have been performed in L. vannamei, most of which focused on the genes involved in the immune response of shrimp during VP AHPND infection (Qi et al., 2017;Maralit et al., 2018;Ng et al., 2018;Qin et al., 2018;Zheng et al., 2018). However, few works have focused on the correlation between the gene expression level and resistance phenotype of shrimp. In this study, we performed RNA-seq analyses to investigate the comparative expression profiles between resistant and susceptible families of shrimp. Gene expression profiles can underlie complex phenotype variations (Crawford and Powers, 1992). The gene expression variation is influenced by various genetic and environmental factors (Leder et al., 2014). Many studies proved that genetic factors influence gene expression in humans (Cheung and Spielman, 2002), mice (Sandberg et al., 2000), Drosophila (Jin et al., 2001), and yeast (Brem et al., 2002). To our knowledge, this is the first report about the comparison of the basal mRNA expression profiles from the omics level of AHPND-resistant and AHPND-susceptible families of L. vannamei. The gene expression levels influencing the phenotype variations could also be considered as molecular markers and could be used for genetic selection (Gilad et al., 2008). For aquaculture species, genes related to phenotype variations have been reported. The expression of BIRC7 was correlated with the survival in Vibrio challenge tests in clam M. petechialis, and gene expression variation of BIRC7 gene was heritable, indicating the feasibility of selective breeding by reliable genetically based markers (Jiang et al., 2018;Jiang et al., 2019).
In the present study, we found that DEGs highly expressed in the susceptible families of shrimp were mainly related to the immunity, which were in agreement with previous research results. It was reported that the basal expression levels of immune-related genes BIRC7 and NFIL3 were higher in Vibrio-susceptible clams (Jiang et al., 2018). In this study, several genes in proPO activation system showed a higher expression level in the susceptible families. A previous report showed that PPAE2 presented a higher expression level in the susceptible line of L. vannamei after AHPND challenge in comparison with the resistant line (Mai et al., 2021). This system is important for fighting against bacteria pathogens in penaeid shrimp (Amparyup, et al., 2013), and the increased activity of proPO system against Vibrio has been reported in Fenneropenaeus indicus (Sarathi et al., 2007), L. vannamei (Boonchuen et al., 2021), and Macrobrachium rosenbergii (Rao et al., 2015). The C-type lectin was also upregulated in the susceptible families of shrimp; the C-type lectin plays an important role in phagocytosis, melanization, respiratory burst, and coagulation; and it can also activate the proPO system (Junkunlo et al., 2012). Interestingly, C-type lectin 1-like and Crustin-P had significant higher expression levels in the AHPND-susceptible line of L. vannamei than resistant line at 24 h post-infection (Mai et al., 2021). In Litopenaeus stylirostris, the basal expression level of antimicrobial peptide was significantly higher in Vibrio-susceptible shrimp line (Lorgeril et al., 2008). In addition, β-arrestin also played a vital role in the antibacterial immunity of shrimp (Jiang et al., 2013;Sun et al., 2016). The GSEA also illustrated that partial immune-related genes were upregulated in the susceptible families; the JAK/STAT signaling pathway, which was an important signal transduction pathway regulating the immune response in invertebrates, was significantly enriched in the susceptible families (Yu et al., 2017).
Besides, the genes in immunity showed differential expression patterns between susceptible and resistant families; another major part of the DEGs was related to myosin and metabolism. In this study, myosin and many other genes related to metabolism were more active in resistant families. Myosin and actin play a diverse role in a wide range of functions such as cytoskeleton, muscle contraction, and immune response (Marston 2018;Sitbon et al., 2019). It was already reported that myosin light chain was related to the phagocytosis against invading pathogens in Penaeus japonicus, and the transcription level of myosin in WSSVresistant shrimp was nearly two times higher than that in the control shrimp (Han et al., 2010). After pathogen infection, myosin and actin were significantly upregulated in shrimp (Shi et al., 2018;Ren et al., 2019). As for metabolism, we found that enzymes like trypsin, ChyA, PL, and SP and glycolysis pathway including GAPDH and TPI were highly expressed in the resistant families (Table 4). Meanwhile, the top 10 gene sets enriched in resistant families in the GSEA were all involved in energy metabolism process. The finding was consistent with the previous report that SP and ChyB had a significantly higher expression level in resistant shrimp line during the AHPND infection (Mai et al., 2021). Several studies have indicated that metabolic processes such as lipid metabolism and carbohydrate metabolic process in shrimp were greatly affected during AHPND infection (Velazquez-Lizarraga et al., 2019;Kumar et al., 2021). Taken together, the activated myosin, actin, and energy metabolism might indicate that the shrimp were healthier, which led to higher resistance of shrimp to disease.
After validation, several genes showed the same expression pattern in the offspring of susceptible and resistant families. These genes are possible to be developed as biomarkers for disease resistance of shrimp. It was already reported that gene expression profile could be used as an indicator for disease resistance trait. For example, Bsr-d1 RNA level in susceptible rice strain was twofold higher than that in resistant rice strain, and it has been further proved that inhibiting the expression of bsr-d1 could increase the rice resistance against blast infection . In E. carinicauda, eight immune-related genes were suggested as potential disease-resistant parameters for evaluating the physiological status and disease-resistant capability of shrimp when infected with VP AHPND (Ge et al., 2018). In this study, the 19 genes successfully verified in their descendant families were expected to be developed as biomarkers of shrimp resistance against Vibrio. Therefore, apart from sib challenge testing and molecular marker-assisted breeding, the gene expression level of these 19 genes could also be used as molecular markers for accelerating the breeding of disease-resistant varieties in L. vannamei.
In summary, this study integrated the VP AHPND -resistance phenotype variation and gene expression profiles to identify the genes related to disease resistance of shrimp. A total of 489 DEGs were identified between the resistant and susceptible families, and they were considered to be associated with the ability of VP AHPND resistance in L. vannamei. Gene annotation and enrichment analysis revealed that the immune response and energy metabolism could influence resistance of shrimp against VP AHPND . The obtained data provide a fundamental basis for clarifying the genetic mechanism of resistance to bacterial pathogen, and the identified disease resistance genes of shrimp could accelerate the genetic breeding in shrimp aquaculture.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
QZ, YY, and FL contributed to the conception and design of the study. ZL performed the statistical analysis. QZ wrote the first draft of the article. QZ, YY, JX, and FL wrote sections of the article. All authors contributed to article revision and read and approved the submitted version.